（流体力学专业论文）格子Boltzmann方法的细观渗流数值模拟.pdf

摘要 摘要 随着科学技术的不断进步和工业生产的迫切需要，从孔隙尺度研究多孔介质 中的渗流已成为国际上渗流力学研究的热点，而细观渗流数值模拟无疑是细观渗 流重要手段和方法之一。 基于分子动理论的格子b o l t z m a n n 方法，是一种典型的微观方法。正是由于 格子b o l t z r n a n n 方法从微观角度出发，使得其在研究细观渗流方面有其它方法不 可比拟的优点。 本文的工作主要包括： 1 基本模型：介绍了几种常见的定常不可压模型和非定常不可压模型。给 出了这些模型的具体格式，并给出了由这些格式利用多尺度展开方法得 到的宏观方程。 2 边界处理：在非平衡外推方法的基础上，对其进行了改进，提出了一种 新的边界处理方法。该方法的基本思想是：边界处的未知宏观量的值采 用邻近的流场内部点二阶外推得到，并用该边界处理方法模拟了 p o i s e u i l l e 流、c o u e t t e 流、c a v i t y 流和柱体绕流等经典问题，数值结果 表明该方法为二阶精度。 3 网格技术：介绍了插值格子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 方法，在此基础上，引入了垂直平分网格划分( p e b i 网格) 方法，利用该方法对圆柱绕流和c a v i t y 流进行了网格划分。并对c a v i t y 流进行了流动模拟。 4 细观渗流：采用格子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 方法有限体积边界处理多点地 质统计三维重构 a b s t r a c t a b s t r a c t w i t ht h ed e v e l o p m e n to fs c i e n c ea n dt e c h n o l o g ya n du r g e n tn e e d so fi n d u s t r i a l p r o d u c t i o n ，r e s e a r c h i n gf l u i d sf l o wi np o r o u sm e d i a o nt h ep o r e s c a l eh a sb e c o m eo n e o fh o ts p o t si n t e r n a t i o n a l l yb yn o w f u r t h e r m o r e ，n u m e r i c a ls i m u l a t i o no ff l u i d sf l o w i np o r o u sm e d i ao nt h em i c r o s c a l ei so n eo fi m p o r t a n tm e t h o d sf o rs t u d y i n gf l u i d s f l o wi np o r o u sm e d i an o d o u b t l y b a s e do ng a sk i n e t i ct h e o r y , t h el a t t i c eb o l t z m a n nm e t h o di sat y p i c a ls c h e m eo n t h em i c r o s c a l e b e c a u s eo ft h eo r i g i nt h a ti sf r o mm i c r o s c a l e ，t h el a t t i c eb o l t z m a n n m e t h o dh a sm o r em e r i t st h a no t h e rm e t h o d si nf l u i d sf l o wi np o r o u sm e d i ao nt h e m i c r o s c a l e t h i sp a p e rm a i n l yi n c l u d e sa sf o l l o w s ： 1 l a t t i c eb o l t z m a n n m o d e l s s e v e r a lk i n d so fm o d e l sf o r s t e a d y i n c o m p r e s s i b l ef l o w sa n du n s t e a d yi n c o m p r e s s i b l ef l o w sa r ed e s c r i b e da n d e x p l i c i tf o r m so ft h e s em o d e l sa r eg i v e n i na d d i t i o n ，t h em a c r oe q u a t i o n s a r ed e r i v e df r o mt h e s em o d e l su s i n gm u l t i s c a l ee x p a n s i o n 2 b o u n d a r yc o n d i t i o n s b a s e do nt h en o n - e q u i l i b r i u me x t r a p o l a t i o nm e t h o d ，a n e wm e t h o di sp r o p o s e d t h em a i np r o c e s so ft h em e t h o da r et h a tu n k n o w n m a c r o s c o p i cv a r i a t i o n sc a nb ed e r i v e du s i n gt h en e a rn e i g h b o r i n gf l u i dn o d e t oe x t r a p o l a t ew i t ht h es e c o n d o r d e r t h en e wm e t h o di sa p p l i e dt os i m u l a t e s e v e r a lc l a s s i c a lp r o b l e m ss u c ha s ：p o i s e u i l l ef l o w s ，c o u e t t ef l o w s ，c a v i t y f l o w sa n ds oo n n u m e r i c a lr e s u l t si n d i c a t et h en e wm e t h o df o rb o u n d a r y c o n d i t i o ni ss e c o n d o r d e ra c c u r a c y 3 g r i dg e n e r a t i o ns c h e m e s i n t e r p o l a t i o ns u p p l e m e n t e dl a t t i c eb o l t z m a n n m e t h o d ，m y l o v s e r i e se x p a n s i o na n dl e a s ts q u a r e b a s e dl a t t i c eb o l t z m a n n m e t h o da n df i n i t ev o l u m el a t t i c eb o l t z m a n nm e t h o da r ei n t r o d u c e d t h e f i n i t ev o l u m el a t t i c eb o l t z m a r mm e t h o di sp r o v e da n d ，an e wm e t h o dt h a t p e b i 鲥dg e n e r a n t i o ni sb r o u g h ti n t ot h ef i n i t ev o l u m el a t t i c eb o l t z m a n n m e t h o di sp r o p o s e d 1 1 1 en e wm e t h o di su s e dt og e n e r a t eg r i d so ff l o w a r o u n dac i r c u l a rc y l i n d e ra n dc a v i t yf l o w f i n a l l y , c a v i t yf l o wi s s i m u l a t e dn u m e r i c a l l yb yt h en e wm e t h o d 4 f l u i d sf l o wi np o r o u sm e d i ao nm i c r o s c a l e t h el a t t i c eb o l t z m a n nm e t h o di s a p p l i e d t os i m u l a t ef l u i d sf l o wi np o r o u sm e d i ao np o r e - s c a l e f i r s t ， r e c o n s t r u c t i o no f p o r o u s m e d i ai s r e s e a r c h e da n d m u l t i p l e p o i n t g e o s t a t i s t i c si sa d o p t e dt or e c o n s t r u c tp o r o u sm e d i ai nt h r e ed i m e n s i o n s i i i a b s t r a c t s e c o n d ，t h ei a t t i c eb o l t z m a n nm e t h o di si m p l e m e n t e dt os i m u l a t et h e r e c o n s t r u c t i o nm o d e lb ym u l t i p l e - p o i n tg e o s t a t i s t i c sa n dt h r e e - d i m e n s i o n a l v e t o r sa n ds t r e a m l i n e sf o r t h ew h o l ed o m a i na r eg i v e n f i n a l l y , m u l t i g r o u p r e c o n s t r u c t i o nm o d e l so fd i f f e r e n tp o r o s i t ya r es i m u l a t e db yt h el a t t i c e b o l t z m a n nm e t h o d f u r t h e r m o r e ，p e r m e a b i l i t yo f p o r o u sm e d i a i sp r e d i c t e d k e yw o r d s ：p o r o u sm e d i a ，f l u i d sf l o wi np o r o u sm e d i ao nm i c r o s c a l e ，l a t t i c e b o l t z m a n nm e t h o d ，f i n i t ev o l u m em e t h o d ，b o u n d a r y c o n d i t i o n ，m u l t i p l e - p o i n t g e o s t a t i s t i c s ，t h r e e d i m e n s i o n a lr e c o n s t r u c t i o n i v 符号说明 分布函数 平衡态分布函数 碰撞算子 粒子速度 非平衡态分布函数 粒子速度大小 模型声速 密度 粒子模型系数 雷诺数 运动粘性系数 时间步长 松弛时间 压力 x 方向速度 y 方向速度 z 方向速度 符号说明 可 口 i ， z 矿q q 矿c q p q & d 缸 f p “ v w 论文原创性和授权使用声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工作 所取得的成果。除已特别加以标注和致谢的地方外，论文中不包含任 何他人已经发表或撰写过的研究成果。与我一同工作的同志对本研究 所做的贡献均已在论文中作了明确的说明。 本人授权中国科学技术大学拥有学位论文的部分使用权，即：学 校有权按有关规定向国家有关部门或机构送交论文的复印件和电子 版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进行检 索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 保密的学位论文在解密后也遵守此规定。 作者签名：l 蛐 沙。年6 月7 日 第1 章绪论 1 1 引言 第1 章绪论 流体通过多孔介质的流动称为渗流【l l ，这种流动现象在自然界中十分广泛， 其涉及到多个领域，如石油天然气开发、水利水电工程和岩土矿山等。所以，研 究渗流力学无论是基础研究还是工程应用，其意义都十分重大。 多孔介质内部结构非常复杂，当流体在多孔介质内流动时，也会呈现出非常 复杂的流动状态。由于多孔介质本身结构的复杂以及在其内流动状态复杂，使得 人们对渗流机理及孔隙内输运现象的认识远远落后于其他领域。 近年来，计算机的快速发展，使得大规模的科学工程计算成为可能，并逐渐 成为一个新兴的研究领域。随着计算机科学的进一步发展，大规模科学与工程计 算在科学研究和工程应用中的地位越来越重要，数值计算被公认为与理论、实验 相并行的第三种科学方法。 自从计算机问世以来，人们就进入了信息时代。随着计算机的不断发展，计 算与存储速度越来越快，计算机所能解决的问题也越来越广泛，但所面临的问题 也越来越复杂。在计算机产生之前，人们对流体力学的研究主要依赖于理论分析 和科学实验两种手段，但所能解决的实际问题非常有限，因为实际中的问题往往 非常复杂。对于这些复杂的实际流动问题，实验和理论分析手段往往无能为力。 计算机及计算技术的发展改变了这一状况，目前，数值计算方法已成为研究流体 力学的主要方法之一。流体的运动一般由一组偏微分方程组来描述的，大多数情 况下，这些方程都是非线性的( 如n a v i e r - s t o k e s 方程) ，采用解析求解几乎不可 能的，只有在非常简单的情况下才能得到解析解。正是这些问题亟待解决，促使 了计算机及数值模拟方法得到了迅速发展。二十世纪六十年代，用数值方法来 求解分析流动问题一计算流体力学( c o m p u t a t i o n a lf l u i dd y n a m i c s ，c f d ) 应运而 生，并经过几十年的发展已在各个领域得到了广泛的应用，如航空航天、大气海 洋工程、环境工程及生命科学等众多领埘2 1 。之后，随着巨型计算机的出现和计 算方法的改进，计算流体力学在理论研究和工程实际中显示了巨大的作用，并在 航空等各个领域起了强大的推动作用。在计算流体力学中，传统的数值方法一般 可以分为两种；一、从宏观角度出发，基于连续性假设，流体的控制方程为 n a v i e r - s t o k e s 或e u l e r 方程，这种方法可以用来解决宏观尺度的流动问题，在工 程中应用非常广泛；二、从微观角度出发，基于分子动力学方法，对流体在微尺 度流动进行模拟，这种方法在纳米等领域应用十分广泛。 第l 章绪论 这两种方法都有各自的优点，并在不同的领域里发挥着重要的作用。随着计 算流体力学的迅速发展，数值方法已经被很多人引入到了渗流力学中，但是，用 传统计算流体力学的数值方法( 如有限元方法、有限差分方法和有限体积方法等) 来模拟多孔介质内流体流动问题可能会遇到一些困难，如孔隙结构复杂，边界难 以描述，边界处理将十分困难。 因此，发展新的数值模拟方法将有力的推动渗流力学的发展。格子b o l t z m a n n 方法( l a t t i c eb o l t z m a n nm e t h o d ) 就是一种全新的数值模拟方法，其从分子动理 论出发，是一种典型的介观方法【3 1 。 a u t o m a t a ) ，后来被一些学者证实可由 它来源于格子气自动机( l a t t i c eg a s b o l t z m a n n 方程进行差分得到【4 1 。格子 b o l t z m a n n 方法是空间、时间和速度空间都离散的一种虚拟微观模型，与以连续 性方程为基础的传统计算流体力学有本质的不刚5 6 1 。它是基于微观模型和分子 动理论，因此不受连续性假设的限制，所以在微尺度流动领域方面( 连续性假设 不成立) 的模拟是传统计算流体力学所不能比拟的。尽管格子b o l t z m a n n 方法也 有很多不足之处，很多模型需要进一步的改进。但是其与传统计算流体力学相比 也有很多优点：边界条件容易实现、具有天然并行性和计算效率高等。 1 2 格子b o l t z m a n n 方法的发展历程 格子b o l t z m a n n 方法是源于8 0 年代提出的格子气自动机( l a t t i c eg a s a u t o m a t a ，l g a ) 方法，所以我们将介绍l g a 及l g a 发展到格子b o l t z m a n n 方 法的过程。 1 2 1 格子气自动机 一些复杂的流体系统，用n a v i e r - s t o k e s 方程进行描述和求解都十分困难。 例如：多孔介质中的流体流动，由于孔隙的大小，分布，走向都极为复杂，很难 对其进行精确描述，用n a v i e r - s t o k e s 方程求解该类问题很难处理这种复杂的边 界，还有一些复杂的流体系统根本就不存在宏观控制方程，更无法进行方程的数 值求解。 对于无法用n a v i e r - s t o k e s 方程进行描述的复杂流体系统，从微观或介观角 度构建恰当的模型是一种可行的方法。格子气自动机实际上是一种简化的介观流 体模型，基本思想如下：流体的宏观运动是流体分子做微观热运动的统计平均结 果，宏观行为对分子的运动细节不敏感。因此，可以构造一个尽可能简单的微观 或介观模型，使之能够遵守基本守恒定律。 在l g a 中，流场被看成由大量的粒子构成，每个粒子驻留在一个规则晶格 2 第1 章绪论 上，并按照一定的规则进行碰撞和迁移。粒子只能沿网格线运动，下个时刻粒子 正好运动到邻近的节点，又与其他方向到达该节点的粒子发生碰撞。有了以上规 定，可建立格子气自动机的演化方程： n i ( x + e i9 f + 1 ) = 惕( x ，) + q ( 刀) ，( f = o ，1 ，m ) ( 1 1 ) 式中：n ( x ，f ) 表示在t 时刻，x 节点上速度为e ，的粒子数；q ( 门) 是碰撞函数。 2 0 世纪6 0 年代，b r o a d w e l l 等人首次提出了离散速度模型，进一步对空间 进行了离散处理，但时间和空间仍然保持连续。7 0 年代，法国三位科学家h a r d y 、 p o m e a u 和p a z z i s 提出了第一个l g a 模型，并根据三位科学家的名字命名为h p p 模型。在该模型中，流体被离散成一系列粒子，空间和时间也被离散到二维正方 形格子上。由于正方形格子缺乏足够的对称性，h p p 模型中的应力张量不能满足 各向同性，这导致了h p p 模型不能反映n a v i e r - s t o k e s 方程的非线性项。 1 9 8 6 年，三位法国科学家f r i s c h 、h a s s l a c h e r 和p o m e a u 提出了具有更好对 称性的二维l g a 模型 7 1 ，即f h p 模型。该模型中，空间和时间是被离散到二维 正六边形格子上，这比h p p 模型具有更好的对称性，碰撞过程也更为丰富。同 年，法国的d h u m i e r e s 、l a l l e m a n d 和f r i s c h 提出四维面心立方( f c h c ) 模型及 在三维空间的投影。 在l g a 中，流体粒子在离散的格子节点上，沿着网格线进行迁移和碰撞， 在固体边界上，l g a 只需让进入边界的粒子做反弹处理，边界处理简单，可适 用于多孔介质等复杂边界条件的模拟。因为采用l 或0 来描述网格线上粒子的有 无( 即布尔运算) ，所以l g a 可以无条件稳定。但是这一方法也存在一些缺陷： 1 、由l g a 演化方程推导出来的动量方程不满足伽利略不变性； 2 、流体状态方程不仅与密度和温度有关，还依赖于流体的宏观速度； 3 、由于采用b o o l 运算，这给计算带来了数值噪音； 4 、碰撞算子极其复杂，给求解带来困难。 为了克服上述l g a 的不足，科学家们一直在寻找解决的方法，这就促使了 格子b o l t z m a n n 方法的诞生。 1 2 2 格子b o l t z m a n n 方法 1 9 8 8 年，m c n a m a r a 和z a n e t t i 把l g a 中的b o o l 运算变成实数运算，提出 用实数厂来表示系综平均后的局部粒子分布函数，这样就用b o l t z m a n n 方程代替 了l g a 的演化方程，这就是最早的格子b o l t z m a n n 方程，并用于了流体的数值 计算，开创了流体计算的一个新方向。 由于m c n a m a r a 和z a n e t t i 的格子b o l t z m a n n 模型( m z 模型) 使用实数型的 第l 章绪论 粒子分布函数代替布尔型的粒子数进行演化，所以m z 模型消除了l g a 系统的 大部分噪音。但是，这一模型采用的碰撞算子仍然具有指数复杂性，也没有克服 l g a 模型的其他缺点。 1 9 8 9 年，h i g u e r a 和j i m e n e z 对m z 模型进行了改进，提出了一种简化模型 ( h j 模型) ，将碰撞算子线性化，引入了平衡态分布函数尸。用一个碰撞矩阵 来代替碰撞算子，矩阵中个元素满足质量守恒和动量守恒，这样该模型不需要碰 撞规则表，各粒子间的碰撞细节也可以忽略，比多粒子碰撞模型容易构造。但当 粒子种类增加时，矩阵将变得很大，碰撞演算也将变得复杂。同年，h i g u e r a 等 提出了强化碰撞算子方法，增加模型的数值稳定性。 上述两类模型消除了统计噪声，克服了碰撞算子的指数复杂性，但由于使用 f e r m i d i r a c 平衡分布函数，仍然存在一些l g a 的缺点。 1 9 9 1 年1 9 9 2 年，q i a n 和c h e n 等人提出了单松模型【8 9 1 ，用一个时间松弛 系数来控制不同粒子靠近平衡态的快慢，进一步简化了碰撞算子，该模型称为格 子b g k 模型( l b g k ) 。该模型的碰撞算子非常简单，可表示如下： 1p1 q ( f ) = 一二iz z 钾l ( 1 2 ) f l 一 其中f 是松弛时间，f 钾是一个待定的平衡分布函数，l b g k 模型极大的提高了 计算效率，只要选择适当的平衡分布函数，就能从该模型推导出正确的 n a v i e r - s t o k e s 方程。l b g k 模型彻底摒弃了f e r m i d i r a c 平衡分布函数，采用了 m a x w e l l b o l t z m a n n 平衡分布函数。计算过程中平衡态分布函数的计算公式与流 场的维数和具体模型有关。选择合适的f 叼，可以从格子b o l t z m a n n 方程恢复到 n a v i e r - s t o k e s 方程具有g a l i l e i 不变性。 1 9 9 3 年1 9 9 8 年间，格子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 方程通过有限差分离散获得【4 ，1 0 , 1 1 】。这样，格子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 方法的稳定性进行了分析 1 2 1 。基本 模型方面，针对多相多组分系统，学者们提出了颜色模型【1 3 1 、伪势模型【1 4 ，1 5 1 和 自由能模型【1 6 1 等；在化学反应系统中，反应扩散和燃烧模型被提出【1 7 1 8 】；在非 规则网格方面，提出了局部插值【1 9 , 2 0 】、局部加密【2 1 , 2 2 】、有限体积格子b o l t z m a n n 方法【2 3 】等模型。 在此期间，格子b o l t z m a n n 方法在应用方面也取得了很大成就。很多新的边 界格式被提出了，如半步长反弹、外推格式等等【2 4 捌】。格子b o l t z m a r m 方法也被 用于了湍流的大涡模拟1 3 。微观渗流领域也是个很重要的应用，1 9 8 9 年，s u c c i 等人用格子b o l t z r n a n n 方法模拟了多孔介质流动，并验证了d a r e y 定律【3 2 1 。也有 4 第l 章绪论 些学者将l b g k 模型用于模拟渗流流动【3 3 1 。后来，一些学者提出了直接模拟宏 观渗流的模型。 1 9 9 9 年至今，格子b o l t z m a n n 方法受到越多越多的重视，众多学者加入到 这个研究领域，其发展十分迅速【3 4 1 。g u o 等人提出了相互作用力处理的新方法1 3 5 】； 边界处理方法也有了新的发展【3 6 4 7 】；网格技术方面，有限差分4 8 2 1 、有限体积【4 3 5 3 铷】、有限元【5 9 石2 】i 分块网格技术【6 3 】迅速发展；在多相多组分流模型方面，几个 不同的研究小组分别提出了大密度比多相多组分模型【6 4 柳】；渗流方面，用于大规 模工程实际问题和r e v 尺度的模型被提出【6 8 1 ；另外，多松弛模型也重新受到了 人们的重视【6 9 。7 4 1 。 1 3 细观渗流数值模拟研究进展 近年来，用格子b o l t z m a n n 方法研究多孔介质孔隙尺度的流动越来越受到学 者们的重视。 1 9 9 5 年，h e j i s 等人【7 5 】用格子b o l t z m a n n 方法研究了两类多孔介质的k o z e n y 常数，并对多孔介质的渗透率进行了预测。 1 9 9 6 年，a d r o v e r 等人【7 6 】通过数值方法对多孔介质的渗透率进行了预测，提 出了一个新模型。 1 9 9 7 年，k o c h 等人旧研究了r e 数与渗透率之间的关系。同年，m i c h e a l 等人【7 8 】研究了流体在非均质多孔介质内的流动情况。 k o p o n e n 等人【7 9 1 于1 9 9 8 年对三维随机生成的多孔介质进了数值模拟，计算 出了渗透率。 2 0 0 0 年，m a i e r 等人【8 0 】模拟了用立方体和球体随机堆积成的多孔介质内的流 动，研究了p e c l e t 数与耗散之间的关系。同年，s i n g h 等人【8 l 】研究了多孔介质空 间相关长度与渗透率之间的关系。 2 0 0 1 年，p a n 等人【8 2 】用格子b o l t z m a n n 方法和孔隙网络模型方法模拟了单相 流在球体构成的多孔介质内流动情况，用格子b o l t z m a n n 方法计算得到的渗透率 在小r e 数情况下和实验得到的相吻合。k i m 等人f 8 3 】对d a r c y 定律进行了非线性 修正。同年，k o c h 等人1 8 4 1 研究了多孑l 介质内流动的相互作用。 2 0 0 4 年，p a n 等人【8 5 1 用格子b o l t z m a n n 方法模拟了单相和多相流在多孔介质 内流动。 2 0 0 5 年，t a n g 等人【8 6 】用三维的格子b o l t z m a n n 模型模拟了气体在多孔介质 的孑l 隙尺度内的流动情况，主要研究了滑脱效应。同年，q u i s p e 等人【8 7 1 用格子 b o l t z m a n n 方法和m o n t ec a r l o 方法模拟了二维多孔介质内流动，研究了渗透率 第l 章绪论 与孔隙度之间的关系。 2 0 0 6 年，j e o n g 等人i 髓】用格子b o l t z m a n n 方法模拟了二维和三维多孔介质的 微流动，对d a r c y f o r c h h e i m e r 阻力进行了预测。 2 0 0 7 年，z e i s e r 等人【8 9 】用格子b o l t z m a n n 方法研究了牛顿流体在多孔介质孔 隙尺度的流动情况。 2 0 0 9 年，n a b o v a t i 等人【蚓分析了流体在多孔介质内的流动，给出了一个预 测渗透率的通用模型。 1 4 论文研究主要内容 2 1 世纪，渗流力学的某些领域将成为力学学科的前沿领域l l 】。格子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 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 a n n 方法网格划分，并引入了垂直平分网格( p e b i 网格) 。第五章为用 格子b o l t z m a n n 方法对细观渗流的研究，研究了多孔介质三维重构方法，并用格 子b o l t z m a n n 方法对重构的模型进行了流体流动的数值模拟，对不同孔隙度的多 孔介质的渗透率进行了预测。模拟了气水两相流在多孔介质内的分离过程。 6 第2 章格子b o l t z m a n n 方法原理 第2 章格子b o l t z m a n n 方法原理 格子b o l t z m a n n 方法发展至今已有2 0 多年，在此期间，理论和应用研究方 面都得到了长足的发展，并逐渐成为计算流体力学领域研究的热点之一。传统的 计算流体力学方法是基于宏观的微分控制方程，而格子b o l t z m a n n 方法是从微观 角度出发，基于微观模型和分子动理论，其基本思想是构造简化的粒子运动规律 来再现微观或细观运动的物理本质，其运动的宏观表现是满足宏观的运动控制方 程的。 2 1 格子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 模型导出n a v i e r - s t o k e s 方程组得到的是一个可压 缩方程。在推导过程中，假设流动的m a c h 数非常小，采用了多尺度展开技术。 我们知道，低m a t h 数的流动等价于不可压流动，因此格子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 方法在推导n a v i e r - s t o k e s 方程时使用了小速度展开，因此 格子b o l t z m a n n 方法只适用于小m a c h 数的流动。又因为小m a c h 数与不可压缩 等价，所以，在小m a c h 数的假设前提下，能够设计出模拟一般不可压 n a v i e r - s t o k e s 方程的格子b o l t z m a n n 模型。 我们以d 2 q 9 模型【8 】为例，如图1 所示，图中每个格点的速度离散成9 个方 向， 7 第2 章格子b o l t z m a n n 方法原理 6 3 2 5 。 e s ”e 、e l 7 夕e ? r 弓 7 4 8 图2 1d 2 q 9 模型 e o = 0 q = ( c o s ( i - 1 ) t r 2 ，s i n ( i - 1 ) j r 2 ) ， i = 1 ，2 ，3 ，4 e i = 压( c 。s ( i - 5 ) r c 2 + ，r 4 ，s i n ( f - 5 ) 万2 + 州4 ) ， i = 5 ，6 ，7 ，8 粒子的演化方程为： 碰撞： f ( x ，) = z ( x ，f ) 一l f j ( x , t ) 一z 叼( x ，f ) ( 2 1 ) 流动： z ( x + c e f a t ，t + a t ) = z + ( x ，) ( 2 2 ) 其中c = a x a t 为粒子速度，缸和缸分别为格子步长和时间步长，f 为无量 纲松弛时间，z ( x ，f ) 是在点( x ，f ) 沿e 1 方向的分布函数，z 叼( x ，) 是局部平衡分布 函数。平衡分布函数可写为如下形式： f q = p n + 3 e u + 兰c e ，u ，2 一三u 2 c 2 3 ， = 4 9 ；c o f = l 9 ，i = 1 ，2 ，3 ，4 ；c o j = 1 3 6 ，i = 5 ，6 ，7 ，8 宏观密度和宏观速度为： p = 彳，p u = e 。z ( 2 4 ) 对演化方程进行多尺度展开，可得到宏观流体运动的控制方程 a ，p + v ( 户u ) = 0 姒川冉( 删) = 一跏冉 p u ( v u + ( v u ) 7 ) 仁5 其中压力p = p ，g = l ；是声速，u = ( 2 r 一1 ) a x 2 6 z l t 是运动粘性。 8 第2 章格子b o l t z m a n n 方法原理 2 1 1z o u h o u 模型 1 9 9 5 年，z o u 等人提出了一种定常不司压n a v i e r - s t o k e s 方程的格子 b o l t z m a n n 模型，称之为z o u h o u 模型【9 1 1 。 对定常流动，方程( 2 5 ) 变为： v ( p u ) = 0 v ( p u u ) = 一跏丹v u + ( v u ) r ) 仁6 为了方便，我们假定密度为常数p = p o ，定常不可压n a v i e r - s t o k e s 方程为： u v u - - - 、t p + 胛2 u ( 2 7 ) z o u 等人选取了如下的平衡态分布函数： 肛q 卜警+ 曙+ 矧 亿8 ， 并定义流动的宏观密度和宏观速度： p = z ，p u = e ，z ( 2 9 ) 则利用多尺度展开法可得到下面的宏观方程： 望+ v u ：o o t ( 2 1 0 ) 罢冉( u u ) = 一印丹u + ( v 例 其中压力p 和粘性系数u 的定义和前面相同，在定常情况下，方程组( 2 1 0 ) 和标 准的常密度定常不可压n a v i e r - s t o k e s 方程( 2 6 ) 完全一致，消除了可压缩性。 2 1 2 l i n f a n g 模型 z o u - h o u 模型在计算压力p 时假定密度是变化的，得到的是常密度定常不可 压n a v i e r - s t o k e s 方程。这一矛盾被我国复旦大学学者l i n 等人注意到了，他们于 1 9 9 6 年提出了一种新的定常不可压的格子b o l t z m a n n 模型，称为l i n f a n g 模型 【9 2 】 o 该模型假设流体的密度为常数p = l ，平衡态分布函数与压力无关，如下： 肛q 卜等+ 譬+ 专i 亿m 9 第2 章格子b o l t z m a n n 方法原理 其中d 是一个与压力无关的参数，定义如下： d = z ( 2 1 2 ) i 流动的宏观速度为： p u = e ，z ( 2 1 3 ) i 利用多尺度展开方法得到的宏观方程与( 2 1 0 ) 式相同，压力为p = d 。定常 情况下与常密度定常不可压n a v i e r - s t o k e s 方程组一致，压力的计算不再与密度 有关。 2 1 3c h e n o h a s h i 模型 以上两个模型都是解定常不可压n a v i e r - s t o k e s 方程的，不能用来解非定常 不可压流动问题。为了用格子b o l t z m a n n 方法求解非定常不可压流动问题，c h e n 和o h a s h i 于1 9 9 7 年对z o u h o u 模型进行了改进，提出了一种新的不可压格子 b o l t z m a n n 模型，称之为c h e n o h a s h i 模型【9 3 】。 z o u h o u 模型利用分布函数( 2 8 ) 和宏观物理量定义( 2 9 ) ，再利用多尺度展开 技术可得到宏观方程( 2 1 0 ) ，但是仍是人工压缩方法。c h e n 等人提出了一种新方 法，将宏观速度u 看做一个临时速度，真实速度w 对u 进行修正得到： w = 1 1 + r f ( 2 1 4 ) 其中f 为外力项，假设为某个势函数矽的梯度： f = v ( 2 1 5 ) 对上面两式求散度得： v 2 ：一三v u( 2 1 6 ) 这里假定了真实速度w 满足不可压方程。在c h e n o h a s h i 模型中，平衡态分 布函数的形式与z o u h o u 模型的分布函数相似，只是用w 代替了1 1 ，即 肛q 卜等+ 簪+ 蓦i 亿聊 宏观密度和速度定义如下： p = zw = e ，z + 刃 ( 2 1 8 ) 因为p = e qw = e ，z q ，所以 z 叼= z ，q 护q z ( 2 1 9 ) 从上式可以看出，碰撞算子可以保持局部质量守恒，但不能保持局部动量守 1 0 第2 章格子b o l t z m a n n 方法原理 恒。采用多尺度展开技术，从c h e n o h a s h i 模型推导出宏观方程为： v w = 0 w v w - 唧+ 闪2 w ( 2 2 0 ) 其中压力p = c ；p 一矽，粘性系数如前所述。 2 1 4h e l u o 模型及d 2 g 9 模型 在c h e n o h a s h i 模型提出的同时，1 9 9 7 年h e 和l u o 也提出了一个不可压格 子b o l t z m a n n 模型，称之为h e l u o 模型【9 4 】。其基本思想是在平衡态分布函数中 消除由密度变化引起的高阶m a c h 数项。在该模型中，独立的变量是压力p 而不 是密度p ，这与不可压n a v i e r - s t o k e s 方程是一致的。 低m a t h 数的流动或不可压流体，流体密度近似为常数风，于是密度p 可分 解为p = p o + 印，其中印为密度的波动。在研究流体流动时，如果考虑为等温 过程，则由伯努利方程可知：p + p u 2 2 = p c ，p c 为驻点压力，有 8 p = p c p = p u 2 2 。绝热过程中，音速c 。可定义为： q = ( 2 2 1 ) 式中y 为比热容比，因为r a p = e 印，于是有： 延p 兰2 22 2 1 肌v a “1 2 ig 卜7 因此，平衡态分布函数为 z 叼= 哆p + p o s t ( u ) + o ( 胁3 ) ( 2 2 3 ) 式中，t ( u ) 为平衡态分布函数中与速度有关的部分： 咖h 睁簪+ 矧 亿2 4 ， 在平衡态分布函数( 2 2 3 ) 中略去m a 的高阶项，则可得到一个新的平衡态分布函 数：