（流体力学专业论文）利用格子博尔兹曼方法数值模拟水波传播.pdf

中山大学硕士学位论文 论文题目：利用格子博尔兹曼方法数值模拟水波传播 专业：流体力学 硕士生：黄广开 指导教师：詹杰民 摘要 传统的流体力学数值模拟计算方法，如有限差分法、有限体积法、有限单元法等， 均采用基于流体力学宏观方程，自上而下离散微分方程，进而使用标准数值方法进行求 解。格子博尔兹曼方法则是基于流体微观粒子运动特征，从粒子的碰撞和传播出发，设 定一系列的演化规则，经过数理统计，得到宏观上实际流体的动力特征。从微观粒子特 征的统计角度出发，模拟得到流体的运动规律。 本文回顾了水波数值模拟、格子博尔兹曼方法研究的历史，介绍了水波运动方程、 格子博尔兹曼方法的推导原理。通过比拟的方法将格子博尔兹曼方法对普通流体的模拟 移植到对水波的数值模拟中，推导得到了计算过程中所需要的水波模拟平衡态分布函 数。 最后，文章通过对方腔驱动流动模拟、全溃坝波模拟、部分溃坝波模拟这三个数值 编程实验的计算，实现了利用格子博尔兹曼方法对水波的数值模拟，通过与理论解或其 他算法的计算解的比较，验证了该方法在水波模拟中的准确性。 关键词：格子博尔兹曼；数值模拟；水波；浅水方程 利用格子博尔兹曼方法数值模拟水波传播 t i t l e ：n u m e r i c a ls i m u l a t i o no fw a t e rw a v e w i t ht 1 1 el a t t i c eb o l t z m a i u lm e m o d m 面o r ：f l u i dm e c h a n i c s n a m e ：g u a n 然a ih u a n g s u p e i s o r ：j i e i i l i nz h a n a b s t r a c t c o n 、，a l t i o n a l l y ，t h ef l u i ds i i l l u l a t i n gm e m o d sw 1 1 i c hr e s e a r c h e r su s et 0s 0 1 v en u i d p r o b l e m ss 1 1 d la sf i i l i t ed i f j f h e l l c em e t l l o d ，f i n i t ev 0 1 u m em e t l l o d ，f i l l i t ee l 锄e i l tm e m o d ， a r eb a s e do nt l l em a c r o s c o p i cna _ v i e r - s t o k e se q u a t i o n t h a t sb e c a u s en a v i e r - s t o k e se q u a t i o n i sm e 矗m d 锄e n t a le q u a t i o nt l l a t 目r o v 锄sj e l u i dn o w w e 昀n s f o n nn l ec o n t i n u u me q u a d o n i n t oas e r i e so fd i s c r e t ea 1g r e b r ae q u a t i o n s h o w e v e r b yu s i n gl a t t i c eb 0 1 t z m 锄m a b o d ，w ec a nl o o ki n t ot 1 1 es e c r e to fn u i dn o w i i lav e 拶i i l i c r o s c o p i cw a y a r e ral o n 哥t i m ee v o l u t i o no fl a t t i c e s c o l l i s i o na n ds 讹a l i l i i 培，w e c a l lh a v ea 血n ym a c r o s c o p i ci m a g eo fm en o w u s i n gs t a t i s t i c a lc a l c u l a t i o n 1 1 1t l l i sa r t i c l e ，l el l i s t o r yo fw 栅w a v er e s e a l 曲a n dl b m s t ：u d yi sp r e s e n t e d a n d 硪e r t h em a t h 锄a t i c a lc a l c u l a t i o n ，w eg 戍ag e l l e r a lv i e wo fm el 】mt l l e 0 巧吼dw a t e rw 打e6 e l d h lt l l ee ：l l d ，w eu s et l l el a t t i c eb 0 1 咖a n nl e m o dt os i n m l a t e 血nd 锄- b r e a l ( i n gw a v ea 1 1 d p 枷a ld 锄一b r e a l 【i n gw a v e c o m p a r i i 培w i t l ln l e 锄a l 如c a la i l do t l l e rm e m o d s s o l u t i o n ，w e c a n s e em a tt l l er g 巯l tc o m i n gf b ml a t t i c eb o l t 刁n 锄m e h o di sa c c l l r a t e k e yw o r d s ：l a t t i c eb o l 乜m a n n ；n u m e r i c a ls i m u l a t i o n ；w a t e rw a v e ； s h a n o ww a t e re q u a 6 明 一i i 中山大学硕士学位 论文原创性声明内容 本人郑重声明：所呈交的学位论文，是本人在导师的指导下， 独立进行研究工作所取得的成果。除文中已经注明引用的内容外， 本论文不包含任何其他个人或集体已经发表或撰写过的作品成 果。对本文的研究作出重要贡献的个人和集体，均已在文中以明 确方式标明。本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名：受广彳 日期：矿7 年6 月户日 学位论文使用授权声明 本人完全了解中山大学有关保留、使用学位论文的规定，即： 学校有权保留学位论文并向国家主管部门或其指定机构送交论文 的电子版和纸质版，有权将学位论文用于非赢利目的的少量复制 并允许论文进入学校图书馆、院系资料室被查阅，有权将学位论 文的内容编入有关数据库进行检索，可以采用复印、缩印或其他 方法保存学位论文。 学位论文作者签名：谈广哥导师签名：垮煮鼬 日期：旷7 年厶月f 。日日期：p 7 年6 月f 口日 中山大学硕士研究生学位论文 第一章绪论 1 1 选题背景 我们知道，物质是由分子构成的，物质的宏观表现物理量，在本质上均来自其内部 的微观分子运动的整体组合平均。为此，只要我们知道宏观流体内部的分子运动特征， 通过统计学的方法，我们即可得到流体的宏观特征如密度、速度等。这种计算方法在物 理学上被称为分子动力模拟方法( s ) 。然而由于分子数量的庞大，若然通过s 来进行计算，将会损耗很大数量级别的存储量和计算量。于是，一种简化型的m d s 应 运而生，这便是本文要研究的格子b o h 嬲方法。 格子b o l t 猢方法通过粗化m d s 模型，人工虚拟较大尺寸( 相对于分子) 的流 体微粒的存在。并在满足整体运动规律不变的约束下，简化流体微粒的演化规则，大大 减少了计算量和存储量，从而为这种简化版的如s 计算机数值模拟提供了推广的可能 【l 】 o 目前的数值模拟流体运动的方法有两种。其一是基于流体运动宏观方程的自顶而下 方法，其二则是基于微观离散模型的自下而上方法。传统的计算流体力学的数值方法大 多采用前者来进行数值模拟设计。这类方法从流体运动的n a v i 昏s t o k e s 方程组出发，采 用有限差分、有限体积、有限单元等离散方法对微分方程进行离散，得到代数方程组或 者常微分方程系统，进而使用标准的数值方法求解。 与常规的模拟方法不同的第二类方法，构造了人工的微观粒子模型，保持真实流体 的基本特征，从流体微观粒子运动特性出发，经过数理统计，得到宏观流体流动特性。 利用这类方法，可以简化流体运动，通过对不同的流体微观粒子设定相应的传播与碰撞 规则，经过统计，得到符合我们研究要求的宏观流体动力特征。本文所提到的格子 b o l t 蚴方法( l a t t i c eb o h 嬲m e t l l o d ，简称u m ) 就属于此类方法。 使用l b m ，我们将流体划分成众多不同流体微粒，同时将流体运动区域分割成不 同的网格。流体粒子在网格上运动，在设定的时间间隔，将保持原位置或者运动到它的 相邻格点上。经过上一时刻传播的粒子到达同一格点之后，将遵循流体运动的规则发生 碰撞，各粒子改变运动速度，成为下一时刻粒子传播的依据，在格点上的微观粒子运动 的统计，将成为我们所观察到的流体宏观特征如密度、速度、压力等。而粒子运动的时 利用格子博尔兹曼方法数值模拟水波传播 间间隔推进，为我们得到宏观流体运动的连续性。粒子碰撞规则的设定，源于连续流体 系统中物理量中质量、动量的守恒【l 】o 并且，l b m 数值模拟算法依据于格子b o h z m a n n 方程( l a t t i c eb o h 珊黜e q u a t i o n ， 简称l b e ) 而实现，我们可以从该模型的碰撞过程中的质量守恒、动量守恒出发，通过 数学推导，得出宏观流体运动n a v i * s t o k e s 方程组，实现微观模型与宏观方法方程的统 一。 格子b o h 础方法算法的微观粒子背景使得它可以提供了诸多分子动力学的优 点。从物理研究的角度，它可以获得清晰的物理图案，更好地处理流体与外界边界、流 体组分之间、流体界面之间较为复杂的作用。从计算的角度，它简单清晰的演化过程使 得它的程序编写天然地具备较为简洁的特点，同时它局部性的计算使得它可以提供并行 编程的便利。基于此，它在进行复杂边界流体、多种接触界面流体、大型流体并行计算 过程中，具有很成功的应用。在流体力学数值模拟的实践中，格子b o h 础方法也得 到许多领域科学家的关注，被广泛应用于多相流、化学反应扩散、粒子悬浮流、磁流体 力学等相关领域。【2 】 控制水波运动的浅水波是一种不可压的流体，同时受到重力等源项作用影响， n a v i * s t o k e s 方程会有一定的特殊化，这样我们的研究未知数将从密度转化到水波波 高。 。 水波的格子b o n 删方法数值模拟，是从微观粒子模型化的角度出发，对水波的 传播运动进行模拟，用此方法进行数学验证和实践计算，计算得到水波传播过程的动力 特征。 作为一种从极其微观化的m d s 模型到宏观方程计算模型的中介物，格子b o h 猢 方法架构起了一种联系微观和宏观的桥梁，通过这个方法的研究，我们可以更清晰地得 到流体的运动实质。【3 】 1 2 格子b o l t z m a n n 方法研究回顾 l b m 方法源自于2 0 世纪7 0 年代的流体计算l g a ( l a 钍i c eg a sa u t o m a t a ) 模型。 j e a nh 锄d y ，o l i v i e rd ep a z z i s ，y v e sp o 瞄踟在1 9 7 3 年和1 9 7 6 年发表的两篇论文创建了 第一个l g a 模型，并以三位作者的名字命名为 玎咿模型【4 】【5 】。 h p p 是以正方形格子为基础的最简单的二维l g a 。在该模型中提及，每个网格上 的流体粒子具有相同的质量，只能沿四个网格线的方向运动，粒子运动的速度刚好使它 在下一计算步长的时间运动到其相邻格点。并且每个格点上最多只能有一个粒子沿着某 条网格线运动，以达到遵循粒子运动的p 踟l i 不相容原则的目的。在h p p 模型的模拟碰 一2 一 中山大学硕士研究生学位论文 撞方式上，将遵循碰撞前后的质量和动量守恒，因此，只发生一种方式的有效碰撞即对 头碰撞，两个粒子沿相反方向到达同一个网格点而另外两个方向上没有粒子。这种碰撞 的结果将使得两个粒子分别旋转9 0 度沿另外两个方向运动。然而，由于对称性的不足， h p p 模型并不能推导出宏观动力学的n a v 衙s t o k 骼方程组。皿p p 模型的格子和离散速 度图如图2 1 所示。 ! ! ! ! 呻卜l i i 令一 iiiiii iiiiii 图2 1 玎叩模型的格子和离散速度示意图 f i 吕2 - 1n 聆l a n i c 嚣a n dd i s c r e t ev e l o c i t i 豁o f h p pm o d e l 在1 9 8 6 年和1 9 8 7 年的两篇论文中【6 l ，网，命名于论文三位作者的f h p 模型被引进。 f r i s c h ，h a 豁1 a c h e r 和p o i 豫a u 构造了基于六边形格子的l g a 模型，提出了比h h p 模型 具有较高对称性的f 肿模型。该模型的流体格子由正六边形构成，每个格点拥有六个 运动方向。每个方向上最多允许有一个运动粒子，该模型具有很强的对称性，能够在宏 观局限下得到n a v i * s t o k e s 方程。该模型的提出让世人开始瞩目l g a ，华盛顿邮报 用头条指出l g a 模型具有比以往的计算方法快一千到一百万倍速度的潜力。这三位作 者的工作为l g a 方法的发展奠定了理论基础【8 j 。 1 9 8 8 年，美国l o sa 1 锄s 国家重点实验室的m c n a m a r a 和z a 瑚n i 提出在l g a 中 直接使用布尔变量豫的统计平均量f ( 单粒子分布函数) 代替进行演化，即使用格子 b o l t z 贼m n 方程代替l g a 的演化方程进行计算【9 】。19 8 9 年，h i g u e r a 和j i i n e n e z 对该模 型进行了进一步的简化，提出线性化碰撞算子模型【1 0 】【1 1 】。该模型引入平衡态分布函数 声对碰撞算子做线性化处理，简化了计算，使得计算储存量大大下降。此后，c h e ns 、 c h e nh 、q i a ny 等人提出了更简单的模型b g k ( b h a t n a g a r - c 衲s s k 0 0 k ) 模型，又称 一3 一 利用格子博尔兹曼方法数值模拟水波传播 为单松弛模型。b g k 模型大大简化了模型的计算，并在一定条件下可以从该模型到处 n a ：v i * s t o k e s 方程【1 2 1 【1 3 】【1 4 】【悦。目前该类模型是格子b o h 舢方法中最主要也是应用 最广的模型。 格子b o h z 嘲n n 方法是一种与传统数值方法完全不同的方法。在该方法中，流体的 宏观力学特性完全由微观粒子的运动所决定，这也印证了物体运动的本贡一宏观现象皆 是大量微观现象的综合体现。流场是运动粒子的宏观效应，将流场划分为网格，让 b o l t 皿姗m 方程中的粒子分布函数沿网格线运动，在网格点上根据一定的定律相互碰撞。 并由分布函数的演化在宏观上得到流体的运动特性和规律，得到流场中的密度、速度、 压力等等各种力学函数。基于它的微观动力学特点，它天生具有其它常规数值方法( 如 有限单元法、有限差分法、有限体积法、谱方法等) 所没有的优点： 1 ) 该方法易以编程实现。由于微观粒子之间的碰撞规则简单，在l b m 计算的主体 算法编程简单而清晰，编程方便； 2 ) 该方法物理表达非常简洁。在整个运算过程中，只有一个微观量，即分布函数， 是未知而需要计算来确定的，单一的微观未知数，使得物理规律和数学表达变 得简洁； 3 ) 该方法适以进行并行编程。运算中某一时刻的分布函数值，完全由碰撞规则和 上一时刻的分布函数值所确定，这一特点使得l b m 程序很容易实现并行运算； 4 ) 该方法易以处理复杂边界。l b m 方法对边界条件的处理比较简单，使它可以处 理复杂几何边界如多孔介质流动的数值模拟计算； 一 5 ) 该方法处理复杂流体较为简单。它比较适宜处理如多相流、边界变动流体等复 杂流体的数值模拟计算。 由于以上的优点，格子b o l t 朋珈m 方法正日益成为一个充满希望的计算方法，在多 个不同领域中得到了成功的应用，尤其是在复杂边界和复杂流场的模拟中获得了巨大的 成功，比如复杂几何结构的湍流绕流、多空介质流、多相流、悬浮粒子流动、化学反应 扩散流动、渗流、磁流体力学、燃烧数值模拟、结晶过程模拟等等。并由纯理论试验研 究向工程实际应用迈进，形成国际计算流体力学的研究热点，并已经出现了专业商业软 件使用格子b o h 舢方法作为软件内核( 如美国e x a 公司的p 0 w c 棵1 0 w 系列) 。目 前，它的应用研究领域已经渗透到了生物流体、交通流、微尺度流、纳米流体等众多领 域，加速了向其他领域的渗透。 1 3 水波模拟与格子b o i t z m a n n 方法的结合研究回顾 在水波数值模拟的应用领域中，格子b 0 1 t 删方法的使用正处于快速发展阶段。 一4 一 中山大学硕士研究生学位论文 国内方面，在1 9 9 6 年和1 9 9 7 年的两篇论文中，吉林大学数学系的阎广武、胡守信、 施卫平等人，利用格子b o h 姗模型，构建了正六边形晶格，研究了波与其粘性实验， 并计算了一维浅水波模型( 水流过隆起物算例) ，得到了与该问题解析解接近的计算值 【1 6 1 7 lo 2 0 0 0 年和2 0 0 4 年，程永光和程冰等人，分别对两种格子b o h 删模型进行研究， 建立了统一的浅水方程l b 模型【1 3 】f 1 9 1 ，并对一维溃坝波进行了数值模拟计算，得到了良 好的计算结果。 同时，中国科学院大气研究所的刘峰和胡非在2 0 0 3 年利用格子b o h 础n n 方法数 值模拟了正压大气方程，并在2 0 0 4 年对该模型进行修改【2 0 1 ，应用于浅水动力学模拟， 准确模拟得到了二维溃坝流动的计算值，与理论解基本吻合。 而他们的研究，均受益于陈十一、郭照立、郑楚光等人对格子b o l t z 胧衄的引进和 推动。 国际方面，j g z l l o u 在2 0 0 1 年建立了格子b o h 础方法模拟浅水波方程中的弹 性碰撞方法来处理流体复杂边界条件，处理了比较复杂的曲面条件【2 1 1 。他还在2 0 0 2 年 利用格子b o h z 瑚n n 方法研究了浅水波方程，并数值模拟了经过隆起物流动、潮汐传播 流、圆柱绕流、明渠扩大流，均得到了良好的计算结果圈。2 0 0 7 年，他还利用格子 b o l t z 胧啪方法数值模拟了非连续流体，他所模拟的水面瞬间扰动流、二维溃坝、波浪 绕圆柱流动，均有良好的计算结果。1 2 3 1 2 0 0 7 年，g r i l i d o1 删n e s ，m o l l a m m e ds e a i d ，m a p u n d ik b 锄d a 等人发表的文章中， 利用格子b o l t 础方法模拟了浅水波流动【2 4 1 。并计算了直布罗陀海峡的实际流动问题， 证明了利用该方法计算水波运动并获取到主要流体参数的能力。 1 4 本文工作 本文的工作，收集和整理了格子博尔兹曼方法、水波数值模拟方面的研究资料，并 掌握了这两方面研究的数学推导原理。在对两者研究掌握熟练的基础上，采用比拟的方 法将格子b o h z 艘m 方法应用到水波的数值模拟计算中，成功模拟了全溃坝波、部分溃 坝波的传播问题。 。一本文的研究结果，若通过改进，可应用到实际水工设计和预防的控制管理上。目前 采用格子b o h z 胧m 方法在水波模拟领域的研究还处在起步阶段，理论和应用的研究尚 待完善，本文对该类研究的发展有一定程度的促进作用。 一5 一 利用格子博尔兹曼方法数值模拟水波传播 第二章计算模型原理 2 1b oit z m a n n 方法流体运动数学模型 2 1 1 分子动力学模型 分子动力学模型研究流体微观分子的空间位置和速度等物理量的时间演化，并利用 统计方法通过这些演化信息推导出宏观流动的物理特征。 众所周知，基于连续微分方程的数值模拟方法不考虑流体的微观分子结构，直接通 过对流体的宏观物理量如密度、速度等进行数值描述。但是在物理上，宏观流体是由流 体分子构成的，流体的宏观运动也是微观的流体分子热运动的统计结果。因此，如果能 够知道流体分子的运动情况，从理论上，是可以通过统计的方法得到流体的宏观描述。 这是分子动力学模型的基本思想。 在分子动力学模型中，通常认为分子的运动遵循经典运动方程如牛顿运动定律方 程： 查：v ，立：a ：里( 2 1 ) = v ，一= a = 一 l z 一上， 防戤聊 其中r ( ，) 、v ( f ) 和a ( f ) 分别是分子在时刻f 的空间位置、速度和加速度，f o ) 是作用 在分子上的外力，肌是分子质量。分子之间存在相互作用势矽( r ) ，这引起相互之间的作 用力： l ；，= v 痧( m ) ( 2 2 ) 其中弓，2 i i = ，是分子i 和分子歹之间的距离。再加上作用其他外力如重力、电磁力 等，作用在分子f 上的总的合力为： e = f ：! ，+ g ， ( 2 3 ) 其中，g ，为其他外力。目前应用最广泛的分子相互作用势为l e n n a r d j 0 m s 势。可 表达为： 中山大学硕士研究生学位论文 ( ，) = 4 占f 】6 【( ! - ) 2 1 】 ( 2 4 ) 厂，- 其中，仃为分子直径。 一旦确定了分子间的相互作用势，即可积分得到分子的运动速度和方位，但是，由 于流体中分子数量的庞大，用解析方法难以求解分子运动。目前该方法主要通过计算机 来实现，这边是分子动力学模拟( m o l e c u k 功m a “cs i n m l a t i o n ，简称m d s ) 方法。 加s 是基于物理上最基础的运动规律，所研究对象为物理存在的分子，所以，该 方法在原则上可以模拟任意流体系统。目前，分子动力学模拟方法在化学、生物学、物 理学和材料科学等领域均有良好应用。 但是，在数值模拟的过程中，作为一个流体系统，往往含有非常庞大的分子数目； 并且，如果使用s 方法，每一个演化的时间步都很小。这两个原因，造成泐s 在 模拟一个流体系统的时候，需要耗费巨大的存储量和计算量。即使模拟一个很短时间的 分子演化过程，也要花费很长的计算时间。不适合于计算方法的通用推广和实际操作。 这个背景下，格子气自动机模拟方法便以一种简化的s 方法出现。在这种模型 中，研究对象不再是实际存在的、数目庞大的流体分子个体，而是研究者自由构造的， 数量和大小可以自由控制的流体分子微粒。这样，模型对存储量的要求可以大幅减小。 另外，模型还严格控制了微观粒子的运动轨迹，它们被规定只能在格线上以相同的离散 时间步长运动，这样一模型对计算量的要求也得以改进。格子b o h 础方法便是以这 种方法为基础发展起来的流体力学模拟方法。 2 1 2b 0 1 t z 髓n n 方程 一 多体离散系统可以从不同层次来描述。从微观的层次，可以用h a i l l i h o n 方程来进行 描述。因此，对于由大量流体分子构成的流体系统，也可以用h a i i l i h o n 方程对所有流体 分子进行描述。但是，即使是一个具有很小体积的流体系统，它所包含的分子数量依然 很庞大，例如，o 摄氏度时，一个大气压下的1 立方厘米的空气里大约含有2 6 9 1 0 1 9 个 分子，所以用h 雅l i l t o n 方程来描述这样的系统几近不可能。 一种更高层次的描述方法是用分布函数厶( q 。，p l ，q ，p ，f ) 来描述由n 个粒子构 成的离散系统( q ，n 分别为第i 个粒子的广义坐标和动量) 网。显然，分布函数厶包 含了所有动力学的全部统计信息。厶( q 。，p l ，q ，p ，r ) 遵循l 幻u v i u e 方程： 煦一挚i 盟盟一盟盟l ：o ( 2 5 ) 钟匀la q ，a p ，卸，a q ，l 。 一7 一 利用格子博尔兹曼方法数值模拟水波传播 其中，巩是系统的h a i n i h o n 量。对厶在相空间中进行部分积分可以定义一个新的 分布函数： e ( q l ，p l ，q 。，p 。，f ) = 圪i _ ( q l ，p l ，q ，p ，r ) d q ，+ ld p ，+ l ，d q d p ( 2 6 ) 圪为归一化因子。可以证明，一个关于只( 1 s ) 的耦合微分方程组与l i o u v i c 方程等价。该方程组称为b b g k y ( b o g u l j u b o v _ b o m g r e e n - k i r k w o o d - y v o n ) 方程组。 b o h 舢方程就是在一些假设条件下从b b g k y 层级方程组中得到的一个结果。这些 假设条件有：1 只考虑二体碰撞，如稀薄气体系统；2 两个发生碰撞的粒子在碰撞前的 速度不相关( 这假设通常称为分子混沌假设) ；3 外力不影响局部碰撞的动力学行为。 b o h 姗方程描述了单个粒子分布函数厂( x ，c ，f ) 虻f l ( q 。，p l ，) 的演化过程： 罢+ 四。+ a 田。厂= q ( 力 ( 2 7 ) 其中，x ：q 。为粒子的空间坐标，c 卫为粒子的速度( 所为粒子的质量) ，a 为 m l n 外力f 引起的加速度，q ( 力为碰撞引起的变化。在这里，单粒子分布函数厂( x ，c ，f ) 的物 理意义非常明确，( x ，c ，f ) d x d c 表示时刻，位于空间【x ，x + d x 】中且速度在【c ，c + d c 】范围 内的流体粒子数。 得到了单粒子分布函数之后流体的宏观物理量如密度p 、速度u 和内能p 可以通 过如下方程统计确定： 一 p ( x ，r ) = 朋i 厂( x ，c ，f ) d c 胛( x ，) = 聊i 旷伍c ，f ) d c ( 2 8 ) 胪( x ，) = 等月c u 陬x ，c ，r ) d c 求解b 0 1 t 舢方程需要知道碰撞算子q ( 力的表达式。不管采用何种形式，q ( 力 都需要满足守恒条件( 质量、动量和能量守恒) ： 恢( c ( 力d c = o ( 2 9 ) 其中纯( c ) ( 七= 0 ，1 ，2 ，3 ，4 ) 为五个基本碰撞不变量。= 1 ，( 鲲，仍，伤) = c ，及伤= c 2 。 般的碰撞不变量矽( c ) 都可以表示成饩的线性组合，即有： 一8 一 中山大学硕士研究生学位论文 ( c ) = 么+ 曰匝+ q 2 ( 2 1 0 ) 2 1 3b 0 1 t z 吼i l i lh _ 定理和平衡态分布函数 b o h 皿姗n 方程描述了流体系统的微观演化过程。可以证明，经过充分长的时间后， 系统最终会达到一个平衡状态。在该平衡状态下，粒子分布函数具有一个特殊的形式， 称为平衡态分布函数。 1 8 7 2 年，l u d w i gb o h 黝n n 定义了h 一函数： 日( f ) = j ( x ，c ，f ) l i l 厂d c d x ( 2 1 1 ) 其中，( x ，c ，) 是满足b o n 砌a n n 方程的任一正值函数。可以证明，日( f ) 是非增函 数，数学表达为： 塑! o ( 2 1 2 ) 讲 此即b o l t z m a n nh 定理。 b o l t z n 瑚mh - 定理表明，经过充分长的时间后，系统最终会达到平衡状态，此时有 望娶= o 。这等价于厂( c ：) ( c ：) 一( c ：) ( c 。) ：o 或者l l l ( c ：) + 厂( c ：) ：l n 厂( c ：) + l n ( c 。) 。 讲 ( 这里，c 是碰撞后c 的速度。) 因此，如果尸是平衡态分布，则有l i l 厂哪是一个碰撞 不变量，所以可以将其表示为基本碰撞不变量的组合形式： l l l 叼= 嘞+ q 睦+ 锡c 2 ( 2 1 3 ) 一! ! ! = ! z 有尸= 2 妒，将其代入到( 2 1 ) 式，可得到 坍f c 一_ _ 1 2 严以南几百 ( 2 _ 1 4 ) 这就是m a x w e l l - b o l t 删分布，为b 0 h 黝常数，r 为温度。 2 1 4b g k 模型逼近 b o n z 衄咖方程描述了流体系统的粒子分布函数的演化过程。通过积分我们可以得 到分布函数得到流体的宏观密度、动量和能量。但是在实际计算中，由于粒子之间的大 量碰撞作用，使得该方程很难求解。 一9 一 利用格子博尔兹曼方法数值模拟水波传播 并且我们知道，流体的宏观性质是流体分子微观运动的平均结果，我们无需对微观 细节过于关注，即可得到我们想要知道的宏观特征量。因此，可以对碰撞算子做某种简 化，只需要使得简化算子满足原有碰撞算子的平均性质即可。一个简化碰撞算子q 应 当满足以下条件： a ) 碰撞不变量依( c ) ( 七= 0 ，1 ，2 ，3 ，4 ) 对q ( 力依然成立。即： i 依( c p d c = 0 = o ，l ，2 ，3 ，4 ) ( 2 1 5 ) b ) q ( 力能反映m a ) ( w e l l b o h 珊锄分布趋势。 基于此限制，b 1 1 a ：t n a g 孤，( 衲豁和k r o o k 引入了一个简化的碰撞模型，即b g k 模型 或单松弛模型。在该模型中，碰撞算子表示为： q ( 力= 一二【( x ，c ，) 一叼( x ，c ，f ) 】 ( 2 1 6 ) 其中，广( x ，c ，f ) 为局部m 删e 1 1 b o h 删平衡态分布函数，由计算式可以计算得 到，是松弛时间。b g k 模型碰撞算子的物理意义是，每次碰撞引起的分布函数变化 量与它偏离局部分布函数平衡态的量成正比，因此彩= 二也叫做碰撞频率。厂叼( x ，c ，f ) 随 时间和空间变化，所以称为局部平衡态分布函数。 2 1 5c h a p m a n e n s k o g 展开方法与宏观方程推导 描述流动现象的n a v 冲s t o k e s 方程是建立在连续介质假设基础上的，同样做为描述 流体运动规律方程，b o h 舢方程在一定条件下与n a 访* s t o k e s 方程应该是一致的。 从b o h 殂瑚彻方程到n 2 l v i 小s t o k e s 方程的推导，需要借助c 1 1 a p m a 玎e i l s k o g 于1 9 1 0 - 1 9 2 0 年间发展起来的多尺度展开方法【1 2 】。 在利用c 1 1 a p m a n e n 灿g 展开方法推导与b o n 蝴方程相应的宏观方程过程中， 为了简便需要，假设外力为0 。由于一个流体系统的密度、速度和能量等宏观物理量是 有限的，因此方程( 2 8 ) 的积分总是有意义的，从而有1 i m 仍矿= o 。利用分部积分有： 巨d c ：o ( 2 1 7 ) 。8 c 馐 善“一詈 ( 2 _ 1 8 ) 中山大学硕士研究生学位论文 b c 2 要d c ：一旦 ( 2 1 9 ) j2 阪 朋8 在b o h 删方程中，空间位置x ，粒子速度c 和时间，是相互独立的变量，因此关 于三者的导数可以交换次序。将基本碰撞不变量仍乘于b 0 h 舢方程式子两端并积分 可得到： 昙n ，d c + v ，口卜够d c + a 口协v 。厂d c = o ( 2 屹o ) 将方程( 2 7 ) 乘于所并取汪o 可得到连续方程： 誓+ 俨) = o ( 2 2 1 ) 国 ” 取i = 1 口3 可得动量方程： 挈+ 伊u ) = 册+ 伊 ( 2 _ 2 2 ) 其中， p = 聊f u o u o 厂d c ( 2 2 3 ) 取f = 4 可得能量方程为： 掣+ 胛p ) ：_ p v u + 河 ( 2 - 2 4 ) 讲 其中，j 是变形速度张量，为： j = 罟 2 u 。d c ( 2 嗡) 称为热通量，方程( 2 2 1 ) - ( 2 - 2 5 ) 一起构成了流体力学的运动方程组。 c h a p m a n e n s k o g 展开方法本质上是一种多尺度方法，其基本思想是将时空变量用 多层次时空尺度表示，而在各级尺度上物理量的量级一致。在c h a p m a i 卜e n s k o g 方法中， 假设系统虽没达到平衡态，但在每个局部小体积中，已离局部平衡态不远，因此可以假 设分布函数可以展开为如下形式： 厂= ”+ 手厂1 + 孝2 2 + ( 2 2 6 ) 利用格子博尔兹曼方法数值模拟水波传播 其中孝为m d n 数，定义为孝= 量。这里名是分子平均自由程，三为宏观长度。 由于流体宏观运动的时间尺度和空间尺度远远大于微观流体分子运动的时间和空间尺 度。所以引入两个宏观时间尺度 = 争，2 = 孝2 ，与及一个宏观空间尺度而= 豇。与之 对应，时间和空间导数可以表示为： 昙= 孝毒+ 孝2 毒，去= 孝毒c v = 尹- ， c 2 屯7 ， 对碰撞算子也做类似展开： q ( 力= q ( o ) + 髀1 厂1 + ( 2 2 8 ) 将( 2 2 6 ) - ( 2 - 2 8 ) 三式代入b o l t z m a n n 方程( 2 - 7 ) 中，由孝。阶的项可以得到： q ( 厂o ) = o ( 2 2 9 ) 该方程的解就是局部平衡态分布函数厂钾，即有o ) = 厂叼。 以o 代入( 2 2 3 ) 和( 2 2 5 ) 可得： 蟛= u 。五u 。，d c = 丁呜 ( 2 3 0 ) j 妒= 詈k 2 厂柳d c = o ( 2 _ 3 1 ) 将( 2 2 7 ) 、( 2 - 3 0 ) 和( 2 3 1 ) 代入到方程组( 2 2 1 ) 、( 2 - 2 2 ) 和( 2 2 4 ) 可得到 和而尺度上的 方程：。 竽+ v l 口 伊) ：o ( 2 3 2 ) 跣 一7 笔掣+ v l 胛u ) ：_ v l 田( o ) + 伊 ( 2 3 3 ) 巩 等+ v 。q 用d p ( o ) v 。u ( 2 - 3 4 ) 这与e u h 方程组是一致的。 对于r 尺度上的计算，经过推导最终可以得到： 中山大学硕士研究生学位论文 蚴= u 。 ，1 ) d c = _ 2 【踢一丢跗】 ( 2 3 5 ) j ! ) - 罟n 2 厂1 ) d c = 一七v l 口r ( 2 3 6 ) 础：丢( v l 口坳+ v l 户坳) ( 2 3 7 ) 其中和七分别是粘性系数和热传导系数； ：一兰掣) 抛 ( 2 _ 3 8 ) = 一一i - 兰一l l z j 石， 。 2 、聊7 七：一粤髑抛 ( 2 - 3 9 ) 4 、m 。 其中盯，白和岛是与分子碰撞相关的参数，并且有炙：筝。由此可以得到： 譬= 等= 主g ( 2 - 4 0 ) 七4 朋2 7 其中g ：警为定容比热。将p ( ，p ( 1 ) ，j ( ，j 代入到( 2 2 1 ) 、乜。2 2 ) 和( 2 2 4 ) z 力l 望：o 钙 旦逊：_ v 。口( 1 ) 一= 一v i r 、 掣：_ p ( 1 ) v l u _ v l ( 1 ) 础 1 ( 2 4 1 ) 、 ( 2 4 2 ) ( 2 4 3 ) 结合，l 和f 2 两个尺度上的宏观方程，可以得到精确到孝2 的整体宏观方程如下： 粤+ 伊) ：o 岔 ”7 掣+ 伽u ) ：一跏+ v 口) + v ( 譬：u ) + 伪 讲 j ( 2 4 4 ) ( 2 4 5 ) 利用格子博尔兹曼方法数值模拟水波传播 里墨呈堕+ v u p ) ：v 七v d 一册+ ( 2 4 6 ) 研 其中，为耗散函数。如此，便从b o h 舢方程导出了流体力学的基本方程，并 得到了粘性系数和热传导系数及定容比热。 2 2 格子b oit z m a n n 方法的原理和模型 2 2 1 格子b 0 1 t z 髓i l i l 方法原理 b o h z 翻m 方程用细观的角度描述了流体分子分布函数的变化规律，再加上可以通 过该方程推导得到流体的宏观运动方程，因此可以通过数值求解b o h z 胧咖方程的方法 来模拟流体的宏观运动。 格子b o h z 胧衄方法是b o l t 黝方程的一种特殊离散格式，该方法由一个规则的 网格构成，并使用该网格对b o l t z m a n n 方程进行离散化操作。离散时粒子的速度由格子 确定，碰撞算子较多地采用了b g k 算子。该方法的模型演化与l g a 相似，由碰撞和流 动两个子步构成。但与l g a 不同的是，在格子b o n 删方法的模型中，运动的不再 是粒子，而是粒子的分布函数。目前该方法正逐渐被更多学者研究，积极应用于流体研 究中。、 通过格子b o h 册姗n 方法，将所要求解的流体系统离散成规则的网格，同时立足于 粒子分布函数厂，将b o h 黝演化方程离散成如下的形式： z ( x + q ，r + 址) 一z ( x ，d = q 【( x ，明 ( 2 4 7 ) 该式子包含了离散网格上演化过程中的碰撞和流动量。在实际的计算中，对于碰撞 算子，基于b o n z m a n n 方程的b g k 算子模型被使用得较多，成为该方法研究中的最主 要算子模型。其形式如下： q ( 力= 一二( z 一铲)( 2 4 8 ) i 在格子b o h 舢方法中，该算子模型被称为l b g k 模型或者单松弛算子模型。它 大大提高了计算效率，只要选择恰当的平衡态分布函数，从该模型可以推导出正确的 n a 、，i e r s t o k e s 方程。粒子分布函数石和平衡态分布函数z 叼需要满足质量和动量守恒条 件如下： p = z = 铲 ( 2 4 9 ) 中山大学硕士研究生学位论文 胂= “= q 护 ( 2 5 0 ) 2 2 2 格子b 0 1 t z 髓n n 方法的离散速度模型 在各种b 0 h 蝴方法的各种模型中，q i 髓等人在1 9 9 2 年提出的隗q m ( n 维空 间，m 个离散速度) 系列是应用最广泛的离散模型。目前在二维流动问题的数值模拟中 所采用的格子b o h 舢方法，较多地使用了采用以正六边形为网格的d 2 q 7 模型和以 采用正方形为网格的d 2 q 9 模型。在三维流动问题的模拟中，较常使用的是d 3 q 1 9 和 d 3 q 2 7 模型。 对于d 2 q 7 模型( 速度结构如图2 - 2 所示) ，由正六边形构成格子，每个格点拥有 七个运动方向： q = ( c o s 器醐啦3 毒o ，6 ( 2 - 5 1 ) q 2 t ( c o s 谚，s i i l g 弦，q ：o l 彷3 ，矧，6 喵曲u 其中c = 缸垃，缸和址分别是网格常数和时间步长。 d 2 q 7 模型中格点上的平衡分布函数有f 2 5 j ： 叼| ，+ 竽+ 竽一学l 蚴 在这里，当扛。时，心= 三；当舀1 ；2 ，6 时，m = 击。该模型的声速g 有：乞2 = 手。 3 4 2 1 6 s 图2 - 2 d 2 q 7 离散速度模型 f 远2 - 2t l l ed i s c r e t ev e l o c 时i i l o d e lo f d 2 q 7 利用格子博尔兹曼方法数值模拟水波传播 d 2 q 7 模型的一些粒子碰撞规则如图2 3 所示： 。v 矽* * 函* 公v j ， * 一* * 南* * 一* 图2 - 3d 2 q 7 模型的粒子碰撞规则图 f i 吕2 3t h er m 鹤o f l a n i c 嚣砌1 i s i i n 龇m o d e l0 f d 2 q 7 对于d 2 q 9 模型( 速度结构如图2 _ 4 所示) ，它在一个二维正方形网格空间上构造 而成，在每个格点上，流体粒子拥有9 个速度方向： i ( o ，o ) ， f = o q = ( c o s 只，s i n b ) c ， 只= o 一1 ) 万2 f = 1 ，2 ，3 ，4 ( 2 5 3 ) 【互( c o s 瞑，s i i l 谚) c ，包= ( f 一5 ) 万2 + 万4 f = 5 ，6 ，7 ，8 其中c = 缸，止和缸分别是网格常数和时间步长。 同时在d 2 q 9 模型中，格点上的平衡分布函数有【2 5 】： 肛嵋p - + 竽+ 警一期 5 4 ， 中山大学硕士研究生学位论文 在这里，当扛。时，心= 罟；当i = l ，2 ，3 ，4 时，坼= 吉；当扛5 ，6 ，7 ，8 时，咩= 妻。 该模型的声速q 有： q 2 = 手。 s 1 8 图2 - 4d 2 q 9 离散速度模型 f i 吕2 - 4d i 鲫e t ev e l o c 时m o d e lo f d 2 q 9 2 2 3 边界条件处理 在实际的问题处理中，我们遇到的边界条件总是根据宏观物理量如密度、速度等的 条件，而在格子b o h 糊方法中，在每一个时刻所进行演化的却是不同类型的分布函 数。这需要我们设法根据流动情况的这些宏观条件，设计出相应的分布函数的边界条件。 边界条件的处理对于格子b o h 嬲方法的实现起到了重要的作用，对模型的精度和稳 定性都有很大影响。 目前格子b o h 删方法对于边界条件的处理，大体可以分为三类：启发式格式、 动