（理论物理专业论文）格子玻尔兹曼方法在表面活性剂体系中的应用.pdf

摘要 本文介绍了格子玻尔兹曼方法( l b m ) 及从格子玻尔兹曼方程出发到n a v i e r s t o k e s 方程推导过程，并应用l b m 模拟了多元不相容体系的相态结构，重点研究了表面活性剂 的加入对体系形态结构及其演化过程的影响，得到如下的结果： 1 针对多元不相容体系，建立了剪切流动下l b m 模型、完成了l b m 模拟程序的编 写、调试和验证。 2 模拟研究了剪切速率、体系粘度等对二元、三元不相容体系的相态结构的影响。 研究结果表明随着剪切速率的提高体系的相态结构出现了明显的取向性，并最 终可形成层状结构。 3 表面活性剂的加入明显改变了不相容体系的相分离动力学过程。在油、水、表 面活性剂的体系中，得到了多种相态结构及其对组成、剪切速率的关系。 关键词：格子玻尔兹曼方法；不相容混合物；表面活性剂；相分离；剪切流动 a b s t r a c t i i ln l i s t i l e s i s ，l a t t i c eb o l t z m 锄m e t l l o d ( l 】3 m ) a 1 1 dt l l ed e r i v a t i o nf 0 皿l a t t i c e b o l t 珊a 1 1 nt 0t h ena _ v i e r - s t o k e se q u a t i o n sa r ef i r s t l yi n 仃o d u c e d ，a i l dt h e nt h el b mi su s e dt 0 s t u d ym ep l l a s es t n j c t u r e sf o ri i l c o m p a t i b l em u l t i c o m p o n e n ts y s t e m t h es t u d yi sf o c u s e do n t h ee 铱；c to ft l l ea d m t i o no ft h es u r f a c t a n to nm ep h a s es t 九l c t i 】r eo fo i l i a t e rs y s t e m t h em a i n r e s u l t so b t a i n e da r e 嬲f o l l o w i n g ： 1 t h el b mn 1 0 d e la i l di t sc o n l p u t e rs i m u l a t i o np r o g r a mh a v e b e e nb u i l ta n de d i t e d t h e r e s u l t si 1 1 d i c a t et h a t t l l e y a r e a p p l i c a b l e t 0t l l e s t u d y f o rm e i n c o m p a t i b l e m u l t i c o m p o n e n ts y s t e mu 1 1 d e rs h e a rn o w 2 t h ee 日e c to fs h e a rr a t ea i l dt h e 、，i s c o s i 够o nm ep h a s es t r u c t u r ew 嬲s i m u l a t e df o r t h e b i i l a 叫觚dt e l m u yi 1 1 c o m p a t i b l es y s t e mu n d e rs h e a rn o w t h es i m u l a t i o nr e s u l t ss h o w m a tm ep h 弱es t l l l c t u r ec a nb eo b v i o u s l yo r i e n t a t e d ，a n df i n a l l y1 a m e l l a rs 廿u c t u r ec a nb e f 0 玎n e dw i t l lt l l ei n c r e a s eo fs h e a rr a t e 3 t h es i l n u l a t i o nr e s u l t s 证d i c a t et 胁t h ea d d i t i o no f 也es u m c t a l l ts i 鲥矗c a l l t l yc h a n g e d t 1 1 ep h a s es t n l c t i l r e m a n yp h 私es 仃l i c 恤e s 锄d 也e i fr e l a t i o 璐w 岫m ec o i n p o s i t i o na i l d s h e a rr a t eh a v e b e e no b t a i n e di nt h i st b e s i s k e yw o r d s ：协t i c eb o l t z m 猢m e 廿1 0 d ；i i l c o m p a t i b l em u l t i - c o m p o n e n tm 议t u r e s ；s u r f a c t a n t ； p h a s es 印a r a t i o n ；s h e a rn o w i l 独创性声明 本人郑重声明：所提交的学位论文是本人在导师指导下独立进行研究工作所 取得的成果。据我所知，除了特别加以标注和致谢的地方外，论文中不包含其他 人已经发表或撰写过的研究成果。对本人的研究做出重要贡献的个人和集体，均 已在文中作了明确的说明。本声明的法律结果由本人承担。 学位论文作者签名：凇 日期：湖g 6 学位论文使用授权书 本学位论文作者完全了解东北师范大学有关保留、使用学位论文的规定，即： 东北师范大学有权保留并向国家有关部门或机构送交学位论文的复印件和电子 版，允许论文被查阅和借阅。本人授权东北师范大学可以采用影印、缩印或其它 复制手段保存、汇编本学位论文。同意将本学位论文收录到中国优秀博硕士学 位论文全文数据库( 中国学术期刊( 光盘版) 电子杂志社) 、中国学位论文全文 数据库( 中国科学技术信息研究所) 等数据库中，并以电子出版物形式出版发行 和提供信息服务。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：燃 日 期：塑星161 e 学位论文作者毕业后去向： 工作单位： 通讯地址： 誓导教师签翥! 毪 日期：二金一力- ， 电话：( 丝盐暨矽9 邮编：弓，也主 7 东北师范大学硕士论文 引言 油、水、表面活性剂的混合物体系可以展示出丰富的微观结构，而当其处于流场中 时则会呈现出更为复杂的形态，对其相行为、自由能变化等特性的研究，可以为生产加 工过程以及材料的应用提供依据。 目前，研究混合溶液相分离问题的方法有实验、理论分析与计算机模拟三种。实验 方法能发现新现象和新原理，它的结果可以检验其它方法的正确性，但由于受到实验条 件和实验手段的限制，对于很多体系尚不能给出完整的描述。理论分析方法虽能获得精 确的解析解，但这样的体系十分有限。计算机模拟方法能够有效的弥补实验和理论分析 方法的不足，在解决复杂溶液问题方面表现出了很大的优势。 随着计算机硬件和软件的不断发展，许多计算机模拟方法已经用来研究流体动力 学。格子气自动机( l g a ) 理论是解决复杂流体问题的有力工具。h a r d y 、p 伽e a u 和p a z z i s 于2 0 世纪7 0 年代建立的所谓h p p 格子气模型实际上是第一个完全的格子气自动机模型。 这个模型是将平面流场划分成正方形网格，每个节点上的粒子只能向四个方向之一运 动，且只有两个对头碰撞才是有效的。这个模型过于简单，不能充分描述流体的特征。 1 9 8 6 年，由法国u f r i s h ，y p o m e a u 以及美国的b h a s s l a c h e r 提出了一个对称度 更高的正六边形格子模型，即f h p 模型。他们用此模型成功地模拟了一些典型的流体力 学问题，并证明了该模型的宏观行为符合标准的n a v i e r s t o k e s 方程。格子玻尔兹曼方 法( l b m ) 则是在格子气自动机的基础上发展起来的新的计算方法，它克服了格子气自动 机模型在描述更广泛的物理系统时不能进行灵活的参数调整等许多缺点，而且在数值上 比布尔动力学有效得多，已被广泛应用到两相流、磁流体动力学、孔隙介质中的渗流和 热流体动力学等各种流体的流动问题中。 格子玻尔兹曼方法作为介观模拟方法的一种，虽然已经成功地对多种复杂流体进行 了模拟计算，但在表面活性剂方面的应用仍不多见。对这一问题的解决，可以对多种模 拟方法的贯通做有力的补充，而且其本身也将具有广阔的应用前景。 东北师范大学硕士论文 第一章l a t t i c eb o l t z m a n n 理论基础 在本章节中，首先简单介绍了在l b 方法中涉及到的有关非平衡态统计物理的知识以 及流体力学方程组。这些内容会对我们后面的理论分析及计算模拟有很大帮助。接下来 介绍如何推导l b g k 方程，这个方程式是模拟过程的出发点。最后列举在文献中给出的应 用实例。 1 1l g a 和l b m 方法简介 从物理学和流体力学的发展来看，描述流体流动主要在如下两个层次上：第一，连 续介质层次上；第二，运动论层次的描述。后者是从非平衡统计力学的观点出发，想象 流体是由大量的微观粒子所组成。这些粒子运动遵守力学定律，同时服从统计规律。这 是一个微观动力学模型。运用统计物理的方法来讨论流体的宏观性质。格子气方法就是 以这个层次的描述为背景的。流体运动论层次的描述所遵循的基本原则就是物理守恒定 律，它要求所有微观过程保持质量、动量和能量守恒。l g a 方法的理论分析手段就是从 此出发，从微观动力学方程出发，推导出连续介质描述的宏观流体力学方程，其微观动 力学是基于布尔运算的。推导中一个重要的假设就是分子混沌条件，即相互碰撞的分子 间统计独立。从另一个观点看，他们的相互作用是局域性的，作用之前和作用之后都彼 此无关。借助这一条件得到l a t t i c eb 0 1 t z m a n n 输运方程，再通过c h a p m a n e n s k o g 解 法得到l a t t i c eb 0 1 t z m a n n 输运方程的近似解，最后可得到n a v i e r s t o k e s 方程。运动 论层次虽然从分子的观点对流体作了离散化处理，但空间和时间依然是连续的。l g a 方 法则是一个更大胆的尝试，它不仅对流体和空间离散，而且时间也是离散的。 l b m 方法的理论基础是l a t t i c eb 0 1 t z 帆n n 方程。l b m 方法是对l a t t i c eb o l t z m a n n 方程的直接模拟。由于l b m 方法是在l g a 的基础上发展起来的，所以下小节将主要介绍 如何由l g a 的微观动力学方程得到l a t t i c eb o l t z m a n n 方程，然后讨论其碰撞项的含义 及如何简化，最后进一步简化得到l b m 方法中的很有前景的l b g k 模型。l b m 方法的基本 思想是：在简化连续b o l t z m a n n 方程的同时，保留流体力学的本质特征。基于对 b o l t z m a n n 方程考虑单一松弛时间的碰撞相( b h a t n a g e r ，g r o s sa n dk r o o k 于1 9 5 4 提出) ， 提出了l b g k 模型。它有处理高雷诺数的能力，但在某些环境下，这种方法在模拟过程 中可能产生数值不稳定性。 格子玻尔兹曼方法( l a t t i c eb 0 1 t z m a n nm o d e l s l b m ) 的发展可分为如下几个阶段： ( 1 ) 在格子气自动机中，l a t t i c eb o l t z m a n n 方程已经被用来计算流体的粘度。 ( 2 ) 为了消除格子气自动机中的噪声( 由布尔变量引起) l b m 作为一种独立的数值方法来 模拟流体流动。在原有的格子框架下用连续的分布函数代替布尔场分布。 2 东北师范大学硕士论文 ( 3 ) l a t t i c eb 0 1 t z m 锄方程碰撞项的线性化。 ( 4 ) b 0 1 t z i i l a n n 平衡分布代替f e r m i d i r a c 单粒子分布函数。 ( 5 ) 基于格子气自动机中的碰撞项被基于单弛豫时间下的b g k 近似所替代，即l b g k 模 型，碰撞过程不再被明确的定义。 两种方法在网格意义上是一致的。在二维空间有h p p 模型、f h pi 、i i 模型及九位 正方形模型，前两者为单速模型，后两者多速模型。对于单速格子气模型，粒子间的碰 撞没有能量交换，因此这种模型没有独立的能量方程( 包含于连续性方程之中) 。无论是 理论分析还是数值模拟都不能反映温度梯度的影响，它们只能用于模拟那些定常温度的 不可压流体的流动问题。多速模型可用来解决单速模型无法解决的传热问题。三维空间 中没有如平面上的正三角形那样的单纯形存在，不能得到具有足够对称的三维格子气自 动机。四维空间的f c h c 模型则具有足够的对称性，它在三维空间的投影，可得到一个 伪四维f c h c 模型。另外，在九位正方形模型的的启示下，又得到了正方体格子气自动 机模型。l 蹦方法的模型则可通过f c h c 模型投影到二位正方形网格上，也可应用六边形 网格进行数值模拟。 两种方法的优缺点比较：l g a 适合于不可压流体和低雷诺数流体但布尔变量使运算 数值稳定，对复杂边界有很强的处理能力。l 蹦可处理高雷诺数流体，对系统参数调节 更灵活，但忽略了多体相关性。并行算法则是此两种方法的共同优点。在两种方法之间 仍有相当太的发展空间。 1 2b 0 1 t z m a n n 方程及c h a p m 肌_ e n s k o g 方法 在导出l a t t i c eb 0 1 t z m a n n 方程之前我们先给出连续的b o l t z m a n n 方程，它是非平 衡态统计物理中描述单粒子分布函数的时间演化行为的方程。 a ，a h 芒耵= n ( 1 1 ) q 是碰撞积分，这里并没有考虑外力。尽管在导出过程中引入了许多假设，如假设两粒 子碰撞前是互不相关的( 这样的假设不可能总是成立的) ，已经不是严格的微观方程，但 和通常的宏观运动方程相比它提供的信息仍是比较微观的。 有了玻尔兹曼方程，通过求解该方程，可得到单粒子约化分布函数，再通过统计平 均，便可得到各种宏观量的统计表达式，但是玻尔兹曼方程中包含七个自变量，并且涉 及到偏导和多重积分，尤其是方程对于未知的分布函数是非线性的，这些因素使得方程 的求解十分困难。 c h a p m a n e n s k o g 方法常用来求解b o l t z 雎加方程。为了尽可能简化，假定体系是 接近于平衡的，为此假定在局域范围内单粒子的速度分布接近于平衡态的m a x w e l l 分布： ，o ( 芦，云，f ) = 玎( 删2 砬r ) 3 7 2 唧 _ ( m 2 足r ) ( 五一厅) 2 】 ( 1 2 ) 厅是气体的局部平均速度，和完全平衡态的m a x w e l l 分布有所不同，粒子数密度n ，温 度t ，速度于和宏观速度厅都是位置和时间的函数。 3 东北师范大学硕士论文 在体积、粒子数和宏观温度相同的条件下，气体在平衡态和非平衡态有相同的平均 密度、平均速度和平均动能，因而有： l 簿v = l ，删赢 肛击= ，5 击 ( 1 3 ) r 晰2 ) ( 弓一厅) 2 正拓= 且坍2 ) 仁一打) 2 ，。击 因为已经假定流体是接近平衡，可以把实际平衡分布f 近似表示为： ，= ，o ( 1 + ) ( 1 4 ) 其中妒1 “1 。将上式带入b o l t z l l l a n n 方程并作泰勒展开，保留矿1 的线性项，并利用 ( 1 3 ) 得到一级近似妒o 满足的方程，该方程不仅是线性的，碰撞积分中庐n 是唯一的未 知函数。一旦求得妒n ，则一级近似分布函数就是已知的，将( 1 4 ) 保留到二级近似庐扪， 在带入b 0 1 t z i a n n 方程，用同样的方法得到二级近似驴2 满足的方程，通过求解该方程 可得到二级近似分布函数，如此重复可得b o l t z 腧n n 方程的级数解： ，= ，o ( 1 + 妒1 + 矿2 + - ) ( 1 5 ) 这种方法在实际应用中相当繁琐，主要限于求解一级近似。 1 3b g k 近似 上述b o l t z 腿肌方程描述了单粒子分布函数随时间演化。流体密度、动量及能量可 通过适当的积分得到。理论上来讲很简单，然而在实施上由于碰撞项的复杂变得很困难。 因此很自然的要求对碰撞项进行简化，b g k 模型f 1 3 通常用来简化碰撞项，即单松弛时间 近似： q = 一二( 厂一，) ( 1 6 ) z ，o 是局域m a x w e l l 平衡分布，彳是松弛时问。这种碰撞项的替代必须满足质量守恒( 粒 子数守恒) 、动量守恒及能量守恒： q ( ，) 击= o p q ( ，) 击= o p 2 q ( ，) 斑= o ( 1 7 ) 东北师范大学硕士论文 1 4 流体力学方程 在非平衡态统计物理中，一般来说宏观量不仅与时问有关，还与空间位置有关。要 表征一个系统的状态需要在物理空间的每个点得到相关物理量的值，它们是随时间变化 的。例如质量密度p ( 产，f ) ，宏观速度 ，这些物理量构成了宏观场，如密度场，速 度场。非平衡态统计物理的一个重要任务是从单粒子分布函数随时间演化方程导出宏观 运动方程。 碰撞过程中守恒量因其值不受碰撞过程的影响，它们衰变到平衡值的速度要比非守 恒量的衰变速度慢得多。因此系统的非平衡态行为往往可以用守恒量的场来表征。描述 碰撞守恒量的场的运动方程通常称为流体力学方程。 首先给出几个重要物理量的定义式： 密度：p ( f ，0 = i ( f ，石，f ) 应( 1 8 ) 宏观流速： = i 蛳 ( 1 9 ) 热速度： 吾( ，r ) = 云一面 ( 1 1 0 ) 应力张量： 能量密度： 温度： 热通量： 利用( 1 1 ) 、( 1 7 ) 及( 1 8 ) 一 = m j 印肺 声( 产，f ) = 片州已2 弦 r ( 尹，f ) = ，( 二足丑) 彳= b 榭2 准 ( 1 1 3 ) 可以得到流体力学方程： 詈w ( 历瑚讲 a 。( 肛，) + a 口( 肛群即+ ) = o m p ( + ) 】+ 伊口( 圭“2 删+ + = o 1 5l b g k 方程的导出 首先根据b g k 近似和方程( 1 1 ) 我们得到： 矽a h 石w = 一 ( ，一，蚋) ( 1 1 1 ) ( 1 1 2 ) ( 1 。1 - 3 ) ( 1 1 4 ) ( 1 1 5 ) ( 1 - 1 6 ) ( 1 1 7 ) ( 1 1 8 ) 分布函数是位置、速度和时间的函数，我们首先在速度空间选取有限且必要的方向，完 成对速度空间的离散化：亭j 己；f = l ，6 ，五为格子速度。b 为选取方向的个数。对 5 东北师范大学硕士论文 应的分布函数为石( 尹，f ) 。这样我们就得到速度空间离散的b o l t z l i l a n n 方程： 奶a f + 己= 一圭( z z ) ( 1 2 0 ) 我们再进一步把时间和空间都离散化，空间间隔、时间间隔及格子速度有如下关系： 氲= 己f ( 1 2 1 ) 再根据有限差分的方法得到l b g k 方程： z ( i + 缸，f + f ) 一z ( 手，f ) = 一! 兰【z ( i ，) 一z 。( i ，f ) 】 ( 1 2 2 ) 密度和宏观速度公式变为： p ( 训= z ( 训 ( 1 2 3 ) 口( j ，f ) = z ( i ，f ) c 胁 ( 1 2 4 ) f 根据l b g k 方程及格子空间算法过程可描述为如下几步： ( 1 ) 设定初始密度和初始动量来计算平衡分布函数。 ( 2 ) 选定松弛时间由l b g k 方程得出碰撞后的分布函数，传播分布函数。 ( 3 ) 重新计算各格点密度和动量。 ( 4 ) 计算此时的平衡分布，重复第二步。 ( 5 ) 这种过程的循环最终实现了由局部平衡到达整体平衡。 1 6l b m 方法应用举例 ( 1 ) 微观流动模拟：已经成功应用到了带压力的微管流动、微管【2 】中的两相混合等。动 力学理论( b o l t z m a n n 方程) 能够处理微观流动效应。 ( 2 ) l b m 方法成功应用于多相和多组分流动问题的模拟。多相和多组分l b m 模型涉及到 流体相变问题及流体包含不同特性的分子两相混合问题。n 6 l i d 0g o n z 矗l e z s e g r e d 0 等人 3 】利用三维极性l b m 方法研究了二元不相容流体相分离动力学过程。在以下几 方面的应用中体现出了它的计算优势：悬浮粒子碰撞问题【4 】、高分子溶液体系嘲及 不同流体表面动力学模拟。 ( 3 ) 此例【6 】是l b 方法在不均匀空穴体系中的流场示意图，体系左，右，下三边为固定墙， 上边加入从左至右的剪切。当系统达到稳定时，体系内会形成一个不均匀的流场。 特别是当雷诺数达到一定数值时，体系的左下方会形成一个涡旋。 东北师范大学硕士论文 图卜1 不同r e 数下的涡流函数等值线图：r e = 5 0 0 0 ，7 5 0 0 ，1 0 0 0 0 ( 上部分：从左至右) ， r e = 1 5 0 0 0 ，2 0 0 0 0 ，5 0 0 0 0 ( 下部分：从左至右) 示意图 东北师范大学硕士论文 第二章l b m 的基本模型 本章中介绍了二维空间中的两种模型：d 2 q 7 模型和d 2 q 9 模型。以及在l b g k 方程中 至关重要的平衡分布函数得出。通过各向同性和伽利略不变性要求可以得到它在具体模 型中的形式。接下来将通过c h a p m a n e n s k o ge x p a n s i o n 由l b g k 方程得到流体力学方程， 从而得到粘度系数与模型参数之间的关系。 2 1d 2 q 7 模型和d 2 q 9 模型 在模型中，我们将整个体系划分成若干个间距为一的正三角形网格，每个节点有六 个邻位，如图2 1 所示： 对应的格子速度为： 32 量 图2 1d 2 q 7 空间格子 c o = ( 0 ，o ) c ，= ( c o s 孚舢等) c f = l ，6 ( 2 1 ) c 是传播速度，表示粒子在两个节点之间运动速度。c 的值与空间位移的关系为： 6 s = c & ( 2 2 ) 在很多模拟中，c 和厨被设定为1 。 这种的划分并不是任意的，为了能够得到较满意的宏观流体力学方程，格子张量必 8 东北师范大学硕士论文 须满足各向同性要求。格子张量定义如下【7 】： ( 2 3 ) 根据( 2 1 ) 和( 2 3 ) 式可以得到d 2 q 7 模型中格子张量的结果： 三口= c ，口= o ( 2 4 ) j = c 以c 妒= 3 ， 卿= c 垴c 妒c 纱= o f = 军q 口c 猡q 肋= 号( 如+ 成) ( 2 5 ) ( 2 6 ) ( 2 7 ) 在许多后来提出的多速l b m 模型中，四阶格子张量通常是非各向同性的，但是通过 引入权重因子，可以使四阶张量具有各向同性，定义一般性的格子张量： gm 2 口一= wf c 川c 川c 胁 ( 2 8 ) f 对于d 2 q 7 模型，一般性的格子张量具体表达式为： c 地= o m c ，口c 猡= 刀2 f c f 口c 班c 秒= o j ( 2 9 ) ( 2 1 0 ) ( 2 1 1 ) 心q ，= 心( 知+ 妨成) ( 2 1 2 ) ， ( 2 9 ) 一( 2 1 2 ) 也就是各向同性条件，权重因子及未知参数可在求解平衡分布函数过程 中求得。d 2 q 7 模型的格子张量已经满足了到四阶量的各向同性，一般性的格子张量的 引入并无实际意义，只是为了说明在格子本身对称性不足时所要引入的张量形式即一般 性的格子张量，这在d 2 q 9 模型中就能够体现出来。 - 9 - 铆 c 2口 c 口 c ， = 铆 2口口 乙 东北师范大学硕士论文 下面给出d 2 q 9 模型及其对称关系【8 】： 对应的格子速度为： 6 3 t 2 4 图2 2d 2 q 9 空间格子 c o = ( o ，o ) 5 1 8 c l ，3 ，乞，4 = ( 1 ，o ) ，( o ，1 ) ( 2 1 3 ) g 6 7 ，8 = ( 1 ，1 ) d 2 q 9 模型，一般性的格子张量具体表达式为： 一嵋c ，口= o f 权重因子为： 军铲吾c 2 c 地c 谚c f ，= o f 军心c 舻泸q 一占= 吉c 4 ( + 妨蛾) 4 2 万 l 2 万 1 嵋= 丽 f = 0 f = 1 ，4 f = 5 ，8 ( 2 1 4 ) ( 2 1 5 ) ( 2 1 6 ) ( 2 1 7 ) ( 2 1 8 ) 东北师范大学硕士论文 为了能够得到n a v i e r s t o k e s 方程，必须选择适当的平衡分布函刿9 1 ： 厂。) = p w ， 1 + 3 c ，口“，口+ 昙( c ，口甜，口) 2 一导“2 】 ( 2 1 9 ) 厂= p w f 1 + 3 c f 口“j 口+ 詈( c j 口甜f 口) 2 一号“2 】 ( 2 1 9 二 ， 么 2 2 平衡分布函数 我们所要选取的平衡分布函数既要满足质量守恒、动量守恒条件，又要满足动量通 量张量具有流体力学中的关系。速度各阶量的伽利略不变性要求也很重要，会使我们在 后面得到的宏观方程具有伽利略不变性。这里只给出求解平衡分布函数所需要的条件， 不做详细的计算，但过程具有一定的普遍性。 在低m a c h 数即矧很小时1 们，若对m a x w e l l 平衡分布进行t y l o r 展开并略去高阶项， 则有： z o = w f ( 1 + 彳甜口c 衄+ b 甜口“卢+ c h 口甜c f 口c 谚) ( 2 2 0 ) z o 物理意义是当宏观速度为零时的平衡分布函数。在d 2 q 7 模型中，因为i = 1 ，6 速度模相同，具有相同的权重w ，i = o 具有单独的权重w o 。d 2 q 9 模型中的速度权重比 d 2 q 7 模型多一种。 为了求得各系数，必须考虑伽利略不变性1 1 1 ，首先定义m j ： m ，= ( 历一己) ，z o ( 2 2 1 ) 伽利略不变性要求 m ， ) = m ，( o ) ，j f = o ，1 ，2 ，3 。 ( 2 2 2 ) 将平衡分布函数代入( 2 9 ) 一( 2 1 2 ) 和( 2 2 2 ) 式得到方程组，最后解出系数。 利用这些条件就构建了一个具有各向同性和伽利略不变性的二维模型： z o = 嘶( 1 + 4 “口c 衄一2 甜口，+ 8 甜口“，c 胁c 泸) ( 2 2 3 ) 并且有= 詈，w = 告。 质量守恒、动量守恒条件可表示如下： p = z = z 。 ( 2 2 4 ) 东北师范大学硕士论文 肛口= z c ，口= z 。c 以 ( 2 2 5 ) ， 2 3 从l b g k 方程到流体力学方程 使用c h a p m a n e n s k o g 多尺度展开方法，在确定的碰撞模型下可以由b o l t z m a n n 方 程得到n a v i e r s t o k e s 方程和输运系数。在格子气自动机时代就己被用来推出e u le r 方 程。l 蹦方法的出发点就是l b g k 方程，所以由此得到宏观方程是非常必要的，同时将得 到模型的粘度与基本参数关系式。展开参数被称为k n u d s e nn u m b e r ( k n ) 记为，它的 物理意义可解释为碰撞时间和流动时间之比或解释为粒子运动平均自由程与系统特征 尺度之比。它是一个小参数。 求解b 0 1 t z m a n n 方程的c h a p 腿n e n s k o g 方法，为了使问题尽可能简化，首先要假 定体系是接近于平衡的，这样便可得到： 质量守恒、动量守恒条件变成： z = z o + 锐1 + 2 z 2 + o ( 3 ) ( 2 2 6 ) = c 妇= o f = c 衄= 0 jj ( 2 2 7 ) ( 2 2 8 ) 矿= 卜学步+ 等2 一筠吼山亿3 。， ip 幌一等) 待o 、。 三口= c 蛔= o = 弘旷害 ( 2 3 1 ) ( 2 3 2 ) 东北师范大学硕士论文 锄= 邻= o ( 2 3 3 ) = 军q 一万= 瓦宗酉( 勃成勃) ( 2 3 4 ) 利用( 2 2 9 ) 一( 2 3 4 ) ，我们得到： z o = p ( 2 3 5 ) z 0 1 = 肛口 ( 2 3 6 ) 军z ) 印= 脓坳+ 华妨 ( 2 3 7 ) ；z 锄= 彘( z ，+ 妨+ 坳) ( 2 3 8 ) 现在我们可从l b g k 方程出发来利用上面给出的关系式得到不同k n 阶所对应的方程 了，首先对方程左侧进行t y l o r 级数展开( 取址= 1 ) 并略去高阶项( 保留到二次导数) ： z ( i + 五，f + 1 ) 一z ( i ，f ) 圭【( a ，+ c 船a 口) + 三( a ，+ c ，口a 口) 2 】z ( j ，f ) ( 2 3 9 ) 将( 2 2 6 ) ，( 2 2 7 ) 代入上式得： 【( 国l ，+ 2 a 2 r ) + 岔j 口a l 口+ 去( 国l ，+ 2 a 2 ，+ 岔f 口a l 口) 2 】( 石o + 玩d + 2 彳2 ) = 一圭( z + 玩o + 2 z 2 一z ) 二 ( 2 4 0 ) 展开( 2 4 0 ) 左侧并忽略o 忙3 ) 项： ( 国l ，+ 占2 a 2 ，) ( z 。+ 锐1 + 2 z 2 ) + 衄a l 口( z 。+ 玩1 + 占2 z 2 ) + 争2 c 地c 泸a 。口a 。，+ 2 岔衄a 。口( 国“+ 2 a ：，) + ( 国。，+ 2 a ：，) 2 ( z ) + 矾m + 9 2 z 但) 圭( a 。，正 ) + q 口a 。口z ) + ? 2 ( a 。，z 。) + a ：，z o ) + c 船a 。口z a ) + 三( c 蛔c 泸a 。口a 。，z + 2 c f 口a l 口a lr z o + a l f a “z o 】 右边= 一吾( 阢+ 2 彳但) 比较系数得到： a 。，z o ) + c 衄a 口z ) = 一z 。) 求和i 并利用( 2 2 7 ) 、( 2 3 5 ) 和( 2 3 6 ) ， ( 2 4 1 ) 东北师范大学硕士论文 ( 2 4 1 ) 宰钿得， a 户+ a l 口( 肛口) = o ( 2 4 2 ) a 。，嘶o ) + 印。0 ) = 扫z m ( 2 4 3 ) 求和i 并利用( 2 1 9 ) 、( 2 3 6 ) 和( 2 3 7 ) ， a “刚h 肛卿) _ - a d 华】 ( 2 4 4 ) 比较2 项系数： ( a 。，z m + a ：，z ) + a 。口z m ) + 三( a 。口a 。卢z ) + 2 a 。口a 。，) + a 。a ，z ) = 一昙z 。) ，( 2 4 5 ) 求和i ，我们求出( 2 4 5 ) 第二个括号中的如下几项： 军a 。口a 。，z 。丝a ，口a 。，【丛与笋星】+ a 一从肌坳) ( 2 4 6 ) 利用( 2 4 4 ) 有 代入( 2 4 6 ) 式 a 一卿) - - k 掣h ( 刚 ( 2 4 7 ) a 。口a 文肛口) = 一a 。口a 。，【丛与盟】一a 。，a 。，( 肛，) ( 2 4 8 ) c 船c 泸a l 口a l ，z 。= - a l 卢a l ，( 肛，) f 2 a 。口a l ，彳o = 2 a l 口a 。，z o ff ( 2 4 9 ) 2 a l 口a l ，( 眺) ( 2 5 0 ) a 。，a 。，z 。= a 。，a 。，p 坚竺一。口( 纵) ( 2 5 1 ) j 一 求和( 2 4 9 ) 、( 2 5 0 ) 及( 2 5 1 ) c f 口c 泸a l 口a l ，z 。+ 2 c 橱a i 口a i ，z o + a a l ，z o = o ( 2 5 2 ) fff 这样( 2 4 5 ) 简化为 东北师范大学硕士论文 a 。+ a 2 + 啊一 利用( 2 2 7 ) 、( 2 2 8 ) 和( 2 3 5 ) 得( 2 4 5 ) a 2 ，p = o ( 2 5 3 ) ( 2 5 4 ) ( 2 4 1 ) 木得， p 。，c ：d + 九c ： + q 廊a 。z d ) + 去心肠印。伊伊+ 狮。a ，) + 扛a 。，“0 ) ) = 弓嘶圆 q _ 5 由( 2 4 1 ) ，得 z 1 = 一彳( a l ，z + c 蛔a l 口z o ) ( 2 5 6 ) 代入上式并求和( 利用了( 2 2 7 ) 和( 2 2 8 ) ) ： d 各c f + c i 脚k m k + 印。) 一勘p 。一脚扩+ a a ，c f 删) + 残a 。，( c f 一。+ a 。，q ) = o ( 2 5 7 ) 由( 2 4 2 ) 可知左边第三项为零， 利用( 2 3 5 ) 、( 2 3 6 ) 和( 2 3 7 ) 得 a ：，油y ) + 睦一力擘学a 。a 口( p 妨) + k a 。【彘( ”口勤+ 即+ “，) 】) i 。 ( 2 5 8 ) 展开计算： a 。，a 。口( 畴) = a 。，( a 。，p ) = - a 。，a ，口( 口) a ，口a ，( 口勤) = a ，口a l ，( 肛口) = a ，a ，口( 口) a l 口a l ，( 卢) = a l ，a i ，( 肛，) = a l ，a l 口( 肛口) a 。口a 。，( 肛，) = a 。口a 。口( 肛，) 代入上式整理，得 九( ，) + 哇一力 尘号显 一a 。户。口慨) 】+ 吉缶吼户。口( 以) + 哇一力苫a 。口( ，) = c 东北师范大学硕士论文 令 则有： 州而彘一半，刀：萼 ( 2 5 9 ) a 2 ，( p “厂) = al y al 口( 尸甜口) + 刀al 口al 口( p “，) ( 2 6 0 ) 对一阶和二阶密度方程求和( 2 4 2 ) + ( 2 5 4 ) ，得到连续性方程： a ，p + a 口( 口) = o ( 2 6 1 ) 对一阶和二阶动量方程求和( 2 4 3 ) + ( 2 5 5 ) ，得到n a v i e r s t o k e s 方程： a 舳小a ( 嗍) = 一【世笋】+ 闰口a 卢( 脚) + 稠户，( 纵) ( 2 6 2 ) 可知刀为动力粘性系数，为本体粘性系数。 小结： 本章主要介绍了l b 方法重点中的d 2 q 7 和d 2 q 9 模型，本文主要采用d 2 q 9 模型。给 出了如何从l b g k 近似出发得到流体动力学方程，从而得到平衡分布函数的具体形式。 由于l b 方法中参数调整很灵活，因此它可在模拟计算中可以发挥更大的优势。 东北师范大学硕士论文 第三章不相容多组分体系的相分离模拟 自由能l b m 能够模拟多种体系相变。通过引入正确的非理想流体应力张量并改写基 本模型中的平衡分布函数，改变了b 0 1 t z l a n n 方程中的碰撞项，使得体系最终平衡与我 们选择的自由能相对应。它与基本模型相比，包含了更多的物理真实性，减少了离散化 带来的非物理问题。边界条件在l b m 中直接影响模拟的结果，无论是模拟精确度，还是 体系内各组分分布，都与边界条件的设置相联系。本章采用剪切边界条件，对不相容多 组分体系进行了相分离模拟。 3 1 自由能l b m 理论框架 在l b m 模型中，在同样的离散格子中用粒子速度分布函数代替流体粒子以离散速度 迁移，在碰撞过程各分部函数松弛到相应的平衡值。在不同方向上粒子速度分布函数的 平衡值由局部平衡分布函数确定，通过选择恰当的平衡分布函数可以得到相应的 n a v i e r s t o k e s 方程，从而体现不同的流体特性。l b m 模拟流体流动是从流体的非平衡 态出发依据格子b o l t z l 玎a n n 方程向平衡态演化。在演化过程中，流体的单粒子速度分布 函数满足格子b o l t z 衄n n 方程，流体的宏观量( 如密度、速度) 由单粒子速度分布函数定 义。 l b 方法最关键处在于平衡分布函数的选取，使流体粒子在每个格点上满足质量、动 量守恒律，然后使用多尺度展开技术( c h 锄p 腿n e n s k o g 展开) 可以恢复流体动力学的连 续方程和n a v i e r s t o k e s 方程。l b m 演化过程主要分两个步骤：( a ) 迁移，粒子从一个节 点在一个时步内，以恒定的速度运动到相邻节点：( b ) 碰撞，在一个节点上与相邻节点运 动来的粒子发生碰撞，根据质量、动量和能量守恒规则改变粒子的速度，然后各个粒子又 以改变后的速度迁移。这两个步骤交替循环，直到整个体系达到收敛。 l b m 能够正确模拟单组分流体的流动【”】，在处理复杂边界问题时也非常适用。随着 l b m 的不断发展，它也被应用到多相流体【1 4 】问题研究中。自由能l b m 理论最早是由s w i f t 等人提出【1 5 】，后来又扩展到模拟二元流体【1 6 】，它的基本物理思想是流体系统演化到一个 与给定自由能相对应的平衡态。 3 2 二元体系模型 二元模型中我们不分别考虑两种组分的密度，而是考虑格点处两种组分的总密度p 和密度差伊，这样就要有两套分布函数 和堙， 。两套分布函数的演化方程为仉1 乃1 8 1 ： - 1 7 - 东北师范大学硕士论文 厢+ c f ，f “) 一胛力一去咪础) 叫o ) ( 酬 ( 3 1 ) “n h 1 ) - “如) 一去噱( 刁) _ g f ( 例 ( 3 2 ) 质量守恒和动量守恒条件为： z k p 一c i 口= p ”口 l g j = 矿 j 为了得到正确的流体力学方程需作如下的定义： z c j 口c 泸= 刖口“声+ f g f c f 口= 缈口 f g j 。c i 口c 泸= 缈口材，+ r 万妒 f ( 3 3 ) ( 3 4 ) 通过这些方程便可求得平衡分布函数z 和g ，的系数。( 3 4 ) 中的应力张量和化学 势雏决定了模型热力学性质，所以要选择恰当的自由能。通常选取l a n d a u 自由能1 9 1 ： ，= 每l v 纠2 + y ( p ，仍丁) ) 方 ( 3 5 ) 自由能密度的形式为【2 0 】： 帅们= 等( 1 一争却+ 吾( p + 加呼) + 吾( p 一舢洋) ( 3 6 ) 其中旯为两流体的相互作用强度。由此得到两流体间的化学势为： 雏：要：一竺+ 要l n ( 业) 一胛z 缈 ( 3 7 ) 雏2 面一荔+ i m 茄卜缈 7 压力张量为： 吡丁+ 蝴) + 砥老考一吾v 2 嘞) ( 3 8 ) 根据以上结论，我们从理论上分析平衡态所对应的序参量值伊。相分离过程将最终 使序参量达到平衡值，平衡值由如下平衡条件决定： ( 仍) = ( 仍) ( 3 9 ) 东北师范大学硕士论文 当体系演化到两组分密度几乎不变时，也可以认为达到平衡。 二元体系相分离的模拟中，我们选取的自由能形式的如下【2 1 】： ，二眯矿+ 兰c v 卅 根据热力学原理，对于这种自由能的选取，可得出两种物质间的化学势： 雏= 篙叫州羽2 缈 ( 3 1 0 ) ( 3 1 1 ) ( 3 1 0 ) 形式的自由能常被用来研究相分离。式中级数项与本体特性相关，梯度项与 界面性质相关。当a 0 时，自由能只有一个最小值，相应的缈= o ，当a 0 时，自 广= 由能有两个最小值，相应的缈= 1 f 一等。选择a _ 一1 ，b = 1 ，平衡时缈= 1 。 y d 3 3 边界条件的引入 在体系相分离的过程中，可受不同条件的影响，如温度、粘度、初始浓度、各 组分比例等，其中，边界条件在体系演化过程中起到的作用是不可忽略的。接下来 我们介绍两种本文用到的边界条件：周期性边界条件和剪切性边界条件。一 由于我们选取的格子空间只是整个空间的一部分，为了能表现出整个空间的情 况，就要用到周期边界条件。a 为我们所要模拟的格子空间，其上边界对应