动量守恒和粒子数守恒: stock mkt!
第22 卷第3 期
2007 年3 月
地球科学进展
ADVANCES IN EARTH SCIENCE
Vol.22 No.3
Mar.，2007
文章编号：1001-8166 （2007 ）03- -12
格子玻尔兹曼方法及其在大气湍流研究中的应用
*
程雪玲，胡 非，赵松年，姜金华
（中国科学院大气物理研究所大气边界层物理和大气化学国家重点实验室，北京 100029 ）
摘 要：文章的目的是对格子玻尔兹曼方法进行系统的介绍，格子玻尔兹曼方法（ Lattice Boltzm ann
Method ）的出现直接来源于20 世纪60 年代的元胞自动机（ Cellular Autom ata ）思想，而这一方法用
于解决流动现象时，又可以追溯到19 世纪的分子运动论，求解的是 Boltzm ann 提出的玻尔兹曼输运
方程，因此将这一方法称为格子玻尔兹曼方法，之前也被称为格子气自动机（ Lattice Gas Autom aton
）。该方法多用于研究复杂现象，如材料晶体凝聚时的生长过程、城市土地利用的演化等方面。
在20 世纪70 年代由 Hardy 、 Pom eau 和 Pazzis 建立了第一个用于研究流体运动的格子气自动机，此
后，这一方法被广泛用来模拟各种流动问题，诸如二相流、孔隙介质中的渗流等，并根据这一方法开
发了相应的商业软件 Pow er Flow 。同时，格子玻尔兹曼方法由于其在微观水平描述运动的特点，成
为研究湍流的一个很好的数值计算工具，特别是用其进行直接数值模拟（ DNS ）计算，成为继传统的
差分法、有限体积法和谱方法之后的又一有力的手段。而作为大气运动的一个主要现象的大气湍
流，比普通湍流更加复杂，在这里着重介绍了大气湍流的特点和应用格子玻尔兹曼方法模拟湍流的
发展过程。
关 键 词：格子玻尔兹曼；元胞自动机；格子气；直接数值模拟；大气湍流
中图分类号： P425.2 文献标识码：A
1 引 言
元胞自动机是空间和时间都离散、物理参量只
取有限数值集的物理系统的理想化模型。元胞自动
机的思想首先是由J ohn von Neum ann 提出的，尽管
由于他在20 世纪40 年代末参与了第一台数字计算
机的设计，而使他的名字与当代串行计算机联系起
来，但他通过模仿人脑的并行化特点而提出的元胞
自动机思想，第一次人工构造出动力演化系统，为解
决复杂问题开创了新的方向。在20 世纪60 年代他
和 S. Ulam 构造了一个自复制系统，这个系统中元
胞的活动是同时进行的，所有元胞的演化规则是相
同的，随邻近元胞的状态而变化，这一系统用来模拟
生物的自复制过程。尽管这一结果连极原始的生命
形式也算不上，但却具有非常重大的意义。因为通
常认为，机器只能构造比自身简单的结构，而采用自
复制元胞自动机，却可以获得一种能产生新的、具有
同样复杂性和功能的“机器”。继 Von Neum ann 的
研究之后，另一些学者也遵循同样的研究路线，继续
研究这个问题。在1970 年，这方面的工作被 Conw -
ay 等加以发展，提出了著名的“生命游戏”（ Gam e of
Life ）模型，这是一个二维元胞自动机的实例。结果
表明，生命游戏机有出乎意料的丰富行为。而元胞
自动机真正被系统地作为动力系统来研究和应用当
推 S. W olfram ，他于1983 年在《现代物理评论》杂志
上发表了奠基性的论文“元胞自动机的统计力学”，
* 收稿日期：2006-11-27 ；修回日期：2007-01-11.
* 基金项目：国家自然科学基金项目“大气湍流能量级串机理及其格子气数值模拟的研究”（编号：40405004 ）；国家自然科学基金重点
项目“非均匀地表通量与大气边界层过程的研究”（编号：40233030 ）；国家自然科学基金项目“非均匀下垫面上中尺度通量
参数化的研究”（编号：40605006 ）资助.
作者简介：程雪玲（1971-），女，辽宁大连人，副研，主要从事大气边界层湍流理论以及湍流扩散的数值模拟等方面的研究.
E - m ail ： chengxl ＠ m ail. lapc. ac. cn
掀起了元胞自动机研究的新高潮［1］。
一般公认， Hardy 等［2］于20 世纪70 年代建立
的所谓 HPP 格子气模型实际上是元胞自动机，这个
模型是由简单、具有全离散动力特性的粒子构成，粒
子以保持动量守恒和粒子数守恒的方式，在二维方
形网格上运动和碰撞。
然而，出乎最初的预想结果，格子气模型并不优
于流体力学的传统数值方法。虽然这种方法在计算
流体运动时，由于其布尔（ Boole ）性质，可以不做什
么近似，如对于多体碰撞项中的多体相关函数可以
直接表达出来，而不用将其分解成单粒子密度函数
之积，因此在计算机上能够精确地进行模型的演化。
但是也正是由于其处理的是布尔量和碰撞规则的布
尔性质，使得一方面在对布尔量求统计平均时产生
噪音，另一方面用布尔性质对模拟的物理系统的参
数的调节欠灵活，如用碰撞规则来确定的粘滞性通
常都很高，且不可调节，因此造成实际研究的困难。
20 世纪80 年代末， Mc Nam ara 等［3］、 Higuerasi 等［4］
证明了把自动机的布尔动力学扩展到直接处理实数
（表示元胞具有指定状态概率）的优点。且 Mc Nam
ara 等［3］在他们的工作中首先研究了将元胞自动
机思想应用于求解玻尔兹曼方程的数值模拟，这也
就是现在广泛用来模拟各种流动的格子玻尔兹曼方
法。
早在18 ～19 世纪，随着欧洲的工业革命带动了
热机的发展，建立了热力学第一定律、热力学第二定
律，而对于后者的研究更是由于玻尔兹曼等人致力
于探究热力学底下的微观层次中的分子机制，而促
使统计物理学的蓬勃发展［5］。1857 年克劳修斯发
表了题为《论热运动形式》一文，首次引进了统计概
念.他指出单系统的宏观性质只是大量分子运动的
平均效果。1860 年麦克斯韦在《气体动力论的说
明》一文中，首先认识到分子的速度各不相同而得
到了速度分布律———麦克斯韦速度分布律。奥地利
物理学家玻尔兹曼把分子运动论理论的研究推向了
高峰。玻尔兹曼在接受了麦克斯韦速度分布律的规
律后，进而考虑到麦克斯韦的分布律，只反映了气体
平衡态下的情况，其中没有时间因子，他试图建立起
非平衡态分布的运动方程，在1872 年，他得出了速
度分布函数必须满足的非线性积分微分方程，即著
名的玻尔兹曼（ Boltzm ann ）方程。而玻尔兹曼基于
分子运动的统计特性研究宏观现象的思想，随着计
算机技术的发展，新的计算方法的出现，更是从热力
学研究领域向流体力学等更广阔的物理学科扩展，
我们这里所介绍的格子玻尔兹曼方法即是用于研究
流动现象的一例。
2 格子玻尔兹曼方法
2.1 玻尔兹曼方程
假设所考虑的系统是稀薄气体，可忽略三体及
三体以上的碰撞，而只考虑二体碰撞，且假设是弹性
碰撞。再假定气体是单一成分的单原子分子，即没
有源项，则粒子在相空间（珒r，珒v）附近的小体积元d珒rd珒v
中的粒子数fd珒rd珒v，它随时间的变化率5
5t
（fd珒rd珒v）由
以下部分组成：
（1） 粒子在坐标空间中运动造成的变化率［ 5
5t
（fd珒r d珒v］1，由于速度在珒v 附近d珒v 范围内分子的运动
所造成d珒r 体积元中粒子数的增加率为：－ dv ∮
σ
f珒v·
dσ，其中σ表示小体积的表面，dσ是小面积元，方
向沿表面外法向。利用 Gauss 定理有：
－ d珒vo ∮
σ
f珒v·dσ ＝－ dv ∫5
5r
·（f珒v）d珒r
由于珒v 和珒r 都是独立变量5
5r
·（f珒v） ＝珒v· 5
5r
，又考虑
到整个体积就是d珒r，便有：
［ 5
5t
（fd珒r d珒v）］1 ＝d珒vd珒r·珒v· 5
5珒r
（1）
（2） 粒子在速度空间中迁移造成的变化率［ 5
5t
（fd珒r d珒v）］2，这是因为在相空间中，坐标中包含有速
度坐标珒v，所以才有这一项，和上一项同样方法可以
简化为，
［ 5
5t
（fd珒rd珒v）］2 ＝d珒vd珒r 5
5珒v
·（f
珗F
m
） （2）
其中珗F 是作用于粒子的外力，m 是粒子质量。如果
5
5v
· 珗F ，则有，
［ 5
5t
（fd珒rd珒v）］2 ＝d珒vd珒r
珗F
m
· 5f
5v
（3）
（3） 粒子间相互碰撞造成的变化率5f
5tc
d珒rd珒v，设珒v 与
珒v2 是
碰
撞
前
两
分
子
的
速
度
，
碰
撞
后
速
度
分
别
为
珒
v′，
珒v′2
，则碰撞项可以写成：
（5f
5t
）c ＝∫dΩσ（｜ g ｜，θ） ｜ g ｜
［f（珒v′）f（珒v′2）－ f（珒v）f（珒v2）d珒v2 （4）
30 地球科学进展 第22 卷
其中，f（珒r，珒v，t）称为单粒子分布函数，g 是两个分子
的相对速度，θ为g 与质点系中心线的夹角，σ为微
分散射截面，Ω为固体角。由于有：
5
5t
（fd珒r d珒v） ＝［ 5
5t
（fd珒r d珒v］1 ＋［5
5t
（fd珒r d珒v］2
＋（5f
5t
）cd珒r d珒v （5）
将式（1）、（3）和（4）代入到式（5）中，得到：
5f
5t
＋珒v·5f
5珒r
＋
珗F
m
·5f
5珒v ＝∫dΩσ（｜ g ｜，θ）
｜ g ｜（ f′ f′2 － ff 2）d珒v2 （6）
这就是稀薄气体的输运方程，也称为玻尔兹曼输运
方程［6，7
］。
考虑系统不受外力作用的简单情形，此时玻尔
兹曼方程简化为：
5f
5t
＋珒v·èf ＝（5f
5t
）c （7）
但该方程当时无法求解，直到40 年后（1916 —
1917 ） S. Chapm an 和 D . Enskog 才第一次得出在一
般情况下求解的一种方法。这是一种将分布函数展
开成努森（ Knudsen ）数的幂级数的逐级近似方法。
他们假设系统虽未达到平衡，但在每一小局域的体
积中已离局部平衡不远。所谓局部平衡，是指分子
的速度服从 Maxwell- Boltzm ann 分布：
f0（珒r，珒v，t） ＝n（珒r，t） m
2πkB T（珒r，t
( )
）
3/2
exp －
m v2
2kB T（珒r，t
( )
）
（8）
其中n（珒r，珒v，t）＝ρ（珒r，t）/m 。
假设f可在f0 附近展开，
f ＝ f 0 ＋Knf 1 ＋Kn 2f2 ＋⋯ （9）
Kn ＝τ
Df/Dt
f n 1，为努森（ Knudsen ）数，表示2 个基
本的时间尺度（当以平均自由程和流体空间尺度表
示时为长度尺度）之比，也就是说：分子碰撞时间尺
度τ（或者与湍流涡的相互作用相关的特征碰撞时
间尺度τ turb ），和由于平均流的非均匀性形成的偏离
当地值的对流时间尺度之比， D / D t是实导数。
根据其一级近似所得的流体动力学方程，就是
通常适用于连续介质的 Navier- Stokes （NS）流体动力
学方程组。但是这之后，玻尔兹曼方程并未能象NS
方程那样，随着计算机技术的发展，在数值模拟方面
得到广泛的应用，直到20 世纪60 年代元胞自动机
理论的出现，在此基础上，科学家们成功地进行了首
批流体力学格子气模拟，为玻尔兹曼方程的格子化
提供了数学手段。
2.2 格子玻尔兹曼方法
Hardy 等所建立的 HPP 格子气模型的物理过程
实际上与玻尔兹曼方程所反映的物理过程相似，其
基本要素是：①空间离散，将流体存在的空间划分成
间距为1 的正方形网格（图1a）；②流体离散，将流
体想象成由许多只有质量（单位质量），没有体积的
微小粒子组成，每个粒子可向4 个方向之一运动；③
时间离散，有了流体及其存在的空间，它们的相互作
用就称为时间作用步，分为粒子在每个节点上的相
互碰撞和从一点到邻近点的移动（图1b）；④碰撞规
则，当且仅当只有两个粒子沿相反的方向到达某节
点时（对头碰撞），它们沿另外的两个方向离开该节
点，其他情形则直接穿透（图1c）；⑤ Pauli 不相容原
理，在同一时刻同一网格点上，每一个速度方向最多
只允许有一个粒子，它为 Boole 型表示法提供了可
能，可以用“1”表示有粒子，“0”表示没有粒子。
有了以上5 条要素，若再给定边界上的碰撞规
则和初始状态，系统就变成了一个离散化 Boole 化
（a） HPP 模型网格（b）时间步作用过程
（c）碰撞规则
图1 HPP 模型
Fig.1 HPP m odel
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 31
动力系统———格子气自动机［7］（ Lattice Gas Autom ata，
LGA ）；它可以随着分离的时间步而运行。
如前所述，如果我们考虑系统不受外力作用，简
化的玻尔兹曼方程为：
5f
5t
＋珒v·èf ＝（5f
5t
）c （10）
考虑我们所讨论的格子气模型，可以知道有2
个因素导致f 随时间变化：① 粒子从一个节点到另
一个节点的运动；②在网格节点上粒子间的相互碰
撞。实际上粒子间不断发生的碰撞，对f 的影响非
常复杂，我们只引入一个简单的碰撞项Ω，得到格子
玻尔兹曼方程：
f1（ x ＋ e i，t ＋1）－ fi（x，t） ＝Ωi （11）
下标i表示粒子的某个运动方向。
在格子玻尔兹曼方程中，如果直接计算碰撞项，
计算量相当巨大，对这个问题的解决是由 Higuera
及其合作者［4，8
］提出来的。早在1954 年， Bhatnager
等［9，1
0］曾假设，那些描述涡与涡之间相互作用的真
实的（也是未知的）碰撞作用的具体细节并不重要，
他们提出了一个简化的碰撞模式，将玻尔兹曼方程
简化为：
5f
5t
＋珒v· 5f
5x
＝－
1
τ
（ f－ feq） （12）
在这种形式的玻尔兹曼方程中，碰撞表达式中包含
了在特征时间尺度τ内达到局地平衡分布τ这一松
弛项，该模式被称为 BGK 模式。它描述了分子相互
作用的物理本质，τ为分子碰撞时间，又可称为松弛
时间。尽管 BGK 模式似乎只适用于局地近平衡态，
但人们认识到只要松弛时间的长短能包含相关的物
理特点，这样一种近似将超越其理论局限得到应用。
在湍流的 BGK 模式中，τ可用一个典型的湍流松弛
时间τ turb 代替。 Higuera 等人根据 BGK 模式将格子
玻尔兹曼方程写为如下形式：
fi （ x ＋ e i，t ＋1）－ fi （x，t） ＝ 1
τ turb
［feq
i（x，t）－ fi （x，t）］
（13）
其中，局部平衡解feq
i 是密度Σfi 和动量Σfivi 的给
定函数。可以选择函数feq
i ，以便产生相应的物理现
象，更一般意义上feq
i 可以包括其他物理性质，例如
局部温度等。此方程被称为 LBGK 方程，也有人将
格子玻尔兹曼方法称作 LBGK 模式。
3 格子玻尔兹曼方法的应用
格子玻尔兹曼方法是基于元胞自动机的并行思
想得以实现的，一方面具有元胞自动机的由简单演
化规则显现复杂动力特性的特点，在研究动力系统
时得到丰富的物理图象，从而带来对物理问题深刻
的认识；另一方面由于以玻尔兹曼方程为基础，使其
在模拟流动现象时比元胞自动机更符合实际情况。
格子玻尔兹曼方法被广泛用来模拟各种流体流
动，并被确信是胜过计算流体动力学传统数值技术
的极好方法，其微观水平描述给出数值方法的自然
解释，使之推广到复杂的流动问题，诸如二相流、磁
流体动力学、孔隙介质中的渗流或热流体动力学，以
及湍流运动的数值模拟等等。
3.1 格子玻尔兹曼方法在各种物理领域的应用
元胞自动机最初是被用来模拟非生物过程的自
繁殖系统［11，1
2］，由元胞自动机规则显现出非常复杂
且意想不到的行为。而这之后，又有许许多多的规
则被提出来［1］，如模拟表面生长过程的规则，模拟
森林火灾的概率元胞自动机规则，以及模拟像砂粒
一样的颗粒的基本堆积和倒塌现象的砂堆规则。S.
W olfram 更是将规则按动力学分类，进行了系统的
研究［13－14 ］。而这些规则也从单纯的理论研究到应
用于实际问题中，如 Vichniac 在20 世纪80 年代提
出的相变统计模型（ Ising 模型）的自旋动力 Q2 R 规
则［15，1
6］。
元胞自动机的应用范围越来越广泛，比如扩散
现象、材料晶体凝聚时的生长过程、岩石等非均质材
料的弹、脆、塑性等性质的解释，城市土地利用的演
化、道路交通流的模拟、对地震的模拟、林火的蔓延、
流行病毒的传播、甚至舆论传播等等方面。
扩散现象在物理学、化学和生物学等许多领域
中起着重要作用，其中有些过程非常复杂，难以用扩
散方程的解来描述，而元胞自动机方法非常自然、很
好地适合于进行数值模拟。如产生两种介质接触间
的粗糙界面，对扩散前锋的性质进行研究。生长也
是一种扩散现象，一般是在远离平衡态的时候形成
的分形结构，比如由于曲率原因驱动扩散进行进而
引起晶粒长大的现象、过冷介质中的树枝状凝固现
象等（图2）。现在，这一方法在研究金属材料的再
结晶［17］、金属凝固过程的枝晶生长［18］等方面得到
了广泛的应用。
由此，我们也可以看到，元胞自动机的优势在于
能够通过规则对各种复杂的传播过程进行模拟，除
了前面介绍的扩散、生长现象之外，还有象岩石的破
裂［19］（图3）、地震波传播过程的模拟［20，2
1］。
32 地球科学进展 第22 卷
图2 二维元胞自动机模拟枝晶的生长（黑色）［1］
灰色点代表尚未凝聚的扩散粒子
Fig. 2 Num erical sim ulation of dendrite deve lopm ent
（ black ） by the 2 D cellular autom aton ［1］
In figure ， grays are un- agglom eration grains
除此之外，利用元胞自动机具有模拟复杂系统
的优势，人们还进行了各种社会现象的研究，如通过
分析在人群集中的地方人员的自适应流动，找到出
现瓶径现象的规律，从而为指导城市交通建设提供
理论依据［22］（图4）；对舆论传播这一社会复杂系统
的演化过程进行模拟，利用不同元胞给定值代表对
舆论持不同态度的人，通过一定规则演化之后，得到
给定值接近的元胞聚集在一起，与人们常说的“物
以类聚，人以群分”的社会现象完全吻合［23］（图5）。
在另一方面，元胞自动机的元胞空间与基于几
何定位的地理空间概念有类似的意义，用元胞状态
的改变就可反映出地理特征的改变情况，在土地状
况日益受到关注的今天，利用元胞自动机研究土地
覆盖的变化［24］、人地关系［25 ～27 ］、城市发展规划［28］
等问题取得了众多有价值的成果。
图3 用模型获得的岩石裂隙的不同形式
Fig.3 Num erical sim ulation of the different form s of the rock cranny
图4 人员流动过程中的成拱现象［22］
Fig.4 The arching phenom enon in crow d
pedestrian flow ［22］
3.2 格子玻尔兹曼方法在流体计算中的应用
20 世纪70 年代，基于元胞自动机思想提出的
HPP 模型是能够进行流体运动计算的第一个格子
气模型，但实际上它最初是用于研究气体交互作用
粒子的基本统计性质的。而随着元胞自动机理论的
发展，到80 年代末，人们已预计全离散计算模型将
替代风洞试验，使人们看到这一前景的格子气模型
图5 模拟得到“物以类聚，人以群分”［23］
Fig.5 Sim ulation result on things of one
kind com e together ［2 3］
则是1986 年由 U . Frisch 等［29］以及 S. W olfram ［30］
几乎在同时提出的著名的 FHP 模型。此后，利用格
子气模型进行了一系列的流动现象的数值模拟［31］，
比如，卡门涡阶和有台阶的回流区计算、自由边界问
题、格子气的 MHD （磁流体动力）模型、表面张力问
题、多孔介质问题、悬浮体模型、传热问题以及湍流
问题等。进一步，为了消除格子气模型布尔量求统
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 33
计平均时产生噪音，以及调节模拟的物理系统的参
数欠灵活等问题，研究人员引入了 Maxwell- Boltzm
ann 速度分布函数代替格子气中的布尔量，这种格
子玻尔兹曼方法能够避免数值的不稳定性，并可以
模拟高 Reynolds 数流体的流动，目前被广泛用来模
拟各种流动问题［32，3
3］，如多相流［34 ～36 ］、孔隙介质中
的渗流［37 ～40 ］，模拟风传输雪的过程等（图6），格子
玻尔兹曼方法的微观水平描述在计算这种复杂流动
现象方面具有传统计算方法不可比拟的优势。
除此之外，利用格子玻尔兹曼方法计算的流动
图6 沟道中雪沉积的三个连续阶段［1］（风从左向右吹）
Fig.6 Three continuous snow aggradation m om ents in street canyon ［1］（ W ind blow s from left ）
图7 不同排列的双圆柱绕流涡线图［41］
Fig.7 Vorticity contours in the w ake flow fie lds of various array cylinders ［41］
图8 建筑物内外的流线［44］
Fig.8 Stream lines around the building ［44］
图9 Pow er Flow 计算汽车外部流线①
Fig.9 Stream lines around the car by Pow er Flo w *
问题也越来越广泛，比如，绕流运动［41 ～43 ］、室内流
动［44］（图8）、可压缩的高速流动［45］、浅水涌浪在固
体表面的反射和绕射现象［46］等。
近年来，格子玻尔兹曼方法在数值计算领域越
来越受到重视，并被用于研发形成商业 CFD 软件，
Power FLOW ′就是由美国 EXA 公司开发的基于格子
玻尔兹曼方法的 CFD 软件，在短短不到20 年的时
间里，实现了人们的期望：用全离散计算模型将替代
风洞试验（图9）。
除了这些工程上的应用之外，在大气领域的计
算也初步取得了成果，比如利用旋转坐标系计算
Ros*****y 涡［47］、计算正压大气运动［48，4
9］以及模拟重
力波［50］等。
3.3 格子玻尔兹曼方法研究湍流
湍流也是一种流动现象，即流体的紊乱流动，在
各种不同的流动过程中都能见到，如汽车或飞机周
围的空气流动，大气运动等，我们之所以将其作为一
个单独的题目介绍，在于它的重要性和困难性。
所有的湍流运动可用一组非线性偏微分方程来
描述，这些方程大约200 年前由 Navier 和 Stokes 首
先提出来。 Navier- Stokes 公式的数学求解非常复
杂，一个多世纪来，人们无法得到方程的解，特别是
湍流运动的解［51 ～53 ］。
3.3.1 格子玻尔兹曼方法封闭NS 方程［54］
当湍流发生时，随着能量的级串，能量由外界传
34 地球科学进展 第22 卷
① http ：∥ www . exa. com /
给大尺度涡L ，然后向小尺度（大约η≈Re－3/4 L 量
级，η称为 Kolm ogorov 长度）涡传递，并最终耗散掉。
采用直接数值模拟（ DNS ）求解湍流问题时，必须对
所有被激发的自由度同时求解。但对于多数实际情
况，这个计算量相当巨大，如计算经过汽车的流动，
雷诺数大约是106 ，至少需要计算1014 个自由度，远
远超出现有的最快的计算机的计算能力；因此，对于
高雷诺数流动必须建立模式。
Boussinesq 提出可以将小尺度场u′对大尺度场
U 的平均作用以类似分子撞击大尺度流场的方式来
描述；也就是说，湍流的小尺度是通过涡粘性系数
ve 对大尺度作用的。按照这种思路，湍流流动的大
尺度常用 Navier- Stokes 方程描述，不过这里v 则用
修正的有效粘性系数vt ＝v ＋ve 所代替。
1941 年， Kolm ogorov ［55］提出当雷诺数很大时，
湍流的统计特性具有普适性，也就是说，它们只依赖
湍能的耗散率。因此，在湍流运动中，湍能从大尺度
（外部作用如加热，可认为是原动力）向非常小的尺
度传递，普适性意味着存在一个相当宽的尺度范围，
其中能量传递的统计特性不依赖（或稍微依赖）于
大尺度流动特性。因此，对于Δn L，人们可以将小
于Δ的尺度的动力学特性参数化，得到有效粘性系
数vt（Δ），就可以正确反映能量向小尺度的传递；这
简化了计算，对计算机的计算能力要求也降低了。
然而，对涡粘性系数的基本假设及修正，依赖于
已存在的尺度分离，即可解的大尺度的运动与那些
不可解的脉动涡的分离。但是，由于缺乏尺度间的
明确的划分，使得涡粘性模型的有效性受到质
疑［56，57］。
格子玻尔兹曼方法为解决湍流问题提供了一个
不同的途径，在此方法中，尽管大的湍流涡仍被看作
存在于涡的海洋中，但却不再依赖于尺度分离的假
设。如前所述， Navier- Stokes 方程可通过 LBGK 方
程以幂Kn 展开的低阶截断得到，当Knn 1 ，在两个
基本尺度间存在尺度分离，格子玻尔兹曼方程（式
13）包含了涡粘性模式，湍流涡粘性系数正比于
τ turb ；当Kn 不很小且不存在尺度分离时，τ turb 表示各
尺度湍流涡粘性的重整化求和，本质上优于 Navier-
Stokes 方程中依赖于尺度分离的涡粘性系数。下
面，用此方程导出更普适的涡粘性模式，即确定湍流
脉动的自松弛时间τ turb
［54］。
τ turb 应依赖于反映不同湍流机制的时间尺度，
τ turb ＝τ0 ＋Ψ（k /ε，S－1，⋯） （ 14）
其中，τ0 是基本的分子松弛时间，k 是湍流动能，ε
是湍流耗散率，S 是当地速度梯度的度量。这里，k /
ε和S－1 是湍流特征时间尺度。由于存在其他的湍
流作用，比如浮力和旋转，方程14 中可能还要包括
其他的固有时间。因此，方程中Ψ的典型的函数形
式应当是湍流时间尺度的谐波平均（也就是说是湍
流频率的叠加）。通过推导得到有效松弛时间τ turb ，
τ turb ＝τ0
＋C μ
k2 ε
T（1 ＋η2 ）1 /2 （15）
其中，C μ＝0.085 ，η＝ Sk / ε，T 是绝对温度。在
Kn→0 时，即两个基本尺度间存在尺度分离，可导出
基于重整化群的k－ ε涡粘性模型。
应用有效松弛时间τ turb 进行计算时，不必调整
模式参数。图9 是应用此方法计算流动经过一个实
际的小轿车。
3.3.2 格子玻尔兹曼方法直接数值模拟
在上一节，我们介绍了 Kolm ogorov ［58］提出当雷
诺数很大时，能量由外界输入给大尺度涡，然后向小
尺度涡传递，并耗散掉。由此可见，涡旋是流体运动
的肌腱，它的生成、发展和演化，以及它与外部流动
和物体之间的相互作用，影响了整个流场的特性。
涡旋研究的对象可以来自于解析解，如 Ossen 涡、
Taylor 涡、 Burgers 涡等，这些具有解析解的涡结构都
比较简单，对这些简单涡旋结构的研究是分析复杂
涡旋流动的基础［59］。而在20 世纪50 ～60 年代，人
们在湍流场中又发现了拟序结构（ Coherent structure
）［60］。拟序结构与涡旋有着极为密切的关系，一
般认为拟序结构是存在于湍流场中的具有较大横向
尺度的有规律的涡旋结构，但严格的说，只要是有规
律的涡旋结构，不论尺度多大，都可以看成是拟序结
构。湍流就是由拟序结构和随机小涡背景流场共同
构成的。
面对湍流这一复杂的现象，人们最早采用的是
实验的方法，如雷诺实验，但这种流动显示实验很难
进行定量的分析。之后的热线风速仪以及超声风速
仪，可以得到流场中单点速度的定量结果，但对流场
有一定影响。近年来激光多普勒测速技术可以对流
场进行非接触性的测量，但成本较高。20 世纪80
年代初，一种名为 PIV （图像粒子测速）的流场测量
技术逐渐发展起来，它的原理是通过向流场中投放
示踪粒子，然后测量和流体混合均匀的粒子的运动
速度来间接测量粒子周围流体的速度。 PIV 技术的
特点是它不再是单点速度的测量，它可以测量一个
平面内二维的流场，使人们可以得到流动的更多的
信息。
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 35
除了这些实验手段，随着计算机存储容量和计
算速度的飞速提高，数值模拟作为一种研究手段也
相应发展起来，利用数值模拟的方法，人们可以得到
整个三维流动速度场的时间序列，而且其成本比实
验要低很多，人们已经利用数值模拟的方法对湍流
流动的各种统计特性进行了许多的研究［61，62］。但
以往的直接数值模拟都是基于求解NS 方程，常用
的方法有差分法、有限元法和谱方法等，由前面的分
析我们知道，其计算量是相当巨大的。而格子玻尔
兹曼方法由于其基于分子运动方程求解，程序设计
非常简单，还有一个优点是易于并行化，是实现这样
庞大的计算问题的一个重要保证。
用格子玻尔兹曼方法进行直接数值模拟，最开
始也是计算周期边界条件的湍流演化过程，并与传
统的谱方法进行比较，比较的结果两者在计算精度
上并无太大差别［32，3
3，6
3，6
4］。之后，还用这种方法计
算槽道流，研究其近壁面的特性，以及直接数值模拟
Rayleigh- Taylor 流［65］、两种流体界面上 Kelvin －
Helm holtz 波［66］等图10。最近，还出现了一种熵格
子玻尔兹曼方法［67］，模拟等温条件下的运动，其中
可调节直接数值模拟的截断尺度，只计算大尺度的
运动。
图10 直接数值模拟 Kelvin - Helm holtz 波［6 6］
Fig.10 Direct num erical sim ulation of Kelv i n－ Helm holtz w av e［66］
3.3.3 格子玻尔兹曼方法研究大气湍流
大气湍流通常比实验室里的湍流要更加复杂，
在大气运动中，除了由于摩擦产生的机械湍流外，还
有由于浮力产生的对流热泡造成的湍流运动，由于
地面吸收太阳辐射后升温导致热气上升，形成对流，
产生湍流运动，此时的湍流运动是硬湍流运动，人们
发现在硬湍流条件下，温度场级串有别于速度
场［68］。
湍流能量级串模式一直是湍流研究中的前沿问
题，对此问题的认识将最终影响湍流模式的参数化。
目前，已经出现了各种级串模型，但对于能量级串这
一问题的研究还有待于进行多方面深入的研究。
Kolm ogorov 理论认为从大涡向小涡的逐级级串
具有自相似结构（由 Richardson 首先提出），它们充
满整个空间，该理论（又称 K41 理论）给出的结构函
数及能量耗散率的标度律在高阶矩时与实验结果明
显不符。 Landau ［69］首先向 K41 理论提出了质疑，指
出将能量耗散律ε看作常数缺乏物理根据，应当考
虑ε的涨落；不久之后，实验发现了湍流的间歇性，
湍流研究者才认识到间歇性对建立速度结构函数标
度律的重要意义，并出现了一系列计入间歇性的一
些重要模型［70 ～73 ］，如β－模型，分形级串模型， GOY
壳模型，对数－Poisson 模型和 S- L 模型等，理论预
测和实验结果在较高阶矩时相当一致。这是由于
时—空间歇性的存在意味着湍能从大涡（大尺度）
向小涡（小尺度）的逐级传递过程，既不充满空间，
也不能看作是平稳随机过程。
最近，一种新的级串过程———同步级串模式被
提出来［74］图11，该模式认为从外界向含能区注入
能量到形成发达湍流其实并不需要经过n 次有序的
图11 湍流能量级串模式
Fig.11 Turbulent energy cascade m odel
36 地球科学进展 第22 卷
级串过程，正如一块玻璃被摔碎时，最可能的情形是
破碎成各种不同尺寸的碎块，含能涡一次散裂同时
形成大、中、小、微各种尺度的湍涡的概率也应远大
于 Kolm ogorov 级串模式的概率。
综上所述，经过湍流研究工作者几十年的积累，
关于湍流速度场的统计特征已经了解得较为深入，
但是对作为被动保守物的标量场（例如温度场）了
解得还很不够［75 ～77 ］。人们过去一度认为标量场应
具有与速度场相同概率分布，但是后来从观测事实
中发现二者之间有很多的不同，标量场的概率分布
形态更显复杂。20 世纪80 年代美国芝加哥大学在
封闭容器中利用氦气做热对流实验（文献中称为
“芝加哥实验”）［78 ～80 ］，结果表明湍流温度脉动场的
概率分布类型与外参数 Rayleigh 数（Ra）有密切关
系。存在一临界Ra 值（≈108），即当 Ra ＜10 8 时，湍
流温度场中两点温度差的概率密度分布非常接近于
Gauss 分布（软湍流）；但当 Ra ＞10 8 时，概率密度分
布变成非 Gauss 分布，而更接近于指数分布（硬湍
流）。Ra 是一个反映浮力作用与热扩散作用相对的
无量纲数，定义为：
Ra ＝
αgδL3
κν
（16）
式中α表示热膨胀系数，g 是重力加速度，δ即驱动
对流的外部温差，L 代表直径或高，v 和κ分别为分
子粘性系数和热扩散系数。到目前为止，还有其它
很多实验室模拟和数值模拟工作也进一步证实了芝
加哥实验的结果［76，7
7］。
由此可见，对于包含热运动的大气湍流，能量级
串过程将会更加复杂，而通过格子玻尔兹曼方法直
接数值模拟这种湍流运动，进行机制的研究应是非
常有意义的。利用格子玻尔兹曼方法对温度场的研
究已有很多［81，8
4］，但都还限于平均场的模拟，可以
先用格子玻尔兹曼方法直接数值模拟含有热对流运
动的 Rayleigh- Bernard 流，了解热对流的湍流运动，
然后模拟更加真实的大气湍流运动。
4 结 语
通过前面的介绍，我们对格子玻尔兹曼方法的
产生和应用有了比较全面的认识：格子玻尔兹曼方
法是基于20 世纪60 年代的并行思想及由此产生的
元胞自动机，同时建立在分子运动论思想上的玻尔
兹曼方程为其物理基础。除了用于模拟各种复杂物
理、社会现象之外，在对流动现象的模拟方面，格子
玻尔兹曼方法一方面由于其分子运动论的物理背
景，使得其在模拟流体流动时不会出现元胞自动机
那样由于布尔运算而产生的噪音，所以广泛应用于
流动现象的数值模拟；另一方面由于其微观水平描
述给出数值方法的自然解释，对于模拟流动中的复
杂现象，如湍流问题要优越于NS 方程，目前在这方
面的探索已取得了许多成果，使得对湍流这一古老
的问题的研究又增添了一个有力的工具；并尝试分
析了利用这一方法研究大气湍流的初步方案。
参考文献（ Referem ces ）：
［1］ Chopard B ， Droz M . Cellular Autom ata Modeling of Physica l System
s ［M］. London ： Cam bridge University Press ，1998.［祝玉学，
赵学龙译. 物理系统的元胞自动机模拟［M］. 北京：清华大学
出版社，2003.］
［2］ Hardy J， Pom eau Y ， Pazzis O De. Molecular Dynam ics of a classical
lattice gas ： Transportproperties and tim e correlation fu nctions
［J］. Phys. Rev. A ，1976 ，13：1 949-1 960.
［3］ Mc Nam ara G G ， Zanetti G . Use ofthe Boltzm ann Equation to Sim -
ulate Lattice Gas Autom ata ［J］. Physical Review Letters ，1988 ，
61：2 332-2 335.
［4］ Higuera F ，J im enez J， Succi S. Boltzm ann approach to lattice gas
sim ulations ［J］. Europhysics Letters ，1989 ，9：663 .
［5］ Feng Duan ， Feng Shaotong. The W orld of Entropy ［M］. Bei jing ：
Science Press ，2005.［冯端，冯少彤. 溯源探幽———熵的世界
［M］. 北京：科学出版社，2005.］
［6］ Harris ， Stewart. An Introduction to the theory ofthe B oltzm ann E -
quation ［M］. New York ： Holt ， Rinehart and W inston ，1971.
［7］ Li Yuanxiang ， Kang Lishan ， Chen Yuping. Lattice Gas Autom aton
［M］. Bei jing ： Tsinghua University Press ，1994.［李元香，康
立山，陈毓屏. 格子气自动机［M ］. 北京：清华大学出版社，
1994.］
［8］ Higuera F ，J im enez J， Succi S. Lattice gas dynam ics with enhanced
collision ［J］. Europhysics Letters ，1989 ，9：345.
［9］ Quian Y H ， Hum ieres D ， Lallem and P. Lattice BGK Models for
Navier- Stokes Equation ［J］. Europhysics Letters ，1992 ，17（6）：
470-484.
［10］ Bhatnager P ， Gross E P ， Krook M K . A m odel for collision
process in gases ［J］. Physics Review ，1954 ，94：511.
［11］ Pesavento U . An im plem entation of von Neum an n ＇ s self- reproducing
m achine ［J］. Artificial Life ，1995 ，2：337-354.
［12］ Yan Guangwu. Advances of the Studies on the Ce llular Autom ata
and Artificial Life ［J］. J ournalofJ ilin University ，2003 ，41（1）：
40-44.［闫广武. 元胞自动机与人工生命研究进展［J］. 吉林
大学学报，2003 ，41（1）：40-44.］
［13］ W olfram S. Statisticalm echanics ofcellula r autom ata ［J］. Review
of Modern Physics ，1983 ，55（3）：601-644.
［14］ W olfram S. Com putation theory of cellular au tom ata ［J］. Com -
m unication in Mathem aticat Physics ，1984 ，95：15-57.
［15］ Vichniac G . Sim ulating physics with cellula r autom ata ［J］. Physica
D ，1984 ，10：96-115 .
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 37
［16］ Yu Lu ， Hao Bolin ， Chen Xiaosong. Phase Transitions and Critical
Phenom ena ［M ］. Bei jing ： Science Press ，2005.［于禄，郝柏
林，陈晓松. 边缘奇迹———相变和临界现象［M］. 北京：科学
出版社，2005.］
［17］ Zheng Chengw u ， Lan Yong jun ， Xiao Nam in ， et al . Mesoscale
sim ulation of Static recrystallization ofho tdeform ed austenite in a
low carbon steel ［J］. Acta Metallurgica Sinica ，2006 ，42（5）：
474-480.［郑成武，兰勇军，肖纳敏，等.热变形低碳钢中奥氏
体静态再结晶介观尺度模拟［J］. 金属学报，2006 ，42（5）：
474-480.］
［18］ Guo Dayong ， Yang Yuansheng. Num erical Sim ulation of Dend rite
Developm ent in Ternary Alloy Solidification ［J］. Foundry ，2006 ，
55（6）：601-603 .［郭大勇，杨院生.三元合金凝固过程枝晶生
长数值模拟［J］.铸造，2006 ，55（6）：601-603.］
［19］ W ang Shim in ， Zhu Hehua ， Feng Xiating ， etal . Influence ofheterogeneity
on m acroscopical crack form of the brittle rock ［J］.
Rock and Soil Mechanics ，2006 ，27（2）：224-227.［王士民，朱
合华，冯夏庭，等.细观非均匀性对脆性岩石材料宏观破坏形
式的影响［J］.岩土力学，2006 ，27（2）：224-227.］
［20］ Dong Qiaoliang ， Zhou Yubing ， Yao Yao. Seism ic Sim ulating by
Using Lattice Boltzm ann Method ［J］. J ournal of China University
of Geosciences （ Earth Science ），1997 ，22（6）：638-642. ［董桥
梁，周玉冰，姚姚.格子 Boltzm ann 方法地震波场模拟［J］.中
国地质大学学报：地球科学，1997 ，22（6）：638-642.］
［21］ W ang Zhenli ， Li Youm ing. Parallel Algorithm for Sim ulatin g
Seism ic W ave Propagation by Cellular Autom at a［J］. Chinese
J ournal of Geophysics ，1999 ，42（3）：410 -415.［王真理，李幼
铭，细胞自动机地震波模拟的并行化算法［J］.地球物理学
报，1999 ，42（3）：410-415.］
［22］ Zhang Peihong ， Huang Xiaoyan ， W an Huanhuan ， et al . Study
on Adaptive Cellular Autom ation Model on Crow d Pedestrian
Flow ［ J］. J ournal of Shenyang J ianzhu University ，2006 ，22
（2）：289-293 .［张培红，黄晓燕，万欢欢，等.人员群集流动
自适应元胞自动机模型研究［J］.沈阳建筑大学学报，2006 ，
22（2）：289-293.］
［23］ Liu Muren ， Deng Minyi ， Kong Ling jiang. The Cellular Autom ata
Model for Public Opinion Evolution ［J］. J ournalof Guangxi Norm
al University ，2002 ，20（2）：1.［刘慕仁，邓敏艺，孔令江.舆
论传播的元胞自动机模型［J］. 广西师范大学学报，2002 ，20
（2）：1.］
［24］ Yan J infeng ， Chen Xi. Analysis and Stim ulation Method Disc ussion
on LUCC of Arid Region Based on GIS ［J］. Arid Land Geography
，2003 ，26（2）：185 -191.［阎金凤、陈曦.基于 GIS 的干
旱区 LUCC 分析和模拟方法探讨［J］.干旱区地理，2003 ，26
（2）：185-191.］
［25］ Luo Ping ， Du Qingyuan ， Lei Yuanxin ， et al . Cellular Autom ata
Based on Geographic Feature and Urban Land Use Evolvem ent
［J］. Geom atics and Inform ation Science of Wuhan Un iversity ，
2004 ，29（6）：504 -507.［罗平、杜清远、雷元新，等.地理特征
元胞自动机及城市土地利用演化研究［J］.武汉大学学报，
2004 ，29（6）：504-507 .］
［26］ Hou Xiyong ， Chang Bin ， Yu Xinfang. Land use change in Hexi
corridor based on A - Markov m ethods ［ J］. Transactions of the
CSAE ，2004 ，20（5）：286 -291.［侯西勇、常斌、于信芳. 基于
CA - Markov 的河西走廊土地利用变化研究［J］. 农业工程学
报，2004 ，20（5）：286-291.］
［27］ Chen J ianping ， Ding Huoping ， W ang Gongwen ， et al . Desertification
Evolution Modeling through the Integr ation of GISand Cellular
Autom ata ［J］. J ournal of Rem ote Sensing ，2004 ，6（3）：
254-260.［陈建平，丁火平，王功文，等.基于 GIS 和元胞自动
机的荒漠化演化预测模型［J］. 遥感学报，2004 ，6（3）：254-
260.］
［28］ Cui W eihong ， Zhang Xianfeng. Study on Dynam ic Monitoring
and Dynam ic Sim ulating of Land Resource ［J］. Geo- Inform ation
Science ，2002 ，1：79-85.［崔伟宏，张显峰.土地资源的动态
监测和动态模拟研究［J］.地球信息科学，2002 ，1：79-85.］
［29］ Frisch U ， Hasslacher B ， Pom eau Y . Lattice gas autom ata for the
Navier- Stokes equation ［J］. Phys. Rev. Lett ，1986 ，56：1505.
［30］ W olfram S. Cellular autom ation fluid ： Basic theory ［J］.J. Stat.
Phys ，1986 ，45：471 .
［31］ Qian Yuehong ， D d ＇ Hum ieres ， Y Pom eau ， et al . The newest developm
ent of Lattice Gas Autom aton Dynam ics ［J］. Mechanics in
Engineering ，1990 ，12（1）：7-16.［钱跃宏，D d’ Hum ieres ，Y .
Pom eau ，等.格子气流体动力学及其最新进展［J］.力学与实
践，1990 ，12（1）：7-16.］
［32］ Benzi R ， Succi S ， Vergassola M . The lattice Boltzm ann equation
： theory and applications ［J］. Physics Reports ，1992 ，222
（3）：145-197.
［33］ Shiyi Chen ， Gary D Doolen. Lattice Boltzm ann Method for Fl uid
Flow s ［J］. Annu. Rev. Fluid Mech ，1998 ，30：329-364.
［34］ Shan Xiaowen ， Chen Hudong. Lattice Boltzm ann m odelfor sim u -
lating flows with m ultiple phases and com pone nts ［J］. Physical
Review E ，1993 ，47（3）：1 815-1 819.
［35］ Martys N S ， Chen Hudong. Sim ulation of m ulticom ponent fl uids
in com plex three- dim ensional geom etries by t he lattice Boltzm ann
m ethod ［J］. Physical Review E ，1996 ，53（1）：743-750 .
［36］ Zhang Raoyang ， Chen Hudong. Lattice Boltzm ann m ethod for
sim ulations ofliquid- vapor therm alflow s ［J］. Physical Review E ，
2003 ，67：066711 .
［37］ Zhan Hongren ， Li Dongyan ， W ang Guoheng. Lattice gas autom ata
m odeling on evaporation process ofporous m edia ［J］. Chem ical
Equipm ent Technology ，2005 ，26（6）：1-5.［战洪仁，李冬艳，王
国恒.多孔介质干燥过程的格子气自动机法模拟［J］.化工装
备技术，2005 ，26（6）：1-5.］
［38］ Zheng Zhong ， Gao Xiaoqiang ， Shi W anyuan. Lattice Gas Autom -
ata Modeling of Fluid Flow Through Porous Medi a［J］. J ournalof
Chem ical Industry and Engineering ，2001 ，52（5）：406-409.
［郑忠，高小强，石万元. 多孔介质中流体流动的格子气自动
机模拟［J］.化工学报，2001 ，52（5）：406-409.］
［39］ Zheng Zhong ， Gao Xiaoqiang. Lattice Gas Autom ata Method fo r
Modeling Fluid Flow and Heat Transfer in Metal lurgical Porous
Media ［J］. Acta Metallurgica Sinica ，2000 ，36（4）：433-437.
［郑忠，高小强.格子气自动机法模拟冶金多孔介质传输现象
［J］.金属学报，2000 ，36（4）：433-437.］
38 地球科学进展 第22 卷
［40］ Yue W enzheng ， Tao Guo ， Zhu Keqin. Sim ulation of electrical
transport properties in oil- water saturated porous m edia with 2 D
lattice gas autom ata ［J］. Chinese J ournal of Geophysics ，2005 ，48
（1）：189-195.［岳文正，陶果，朱克勤.二维格子气自动机模
拟孔隙介质的电传输特性［J］.地球物理学报，2005 ，48（1）：
189-195.］
［41］ Surm as R ， Santos L ， Philippi P. Lattice Boltzm ann sim ulation of
the flow interference in bluff body w akes ［J］. Future Generation
Com puter System s ，2004 ，20 ：951-958.
［42］ Takada ， Naoki ， Tsutahara ， Michihisa. Evolution of viscous flow
around a suddenly rotating circular cylinder in the Lattice Boltzm
ann Method ［J］. Com puters ＆ Fluids ，1998 ，27（7）：807-828.
［43］ Feng Shide ， Zhao Ying ， Gao Xianlin ， et al . The sim ulation of
2 D flow around a cylinder by lattice BGK m odel ［J］. Acta Mechanica
Sinica ，2002 ，34：254-258.［冯士德，赵颖，郜宪林，等.
用格子 BGK 模型模拟两维圆柱绕流［J］. 力学学报，2002 ，
34：254-258.］
［44］ Crouse B ， Krafczyk M ， Kuhner S ， et al . Indoor air flow analysis
based on lattice Boltzm ann m ethods ［J］. Energy and Buildings ，
2002 ，34：941-949.
［45］ Yu Huidan ， Zhao Kaihua. Lattice Boltzm ann m ethod for com -
pressible flow s w ith high Mach num bers ［J］. Physical Review E ，
2000 ，61（4）：3 867 -3 870 .
［46］ Liu Feng ， Hu Fei. Application of lattice Boltzm ann m eth od in
sim ulating reflection and diffraction ofsur ge waves ［J］.J ournalof
Hydraulic Engineering ，2004 ，3：22-30.［刘蜂，胡非.用格子
Boltzm ann 方法模拟涌波的反射与绕射［J］.水利学报，2004 ，
3：22-30.］
［47］ Yu huidan ， Zhao Kaihua. Ros*****y vortex sim ulation on a par aboloidal
coordinate system using the lattice Bo ltzm ann m ethod ［J］.
Physical Review E ，2001 ，64：056703.
［48］ Liu Feng ， Hu Fei. A prelim inary study on the constructio n of atm
ospheric dynam ic m odels w ith lattice Boltzm ann m ethod ［J］. Acta
Meteorologica Sinica ，2003 ，61（3）：267-274.［刘峰，胡非.
用格子 Boltzm ann 方法构建大气动力学模式的初步研究［J］.
气象学报，2003 ，61（3）：267-274.］
［49］ Liu Feng ， Hu Fei. Lattice Boltzm ann m ethod to sim ulate b arotropic
atm osphere ［J］. J ournalofthe Graduate Schoolofthe Chinese
Academ y of Sciences ，2003 ，20（2）：168-171.［刘峰，胡非，
用格子 Boltzm ann 方法模拟正压大气运动［J］.中国科学院研
究生院报，2003 ，20（2）：168-171.］
［50］ Buick J， Greated C . Lattice Boltzm ann m odeling of inte rfacial
gravity waves ［J］. Physics of Fluids ，1998 ，10（6）：1 490 -1 511.
［51］ Frisch U . Turbulence- The Legacy of a. Kolm og orov ［M ］. Cam -
bridge London ： Cam bridge University Press ，1996.
［52］ Frisch U ， Orszag S. Turbulence- challenges for theory a nd experim
ents ［J］. Physics Today ，1990 ，43：24-32.
［53］ Sreenivasan K R . Fluid Turbulence ［J］. Review of Modern Physics
，1999 ，71：383-395.
［54］ Chen Hudong ， Satheesh Kandasam y ， Steven Orszag ， et al . Extended
Boltzm ann kinetic equation for turbul ent flow s ［J］. Science
，2003 ，301 ：633-636 .
［55］ Kolm ogorov A N . Localstructure ofturbulenc e in an incom pressible
viscous fluid at very high Reynolds num ber s［J］. Dokl Akad
Nauk SSSR ，1941 ，30：301-305 （英文重印： Proceeding of the
Royal Society of London A ，1991 ，434 ：9-13.）
［56］ Yakhot V ， Orszag S. Renorm alization group analysis of t urbulence
I. Basic theory ［J］. J ournal of Science Com put ，1986 ，1：
3-51.
［57］ Yakhot V ， Orszag S ， Thangam S ， et al . Developm ent of turbulence
m odels for shear flow s by a double expans ion technique
［J］. Physics of Fluids A ，1992 ，4：1 510 -1 520.
［58］ Kuchem ann D . Reporton the IUTAM sym posium on concentrated
vortex m otion in fluids ［J］. J ournal of Fluid Mechanics ，1965 ，
21：1.
［59］ Ma Bin. The identification and visualizatio n of vortex ［D］. Beijing
： Tsinghua University ，1997.［马兵. 涡旋的识别与可视化
［D］. 北京：清华大学工程力学系，1997.］
［60］ Kline S J， Reynolds W C ， Schraub F A ， et al . The structure of
turbulent boundary layers ［J］. J ournalof Fluid Mechanics ，1967 ，
30：741.
［61］ Chen J ian. Direct Num erical Sim ulation of Sp atial Sm ooth Pipe
Flow Transition ［D］. Bei jing ： Tsinghua University ，2002.［陈
健. 空间发展光滑圆管转捩的直接数值模拟［D］. 北京：清
华大学工程力学系，2002.］
［62］ W u Minw ei. Num erical Research on Turbulent P ipe Flow ［D］.
Bei jing ： Tsinghua University ，2002.［伍敏伟.圆管湍流的数值
研究［D］. 北京：清华大学工程力学系，2002.］
［63］ Martinez D ， Matthaeus W ， Chen S. Com parison ofspectralm ethod
and lattice Boltzm ann sim ulations oftw o- d im ensional hydrodynam
ics ［J］. Physics Fluids ，1994 ，6（3）：1 285-1 298.
［64］ Yu H ， Girim a ji S ， Luo L. Lattice Boltzm ann sim ulations ofdecaying
hom ogeneous isotropic turbulence ［ J］. Physical Review E ，
2005 ，71：016708 .
［65］ He Xiaoyi ， Chen Shiyi and Zhang Raoyang. A lattice Boltzm ann
schem e for incom pressible m ultiphase flow an d its application in
sim ulation of rayleigh- taylor Instability ［J］. J ournal of Com putational
Physics ，1999 ，152 ：642-663 .
［66］ Zhang Raoyang ， He Xiaoyi ， Doolen Gary ， et al . Surface tension
effects on two- dim ensional tw o- phase Kelvin - Helm holtz instabilities
［J］. Advances in Water Resources ，2001 ，24 ：461-478.
［67］ Boghosian B ， Yepez J， Coveney P ， et al . Entropic lattice Boltzm
ann m ethods ［J］. Proceeding ofthe Royal Society of Lond on A ，
2001 ，457 ：717-766 .
［68］ Hu Fei ， Cheng Xueling ， Zhao Songnian ， etal . Hard state ofthe
urban canopy layer turbulence and its self- si m ilar m ultiplicative
cascade m odels ［J］. Science in China ，（ Series D ），2005 ，35
（ Supp. I ）.［胡非，程雪玲，赵松年，等.城市冠层中温度脉动
的硬湍流特性和相似性级串模型［J］.中国科学D 辑，2005 ，
35（增刊1）.］
［69］ Landau L D ， Lifshitz E M . Fluid Mechanics ［M ］. Oxford ： Pergam
on Press ，1987.
［70］ She ZS ，J ackson E ， Orszag SA . Interm ittentvortex structures in
hom ogeneous isotropic turbulence ［J］. Nature ，1990 ，344 ：226-
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 39
228.
［71］ She Z S ， Leveque E . Universal scaling law s in fully dev eloped
turbulence ［J］. Phys. Rev. Lett ，1994 ，72：336-338.
［72］ Ruelle D ， Takens F . On the nature ofturbulence ［J］. Com m unication
in Mathem atical Physics ，1971 ，20：167 -192.
［73］ Benzi R ， Biferale L ， Ciliberto S ， et al . Generalized scaling in
fully developed turbulence ［J］. Physica D ，1996 ，96：162-181.
［74］ Zhao Songnian. Synchrocascade pattern in th e atm ospheric turbulence
［J］.J ournal of Geophysical Research ，2003 ，108.
［75］ She ZS. On the scaling laws oftherm alturbule ntconvection ［J］.
Physics Fluids A ，1989 ，1（6）：911 -913.
［76］ Ching E S C . Probability densities of turbule nt tem perature fluctuations
［J］. Physics Rewiew Lett ，1993 ，70（3）：283-286.
［77］ Pope SB ， Ching E SC . Stationary probability density fu nctions ：
an exact result ［J］. Physics Fluids A ，1993 ，5（7）：1 529-1 531.
［78］ Heslot F ， Castaing B ， Libchaber A . Transition to turbulence in
helium gas ［J］. Physics Review A ，1987 ，36：5 870-5 873.
［79］ Castaing B ， Gunaratne G ， Heslot F ， et al . Scaling of hard therm
al turbulence in Rayleigh- Bernard convecti on［J］. J ournal of
Fluid Mechanies ，1989 ，204 ：1-30.
［80］ W u X Z ， Kadannoff L ， Libchaber A ， et al . Frequency power
spectrum oftem perature fluctuations in free convection ［J］. Physical
Review Letters ，1990 ，64：2 140-2 143.
［81］ Feng Shide ， Zhao Ying ， Gao Xianlin ， et al . Som e progress in
the lattice Boltzm ann m odel and hydrodynam ic equations ［J］. Acta
Mechanic Sinica ，2002 ，34：254-258 .［冯士德，赵颖，郜宪
林，等. 用格子 Boltzm ann 模型的改进与流体力学方程［J］.
力学学报，2002 ，34：254-258.］
［82］ Zhang Raoyang ， Chen Hudong. Lattice Boltzm ann m ethod for
sim ulations ofliquid- vapor therm alflow s ［J］. Physical Review E ，
2003 ，67：066711 .
［83］ He Xiaoyi ， Chen Shiyi ， Doolen Gary D . A novel therm al m odel
for the lattice Boltzm ann m ethod in incom pres sible lim it ［J］.J ournal
of Com putational Physics ，1998 ，146 ：282-300.
［84］ Shi Yong ， Zhao T S ， Guo Z L. Therm allattice Bhatnagar- Gross-
Krook m odel for flows w ith viscous heat dissip ation in the incom -
pressible lim it ［J］. Physical Review E ，2004 ，70：066310.
The Applic atio n of Lattic e Boltzm ann M eth od i n
th e Atm ospheric Turbule nce Stu dy
CHENG Xue- lin g ， HU Fei ， ZHAO Song- nia n ，J IA NG J in - hua
（ State Key Laboratory of Atm ospheric Boundary Layer Physics and Atm ospheric Chem istry
Institute of Atm ospheric Physics ， Chinese Academ y of Science ， Bei jing 100029 ， China ）
Abstract ： The ob jective of this paper is to system ically introduce the lattice Boltzm ann m ethod from i ts m athem
atic and physical bases ， and investigate broadly its application in ev ery field. The lattice Boltzm ann m ethod
based on the idea ofcellular autom aton w as pro posed by J ohn von Neum ann in the 60 ＇s oflast cen tury. The cellular
autom aton sim ulates the parallelizable char acter of the brain and constructs the dynam ic e volvem ent system . That
starts a new direction to solve the com plex pro blem s. W hen it is used to solve the flow problem ， it com bines w ith
the m olecule kinetic theory of 19 century to so lve Boltzm ann transportequation. Itm eans to study the flow by sim ulating
the m olecules m ovem ent in the flow . So i t is also called lattice Boltzm ann m ethod ， or lattice gas autom aton.
From investigation ， itcan be seen thatthe application fields are v ery broad. For exam ple ， itis used to sim ulate the
crystal agglom eration ， the evolvem ent of city land use ， the traffic flow ， the seism ic w ave ， the spread of fire in the
forest ， the prevalence of virus and the spread of publi c opinion. And by this m ethod ， good results are get. In the
70＇s of last century ， Hardy ， Pom eau and Pazzis builtthe first lattice gas a utom aton to sim ulate flow . Then ， because
it can explain flow from m icrocosm ic level ， itis very suitto solve com plex flow such as m ul tiphase flow ， porous flow
and snow grain transportin the w ind. Moreover ， the com m erce softw are — Pow er Flow is developed based on the lattice
Boltzm ann m ethod. For the old problem — turbulence ， m any m ethods are used and LBM now gradually bec om es
a new w ay to directly sim ulate it because LBM ca n calculate in the m olecule level. The atm osph eric turbulence is
m ore com plex then usual turbulence. Here ， w e especially introduce the atm ospheric turb ulence character and give
an elem entary pro ject to sim ulate by LBM .
Key w ords ： Lattice Boltzm ann m ethod ； Cellular autom ata ； Lattice gas ； Direct num erical sim ulation ； Atm ospheric
turbulence.
40 地球科学进展 第22 卷
第22 卷第3 期
2007 年3 月
地球科学进展
ADVANCES IN EARTH SCIENCE
Vol.22 No.3
Mar.，2007
文章编号：1001-8166 （2007 ）03- -12
格子玻尔兹曼方法及其在大气湍流研究中的应用
*
程雪玲，胡 非，赵松年，姜金华
（中国科学院大气物理研究所大气边界层物理和大气化学国家重点实验室，北京 100029 ）
摘 要：文章的目的是对格子玻尔兹曼方法进行系统的介绍，格子玻尔兹曼方法（ Lattice Boltzm ann
Method ）的出现直接来源于20 世纪60 年代的元胞自动机（ Cellular Autom ata ）思想，而这一方法用
于解决流动现象时，又可以追溯到19 世纪的分子运动论，求解的是 Boltzm ann 提出的玻尔兹曼输运
方程，因此将这一方法称为格子玻尔兹曼方法，之前也被称为格子气自动机（ Lattice Gas Autom aton
）。该方法多用于研究复杂现象，如材料晶体凝聚时的生长过程、城市土地利用的演化等方面。
在20 世纪70 年代由 Hardy 、 Pom eau 和 Pazzis 建立了第一个用于研究流体运动的格子气自动机，此
后，这一方法被广泛用来模拟各种流动问题，诸如二相流、孔隙介质中的渗流等，并根据这一方法开
发了相应的商业软件 Pow er Flow 。同时，格子玻尔兹曼方法由于其在微观水平描述运动的特点，成
为研究湍流的一个很好的数值计算工具，特别是用其进行直接数值模拟（ DNS ）计算，成为继传统的
差分法、有限体积法和谱方法之后的又一有力的手段。而作为大气运动的一个主要现象的大气湍
流，比普通湍流更加复杂，在这里着重介绍了大气湍流的特点和应用格子玻尔兹曼方法模拟湍流的
发展过程。
关 键 词：格子玻尔兹曼；元胞自动机；格子气；直接数值模拟；大气湍流
中图分类号： P425.2 文献标识码：A
1 引 言
元胞自动机是空间和时间都离散、物理参量只
取有限数值集的物理系统的理想化模型。元胞自动
机的思想首先是由J ohn von Neum ann 提出的，尽管
由于他在20 世纪40 年代末参与了第一台数字计算
机的设计，而使他的名字与当代串行计算机联系起
来，但他通过模仿人脑的并行化特点而提出的元胞
自动机思想，第一次人工构造出动力演化系统，为解
决复杂问题开创了新的方向。在20 世纪60 年代他
和 S. Ulam 构造了一个自复制系统，这个系统中元
胞的活动是同时进行的，所有元胞的演化规则是相
同的，随邻近元胞的状态而变化，这一系统用来模拟
生物的自复制过程。尽管这一结果连极原始的生命
形式也算不上，但却具有非常重大的意义。因为通
常认为，机器只能构造比自身简单的结构，而采用自
复制元胞自动机，却可以获得一种能产生新的、具有
同样复杂性和功能的“机器”。继 Von Neum ann 的
研究之后，另一些学者也遵循同样的研究路线，继续
研究这个问题。在1970 年，这方面的工作被 Conw -
ay 等加以发展，提出了著名的“生命游戏”（ Gam e of
Life ）模型，这是一个二维元胞自动机的实例。结果
表明，生命游戏机有出乎意料的丰富行为。而元胞
自动机真正被系统地作为动力系统来研究和应用当
推 S. W olfram ，他于1983 年在《现代物理评论》杂志
上发表了奠基性的论文“元胞自动机的统计力学”，
* 收稿日期：2006-11-27 ；修回日期：2007-01-11.
* 基金项目：国家自然科学基金项目“大气湍流能量级串机理及其格子气数值模拟的研究”（编号：40405004 ）；国家自然科学基金重点
项目“非均匀地表通量与大气边界层过程的研究”（编号：40233030 ）；国家自然科学基金项目“非均匀下垫面上中尺度通量
参数化的研究”（编号：40605006 ）资助.
作者简介：程雪玲（1971-），女，辽宁大连人，副研，主要从事大气边界层湍流理论以及湍流扩散的数值模拟等方面的研究.
E - m ail ： chengxl ＠ m ail. lapc. ac. cn
掀起了元胞自动机研究的新高潮［1］。
一般公认， Hardy 等［2］于20 世纪70 年代建立
的所谓 HPP 格子气模型实际上是元胞自动机，这个
模型是由简单、具有全离散动力特性的粒子构成，粒
子以保持动量守恒和粒子数守恒的方式，在二维方
形网格上运动和碰撞。
然而，出乎最初的预想结果，格子气模型并不优
于流体力学的传统数值方法。虽然这种方法在计算
流体运动时，由于其布尔（ Boole ）性质，可以不做什
么近似，如对于多体碰撞项中的多体相关函数可以
直接表达出来，而不用将其分解成单粒子密度函数
之积，因此在计算机上能够精确地进行模型的演化。
但是也正是由于其处理的是布尔量和碰撞规则的布
尔性质，使得一方面在对布尔量求统计平均时产生
噪音，另一方面用布尔性质对模拟的物理系统的参
数的调节欠灵活，如用碰撞规则来确定的粘滞性通
常都很高，且不可调节，因此造成实际研究的困难。
20 世纪80 年代末， Mc Nam ara 等［3］、 Higuerasi 等［4］
证明了把自动机的布尔动力学扩展到直接处理实数
（表示元胞具有指定状态概率）的优点。且 Mc Nam
ara 等［3］在他们的工作中首先研究了将元胞自动
机思想应用于求解玻尔兹曼方程的数值模拟，这也
就是现在广泛用来模拟各种流动的格子玻尔兹曼方
法。
早在18 ～19 世纪，随着欧洲的工业革命带动了
热机的发展，建立了热力学第一定律、热力学第二定
律，而对于后者的研究更是由于玻尔兹曼等人致力
于探究热力学底下的微观层次中的分子机制，而促
使统计物理学的蓬勃发展［5］。1857 年克劳修斯发
表了题为《论热运动形式》一文，首次引进了统计概
念.他指出单系统的宏观性质只是大量分子运动的
平均效果。1860 年麦克斯韦在《气体动力论的说
明》一文中，首先认识到分子的速度各不相同而得
到了速度分布律———麦克斯韦速度分布律。奥地利
物理学家玻尔兹曼把分子运动论理论的研究推向了
高峰。玻尔兹曼在接受了麦克斯韦速度分布律的规
律后，进而考虑到麦克斯韦的分布律，只反映了气体
平衡态下的情况，其中没有时间因子，他试图建立起
非平衡态分布的运动方程，在1872 年，他得出了速
度分布函数必须满足的非线性积分微分方程，即著
名的玻尔兹曼（ Boltzm ann ）方程。而玻尔兹曼基于
分子运动的统计特性研究宏观现象的思想，随着计
算机技术的发展，新的计算方法的出现，更是从热力
学研究领域向流体力学等更广阔的物理学科扩展，
我们这里所介绍的格子玻尔兹曼方法即是用于研究
流动现象的一例。
2 格子玻尔兹曼方法
2.1 玻尔兹曼方程
假设所考虑的系统是稀薄气体，可忽略三体及
三体以上的碰撞，而只考虑二体碰撞，且假设是弹性
碰撞。再假定气体是单一成分的单原子分子，即没
有源项，则粒子在相空间（珒r，珒v）附近的小体积元d珒rd珒v
中的粒子数fd珒rd珒v，它随时间的变化率5
5t
（fd珒rd珒v）由
以下部分组成：
（1） 粒子在坐标空间中运动造成的变化率［ 5
5t
（fd珒r d珒v］1，由于速度在珒v 附近d珒v 范围内分子的运动
所造成d珒r 体积元中粒子数的增加率为：－ dv ∮
σ
f珒v·
dσ，其中σ表示小体积的表面，dσ是小面积元，方
向沿表面外法向。利用 Gauss 定理有：
－ d珒vo ∮
σ
f珒v·dσ ＝－ dv ∫5
5r
·（f珒v）d珒r
由于珒v 和珒r 都是独立变量5
5r
·（f珒v） ＝珒v· 5
5r
，又考虑
到整个体积就是d珒r，便有：
［ 5
5t
（fd珒r d珒v）］1 ＝d珒vd珒r·珒v· 5
5珒r
（1）
（2） 粒子在速度空间中迁移造成的变化率［ 5
5t
（fd珒r d珒v）］2，这是因为在相空间中，坐标中包含有速
度坐标珒v，所以才有这一项，和上一项同样方法可以
简化为，
［ 5
5t
（fd珒rd珒v）］2 ＝d珒vd珒r 5
5珒v
·（f
珗F
m
） （2）
其中珗F 是作用于粒子的外力，m 是粒子质量。如果
5
5v
· 珗F ，则有，
［ 5
5t
（fd珒rd珒v）］2 ＝d珒vd珒r
珗F
m
· 5f
5v
（3）
（3） 粒子间相互碰撞造成的变化率5f
5tc
d珒rd珒v，设珒v 与
珒v2 是
碰
撞
前
两
分
子
的
速
度
，
碰
撞
后
速
度
分
别
为
珒
v′，
珒v′2
，则碰撞项可以写成：
（5f
5t
）c ＝∫dΩσ（｜ g ｜，θ） ｜ g ｜
［f（珒v′）f（珒v′2）－ f（珒v）f（珒v2）d珒v2 （4）
30 地球科学进展 第22 卷
其中，f（珒r，珒v，t）称为单粒子分布函数，g 是两个分子
的相对速度，θ为g 与质点系中心线的夹角，σ为微
分散射截面，Ω为固体角。由于有：
5
5t
（fd珒r d珒v） ＝［ 5
5t
（fd珒r d珒v］1 ＋［5
5t
（fd珒r d珒v］2
＋（5f
5t
）cd珒r d珒v （5）
将式（1）、（3）和（4）代入到式（5）中，得到：
5f
5t
＋珒v·5f
5珒r
＋
珗F
m
·5f
5珒v ＝∫dΩσ（｜ g ｜，θ）
｜ g ｜（ f′ f′2 － ff 2）d珒v2 （6）
这就是稀薄气体的输运方程，也称为玻尔兹曼输运
方程［6，7
］。
考虑系统不受外力作用的简单情形，此时玻尔
兹曼方程简化为：
5f
5t
＋珒v·èf ＝（5f
5t
）c （7）
但该方程当时无法求解，直到40 年后（1916 —
1917 ） S. Chapm an 和 D . Enskog 才第一次得出在一
般情况下求解的一种方法。这是一种将分布函数展
开成努森（ Knudsen ）数的幂级数的逐级近似方法。
他们假设系统虽未达到平衡，但在每一小局域的体
积中已离局部平衡不远。所谓局部平衡，是指分子
的速度服从 Maxwell- Boltzm ann 分布：
f0（珒r，珒v，t） ＝n（珒r，t） m
2πkB T（珒r，t
( )
）
3/2
exp －
m v2
2kB T（珒r，t
( )
）
（8）
其中n（珒r，珒v，t）＝ρ（珒r，t）/m 。
假设f可在f0 附近展开，
f ＝ f 0 ＋Knf 1 ＋Kn 2f2 ＋⋯ （9）
Kn ＝τ
Df/Dt
f n 1，为努森（ Knudsen ）数，表示2 个基
本的时间尺度（当以平均自由程和流体空间尺度表
示时为长度尺度）之比，也就是说：分子碰撞时间尺
度τ（或者与湍流涡的相互作用相关的特征碰撞时
间尺度τ turb ），和由于平均流的非均匀性形成的偏离
当地值的对流时间尺度之比， D / D t是实导数。
根据其一级近似所得的流体动力学方程，就是
通常适用于连续介质的 Navier- Stokes （NS）流体动力
学方程组。但是这之后，玻尔兹曼方程并未能象NS
方程那样，随着计算机技术的发展，在数值模拟方面
得到广泛的应用，直到20 世纪60 年代元胞自动机
理论的出现，在此基础上，科学家们成功地进行了首
批流体力学格子气模拟，为玻尔兹曼方程的格子化
提供了数学手段。
2.2 格子玻尔兹曼方法
Hardy 等所建立的 HPP 格子气模型的物理过程
实际上与玻尔兹曼方程所反映的物理过程相似，其
基本要素是：①空间离散，将流体存在的空间划分成
间距为1 的正方形网格（图1a）；②流体离散，将流
体想象成由许多只有质量（单位质量），没有体积的
微小粒子组成，每个粒子可向4 个方向之一运动；③
时间离散，有了流体及其存在的空间，它们的相互作
用就称为时间作用步，分为粒子在每个节点上的相
互碰撞和从一点到邻近点的移动（图1b）；④碰撞规
则，当且仅当只有两个粒子沿相反的方向到达某节
点时（对头碰撞），它们沿另外的两个方向离开该节
点，其他情形则直接穿透（图1c）；⑤ Pauli 不相容原
理，在同一时刻同一网格点上，每一个速度方向最多
只允许有一个粒子，它为 Boole 型表示法提供了可
能，可以用“1”表示有粒子，“0”表示没有粒子。
有了以上5 条要素，若再给定边界上的碰撞规
则和初始状态，系统就变成了一个离散化 Boole 化
（a） HPP 模型网格（b）时间步作用过程
（c）碰撞规则
图1 HPP 模型
Fig.1 HPP m odel
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 31
动力系统———格子气自动机［7］（ Lattice Gas Autom ata，
LGA ）；它可以随着分离的时间步而运行。
如前所述，如果我们考虑系统不受外力作用，简
化的玻尔兹曼方程为：
5f
5t
＋珒v·èf ＝（5f
5t
）c （10）
考虑我们所讨论的格子气模型，可以知道有2
个因素导致f 随时间变化：① 粒子从一个节点到另
一个节点的运动；②在网格节点上粒子间的相互碰
撞。实际上粒子间不断发生的碰撞，对f 的影响非
常复杂，我们只引入一个简单的碰撞项Ω，得到格子
玻尔兹曼方程：
f1（ x ＋ e i，t ＋1）－ fi（x，t） ＝Ωi （11）
下标i表示粒子的某个运动方向。
在格子玻尔兹曼方程中，如果直接计算碰撞项，
计算量相当巨大，对这个问题的解决是由 Higuera
及其合作者［4，8
］提出来的。早在1954 年， Bhatnager
等［9，1
0］曾假设，那些描述涡与涡之间相互作用的真
实的（也是未知的）碰撞作用的具体细节并不重要，
他们提出了一个简化的碰撞模式，将玻尔兹曼方程
简化为：
5f
5t
＋珒v· 5f
5x
＝－
1
τ
（ f－ feq） （12）
在这种形式的玻尔兹曼方程中，碰撞表达式中包含
了在特征时间尺度τ内达到局地平衡分布τ这一松
弛项，该模式被称为 BGK 模式。它描述了分子相互
作用的物理本质，τ为分子碰撞时间，又可称为松弛
时间。尽管 BGK 模式似乎只适用于局地近平衡态，
但人们认识到只要松弛时间的长短能包含相关的物
理特点，这样一种近似将超越其理论局限得到应用。
在湍流的 BGK 模式中，τ可用一个典型的湍流松弛
时间τ turb 代替。 Higuera 等人根据 BGK 模式将格子
玻尔兹曼方程写为如下形式：
fi （ x ＋ e i，t ＋1）－ fi （x，t） ＝ 1
τ turb
［feq
i（x，t）－ fi （x，t）］
（13）
其中，局部平衡解feq
i 是密度Σfi 和动量Σfivi 的给
定函数。可以选择函数feq
i ，以便产生相应的物理现
象，更一般意义上feq
i 可以包括其他物理性质，例如
局部温度等。此方程被称为 LBGK 方程，也有人将
格子玻尔兹曼方法称作 LBGK 模式。
3 格子玻尔兹曼方法的应用
格子玻尔兹曼方法是基于元胞自动机的并行思
想得以实现的，一方面具有元胞自动机的由简单演
化规则显现复杂动力特性的特点，在研究动力系统
时得到丰富的物理图象，从而带来对物理问题深刻
的认识；另一方面由于以玻尔兹曼方程为基础，使其
在模拟流动现象时比元胞自动机更符合实际情况。
格子玻尔兹曼方法被广泛用来模拟各种流体流
动，并被确信是胜过计算流体动力学传统数值技术
的极好方法，其微观水平描述给出数值方法的自然
解释，使之推广到复杂的流动问题，诸如二相流、磁
流体动力学、孔隙介质中的渗流或热流体动力学，以
及湍流运动的数值模拟等等。
3.1 格子玻尔兹曼方法在各种物理领域的应用
元胞自动机最初是被用来模拟非生物过程的自
繁殖系统［11，1
2］，由元胞自动机规则显现出非常复杂
且意想不到的行为。而这之后，又有许许多多的规
则被提出来［1］，如模拟表面生长过程的规则，模拟
森林火灾的概率元胞自动机规则，以及模拟像砂粒
一样的颗粒的基本堆积和倒塌现象的砂堆规则。S.
W olfram 更是将规则按动力学分类，进行了系统的
研究［13－14 ］。而这些规则也从单纯的理论研究到应
用于实际问题中，如 Vichniac 在20 世纪80 年代提
出的相变统计模型（ Ising 模型）的自旋动力 Q2 R 规
则［15，1
6］。
元胞自动机的应用范围越来越广泛，比如扩散
现象、材料晶体凝聚时的生长过程、岩石等非均质材
料的弹、脆、塑性等性质的解释，城市土地利用的演
化、道路交通流的模拟、对地震的模拟、林火的蔓延、
流行病毒的传播、甚至舆论传播等等方面。
扩散现象在物理学、化学和生物学等许多领域
中起着重要作用，其中有些过程非常复杂，难以用扩
散方程的解来描述，而元胞自动机方法非常自然、很
好地适合于进行数值模拟。如产生两种介质接触间
的粗糙界面，对扩散前锋的性质进行研究。生长也
是一种扩散现象，一般是在远离平衡态的时候形成
的分形结构，比如由于曲率原因驱动扩散进行进而
引起晶粒长大的现象、过冷介质中的树枝状凝固现
象等（图2）。现在，这一方法在研究金属材料的再
结晶［17］、金属凝固过程的枝晶生长［18］等方面得到
了广泛的应用。
由此，我们也可以看到，元胞自动机的优势在于
能够通过规则对各种复杂的传播过程进行模拟，除
了前面介绍的扩散、生长现象之外，还有象岩石的破
裂［19］（图3）、地震波传播过程的模拟［20，2
1］。
32 地球科学进展 第22 卷
图2 二维元胞自动机模拟枝晶的生长（黑色）［1］
灰色点代表尚未凝聚的扩散粒子
Fig. 2 Num erical sim ulation of dendrite deve lopm ent
（ black ） by the 2 D cellular autom aton ［1］
In figure ， grays are un- agglom eration grains
除此之外，利用元胞自动机具有模拟复杂系统
的优势，人们还进行了各种社会现象的研究，如通过
分析在人群集中的地方人员的自适应流动，找到出
现瓶径现象的规律，从而为指导城市交通建设提供
理论依据［22］（图4）；对舆论传播这一社会复杂系统
的演化过程进行模拟，利用不同元胞给定值代表对
舆论持不同态度的人，通过一定规则演化之后，得到
给定值接近的元胞聚集在一起，与人们常说的“物
以类聚，人以群分”的社会现象完全吻合［23］（图5）。
在另一方面，元胞自动机的元胞空间与基于几
何定位的地理空间概念有类似的意义，用元胞状态
的改变就可反映出地理特征的改变情况，在土地状
况日益受到关注的今天，利用元胞自动机研究土地
覆盖的变化［24］、人地关系［25 ～27 ］、城市发展规划［28］
等问题取得了众多有价值的成果。
图3 用模型获得的岩石裂隙的不同形式
Fig.3 Num erical sim ulation of the different form s of the rock cranny
图4 人员流动过程中的成拱现象［22］
Fig.4 The arching phenom enon in crow d
pedestrian flow ［22］
3.2 格子玻尔兹曼方法在流体计算中的应用
20 世纪70 年代，基于元胞自动机思想提出的
HPP 模型是能够进行流体运动计算的第一个格子
气模型，但实际上它最初是用于研究气体交互作用
粒子的基本统计性质的。而随着元胞自动机理论的
发展，到80 年代末，人们已预计全离散计算模型将
替代风洞试验，使人们看到这一前景的格子气模型
图5 模拟得到“物以类聚，人以群分”［23］
Fig.5 Sim ulation result on things of one
kind com e together ［2 3］
则是1986 年由 U . Frisch 等［29］以及 S. W olfram ［30］
几乎在同时提出的著名的 FHP 模型。此后，利用格
子气模型进行了一系列的流动现象的数值模拟［31］，
比如，卡门涡阶和有台阶的回流区计算、自由边界问
题、格子气的 MHD （磁流体动力）模型、表面张力问
题、多孔介质问题、悬浮体模型、传热问题以及湍流
问题等。进一步，为了消除格子气模型布尔量求统
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 33
计平均时产生噪音，以及调节模拟的物理系统的参
数欠灵活等问题，研究人员引入了 Maxwell- Boltzm
ann 速度分布函数代替格子气中的布尔量，这种格
子玻尔兹曼方法能够避免数值的不稳定性，并可以
模拟高 Reynolds 数流体的流动，目前被广泛用来模
拟各种流动问题［32，3
3］，如多相流［34 ～36 ］、孔隙介质中
的渗流［37 ～40 ］，模拟风传输雪的过程等（图6），格子
玻尔兹曼方法的微观水平描述在计算这种复杂流动
现象方面具有传统计算方法不可比拟的优势。
除此之外，利用格子玻尔兹曼方法计算的流动
图6 沟道中雪沉积的三个连续阶段［1］（风从左向右吹）
Fig.6 Three continuous snow aggradation m om ents in street canyon ［1］（ W ind blow s from left ）
图7 不同排列的双圆柱绕流涡线图［41］
Fig.7 Vorticity contours in the w ake flow fie lds of various array cylinders ［41］
图8 建筑物内外的流线［44］
Fig.8 Stream lines around the building ［44］
图9 Pow er Flow 计算汽车外部流线①
Fig.9 Stream lines around the car by Pow er Flo w *
问题也越来越广泛，比如，绕流运动［41 ～43 ］、室内流
动［44］（图8）、可压缩的高速流动［45］、浅水涌浪在固
体表面的反射和绕射现象［46］等。
近年来，格子玻尔兹曼方法在数值计算领域越
来越受到重视，并被用于研发形成商业 CFD 软件，
Power FLOW ′就是由美国 EXA 公司开发的基于格子
玻尔兹曼方法的 CFD 软件，在短短不到20 年的时
间里，实现了人们的期望：用全离散计算模型将替代
风洞试验（图9）。
除了这些工程上的应用之外，在大气领域的计
算也初步取得了成果，比如利用旋转坐标系计算
Ros*****y 涡［47］、计算正压大气运动［48，4
9］以及模拟重
力波［50］等。
3.3 格子玻尔兹曼方法研究湍流
湍流也是一种流动现象，即流体的紊乱流动，在
各种不同的流动过程中都能见到，如汽车或飞机周
围的空气流动，大气运动等，我们之所以将其作为一
个单独的题目介绍，在于它的重要性和困难性。
所有的湍流运动可用一组非线性偏微分方程来
描述，这些方程大约200 年前由 Navier 和 Stokes 首
先提出来。 Navier- Stokes 公式的数学求解非常复
杂，一个多世纪来，人们无法得到方程的解，特别是
湍流运动的解［51 ～53 ］。
3.3.1 格子玻尔兹曼方法封闭NS 方程［54］
当湍流发生时，随着能量的级串，能量由外界传
34 地球科学进展 第22 卷
① http ：∥ www . exa. com /
给大尺度涡L ，然后向小尺度（大约η≈Re－3/4 L 量
级，η称为 Kolm ogorov 长度）涡传递，并最终耗散掉。
采用直接数值模拟（ DNS ）求解湍流问题时，必须对
所有被激发的自由度同时求解。但对于多数实际情
况，这个计算量相当巨大，如计算经过汽车的流动，
雷诺数大约是106 ，至少需要计算1014 个自由度，远
远超出现有的最快的计算机的计算能力；因此，对于
高雷诺数流动必须建立模式。
Boussinesq 提出可以将小尺度场u′对大尺度场
U 的平均作用以类似分子撞击大尺度流场的方式来
描述；也就是说，湍流的小尺度是通过涡粘性系数
ve 对大尺度作用的。按照这种思路，湍流流动的大
尺度常用 Navier- Stokes 方程描述，不过这里v 则用
修正的有效粘性系数vt ＝v ＋ve 所代替。
1941 年， Kolm ogorov ［55］提出当雷诺数很大时，
湍流的统计特性具有普适性，也就是说，它们只依赖
湍能的耗散率。因此，在湍流运动中，湍能从大尺度
（外部作用如加热，可认为是原动力）向非常小的尺
度传递，普适性意味着存在一个相当宽的尺度范围，
其中能量传递的统计特性不依赖（或稍微依赖）于
大尺度流动特性。因此，对于Δn L，人们可以将小
于Δ的尺度的动力学特性参数化，得到有效粘性系
数vt（Δ），就可以正确反映能量向小尺度的传递；这
简化了计算，对计算机的计算能力要求也降低了。
然而，对涡粘性系数的基本假设及修正，依赖于
已存在的尺度分离，即可解的大尺度的运动与那些
不可解的脉动涡的分离。但是，由于缺乏尺度间的
明确的划分，使得涡粘性模型的有效性受到质
疑［56，57］。
格子玻尔兹曼方法为解决湍流问题提供了一个
不同的途径，在此方法中，尽管大的湍流涡仍被看作
存在于涡的海洋中，但却不再依赖于尺度分离的假
设。如前所述， Navier- Stokes 方程可通过 LBGK 方
程以幂Kn 展开的低阶截断得到，当Knn 1 ，在两个
基本尺度间存在尺度分离，格子玻尔兹曼方程（式
13）包含了涡粘性模式，湍流涡粘性系数正比于
τ turb ；当Kn 不很小且不存在尺度分离时，τ turb 表示各
尺度湍流涡粘性的重整化求和，本质上优于 Navier-
Stokes 方程中依赖于尺度分离的涡粘性系数。下
面，用此方程导出更普适的涡粘性模式，即确定湍流
脉动的自松弛时间τ turb
［54］。
τ turb 应依赖于反映不同湍流机制的时间尺度，
τ turb ＝τ0 ＋Ψ（k /ε，S－1，⋯） （ 14）
其中，τ0 是基本的分子松弛时间，k 是湍流动能，ε
是湍流耗散率，S 是当地速度梯度的度量。这里，k /
ε和S－1 是湍流特征时间尺度。由于存在其他的湍
流作用，比如浮力和旋转，方程14 中可能还要包括
其他的固有时间。因此，方程中Ψ的典型的函数形
式应当是湍流时间尺度的谐波平均（也就是说是湍
流频率的叠加）。通过推导得到有效松弛时间τ turb ，
τ turb ＝τ0
＋C μ
k2 ε
T（1 ＋η2 ）1 /2 （15）
其中，C μ＝0.085 ，η＝ Sk / ε，T 是绝对温度。在
Kn→0 时，即两个基本尺度间存在尺度分离，可导出
基于重整化群的k－ ε涡粘性模型。
应用有效松弛时间τ turb 进行计算时，不必调整
模式参数。图9 是应用此方法计算流动经过一个实
际的小轿车。
3.3.2 格子玻尔兹曼方法直接数值模拟
在上一节，我们介绍了 Kolm ogorov ［58］提出当雷
诺数很大时，能量由外界输入给大尺度涡，然后向小
尺度涡传递，并耗散掉。由此可见，涡旋是流体运动
的肌腱，它的生成、发展和演化，以及它与外部流动
和物体之间的相互作用，影响了整个流场的特性。
涡旋研究的对象可以来自于解析解，如 Ossen 涡、
Taylor 涡、 Burgers 涡等，这些具有解析解的涡结构都
比较简单，对这些简单涡旋结构的研究是分析复杂
涡旋流动的基础［59］。而在20 世纪50 ～60 年代，人
们在湍流场中又发现了拟序结构（ Coherent structure
）［60］。拟序结构与涡旋有着极为密切的关系，一
般认为拟序结构是存在于湍流场中的具有较大横向
尺度的有规律的涡旋结构，但严格的说，只要是有规
律的涡旋结构，不论尺度多大，都可以看成是拟序结
构。湍流就是由拟序结构和随机小涡背景流场共同
构成的。
面对湍流这一复杂的现象，人们最早采用的是
实验的方法，如雷诺实验，但这种流动显示实验很难
进行定量的分析。之后的热线风速仪以及超声风速
仪，可以得到流场中单点速度的定量结果，但对流场
有一定影响。近年来激光多普勒测速技术可以对流
场进行非接触性的测量，但成本较高。20 世纪80
年代初，一种名为 PIV （图像粒子测速）的流场测量
技术逐渐发展起来，它的原理是通过向流场中投放
示踪粒子，然后测量和流体混合均匀的粒子的运动
速度来间接测量粒子周围流体的速度。 PIV 技术的
特点是它不再是单点速度的测量，它可以测量一个
平面内二维的流场，使人们可以得到流动的更多的
信息。
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 35
除了这些实验手段，随着计算机存储容量和计
算速度的飞速提高，数值模拟作为一种研究手段也
相应发展起来，利用数值模拟的方法，人们可以得到
整个三维流动速度场的时间序列，而且其成本比实
验要低很多，人们已经利用数值模拟的方法对湍流
流动的各种统计特性进行了许多的研究［61，62］。但
以往的直接数值模拟都是基于求解NS 方程，常用
的方法有差分法、有限元法和谱方法等，由前面的分
析我们知道，其计算量是相当巨大的。而格子玻尔
兹曼方法由于其基于分子运动方程求解，程序设计
非常简单，还有一个优点是易于并行化，是实现这样
庞大的计算问题的一个重要保证。
用格子玻尔兹曼方法进行直接数值模拟，最开
始也是计算周期边界条件的湍流演化过程，并与传
统的谱方法进行比较，比较的结果两者在计算精度
上并无太大差别［32，3
3，6
3，6
4］。之后，还用这种方法计
算槽道流，研究其近壁面的特性，以及直接数值模拟
Rayleigh- Taylor 流［65］、两种流体界面上 Kelvin －
Helm holtz 波［66］等图10。最近，还出现了一种熵格
子玻尔兹曼方法［67］，模拟等温条件下的运动，其中
可调节直接数值模拟的截断尺度，只计算大尺度的
运动。
图10 直接数值模拟 Kelvin - Helm holtz 波［6 6］
Fig.10 Direct num erical sim ulation of Kelv i n－ Helm holtz w av e［66］
3.3.3 格子玻尔兹曼方法研究大气湍流
大气湍流通常比实验室里的湍流要更加复杂，
在大气运动中，除了由于摩擦产生的机械湍流外，还
有由于浮力产生的对流热泡造成的湍流运动，由于
地面吸收太阳辐射后升温导致热气上升，形成对流，
产生湍流运动，此时的湍流运动是硬湍流运动，人们
发现在硬湍流条件下，温度场级串有别于速度
场［68］。
湍流能量级串模式一直是湍流研究中的前沿问
题，对此问题的认识将最终影响湍流模式的参数化。
目前，已经出现了各种级串模型，但对于能量级串这
一问题的研究还有待于进行多方面深入的研究。
Kolm ogorov 理论认为从大涡向小涡的逐级级串
具有自相似结构（由 Richardson 首先提出），它们充
满整个空间，该理论（又称 K41 理论）给出的结构函
数及能量耗散率的标度律在高阶矩时与实验结果明
显不符。 Landau ［69］首先向 K41 理论提出了质疑，指
出将能量耗散律ε看作常数缺乏物理根据，应当考
虑ε的涨落；不久之后，实验发现了湍流的间歇性，
湍流研究者才认识到间歇性对建立速度结构函数标
度律的重要意义，并出现了一系列计入间歇性的一
些重要模型［70 ～73 ］，如β－模型，分形级串模型， GOY
壳模型，对数－Poisson 模型和 S- L 模型等，理论预
测和实验结果在较高阶矩时相当一致。这是由于
时—空间歇性的存在意味着湍能从大涡（大尺度）
向小涡（小尺度）的逐级传递过程，既不充满空间，
也不能看作是平稳随机过程。
最近，一种新的级串过程———同步级串模式被
提出来［74］图11，该模式认为从外界向含能区注入
能量到形成发达湍流其实并不需要经过n 次有序的
图11 湍流能量级串模式
Fig.11 Turbulent energy cascade m odel
36 地球科学进展 第22 卷
级串过程，正如一块玻璃被摔碎时，最可能的情形是
破碎成各种不同尺寸的碎块，含能涡一次散裂同时
形成大、中、小、微各种尺度的湍涡的概率也应远大
于 Kolm ogorov 级串模式的概率。
综上所述，经过湍流研究工作者几十年的积累，
关于湍流速度场的统计特征已经了解得较为深入，
但是对作为被动保守物的标量场（例如温度场）了
解得还很不够［75 ～77 ］。人们过去一度认为标量场应
具有与速度场相同概率分布，但是后来从观测事实
中发现二者之间有很多的不同，标量场的概率分布
形态更显复杂。20 世纪80 年代美国芝加哥大学在
封闭容器中利用氦气做热对流实验（文献中称为
“芝加哥实验”）［78 ～80 ］，结果表明湍流温度脉动场的
概率分布类型与外参数 Rayleigh 数（Ra）有密切关
系。存在一临界Ra 值（≈108），即当 Ra ＜10 8 时，湍
流温度场中两点温度差的概率密度分布非常接近于
Gauss 分布（软湍流）；但当 Ra ＞10 8 时，概率密度分
布变成非 Gauss 分布，而更接近于指数分布（硬湍
流）。Ra 是一个反映浮力作用与热扩散作用相对的
无量纲数，定义为：
Ra ＝
αgδL3
κν
（16）
式中α表示热膨胀系数，g 是重力加速度，δ即驱动
对流的外部温差，L 代表直径或高，v 和κ分别为分
子粘性系数和热扩散系数。到目前为止，还有其它
很多实验室模拟和数值模拟工作也进一步证实了芝
加哥实验的结果［76，7
7］。
由此可见，对于包含热运动的大气湍流，能量级
串过程将会更加复杂，而通过格子玻尔兹曼方法直
接数值模拟这种湍流运动，进行机制的研究应是非
常有意义的。利用格子玻尔兹曼方法对温度场的研
究已有很多［81，8
4］，但都还限于平均场的模拟，可以
先用格子玻尔兹曼方法直接数值模拟含有热对流运
动的 Rayleigh- Bernard 流，了解热对流的湍流运动，
然后模拟更加真实的大气湍流运动。
4 结 语
通过前面的介绍，我们对格子玻尔兹曼方法的
产生和应用有了比较全面的认识：格子玻尔兹曼方
法是基于20 世纪60 年代的并行思想及由此产生的
元胞自动机，同时建立在分子运动论思想上的玻尔
兹曼方程为其物理基础。除了用于模拟各种复杂物
理、社会现象之外，在对流动现象的模拟方面，格子
玻尔兹曼方法一方面由于其分子运动论的物理背
景，使得其在模拟流体流动时不会出现元胞自动机
那样由于布尔运算而产生的噪音，所以广泛应用于
流动现象的数值模拟；另一方面由于其微观水平描
述给出数值方法的自然解释，对于模拟流动中的复
杂现象，如湍流问题要优越于NS 方程，目前在这方
面的探索已取得了许多成果，使得对湍流这一古老
的问题的研究又增添了一个有力的工具；并尝试分
析了利用这一方法研究大气湍流的初步方案。
参考文献（ Referem ces ）：
［1］ Chopard B ， Droz M . Cellular Autom ata Modeling of Physica l System
s ［M］. London ： Cam bridge University Press ，1998.［祝玉学，
赵学龙译. 物理系统的元胞自动机模拟［M］. 北京：清华大学
出版社，2003.］
［2］ Hardy J， Pom eau Y ， Pazzis O De. Molecular Dynam ics of a classical
lattice gas ： Transportproperties and tim e correlation fu nctions
［J］. Phys. Rev. A ，1976 ，13：1 949-1 960.
［3］ Mc Nam ara G G ， Zanetti G . Use ofthe Boltzm ann Equation to Sim -
ulate Lattice Gas Autom ata ［J］. Physical Review Letters ，1988 ，
61：2 332-2 335.
［4］ Higuera F ，J im enez J， Succi S. Boltzm ann approach to lattice gas
sim ulations ［J］. Europhysics Letters ，1989 ，9：663 .
［5］ Feng Duan ， Feng Shaotong. The W orld of Entropy ［M］. Bei jing ：
Science Press ，2005.［冯端，冯少彤. 溯源探幽———熵的世界
［M］. 北京：科学出版社，2005.］
［6］ Harris ， Stewart. An Introduction to the theory ofthe B oltzm ann E -
quation ［M］. New York ： Holt ， Rinehart and W inston ，1971.
［7］ Li Yuanxiang ， Kang Lishan ， Chen Yuping. Lattice Gas Autom aton
［M］. Bei jing ： Tsinghua University Press ，1994.［李元香，康
立山，陈毓屏. 格子气自动机［M ］. 北京：清华大学出版社，
1994.］
［8］ Higuera F ，J im enez J， Succi S. Lattice gas dynam ics with enhanced
collision ［J］. Europhysics Letters ，1989 ，9：345.
［9］ Quian Y H ， Hum ieres D ， Lallem and P. Lattice BGK Models for
Navier- Stokes Equation ［J］. Europhysics Letters ，1992 ，17（6）：
470-484.
［10］ Bhatnager P ， Gross E P ， Krook M K . A m odel for collision
process in gases ［J］. Physics Review ，1954 ，94：511.
［11］ Pesavento U . An im plem entation of von Neum an n ＇ s self- reproducing
m achine ［J］. Artificial Life ，1995 ，2：337-354.
［12］ Yan Guangwu. Advances of the Studies on the Ce llular Autom ata
and Artificial Life ［J］. J ournalofJ ilin University ，2003 ，41（1）：
40-44.［闫广武. 元胞自动机与人工生命研究进展［J］. 吉林
大学学报，2003 ，41（1）：40-44.］
［13］ W olfram S. Statisticalm echanics ofcellula r autom ata ［J］. Review
of Modern Physics ，1983 ，55（3）：601-644.
［14］ W olfram S. Com putation theory of cellular au tom ata ［J］. Com -
m unication in Mathem aticat Physics ，1984 ，95：15-57.
［15］ Vichniac G . Sim ulating physics with cellula r autom ata ［J］. Physica
D ，1984 ，10：96-115 .
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 37
［16］ Yu Lu ， Hao Bolin ， Chen Xiaosong. Phase Transitions and Critical
Phenom ena ［M ］. Bei jing ： Science Press ，2005.［于禄，郝柏
林，陈晓松. 边缘奇迹———相变和临界现象［M］. 北京：科学
出版社，2005.］
［17］ Zheng Chengw u ， Lan Yong jun ， Xiao Nam in ， et al . Mesoscale
sim ulation of Static recrystallization ofho tdeform ed austenite in a
low carbon steel ［J］. Acta Metallurgica Sinica ，2006 ，42（5）：
474-480.［郑成武，兰勇军，肖纳敏，等.热变形低碳钢中奥氏
体静态再结晶介观尺度模拟［J］. 金属学报，2006 ，42（5）：
474-480.］
［18］ Guo Dayong ， Yang Yuansheng. Num erical Sim ulation of Dend rite
Developm ent in Ternary Alloy Solidification ［J］. Foundry ，2006 ，
55（6）：601-603 .［郭大勇，杨院生.三元合金凝固过程枝晶生
长数值模拟［J］.铸造，2006 ，55（6）：601-603.］
［19］ W ang Shim in ， Zhu Hehua ， Feng Xiating ， etal . Influence ofheterogeneity
on m acroscopical crack form of the brittle rock ［J］.
Rock and Soil Mechanics ，2006 ，27（2）：224-227.［王士民，朱
合华，冯夏庭，等.细观非均匀性对脆性岩石材料宏观破坏形
式的影响［J］.岩土力学，2006 ，27（2）：224-227.］
［20］ Dong Qiaoliang ， Zhou Yubing ， Yao Yao. Seism ic Sim ulating by
Using Lattice Boltzm ann Method ［J］. J ournal of China University
of Geosciences （ Earth Science ），1997 ，22（6）：638-642. ［董桥
梁，周玉冰，姚姚.格子 Boltzm ann 方法地震波场模拟［J］.中
国地质大学学报：地球科学，1997 ，22（6）：638-642.］
［21］ W ang Zhenli ， Li Youm ing. Parallel Algorithm for Sim ulatin g
Seism ic W ave Propagation by Cellular Autom at a［J］. Chinese
J ournal of Geophysics ，1999 ，42（3）：410 -415.［王真理，李幼
铭，细胞自动机地震波模拟的并行化算法［J］.地球物理学
报，1999 ，42（3）：410-415.］
［22］ Zhang Peihong ， Huang Xiaoyan ， W an Huanhuan ， et al . Study
on Adaptive Cellular Autom ation Model on Crow d Pedestrian
Flow ［ J］. J ournal of Shenyang J ianzhu University ，2006 ，22
（2）：289-293 .［张培红，黄晓燕，万欢欢，等.人员群集流动
自适应元胞自动机模型研究［J］.沈阳建筑大学学报，2006 ，
22（2）：289-293.］
［23］ Liu Muren ， Deng Minyi ， Kong Ling jiang. The Cellular Autom ata
Model for Public Opinion Evolution ［J］. J ournalof Guangxi Norm
al University ，2002 ，20（2）：1.［刘慕仁，邓敏艺，孔令江.舆
论传播的元胞自动机模型［J］. 广西师范大学学报，2002 ，20
（2）：1.］
［24］ Yan J infeng ， Chen Xi. Analysis and Stim ulation Method Disc ussion
on LUCC of Arid Region Based on GIS ［J］. Arid Land Geography
，2003 ，26（2）：185 -191.［阎金凤、陈曦.基于 GIS 的干
旱区 LUCC 分析和模拟方法探讨［J］.干旱区地理，2003 ，26
（2）：185-191.］
［25］ Luo Ping ， Du Qingyuan ， Lei Yuanxin ， et al . Cellular Autom ata
Based on Geographic Feature and Urban Land Use Evolvem ent
［J］. Geom atics and Inform ation Science of Wuhan Un iversity ，
2004 ，29（6）：504 -507.［罗平、杜清远、雷元新，等.地理特征
元胞自动机及城市土地利用演化研究［J］.武汉大学学报，
2004 ，29（6）：504-507 .］
［26］ Hou Xiyong ， Chang Bin ， Yu Xinfang. Land use change in Hexi
corridor based on A - Markov m ethods ［ J］. Transactions of the
CSAE ，2004 ，20（5）：286 -291.［侯西勇、常斌、于信芳. 基于
CA - Markov 的河西走廊土地利用变化研究［J］. 农业工程学
报，2004 ，20（5）：286-291.］
［27］ Chen J ianping ， Ding Huoping ， W ang Gongwen ， et al . Desertification
Evolution Modeling through the Integr ation of GISand Cellular
Autom ata ［J］. J ournal of Rem ote Sensing ，2004 ，6（3）：
254-260.［陈建平，丁火平，王功文，等.基于 GIS 和元胞自动
机的荒漠化演化预测模型［J］. 遥感学报，2004 ，6（3）：254-
260.］
［28］ Cui W eihong ， Zhang Xianfeng. Study on Dynam ic Monitoring
and Dynam ic Sim ulating of Land Resource ［J］. Geo- Inform ation
Science ，2002 ，1：79-85.［崔伟宏，张显峰.土地资源的动态
监测和动态模拟研究［J］.地球信息科学，2002 ，1：79-85.］
［29］ Frisch U ， Hasslacher B ， Pom eau Y . Lattice gas autom ata for the
Navier- Stokes equation ［J］. Phys. Rev. Lett ，1986 ，56：1505.
［30］ W olfram S. Cellular autom ation fluid ： Basic theory ［J］.J. Stat.
Phys ，1986 ，45：471 .
［31］ Qian Yuehong ， D d ＇ Hum ieres ， Y Pom eau ， et al . The newest developm
ent of Lattice Gas Autom aton Dynam ics ［J］. Mechanics in
Engineering ，1990 ，12（1）：7-16.［钱跃宏，D d’ Hum ieres ，Y .
Pom eau ，等.格子气流体动力学及其最新进展［J］.力学与实
践，1990 ，12（1）：7-16.］
［32］ Benzi R ， Succi S ， Vergassola M . The lattice Boltzm ann equation
： theory and applications ［J］. Physics Reports ，1992 ，222
（3）：145-197.
［33］ Shiyi Chen ， Gary D Doolen. Lattice Boltzm ann Method for Fl uid
Flow s ［J］. Annu. Rev. Fluid Mech ，1998 ，30：329-364.
［34］ Shan Xiaowen ， Chen Hudong. Lattice Boltzm ann m odelfor sim u -
lating flows with m ultiple phases and com pone nts ［J］. Physical
Review E ，1993 ，47（3）：1 815-1 819.
［35］ Martys N S ， Chen Hudong. Sim ulation of m ulticom ponent fl uids
in com plex three- dim ensional geom etries by t he lattice Boltzm ann
m ethod ［J］. Physical Review E ，1996 ，53（1）：743-750 .
［36］ Zhang Raoyang ， Chen Hudong. Lattice Boltzm ann m ethod for
sim ulations ofliquid- vapor therm alflow s ［J］. Physical Review E ，
2003 ，67：066711 .
［37］ Zhan Hongren ， Li Dongyan ， W ang Guoheng. Lattice gas autom ata
m odeling on evaporation process ofporous m edia ［J］. Chem ical
Equipm ent Technology ，2005 ，26（6）：1-5.［战洪仁，李冬艳，王
国恒.多孔介质干燥过程的格子气自动机法模拟［J］.化工装
备技术，2005 ，26（6）：1-5.］
［38］ Zheng Zhong ， Gao Xiaoqiang ， Shi W anyuan. Lattice Gas Autom -
ata Modeling of Fluid Flow Through Porous Medi a［J］. J ournalof
Chem ical Industry and Engineering ，2001 ，52（5）：406-409.
［郑忠，高小强，石万元. 多孔介质中流体流动的格子气自动
机模拟［J］.化工学报，2001 ，52（5）：406-409.］
［39］ Zheng Zhong ， Gao Xiaoqiang. Lattice Gas Autom ata Method fo r
Modeling Fluid Flow and Heat Transfer in Metal lurgical Porous
Media ［J］. Acta Metallurgica Sinica ，2000 ，36（4）：433-437.
［郑忠，高小强.格子气自动机法模拟冶金多孔介质传输现象
［J］.金属学报，2000 ，36（4）：433-437.］
38 地球科学进展 第22 卷
［40］ Yue W enzheng ， Tao Guo ， Zhu Keqin. Sim ulation of electrical
transport properties in oil- water saturated porous m edia with 2 D
lattice gas autom ata ［J］. Chinese J ournal of Geophysics ，2005 ，48
（1）：189-195.［岳文正，陶果，朱克勤.二维格子气自动机模
拟孔隙介质的电传输特性［J］.地球物理学报，2005 ，48（1）：
189-195.］
［41］ Surm as R ， Santos L ， Philippi P. Lattice Boltzm ann sim ulation of
the flow interference in bluff body w akes ［J］. Future Generation
Com puter System s ，2004 ，20 ：951-958.
［42］ Takada ， Naoki ， Tsutahara ， Michihisa. Evolution of viscous flow
around a suddenly rotating circular cylinder in the Lattice Boltzm
ann Method ［J］. Com puters ＆ Fluids ，1998 ，27（7）：807-828.
［43］ Feng Shide ， Zhao Ying ， Gao Xianlin ， et al . The sim ulation of
2 D flow around a cylinder by lattice BGK m odel ［J］. Acta Mechanica
Sinica ，2002 ，34：254-258.［冯士德，赵颖，郜宪林，等.
用格子 BGK 模型模拟两维圆柱绕流［J］. 力学学报，2002 ，
34：254-258.］
［44］ Crouse B ， Krafczyk M ， Kuhner S ， et al . Indoor air flow analysis
based on lattice Boltzm ann m ethods ［J］. Energy and Buildings ，
2002 ，34：941-949.
［45］ Yu Huidan ， Zhao Kaihua. Lattice Boltzm ann m ethod for com -
pressible flow s w ith high Mach num bers ［J］. Physical Review E ，
2000 ，61（4）：3 867 -3 870 .
［46］ Liu Feng ， Hu Fei. Application of lattice Boltzm ann m eth od in
sim ulating reflection and diffraction ofsur ge waves ［J］.J ournalof
Hydraulic Engineering ，2004 ，3：22-30.［刘蜂，胡非.用格子
Boltzm ann 方法模拟涌波的反射与绕射［J］.水利学报，2004 ，
3：22-30.］
［47］ Yu huidan ， Zhao Kaihua. Ros*****y vortex sim ulation on a par aboloidal
coordinate system using the lattice Bo ltzm ann m ethod ［J］.
Physical Review E ，2001 ，64：056703.
［48］ Liu Feng ， Hu Fei. A prelim inary study on the constructio n of atm
ospheric dynam ic m odels w ith lattice Boltzm ann m ethod ［J］. Acta
Meteorologica Sinica ，2003 ，61（3）：267-274.［刘峰，胡非.
用格子 Boltzm ann 方法构建大气动力学模式的初步研究［J］.
气象学报，2003 ，61（3）：267-274.］
［49］ Liu Feng ， Hu Fei. Lattice Boltzm ann m ethod to sim ulate b arotropic
atm osphere ［J］. J ournalofthe Graduate Schoolofthe Chinese
Academ y of Sciences ，2003 ，20（2）：168-171.［刘峰，胡非，
用格子 Boltzm ann 方法模拟正压大气运动［J］.中国科学院研
究生院报，2003 ，20（2）：168-171.］
［50］ Buick J， Greated C . Lattice Boltzm ann m odeling of inte rfacial
gravity waves ［J］. Physics of Fluids ，1998 ，10（6）：1 490 -1 511.
［51］ Frisch U . Turbulence- The Legacy of a. Kolm og orov ［M ］. Cam -
bridge London ： Cam bridge University Press ，1996.
［52］ Frisch U ， Orszag S. Turbulence- challenges for theory a nd experim
ents ［J］. Physics Today ，1990 ，43：24-32.
［53］ Sreenivasan K R . Fluid Turbulence ［J］. Review of Modern Physics
，1999 ，71：383-395.
［54］ Chen Hudong ， Satheesh Kandasam y ， Steven Orszag ， et al . Extended
Boltzm ann kinetic equation for turbul ent flow s ［J］. Science
，2003 ，301 ：633-636 .
［55］ Kolm ogorov A N . Localstructure ofturbulenc e in an incom pressible
viscous fluid at very high Reynolds num ber s［J］. Dokl Akad
Nauk SSSR ，1941 ，30：301-305 （英文重印： Proceeding of the
Royal Society of London A ，1991 ，434 ：9-13.）
［56］ Yakhot V ， Orszag S. Renorm alization group analysis of t urbulence
I. Basic theory ［J］. J ournal of Science Com put ，1986 ，1：
3-51.
［57］ Yakhot V ， Orszag S ， Thangam S ， et al . Developm ent of turbulence
m odels for shear flow s by a double expans ion technique
［J］. Physics of Fluids A ，1992 ，4：1 510 -1 520.
［58］ Kuchem ann D . Reporton the IUTAM sym posium on concentrated
vortex m otion in fluids ［J］. J ournal of Fluid Mechanics ，1965 ，
21：1.
［59］ Ma Bin. The identification and visualizatio n of vortex ［D］. Beijing
： Tsinghua University ，1997.［马兵. 涡旋的识别与可视化
［D］. 北京：清华大学工程力学系，1997.］
［60］ Kline S J， Reynolds W C ， Schraub F A ， et al . The structure of
turbulent boundary layers ［J］. J ournalof Fluid Mechanics ，1967 ，
30：741.
［61］ Chen J ian. Direct Num erical Sim ulation of Sp atial Sm ooth Pipe
Flow Transition ［D］. Bei jing ： Tsinghua University ，2002.［陈
健. 空间发展光滑圆管转捩的直接数值模拟［D］. 北京：清
华大学工程力学系，2002.］
［62］ W u Minw ei. Num erical Research on Turbulent P ipe Flow ［D］.
Bei jing ： Tsinghua University ，2002.［伍敏伟.圆管湍流的数值
研究［D］. 北京：清华大学工程力学系，2002.］
［63］ Martinez D ， Matthaeus W ， Chen S. Com parison ofspectralm ethod
and lattice Boltzm ann sim ulations oftw o- d im ensional hydrodynam
ics ［J］. Physics Fluids ，1994 ，6（3）：1 285-1 298.
［64］ Yu H ， Girim a ji S ， Luo L. Lattice Boltzm ann sim ulations ofdecaying
hom ogeneous isotropic turbulence ［ J］. Physical Review E ，
2005 ，71：016708 .
［65］ He Xiaoyi ， Chen Shiyi and Zhang Raoyang. A lattice Boltzm ann
schem e for incom pressible m ultiphase flow an d its application in
sim ulation of rayleigh- taylor Instability ［J］. J ournal of Com putational
Physics ，1999 ，152 ：642-663 .
［66］ Zhang Raoyang ， He Xiaoyi ， Doolen Gary ， et al . Surface tension
effects on two- dim ensional tw o- phase Kelvin - Helm holtz instabilities
［J］. Advances in Water Resources ，2001 ，24 ：461-478.
［67］ Boghosian B ， Yepez J， Coveney P ， et al . Entropic lattice Boltzm
ann m ethods ［J］. Proceeding ofthe Royal Society of Lond on A ，
2001 ，457 ：717-766 .
［68］ Hu Fei ， Cheng Xueling ， Zhao Songnian ， etal . Hard state ofthe
urban canopy layer turbulence and its self- si m ilar m ultiplicative
cascade m odels ［J］. Science in China ，（ Series D ），2005 ，35
（ Supp. I ）.［胡非，程雪玲，赵松年，等.城市冠层中温度脉动
的硬湍流特性和相似性级串模型［J］.中国科学D 辑，2005 ，
35（增刊1）.］
［69］ Landau L D ， Lifshitz E M . Fluid Mechanics ［M ］. Oxford ： Pergam
on Press ，1987.
［70］ She ZS ，J ackson E ， Orszag SA . Interm ittentvortex structures in
hom ogeneous isotropic turbulence ［J］. Nature ，1990 ，344 ：226-
第3 期 程雪玲等：格子玻尔兹曼方法及其在大气湍流研究中的应用 39
228.
［71］ She Z S ， Leveque E . Universal scaling law s in fully dev eloped
turbulence ［J］. Phys. Rev. Lett ，1994 ，72：336-338.
［72］ Ruelle D ， Takens F . On the nature ofturbulence ［J］. Com m unication
in Mathem atical Physics ，1971 ，20：167 -192.
［73］ Benzi R ， Biferale L ， Ciliberto S ， et al . Generalized scaling in
fully developed turbulence ［J］. Physica D ，1996 ，96：162-181.
［74］ Zhao Songnian. Synchrocascade pattern in th e atm ospheric turbulence
［J］.J ournal of Geophysical Research ，2003 ，108.
［75］ She ZS. On the scaling laws oftherm alturbule ntconvection ［J］.
Physics Fluids A ，1989 ，1（6）：911 -913.
［76］ Ching E S C . Probability densities of turbule nt tem perature fluctuations
［J］. Physics Rewiew Lett ，1993 ，70（3）：283-286.
［77］ Pope SB ， Ching E SC . Stationary probability density fu nctions ：
an exact result ［J］. Physics Fluids A ，1993 ，5（7）：1 529-1 531.
［78］ Heslot F ， Castaing B ， Libchaber A . Transition to turbulence in
helium gas ［J］. Physics Review A ，1987 ，36：5 870-5 873.
［79］ Castaing B ， Gunaratne G ， Heslot F ， et al . Scaling of hard therm
al turbulence in Rayleigh- Bernard convecti on［J］. J ournal of
Fluid Mechanies ，1989 ，204 ：1-30.
［80］ W u X Z ， Kadannoff L ， Libchaber A ， et al . Frequency power
spectrum oftem perature fluctuations in free convection ［J］. Physical
Review Letters ，1990 ，64：2 140-2 143.
［81］ Feng Shide ， Zhao Ying ， Gao Xianlin ， et al . Som e progress in
the lattice Boltzm ann m odel and hydrodynam ic equations ［J］. Acta
Mechanic Sinica ，2002 ，34：254-258 .［冯士德，赵颖，郜宪
林，等. 用格子 Boltzm ann 模型的改进与流体力学方程［J］.
力学学报，2002 ，34：254-258.］
［82］ Zhang Raoyang ， Chen Hudong. Lattice Boltzm ann m ethod for
sim ulations ofliquid- vapor therm alflow s ［J］. Physical Review E ，
2003 ，67：066711 .
［83］ He Xiaoyi ， Chen Shiyi ， Doolen Gary D . A novel therm al m odel
for the lattice Boltzm ann m ethod in incom pres sible lim it ［J］.J ournal
of Com putational Physics ，1998 ，146 ：282-300.
［84］ Shi Yong ， Zhao T S ， Guo Z L. Therm allattice Bhatnagar- Gross-
Krook m odel for flows w ith viscous heat dissip ation in the incom -
pressible lim it ［J］. Physical Review E ，2004 ，70：066310.
The Applic atio n of Lattic e Boltzm ann M eth od i n
th e Atm ospheric Turbule nce Stu dy
CHENG Xue- lin g ， HU Fei ， ZHAO Song- nia n ，J IA NG J in - hua
（ State Key Laboratory of Atm ospheric Boundary Layer Physics and Atm ospheric Chem istry
Institute of Atm ospheric Physics ， Chinese Academ y of Science ， Bei jing 100029 ， China ）
Abstract ： The ob jective of this paper is to system ically introduce the lattice Boltzm ann m ethod from i ts m athem
atic and physical bases ， and investigate broadly its application in ev ery field. The lattice Boltzm ann m ethod
based on the idea ofcellular autom aton w as pro posed by J ohn von Neum ann in the 60 ＇s oflast cen tury. The cellular
autom aton sim ulates the parallelizable char acter of the brain and constructs the dynam ic e volvem ent system . That
starts a new direction to solve the com plex pro blem s. W hen it is used to solve the flow problem ， it com bines w ith
the m olecule kinetic theory of 19 century to so lve Boltzm ann transportequation. Itm eans to study the flow by sim ulating
the m olecules m ovem ent in the flow . So i t is also called lattice Boltzm ann m ethod ， or lattice gas autom aton.
From investigation ， itcan be seen thatthe application fields are v ery broad. For exam ple ， itis used to sim ulate the
crystal agglom eration ， the evolvem ent of city land use ， the traffic flow ， the seism ic w ave ， the spread of fire in the
forest ， the prevalence of virus and the spread of publi c opinion. And by this m ethod ， good results are get. In the
70＇s of last century ， Hardy ， Pom eau and Pazzis builtthe first lattice gas a utom aton to sim ulate flow . Then ， because
it can explain flow from m icrocosm ic level ， itis very suitto solve com plex flow such as m ul tiphase flow ， porous flow
and snow grain transportin the w ind. Moreover ， the com m erce softw are — Pow er Flow is developed based on the lattice
Boltzm ann m ethod. For the old problem — turbulence ， m any m ethods are used and LBM now gradually bec om es
a new w ay to directly sim ulate it because LBM ca n calculate in the m olecule level. The atm osph eric turbulence is
m ore com plex then usual turbulence. Here ， w e especially introduce the atm ospheric turb ulence character and give
an elem entary pro ject to sim ulate by LBM .
Key w ords ： Lattice Boltzm ann m ethod ； Cellular autom ata ； Lattice gas ； Direct num erical sim ulation ； Atm ospheric
turbulence.
40 地球科学进展 第22 卷
选择“Disable on www.wenxuecity.com”
选择“don't run on pages on this domain”
