（计算数学专业论文）求解含源项对流扩散方程的格子bgk模型.pdf

华中科技大学硕士学位论文 摘要 建立在微观模型上的格子b o l t z m a r m 方法是近年来发展起来的种模拟流体 流动新的计算方法。与传统算法相比较，格子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 t t z m a n n 模型中，现有的模型中， 都有对源项做不同的假设，这样的话，使得这类模型的适用范围变得狭窄。 本文给出了三个求解含源项对流扩散方程的格子b g k 模型( 最常见的格子 b o l t z m a n n 模型) ，这三个模型是针对原有模型都有对源项的不同的假设这一不足 而提出的，通过将演化方程重写，在演化方程的右端加上一个小量，使得本文的 三个模型都不需要对源项做任何假设，从理论上说明了改进模型较原有的模型有更 广的使用范围。第三章中我们使用改进模型对含源项的扩散方程、反应扩散方程、 含源项的对流扩散方程进行了模拟，实验结果和解析解吻合的很好，和原有模型 模拟的结果比较进一步说明了改进模型较原有的模型有更广的适用范围和更高的 数值精度。第四章研究了改进模型i 对c i m a 系统的模拟，并给出了二维空间中 图灵斑图的形成过程。 关键词： 格子b o l t z m a n n 方法，含源项的对流扩散方程，改进模型，c i m a 系统，图灵斑图。 华中科技大学硕士学位论文 a b s t r a c t t h el a t t i c eb o l t z m a n nm e t h o di san e wn u m e r i c a lt e c h n i q u eb a s e do nk i n e t i c t h e o r yt os i m u l a t ef l u i df l o w s i nr e c e n ty e a r s c o m p a r e dw i t ho t h e rc o m p u t a t i o n a lf l u i d d y n a m i c sa p p r o a c h e s ，t h e m e t h o di s s i m p l e ，i n t r i n s i c a l l yp a r a l l e l ，m a de a s y t o i n c o r p o r a t ec o m p l e xb o u n d a r yc o n d i t i o n s a l t h o u g h ，t h e r eh a sap r o g r e s si ne m p l o y i n g t h el a t t i c eb o l t z m a n ne q u a t i o nf o rf l u i dd y n m l a i c sa n dm o d e l i n gc o m p l e xp h y s i c si n f l u i d sr e c e n t l y , t h e r ea r es o m ep r o b l e m sn e e dt ob es o l v e d s u c ha st h e r ea r ed i f f e r e n t a s s u m p t i o n s a b o u ts o u r c et e r mi nt h e o r i g i n a l l b g km o d e l sf o rt h e c o n v e c t i o n d i f f u s i o n e q u a t i o nw i t hs o u r c e t e r m i nt h i sp a p e rt h r e ei m p r o v e dl b g km o d e l sa r cg i v e n t h e s et h r e em o d e l sr e m o v e a l la s s u m p t i o n sb yr e w r i t et h ee v o l u t i o ne q u a t i o n t ot e s tt h ei m p r o v e dl b g km o d e l s p r o p o s e di nt h ep a p e r , n u m e r i c a ls i m u l a t i o n so f d i f f u s i o ne q u a t i o nw i t hs o n r c et e r m ， r e a c t i o n - d i f f u s i o n e q u a t i o na n dc o n v e c t i o n d i f f u s i o ne q u a t i o nw i t h s o u r c et e r ma r e c a r r i e do u ta t c h a p t e r3 i t i sf o u n dt h a tt h es i m u l a t i o nr e s u k so fi m p r o v e dl b g k m o d e l sa r ei ne x c e l l e n ta g r e e m e n tw i t ht h ea n a l y t i c a ls o l u t i o n s ，ni sa l s os h o w nt h a tt h e n u m e r i c a la c c u r a c yo ft h e i m p r o v e dl b g km o d e li s m u c hb e t t e rt h a nt h a to ft h e e x i s t i n gm o d e l s a tc h a p t e r4w e u s et h ei m p r o v e dl b g km o d e lit os i m u l a t et h e c i m as y s t e m t h ep r o c e s so ft a r i n gp a t t e r n sf o r m si nt w o - d i m e n s i o n a ls p a c ea r e g i v e n k e y w o r d s ：l a t t i c eb o l t z r n a n nm e t h o d s ，c o n v e c t i o n d i f f u s i o ne q u a t i o nw i t hs o u r c e t e r m ，i m p r o v e dl b g km o d e l s ，c i m as y s t e m ，t u r i n gp a t t e r n s 独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究 成果。尽我所知，除文中已经标明引用的内容外，本论文不包含任何其他个人或集体己 经发表或撰写过的研究成果。对本文的研究做出贡献的个人和集体，均己在文中以明确 方式标明。本人完全意识到，本声明的法律结果由本人承担。 学位论文作者签名：商 照 日期：噼年 月3 日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即：学校有权保留 并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授 权华中科技大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以采 用影印、缩印或扫描等复制手段保存和汇编本学位论文。 保密口，在年解密后适用本授权书。 本论文属于不保密口。 ( 请在以上方框内打“”) 学位论文作者签名：覃p 癌 日期：砰年占月8 日 指导教师签名：匆绯启 日期；口弘年钥彦日 华中科技大学硕士学位论文 1 1 引言 1 绪论 流体力学是研究流体运动规律的门学科，其范围非常广泛。经过多年的发展， 流体力学已经取得了丰硕的成果。但是由于流体力学的复杂性，流体力学还面临 着巨大的挑战。数学上，流体力学的复杂性反映在控制方程的非线性上。除了某 些简单的情况外，很难获得这些偏微分方程的精确解。对于大多数的实际问题， 必须采用试验方法或者是数值解法。 随着计算机技术的不断发展，数值模拟方法已经成为流体力学研究中的一种 重要的手段和工具，并且日益受到人们的重视。这一领域已经反展成为流体力学 的一个分支，即计算流体力学( c f d ) 。 设计模拟流体运动的数值方法有两种途径，即基于宏观连续模型的自顶向下 方法和基于微观离散模型的自底向上方法。传统的计算流体力学中的数值方法大 多是利用自顶向上方法设计的。这类方法以非线性的微分方程为出发点，采用有 限差分、有限体积、有限元或有限谱等离散方法对微分方程进行离散，得到代数 方程组或者是常微分方程系统，然后再用标准的数值方法求解。虽然这类自顶向 下方法比较直观，但仍然存在许多不足。例如，在这类方法中人们往往看重分析 这类从连续微分方程到离散代数方程的截断误差，而忽视了离散过程中某些物理 量的守恒性。对某些系统而言，为了得到合理的结果，这种守恒性的要求是非常 重要的。另外，数值稳定性也是这类方法的一个重要问题。 与上述自顶向下的方法完全不同，自底向上的方法的出发点是流体的微观离 散模型。我们知道，流体的宏观运动是大量流体分子微观运动的统计平均结果， 单个分子的运动细节并不影响宏观运动的特性。例如，气体和液体有完全不同的 微观分子结构，分子的运动方式也大不相同。但是，对同一类流体( 圆柱绕流) ，只 要流动的r e y n o l d s 数相同，气体和液体的运动特性是完全相同的。因此，我们可 以构造这一种人工微观模型，使其在保持真实流体的基本特征的前提下，结构尽 华中科技大学硕士学位论文 可能的简单，粒子运动的细节尽可能的简化，且其宏观的统计特性符合客观的运 动规律。这样，我们就可以借助这种人工的微观模型模拟现实的流体系统。近年 来备受人们关注的格子方法( 包括格子b o l t z a m n n 方法及其前身格子气自动机) 就属 于这类模型。 在格子方法中，流体被抽象为大量的微观粒子，并且这些粒子根据某些简单 的方式在规则的离散格子上碰撞和迁移。通过对粒子的运动进行统计，就可以得 到流体的宏观运动特性。从离散的网格来说，格子方法具有e u l e r 方法的属性；从 离散的粒子观点来说，格子方法又具有l a n g r a n g e 方法的属性。格子方法的这种粒 子属性也使其具有许多常规数值方法没有的独特优点，如物理图像清晰，边界处 理容易和本质并行等。格子方法还提供了联系宏观和微观的可能性和现实性。它 既能直接计算流体的粘性，又可以在一定的条件下逼近n a v i e r s t o k e s 方程。同时， 格子方法这种用简单模型实现复杂系统的数学建模方法，也打破了传统的建模观 念，为其他复杂系统的建模提供了新的途径。从这种意义上来说，格予方法不应 当仅仅被认为使一种计算方法。 1 2 格子自动机模型概述 2 0 世纪5 0 年代初，j y o n n e u m a n n 提出了元胞自动机( c e l l u l a ra u t o m a t a 简记c a ) 。在c a 这个重要的学科生长点上，人们提出了格子自动机( l a t t i c eg a s a u t o m a t a ，简记l g a ) 和格子b o l t z m a n n 方法( l a t t i c eb o l t z m a n n m e t h o d ，简记l b m ) 。 它们都是离散演化动力系统的计算模型经典，这些模型均采用离散的空间布局( 格 点) 和离散的时矧间隔( 时间步) ，在微观或细观层次| 二，基于局部动态平衡原理， 着眼于个体自适应的演化运动描述，模拟复杂系统的复杂现象。它们在研究、探 索复杂现象的形成机理、发展趋势等方面取得了很大的成功。 1 2 1 元胞自动机 元胞自动机【1 】是j v o n n e u m a n n 了模拟生物的自复制功能而设计的哟中简单 的计算模型。他曾构造了一个2 9 种状态的c a 具有自复制功能。1 9 8 4 年，c g _ 2 华中科技大学硕士学位论文 l a n g t o n 构造了更为简单的具有自复制功能的c a ”，并迅速的扩展唯一个学科分 支人工生命。另有一些理论工作证明了c a 与图灵机之间的计算等价性。1 9 7 0 年，j c o n w a y 的“生命游戏”形象地表明，随机的初始布局是如何经过元胞间的 局部作用、微观决策、同步协调，丽演化成一定时空结构的。1 9 8 4 年前后，s w o l f r a m 的杰出工作 卜4 1 将c a 的研究与应用推向新的高潮。c a 在斑图形成、枝晶生成、 交通流、地理信息系统、图像处理、曲面拟合造型【5 1 等方面取得了良好的应用效果。 c a 还是一个重要的学科生长点，l g a 与l b m 均由c a 反展而来，并在定量计算 与分析方面取得了重大的突破。 一维c a 的基本概念 6 1 ：设在一条直线上按等间隔方式分布着完全相同的一系 列元胞( c e l l ) ，抽象地认为这里的每个元胞是具有一定功能的信息处理机或自动 机，每个元胞的状态只有有限种。假设直线在两个方向上没有限制，因此有无限 多个元胞。某个时刻全体元胞的状态集合可以用双侧无限的符号序列表示，这样 的序列称为元胞自动机的一个构型。例如： a o k ( a ! ? a ：a l ) 表示c a 在时刻l 的构型，a 为时刻f 位于空间位置i 处的元胞的状态。位于f 处 的元胞在下一时刻o + l 时亥6 ) 的状态a ，由其t 时刻某个邻域f ( i ，t ) 中的有限构 型所决定，简记为 a l “一十( n r ( f ，f ) ) = 厂( a 婴，a 葛，a i ) ，a 黯a “i + ，) 其中，称为领域半径，称厂为元胞自动机的局部演化规则，且，与i ，f 无关。 因此，c a 可用空间c 的、有限状态集( s ) 、领域( ) 和演化规则( 十) 四要素描述。 初等元胞自动机是最简单的c a 。其状态集海 0 ，1 ) ，邻域半径r = l 。这时，邻 域有限子构型( a ! ) a a ：i ) ) 仅有8 种，即( ( 1 ，l ，1 ) ，( 1 ，1 ，o ) ，( 1 ，0 ，1 ) ，( 1 ，0 ，o ) ，( o ，l ，1 ) ，( o ，】，o ) ， ( 0 , 0 ，1 ) ，( 0 ，0 ，o ) ) 。只要给出这8 种有限子构型到状态集s 的一种映射，元胞自动机 也就确定了。例如；可以定义如下的局部演化规则 华中科技大学硕士学位论文 1 1 11 1 01 0 11 0 00 1 10 1 00 0 10 0 0 j ，上j rj山j ，j ，山 o1o0l1o0 显然，这样的演化规则共有2 8 = 2 5 6 种，这2 5 6 种c a 统称为初等元胞自动机。 初等元胞自动机按其演化规则所对应的二进制数( ( 0 1 0 0 1 1 0 0 ) 2 = ( 7 6 h o ) 进行编号。上 例便是第7 6 号元胞子定级的演化规则。 给定系统的初始构型( 例如，可给每个元胞赋一随机状态值) ，c a 便可根据演 化法则按时间步同步演化。 二维及高维c a 是一维c a 的简单推广。二维c a 的元胞分布在平面上，常用 的有正方形网格、三角形网格等，适当定义领域概念，演化规则即可。 元胞自动机的主要特点可概括如下：空间离散、时间离散、状态有限、演化规 则局部、操作同步并行。 1 2 2 格子自动机 1 9 7 2 年，j h a r d y 、yp o m e a u 和o ，d ep a z z i s 将c a 应用于流体运动模拟，提 出了第一个l g a 模型( 现称为h p p 模型) 【”。这个模型中，将流体视为由微观粒子 所组成。为了模拟粒子的运动，引入了元胞的方向属性，碰撞、迁移等操作。该 模型将平面划分为正方形网格，每个网格点由上、下、左、右四个方向与其他网 格点相连。所有粒子根据其方向属性在下一时间步沿网格连线并发的运动到下一 网格点。当两个相对运动的粒子碰撞后，它们分别沿另外两个方向运动。一般的 l g a 采用互斥的原则( 或称为不相容原则) ，即各网格点的每个方向上( 包括静止的， 即0 方向) 的粒子数不超过1 。因而，在h p p 模型中，每个网格点上，包括静止粒 子在内，允许最多5 个粒子。 1 9 8 6 年，u f r i s e h 和yp o m e a u 提出了一个对称度更高的正六边形的l g a 模 型( 即f h p 模型 码。他们用此模型成功的模拟了一些典型的流体运动问题，并证明 了该模型的宏观行为基本上符合流体力学的n a v i e r s t o k e s 方程。 f h p 模型是将流场划分为间距为l 的正三角形网格，每个结点有六个邻居，正 4 华中科技大学硕士学位论文 好是一个正六边形的六个顶点。因此，f h p 模型又叫做正六边形格子气模型。结 点上的粒子可以具有六个不同的速度中的一个或者静止。 根据结点上有无静止的粒子，将f h p 模型分为三类：1 无静止粒子；2 一 一允许至多一个静止粒子；3 允许多个几个静止粒子。f h p 模型的碰撞规则设 计主要考虑到粒子数守恒、动量守恒。在保持粒子数守恒、动量守恒的前提下， 碰撞规则有多种设计方案。当一个入射状态对应两个出射状态时，按等概率选出 其中的一个出射状态 8 1 。 l g a 模型在微观层次上，从非平衡统计力学的观点出发，将流体视为由大量 的微观粒子所组成。这些粒子的运动时局部的，它们遵守力学定律、服从统计规 律，同时也能和宏观层次上的模型沟通。l g a 具有：l 绝对的数值稳定性；2 边界 处理简单；3 程序容易实现、并行性强等优点。然而，1 存在高统计噪声：2 存在 一些非物理的影响( 如，虚拟的不变量、不满足伽利略不变性谗 其不足之处。 1 2 3 格子b o l t z m a n n 方法 l b m 是为了保留l g a 的优点，尽可能地克服l g a 地缺点而提出来的。1 9 8 8 年，m c n a m a r a 和z a n e t t i 从分子混沌的假设出发，提出了格子b o l t z m a n n 方程模 型9 1 ，有效的消除了l g a 模型中的统计噪声，但该模型使用了规模很大的碰撞矩 阵，其计算量大且仍然不满足伽利略不变性。 1 9 9 2 年，q i a n 、c h e n 等人分别独立地采用b h a t n a g a r - g r o s s k r o o k ( 简记b g k ) 碰撞松弛模型1 0 1 提出了格子b g k ( l a t t i c eb g k ，简记l b g k ) 模型i u 2 ，将复杂的 碰撞操作转化成一个简单地松弛过程： ，( x + e 。，t + 1 ) = ( j 一) z ( x ，t ) + 吐矿唧( x ，t ) 其中，f a x ，t ) 是定义在离散速度集e j 上，时刻t 位于空间网格点x 处粒子密度 分布函数：，”( x ，t ) 是由系统当前的宏观量构造出来的局部动态平衡分布函数；松 弛因子6 0 与物性系数( 如：流体的粘性系数、热传导系数、质量扩散系数1 有关。微 观的粒子密度分布函数， ，t ) 与宏观密度p ( x ，t ) 、宏观速度u ( x ，t ) 的关系为： 华中科技大学硕士学位论文 z ( x ，t ) = p ( x ，t ) ，q ，( x ，t ) = p ( x ，t ) u ( x ，t ) 在l b g k 模型中，宏观的系统状态由微观粒子团的状态构成：反过来，微观 粒子团的运动又受控于系统动态的局部宏观状态。这体现了个体局部的自律性和 自适应性。局部动态平衡分布函数f e q ( x ，t ) 犹如一位“向导”动态地将系统引向符 合客观规律的状态。这种松弛技术还使得计算效率非常高。 自l b g k 模型提出以来，l b g k 的理论和应用都得到了较大的发展f 1 3 - 15 1 。一 些基本问题得到了解决或有了进展1 6 - 2 6 】；应用领域不断扩展【2 7 _ 2 9 。l b g k 模型的 主要优点有：l 与传统的数值计算方法相比计算效率高，程序容易实现、并行性强； 2 边界处理简单；3 无统计噪声；4 模型的宏观行为符合所研究的宏观方程( 如，流 体力学中n a v i e r - s t o k e s 方程等) 。 建立l b g k 模型的一般步骤是：l 选定离散速度集 ；2 根据所讨论的宏观 控制方程选定平衡态分布函数的类型( 平衡态分布函数仅依赖于宏观量) ；3 运用多 尺度展开以及t a y l o r 展开，根据控制方程确定松弛因子及其他待定系数。 l b g k 模型的计算步骤：1 选定网格大小确定血，选定出，根据所讨论的问 题模型，确定松弛因子甜；2 根据演化方程进行迭代计算；3 根据一定的要求结束 运算，其中根据演化方程进行迭代计算又细分为；( i ) 根据宏观初值或边值条件， 计算初始的平衡态分布函数，”( x ，t ) 作为分布函数，( x ，t ) 的初值；( i i ) 选用合适的 边界处理方法( 简单的有周期边界条件、反弹边界处理、边界外推插值等) ；( i i i l 计 算宏观量、计算平衡态分布函数、根据演化方程计算下一时刻的分布函数数值。 1 2 4 格子b o l t z m a n n 方法的应用 ( 1 ) 多孔介质流动 格子b o l t z m a r m 方法作为一种有效的计算方法目前在多孔介质流动方面得到 了比较成功的应用。与多孔介质流动的尺度相对应，利用格子b o l t z r a a n n 方法研究 多孔介质流动的方法有两类：孔隙尺度模拟和p , j g v 尺度模拟。早在1 9 8 8 年， 6 华中科技大学硕士学位论文 r o t h m a n 就利用格子气自动机来模拟多孔介质内流体流动，对d a r c y 定律进行了验 证 30 1 。s u c c i 等人利用格子b o l t z m a n n 方法对一个三维随机介质的渗透率进行了测 量，证实了d a r c y 定律1 3 1 1 oa d r o v e r 和g i o n a 研究了高空隙率情况下的二维随机介 质流体流动，测定了k o z e n y - - c a r m e n 公式中经验常数【3 “。此后，s i n g h 和m o h a n t y 对一般情况下的三维随机介质流动进行了模拟，研究了空隙率和空间相关性对渗 透率的影响【3 3 】。k i m 的研究小组研究了d a r c y 定律的修正形式t 3 4 1 。d a r d i s 和 m c c l o s k e y 直接在标准格子b o l t z m a n n e 中增加一阻力项，用以反映反弹格式的影 响【3 5 】。s p a i d 和p h e l a n 基于b r i n k m a n 模型，提出了一类模拟多孔介质流动的格子 b o l t z m a r m 模型 3 6 , 3 7 。后来，m a r t y s 对该模型进行了改进 3 8 1 。他们的计算结果表明 这些r e v 尺度模型的确可以描述宏观的多孔介质流动。 ( 2 ) 湍流的模拟 从数值计算角度看，格子b o l t z m a r m 方法是可以看作是求解n a v i e r s t o k e s 方程的一种数值方法，因此也能够用于湍流的直接模拟。事实上，早在1 9 9 0 年， b e n z i 和s u c c i 就对二维剪切湍流进行了模拟p ，并与谱方法的结果进行了对比， 使用的各自的规模最大达到5 1 2 x5 1 2 。他们发现格子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 方法的模拟结果仍然反映了二维湍流 的基本现象。s u c c i 等人使用b g k 模型和强化碰撞模型模拟了二维剪切湍流【4 刚， 比较了二者的差别。当r e y r l o l d s 数较小的时候，两种模拟结果的能谱基本相同， 但是强化碰撞模型的小尺度动力学的噪声较大。他们还发现，格子b o l t z m a n 模型 的数值稳定性和适应性也优于强化碰撞模型。c h e n 等人也做了相似的验证，将 利用b g k 模型和谱方法对三维各向同性湍流进行了模拟。在其研究工作中，他们 模拟了三维的b e l t r a m i 流动、衰减t a l o r - - g r e e n 涡和衰减的三维湍流。两科t 速度 场和漩涡场的空间和时间分部较好。这个评价得到了t r e v i n o 和h i g u e r a 的支持。 他们用伪谱法和格子b o l t a m a n n 方法研究了不同r e 数下的k o l m o g o r o v 流动的非 线性稳定性。为了进一步验证格子b o l t z m a r m 方法对湍流的模拟能力，m a r t i n e z 等 人用谱方法和格子b o l t z m a n n 方法研究了衰减剪切湍流【4 2 1 ，模拟的初始r e y n o l d s 7 华中科技大学硕士学位论文 数高达1 0 0 0 0 。格子b o l t z r n a r m 方法使用5 1 2 5 1 2 的格子，运行了8 0 个涡旋转时 间。他们仔细地比较了两种模拟方法的流函数和涡量场的空间分布和时间演化过 程，也研究了能量谱的时间演化。他们的结论是，格子b o l t z m a n n 方法模拟的总体 结果是精确的，波谱和涡量场与谱方法的结果都非常相似。格子b o l t z m a n n 方法同 样也可用于模拟非均匀各向同性的湍流。例如，b e n z i 等人利用格子b o l t z m a n n 模 型模拟了三维的剪切流动 4 3 1 ，用以验证各向异性流动中的自相似性。基于对流动 的分析，他们认为广义扩展自相似性对各向异性流动是成立的。a m a t i 等人则用格 子b o l t z m a n n 方法研究了湍流管道流【4 4 4 5 1 ，e g g e l s 也利用基于f c h c 格子的 b o l t z m a n n 方法模拟了三维湍流管道流 4 6 1 。 ( 3 ) 多相流 由于格子b o l t z m a n n 方法的微观粒子特性，它可以方便的描述不同相之间的相 互作用，使之在模拟这类复杂流动问题上具有常规方法所没有的优势。从而提供 了研究多相流系统的一种有效途径。按照设计方法，模拟多耜流的格子b o i t z m a n n 方法有着色模型、伪势模型、自由能模型和其他模型。1 9 8 8 年，r o t h m a n 和k e l l e r 提出了两相格子气自动机模型【4 7j 。g u n s t e n s e n 等人在此基础上提出了多相格子 b o k z m a n n 模型【4 8 1 ，g r u n a u 等人将此模型推广到密度和粘性变化的多相流系统4 9 1 。 k o n o 等将相变过程引入到这类模型中，并模拟了二维相变问题f 5 。】。这类模型的缺 点是在相界面附近会产生非物理现象，并且“颜色能量”的极小化过程的计算量 比较大。d a r t o n a 等对这类模型进行了改进，使计算效率得以提高。另外一种处 理相间相互作用的方法是基于均场理论的伪势方法，即通过一个伪势函数反映不 同相间的斥力或引力，从而导出非理想状态方程。伪势函数可以是指数型的( 旺1 ， 也可以是分数型的田 。上述两种模型都是基于界面动力学的唯象模型，并且只可 用于等温系统。s w i f t 等人通过修改平衡态分布函数，提出了基于自由能函数的格 子b o l t z m a n n 模型i s 4 1 。该模型与热力学理论一致，且满足局部动量守恒。虽然自 由能模型考虑了热力学效应，但是不满足伽利略不变性。最进，i n a m u r o 等人对此 模型作了改进口”，使伽利略不变性得以满足。上述各类模型中，湿度只是一个不 变的参数，并且都或多或少地与热力学理论矛盾 5 6 , 5 7 。最近的研究表明，可以直 8 华中科技大学硕士学位论文 接从绸密流体的动力学方程构造多相流格子b o l t z m a r m 模型陋5 9 。l u o 直接从绸密 气体的e n s k o g 方程出发，导出了一个与热力学相容的模型 5 6 , 5 7 ；h e 等人基于动力 学方程提出了一种模拟非理想气体的模型【5 8 1 ：从连续b o l t z m a n n 方程出发，h e 等 人还设计了类基于指标函数的不可压多相流模型【5 9 】，并模拟了二维和三维 r a y l e i g h - - t a y l o r 不稳定现象l ”“j 。 ( 4 ) 粒子悬浮流 为了模拟悬浮与流体中的颗粒，需要对流体一颗粒间相互作用进行处理。l a d d 在文献 6 2 , 6 3 1 中的开创性工作为使用格子b o l t z m a r m 方法模拟这类流动奠定了基础。 在l a d d 的模型中，颗粒内部和外部的流动可以穿透颗粒的边界进出颗粒，因而在 颗粒运动的过程中颗粒内部的流体质量并不守恒。a i d u n 和l u 修正了边界处理方 法【6 ”，使颗粒内流体满足质量守恒条件。此后人们还提出了消除或降低颗粒内部 流体影响的其他方法 6 5 - 6 7 】。目前人们已经用l a d d 的模型模拟了二维的和三维的球 形颗粒、椭球形颗粒、柱状颗粒的悬浮流动k 5 。6 8 1 。与实验结果和其他数值模拟结 果的比较表明这类模型是可行的。上述格子b o l t z m a n n 模型不考虑颗粒的布朗运 动。l a d d 在1 9 9 3 年构造了一类波动格子b o l t z m a n n 模型可以模拟这类运动【6 9 1 。在 该模型中，分布函数在演化过程中受一个随机波动的影响。该随机波动满足质量 和动量守恒，但对流体的应力张量有影响。s e g r e 等比较了这一模型的模拟结果和 实验结果，发现二者定量相符【70 】。l a d d 使用该模型模拟了3 2 0 0 个颗粒的悬浮问题 7 1 ，其结果提供了在实验中观察到但缺乏理论依据的一些现象。 除了前面介绍的几个方面外，格子b o l t z m a n n 方法的应用目前已经涉及到生物 流体、磁流体、交通流、燃烧、微尺度流动、化工、图像处理、量子力学、纳米 流体以及光学、声学等众多领域，并且这种向其它领域渗透的趋势越来越快。 1 3 本文研究工作介绍 本文对求解含源项的对流一扩散方程原有的格子b o r z m a n n 模型进行了改进， 消去了原来模型中对源项做的假设，使得模型的适用范围更大。在第二章中，将 给出三种不同的改进模型具体推导，三种模型的出发点都是为了消去了原来模型 9 华中科技大学硕士学位论文 中对源项做的假设。在第三章中，我们将使用第种改进模型对对流一扩散方程、 反应一扩散方程、含源项的对流扩散方程进行模拟，并和原有的模型模拟的结果 进行比较，这些模拟都是在二维空间中进行的：我们还将使用第二种改进模型对 三维的含源项的对流一扩散方程进行模拟。另外我们还将在第四章中使用改进模 型l 刺c i m a 模型进行模拟，并给出了二维空间中图灵班图的形成过程。 华中科技大字硕士学位论文 2 求解含源项对流扩散方程的l b g k 模型 2 。1 含源项对流扩散方程的介绍 含源项对流扩散方程的一般形式如下： o p ：( x , 一o 十“v p ( x , t ) = 口v 2 户( 工，f ) + f ( x ，f ) ( 2 1 1 ) 0 f 其中4 = ( u ，u 。，u n ) 1 是一个”维的常向量，代表对流速度；口是扩散系数，是一 常数；x = ( x ，x ：，x 。) 是空f 刮的坐标，r 代表时间：户( x ，f ) 是是物质在空削点x 和时刻，的密度；v 是一个关于x 的梯度算子；v 2 是一个关于x 的拉普拉斯箨予； f ( x ，f ) 是在空间点x 和时刻，的源项。 当脚= 0 时，方程( 2 1 ) 就变为含源项的扩散方程；多个这样的方程联立在一 起就形成了反应扩散方程；当f ( x ，r ) = 0 时，方程( 2 1 ) 就演变为对流扩散方程； 当“= 0 ，f ( x ，r ) = 0 时，方程( 2 1 ) 就变成了纯扩散方程。 2 2 求解含源项对流一扩散方程的格子b g k 模型的一些研究工作介 绍 1 9 9 5 年w o l f - - g l a d r o w 研究了高维空间中纯扩散问题的格子b o l t z m a n n 模型 7 2 , 7 3 ；1 9 9 9 年刘慕仁等人给出了求解二维古源项对流一扩散方程的l b g k 模型1 7 4 】； 同样是1 9 9 9 年，郭照立等人提出了空间任意维数对流扩散方程的通用l b o k 模 型7 5 j ；d a w s o n 等人于1 9 9 3 年提出了计算反应一扩散系统的格子b o l t z m a n n 模型 【7 6 】，并研究了s e l k o v 反应一扩散模型中的不动点、极限环以及斑图形成问题，他 们在模型中假设源项f ( 工，) 。ca t 2 ，其中出是时间步长。b l a a k 等人则研究了l b g k 模型模拟反应扩散系统时格子对称性于尺寸对模拟结果的影响。2 0 0 1 年李青 华中科技大学硕士学位论父 【7 8 1 等人利用了d a w s o n 的模型研究了c i m a 反应扩散系统中斑图的形成，他们指 出，模型需要假设f ( x ，) za t 。阎广武| 7 9 1 于2 0 0 1 年提出了- - , e e 新的l b o k 模型 来求解l o r e n z 方程，他的模型中也有假设f ( x ，f ) 。c a t ，而且他还指出这个假设是 非常重要的，在目前的所有工作中，这种假设无法消去。 2 3 求解含源项对流一扩散方程的格子b g k 模型 本节将针对原有求解含源项对流一扩散方程的格子b g k 模型都有对源项的不 同假设而提出三个改进模型，这三个模型消除了对源项的假设。 2 3 1 改进模型i 不币中米用的溟化方程为： ，j ( x + q a f ，+ r ) 一，( x ，f ) = 一当( ，( x ，r ) 一f ”( x ，f ) ) + f f ( x ，r ) + i _ a f 2 掣， z z o t 0 i q - i ，f 2 3 ，1 ，1 ) 不同于一般的格子b o l t z m a n 模型的演化方程，我们在演化方程的右边加上了 ；r8 f i 。( x ，, t ) 。演化方程中，血和m 分别是空间格子氏度和时间步长，c ：血出 z删 为粒子迁移的速率， c ，= c e ，是粒子运动的速度，e = 缸，0 兰i q - 1 是离散速度集 合：f 是松弛时间，其中f = ；f e q ( x ，f ) 是平衡态分布函数，这里我们取平衡态 函数为如下形式： 舷归州h 半+ 等一黔 眨，小：， c iz c iz c 这里q 是声速，由所取的格子模型决定；权系数q 满足q = 1 。 z ( x ，) 和，”( x ，f ) 满足： ；、f = 2 p ( 2 3 1 3 ) 1 2 华中科技大学硕士学但论文 q 护= p u ( 2 3 ，1 4 ) c ，q ，”= p u u + c ；p i ( 2 3 1 5 ) 演化方程中f 取为： 嘲用+ 半半) ( 2 3 1 6 ) 下面推导松弛时间f 和扩散系数口之间的关系。首先对演化方程中分布函数， 采用多尺度技术： ：= ，州+ 够1 + 占2 ：2 + 这里我们取占= a t 。上式两边对i 求和，然后利用方程( 2 3 1 3 ) ，我们可以得到： = 0 ，k 1 ( 2 3 1 7 ) 对演化方程应用多尺度展开和t a y l o r 展开，就可以得到： d f f , q + a t d ”+ 擎肛一妙) + t z m + 等等+ 0 ( 舻) ( 2 3 瑚) 其中皿2 昙+ e v 。对上述方程取出g o 一阶量，消去缸的1 一d r y2 一阶 量，我们得到 d j , ”= 一圭，。+ f + o ( f ) ( 2 3 1 9 ) 将方程( 2 3 1 9 ) 应用到方程( 2 3 ，1 8 ) 中，我们可以将( 2 3 1 8 ) 重写为： d ，z + a t ( 1 一扣肚a _ 坐2 td 扩e 一抄蟛锄m + 等等+ 0 ( f 2 ) 方程两边同时消去竽墨后得到： zo t q ，”+ & ( 1 一去) d j ，m + 了a t c ，- v f = 一( z m + 斫2 ) + 只+ d ( 血2 ) ( 2 3 1 1 0 ) 方程( 2 3 1 9 ) 两边对i 求和，然后应用方程( 2 3 i 3 ) ，( 2 ，3 1 4 ) ，( 2 3 1 5 ) ，( 2 ，3 1 7 ) ， 我们可以得到： 竺栅即= e + o ( a r ) (2肌11)3t 。j j ， 、。1 华中科技大学硕士学位论文 方程( 231 1 0 ) 两边对i 求和，然后应用方程( 2 3 ，13 ) ，( 2 3 1 ，4 ) ，( 2 3 1 5 ) ，( 2 3 1 7 ) ， 我们可以得到： o d p ，+ u v p + a t ( 1 一去) ；口，+ 了a t v ( ；t f ) = ( + 。( ) ( 2 3 川2 ) 应用方程( 2 3 1 9 ) ，( 23 13 ) ，( 231 4 ) ，( 2 3 ，15 ) ，( 2 3 1 7 ) 和( 2 3 1 1 1 ) ，有 莩口2 鲁军+ 甲( 军q ) = 。一刃；呱p ，“一f ) + 。( 出) 一刃( 莩昙c q r 。+ v t 莩c 一莩c 。只 + 。e 力 = 一刃( 昙( p h ) + v ( p “ + c ；p ，) 一；c f + 。( r ) = 一r c ! v 2 p - v ( 咧昙p 协跏) ) 舟( r c 朋+ 。( f ) 一2 c j v 2 p v ( r u f ) + v o c ，只) + o ( a t )( 2 3 1 1 3 ) 将方程( 23 1 1 3 ) 代入方程( 2 3 1 1 2 ) ，得到： o 洲p + u v p + v ( 折莩。，f ) 冉( “1 2 - r ) u ，2 f , ) = c 妒1 2 ) 掷2 p + z ，卢：+ 。( f 2 ) ( 2 3 1 1 4 ) 由f 的定义，有： 只= f 孕巧= 孚f 将上面的式子代x 2 h - i ；呈( 2 3 1 1 4 ) ，我们得到 詈协即刮卜i ) a t v p + f 蜊) 将上述方程和方程( 2 1 1 ) 比较，我们得到了松弛时间r 和扩散系数口之间的 关系： 华中科技大学硕五学位论文 2 ，3 2 改进模型i i a = ( r 一告) 血 本节中采用的演化方程为： 似托，h ) 一似，忙一号( 似力吖杈朋) + 幔( 圳) + j 1 f 2 口f ( 列) ， 0 i q 一1 ， ( 2 3 2 ，1 ) 不同于一般的格子b o l t z m a n 模型的演化方程，我们在演化方程的右边加上了 去础2 d f ( x ，f ) 。演化方程中，相关的量的意义和2 3 1 节中的一样，”( x ，r ) 是平 衡态分布函数，这里我们取平衡态函数为如下形式： f q ( x , t ) = c o p ( 1 + 半+ 等一蜘 皿，2 z ， 这里c s 是声速，由所取的格子模型决定；权g g z o , ，满足q = 1 。 a ( x ，f ) 和”( x ，f ) 满足： z = z ”= 尸( 2 3 2 3 ) c f 矿= p u( 2 324 ) c ，c ；矿= p u u + c ；p l( 2 3 2 5 ) 演化方程中f 取为： f = q f ( 1 + 掣)( 2 3 c ： 、 。 下面推导松弛时间r 和扩散系数口之间的关系。首先对演化方程中分布函数z 采用多尺度技术： ：= ，呵+ 锐+ 占2 ：2 + 这里我们取5 = a t 。上式两边对i 求和，然后利用方程( 2 3 2 3 ) ，我们可以得到： 华中科技大学硕士学位论文 = o ，t 1( 2 3 2 7 ) 对演化方程应用多尺度展开和t a y l o r 展开，就可以得到： 。，f