（应用数学专业论文）若干非线性偏微分方程的格子bgk模拟.pdf

中文摘要 格子b o l t z m a 丑n 方法( l b m ) 是一种新兴的模拟流体和复杂物理系统的数值计算方 法。不像基于宏观连续方程的传统数值方法，l b m 是起源于微观模型和细观运动论的介 观方法，它具有许多分子动力学的优点，如物理图像清晰、容易处理复杂边界、编程容 易实现等。近年来，l b m 在模拟线性和非线性偏微分方程方面取得了重要进展，但是理 论部分仍有许多问题有待完善，例如如何构造出精度较高的模型和如何模拟更复杂的非 线性偏微分方程。 本文首先在绪论部分简要介绍了l b m 的发展历史及其应用，然后在接下来的 四章中分别针对几类非线性偏微分方程，利用多尺度技术，建立了相对应的几 种格子b g k 模型。第一章中针对二维对流扩散方程建立d 2 q 4 模型；第二章中针 对s i n e - g o r d o n 方程建立隐式格子b o l t z m a n n 模型；第三章中针对广义k d v 方程，k d v - b u r g e r s 方程，组合k d v - m k d v 方程和广义b u r g e r s - h u x l e y 方程建立统一的具有五阶 精度的格子b g k 模型；第四章中针对广义k u r a m o t o - s i v a s h i n s k y 方程建立d 1 q 5 格 子b g k 模型，由于宏观方程含有四阶导数，标准的l b m 无法恢复出来，因此，本模 型的提出填补了这一方面的空白，拓展了l b m 在模拟复杂非线性偏微分方程方面的领 域。 数值结果表明所建模型均十分有效，为以后更复杂和更高维的复杂非线性偏微分方 程的数值模拟积累了经验。 关键词：格子b o l t z m a n n 方法，c h a p m a n e n s k o g 展开，非线性偏微分方程， 对流扩散方程，s i n e - g o r d o n 方程，广义k u r a m o t o - s i v a s h i n s k y 方程，k d v - b u r g e r s - k u r a m o t o 方程 中文文摘 非线性偏微分方程( n p d e s ) 在数学、物理学和化学等不同领域里扮演着非常重要 的角色。对这些非线性方程的研究成为广大物理学、力学、地球科学、生命科学、应用 数学和工程技术等科学工作者的重要课题。在物理系统中，许多很有趣且有用的特性隐 含在它们的非线性行为里，如果它们的非线性方程的解析解存在，则通过这些解析解我 们可以更好地了解其复杂物理现象和动力学过程的机理。 虽然人们已建立和发展了一些求解非线性偏微分方程的方法，但对于非线性演化方 程目前仍无统一的求解方法。一般的非线性偏微分方程只在特殊的初边值条件下才有解 析解，而在通常情况下只能进行数值求解。大部分n p d e s 的研究是用近似的数值方法 来处理方程中的非线性项。许多关于这些n p d e s 的数值模拟方法发展了起来，包括有 限差方法、热平衡积分法、有限元法、谱方法、变分迭代法等。 作为一种新兴的数值方法，格子b o l t z m a n n 方法( l b m ) 不同于传统的数值方法，它 是基于微观模型和细观运动论的介观方法。l b m 在求解非线性方程以及复杂系统的演 化，特别在流体力学的研究中取得了很大成果。这是由于l b m 具有物理图像清晰、边 界处理容易、编程实现简单等优点。从计算量的层次面来看，l b m 属于显示时间推进 方法，每个时间步的计算量为o ( m n ) ( m 为离散速度数，n 为计算格点数) ，其计算效 率要高于一般的数值方法。由于模型所涉及的计算都是具有局部性，所需局部平衡态分 布函数是同时进行计算的，具有天然的本质并行性，非常适合在大规模并行计算机上运 行。再加上“迁移 步在实际编程计算中只是一个赋值过程，并不占计算时间，所以其 效率是比较高的，这也是这个方法蓬勃发展的原因之一。l b m 提供了联系宏观和微观的 可能性和现实性，除了在一般的流体力学问题中得到了成功的验证之外，在湍流、多相 流、多组分流、粒子悬浮流、量子力学以及磁流体力学等相关领域也具有广阔的应用前 景。 近2 0 年来，l b m 在某些复杂的演化方程模拟中取得了重要的进展，这些复杂系 统包括如对流扩散方程、反应扩散方程、b u r g e r s 方程、m k d v 方程，k d v - b u r g e r s 方 程，s c h r s d i n g e r 方程等。但是理论部分仍有许多问题有待完善，例如如何构造出精度较 高的模型和如何模拟更复杂的非线性偏微分方程。本文针对各类非线性偏微分方程建立 了比较有效的格子b g k 模型，提高原有模型的精度，扩展了l b m 模拟非线性偏微分方 程的应用领域。 i i i 福建师范大学赖惠林硕士学位论文 我们首先在绪论中简要介绍了l b m 的发展历史及其应用。 在第一章中，我们针对二维对流扩散方程，建立了带修正项的d 2 q 4 的格子b g k 模 型，并用数值算例进行了验证比较。与之前的七速和九速格子模型相比，在同样精度的 情况下，我们的模型少了几个速度数，这样在计算量方面就节省了许多。特别是第一个 算例，我们可以发现经过相当长时间的演化，当t = 5 0 0 0 0 时，模拟解与解析解仍然相 当吻合，可见本模型的数值有效性和稳定性。 在第二章中，我们针对一维s i n e - g o r d o n 方程，建立了一维隐式格子b o l t z m a n n 模 型，此模型可以很容易转换成显格式，其中的自由参数口提高了模拟s i n e - g o r d o n 方程 的精度和灵活性。而且相比之前的显格式模型，此模型在数学推导上更为简洁。根据数 值模拟结果，我们可以得知数值解随着时间的变化而呈周期性变化，即具有孤立波和周 期波的特点，这有效地支持了b r e z i s 的理论分析。 在第三章中，我们针对一大类非线性偏微分方程，在现有模型的基础上，建立了 统一的一维带修正项d 1 q 5 的高阶格子b g k 模型，可以有效模拟广义k d v 方程，k d v - b u r g e m 方程，组合k d v - m k d v 方程，广义b u r g e r s - h u x l e y 方程等，不仅提高了模型精 度，使其具有五阶精度，而且其中的参数入与，7 的选择保留了原先模型的灵活性。本模 型的提出为以后高阶格子b g k 模型的建立积累了经验。 在第四章中，我们针对广义k u r a m o t o - s i v a s h i n s k y ( k - s ) 方程，建立了一维带修正 项的d 1 q 5 格子b g k 模型，通过增加一个修正函数，可以正确地恢复出含有四阶导数 的广义k s 方程，由于广义k s 方程中含有四阶导数，而标准l b m 无法恢复出来的， 因此，本模型的提出填补了l b m 对于宏观方程具有四阶导数的数值模拟这一方面的空 白，使得l b m 的应用领域扩展到含四阶导数的非线性偏微分方程。这为以后模拟更复 杂的非线性偏微分方程和更高维数的复杂系统提供经验。 最后我们对全文作了总结和展望。 i v a b s t r a c t l a t t i c eb o l t z m a n nm e t h o d ( l b m ) i san e wt e c h n i q u ef o rs i m u l a t i n gf l u i d sa n dr o o d e l i n gc o m p l e xp h y s i c si nf l u i d s u n l i k ec o n v e n t i o n a ln u m e r i c a lm e t h o d sb a s e do nam a c r o - s c o p i cc o n t i n u u me q u a t i o n ，t h el b ms t a r t sf r o mm i c r o s c o p i cm o d e l sa n dm e s o s c o p i ck i n e t i ce q u a t i o n s i tp r o v i d e sm a n yo ft h ea d v a n t a g e so fm o l e c u l a rd y n a m i c s ，i n c l u d i n gc l e a r p h y s i c a lp i c t u r e s ，e a s ei ni n c o r p o r a t i n gc o m p l e xb o u n d a r yc o n d i t i o n sa n ds i m p l i c i t yo f p r o g r a m m i n g r e c e n t l y , t h el b mh a v eb e e nd e v e l o p e dt os i m u l a t el i n e a ra n dn o n l i n e a r p a r t i a ld i f f e r e n t i a le q u a t i o n s ( n p d e s ) h o w e v e r ，t h e r ei sat r o u b l e s o m ep r o b l e mf o rs o l v - i n gn p d e 8 i nm a n ye x i s t i n gl a t t i c eb o l t z m a n nm o d e l s ，i e ，h o wt od e r i v eh i g h e ra c c u r a c y a n dm o r ec o m p l e xn o n l i n e a rt e r m si nn p d e s a tt h eb e g i n n i n go ft h i st h e s i s ，s o m eb a s i ss u m m a r i e so ft h el b ma r ei n t r o d u c e d i n c h a p t e 1 ，w ec o n s t r u c tad 2 q 4m o d e lf o rs i m u l a t i n gt w o - d i m e n s i o n a lc 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 n ；i nc h a p t e r2 ，w ee s t a b l i s ha ni m p l i c i ts c h e m eo fl a t t i c eb g k m o d e lf o rs i m u l a t i n g s i n e g o r d o ne q u a t i o n ；i nc h a p t e r3 ，w ep r o p o s eau n i f i e dl a t t i c eb g km o d e l ，w h i c h h a sf i v eo r d e ra c c u r a c y , f o rs i m u l a t i n gg e n e r a l i z e dk ( e q u a t i o n ，k d v - b u r g e r se q u a t i o n ， c o m b i n e dk d v - m k d ve q u a t i o na n dg e n e r a l i z e db u r g e r s - h u x l e ye q u a t i o n ；i nc h a p t e r4 ， w ep r o p o s ead i q 5l a t t i c eb o l t z m a n nm o d e lf o rt h eg e n e r a l i z e dk u r a m o t o - s i v a s h i n s k y ( g k s ) 吾q u a t i o nw h i c hh a sf o u r - o r d e rd e r i v a t i v e d u et ot h ef o u r - o r d e rd e r i v a t i v ee x i s t s - i nt h eg k se q u a t i o n 一t h e nt h e s t a n d a r d l b m c a n n o t r e c o v e r t h e g o v e r n i n g 倒o l u t i o n _ _ e q u a t i o n s o u rp r e s e n t e dm e t h o d 丘ut h i sg a pa n de x t e n dt h ea p p l i c a t i o no ft h el b mt o t h eh i g h e rd e r i v a t i v et e r m sp r o b l e m f r o mt h es i m u l a t i o n s ，w ef i n dt h a tt h es i m u l a t i n gr e s u l t sa r ei ne x c e l l e n ta g r e e m e n t w i t ht h ea n a l y t i c a ls o l u t i o n s t h i si n d i c a t e st h a tt h ep r e s e n tm e t h o d sa r ee f f i c i e n ta n d f l e x i b l ea p p r o a c h sf o rp r a c t i c a la p p l i c a t i o n t h ep r e s e n tm o d e l sc a nb ee x t e n d e dt om o r e c o m p l e xs y s t e m s k e yw o r d s ：l a t t i c eb o t z m a n nm e t h o d ，c h a p m a n - e n s k o ge x p a n s i o n ，n o n l i n e a rp a r - t i a ld i f f e r e n t i a le q u a t i o n 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 n ，s i n e - g o r d o ne q u a t i o n ，g e n e r a l i z e d k u r a m o t o - s i v a s h i n s k ye q u a t i o n ，k d v - b u r g e r s - k u r a m o t oe q u a t i o n i i 福建师范大学硕士学位论文独创性和使用授权声明 专业应用数学所呈交的 论文( 论文题目：若干非线性偏微分方程的格子b g k 模拟) 是本人在导 师指导下，独立进行的研究工作及取得的研究成果。尽我所知，除论文中 已特别标明引用和致谢的内容外，本论文不包含任何其他个人或集体已经 发表或撰写过的研究成果。对本论文的研究工作做出贡献的个人或集体， 均已在论文中作了明确说明并表示谢意，由此产生的一切法律结果均由本 人承担。 本人完全了解福建师范大学有关保留、使用学位论文的规定，即：福建 师范大学有权保留学位论文( 含纸质版和电子版) ，并允许论文被查阅和 借阅；本人授权福建师范大学可以将本学位论文的全部或部分内容采用影 印、缩印或扫描等复制手段保存和汇编本学位论文，并按国家有关规定， 向有关部门或机构( 如国家图书馆、中国科学技术信息研究所等) 送交学 位论文( 含纸质版和电子版) 。 ( 保密的学位论文在解密后亦遵守本声明) 学位论文作者签名：始惠杯 签字日期：2 。c i 年月5 日 指导教师签名：善学钐q 签字日期1 年阳厢 绪论 当前，非线性偏微分方程( n p d e s ) 研究领域中的一个重要的研究方向就是运用偏 微分方程来研究物理、化学、生物和经济等领域中的非线性现象【l 】1 。现实生活的许多领 域内数学模型都可以用n p d e s 来描述，很多重要的物理、力学等学科的基本方程本身 就是n p d e s ，所以对n p d e s 的研究有着重要的现实意义【2 ，3 】。基于非线性方程的复杂 性，至今能够求出解析解的方程很少，大部分方程只能进行数值求解。因此，如何发展 比较有效的数值模拟方面就显得十分重要。传统的数值方法大多是基于宏观连续模型的 自顶向下的方法，如有限差分法【4 ，5 l 、有限元法【6 】、有限体积法【7 】、有限谱方法【8 】等，它 们是以非线性的微分方程出发，对微分方程进行离散得到代数方程组或者常微分方程系 统，然后再用标准的数值方法求解。近年来备受人们关注的格子玻耳兹曼方法( l a t t i c e b o l t z m a n nm e t h o d ，l b m ) 9 1 0 】发展了起来，与传统的数值方法不同，l b m 是基于微 观离散模型的自底向上的方法。在l b m 中，流体被抽象为大量的微观虚拟粒子，这些 粒子通过简单的方式在规则的离散网格上迁移和碰撞。通过对这些粒子进行统计平均， 就可以得到流体的宏观运动特性。与常规的数值方法相比，l b m 具有物理图形清晰、边 界处理简单、程序容易实现、并行性较好、模型健壮性较高等优点【1 1 】。因此，它被誉为 最具前途的数值模拟方法之一。 下面我们首先对l b m 的发展历史做一个简要回顾。 格子b o l t z m a n n 方法的发展历史 格子气自动机 格子b o l t z m a u n 方法( l a t t i c eb o l t z m a u nm e t h o d ，l b m ) 历史上源于1 9 7 0 年代提出 并发展的格子气自动机( l a t t i c eg a sa u t o m a t a ，l g a ) 【1 2 】方法，它的出发点是流体微观 离散的物理模型，可以看做是求解连续b o l t z m a n n 方程的一种特殊的离散格式。l g a 是一种计算流体的模型，它是更广泛的元胞自动机( c e l l u l a ra u t o m a t a ) 在流体力学上的 应用。第一个l g a 是由法国的h a r d y 、d ep a z z i s 和p o m e a u 在2 0 世纪7 0 年代提出的，并 以这三位作者的名字命名为h p p 模型【1 3 1 4 】，此模型是空间时间完全离散的数学模型， 它将流体视为离散的大量粒子，流场离散为一个规则正方形的格子，粒子按照一定的规 则在格子上进行碰撞和迁移，其运动速度方向定义如下： e l = c ( i ，o ) ，e l = c ( o ，1 ) ，e l = c ( - 1 ，o ) ，e 4 = c ( o ，- 1 ) ， 1 福建师范大学赖惠林硕士学位论文 这里c = 瓦瓦，瓦和况分别表示格子步长和时间步长。在l g a 中一般取如和况为长度 和时间单位，因此c = l 。这些单位称为格子单位( l a t t i c eu n i t ) 。在l g a 中一般还要求 粒子分布满足p a u l i 不相容原理，即在每个格点处以某速度的粒子最多只能有一个。因 此每个格点处的粒子分布状况可以用一个4 位的布尔变量表示，即n ( x ，t ) = 死l n 2 n 3 n 4 ， 其中m = 1 或者0 表示有或无以速度q 运动的粒子。h p p 模型的状态演化可以分为两 个阶段，即( 1 ) 迁移( s t r e a m i n go rp r o p a g a t i o n ) 过程：粒子从一个节点在一个时间步内， 以恒定的速度运动到相邻的节点；( 2 ) 碰撞( c o l l i s i o n ) 过程：在一个格点上从相邻节点运 动来的粒子发生碰撞，根据质量、动量和能量守恒规则改变粒子的速度，然后各个粒 子又以改变后的速度迁移。这两个步骤交替循环，直到流场达到稳定收敛。碰撞方式 是l g a 的核心部分，必须保证碰撞前后粒子的质量、动量和能量守恒。h p p 模型碰撞 规则如下：当两个速度相反的粒子到达同一个格点而另外两个方向上没有粒子时，发生 对头碰撞，即两个粒子的速度分别旋转9 0 “，而在其他情况下粒子速度不发生改变。见 图0 - 1 。 图0 - 1 ：h p p 正方形网格模型及其碰撞规则 t h el a t t i c es q u a r eg r i da n dc o l l i s i o nr u l eo ft h eh p pm o d e l h p p 的演化方程可以写为： 啦( z + e i 况，t + 瓦) = 啦( z ，t ) + q t ( n ( z ，t ) ) ，i = 1 ，2 ，3 ，4 其中q t 是碰撞算子。满足局部质量、动量和能量守恒： q = 。，e 她= 。，譬继= 。 2 绪论 h p p 演化方程也可以按照粒子运动的物理过程表示为两个步骤： 碰撞： 礼：( z ，t ) = 佗t ( z ，t ) + q i ( n ( z ，t ) ) 迁移： m ( z + e t & ，t + 6 t ) = 仃：( z ，t ) 虽然h p p 模型在微观上满足质量和动量守恒，但其宏观动力学方程并不能恢复 出n a v i e r - s t o k e s 方程的非线性项，也不能恢复耗散项，其原因是h p p 模型中的格子缺 乏足够的对称性，应力张量不能满足各向同性性，且存在伪随机守恒量，因此，为了得 到正确的宏观动力学方程，需要寻找具有足够对称性的格子模型。 1 9 8 6 年，法国的f r i s c h ，p o m e a n 和美国的h a s s l a c h e r 构造出了一个具有足够对称 性的二维l g a 模型，命名为f h p 模型【1 2 l 。该模型使用的格子是规则的正六边形，见 图0 - 2 ，流体粒子具有6 个离散速度： e ：cc o s 仇，s i n 仇) ，仇：垒二! 堕，i ：1 6 岛 岛 弓 岛 图0 - 2 ：f h p 正* 边形网格模型 t h el a t t i c er e g u l a rh e x a g o ng r i do ft h ef h pm o d e l f h p 的碰撞规则比h p p 更为复杂，碰撞类型有对头二体碰撞、对称三体碰撞、非 对称三体碰撞和对称四体碰撞。有些模型还引入静止粒子，此时还发生运动粒子与静止 粒子的碰撞。f h p 模型比h p p 模型有更高的对称性，碰撞过程更为丰富。f h p 在一定 条件下可以恢复到宏观的n a v i e r - s t o k e s 方程。 在l g a 中，粒子存在于离散的格子节点上，并沿着格线迁移。粒子的演化只涉及 到相邻的节点，因此l g a 方法具有物理图形清晰、计算稳定等优点，但是也存在一些 3 福建师范大学赖惠林硕士学位论文 缺点，如动量方程不满足伽利略不变性、局部量存在统计噪声、碰撞算子具有指数复杂 性等。 从格子气自动机到格子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 中直接使用布尔变量的统计平均量粒子分布函数进行演化f i s ，这是最早的格 子b o l t z m a n n 模型。由于该模型仍采用l g a 的碰撞方式，具有指数复杂性。 1 9 8 9 年，h i g u e r a 和j i m e n e z 对上述模型作了进一步简化，提出了线性碰撞算子模 型【1 8 1 ，该模型引入了平衡态分布函数，将碰撞算子用一个碰撞矩阵来代替，将其线性化 处理。这使得计算复杂性大大降低。但是数值稳定性仍较差。同年h i g u e r a 等进一步提 出强化碰撞算子方法【1 7 1 ，克服了指数复杂性，使得计算量和存储量大大降低。 1 9 9 1 1 9 9 2 年，c h e n 、q i a n 等人分别独立采用b h a t n a g a r - g r o s s - k r o o k ( b g k ) 碰撞 松弛模型提出了一种更简单的模型，即单松弛时间( s i n g l er e l a x a t i o nt i m e ，s r t ) 或格 子b g k ( l b g k ) 模型【1 8 t1 9 l ，将复杂的碰撞操作转化为一个简单的松弛过程。在该模型 中，碰撞矩阵由一时间松弛系数来确定，即碰撞算子为： q = 一圭( ，一，( 。) ) 则b o l t z m a n n 方程在不受外力的情况下可以近似表示为： 筹+ e v ，= 一 ( ，产，) 这就是b o l t z m a n n 方程的格子b g k 模型，其中下为驰豫时间，( 叼) 为局部平衡态分布 函数，即m a x w e l l b o l t z m a n n 分布函数。格子b c k 模型继承了格子气自动机中空间离 散、流体离散和时间离散的思想，引入碰撞间隔理论来处理碰撞项，使b o l t z m a n n 方程 在格子的基础上能够顺利的求解。在应用这一模型时关键是如何选择局部平衡态分布函 数，( 叼) 和驰豫时间下。局部平衡态分布函数，( 嘲的选取与所研究的对象有着密切的关 系，必须满足质量守恒、动量守恒和能量守恒。 在实际计算中，可假设上述方程沿特征方向有效，将粒子速度e 简化为有限维的速 度空间( e o ，e 1 ，e 2 ，e ) ，这边代表速度方向数。分布函数，和局部平衡态分布函 数，( 叼) 也做同样处理，则上述连续的b o l t z m a n n 方程可写成离散形式如下： 警怕v = 一季( 五一) ( _ l 2 ，) 4 绪论 根据质量和动量守恒规则，微观粒子分布函数五与宏观密度p 和宏观速度1 1 的关系 为： p = 五，胆= 五e i tt 上式为宏观变量和微观变量间相互转换的关系式。 我们可以根据模型，选择合适的局部平衡态分布函数矗咖，通过多尺度分析可恢 复出相应的宏观物理方程。所用b o l t z m a n n 方程实际上是一个线性偏微分方程，因 此容易在计算机上实现，在一定条件下可恢复出正确的n a v i e r - s t o k e s 方程，有效克服 了l g a 方法的不足。格子b g k 模型的提出使得格子b o l t z m a n n 方法的研究达到了一个 新的水平，并在后面的发展中逐渐成熟，成为目前应用最广泛的模型2 1 】。 建立格子b g k 模型的一般步骤是： ( 1 ) 针对宏观方程，选择适当的离散速度，例如对于一维模型，一般可以采 用d 1 q 2 ，d 1 q 3 ，d 1 q 5 ，d 1 q 7 等模型；二维模型可采用d 2 q 4 ，d 2 q 5 ，d 2 q 7 ，d 2 q 9 等模型；三维模型可采用d 3 q 1 5 ，d 3 q 1 9 等模型；当速度数增加时，所需平衡态分布函 数也需增加，这样就会增加了计算量，因此，在选择模型的时候可以根据需要选择适当 的速度数模型。 ( 2 ) 选择适当的格子b o l t z m a n n 方程，有时需要根据宏观方程进行修正，也就是增 加一修正项。然后利用泰勒展开和多尺度技术进行宏观方程的恢复，根据所讨论的宏观 方程确定适当的局部平衡态分布函数和所需要的参数表达式。 格子b g k 模型的计算步骤一般是： ( 1 ) 选择适当的网格大小，确定时间步长民和空间步长如，然后计算出模型所需要 的各个参数值。 ( 2 ) 选择合适的边界处理方法，大的方法包括启发式格式、动力学格式、外推格 式、复杂边界处理格式等。 ( 3 ) 对计算区域进行初始化，包括分布函数的初始化和平衡态分布函数的初始化。 ( 4 ) 根据演化方程进行迭代计算，其中包括碰撞和迁移。实际上迁移在计算中只是 一个赋值过程，并不占计算量。 ( 5 ) 计算宏观量，计算平衡态分布函数，根据演化方程计算下一时刻的分布函数。 格子b o l t z m a n n 方法的应用 与传统的宏观数值方法相比，l b m 是基于微观模型与细观运动论，具有许多常规 5 福建师范大学赖惠林硕士学位论文 数值方法没有的独特优点，如物理图形清晰，边界处理简单，易处理复杂区域，对于多 项流和多组分流，不需要跟踪不同相间的界面，更重要的是，它的本质并行性使得它 能够适用于并行计算。事实上，在2 0 年的发展过程中，l b m 已成为一个十分活跃和极 具发展前景的模拟方法，并迅速在不可压流【2 2 】，可压缩流卿7 1 、微管道流【2 删、多项 流【3 1 】，多孔介质测3 2 】，悬浮流【3 3 】，湍流叫，磁流体力学【1 8 ，3 5 1 ，生物流体删等相关领域 得到了比较成功的应用。 由于n p d e s 在物理学和数学等不同领域里扮演着非常重要的角色。在物理系统 中，许多很有趣且有用的特性隐含在它们的非线性行为里。如果它们的非线性方程的解 析解存在，则通过这些解析解我们可以更好地了解其复杂物理现象和动力学过程的机 理。因此，寻找和构造n p d e s 的解析解显得越来越重要【3 7 ，3 田。但是由于解析解一般只 存在于某些比较严格的条件下，大部分n p d e s 的研究是用近似的数值方法来处理方程 中的非线性项。 近年来，随着l b m 的发展，l b m 模拟一些线性和非线性的偏微分演化方程 取得了很大的成功，所做的工作包含了波动方程【3 9 】，反应扩散方程阳，4 1 】，对流 扩散方程【4 2 删，b u r g e r s 方程【4 1 1 ，k d v 方程 5 2 - 5 4 1 ，m k d v 方程【5 8 】，s c h r s d i n g e r 方 程 5 9 1 ，k d v - b u r g e r s 方程 e , o - e 2 l ，p o i s s o n 方程【6 3 】，r o s s l e r 方程酬等，取得了一系列成 果，扩展了l b m 在模拟非线性偏微分方程方面的应用领域。但是理论部分仍有许多问 题有待完善，例如如何构造出精度较高的模型和如何模拟更复杂的非线性偏微分方程。 最近一些学者在这方面做了许多重要的工作 5 4 ，6 2 ，弱】。 由于l b m 具有几何灵活、并行本性、易于编程、精度较高等特点，从它诞生之日 起，就倍受物理学家、数学家、计算机专家和其它领域的科学家的关注。在这十几年 里，它得到了充分的发展和广泛的应用。 6 第一章二维对流扩散方程的格子b g k 模拟 1 1引言 二维对流扩散方程是非线性偏微分方程中最重要的方程之一，在科学工程中有着 重要的应用，包括气体、地下水污染物的运输，油层流动等。有许多学者用不同的方 法研究了该方程，例如有限体积法【删、有限差方法【6 7 】、局部不连续g a l e r k i n 法删、不 连续h p 有限元法【6 9 】、外推加速法1 7 0 】、迎风有限体积元法【7 l 】、交替方向差分流线扩散 法 h i 、高阶紧致边界元法【7 3 】等。近年来许多学者尝试应用格子b o l t z m a n n 方法研究扩 散问题并取得了不少成果【州5 7 4 】。对流扩散方程是描述扩散现象的一个重要方程， 本章针对二维对流扩散方程，采用d 2 q 4 速度模型，构造了一类带修正项的b g k 型格 子b o l t z m a n n 方法【7 5 】。计算机模拟结果表明，这一方法是十分有效的。 1 2 格子b g k 模型 设qcr 2 是一个有界区域，z = ( z l ，z 2 ) q 。考虑如下初边值对流扩散方程 f 象幢v “= 山，z 呱 = 0 ，z f = a q ，亡 0 ， ( 1 2 1 ) i “( z ，0 ) = 咖( z ) ， 其中 _ b v = b 1 丢z i + 去，= 昌+ 岛， u 0 是密度， 0 是扩散系数，6 1 ，6 2 为对流系数 计算流体力学研究的格子b o l t z m a n n 方程，在模拟抛物型方程时，要求时间步长与 空间步长的平方相当。在此，本文用小参数s 来表示空间步长( e 可以看作平均自由程与 特征长度的比值) 。采用d 2 q 4 模型使用的离散速度方向如下 【e l 忍忍心】= f1 。 i l 010 1j 7 福建师范大学赖惠林硕士学位论文 图i - i ：d 2 q 4 模型 t h el a t t i c ef o u r - v e l o c i t ym o d e l 带修正项的格子b g k 演化方程为 五( z + c s e t ，亡+ e 2 t ) 一 ( z ，t ) = 一：1 f i c z , t ) 一并o ( z ，亡) + s ( z ，亡) ， ( 1 2 2 ) 其中五( 茁，亡) 和群o ( ，亡) 分别表示密度分布函数和平衡态分布函数，c h i ( x ，t ) 为修正 项，t 为无量纲的特征时间，c 为常数，下为驰豫时间，稳定性要求r 满足r 0 5 。 宏观变量让满足守恒律 u ( z ，亡) = 五( z ，亡) = ( z ，亡) ， ( 1 2 3 ) it 这样，通过选择适当的平衡态分布函数，我们可以恢复对应的宏观方程。 事实上，对方程( 1 2 2 ) 左边泰勒展开并保留三阶项 打警+ 色v + 三c 2 舻慨v ) 2 + 3 ) = 一 五一棚+ 积 ( 1 2 4 ) 为恢复宏观方程( 1 2 1 ) ，对分布函数进行c h a p m a n - e n s k o g 多尺度展开 五= 0 0 + e 鼻1 + e 2 2 + d 3 ) ， ( 1 2 5 ) 其中k n u d s e ne ( 可以取为时间步长a t ) 定义为：s = e l ，这里z 为平均自由程，l 为特征长度。其中非平衡态分布函数z ，矗，g = 1 ，2 ，3 ，4 ) 满足守恒律 z d = = 0 8 ( 1 2 6 ) 第一章二维对流扩散方程的格子b g k 模拟 将( 1 2 5 ) 代入( 1 2 4 ) ，并比较两端小参数e 的同阶项得到 0 ( e ) ： 。自v j c f 0 1 一& = 一享矗， 0 ( ，) ： 警+ 争v 矗1 ，+ 嘉( 州胛：一寺伊 将( 1 2 7 ) 代入( 1 2 8 ) 得 笔+ 争v & 一事一丢) ( e i v ) 2 f i ( 0 ，= 一寺并 对方程( 1 2 9 ) 两边关于t 求和，利用方程( 1 2 3 ) 和( 1 2 6 ) ，我们可以得到 筹+ 等莩( 川& ) 一事( r 一丢) 莩( 棚) 2 卸 我们选择局部平衡态分布函数矗o 任，舌) ，g = l ，2 ，3 ，4 ) 满足 椴z ，亡) = 互1 “( 茁，t ) 同时修正项函数最( z ，亡) = 1 ，2 ，3 ，4 ) 满足 则有 & ( z ，亡) = 士d l u ，i = 1 ，3 ， 士d 2 u ，t = 2 ，4 ， ( 1 2 7 ) ( 1 2 8 ) ( 1 2 9 ) ( 1 2 1 0 ) ( 1 2 1 2 ) ( e 门& ) = ( 如- 怕e 2 - d - e a - d 2 e 4 ) v u 划差+ 2 也是， ( 1 2 1 3 ) 莩( 印v 删。) = 三莩( 龟。v 胤= 三( 毫+ 毫) 把方程( 1 2 1 3 ) 和( 1 2 1 4 ) 代入方程( 1 2 1 0 ) 得 ( 1 2 1 4 ) 象+ 等( d - 差+ d 协o u 2 一丢( 丁一三) ( 蠡+ 毫) = 。 2 朋， 9 。 福建师范大学赖惠林硕士学位论文 为恢复方程( 1 2 1 ) ，只需令 学呐，学吨嘉一三) = 以 ( 1 2 加) 在计算过程中，我们令 z = a y = c g = h ， a t = t e 2 ( 1 2 。1 7 ) 则有 下= 三+ 万2 v a t ，d i e - - 譬，d 2 s - - 而b 2 一a t h ( 1 2 1 8 ) 下。互+ 1 r ，万百，万一 l l 乙埘， 其中b 1 ，b 2 和矿为已知参数。 1 3 数值实验 为验证上述模型的正确性和有效性，本节给出对流扩散方程数值模拟的两个例子。 对于边界处理，我们统一采用郭照立等人提出的非平衡态外推格式【7 6 】。除此之外，为测 定模型误差精度，我们定义总体相对误差( g r e ) 为 g r e = 掣黎若产 3 m 其中u ( x i ，亡) ，让( 戤，t ) 分别为数值解和解析解，在所有格点进行求和。 1 3 1 算例1 我们给出含有初始条件和边界条件的对流扩散方程 4 2 ，删 6 。爱+ 皤= ( 器+ 券) 肛 3 0 ) 2 ，6 蕊+ 6 2 瓦= i 丽+ 万= ( o ，州， ，t ) = 珏扛，3 0 ，t ) = u ( o ，y ，t ) = 牡( 3 0 ，y ，t ) = 0 ， ，o ) = 唧 刍扛+ ) 或n 筹咖嚣， 解析解为 u ( z m 亡) = e x p 一2 优+ 去p + 彰) ) s 泣筹s 组等， 其中 p = 孔孟) 2 + ( 抛 1 0 卜 o y 丝珧出 出 、-i_i_j-il_【 第一章二维对流扩散方程的格子b q k 模拟 模拟中采用如下一族参数：o = 0 1 ，= 0 i ，h = k = 0 0 0 0 4 ，p = 00 0 4 。我们分 别给出了当扛l o o o o 和t = 5 0 0 0 0 的三维可视化模拟结果和在y = 1 5 处的数值解与解 析解对比。见图i - 2 和1 - 3 。我们还给出了不同时刻的全局相对误差见表11 。从模拟 结果我们可以发现，经过长时间的演化，数值解与解析解依然非常吻合。 图1 - 3 ：t = 5 0 0 0 0 时模拟结果一左边为数值解三雏可视化图，右边为在y = 1 5 处数值解与解 析解的对比 脚t f o r t h w h e e n c 。t m = p 5 0 0 0 0 0 警n 怒墓搿d 蛆d v i s u a 蛳ln u i c m d e r s i c 删a ls 。i m n u l a t a i t o h e nr p e r o s u n l t l ；0 f y ：1 5 福建师范太学赖惠林硕士学位论文 表l1 ：算例1 不同时刻的总体相对误差 t h e t h e 出b mr d “i v ee r r o r f o r8 0 l u t i o r b o f e x 8 p b 、1 a td 证e r e n t t l m ， 13 2 算例2 藐们给出舍有初始条件的对流扩散方程州 伟8 ub 。o u + * 舅叫象+ n = 【0 ，吼 iu ( z ，”，o ) ：e x p 一! 二二! ! l ：；：i ! = 。! 壁) ， 解析解为 咖，归南m 卜虹生噪鹋掣 模拟中采用如下一族参数：a x = 0 0 1 ，t = 00 0 0 1 ，b l = b 2 = 0 8 ，p = 00 1 。我 们分别给出了t = 1 0 和t = 12 5 的三维可视化模拟结果和不同时刻最大剖面处数值解 与解析解之间的对比。从解析解可知：当t = 1 0 时，最大剖面为口= 1 3 见图1 - 4 ； 当t = l2 5 时，最大剖面为= 1 5 ，见图1 - 5 。我们还给出了不同时刻的总体相对误 差见表1 2 。可以发现数值解与解析解吻合较好。 h 1 o j o n 梯r b d n u m e 栅a d y 蚓m m p l3 帅10 图1 - 4 ：t = 1 0 时模拟结果：左边为数值解三维可视化图，右边为在最大剖面y = 1 3 处数值 解与解析解的对比 r 培h tf o r t w h 。h e 。n 。t p = 一1 o 。l e 。f f t 璺孟麓i 已盏p 型。：嚣璺i 舞一a t a t i o t n h e r e p s ，u 砌l t 。o f y ：13 一_：翻冒_ 第一章二维对流扩散方程的格子b g k 模拄 图1 - 5 ：t f f i l2 5 时模缸结果：左边为数值解三堆可视化图，右边为在最大剖面y = 15 处数 值解与解析解的对比 脚tf o rw t h e h e n 。t w = 1 。2 i 日5 ：。l o f e f t 怒兰搿狩出嚣忠絮i o t h e nr p e s 删u l t 。o fy ：1 5 表1 _ 2 ：算倒2 不同时刻的总体相对误差 t h e t h eg l o b a lr e l a t i v ee r r o r f o ra o l u t 4 0 n so f e x a m p l e2a td i f f e r e n t t h n e $ 1 4 本章小结 本章通过对标准b o l t z m a n n 方程添