（流体力学专业论文）格子Boltzmann方法及其应用研究.pdf

南京航空航天人学硕十学似j 文 摘要 本文的研究内容主要包括两个部分：格子b o l t z m a n n 方法的理论 研究：格子b o l t z m a n n 方法在工程实际中的应用一一对二维驱动空腔 流、圆柱绕流及二维翼型绕流的模拟计算与结果分析。 理论研究部分系统地介绍了格子b o l t z m a n n 方法的起源、发展和 现状，对该方法的物理及理论背景知识进行了详细的描述，并给出了 标准格子b o l t z m a n n 方法的数学模型。介绍了格子b o l t z m i i n 方程与 宏观n a y i e r - s t o k e s 方程的关系，并具体介绍了格子b o l t z m a n n 方法 的二维粒子速度模型。为了提高格子b o l tz m a n n 方法在工程实际中的 应用，可以采用非均匀网格进行模拟计算。同时介绍了一些目前已有 的方法。在此基础上，提出了一种新的采用非均匀网格计算的格子 b o l t z m a n n 方法，并给出了详细的推导过程。然后，介绍了格子 b o l t z m a n n 方法的不可压模型及一些常用的边界处理方法，并给出了文 中所采用的方法。 最后，针对二维驱动空腔流、圆柱绕流及二维n a c a 0 0 1 2 翼型在不 同雷诺数的条件下，用新提出的方法进行了数值模拟。通过对模拟的 结果的详细分析与比较，显示该方法是行之有效的。 关键词格子b o l t z m a n n 方法，非均匀网格，泰勒级数展开，最小二 乘法，粒子速度模型，n s 方程 格子b o l t z n a n n 方法及其应心研究 a b s t r a c t t h ep r e s e n tt h e s i sc o v e r st w om a i np a r t s ：t h et h e o r e t i c a lr e s e a r c ho f l a t t i c eb o l t z m a n nm e t h o d ，a n da p p l i c a t i o no fl b mi ne n g i n e e r i n g - - f l o w i na2 - dl i dd r i v e n c a v i t y ，f l o w a r o u n dac i r c u l a r c y l i n d e r ，a n d f l o w s i m u l a t i o na r o u n da i r f o i la n da n a l y s i so ft h ef i n a lr e s u l t s a f t e r s y s t e m a t i c a l l yi n t r o d u c i n g t h e o r i g i n a t i o n ，d e v e l o p m e n t a n d c u r r e n ts i t u a t i o no fl b m ，t h e p h y s i c a l a n dt h e o r e t i c a l k n o w l e d g e i s d e s c r i b e da n dm a t h e m a t i c a lr e p r e s e n t a t i o no ft h es t a n d a r dl b ma n dt h e r e l a t i o n s h i p b e t w e e nl a t t i c eb o l t z m a n n e q u a t i o n a n d m a c r o s c o p i c n a v i e r - s t o k e se q u a t i o n si si n t r o d u c e d m o r e o v e r ，t h et w od i m e n s i o n a l p a r t i c l ev e l o c i t y m o d e l so fl b ma r ed e s c r i b e di nd e t a i l i no r d e rt o i m p r o v et h ea p p l i c a t i o no fl b m i ne n g i n e e r i n g ，n o n - u n i f o r mg r i dc o u l d b eu s e di nt h es i m u l a t i o n s o m er e c e n tm e t h o d sa r es h o w n b a s e do nt h e a b o v e d e s c r i p t i o n ，a n e wv e r s i o no fl b mu s i n gn o n u n i f o r m g r i d i s p r o p o s e d a n dt h e d e r i v i n gp r o c e s s i s p r e s e n t e d i nd e t a i l i n a d d i t i o n ， i n c o m p r e s s i b l e m o d e lf o rt h el b ma n ds o m eb o u n d a r yc o n d i t i o n sa r e d i s c u s s e d f i n a l l y ，b yu s i n gt h i sn e wm e t h o d ，f l o wi n a2 - dl i dd r i v e nc a v i t y ， f l o wa r o u n dac i r c u l a rc y l i n d e r ，a n df l o ws i m u l a t i o na r o u n dn a c a 0 012 a i r f o i lw i t hd i f f e r e n t r e y n o l d s n u m b e r sa r e p e r f o r m e d t h e c a r e f u l a n a l y s i sa n dc o m p a r i s o n a r ef i n i s h e d ，w h i c hs h o w st h a tt h en e wm e t h o di s e f f i c j e n ta n dr o b u s t k e yw o r d s l a t t i c eb o l t z m a n nm e t h o d ，n o n u n i f o r mg r i d ，t a y l o rs e r i e s e x p a n s i o n ，l e a s ts q u a r e s ，p a r t i c l ev e l o c i t y m o d e l ，n s e q u a t i o n s 承诺书 本人郑重声明：所呈交的学位论文，是本人在导师指导下， 独立进行研究工作所取得的成果。尽我所知，除文中已经注明 引用的内容外，本学位论文的研究成果不包含任何他人享有著 作权的内容。对本论文所涉及的研究工作做出贡献的其他个人 和集体，均已在文中以明确方式标明。 本人授权南京航空航天大学可以有权保留送交论文的复印 件，允许论文被查阅和借阅，可以将学位论文的全部或部分内容 编入有关数据库进行检索，可以采用影印、缩印或其他复制手 段保存论文。 ( 保密的学位论文在解密后适用本承诺书) 作者签名： 日期： 南京航空航天大学硕 = 学俺论文 1 1 问题背景 第一章绪论 自从二十世纪四十年代出现了第一台电子计算机以来，人们开始 进入了电子信息时代。随着高存储、高速度计算机的出现，人们所能 解决的问题也越来越广泛，同时所面临的问题也越来越复杂。在对流 动现象的研究中，以往人们大部分依靠的是解析方法，但所解决的问 题非常有限。而现实生活中所面临的流动问题往往十分复杂，如航空 航天器的亚跨超音速飞行、舰船的航行等等，依靠解析的方法来解决 这些复杂的流动现象是不可能的。到现今为止，人们对流体运动的研 究主要靠实验方法和数值计算方法。实验方法具有直观、结果基本可 靠的特点。但也存在较大的缺点：耗费大、周期长，并且结果受实验 条件的影响也较大，尤其是如今的航空航天飞行，速度离、飞行条件 复杂，用风洞来模拟困难是相当大的。而流体的运动可以由一组偏微 分方程描述。在大多数情况下，这些方程( 如n s 方程) 都是高度非 线性的，采用解析的求解方法是不实际也是不可行的。随着大型计算 机的出现，使人们可以借助于计算机用数值计算方法来解决复杂的流 动问题。因此，在二十世纪六十年代，用数值方法分析求解流动问题 的学科一一计算流体力学【1 - 3 1 ( c f d ) 逐渐发展起来。 伴随着电子计算机的飞速发展以及各种新颖算法的不断出现，c f d 已经形成了一门独立的学科，并且在航空航天、船舶、大型能源装置 ( 如核电站) 、新型交通工具、海洋工程、环境保护等众多工程技术部 门和领域都得到了广泛的应用，随着计算技术的发展、巨型计算机的 出现、计算方法的不断改进，计算流体力学在解决流动的理论和工程 实际问题中愈加显示出它的巨大作用。目前，计算流体力学已经成为 现代计算科学的最有力的推动力之一。 在计算流体力学中，传统的数值模拟方法可以分为两大类：( 1 ) 从 宏观角度出发，基于连续介质假设，采用数值计算方法，求解全位势 方程或e u l e r 方程或n - s 方程 ”；( 2 ) 从微观角度出发，采用分子动力 学【5 6 】的方法，对流动进行数值模拟。其中，格子b 0 1 t z m a n n 方法【7 】就 是典型的一种。格子b o l t z m a n n 方法( l a t t i c eb o l t z n l a n nm e t h o d ， 格子b o t g m a n n 方法及其席川研究 l b m ) 起源于格子气自动机州( l a t t i c eg a sa u t o m & t a l g a ) 。l g a 方 法是元胞自动机 9 1 ( c e l lu l a ra u t o m a t a ，c a ) 在流体力学中的具体应 用，是空间、时问和速度空间都离散的一个虚拟微观模型，与以连续 微分方程为基础的宏观计算流体力学方法有着本质的不同。l g a 的微观 特性使得它的边界条件非常容易实现，并且计算也很简单。因此，l g a 方法非常适于处理边界复杂的问题。更为重要的是，l g a 的计算具有局 部性和并行性，非常容易在并行机上实现，l g a 的出现不但为并行计算 提供了许多新思想，而且对并行计算机制造技术产生了重要的影响。 但是l g a 方法也有许多不足之处。例如，由于含有随机因素，l g a 的计算结果往往包含很大的统计噪声，l g a 的宏观方程也不是标准的流 体运动宏观方程。格子b o l t z m a n n 方法是为克服l g a 方法的一些内在 不足而发展起来的一种新方法。l b m 不但克服了l g a 的缺点，继承了 l g a 的主要优点，而且还有许多新的优点，如计算量小、计算效率高、 编程简单等。l b m 的产生与发展，不仅在计算流体力学领域中产生了深 远的影响，它所使用的处理方法和观点对其他许多学科也是富有启发 性的。 2 1 本文的工作 虽然格子b 0 1 t z m a n n 方法具有很多优点，而且它在工程实际的应 用中也取得不错的效果，例如多介质流1 t 0 - 13 1 、多孔介质渗透流动 ，”l 、 磁流体动力学 1 6 a 7 1 、传热问题 1 8 19 、m e m s 2 0 1 及化学反应2 1 。2 3 1 等等。但 是，它仍处于发展的阶段，它的理论基础还不是很完善，如系统的建 模方法和理论、模型精度的概念与分析方法、模型应用的定性和定量 分析等都有待研究，边界处理技术、非均匀网格等基本问题尚未得到 完全解决，这些问题在很大程度上影响了l b m 在工程实际中更大范围 内的应用。 由于格子b 0 1 t z m a n n 方法起源予l g a 方法，所以与l g a 一样，它 在计算中所使用的网格仍是均匀的。均匀网格的使用在很大程度上限 制了l b m 的适用范围，因为大量实际问题的数值模拟都需要非均匀网 格。所以，一些用于非均匀网格的l b m 相继出现 2 4 3 ”，这些方法既有 各自的优点，求解一些相应的问题并取得了比较理想的结果，但也存 在着一定的不足。本文的工作就是在标准l b m 的基础上，提出一种采 用非均匀网格计算的新方法。它与标准l b m 相比，不仅保持时间和空 南京航兰；三航天人学硕十学崽沦文 间上的精度不变而且具有实际计算量小、计算形式简单的特点。 本文的安排如下： 第二章介绍了格子b 0 1 t z m a n n 方法的起源一一l g a 方法，以及标准 l b m 的两种不同的导出过程。对从微观角度出发进行数值模拟的分子动 力学方法作出了一些具体描述。第三章介绍了格子b 0 1 t z m a n n 方程与 n - s 方程的关系以及在l b m 中使用的些粒子速度模型。 第四章详细介绍了我们提出的新方法一一基于泰勒级数展开和最 小二乘的格子b 0 1 t z m a n n 方法。从二维标准格子b o l t z m a n n 方程开始， 先后引入了泰勒级数展开方法和最小二乘法，使得l b m 能用于非均匀 网格。最后在确定速度模型的基础上，给出了用于实际计算的代数表 达式。 第五章介绍了格子b o l tz m a n n 方法的不可压模型和一些边界处理 方法，最后给出了在本文中所采用的模型和边界处理方法。第六章提 供了一些用新方法对实际问题进行的数值模拟。最后的结果表明，该 方法在工程实际中是可行的。 格子b o lt z m a n n 方法及其应心研究 第二章格子b o i t z m a n n 方法 2 1 格子b o i t z m a n n 方法的起源 流体是由流体分子构成的，流体的宏观运动是微观分子热运动的 平均结果。因此，如果能够用计算机模拟流体分子的运动，并能准确 地反映出分子之间的相互作用我们就可以模拟出流体的宏观运动。 这就是分子动力学模拟的基本思想。 分子动力学模拟是一种真正的微观模拟方法【3 2 。但是，这种方法 的计算代价是巨大的。模拟中，系统演化的步长必须很小，并且在每 一个时i q j 步，需要根据作用在每个分子上的外力和前一时刻的位置来 计算其新位置和新速度。对在前一时刻发生碰撞的任何分子，都要对 其进行判别并计算其新轨道。显然，这需要非常大的计算量。所以， 分子动力学模拟目前还只能用于简单的二维流动模拟 3 3 】，对三维复杂 流动进行模拟几乎是不可能的。 1 9 8 6 年，法国的f r i s c h 等人 3 4j 在一篇名为“l a t t ic e g a s a u t o m a c af o r t h en a y ie r sc o k e se q u a t i o n ”的文章中提出了格子气 自动机( l g a ) 的概念。可以说，它就是格子b o l t z m a n n 方法( l b m ) 的起源。l g a 是一种简化的分子动力学模型。在这种模型中，模拟的对 象不是数目庞大的流体分子个体，而是数目大大减少的流体粒子，即 微观充分大、宏观充分小的流体分子微团。这样，就大大降低了对计 算机存储量的要求。可以说，l g a 方法是一种离散的粒子动力学方法， 这是因为它采用的计算网格、时间和速度空间都是离散的。由于l g a 方法是元胞自动机( c a ) 在流体力学中的具体应用，因此它经常被称 为格子气元胞自动机。 一般的，一个格子气自动机包含一个均匀网格和网格节点上的粒 子( 即分子微团) 。用x 表示网格节点，每个节点上有m 个粒子运动方 向，可以用一个布尔变量一( x ，f f - l ，m ) 表示时刻，在节点x 上沿某个 速度方向是否有一个粒子在运动，即当有一个粒子运动时，n ，= 1 ，否 则，h = 0 。l c a 的控制方程为 n ，( x + f j ，t + 1 ) = ，( x ，f ) + q ，( ( x ，f ) ) ( i = o ，1 ，一，吁) ( 2 1 ) 这里，p ，是当地粒子速度；q ( 月( x ，) ) 是碰撞算子，它确定粒子接触后如 南京航空航天大学硕士学俺论文 何碰撞。设计l g a 的碰撞规则时要求碰撞算子满足质量和动量守恒 mm q ，( n ( x ，f ) ) = 0 ，e ，q ，( 出，f ) ) = 0 ( 2 2 ) ，1 0 ，_ 0 从初始状态开始，粒子在每个时间步上的变化可以分为连续的两 小步：( a ) 流动；在这一步中，每个粒子沿各自的速度方向运动到离 自己位置最近的节点，在这个过程中，粒子只能沿网格线运动：( b ) 碰 撞；当粒子到达某个节点时，它们之间相互作用并且根据一定的碰撞 规则改变各自的速度方向。在某一时刻某个节点上。最多只能有一个 粒子，并且该粒子的动量方向指向离自己最近的节点。 由f r is c h 等人和w o l f r a m 3 5 1 相互独立提出的l g a 模型( 又叫f h p ) 描述的是粒子在二维等边三角形网格上的运动变化。l g a 模型的碰撞算 子的复杂程度与粒子速度方向个数有关，一般地，使用肘个粒子速度 的碰撞算子的计算复杂度为2 “。图2 1 给出了一些碰撞过程： i n p u t s t 曩“o u t & “ i n p u t 鼬t e o u ts t a t e * t * ， ； 4 * 1* ，* - 。 * - * 1 图2 ，1f h p 模型中的碰撞规则 虽然格子气自动机具有算法简单、边界条件容易处理及完全并行 等优点，但它也有一些不足之处：( a ) 统计噪声：由于l g a 的碰撞算 子含有随机因素，因此不可避免的会有噪声的影响。虽然可以通过时 间平均或空间平均的方法降低噪声成分，但噪声的影响还是很大的， 并且需要更大的计算量。( b ) 碰撞算子的复杂性：l g a 碰撞算子的计算 复杂度与粒子速度方向个数成指数关系，不但增大了l g a 模型的设计 难度，而且不利于l g a 的应用。这一问题对三维情况尤为突出。( c ) 缺乏伽利略不变性；在l g a 宏观方程中，对流项前面有一个非单位的 格子b o l t z m a n n 方法及其应川研究 因子，虽然通过重新标定可以得到正确的n - s 方程 36 l ，但对多相流束 说这种方法是不可行的。 格子b o l tz n l g r l l 3 方法就是为了克服l g a 的这些不足而发展起来的。 正如首先提出用l b m 模拟流体运动的m c n a m a r a 和z a n e t t i t ”l 所言，如 果用单粒子分布函数( 真实变量) ，：= 代替粒子占位变量( 布尔 变量) n ，并且忽略粒子的单个运动特性和粒子间的相关性，最终， 我们就能够得到格子b o l t z m a n n 方程( l b e ) 。这里，表示总体平均。 所以，l b e 为 ，( x + e t a t ，f h a t ) = ，( x ，f ) + q ，杪( x ，f ) ) ，( f = o ，1 ，m )( 2 3 ) 这里，是沿f 方r n 的粒子速度分布函数；q ，驴“，) ) 仍是碰撞算子，它 表示发生碰撞后，的变化率。q ，仅仅与当地分布函数有关。此外，在 l b n 中，压力与密度存在这样的关系：po c p 。 在l b m 中，原始变量是粒子分布的平均值。由于和l g a 有着相同 的运动形式，l b m 保持了动力学方法中的局部性优点。而局部性是并行 性的基础。 与其它数值计算方法相比，格子b o l t z m g r lr l 方法具有这样的特点： 1 l b m 中的流动算子( 或移动过程) 在相空间( 或速度空间) 上 是线性的；这个特性是从分子动力学中继承而来的。在采用宏观量描 述的数值计算方法中，对流项是非线性的。通过多尺度展开技术，可 以把结合了碰撞算子的线性流动项转化为宏观的非线性对流项。在不 可压条件下，l b e 能转化成不可压n s 方程。 2 l b m 中的压力项是采用状态方程计算出来的；相反的，在对不 可压n - s 方程的直接数值模拟过程中，压力满足一个p o is s o n 方程。 其中，速度应变是它的源项。通常在解这样的压力方程时，会遇到一 定的困难。这就需要特殊的处理，例如松弛迭代。 3 l b m 中使用的粒子速度方向个数在相空间上是最少的：在包含 m a x w e l l b o l t z m a n n 平衡分布的传统动力学理论中，相空间是一个彻底 的函数空间。它的平均过程包含整个速度相空间中的所有信息。在l b m 中，由于仅使用了一两个或很少的几个速度方向，使得微观分布函数 和宏观量之间的转换非常简单，因此它的计算算法也是相当容易的。 4 l b m 中存在着最适于并行处理的局部相互关系模型； 5 边界条件容易设定。 南京航空皖大人学硕十学他论文 2 2 标准l b m 前面，我们已经导出了格子b o l t z m a n n 方程( 方程( 2 - 3 ) ) 。但是， 计算其中的碰撞算子q ，仍然是相当困难的。l b e 的碰撞算子最早是直 接根据l g a 的碰撞算子 ”】得到的，但后来对其进行了多种简化。其中 最重要的一个简化是h i g u e r a 提出来的【3 8 1 ，他通过假设分布函数近似 于当地平衡状态，使得碰撞算子达到了线性化 ，= j + ，“”( 2 - 4 ) 其中，”是非平衡态部分，3 是与时间步长同阶的一个小量。这样， q ，吨卜，掣矿+ o ( ) ：) ( 2 - 5 ) 这里采用了张量约定求和方法，即两个指标相同时表示对该指标求和。 注意到q 。o ) = 0 ，于是得到一个线性化的碰撞算子 q ，驴) = q ，一o )( 2 6 ) 其中，q ，：垦穹乒型称为碰撞矩阵。显然，如果使用肠个离散速度，q 的计算复杂度为m x m ，与l g a 的碰撞算子相比，计算复杂度大大降低 了。如果进一步要求q 。满足守恒条件，即 枷。= 0 ，腑，q ，= 0 ( 2 7 ) 则得到了强化碰撞算子 3 9 1 。这时碰撞矩阵是一个对称的循环矩阵，且 其元素仅仅与离散速度的夹角有关。 而最简单的线性化碰撞算子则是采用了一个松弛时间，这个松弛 时间是与使用了单松弛时间的当地平衡状态是接近的。它就是著名的 b h a t n a g a r - g r o s s - k r o o k ( b g k ) 碰撞算子【4 0 j q ( 旭，) ) ：幽二地：尘( 2 8 )q j ( 厂( x ，”= 丑旦尘丑型2 ( 2 8 ) 这里，”是工的平衡状态：f 单松弛时闻。因此，含b g k 模型的l b e ， 即标准的l b m ，可以写为 ，( 。+ 吒f ，f + f ) ：z ( x ，) + 芷堕生丛业，“：o ，1 ，m ) ( 2 9 ) 宏观量中的密度p 和动量密度用定义为分布函数，的粒子速度积率 格子b qlc z m a n n 方法及其应_ l 研究 m p = ： f ；0 刖= ：x s ；e ， ( 2 - 1o ) 在格子b g k 模型中，选用当地平衡分布函数是为了能够导出n s 宏观方程1 。采用格子b g k 模型还可以使计算效率提高并且使得传递 系数的选取具有灵活性。 2 3 用另一种方法导出l b e 虽然l b e 是从格子气自动机方法推导出来的，但是，通过在低马 赫数下的速度离散，l b e 也可以从连续b o l t z m s n n 方程导出 4 2 , 4 3 】。在 整个推导过程中，我们可以从含b g k 模型的b 0 1 t z m a n n 方程开始 瓦a f 垮夥= 一去杪一广) ( 2 - 1 1 ) 这里，孝是粒子速度：厂= s ( x ，孝，) 是在连续相空间( x 孝) 里的单粒子分布 函数；f ”是平衡分布函数( 或m a x w e l l - b o l t z m a n n 分布函数) ；五是松 弛时间。方程的右边项( r h s ) 模拟的是发生碰撞以后流体粘性在分 子级别上的影响。而宏观量( 如密度p 和动量密度胂) 则是，的流体 运动积分 p = 陟( x ，毒，f p 3 掌，p u = 阿( x ，毒，f 3 ( 2 1 2 ) 在不改变守恒率的条件下，速度空间掌可以被离散为一系列有限的点 k 。在离散的速度空间上，b 0 1 t z i l l a n i q 方程( 2 一1 1 ) 可以化为 鲁蝎w = 弓e 竹) ，( f ：0 1 ，m ) ( 2 - 1 3 ) 而离散速度的分布函数则有，( x ，) = f ( x ，孝，f ) 。 在方程( 2 - 1 3 ) 中，如果用一阶时间差分代替时间导数，并且用 一阶逆风格式离散传送项q - v , 而用顺风碰撞项q ( x p ，f ，r ) 代替 q ( x ，r ) ，则可以得到 ，( x ，f + f ) = ，( x ，f ) 一口( x ，f ) 一，( x q f ，f ) 一譬l ( ；_ f f ) 一；”x - e i a t , t ) ( 2 - 1 4 ) 这里，口= a t i e ，l x ：卢= r 五：a t 和a x 分别表示时间步长和网格步长a 取口= 1 ，即使用的网格是均匀的，卢：1 ，方程( 2 1 4 ) 就可以化为和 南京航空航天人学硕十学位论文 方程( 2 9 ) 一致的标准l b e 了。 从整个离散过程我们可以注意到。方程( 2 一1 4 ) 在空间和时间上 仅仅是一阶收敛到方程( 2 13 ) 。但是，正如s t e r l i n g 和c h e n 4 4 1 所展 示的，因为方程( 2 一1 4 ) 在它的空间离散中有一个拉格朗日特性，离 散误差中含有特殊的形式，它是包含在粘性项里的，所以最终的空间 和时间精度都是二阶的。 9 格子b 0 1 z i r a n n 方法及其应刖研究 第三章l b e 和n s 方程 3 1 连续8 0 l t z m a n n 方程及其流体动力性 我们口j 以从笛b g k 模型的惩续b o l tz m & n n 万程升始 盟a t 悖可= 丢杪一g ) ( 3 - 1 ) 这里，厂；，( x ，孝，f ) 是单粒子分布函数；毒是连续速度空间；f 是松弛时 间，它与粒子的碰撞过程有关：g 是当地m a x w e l l i 8 , i 3 平衡分布函数， 定义为 如，u ，班研p 唧 _ 譬 p ：， 这里，d 是速度空间孝的维数：p 、u 和口= h r m 分别是质量密度、宏 观速度及单位质量的标准化温度；、t 和m 分别是b o l t z m a n n 常量、 温度及分子质量。质量密度p 、速度1 1 和标准化温度口( 或内能密度) 是，或g 的流体运动积分 p = l 燃= l g a g ( 3 - 3 曲 加= j 彩善= j c g d # ( 3 - 3 b ) 詈矽= 睦倍一u ) 2 g = 珐倍一u ) 2 告( 3 - 3 c ) 在平衡状态下，厂= g ，连续b o l t z m a n n 方程( 3 - 1 ) 变为 d g ：0 ( 3 4 ) d t 这里， 昙：旦+ 孝v ( 3 5 ) d fa t 。 通过求解方程( 3 - 4 ) 的流体运动积分，可以很容易地导出e u l e r 方程 詈+ v ( 俨) = o ( 3 - 6 a ) 南京航空肮天人学硕士学位论文 掣冉( p u u + p 护) ：o 凹 堕丛+ v ( 旬+ 。+ p ) i l ：o o t 、。 这里，岛= d p o 2 、t = p u2 2 及p 分别是内能密度、 力；“2 = u u 。动量方程( 3 - 6 b ) 也可以写为 p 罢+ 刖v u 柚p 这里， ( 3 6 b ) ( 3 - 6 c ) 平动能密度及压 ( 3 7 ) p = p 【3 。5 j 是理想气体的状态方程。 通过使用c h a p m a n e n s k 。9 1 45 1 方法，厂的解可以展开h 厂= 厂f 0 + 厂1 + - ( 3 - 9 ) 这里，方程( 3 - 9 ) 中的户要比厂( o 小得多。由于，( o ) = g ，把方程( 3 - 9 ) 代入方程( 3 - 1 ) 可得 f ( i ) = _ r d 肼g ( 3 - 1 0 ) 因此， 厂+ ，( i ) - - - 鲁( 3 - 1 1 ) 由上述对f 的求解，方程( 3 一i ) 可以化为 害一r 粤d t = 一 仔，z ， d f 。 f 。 对方程( 3 - i 2 ) 求积分，就可以导出n - s 方程 警十v 恤) = o ( 3 - 1 3 a ) 掣十v 汹u + p 口) ：一v 1 7 ( 3 - 1 3 b ) o t 这里，n ( 1 ) 是一阶切应力张量 n 牡胁，i 1 伽一护( 等寺) m a ， 格子b o l t z m a n n 方法及其应用研究 很显然的是，在不可压条件下，即p z c o n s t a n t ，方程( 3 13 b ) 中的运 动粘性系数为 u = 胡 ( 3 15 ) 3 2 二维速度模型 在l b m 中，有着各种各样的粒子速度模型，如二维空间里的7 - b i t 模型、9 - b i t 模型 4 6 1 ；三维空间里的1 5 一b i t 模型、1 9 一b i t 模型、2 7 一b i t 模型 4 7 1 。由于在接下来的叙述中我们将用到二维速度模型，所以这里 给出了它们的具体描述。 3 2 17 - bi t 模型 7 - b i t 模型是在一个二维等三角形网格空间上构造的，所以它也被 称为d 2 q 7 模型。如图3 1 所示。 图3 1d 2 q 7 速度模型 它的粒子速度为 e 。= 1 ( ( c o 。, o 。t 以，。i n 铭k ，吼；仁一1 k 3 ，：! ：2 ，6 ( 3 1 6 ) 。口2 1 ( c o s 以，sn 铭，吼= 仁一1 k 3 ，口= 1 ，2 ， 。卜 这里，c = 蠡协，疵和国分别是网格常数和时i 到步长ad 2 q 7 模型的平衡 分布函数为 查重塾兰堕奎_ 人兰婴主堂垡堡奎 口= k p f ，+ 掣+ 掣等lc c c l 这里，当口= 0 时，w a = 1 2 ： 的声速c ，为：c 。2 = 臼= c 2 4 。 32 29 - hi t 模型 当口= 1 , 2 ，6 时，w a = 1 1 2 。 这里口为内能密度( 见方程 ( 3 1 7 ) 同时，该模型 ( 3 2 ) ) 。 9 一b i t 模型是在一个二维正方形网格空间上构造的，所以它也被称 为d 2 q 9 模型。如图3 2 所示。它的粒子速度为 f ( o ，o l 口= 0 p 。= ( e o s e ，s i n 口。k ，巳= q l k 2 ，g = 1 , 2 ，3 ，4( 3 1 8 ) l 2 ( c o s 眈，s i n 口a k ，以= g 一5 k 2 + 吖4 ，口= 5 , 6 ，7 ，8 d 2 0 7 模型的平衡分布函数为 口l + 掣+ 掣一蒡l 这里，当口= 0 时，k = 4 9 ；当口= 1 ，4 时，w a = t 9 ；当口= 5 ，8 时， w a = 1 3 6 。同时，该模型的声速为：f 。2 = 口= c2 3 。 图3 2d 2 0 9 速度模型 格子b o l t z m a n n 方法及其应用研究 3 3 格子b o it z m a n n 方程及其流体动力性 在上一章里，我们已经导出了含b g k 模型的l b e 以( x + 国，r + 西) = 厶( x ，r ) + ”( i ，r ) 一兀( x ，r ” 和它的流体运动积率 p = 正= 口 p u = e 。厶= 气口 口a r 3 2 0 ) ( 3 - 2 1 a ) ( 3 - 2 1 b ) 詈p 臼= ；圭o 。一u ) 2 丘= 莓圭o 。一u ) 2 尼9 ( 3 - 2 1 c ) 通过采用c h a p m a n e n s k o g 方法，可以导出l b e 的流体运动特性。 由下面的展开式【4 8 】 聃坞却训= 砉等 掣 ( 3 - 2 2 a ) 厶= ，“一” ( 3 2 2 b ) 这里，昙= ( 昙+ 气- v ) ；筘国。对上述方程取参数，的一阶精度得 臂) = 口 ( 3 2 3 a ) ，f 。0 ) = - rd f皿(口。l(3-23b) 由上述方程，我们可以导出l b e 的控制方程 d f 2 ：o( 3 2 4 a ) 型一1 ( 2 f - l 坶孥：一切) ( 3 - 2 4 b ) d t2 d r 2 f 。“ 在离散的动量空间( 由方程( 3 - 2 1 ) 定义) 上，求零阶方程( 3 - 2 4 a ) 的积率可以得到e u l e r 方程 南京航空航天大学硕士学位论文 鲁丹油) = o 掣冉+ 肌u ) = 。 堕坐+驴+尸)u：oot 、。7 ( 3 2 5 a ) ( 3 - 2 5 b ) ( 3 2 5 c ) 可以发现，方程( 3 - 2 5 c ) 与方程( 3 - 6 c ) 存在着一定的差距。由 流体的平动而引起的能量通量( 1 2 ) 2 u ( 即平动能) 消失了。只有在如 下的条件下，平动能才可能消失 y p o u 这里，= + 2 ) 2 ；“= 1 1 = 1 l 。因为0 = c ，2 ，所以，由上述不等式，我们 可以导出 m l c ，。其中，r 和分别是流动的特征时间和特征长度。 使用r 和对方程( 5 - 9 a ) 进行无量纲化可得到 土竺+ v “：0 ( 5 - 1 0 ) c 、i 墩 其中，p = p c ；，= t r ，v = l v ，u = u c ，。 从方程( 5 一l o ) 中可以发现，为了尽量降低人工压缩项寺筹的 影响，必须满足t l c 。 5 2 l b m 的边界处理 在使用l b 进行计算模拟的过程中，不断变化着的是微观的密度 分布函数或其它类型的分布函数，雨实际阊题的边界条件总是关于宏 观物理蠡的条件。因此，必须根据这些宏观边界条件确定出相应的关 于分布函数的微观边界条件，这就是l b m 的边界处理方法。 l b m 中的边界处理方法是从l g a 方法中继承而来的。例如，为了得 到无滑移速度边界条件。可以采用粒子分布函数反弹方法。之所以称 反弹方法，是指当一个节点a 上的密度分布函数尼沿p 。方向流动到边 界节点8 时，则下步该密度分布函数疋将沿e 。的反方向一。反弹回原 来的节点4 。如图5 1 所示。 反弹方法适用范围广，能够处理边界复杂的流动问题。但是这种 方法也有其缺点：该方法的精度仅为一阶，而l b m 在流场内部节点上 的精度为二阶。这样，反弹方法就降低了l b m 的整体精度垆驯；其次， 反弹方法不能处理运动边界条件和n u m a n n 边界条件以及其它的复杂边 界条件【5 9 j 。为此，人们提出了多种其它的边界处理方法。s k o r d o s l 6 0 i 格fb o l t z m a n n 方法及其应川研究 提出可以在边界节点上的平衡分布函数中增加速度梯度。n o b e l 等人【5 9 提出可以采用增加了一个压强约束以后的流体力学无滑移物面边界条 件。i n a m u r o 等人【6 i j 发现当采用反弹方法处理物面边界条件时，在物 面节点附近会产生一个滑移速度；他们提议可以用一个反向滑移速度 消除这样的影响。z ie g l e r 【62 1 注意到，如果把物面边界移到流体流动区 域的一半网格步长，即把无滑移条件放在网格节点之间，那么，反弹 方法就能给出二阶数值精度。 幽5 1反弹方法不恿幽 这些边界处理方法都不同程度地改善了l b m 的整体精度，但都有 一定的限制，不能处理一般的边界。1 9 9 6 年，c h e n 等人 6 3 l 提出了一种 简单的外推方法。该方法通过增加一个虚拟边界，把边界上的点看作 内点进行计算，而虚拟边界上的值则由内点和实边界上的值通过线性 外推确定。该方法能够保证l b m 的整体精度，同时能够处理不同的边 界条件。这里，我们就是采用的外推方法来处理边界条件的。下面给 出了外推方法的简单叙述。 含b g k 模型的一般l b e 为 ，( x + 岛a t , t + a t ) ：( x ，f ) 4 芷幽五蚴( 5 - 1 1 ) 可以发现，l b m 是一种特殊的差分格式。因此在设计l b m 的边界条 件处理方法时，可以借鉴有限差分方法中的边界处理方法。基于这一 思想，c h e n 等人提出了所谓的外推方法。以d 2 q 7 速度模型为例，如图 5 2 所示。 南京航空航天火学硕十学位论文 图5 2d 2 q 7 的外推方法不意图 设在流场边界d o a 之夕卜有一虚拟边界e f 。在执行每个时间步的流 动分步之前，使用一种线性外推方法，根据距边界d o a 一个网格步长 的流体层b c 和边界层d o a 上的密度分布函数值计算出虚拟边界e f 上 的密度分布函数值，即在发生流动分步之前对虚拟边界节点设置如下 的密度分布函数 f 2 = 2 r 一爿( 5 1 2 ) 其中，f 2 ，刀和z 分剐是虚拟边界、实际边界和第一层流体节点的 分布函数。 经过上述外推步之后，所以节点( 包括虚拟节点) 上的分布函数 值都是知道的，因此可以对所有的节点执行流动分布( 包括虚拟边界 节点) 。然后执行碰撞分布，但在实际边界上使用真正的边界条件( 速 度或压力) 计算平衡分布函数。 这种外推方法的优点是实酱适性较好，可以根据这种方法容易地 设计出一般边界条件的处理方法。同时，这种方法的计算量也比较小， 容易实现。 格子b o l ? , z n l a n p , 方法及其鹿j _ ；i 研究 第六章t l l b n i 的应用 格子b o l t z m a n n 方法是完全离散的粒子动力学方法，它不需要建 立和求解复杂的偏微分方程。具体到t l l b m 仅仅只需要在每个时间步 上计算一个简单的代数表达式： n + i 厶y ，+ 国) = k = 口蛐g 。 女一l 其中，系数a 。可以在确定了粒子速度模型和计算网格的情况下事先计 算好，并且在每个时间步上它是不变的。 格子b o l t z i i l & n n 方法的计算过程是一个循坏演化的过程，每个演 化过程包括碰撞和流动两个分步，整个过程如图6 1 所示。 图6 1 格子b o l t z i i l & n n 方法计算过程 由某一时刻( 时i n 步r ) 的密度分布函数五( x ，f ) ，可以计算出该时 刻的宏观密度p ( x ，) 和宏观速度u ( x ，t ) ；接着可以得到此时的平衡分布函 数臂( x ，f ) ：根据b g k 碰撞模型可以得到碰撞后的密度分布函数；通过 流动分步，就可以得到下一时刻( 即时间步f + 1 ) 的密度分布函数 南京航空航天人学硕十学位论文 厶 x ，+ 1 ) 。重复上诉过程，可以实现整个过程的不断演化。 此外，必须注意的是计算过程中边界的处理和初始值的设定。对 于边界的处理，在第五章中已经作了详细的介绍。对于初始值，一般 按被模拟的流体运动的物理意义来设定。对定常流动来说，由于最终 的稳定状态与初始状态无关，初始值的设定显得不是那么重要。而流 场是否达到稳定态可以用两个时间步之间的速度之差小于某个小量来 判定。 6 。1 二维驱动空腔的流动模拟 如图6 2 所示，在一个有开口的正方形区域内充满流体。在开口 边界外，有相同的流体以恒定的速度流过。这样，外面的流体就会驱 动正方形区域内的流体产生运动，该运动为有一定规则的漩涡运动， 如图6 3 所示。 图6 2空腔流动示意图图6 3空腔流动的流线 定义雷诺数为r e = u l u ，这里u 和分别是空腔外流体的流速和 空腔的宽度。为了便于比较，这里分别采用均匀网格和非均匀网格进 行了计算，如图6 4 和图6 j 所示。 采用非均匀网格是为了能在相同网格尺度的条件下更加准确地捕 捉边界处的流动状况，所以网格在边界处比较密而在远离边界处比较 稀。对于r e = 1 0 0 和4 0 0 ，网格尺度为6 5 6 5 ：r e = 1 0 0 0 ，网格尺度为 9 5 9 5 ：r e = 5 0 0 0 和t 0 0 0 0 ，网格尺度为1 4 5 x 1 4 5 。 格子b o i t z m a n n 方法及其应 _ j 研究 图6 4均匀网格 图6 5非均匀网格 初始时刻，整个流场的密度为p 。= 1 0 ；在空腔内部，速度为( o ，0 ) ； 在空腔顶部，只有沿z 方向的速度u = 0 1 5 。在每发生了一次流动和碰 撞后，空腔顶部的密度分布函数为厶= ；对于其余的三个固体边 界，采用的是外推边界处理方法。 在不同雷诺数下，附图1 给出了空腔流动的流线，附图2