（计算数学专业论文）流体流动的格子boltzmann方法.pdf

华中科技大学硕士学位论文 摘要 建立在微观模型上的格子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 l 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 l t z m a n n 方法求解，区域边界 的分布函数通过粗细网格转换与非平衡态外推获得。 对于实际流动模拟，本文模拟了三种情况的方柱绕流。在第一种情况中，方 柱位于流场中央，模拟了卡门涡街现象，给出了斯特鲁哈数随雷诺数变化曲线； 在第二种情况中，方柱位于流场壁面，分析了雷诺数对方柱后回流区的影响：在 第三种情况中，两方柱并列在流场中央，考察了方柱间距对流场的影响。 关键词：格子b o l t z m a n n 方法，隐格式，非均匀网格，嵌套边界，曲边界 伯肃流，方腔流，方柱绕流 华中科技大学硕士学位论文 a b s t r a c t t h e1 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 o s i m u l a t ef l u i df l o w si ar 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 l f l u i dd y n a m i c sa p p r o a c h e s ，t h em 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 ，a n de a s yt 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 u ne q u a t i o nf o rf l u i dd y n a m 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 es o m ep r o b l e m sn e e dt ob es o l v e d s h c ha sn l en o n u n i f o r ml a t t i c e a n dc a r v eb o u n d a r yc o n d i t i o n sa r ev e r yd i 箭c u l tf o rt h es t a n d a r dl a t t i c eb o l t z m a n n m e t h o d d e a l i n gw i t h t h ec o m p u t a t i o n a lm o d e l ，b o u n d a r yc o n d i t i o n sa n dn o n - u n i f o r ml a t t i c em e t h o d a r es t u d i e di n d e t a i l i nt h ep a p e r f o rc o m p u t a t i o n a lm o d e l ，a l li m p l i c i ts c h e m ei s p r e s e n t e d ap a r a m e t e r f a c t o ri si n t r o d u c e di n t ot h el a t t i c eb o l t z m a n n e q u a t i o na n dt h e p a r a m e t e r f a c t o ri s i n d e p e n d e n t w i t ht h er e l a x a t i o nt i m ef a c t o r f o r b o u n d a r y c o n d i t i o n s ，an e ws c h e m ei sd e v e l o p e di nw h i c ht h em a c r o s c o p i cq u a n t i t i e so nt h e c o m p u t a t i o n a lb o u n d a r yn o d e sc a n b eo b t a i n e df r o mt h em a c r o s c o p i cq u a n t i t i e so nt h e p h y s i c a lb o u n d a r yn o d e s f o rn o n u n i f o r ml a t t i c em o d e l ，an e s t e db o u n d a r yc o n d i t i o n s c h e m ei sp r o p o s e d i nt h es c h e m e t h ec o m p u t a t i o n a ld o m a l ni sd i v i d e di n t os o m e r e g u l a rs u b d o m a i n s w i t ht h en e s t e db o u n d a r i e s t h ed i s t r i b u f i o nf u n c t i o nc a nb e c o m p u t e d b y u n i f o r ml a t t i c eb o l t z m a n nm e t h o dw i t h i n s u b d o m a i n sa n d b y n o n e q u i l i b r i u me x t r a p o l a t i o ns c h e m e o nt h eb o u n d a r i e s t h r e ee a s e so ff l o w p a s ts q u a r ec y l i n d e r sa r ei n v e s t i g a t e di nt h ep a p e r f o rc a s e1 ， t h e s q u a r ec y l i n d e r i sl o c a t e da tt h ec e n t e ro ft h ec h a n n e l t h ek a r m a nv o r t e xi s s u c c e s s f u l l ys i m u l a t e da n dt h es t r o u h a ln u m b e r sa td i f f e r e n tr e y n o l d sn u m b e r sa r e g i v e n ；f o r c a s e2 ，t h es q u a r ec y l i n d e ri sl o c a t e da taw a l lo ff l o w f i e l d ，t h er e c i r c u l a t i o n l e n g t h sb e h i n dt h es q u a r ec y l i n d e ra td i f f e r e n tr e y n o l d sn u m b e r sa r ec o n s i d e r e d ；f o r c a s e3 ，t w os q u a r ec y l i n d e r sa r r a n g e ds i d eb ys i d ei nt h ec e n t e ro f t h ec h a n n e l ，t h ef l o w f e a t u r e sa td i f i e r e n ts p a c i n gr a t i o sa r es t u d i e d k e y w o r d s ：l a t t i c eb o l t z m a n nm e t h o d ；i m p l i c i ts c h e m e ；n o n u n i f o r mm e s h ； n e s t e db o u n d a r y ；c u r v eb o u n d a r y ；p o i s e u i l l ef l o w ；d r i v e nc a v i t y f l o w ；f l o wp a s tt h es q u a r ec y l i n d e r s i i 独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得 的研究成果。尽我所知，除文中已经标明引用的内容外，本论文不包含任何其他 个人或集体已经发表或撰写过的研究成果。对本文的研究做出贡献的个人和集 体，均已在文中以明确方式标明。本人完全意识到，本声明的法律结果由本人承 担。 学位论文作者签名：王f 袁銎 日期：2 口醉年3 月加同 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即：学校有 权保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和 借阅。本人授权华中科技大学可以将本学位论文的全部或部分内容编入有关数据 库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 保密口，在年解密后适用本授权书。 本论文属于不保密回。 ( 请在以上方框内打“、，”) 学位论文作者签名：王r 趟 同期洳斧岁月p f 1 指导教师签名：旋徉易 同期辨岁月7 0 日 华中科技大学硕士学位论文 1 1 计算流体力学的兴起 1 绪论 流体( 如空气，水等) 是自然界中最常见的物质，与人们的日常生活密切相 关。流体力学是研究流体运动规律的一门学科，其应用范围非常广泛。航空、航 海、宇航、机械、动力、水利、冶金、化工、石油、海洋、气象、天文、生物等 部门无不涉及流体力学问题。经过多年发展，流体力学已经取得了丰硕成果。十 七世纪末到十九世纪初，流体力学理论在实验的基础上得到极大的丰富和发展。 人们建立了流体流动所普遍遵循的n s 方程。由于n s 方程是非线性的，而 复杂流动中仅有极少数问题有解析解，在工程实际中遇到的大多数流体力学问题 求得精确解或解析解是十分困难的。 为了解决工程实际问题和深入研究流体动力学行为，人们利用相对运动原理， 建立了地面实验设备，如水洞、风洞、激波管、电弧加热器等，直接测量流动参 数，获取力、力矩与热力学方面的数据，从而出现了流体力学的一个分支：实验 流体动力学( 如实验空气动力学等) 。实验空气动力学不仅是获取航空航天飞行器 设计所需气动数据的重要手段，而且从实验研究得到的一些新的概念、新的规律， 促进了航空航天飞行器的发展。从流体动力学基本方程所推演出的相似定律、相 似准则是实验的理论基础，但几乎所有的实验设备都不能完全满足所有相似参数 的要求，另外实验研究费用离昂和场地安全等条件等亦对实验做了很大限制。 2 0 世纪5 0 年代以后，计算机发展非常迅速，差不多每两年计算机速度提高 一个数量级，同时计算机的价格和运行费用也在不断下降，使得数值求解流体动 力学成为可能。2 0 世纪6 0 年代起，许多流体力学专家、计算数学专家研究求解 流体力学基本方程数值方法，研制相应的计算机软件，在这种情况下，计算流体 力学继实验流体力学与理论流体力学之后成为解决流体问题的第三种方法，并且 渗透到现代科学的相关学科和工程之中。 华中科技大学硕士学位论文 计算流体力学的兴起促进了实验研究和理论分析方法的发展，为简化流动模 型的建立提供了更多的依据，使很多分析方法得到发展和完善，例如目前在飞机 3 2 , k 中广泛应用面元法。然而，更重要的是计算流体力学采用它独有的研究方法： 实际问题一数学建模一数值模拟一流体运动规律。这种方法的特点如下：( 1 ) 给 出流体运动区域内的离散解，而不是解析解，这区别一般的分析方法；( 2 ) 它的 发展与计算机技术的发展直接相关，这是因为可能模拟的流体运动的复杂程度、 解决问题的广度和所能模拟的物理尺度以及给出解的精度，都与计算机速度、内 存、图形处理直接相关；( 3 ) 若物理问题的数学提法( 如数学方程、边界条件等) 是正确的，则可以在较广泛的流动参数( 如马赫数、雷诺数、飞行高度、气体性 质、模型尺度等) 范围内研究流体力学问题，且能给出流场参数的定量结果，这 常常是风洞实验和理论分析难以做到的。计算流体的兴起促进了流体力学的发展， 它不但可用于研究已知的一些物理问题，而且可用于发现新的物理现象。例如 c a m p b e l l 和m u l l e r 等人在数值实验中，发现了亚声速斜坡绕流的分离现象，在 风洞实验中作了证实；又如k i m 和m o i n 在数值计算中发现了倒马蹄涡，后来被实 验所证实。 事实上，理论分析、实验研究和数值模拟是研究流体运动规律的三种基本方 法，它们的发展是相互促进的。在数值模拟中，实际问题中所求解的偏微分方程 组十分复杂，其数值解的现有数学理论尚不够充分，严格的稳定性分析、误差估 计和收敛性理论的发展还跟不上数值模拟的进展。所以在计算流体力学中，仍必 须依靠一些较简单的、线性化的、与原问题有密切关系的模型方程的严格数学分 析，然后再依靠实验研究结果，验证计算方法的可靠性，从而进一步改进计算方 法。 1 2 格子b o l t z m a n n 方法简介 模拟流体i g z j j 的数值方法有两种途径，即基于宏观连续模型的自顶向下的方 法和基于微观离散模型的自底而上的方法。传统计算流体力学的数值方法大多采 用第一种形式。这些方法以非线性微分方程为出发点，采用有限差分、有限体积、 华中科技大学硕士学位论文 有限元或有限谱等离散方法对微分方程进行离散，得到代数方程组或常微分方程 系统，然后再用标准的数值方法求解。在这类方法中，人们往往着重分析从连续 微分方程到离散代数方程的截断误差，而忽视了一些物理量的守恒性。对一些系 统而言，这种守恒性是十分重要的。另外，计算复杂性和数值稳定性是这类方法 的重要问题。 格子b o l t z m a n n ( l a t t i c e b o l t z m a r m ，l b ) 【l 2 j 方法是近十几年来发展起来的 一种新的数值就算方法。与以宏观连续方程的离散化为基础的传统数值方法不同， l b 方法的基础是微观模型与细观运动论方法。l b 方法的基本思想是构造简化运 动论模型，该模型能够反映微观的物理本质，并且宏观平均上能够满足所遵循的宏 观方程。l b 方法保留了分子运动学的许多优点，如物理条件清晰、边界条件易于 实现、具有完全并行性等。l b 方法与其它方法相比具有以下重要特点：第一，从 宏观方程出发的数值方法的对流项是非线性的，而l b 方法的时间演化方程分成 碰撞和流动过程，方法的对流算子在相空间是线性的，结合碰撞算子，通过多尺 度展开，可以恢复宏观方程的对流项；第二，在宏观n - - s 方程中，压力满足p o i s s o n 方程，数值计算时要特殊处理，而在l b 方法中，在近似不可压的条件下，可以 导出不可压n s 方程，流场压力通过状态方程来描述，并且计算格式是线性显 格式；第三，传统分子运动论中相空间是一个完全的函数空间，平均过程需要整 个空l 、自j 的信息，而l b 在相空间中使用一个最小限度的速度集合，运算由简单代 数运算组成，因而微观分布函数与宏观物理量的计算都非常简单；第四，l b 方 法的边界处理灵活简单，可以根据问题的物理意义和实际需要进行适当选择，计 算效率较高。 l b 方法源于格子气动机( l a t t i c eg a sa u t o m a t a ，l g a ) 。l g a 是一种利用离 散格子和离散时间的离散粒子运动论方法，是元胞自动机( c e l l u l a ra u t o m a t a ， c a ) 在流体力学中的应用。第一个l g a 模型是h p p 模型由h a r d y 、p a z z i s 与 p o m e a u l 3 1 于2 0 世纪7 0 年代提出。该模型采用二维正方形格子，每个格子点上有 粒子驻留，在同一个格点上的流体粒子按照一定的规则进行碰撞，并在一个时间 步内沿网格线运动到相邻的格点。在碰撞和流动过程中，h p p 模型遵循定的物 3 华中科技大学硕士学位论文 理规律，即保持质量和动量守恒，能够反应流体的一些基本特征，如输运性质等。 但是由于正方形格子缺乏足够的对称性，h p p 模型不能反应宏观方程所对应的非 线性和耗散效应。1 9 8 6 年，f r i s h 【4 等人在其关于二维流体动力学的格子气动机方 法研究中开创性的认识到：格子对称性在导出n - - s 方程很重要。他们第一次从 六边形的格子气动机得到j 下确的n s 方程。据此，他们提出了一个新的l g a 模 型，即f h p 模型。w 0 1 f r a m 5 1 和f r i s c h 6 1 等人工作为l g a 方法奠定了理论基础。 在l g a 中，可以定义一组描述粒子有无的布尔变量”，( x ，r ) ( 仁1 ，2 ，m ) 表示空间点x 处在时刻f 是否存在离散速度的粒子，m 表示每个格点处粒子的离 散速度方向的个数。l g a 的演化方程如下： 月。( x + c a t ，f + ，) = h ，( x ，f ) 4 - q ，( n ，( z ，f ) )( 1 1 ) 其中q 。( n ，( x ，f ) ) 称为碰撞算子，表示由于粒子之间的碰撞所引起的速度e ，的粒 子速度变化。在l g a 中，为提高存储效率，常常采用利用不相容原理，即在给 定时刻和节点上具有某个给定速度的粒子不多于1 ，这一原理导出了f e r m i d i r a c 局部平衡态分布。采用布尔型变量，l g a 模型的数值稳定性问题得到解决，但在 运行的过程会引入随机噪声，模拟结果往往含有统计噪声。为消除这种噪声， m c n a m a r a 和z a n e t t i 【7 j 在1 9 8 8 年首次提出在l g a 中使用布尔变量的月，的统计平 均量厂( 或称为单粒子分布函数) 进行演化，即l b 方程代替l g a 方程演化。在 l b 方法中，基本变量为平均粒子分布，他们是细观变量。l g a 模型的另一个缺 点是碰撞算子具有指数复杂性( 2 ”) ，在该模型中仍旧保留。h i g u e r a 和j i m e n e z 8 l 对l b 方法作了一个重要的简化，他们引入平衡态分布函数t ，：”，对碰撞算子作线 性化处理，即q ，( 厂) = 世。( ；一；”) ，其中k 。= 船，阢为碰撞矩阵。该模型使碰 撞算子的复杂性由2 ”降为m 2 。不久，h i g u e r a 9 等人提出了一种线性稳定的强化 碰撞算子方法。此后几个不同的研究小组f 1 0 1 川分别提出了一种更简单的模型，即 单松弛模型( l b g k 模型) 。在这种模型中，碰撞过程用趋于某一平衡态的松弛过 4 华中科技大学硕士学位论文 程代替，矩阵由松弛时间确定，即k 。= 一( 1 r ) j ，。该模型大大简化了计算r 并且 在一定条件下可以恢复n s 方程。该模型的提出使格子方法的研究达到一个新水 平，并且使l b 方法得到广泛应用。 1 3 格子b o l t z m a n n 方法的国内外研究状况 l b 方法正处于不断发展阶段，近年来在计算模型、边界条件、多重网格与区 域分裂、l b 方法的应用等问题上都取得了较大进展。 ( 1 ) 计算模型 1 9 8 6 年，f r i s c ht 4 ，w a l f r a m 【5 1 等分别提出用格子气动机与元胞自动机解决流 体流动问题，格子气动机在流体中的应用研究逐渐引起学者的注意。对b o l t z m a r m 方程线性化，并结合b g k 1 2 1 碰撞算子，q i a n 【l o ，c h e n 【1 1 】分别选择平衡态函数恢 复n - - s 方程，l b 方法的最初实用计算模型建立。z o u 1 3 1 等通过对密度p 与速度 “的选择，得到更实用的两种计算模型。不久z o u 1 4 1 与h e 1 5 】等又提出了l b 分 析解决办法。与习惯的密度函数相比较，l u o 【”1 等采用压强分布函数建立了新的 演化模型。g u o 17 】等通过引入一种新的分布函数，提出了一种更为理想的不可压 计算模型，计算效果有很大的提高。最近，f a n g 1 8 l 等提出松弛密度模型。 ( 2 ) 边界条件 l b 方法的边界处理方法来源于格子气动机方法，伴随着l b 方法的发展而不 断改进。w a l f r a m 5 1 提出反弹规则，即：当一个粒子分布流向壁面节点时，此粒子 分布反射回原节点。反弹规则在l b 方法的数值模拟中得到很好的应用。 c o r n u b e r t 1 9 1 在此基础上提出修f 反弹格式，计算效果在一些流体流动模拟中得 到改进。h e 15 】等提出了一种半途反弹格式，此外反弹规则的另外种改进是镜面 反弹方式。为提高数值精度与稳定性，许多学者提出了多种边界处理方法。n o b l e 2 0 1 根据压力约束，提出无滑移边界动力学格式：c h e n 2 1 】根据插值方法了二阶外推模 式；i n a m u r o 【2 2 j 利用反方向滑移速度，提出了壁面反方向滑移的边界处理方式； d u p i u s 2 3 】根据质量守恒条件，提出质量守恒的边界处理方法。最近，g u o l 2 4 】等提 华中科技大学硕士学位论文 出的非平衡态外推方法在复杂边界与高雷诺数的流体模拟方面有很好的效果。 ( 3 ) 多重网格与区域分裂 在流体模拟中，我们的目标是精确的模拟连续介质，从而要求网格尺度的极 限为零，但直接加细网格，会使计算量剧增，在此意义上讲，计算机的发展远远 跟不上实际问题的需要，探索高效、快速算法成为当前的数值模拟的热点。l b 方 法的一个限制是其离散格式使用的是一一类特殊的均匀和规则格子。为提高计算精 度和效率，解偏微分方程的多重网格与区域分裂的思想开始引入到l b 方法中来。 f i l i p p o v a 2 5 0 6 1 ，y u 2 卅，d u p i u s 2 3 1 等国外学者在这方面取得了较大的进展，但对 实际计算的许多问题有待解决。关于多重网格与区域分裂是l b 方法的充满活力的 研究方向之一。 ( 4 ) l b 方法的应用 l b 方法从上世纪8 0 年代中期提出以来发展迅速，渗透到流体力学的多个领 域。对于许多物理问题建立了相应的格子气模型，其中包括多孔介质、多相流、 湍流、反应扩散，等等。b e n z i 2 9 l 、q i a n 2 9 、c h e n 1 1 等在l b 方法的发展上写了三 篇综述性文章，对l b 方法的应用作了介绍。应当指出的是，l b 方法是一种崭新 的计算方法，许多问题还没有解决，特别是对于一些挑战性领域，如多相流的各 相闻的分离处理、完善的湍流场模拟等方面都在探讨阶段。 1 4 本文研究工作介绍 本文对l b 方法的计算模型、边界处理与非均匀网格方法进行研究。在第二 章中，介绍了l b 方法的基本原理，并对当前流行的计算模型和边界条件处理方 法数值稳定性作了实验对比。在第三章中，提出了一种隐格式方法，该方法中引 入的参数因子不受松弛时间限制。第四章中，设计了一种嵌套边界的非均匀l b 方法，该方法能够灵活处理复杂流场划分问题。第五章中，提出了一种曲边界处 理方法，并对圆柱绕流作了实验对比。第六章中，对三种情况的方柱绕流作了数 值研究，对于一些物理参数，给出了定量数据。 华中科技大学硕士学位论文 2 格子b o l t z m a n n 方法基本原理 近几年来，l b 方法作为一种模拟流体流动的数值方法已经引起人们的普遍 重视。l b 方法是从微观的粒子尺度出发建立离散的速度模型，在满足质量、动 量和能量守恒的条件下，得出粒子的分布函数，然后对分布函数进行计算得到压 力、速度等宏观量。l b 方法的动力学演化方程是： ，( x + c e ，a t ，r + f ) = ：( z ，r ) + q ，( ( z ，f ) ) ( i = o ，1 ，b ) ( 2 1 ) 其中，表示沿速度方向i 的粒子分布函数；q ，( 厂( x ，f ) ) 是碰撞因子，表示 由于碰撞过程引起，( x ，r ) 的改变率，仅依赖于局部分布函数；c = a x a t ，x 和 t 分别为空间和时间步长；e ，表示粒子的运动方向，我们考虑九速二维平面格子 ( 参考图2 1 ) ： e 。= ( o ，o ) ， q = ( c o s ( 一1 ) n 2 】，s i n ( i 一1 ) 7 r 1 2 ) ，i = 1 , 2 ，3 ，4 ， ( 2 2 ) p ，= 压( c o s ( 一5 k 2 + 万4 ls i n ( f 一5 ) l r 2 + r c 4 d ，i = 5 , 6 ，7 ，8 ， 格点处的局部宏观密度、速度由以下各式确定： p = ：，p u = c q ：， ( 2 3 ) f ， 动力学演化方程2 1 把流体运动分成两部分：流动过程与碰撞过程。在流动 过程中， ：( x ，f ) 沿e 方向运动，并且与之邻点交换数值；在碰撞过程中，在个 节点上从相邻节点运动来的粒子发生碰撞，根据质量、动量和能量守恒规则改变 粒子分布函数，具体的讲，采用b g k ”l 松弛模型： q ，= 一( ：一；”)( 2 4 ) 其中，r 是松弛时间，w 局部平衡函数。碰撞函数n ，满足局部质量、动量 守恒： 7 q ，：o ， q 。p ，= o ， 2 - 5 2 1 数学基础 图2 1 d 2 q 9 模型 ( 1 ) 平衡态分布函数 采用b g k 松弛模型( 方程2 4 ) ，l b 演化方程： ，o + c 口，r ，t + a t ) = ，( x ，r ) 一圭( e r ) 一z ”( x ，f ) ) ( i = o ，1 ，b ) ( 2 - 6 ) 根据上式，我们知道确定松弛时间f 条件下，只要确定恰当的平衡态函数；“， l b 方程即可正常演化。平衡态函数z w 有赖于局部宏观变量( p 和u ) ， 现采用 q i a n 1 0 1 等d 2 q 9 平衡态分布： 一1 + 3 半“s 学乩s 料 眨， w 为权系数其数值分别为：w 0 = 4 9 ，w ，= 1 9 “= 1 ：4 ) 和w ，= 1 3 6 ( f 。5 ：8 ) 。 对于离散速度e 满足以下等式：i 亡。一n ( 2 8 a ) e ，。= p ，。= o 2 = ll = s 圭：2 2 8 6 f = 1 一一 8 华中科技大学硕士学位论文 e t 。e | b = 4 6 筇 p ，。e ，口g ，= p ，。p ，日= 0 e j j p e ? i e 扩2 6 铆 8 e j ，e j 乒j ；e j b = 4 。p 驴一栅e f = 5 其却，a ；8 = 6 婶6 + 6 。+ 6 。8 6 畦 通过方程2 7 与2 8 ，我们可以得到下面各式 矿= p c 吃铲= 彤 8 y i = c2 叩，铲= p u u + p i i = o ( 2 1 0 a ) ( 2 i o b ) ( 2 1 0 c ) 嘤= c 3 e i a e i f l p 。f ：”= p ( “f + “+ ) ( 2 i o d ) f _ 0 j 其中p = ；c 2 p ( 关于2 1 0 的详细的推导过程参考附录1 ) 。 ( 2 ) 宏观方程 为了从l b 演化方程导出宏观n s 方程，对l b 演化方程中分布函数，采用 多尺度技术，展开成局部平衡分布函数，”： j = 七? ( 2 1 1 ) 其中f 为k u n d s e n 数，满足占 1 。在计算网格充分细的条件下，可认为占与出同 阶的，即o ( s ) = o ( a t l 。 记，：。= ，“，；”= ：”，则： = 0 ，e ，= 0 ( 2 1 2 ) 令昙+ 钾，v = d j ，对l b 演化方程( 方程2 6 ) 进彳t t a y l 。r 展开： 9 d d 酣 d & 墙 舟 碍 z 2 汜 0 华中科技大学硕士学位论文 蹦扣a 7 t d 2 ，十p 2 ) 一去( ；叫”) 2 1 3 由方程2 1 1 ，2 1 2 ，并且由于。( ) = o ( a t ) ，上式可以转化为： d 。，。) + ：i ) + 等d ? ：0 1 = 一 ，+ o ( ( r ) 2 ) ( 2 1 4 ) 对方程2 1 4 进行0 阶近似，可以得到： d ， 【0 ) = 一圭+ o ( f ) ( 2 1 5 ) 方程2 1 5 可以变形为： ，”= 一，f 0 1 十0 ( f ) ( 2 1 6 ) 从而方程2 1 4 可以转化为： d 。，。) + ( ：1 一f ) 圮f ，。) ：一! ：1 ) + 。( ( f ) 2 ) ( 2 1 7 ) 由于d ，= 昙+ c e ，v ，根据方向i ，对方程2 1 5 两边求和： 壹笔+ 杰崛wc ：一杰+ o ( a r ) 善_ + 丢w k 一否。 ) ( 2 1 8 ) 利用方程2 1 0 a ，2 1 0 b ，2 1 2 ，可得： 孚+ v ( ) = 0 ( o ( a t ) ) ( 2 1 9 ) 对方程2 1 5 两边同时乘以c e ，并根据i 方向求和，根据2 1 0 ，其中方程左端可以 变为： 88 月，【u ) e = ( 一q - c 2 d l e j v f , 姆) j 1 0i = o “ a ( c e i f ，1 ) a ：塑+ v - ( o ) 钟 根据2 1 2 ，方程右端： + 甲( c2 q g ，) ( 2 2 0 ) + d ( z x t ) j ：o ( z x t )( 2 埘) 从而可以得到下面方程： 挈+vn【。1：o(d(f)ol 、， 方程2 1 7 中根据i 方向求和，利用前面相似步骤可以得到 詈w 伽)+ d ( ( 矿) 其中萎珥，。根据方程2 1 9 ，2 2 2 可变为 西a ( o 耐_ + v ( 删) ) + v 3 西a ( o 优冉( 删冉( 詈m n ( 0 】) = o ( a t ) 从而可以二阶近似的连续性方程： 鲁冉( ) = o 方程2 1 7 中两边同时乘以卯，并根据i 方向求和，可以得到 由于 其中 罢再肌拉玉卿o ) + 。) ( 2 2 2 ) ( 22 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) 知删吲静百o f j o ) + 驴8 ，。聊) = # 。( a 。p u + v y i ( o ) ) 、西o 委8 也，钟，+ 喜如局。刖) ( 2 刀， 2 嘉( 警m n m ( 警+ v ( v 一) 1 1 ， 一f 一 ，rl 曾c 。h 0 ， d 。 r ) 一2 一 ( ) o d 。 ( d 1 1 0 z d 。 ) o zq 。 ， 华中科技大学硕士学位论文 v ( vr ( 0 ) ) _ v 小，；( i c 2 p 5 ，+ “，+ ) ) 2 去c2 v 口( v f ( p u ) + v 。( p ) + v 口( 。) ) = ；c2 v p ( v ，( 删。) + v 。( 用劫+ ；c 2 v f ( v f ( 功 ( 2 2 8 ) ：；c ：( v ：( p u ) + v ( v ( ) ) r ) + ；c2 v ( v ( 倒) ，) ：；c 2 ( v2 ( ) + v ( v ( 彤) ) ) + ；e ：审( v ( ) d ( 推导中用到v ( v ( 脚) ) 。= v ( v ( 肛) ) ) 利用2 1 0 c ，2 2 2 ，2 2 8 ，方程2 2 7 可化为： 扣一v 掣+ 掣m c 训， ：， + 土3 。2 ( v 2 ( 膨) + v ( v ( 彤) ) ) 由于，在低马赫数下有下列近似： 0 ( 印) = 0 ( 印) = o ( m2 ) ，0 ( “) = 0 ( 肘) ( 2 3 0 ) n , n2 2 5 ，2 2 9 ，2 3 0 方程2 2 7 可以转化为： c q 研，。= d ( f ) + o ( m 2 ) + 妻c 2 ( v 2 ( ) + v ( v ( 彤) ) ) ( 2 3 1 ) 从而在低马赫数条件下，利用方程2 1 0 c ，2 3 1 ，可以从方程2 2 6 得到二阶近似的动 量方程： ，警+ 口( “) = 一即+ u ( v 2 ( 删) + v ( v ( 倒) ”( 2 3 2 ) 其中，0 为粘性系数，可以表示为： ( 2 r 一1 ) ( 缸) 2 d = = = = o 6f ( 2 3 3 ) 方程2 2 5 ，2 3 2 联立即为n s 方程组，在不可压条件下可咀简化为： v “= 0 ( 2 3 4 ) 宴“( 蚴：一半+ 。 (2-35)ot 口 。7 华中科技大学硕士学位论文 2 2 计算模型 关于计算模型近年来的发展动态，在绪论l3 中已经做了初步介绍，现在选 取当前国际比较流行的三中模型，即q i a n - c h e n 1 0 , 1 1 模型、z o u h o u 1 3 1 模型、不可 压模型。对于离散速度采用二维九速模型( 方程2 2 ) 。关于o i a n c h e n 模型， 即在平衡态分布函数中使用方程2 7 ，下面主要介绍z o u h o u 模型与不可压模型。 2 2 1 模型介绍 ( 1 ) z o u h o u 模型 在流体模拟中，l b 方法所反应的是可压缩方程，在处理不可压问题时，不 可压效应往往对计算结果带来很大影响。事实上从方程2 2 4 到2 3 3 为理想情况， 著考虑到密度p 的空间的导数项，可以得到： v “：一土“v 口 口 ( 2 - 3 6 ) 该方程的右端代表了不可压缩效应引起的误差，这个误差是来源于l b 方法本身 与计算网格或时间离散无关。为此z o u 等人提出了，使用一种新的速度，即令 v = p u( 2 3 7 ) 从向计算模型的平衡态分布函数为： 肛p 半“s 掣c “s 簪l 汜，s ， i c 。 c “i 并定义流动的宏观密度与宏观速度为： p = z ，矿= c e ， ( 2 3 9 ) 通过多尺度展开方法可以得到下面的宏观方程： 窭+ v “= 0( 2 4 0 ) 詈+ v 巾“) 一印州v 2 ”+ v ( v )( 24 1 ) 此模型可模拟定常的不可压模型。此后c h e n 3 0 1 等人推广到非定常流动计算，但计 华中科技大学硕士学位论文 算过程增加了p o i s s o n 方程求解。h e l u o 1 6 1 模型是使用人工压缩办法，即采用压 强分布函数进行计算，模型中附加条件限制了方法的运用。 ( 2 ) 不可压模型 为消除人工压缩效应的影响，郭照立等提出使用一种新的平衡态分布函数进 行计算。在其模型中，平衡态分布函数定义为： g ? = 一4 a 善+ s o ( “) c 一 五尝+ e ( “) ， c 。 y 善+ s 肋) ， c 。 i = 0 i = 5 ，6 ，7 ，8 ， ( 2 4 2 ) 这里叭旷w ，l 。半“s 字“s 粤l ，为权系数其数值捌为：l c c c l 。 w o2 4 9 ，w ，= 1 9 ( f = 1 ：4 ) 和w ，= 1 3 6 ( f _ 5 ：8 ) 盯，z ，y 为模型参数，满足： 五+ y = o - ，五+ 2 f = l 2 该模型的演化方程为： g ，0 + c q r ，t + a t ) 一毋g ，f ) = 一三g ，g ，f ) 一g o ( z ，f ) ) ( 2 4 3 ) 流体的宏观速度和压力可由分布函数计算得到( 对于宏观量的计算推导参考附录 ”= i = 1 朋= 丢瞽蝇 。， 通过多尺度展开技术，可以从( 1 ) 得到不可压条件下的n a v i e r s t o k e s 方程 粘性系数u 由方程2 3 2 确定。 ( 2 4 5 a ) ( 2 4 5 b ) 1 4 “vu+ 即 一 j j 甜 o 呷 i i 甜 “ 十 口抛一西 华中科技大学硕士学位论文 = = = = = = = = = = = = = = = = = = = = = = = = = = = = ；= = = = = = = = = = 一 2 2 2 实验对比 兰三! ：三! 图2 2 方腔流的计算区域与边晃条件 对于l b 模型算法的研究，前人取得了一些重要成果，对于这然而模型的稳 定性研究，当前刊发文章不多。本节方腔流为例，从两个角度以上三种模型进行 数值实验，即：( 1 ) 粒子速度定数条件下，雷诺数的计算对比；( 2 ) 雷诺数定数 条件下，粒子速度计算对比。 对于方腔流的计算区域与边界条件( 如图2 2 ) 。对于边界条件处理方法采用 非平衡态外推方法。取边长为1 的方腔，对于边界条件在方腔项部“= 1 ，v = 0 ， 在其余的壁面上“= v = 0 。方腔内部的初始条件“= v = 0 ，p = 0 1 。对于雷诺数， 定义为r e = u o l u ，其中，“。与分别为方脸顶部速度与长度。对于第一种情况 取粒子速度c = 1 0 ，表2 1 给出了三种计算模型在不同雷诺数下的计算稳定性对 比。对于第二中情况取雷诺数r e = 1 0 0 0 0 ，表2 2 给出了三种计算模型在不同粒 子速度下的计算稳定性对比。通过实验我们可以看出，不可压模型在复杂流体计 算方面有较好的稳定性。 华中科技大学硕士学位论文 表2 1三种计算模型在不同雷诺数下的计算稳定性对比 雷诺数r e 1 0 0 0 01 3 0 0 01 4 0 0 03 0 0 0 0 q i a n 1 。1 稳定稳定不稳定不稳定 z o u j 3 稳定稳定不稳定不稳定 g u o 1 7 】稳定稳定稳定稳定 表2 2 三种计算模型在不同粒子速度下的计算稳定性对比 速度c 1 0764 q i a n e l o 】稳定稳定不稳定不稳定 z o u 【1 3 稳定稳定不稳定不稳定 g u o 【1 7 】 稳定稳定稳定稳定 2 3 边界条件 2 3 1 边界条件处理方法 关于边界条件处理方法，在绪论1 | 3 中已经做了初步介绍，现在对应用比较 广泛的三种边界条件处理方法做介绍，即：反弹格式5 ，伸 、插值外推格式2 1 】和非 平衡态外推方法【2 4 ) 。 ( 1 ) 反弹格式 反弹格式主要用于无滑移壁面。所谓反弹格式，是指如果一个流体节点上的 粒子沿粒子运动方向流动一步达到边界节点，则该粒子沿原方向反弹回原流体节 点。为叙述方便，将l b 演化方程( 方程2 6 ) 分成碰撞和流动两个过程，即： t 碰撞：，+ ( z ，) = ；( x ，f ) 一( 工，f ) 一，w ( x ，r ) ) ( 2 4 6 a ) t 流动：； + c e ，a t ，h a t ) = ，+ ( z ，f )( 2 4 6 b ) 6 华中科技大学硕士学位论文 反弹格式可以表示为 厂_ ，( ，r ) = 厂( x l ，f ) ( 2 4 7 ) 其中为壁面格点，z ，表示关于方向i 离最近的流体格点。在反弹格式中，在 壁面上不发生碰撞，以下壁面为例( 参考图2 3 ) ，关于0 点的分布函数 、兀、 工分别有a 、b 、c 相对应时刻的分布函数确定： ( 0