 |
|
|
| 本站文章非会员原创,如果侵犯了你的权利请联系13978471165 |
| 作者：庆熙 毕业来源：仙剑奇侠博客 点击数： 更新时间：2010-7-7 |
|
摘 要 在上世纪八十年代后期提出的格子Boltzamnn方法克服了格子气方法的缺点，其本身也在不断的发展之中.格子Boltzamnn方法在流体运动计算方面展现了非凡的风采，成功地模拟了包括均相不可压缩湍流和多孔介质中的多相流动在内的流体动力学问题.但和成熟的流体动力学计算方法相比，特别在工程实际应用上，该方法还有许多值得研究的地方. 本文主要介绍工程实际应用时，具体模型的选择问题.首先从理论上对应用最为广泛的几种基本模型进行了详尽的分析和比较.选择了Poiseuille流动，然后从计算精度、数值稳定性和收敛速度这几个方面进行了细致的比较.从理论和实验两个角度验证了D2G9模型的优越性，为工程实际应用上模型的具体选择提供了一定的参考依据.通过研究二阶精确的格子Boltzamnn模型，提出了非牛顿流体.非牛顿流动性是使用幂法则模型实现的.它可以估算出模型的精确程度，同时不会限制这个模型.二阶精度由剪切变稀和剪切增稠液体的幂法则模型参数范围给出.这些结果与Gabbanelli等人的结果相比，精确度更高，并且得到了更快的计算效率.结果表明了格子Boltzamnn方法适用于非牛顿流体模拟. 对于实际流动模拟，本文应用二维9速度模型模拟了四种情况的方柱绕流问题.在第一种情况中，单个方柱位于流场中央，给出了流线图，等涡线图，模拟了卡门涡街现象，并计算了升、阻力系数，Strouhal数等参数；在第二种情况中，计算细长矩板截面柱绕流问题，得到了Strouhal数随着矩形长宽不同的比值下的变化情况；在第三种情况中，两个方柱并列位于流场中央，考察了方柱间距对于流场的影响；在第四种情况中，计算了水平来流为剪切流的方柱绕流问题，比较了速度梯度取不同值下流场的变化情况.所有有关力的求解均采用动量转换法.所得结果，包括流线、等涡线、升/阻力系数曲线等均与已有文献的实验或数值结果基本一致，显示LBM方法及其力的求解方法——动量转换法是有效的，能够精确的模拟各流场.其次，我们还引入一种两相耦合机制对D2G9模型进行了修正，从而使之可以正确处理气固两相流中输运相和颗粒相之间的相互作用.随后，我们模拟了后台阶流动，并和传统CFD方法的模拟结果以及修正其他模型的模拟结果进行了验证，得到了令人满意的结论.从一定程度上验证了两相耦合机制的可行性.通过软件模拟获得了水包油、过渡流型和油包水三种流型的典型模拟图.经分析发现：由软件模拟的流型特点和由探针获得的流型特点具有较好的一致性. 在本文最后，我们介绍了以经典算例一方腔流为例，对格子Boltzamnn方法的核心代码进行了优化的方法，主要讲述对时间和空间上的优化，优化的程序使计算效率提高数倍.在并行的框架下，核心演化的代码换为优化后的程序，计算效率有大幅度的提高.
Abstract In the latter of 80’s，the Lattice Boltzamnn Method（LBM）was introduced mainlyto cope with major drawbacks of its ancestor，the Lattice Gas Automata（LGA）．Eversince，it has undergone a number of refinements and extensions which have taken it tothe point where it can successfully compute a number of non trivial flows，raging fromhomogeneous incompressible turbulence to multiphase flows in porous geometries．Yet，when compared with conventional computational fluids dynamics methods，such as finiteelement，finite difference，it is apparent that there is still a way to go before LBM canachieve full engineering status． In this paper，we mostly focus on the choice of the basic LB models in theengineering application fields．Firstly，we expatiate the basic LB models in theory．Then，we simulate the Poiseuille flow with those basic LB models．And wecompare the simulation results from the computation precision、the numerical stabilityand the convergence rate．Finally，we draw a conclusion that the D2G9 model is the bestchoice in the engineering application fields． Simulation of Flow past square cylinder with LB Method.For the simulation of actual flow，we use D2Q9 investigate fourcases of flow past square cylinders in this paper.For case 1，one singlesquare cylinder is located at the center of the channel，we describe thestreamline contour，vortices contours，simulate the Karman vortex，then compute the lift coefficient，drag coefficient，Strouhal numbersetc.For the case 2，simulate the flow past a cylinder of rectangularcross-section;compute the change of Strouhal numbers varying withthe side ratio.For case 3:two square cylinders arranged side by side inthe center of the channel，the flow features at different spacing ratiosare studied.For case 4:we compute the linear shear flow over a squarecylinder，compare the evolution of flow with different velocitygradient.The results of thesimulation including the streamlines，vorticity contours，lift and drag coefficients etc.are agreed with thoseof available literatures，and show that LB method and itsmomentum-exchange method can achieve accurate results and obtainthe reasonable flow in detail. we employ a two-way coupling mechanisms to modify theD2G9 model．With the modified D2G9 model，we can handle with the interactionsbetween carrier phase and dispersed phase in the model．Then，we simulate abackward-facing step model，and the results are compared qualitatively with the result ofthe traditional CFD method and the other modified LB models．Though the comparison，we can see that the two-way coupling mechanisms can handle with the gas-solid twophases flows successfully． Three kinds of flow pattern，which are oil-in-water flow，transitional flow andwater-in-oil flow，have been got by simulation.According to the result of simulation，theoil-water two-phase flow pattern transition boundary model has been got by.By the analysisof simulation，the characteristic of three kinds of flow pattern of vertical oil which has beengot by analysis of the signals is consistent with results by simulation. We take the classical problem-cavity flow as an example and optimize the kerne codes of the LBM. The optimization include two aspects :time and space .The efficiency of the optimized code increased much more .In the parallel frame，the efficiency also increased if the kernel code is taken the optimized code.
LBM 目 录 第1章 概 述 1 1.1研究格子 Boltzamnn方法的意义 1 1.2 格子 Boltzamnn方法的发展历程 3 1.2.1孕育阶段 3 1.2.2 萌芽到成长阶段 3 1.3 格子 Boltzamnn方法应用概况及优缺点 5 1.3.1格子Boltzamnn方法应用概况 5 1.3.2格子Boltzamnn的优缺点 6 1.4本论文的研究目的 8 1.5 相关研究的综述与专注情况 8 第2章 格子Boltzamnn方法介绍 10 2.1 Boltzamnn方程的产生 10 2.2细胞自动机（CA） 11 2.3格子气自动机（LGA） 12 2.4格子Boltzamnn方法（LBM） 13 2.5 格子Boltzamnn的基本结构 16 2.6本章小结 17 第3章 格子Boltzamnn方法的基本模型比较 18 3.1 格子 Boltzamnn 方法基本模型概述 18 3.2 进行常压力梯度驱动的Poiseuille流动模拟比较几种基本模型 23 3.3本章小结 27 第4章 格子Boltzamnn方法的算法设计 28 4.1格子Boltzamnn方法的算法实现 28 4.2格子Boltzamnn方法的高效算法设计 30 4.2.1优化算法 30 4.2.2优化实验 32 4.3 本章小结 34 第5章 格子Boltzamnn方法的实际应用 35 5.1二阶精确格子Boltzamnn非牛顿流体的流动模拟 35 5.1.1理论背景 35 5.1.2方法和计算结果分析 38 5.1.3 本节小结 40 5.2 格子Boltzamnn方法的方柱绕流模拟 40 5.2.1 单个方柱位于流场中央的绕流问题 40 5.2.2 细长矩形截面住绕流问题 42 5.2.3 两个并列方柱的绕流问题 44 5.2.4来流为剪切流的绕流问题 49 5.3格子Boltzamnn方法模拟气固两相流 51 5.3.1对气固两相流的模拟模拟对象简介 51 5.3.2 计算结果分析 54 5.3.3本节小结 56 5.4 格子Boltzamnn方法模拟油水两相流软件设计 56 5.4.1 LBM油水两相流的关键因素选取 57 5.4.2 软件的设计 60 5.4.3 本节小结 63 5.5 简述格子Boltzamnn方法在其他领域中的应用 64 5.5.1 颗粒悬浮问题的模拟 64 5.5.2 热导和对流—扩散问题的模拟 64 5.5.3 偏微分方程的模拟 65 5.5.4 多相流和多元流的模拟 65 结论及展望 67 参考文献 68 致 谢 69 第1章 概 述
伴随着电子计算机的飞速发展以及各种新颖算法的不断出现，CFD已经形成了一门独立的学科，并且在航空航天、船舶、大型能源装置（如核电站）、新型交通工具、海洋工程、环境保护等众多工程技术部门和领域都得到了广泛的应用.随着计算技术的发展、巨型计算机的出现、计算方法的不断改进，计算流体力学在解决流动的理论和工程实际问题中愈加显示出它的巨大作用.目前，计算流体力学已经成为现代计算科学的最有力的推动力之一. 在计算流体力学中，传统的数值模拟方法可以分为两大类： （1）从宏观角度出发，基于连续介质假设，采用数值计算方法，求解全位势方程或Euler方程或N-S方程； （2）从微观角度出发，采用分子动力学的方法，对流动进行数值模拟.其中，格子Boltzamnn方法就是典型的一种.格子Boltzamnn方法（Lattice Boltzamnn Method，LBM）1.1.2格子Boltzamnn法（lattice Boltzamnn method）起源于格子气自动机（Lattice Gas Automata，LGA）.LGA方法是元胞自动机（Cellular Automata，CA）在流体力学中的具体应用，是空间、时间和速度空间都离散的一个虚拟微观模型，与以连续微分方程为基础的宏观计算流体力学方法有着本质的不同.LGA的微观特性使得它的边界条件非常容易实现，并且计算也很简单.因此，LGA方法非常适于处理边界复杂的问题.更为重要的是，LGA的计算具有局部性和并行性，非常容易在并行机上实现.LGA的出现不但为并行计算提供了许多新思想，而且对并行计算机制造技术产生了重要的影响. 但是，LGA方法也有许多不足之处.例如，由于含有随机因素，LGA的计算结果往往包含很大的统计噪声，LGA的宏观方程也不是标准的流体运动宏观方程.格子Boltzamnn方法是为克服LGA方法的一些内在不足而发展起来的一种新方法.LBM不但克服了LGA的缺点，继承了LGA的主要优点，而且还有许多新的优点，如计算量小、计算效率高、编程简单等.LBM的产生与发展，不仅在计算流体力学领域中产生了深远的影响，它所使用的处理方法和观点对其他许多学科也是富有启发性的. 格子Boltzamnn法是一种应用非连续介质思想研究宏观物理现象，并可平行运行，求解流体力学问题的新方法.它是由格子气自动机（lattice gas automata，简称LGA）方法发展而来的.该法把流体及其存在的时间、空间完全离散，把流体看成由许多只有质量没有体积的微小粒子组成，所有这些粒子同步地随着离散的时间步长，根据给定碰撞规则在网格点上相互碰撞，并沿网格线在节点之间运动.碰撞规则遵循质量、动量和能量守恒定律. 流体运动的宏观特征是由微观流体格子相互碰撞并在整体上表现出来的统计规律.该法是直接从微观模型出发，经过Boole化处理后进行计算，可认为是N-S差分法逼近的一种无限稳定的格式.被广泛应用于复杂几何边界流体流动、多孔介质流、多相流及反应流等. 格子气自动机的基本思想是，把计算区域分成许多均匀的正三角形（或正方形）的网格，而那些只有质量无体积的粒子只能在网格点上存在，并沿着网格线在网格间运动.当某一个粒子从某一网格点到邻近的网格点时，有可能和从其他网格点到达该点的粒子相碰撞.根据Pauli不相容原理，在同一时刻同一点上，沿着每一网格线运动方向最多只有一个粒子，流场中的粒子速度不是0（静止）就是1（设格子边长及时间间隔都为1）.以三角形网格为例，每一个网格上在某一时刻，其周围的6个网格上粒子沿着网格线聚集到该点，加上该点可能还有一个静止粒子，这样，可能有7个粒子在该点发生碰撞，然后根据碰撞规则再散射出去，演化为新的运动粒子流向各节点的邻居，形成格子气自动机. 1986年MeNamaxa和Zaneltti，提出把格子气自动机中的整数运算变成实数运算，建立了格子Boltzamnn模型，克服了格子气自动机的数值噪声的缺点.后来陈十一和钱跃宏采用了单一时间松弛方法，满足了各项同性，GalIean不变性，并得到了独立于速度的压力项．使格子Boltzamnn模型保留了格子气自动机的优点，克服了其不足，并在理论分析和数值模拟方面都具有很大灵活性，而且程序编制简单，计算效率较高. 从格子Boltzamnn方法诞生至今天已有20年，20年间，其在理论和应用研究等方面都取得了迅速发展，并逐渐成为在相关领域研究的国际热点之一，受到国内外众多学者关注.与之传统模拟方法不同，格子Boltzamnn方法基于分子动理论，具有清晰的物理背景.该方法在宏观上是离散方法，微观上是连续方法，因而被称为介观模拟方法.在许多传统模拟方法难以胜任的领域，入微尺度流动与换热、多孔介质、生物流动、磁流体、晶体生长等，格子Boltzamnn方法都可以进行有效的模拟，因此它被用于多种复杂现象的机理研究，推动了相关学科的发展.可以说，格子Boltzamnn方法不仅仅是一种数值模拟方法，而且是一项重要的科学研究手段.此外，格子Boltzamnn方法还具有天生的并行特性，以及边界条件处理简单、程序易于实施等优点.可以预计，随着计算机技术的进一步发展，以及计算方法的逐渐丰富，格子Boltzamnn方法将会取得更多成果，并为科技发展发挥更重要的作用.
1.2.1孕育阶段： 对格子Boltzamnn方法发展使得了解，得先从格子自动机说起.格子气自动机使更广泛的元胞自动机在流体学中的应用.元胞自动机是一个时间和空间离散的数学模型.20世纪60年代，Broadwell等人首先提出了离散速度模型，用以研究流体中的激波结构. 20世纪70年代，为了研究流体的运输性质，法国的Hardy、Pomeau和Pazzis提出了第一个完全离散模型，该模型命名HPP模型.这是历史上的第一个格子气自动机模型. 1986年，法国的Frisch、Pomeau和美国的Hasslacher提出具有足够对称的二维正六变形格子气自动机模型，，命名为FHP模型. 由于这些方法在还处在一些缺点： （1）有格子气自动机演化方程推导出来的动量方程不满足Gaililei不变形； （2）流体状态方程不仅仅依赖于密度和温度，还与宏观流速有关； （3）破装蒜子具有指数复杂性，对计算量和存储量也有较大要求.因而，我们将这一段格子气自动机的发展过程称作格子Boltzamnn方法的孕育期.
1988年，McNamra和Zanetti提出把格子气自动机中的Bool运算变成时数运算，格子点上的粒子数不是用整数0或1来表征，而是用实数f来表示系综平均后的局部粒子分布函数，用Boltzamnn方程代替格子气自动机的演化方程，并将该模型用于流体的数值计算.这是最早的格子Boltzamnn模型，从此开启了格子Boltzamnn方法的历史大门. 1989年，Higuera和Jimenez提出了一种简化模型：通过引入平衡分布函数，将碰撞算子线性化.该模型不需要碰撞模型，并忽略各自粒子间的碰撞细节，相比于多粒子碰撞模型，容易构造.同年，Higuera等进一步提出了强化碰撞算子方法，以增加模型的数值稳定性.这两模型统成为矩阵模型. 经历了上述两类模型，格子Boltzamnn方法消除了统计噪声，克服了碰撞算子指数复杂性，但是由于依然使用Fermi-Dirac平衡态分布函数，格子气自动机的其他缺点仍然存在. 1991年，Chen等提出了单松弛时间法，用同一个时间松弛系数来控制不同例子靠近各自平衡态的快慢，进一步简化了碰撞算子；Qian等人在1992年也提出了类似的方法，称之为格子BGK（LBGK）模型.LBGK模型与矩阵模型类似，但与前面两种模型不同的是，当粒子种类数增加时，碰撞算子本身发生生变化，不会变得复杂. 至此，格子Boltzamnn方法完全克服了格子气自动机的一系列缺点，并逐渐成熟，成为国际研究的热点. 早期的格子Boltzamnn模型只能用于等温不可压缩流动的模拟.但因为存在可压缩效应，会引起一定的误差.为了消除或强敌有可压缩效应引起的误差，许多学者致力于新的格子Boltzamnn模型的研究，并提出了多种等温不可压模型.而后，一些不可压缩热模型成功实现了对有效范围温度变化的热力学和传热学问题的模型.其中，最成功的要数双分布函数模型.他是在密度分布函数的基础上引入了温度分度函数、或内能分布函数、或总能分布函数，并用密度分布函数演化得到速度场，这类模型具有与等温不可压模型相同的数值稳定性，而且可以从根本上解决压缩功和耗热问题. 边界处理方面，经历了20年的发展，格子Boltzamnn方法已逐渐发展出适合不同边界条件、不同模型的边界处理格式. 网格划分方面，最初的格子Boltzamnn方法是基于正六边形或正四边形的均匀对称网格.由于均匀网格在计算效率、计算精度等方面的不足，从而促进了非均匀网格、多快以及多重网格、无网格等多技术出现. 总的来说，这些网格技术延展了格子Boltzamnn方法的应用范围，使得格子Boltzamnn方法主机去年从理论的神殿走向更可能多的实际应用领域.
事实上，在20年的发展过程中，格子Boltzamnn方法的确也已成一个十分活跃极具发展前景的模拟手段.并迅速在微/纳米尺度流、多孔介质流、多相多质流、非牛顿流体、粒子悬隔i浮流、湍流、化学反应流、燃烧问题、磁流体、晶体生长等许多领域得到应用.下面分别以多孔介质流、多相流和非牛顿流体三个方面为例，做较详细说明. 由于格子Boltzamnn方法边界条件易于实施，在模拟具有复杂几何构型的问题具有较大的优势，因而这个方向的发展非常迅速.目前，采用格子Boltzamnn方法对多孔介质流进行模拟主要在空隙尺度和代表单元尺度上进行.在孔隙尺度上，可以直接使用格子Boltzamnn方法描述孔隙内的流体流动，多孔介质则当做固体壁面，流体与介质相互作用使用边界处理格式来描述. 在多相流方面，由于真实的流动问题常常是多相的，因而对其开展研究具有重要的现实意义.由于格子Boltzamnn方法的介质特性，它可以方便地描述数流动中不同相之间的相互作用，因而在多相流领域具有较好的应用前景.按照设计方法的不用，现有模拟多相流的格子Boltzamnn模型可分为四大类：着色模型、伪势模型、自由模型和其他模型. 格子Boltzamnn方法在非牛顿流体领域的应用刚刚起步，主要研究对象是非牛顿幂律流体.Aharonov等最早提出使用矩阵碰撞该算子来计算幂律流问题，即在每一个时步内，调整碰撞算自来该表局部的动力学黏性系数.Boek用该模型模拟了幂律流体在简化多孔介质中模型的流动，模拟结果与达西定律符合良好.最近，Gabbanelli又对上述模型进行了改进，引入分段幂律方程描述剪切率和表现黏度的关系. 以上可看出，到目前为止，格子Boltzamnn方法的研究者主要局限在科学界.尽管如此，随着格子Boltzamnn方法理论体系逐渐完善，以及计算机技术的进一步发展，格子Boltzamnn方法也会走向更加广泛的工业实际应用中.
除了求解的困难外，作为一种对流体物理的描述，与描述经典力学运动的牛顿运动方程，或与描述量子力学运动的薛定谔方程等原理方程不同，纳维--斯托克斯方程是从更根本的原理性方程出发，在合理地假定某些物理机制可以忽略后，经过统计平均得到的.本质上纳维--斯托克斯方程当然不可能描述那些被忽略了的物理机制带来的宏观现象，比如流体系统中的相变、非牛顿的本构关系以及在分子运动自由程尺度上的物理现象，在这些领域，纳维--斯托克斯方程明显的显示出了他的局限性. 从20世纪80年代末开始，一种对于流体力学的全新的理论表相及有效的计算方法初步形成，这就是现在人们通常所谓的格子Boltzamnn方法. 关于格子Boltzamnn方法的早期发展，上文已有较全面的综述，在此仅作简单介绍.从历史角度来讲，格子Boltzamnn方法最初是从所谓的格子气模型演化而来的，而后者是一种抽象简化的分子运动数学模型.格子Boltzamnn方法最初的引入有两个主要原因：一是为了降低模型导致的数值噪音；而是能够克服格子气模型里处在的非物理缺陷.可以证明，格子Boltzamnn系统的宏观表象基本满足纳维--斯托克斯方程.从而，人们可以模拟格子Boltzamnn系统地方法来间接地解纳维--斯托克斯方程. 标准格子Boltzamnn方程一般用一下的数学表达式描述： 式中 ——粒子分布函数； ——碰撞项. 用格子玻尔兹曼模型进行流体的数值模拟有一些明显的优越性．如，它的对流（advection）过程是通过常数值速度实现的．这相应的计算是一项极其简单的操作步骤．当适当的格子网格选定后，该过程通常可以用完全平移的方式实现．用计算数学里的常规有限插值语言来讲，它对应于上风插值．但所不同的是其对应的柯郎数（Courant Number）等于1．相比之下，纳维——斯托克思方程的对流项是一个随时空变化的非线性函数．众所周知，对于它的计算不是一项简单的事，并且，数值稳定性的要求迫使人们在实际问题的计算中只能使用比1小得多的柯朗数．在给定空间分辨度的情况下，小柯朗数意味着小时间步长，从而大大延长了计算时间：同时，小柯朗数也增大了数值扩散误差，迫使人们采用更高精度格式或隐式格式．其后果是，或者算法变得极为复杂，并行效率大大降低；或者计算只限制在处理定常流的情况下．事实上，定常流是对流动情况的极大限制．许多重要的流体力学问题，如分离流，即使我们只关心它的时间平均的结果，也是不能用定常流假设来近似的．在此我们也要提一下格子玻尔兹曼方程的另一个本质特性：所有非线性效应在格子玻尔兹曼方法里都包含在碰撞项中，并且是以纯粹局部信息的方式体现的．这进一步发挥了并行计算的长处．所有这些理由意味着格子玻尔兹曼方法是对非定常流动实行大规模并行模拟计算的一种比较优越的方法． 相比之下，以流体力学方程（纳维一斯托克思方程或Burnett类型方程）宏观描述为基础的传统计算方法对许多这类问题存存基本困难．除边界条件之外，利用各种封闭性假设推导出的超越纳维一斯托克思的宏观方程直至现今仍存在对其数学规范性的疑问和争议，多相流的计算也存存同样问题．众所周知，流体系统中存在多相的物理机制是分子问的长程作用力，这种机制早已超出了流体力学方程所能描述的物理现象范围．以流体力学方程为基础的多相流计算方法必须依赖额外的模型来模拟流体力学方程本身所不包含的物理现象．除了实际数值结果显示的问题之外，这种方法本质上隐含着严重的基本物理缺陷，这种缺陷集中表现在对相交界面的准确描述上面，即在十分尖锐的相界面附近，纳维一斯托克思方程之类近平衡态的近似表象是有相当疑问的．这也反映在相界面和兀滑动（no—slip）固体边界条件的互斥性上面，为了修补这一缺憾，人们不得不引入各种滑动经验模型.反之，以细观（mesoscopic）为表象基础的格子玻尔兹曼方法可容忍更大的非平衡态程度及更广义的严格边界条件．另外，压力的状态方程在细观表象中是由粒子的相互作用自然得出的，而不用直接输入和处理．在相变情况下，物体的宏观特性将产生不连续性，而对应的微观和细观力学机制并无改变．格子玻尔兹曼方法在模拟多相流上有着广泛的使用. 然而，这种为大多数人所熟悉的格子玻尔兹曼方法的理论框架存在本质上的缺陷．由于它运用逆向切普曼一安斯柯格展开的途径来适定平衡态分布函数中的关键参数，以达到复建宏观物理体系的目的，这就使其对比纳维一斯托克思方程更为广泛的流体物理问题变得束手无策．因为后者通常缺乏明确的宏观表述方程.除了这一限制外，大多数通常所知的格子玻尔兹曼模型还存在着其他明显的局限，如它的离散速度集只能达到四阶各项同性的要求，而其平衡态分布也只能适用于极低马赫数下，温度近似为常数的情况，这就限制了格子玻尔兹曼方法在可压缩流动模拟中的应用．在这里我们也顺便提一下和此有关的在学术领域中的思想谬误．这种谬误把格子玻尔兹曼方法和某些通常意义下特定的格子玻尔兹曼模型混为一谈，从而得出该方法只能适用于处理近似不可压及恒温纳维一斯托克思流体问题的错误结论.
近年来，我们在格子Boltzamnn方法的理论和应用两方面都展开了大量的研究工作.本论文除介绍格子Boltzamnn的基本理论原理以外，也包括了我们在格子boltzamnn方法领域所取得的部分研究成果，以及在相关科研实践中积累的经验.我们将它们整理与此，愿与国内外读者分享交流. 本文选择了应用最为广泛的标准模型（D2Q9）、He-Luo模型以及D2G9模型 来进行比较.在理论上的分析以后，具体模拟了二维方腔流和压力梯度驱动的 Poiseuille流动，从理论和实际两个方面来具体分析了这几种模型的优缺点以及使用范围. 因此，本论文的定位是尽可能地全面的介绍格子Boltzamnn方法的发展现状，并系统阐明该方法的流体力学与统计力学背景知识、模型推导、边界处理、数值实施等问题.我们认为，只有牢固掌握相关基础知识、并有着全局的视角，才能做出更出色的研究成果.
2000年，wolf-gladrow所著Lattice-Gas Celluar Automata and Lattice Boltzamnn一书出版，这是一支的国际上第一本格子Boltzamnn方法方面专著.该书由格子气自动机讲起，内容浅显易懂，适合新手热门学习.2001年，succi所著的Lattice Boltzamnn Equation for Fluid Dynamics and Beyond一书出版，该书涉及面广，对格子Boltzamnn方法的理论，在流体力学、热力学以及量子力学领域的应用，格子Boltzamnn方法适用领域等都做了讨论，影响了不少学者.2005年，Sukop和Thorne完成了Lattice Boltzamnn Modeling一书.我国学者郭照立等人于2002年所著的《流体力学的格子Boltzamnn方法》一书，是目前国内少有的相关专著.此外，也有日文所著的相关书籍. 目前，格子Boltzamnn方法领域已有两个专门的国际会议：DSFD以及ICCMMES.并发展出一个专门的商业软件PowerPlow，其应用结果甚至在世界顶级杂志Science上发表在图1.1中，我们给出了这20年来，不同年限里SCI所收录的格子Boltzamnn方法论文数.由此可见，逐年所发表的相关论文，易近似成指数增长.在我国，格子Boltzamnn方法的研究起步稍晚，但最近几年，论文发表数量已经在国际上占有重要分量.
图1.1 我国在论文近年发表状况 第2章 格子Boltzamnn方法介绍
本章从Boltzamnn方程和元胞自动机的理论出发，系统的阐述了格子气自动机和格子Boltzamnn方法的理论思想.并对格子Boltzamnn方法的由来以及国内外研究现状作了详尽的总结.
在对流体系统进行描述的时候，一般的连续建模思想假设流体分布于它所占据的区域，相应的物理量（流速、密度）均是在时间和空间上满足一定光滑性的函数.从最基本的质量、动量和能量守恒律以及空间不变性（平移、转动和伽利略不变性）可以导出流体的运动方程[1]，这些方程就是如下的连续方程： （2.1） 和Navier-Stokes方程： （2.2） 同时我们还知道：Reynolds数：.它影响方程（2.1）的非线性和稳定性.流动在高Reynolds数下发生非定常流和强湍流现象.Mach数：.Mach数决定流体的可压缩性.低Mach数流动通常看作是不可压的. 一般来说，使用Navier-Stokes方程描述流体的运动是比较方便的.但是，在某些情况下，求解这一组偏微方程组是非常困难或者说是根本不可能的[5].还有一些系统，我们甚至不知道其数学模型，更不用说进行数值模拟了. 对于这些用宏观方程难以描述的系统，使用微观的分子动力学方法或者介观的动理论方程进行描述更为恰当和可行.但是，这两类方法又会带来其它问题，分子动力学需要跟踪多个粒子的运动，计算量非常之大，对计算机的存储量和计算速度都有很高的要求. 运动层次论的描述是一个微观动力学模型，从流体的微观结构出发，运用非平衡统计物理的观点，一切宏观特征都看作是流体分子作混乱运动的结果.在这个微观模型中，基本单位是流体分子，它们的运动遵循物理守恒律；基本方法是统计方法.得到的基本方程是Boltzamnn方程．
（2.3） 上面的方程也叫做输运方程，它描述的是稀薄气体分子的动力学规律.值得注意的是在这个模型中流体不再是连续介质，而是由大量离散的分子组成.与连续介质模型比较，运动论描述更多地考虑了流体的物理属性.运用统计物理的方法可从Boltzamnn方程过渡到宏观流体动力学方程.流体的运动论描述是格子Boltzamnn方法诞生的物理背景. 2.2细胞自动机（CA）
一般的，一个细胞自动机可以定义为一个四元组.是空间点阵，称为细胞陈列，它可以是一维或多维空间的点阵，可为有限或无限；是细胞可取状态的集合，称为状态空间；为相邻细胞的集合；是状态的演化函数. 对任意的细胞，在时刻有，在下一时刻，它的值由邻域在时刻的状态确定，即.根据此式，细胞自动机便随着离散的时间演化下去. 细胞自动机有着显著的三个特征：大规模同步并行，局域相互作用和细胞结构简单.同步并行指在各个格点的演化是同步进行的.细胞自动机适用于不同结构的并行计算机系统，尤其适合于分布式计算机上.细胞自动机的计算区域可分为一些独立的子域，在各自的子域里独立的演化，在各子域之间只需交换共同边界上的数据.这对于解决大空间的数学物理问题效果极好.细胞的演化规则一般都比较简单，但其演化能显示复杂的物理行为.这引起了物理学家、化学家和生物学家的广泛兴趣. 细胞自动机已广泛用于模式识别、机器自组织、信息恢复、演化理论、物理、 化学和生物系统的模拟及其并行处理等领域.例如：一个偏微分方程问题的差分算法可把空间划分的格子看作细胞，求解过程则可近似看作一个细胞自动机的演化；细胞自动机可以用于非线性化学系统的反应扩散方程的模拟许多生物系统可用于细胞自动机模拟.
格子气自动机则是一个更大胆的尝试，它不仅对流体的空间离散，而且时间也是离散的.这样就更便于计算机进行处理，为计算流体力学提出了一个全新的方向. 格子气自动机方法模拟流场，是直接从Boltzamnn方程出发作完全离散处理，将流体存在的空间划分成离散的网格（如正方形网格，正三角形网格等）；流体想象成由大量只有质量没有体积的微小粒子组成；时间也离散成整时间步（=0，1，…，n，…）.流体粒子存在于网格节点上并沿网格线运动.所有粒子根据一定的规则（称为碰撞规则）同步地随着时间步相互碰撞和移动（称为时间演化）.这里所谓“同步"，就是要求在每一个整时间步，所有粒子都处于相应的节点上，并且同一个节点上的粒子之间发生碰撞.在下一时刻，粒子正好运动到邻近的节点，又和其他到达该节点的粒子发生碰撞，如此等等.碰撞规则指的是粒子间相互碰撞和如何运动（碰撞后该向哪个方向运动）的准则. 有了以上的一些规定之后，就建立了一个格子气模型，也叫做格子气自动机. （2.4） 式中 ——表示在时刻，节点上以的速度运动的粒子数，它是个整数； ——碰撞算子，表示由于粒子间的碰撞所引起的速度为的粒子数目的变化. 式（2.4）的具体物理意义是，到达同一格点处的流体粒子之间发生碰撞，各个粒子的运动方向可能发生变化.随后，粒子按照其运动速度流动到相邻格点.其中，δ是离散时间步长. 具体地说，它一般应包括以下几个基本要素：速度离散，空间离散，时间离散和碰撞规则. 利用格子气方法求解流体力学问题大致的求解步骤是： （1）建立格子气自动机； （2）运行格子气自动机，即让流体粒子在离散网格上，根据碰撞规则随着整时间步演化； （3）统计平均求数值结果.这里的统计平均一般包括时间平均和空间平均两个过程. 这里我们可以看出格子气方法与传统的数值方法相比较，完全是一个逆向的 求解过程，而且在许多方面优于传统的数值方法.
自从1986年Frisch，Hasslacher和Pomeau（FHP）提出第一个二维的离散格子气模型以来，格子气自动机作为模拟复杂物理系统的重要方法，日益受到人们的关注.尤其是对流体现象的成功模拟，引起国内外专家的高度重视. 由于采用布尔型变量，LGA模型不会遇到数值不稳定问题；同时，所有粒子的碰撞和流动同时发生，且粒子之间的相互作用是局部的，因而LGA具有本质的并行性，非常适于在并行计算机上实现；在固壁边界上，流体粒子反弹回流场内部，因而LGA在处理复杂边界时非常方便.以上这些都是LGA模型不同于传统数值方法的优点. 但是，格子气自动机也有以下几个不足之处： （1）统计噪声.由于LGA的碰撞算子含有随机因素，因此不可避免的有噪声的影响.虽然可以通过时间平均或空间平均的方法降低噪声成分，但噪声的影响还是很大，并且所需要的计算和存贮量也加大了. （2）算子的指数复杂性.LGA的碰撞算子与离散方向数成指数关系，不但增大了LGA模型的设计难度，而且不利于LGA的应用.这一问题对三维情况尤为突出. （3）不满足伽利略不变性.与LGA对应的宏观方程中，对流项前面有一个非单位的因子，因此该方程不满足伽利略不变性的要求.虽然通过重新标度可以得到正确的Navier-Stokes方程，但这种方法只能用于简单系统的LGA模型，对一些复杂的LGA模型（如多相流模型）来说这种方法是不可行的格子Boltzamnn方法就是为克服LGA的这些不足而发展起来的.在格子Boltzamnn方法中，用粒子分布函数的演化代替LGA中粒子的演化，其演化方程直接采用格子Boltzamnn方程，并根据分布函数直接计算流体的密度和速度，因此消除了统计噪声. LGA的统计噪声可以使用空间或时间平均方法加以克服，但这需要较大的计算和存储代价.为了消除这种噪声，McNamara和Zanetti（1988）提出了直接使用平均粒子数或粒子分布函数代替布尔变量进行演化，从而开创了流体计算的一个新方向.McNamara和Zanetti提出的这个模型也就是最早的格子Boltzamnn模型.在该模型中，粒子分布函数按照与LGA类似的方式进行演化：
该方程就是格子Boltzamnn方程.其中碰撞算子与LGA的碰撞算子具有相同形式，但布尔变量用平均粒子数代替.由于直接使用了实数型的粒子分布函数代替布尔型的粒子数进行演化，MZ模型克服了LGA方法的噪声.但是，这一模型采用的碰撞算子仍然具有指数复杂性，也不能克服LGA的其他缺点.MZ模型提出不久，Higuera和Jimenez（1989）就对其做了改进，本文称之为HJ模型.他们证明了上述计算量很大的复杂碰撞算子可以用一个线性算子近似，其方法是引入一个平衡态分布函数，并假设： （2.5） 式中 ——为非平衡态部分. 同时对碰撞算子进行泰勒展开，可以得到： （2.6） 式中 ——碰撞矩阵. 虽然HJ模型的碰撞矩阵和平衡态分布函数仍然依赖原LGA模型，但其计算复杂度比MZ模型的指数复杂度大大降低，大大提高了计算效率.HJ模型首次在演化过程中使用平衡态分布函数，为格子Boltzamnn方法后来的发展奠定了基础.从这一点来说，HJ模型是格子Boltzamnn方法研究历史上的一个重大突破.HJ模型的一个主要缺点是数值稳定性不够好. MZ和HJ模型是LGA模型的直接发展.它们克服了LGA方法的统计噪声，HJ模型的计算复杂度也大大下降.但是，这两种格子Boltzamnn方法都与基本的LGA模型有关，碰撞项都来源于LGA的碰撞规则，平衡态分布函数本质上都是Femi-Dirac分布.这些特点限制了它们的适用范围. 在HJ模型提出不久，Higuera，Succi和Benzi构造了一种新的格子Boltzamnn模型[22]，本文称之为HSB模型.该模型完全抛弃了LGA的Fermi-Dirac分布，而使用Boltzamnn分布，同时根据需要设计适当的碰撞矩阵，以得到正确的宏观方程. 所以HSB又称为强化碰撞模型.设计HSB模型的关键是构造出相应的平衡态分布函数和碰撞矩阵.与HJ模型不同的是，HSB模型的碰撞矩阵与LGA模型无关，矩阵的元素是一些待定参数，根据所要满足的宏观方程的需要确定. HSB模型的碰撞矩阵一般是满矩阵.1991年前后，几个不同的研究小组各自独立的提出了一种单松弛模型.该模型的碰撞矩阵是一个对角矩阵： （2.7） 式中 ——待确定的平衡态分布函数； ——无量纲的参数，表示粒子分布函数达到平衡态的松弛时间. 所以单松弛模型也称为LBGK模型.在LBGK模型中，确定平衡态分布函数最为关键.至少要满足如下的质量和动量守恒条件： （2.8） LBGK模型极大地提高了计算效率，并且只要选择恰当的平衡态分布函数，从该模型可以导出正确的Navier-Stokes方程.LBGK模型至少在目前是格子Boltzamnn方法研究和应用中最主要的模型. LBM方法不但是LGA方法的发展，而且还可以看作是连续Boltzamnn方程的LBM方法不但是LGA方法的发展，而且还可以看作是连续Boltzamnn方程的一种特殊的有限差分格式，甚至可以看作是求解守恒系统的一种Lagrange格式. 与格子气方法相比，格子Boltzamnn方法用并行计算更有效，同时克服了格子气方法的缺点.和其他数值计算方法相比较，格子Boltzamnn方法具有以下特点： （1）与流体宏观描述中非线性对流项相对照，在以流体的分子运动论的描述为基础的格子Boltzamnn方法中，相对于位相空间的对流过程是线性的.使用Chapman-Enskog展开可从格子Boltzamnn方程导出宏观的非线性对流过程； （2）N-S方程中，压力作为速度的来源用运动方程式来表示.在这种情况下，用有限差求解N-S方程时，必需采用反复计算和松弛法处理.格子Boltzamnn方程中压力以状态方程式表示.其解法无需反复计算的特征使其可使用并行计算； （3）在使用Maxwell平衡分布的运动论中，为了使位相空间在连续空间存在，在求解平均量的过程中需要知道所有位相空间的情况.而在格子Boltzamnn方法中，以最小的离散化速度来代表导出N-S方程，并可从局部粒子分布函数统计得出宏观量； （4）格子Boltzamnn方法存在最适于并行处理的局部相互关系模型； （5）与格子气方法相比，统计噪声不复存在； （6）边界条件容易设定. 格子Boltzamnn方法（以及格子气方法）为计算流体力学提供了新的思路.格子气和格子Boltzamnn方法自提出以来得到了飞速的发展.在流体动力学，多孔介质渗透流动，磁流体动力学，传热问题，多相流，MEMS和化学反应等多方面得到了应用.格子气权威G.D.Doolen指出，对于一些新的物理现象，特别是那些尚无法用宏观方程描述的物理现象格子Boltzamnn方法将是一个强有力的工具. 由于格子Boltzamnn方法的影响的扩大，其在国内也受到日益的重视.从上世纪八十年代后期朱照宣、钱跃竑等人在《力学和实践》等杂志上介绍了格子气方法以来，华中科技大学、清华大学、北京应用物理与计算数学研究所、吉林大学和武汉大学等单位都在不同的方向开展了富有成效的工作.国家自然科学基金委的力学发展战略报告中指出，研究多相流问题，格子Boltzamnn方法可以对某些流动图案或机制做出定性的说明，随着计算机技术和并行算法的飞速发展，用这个方法进行定量计算的日子不会太遥远.
（2.9） 式中 ——格子上的一个格点； ——流体粒子的离散速度集合，为离散时间步长，｛=1，2…，｝； ——为当前时间步； ——速度c运动的速度分布函数，是碰撞算子，表示分子间碰撞对速度分布函数的影响. ——一个自封闭系统， 即如果则，所以格子Boltzamnn方法中流体粒子总是在网格线上运动.流体宏观物理量及密度、速度和内能有离散函数的速度矩得到 （2.10） 式中 ——维空间位数. 演化方程又称为格子Boltzamnn方程，它可以看作是LGA演化方程的统计平均，也可以看作是连续Boltzamnn方程的一种特殊离散格式.LBE中的碰撞算子反应了微观流体粒子的相互作用，因此对模型能否准确刻画出流体系统的物理规律起着关键作用.目前LBE一般采用线性化碰撞算子，即 （2.11） 式中 ——一个的矩阵； ——依赖宏观物理量的平衡态分布函数. 因此，碰撞矩阵和平衡态分布函数完全决定了碰撞算子的特性，也决定了一个LBE模型刻画流动问题的能力和精度.当前的LBE模型主要是为解决小Mach数近似不可压Navier-Stokes方程设计的，但也有学者研究可用于大Mach数近似可压流动的LBE模型. 从方程可以发现，格子Boltzamnn方法是一种不同于传统数值方法的流体计算和建模方法.从离散的网格来说，LBE具有Euler方法的属性；从离散的例子观点来看，LBE又具有Lagrange方法的特点.LBE的这种特性使其具有如下的几个显著特点：首先，虽然LBE模型主要用于模拟宏观连续流动，但他不是基于连续模型的Navier-Stokes方程，而是基于介观模型，本身没有连续介质条件的假设.因此，只要设计得当，可以用于描述尺度、稀薄流等连续流动问题.其次，格子Boltzamnn方法的微观例子背景使得它可以比较直观、方便地处理流体内部以及流体与周围环境的相互作用，从在多组分、多相态系统、界面动力学、渗流等复杂流动现象的描述方面，表传统的数值方法更有优势.最后，从计算的角度看，LBE的演化过程物理清晰，计算简单，容易编程，并且计算局部的，具有良好的并行性和可扩展性，对大规模流动问题的计算具有很大优势.
机和格子Boltzamnn方法的理论思想.并对格子Boltzamnn方法的由来以及国内外研究现状作了详尽的总结. 本章是全文的理论基础.通过本章讲述的基本概念，在接下来的一章，我们 首先介绍了实际应用中所要用到的几种基本模型，然后通过对二维方腔流以及压力梯度驱动的Poiseuille流动的模拟，从计算精度、收敛速度、数值稳定性等几个方面对其进行了比较详尽的比较.希望通过我们的工作能够为工程应用时模型的具体选择提供一定的依据. 第3章 格子boltzamnn方法的基本模型比较
前面所提到的这些基本模型都在不同的方面有着自己的优势，从而导致了这样的一个问题：面对实际问题到底应该选择何种基本模型，以及为什么要这样选择.这正是本文的工作重点.在 2002 年，一位日本学者也做了类似的工作，不过他仅仅只是从收敛速度方面对不同模型进行了一个简单的比较，并得出了 D2G9 模型收敛速度最快的一个结论. 本文将在查阅大量文献资料的的基础上，进行了认真的调研后。从计算精度、数值稳定性和收敛速度三个方面出发，系统的对几种应用最为广泛的格子 Boltzamnn 方法基本模型作了理论上的分析和数值上的模拟验证工作.通过我们的讲述，希望能够为工程应用时选择正确的基本模型提供一些依据.
1992 年，Qian 等人提出了 DnQb（n 维空间、b 个离散速度）系列模型，这些模型到目前为止仍然应用的最为广泛.本文将选择 D2Q9 模型来进行介绍和比较. D2Q9 模型的离散速度： （3.1） 当i =0 时，=0 ； 当i =1～4 时， 当i =5～8 时，=. 为网格速度.，分别为网格步长和时间步长.平衡态分布函数： （3.2） 其中 =/3，权系数ω=4/9，=1/9（i=1-4），=1/36（i=5-8）. 演化方程： （3.3） 宏观方程和宏观速度： （3.4） 可以推导得出如下的宏观方程： （3.5） （3.6） 式中 =——压力 ——粘性系数. 当流体密度变化不大时，以上方程就是标准的不可压 Navier-Stokes 方程. 3.1.2 Zou-Hou 模型 1997 年，Zou 及其合作者提出了一个计算定常不可压 Navier-Stokes 方程的一个格子 Boltzamnn 模型，我们称之为 Zou-Hou 模型. 对于定常流动，使用 D2Q9 模型时可以得到如下宏观方程： （3.7a） （3.7b） 为了便于比较，我们给出常密度时定常不可压 Navier-Stokes 方程： （3.8a） （3.8b） 比较以上方程可以发现，如果要使方程组（3.7）逼近方程组（3.8），必须在（3.7）中忽略包含有密度的空间导数的项（压力项除外）. Zou 等人发现，如果仍然采用标准的格子 Boltzamnn 演化方程，但是选取如下的平衡态分布函数： （3.9） 并定义流动的宏观密度和宏观速度为： （3.10） 利用多尺度展开方法可以得到如下的宏观方程： （3.11a） （3.11b） 式中 ——压力，=； ——粘性系数，. 在定常情况下，方程组（3.11）与标准的常密度定常不可压 Navier-Stokes 方程组完全一致，从而消除了可压缩性误差.应当指出的是，Zou-Hou 模型针对的是定常不可压流动，不适用于非定常的情况. 3.1.3 Chen-Ohashi 模型 1997 年，Chen 和 Ohashi 将 Zou-Hou 模型推广到非定常不可压问题，建立了一个新的不可压格子 Boltzamnn 模型，本文称之为 Chen-Ohashi 模型.为了消除由于模型的可压缩效应所导致的人工压缩效应，Chen 等人将方程（3.10）所定义的宏观速度u 看作一个临时速度，流体的真实速度w通过对u 进行校正来得到： （3.12） 其中是一个外力项.从而可以得到如下的平衡态分布函数： （3.13） 宏观密度和宏观速度定义如下： ， （3.14） 利用多尺度展开技术，可以得到 Chen-Ohashi 模型的宏观方程： （3.15a） （3.15b） 式中 p——压力， ——粘性系数， 3.1.4 He-Luo 模型 在 Chen-Ohashi 模型提出的同时，He 和 Luo 提出了另一个针对一般不可压流动的格子 Boltzamnn 模型，本文称之为 He-Luo 模型.这一模型的基本思想是直接在平衡态分布函数中消除由于密度变化引起的高阶 Mach 数项.在 He-Luo 模型中，独立的动力学变量是压力而不是密度. 将密度分解为，从而可以得到 （3.16） 式中 ——平衡态分布函数中与速度有关的部分： （3.17） 我们略去上式中关于的项，可以得到一个全新的平衡态分布函数： （3.18） He-Luo 模型在不可压 Navier-Stokes 方程中使用压力作为一个独立的变量，从 而引入了一个压力分布函数及相应的平衡态分布函数： （3.19） 式中通过的演化方程我们可以得到 的演化方程： （3.20） 基于这种压力分布函数表示方法，流动的宏观压力和速度确定如下： （3.21） 通过多尺度展开方法，同样可以得到 He-Luo 模型的宏观方程： （3.22a） （3.22b） 式中 . 3.1.5 D2G9模型 前面介绍的这几个格子Boltzamnn模型在一定程度上都降低或消除了原格子Boltzamnn模型的由可压缩效应引起的误差.但是我们也注意到，这些不可压格子Boltzamnn模型还存在一定的局限性.Zou-Hou的模型在理论上只能用于定常流动，而不适用于非定常的不可压流动；Chen-Ohashi模型虽然可以模拟一般的不可压流动，但在计算的每一个时间步都需要解一个Poisson方程，从而增加了较大的计算量，并部分地丧失了格子Boltzamnn方法的优势；He-Luo模型本质上仍然是求解不可压Navier-Stokes方程的一个人工压缩方法，要忽略人工压缩性造成的影响，需要增加额外的限制条件，限制了方法的应用范围.最近郭照立等人构造了一个能够模拟一般不可压Navier-Stokes方程组的格子Boltzamnn模型，在一定程度上克服了以上几个模型的不足，我们记为D2G9模型. D2G9模型的基本思想是构造一个格子Boltzamnn模型，使其平衡态分布函数 满足，从而就可能得到一个针对一般不可压Navier-Stokes方程的模型. 在该模型中，引入了一类新的分布函数，其平衡态分布函数定义为： （3.23） 式中 常数——流体的平均密度， ， 、和——一些模型参数，并满足条件 模型的演化方程为： （3.24） 流体的宏观速度和压力定义为： （3.25） 在演化的过程中，要求碰撞过程保持质量和动量守恒条件，即 （3.26） 同样的通过多尺度展开技术，我们可以得到D2G9模型的宏观方程： （3.27a） .（3.27b） 我们可以看出，D2G9模型可以推导出完整的Navier-Stokes方程组. 对这几个基本模型做了一个简单的介绍以后，下一步我们就将利用它们分别来模拟二维方腔流动以及压力梯度驱动的Poiseuille流动. 3.2 进行常压力梯度驱动的Poiseuille流动模拟比较几种基本模型
Poiseuille流动 平板间的粘性层流运动虽然简单，却能够揭示粘性流动的一些本质特征.所以我们首先对有压力梯度驱动的平板流进行了模拟. 如图3.1所示，两平板间充满了粘性系数为ν的流体，在流场的出口、入口处施加一恒定的压力梯度.我们知道理论上可以得到其流场速度的解析解： （3.28） 式中 ——管道宽度； ——流体的粘性系数； ——压力梯度.
图3.1 Pliseuille流流场示意图 模拟时，使用压力边界的二维Poiseuille流动的区域定义为0 ≤x≤2，0 ≤y≤1；流场的格子划分为；流场内的初始化速度为0，初始化密度为1.0.定义流动的雷诺数，.其中；边界处理我们均采用郭照立等人提出的非平衡态外推方法.在上述初始和边界条件下，经过若干次迭代后流动达到了稳定状态. 为了比较几个模型的精度和收敛速度，我们主要是取不同的雷诺数来进行模拟.具体就是模拟了Re=100、400、1000和1500这几种情况，然后把模拟结果与Poiseuille流动的解析解进行比较. 首先，我们分别做出了Re=100、400、1000和1500这几种情况下速度水平分量u在管道中心线x=1.0处的图形.
图3.2 Re=100时速度水平分量
表1：速度水平分量误差比较表
相对应的可以画出相对误差图见图：
另外一个方面，我们利用计算时间来比较了一下这几种模型的收敛速度.我们定模拟的时候流动达到稳定状态的标准是： （3.29） 式中 ——时刻在点处的计算值，. 以上几种模型在相同条件下模拟Poiseuille流动时，计算时间都在2个小时左 右，达到收敛时计算步数都在14万步左右.随着Re数的增大，计算时间和计算步数也有所增大，不过最大也不超过3个小时和20万步（Re=1500）.D2G9模型的计算效率仍然是最快的，不过和精度的比较一样，优势体现的不够明显. 同时，我们在Re=1500的情况下，变化不同的松弛时间来讨论各种模型的数值稳定性.我们可以得到如下这个计算结果：
Re=1500
从上面这个表格可以看出，D2G9模型毫无疑问的拥有最为广泛的数值稳定性.
从模拟结果我们可以看出，在雷诺数比较小的时候，这几个模型都能够得到比较好的模拟结果，无论是计算精度还是收敛速度等各个方面.各自的相对误差最大也不超过5%，都具有工程应用的价值.此时，D2G9模型虽然具有相对更好的计算精度和收敛速度，但是优势体现的并不是特别明显.不过，我们可以看出这样的一个趋势，随着雷诺数的增大，模型的相对误差也随之增大，当雷诺数超过5000（Re=10000）的时候，He-Luo模型和D2Q9模型都出现了数值不稳定和计算结果发散的情况.相对应的，D2G9模型仍然能够得到比较好的模拟结果，相对误差也在工程允许的误差范围之内. 从而，我们可以得出这样的结论：在实际应用时，如果模拟对象的雷诺数比较小，这三种模型基本可以通用，当雷诺数增大到一定的数值时，D2G9模型就体现出了他在计算精度以及收敛速度等方面的优势，可以得到很好的模拟结果. 第4章 格子boltzamnn方法的算法设计
平衡态分布函数满足质量和动量守恒： （4.2a） （4.2b） 宏观的密度和速度为： （4.3a） （4.3b） 给定一个问题的初始条件和边界条件，选用格子Boltzamnn模型和边界处理格式，经过若干次的演化迭代就可以得到该问题的数值解.格子Boltzamnn方法的一个很大的优点就是算法简单，程序代码简洁，核心代码只有数百行，但是可以求解复杂的物理问题.
图4.1:计算流程图 用格子Boltzamnn方法模拟二维间题的程序主结构如下， { alloeate（）；//动态分配数组（f，F）和（u，v，p） init （）;//流场初始化 for（t=1;t<max-t;t++）//max-t为演化步数 { evol（）;//流场内部各个格点各个方向的演化计算 exboulld（）;//流场的边界处理一非平衡态外推方法 } Write-data（）;//存储计算结果:速度和压强 }
点就是以简御繁.正是迅速发展的计算机技术使得复杂流动的大规模模拟计算成为可能.尤其是随着超级计算机的应用，复杂流动模拟的精度和复杂度提高很快，例如对复杂几何形状的三维N一S流的模拟.对于高性能高精度计算，很重要的途径就是网格细化，缩短时间步长，但是这会带来计算量的剧增.一个间题的解决需要数天甚至更长的时间，这样就不能及时得到问题解，影响对数据的后处理，导致对数据的分析延迟.诚然，随着计算机硬件的发展，人们可以选择性能高的并行计算机.但是对算法程序的优化也是至关重要的.优化的算法程序可以提高计算效率数倍，应用优化技术，通过使用程序分析工具，采用算法优化、代码实现优化和编译优化等手段，可以大大优化应用程序的性能.因此对格子Boltzamnn方法程序的优化也至关重要，它直接影响着该方法的应用.目前国外已经开展了相关工作的研究.在本节中，我们在已有的文献基础上，针对不可压模型和非平衡态边界处理格式做了详细的分析和实验数据.另外，我们主要讨论串行程序的优化，对于并行的程序设计，可以直接应用到并行的环境中. 4.2.1优化算法
函数 feq（） evol（） int（） exbound()
从表4.1中可以看出，占用时间最多的就是feq（）函数和evol（）函数，是程序的主体部分.程序的瓶颈（Bottleneck）和热点（Hotpoint）就是演化函数和平衡态分布函数.我们主要做了三个方面的优化. 首先，演化可以分为两步:碰撞步和流动步，在碰撞函数中，平衡态分布函数的计算量很大.注意到计算式，反复调用了数组e[Q][2]，而该数组的值极其简单（0，1，-1），所以在计算中直接代入，减少了数组的调用.另一方面，在平衡态分布函数的计算中利用格子方法的对称性和相关性，比如，
这样可以减少浮点数的运算，而且减少了一层循环（k）.其次，流动函数只是简单的赋值，但是占了36.864%的计算时间，这种简单的赋值显然不优.在流动函数中，注意到格子方法的方向性，分布函数如果沿方向更新赋值，如沿方向的流动为 f[1]][y][x]=f[1][y][x-1] 就不需要F数组存储下一个时间步分布函数的值，换句话说，各个网格点只需要一个分布函数就可以计算和存储下一时刻的分布函数的值，完成演化计算.另外对于存放各个网格点的宏观量数组也占用大量的内存空间，在程序中只是在计算平衡态分布函数的时候用到，但是这样耗费了大量的空间内存，因而在计算中需要宏观量的地方直接用分布函数计算而不是计算好存放在数组中调用，这样减少了空间内存，程序分配的内存空间只是原来分配空间的9/21，从而在同样的硬件条件下，可以模拟更大规模的流体间题.由于流动步按方向流动，为了减少寻址的时间，因此分布函数的动态分配改为f[Q] [max_y] [max_x]. 最后，在边界处理中，针对非平衡态外推的方法[97]做了优化.边界点的分布 函数计算公式为 （4.4） 我们知道平衡态分布函数的计算占用了大量的时间，从边界处理的格式，我们调整顺序： （4.5） 从平衡态分布函数的计算式，在计算时只需计算内点和边界点的平衡态分布函数的差值，由于边界点的压强是由内点的压强近似得到的，密度前后抵消，这样减少了浮点数的运算另外，我们也做了inline、减少catche missing等方法优化程序.格子Boltzamnn方法具有天然的并行性，各个格点的运算只和它邻近的节点 相关，可以设计阻塞通信和非阻塞通信两种并行程序，以及多向划分的并行算法. 但是对于并行算法每个计算节点的程序核心代码和串行程序的相同，因此，基于我们优化的串行程序再并行化，将远远节省并行计算的计算时间.
图4.2方腔流示意图 为了使得我们所得的数据具有稳定性，可靠性，可比拟性，我们串行程序的计算机模拟均是在同一P4 3.0GHz，内存1G的单机上实现的，实验环境如下：
CPU Arehiteeture:IA32 （Family15;Model4:Stepping l） CPU Speed（MHz）:3000 Number of CPUs:2logieal;1physieal Colleetion Duration（sec）: 20 OS: Mierosoft、WindowsXP Professional CPU L2 Caehe Size（KBytes）:1024 Bus Speed（MHz;user seleeted）:533 Bin Size:Funetion Level 全部程序源代码采用ANSI-C编写，变量均采用双精度型.用Intel的Compiler9.0编译工具进行编译，得到可执行文件.exe，参数选取/Od，/O2和/fast.
图4.3 三种算法的运行时间对比
经过我们的算法优化运行速度提高了2.877倍，对原始程序用Intel的编译器选参数/fast，运行时间和我们对算法的优化几乎相当.都采用Intel的编译器选参数/fast，优化的程序仅提高了0.957倍.优化的程序，对于编译器的各种参数运行时间相对比较稳定.总之，我们对算法的手动优化和编译器的优化，使得运行速度比没有优化的算法提高了7.264倍. 对于模拟更复杂的物理问题，如后台阶流、绕流、自然对流、磁流体等，格子方法也同样可以扩展到相关领域，程序的实现也不复杂.对于更复杂的三维问题，只是空间位置由二维变为三维，LBM模型的参数不一样罢了.但是计算量随着网格数的增大而剧增，优化的程序（时间和空间）使得人们探讨三维流体问题成为可能.我们将二维的方腔流扩展到三维，分别用原始算法和优化算法做数值对比实验，不同网格不同编译参数下程序的运行时间（单位:秒）见表4.2.
/Od /O2 /fast =1.8 优化算法 33.547 23.109 14.547 =1.8 优化算法 115.984 80.844 47.41 =1.8 优化算法 338.250 267.219 119.453 =1.8 优化算法 855.265 591.922 361.297
原始算法和优化算法的运行速度比如图4.4所示：
图4.4不同网格三种参数下的运行时间比 从模拟的结果看，对于不同的网格，不同的编译参数，优化后算法的运行速度都有大幅提高，细网格的计算情形下，运行速度提高更为显著，体现了优化步骤中空间内存节省的好处.从表中可以看出，在单机上用 969696的网格计算三维方腔流问题，最优化的程序运行10，000步也只需一个小时.借助优化的程序和编译器的优化作用，在单机上就可以方便地研究三维问题.
第5章 格子Boltzamnn方法的实际应用
格子Boltzamnn方法（Chen 和Doolen 1998）最近已经发展成为对流体流动模拟的一种替代方法.在格子Boltzamnn方法粒子分布函数中，在点x时间t处只限于在正则格上同步移动.分布函数在保持质量，动量格，各项同性和伽利略不变性的同时，对格子有所影响.在这里，表示格子是否与分布函数相链接.在本文中使用的格子是D2Q9，如图5.1.1所示.
（5.1） 如格子D2Q9，见图5.1.1. （5.2） 式中 ——碰撞算子. 在每个节点处，由 和 （5.3）流体密度和速度u可以直接由分布函数计算得到.假设分布函数可以围绕局部均衡分布正式展开，使得 （5.4） 式中 ——一个很小的参数，通常取Knudson数[38]； ——一个均衡分布函数； ——非均衡分布函数， 选择使其满足 且 （5.5） 并且假设非均衡分布函数可以进一步展开为 （5.6） 式中 k=1， （5.7） 碰撞算子由Bhatnagar-Gross-Krook近似给，即 （5.8） 式中 ——松弛时间. D2Q9格子在二维中分布函数的均衡形式由式 （5.9） 其中，当=1，2，3，4时，； 当=5，6，7，8时，； 松弛时间与运动粘度有关，且满足 （5.10） 格子Boltzamnn方法在几乎不可压缩极限条件下与Navier Stokes方程相同，并且在流体的内部满足二阶精确性.压强为p的不可压缩流体的应力张量由式 （5.11） 其中，为Kronecker函数且为应变率张量并满足 （5.12） 可证明在格子Boltzamnn方法中可以在每个节点的局部区域中计算即 : （5.13) 在格子Boltzamnn方法算法中，项通常作为速度计算的一部分计算.因此由于这种方法的计算剪切无需计算速度衍生，它是有效的.此外，在局部计算剪切时，如果并行使用格子Boltzamnn方法将会非常便捷. （2）幂法则模型 在随后的讨论中，我们定义应变率张量的第二个变量为 （5.14） 其中在二维情况下.剪切速率定义为 （5.15） 幂法则模型[39]是对非牛顿流体最简单的概括之一.在此模型中，表观粘度由 （5.16） 其中m和n通常是由拟合方程（5.16）得到的物理粘度参数数据.在刚性二维管中稳定流的情形下，这种模型有一个简单的解析解，如下： （5.17） 式中 L——管径； ——驱动流动的压力梯度. 参数n的值决定了在剪切速率变化时流体的响应，当<1时，流体为剪切变稀，当=1时，流体为牛顿流体，当>1时，流体为剪切增稠. 我们注意到，对剪切变稀流体（即<1）.同时m的量纲为，因此此参数与流体的任何物理性质无关.方程（5.16）可以无量纲地得到下列类似于Reynold数的无量纲数 （5.18） 其中m和为幂法则参数，为在管宽中的最大流量.
为了避免不可压缩，应确保马赫数小于0.03，此时模拟中的G，L和m可以取很多值.当满足如下的标准时，模拟才可进行，即满足 （5.19） 式中，很小，可取其为. 得到的结果与方程（5.17）给出的解析解相比较，计算出全局误差
其中下标a指精确解析解，下标b指Boltzamnn模拟值.全局误差结果见图5.1.2.
图5.1.2.幂法则模型流的全局误差
图5.1.3.当参数n分别取0.25（a），0.50（b），0.75（c），1.25（d）时，用格子Boltzamnn方法和解析解得到幂法则模型流的流速剖面的比较
图5.1.4.正常幂法则模型流的流速剖面.粗黑线代表标准牛顿流的抛物线分布 图5.1.2中的黑线代表了斜率为-2的线，表明了二阶性.可以看出，不同n参数代表的数据与线的斜率相符合.当n较小时，误差较大. 图5.1.3（a）-（d）描述了在（a）=0.25，（b）=0.50，（c）=0.75和（d）=1.25情形下标准化流速剖面分析结果（实线）与格子Boltzamnn方法所得结果（圈）的比较.可以看出，格子Boltzamnn方法精确地模拟出剪切变稀（<1）和剪切增稠（>1）流体的流速剖面. 当<1时，我们看到流速剖面整体比较平，且越小，流速剖面越平.相反，当=1.25时，在中间剖面的高峰值处有较大的曲率.在图5.1.4中可以更明显地看出这点，它描述了当不同时，正常流体剖面间的比较. 图5.1.2表明在二维幂法则流通过一个刚性管的情形下，格子Boltzamnn方法保持了二阶精确性.这种方法是对Gabbanelli等方法的改进，他们用一阶差分方法（Aharanov 和Rothman 1993）来估计剪切，只取得了一阶精确地结果.这里的结果提升到二阶的精确度，并且加快了算法的计算效率. 在图5.1.3中，对所有的，模拟得到的流速剖面与解析方法得到的剖面基本相符.当=0.25时，误差最大，此时格子Boltzamnn方法与解析方法相比，在中间高峰值处不够平.
5.2 格子Boltzamnn方法的方柱绕流模拟
在本章的算例中，对于边界条件，为叙述方便，左端为流场的入口，右端为流场的出口，并都取为充分发展的边界，流场上下边界均取为自由滑移的边界条件，方柱的边界均取为固壁无滑移的边界条件.
， （5.2.1） 式中 和 ——圆柱受到的阻力和升力.
图5.2.1单方柱绕流的流场设置 图5.2.2，图5.2.3分别为流线图和等涡线图，图5.2.4，图5.2.5分别为阻力系数和升力系数随无量纲时间T从第200步到第300步的变化的曲线.
图5.2.4阻力系数随时间变化 图5.2.5升力系数随时间变化 注：图5.2.2至5.2.5皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 从等涡线中我们可以看出漩涡分别从方柱后面的上下顶角处交替产生并在尾流中交替脱落，形成卡门涡街.从图5.2.4，5.2.5中可以看出阻力系数和升力系数在流场稳定后，它们的振幅近似为常数，阻力系数振动的频率约为升力系数的两倍，这些结果都与文献基本一致. Strouhal数是在方柱绕流中的一项重要的数据，其数值可由下式确定： （5.2.2） 式中 ——漩涡脱落的频率 我们通过纪录方柱后尾流8倍边长处一点的垂直方向速度随时间的变化，对其进行快速傅利叶变换得到.D为方柱边长，为流场的入口速度.我们计算所得的=0.1472，与文献[46]中的0.145相差1.5%.
在一个δ的渠道内放置长C=100（即C/L=10）的矩形截面柱，入口处速度是均匀的且保持不变，平板左侧与入口边界之间的距离为100.5，平板下端与渠道之间距离为90.5，雷诺数定义为，取Re=300情况下进行模拟，阻力系数和升力系数分别定义为： ， 图5.2.6给出在T=290时的流场图，所得结果与文献一致.图5.2.7，图5.2.8是平板的阻力系数和升力系数随时间变化的曲线.图中时间为无量纲的时间，结果与文献基本一致.
图5.2.7阻力系数随时间的变化 图5.2.8升力系数随时间的变化 注：图5.2.6至5.2.8皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 依照§4.1节中Strouhal数的定义方法，我们定义基于长度C的Strouhal数为，记录矩形截面柱后15倍L处一点的垂直方向速度随时间变化的情况，利用快速傅利叶变换，得出流场涡脱落的频率f=1.59×，图5.2.9给出当T=290时流场的等涡线图，图5.2.10为流场涡脱落的频谱图.
图5.2.9矩形平板周围的等涡线图
图5.2.10流场涡脱落频谱图 图5.2.11长宽不同比值下的Strouhal数 注：图5.2.9至5.2.11皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 将f代入公式（3.2）得到在Re=300，C/L=10情况下=1.59，与文献中的实验结果基本一致.应用如上方法计算矩形截面柱长宽不同的比值下的Strouhal数，并与2D Navier-Stokes方程的计算结果进行比较，我们发现在LBM方法下的值随着C /L的大小线性变化的.
计算区域如图5.2.12所示，为了更好地与单个方柱绕流比较，我们依然取计算长度，计算高度，方柱边长，方柱左侧边界距流场入口为，取雷诺数为100，流场及方柱周围的边界条件也与单个方柱的相同.
图5.2.12两并列方柱绕流的流场设置 我们计算了方柱之间间距L分别为0.2D，0.4D，0.8D，1.0D，1.2 D，1.4 D，1.6 D，1.8 D，2.0 D，3.0 D，4.0 D的情况.图5.2.13和图5.2.14仅给出了0.2D，0.8D，1.4D，1.6D，2.0D，4.0D时几种具有代表性情况的流线图和等涡线图：
图5.2.13不同间距下两并列方柱的流线图
图5.2.14不同间距下两并列方柱的等涡线图 注：图5.2.13至5.2.14皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 为了便于比较，我们还计算了迎风面长H=20，宽边长的单矩形方柱在相同边界条件下的流场情况，图5.2.15给出其流线图和等涡线图.
图5.2.15迎风面长H=20的方柱的流线图和等涡线图 注：图5.2.15皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 通过比较我们可以看出，当柱面间距较小时，流场更接近于迎风面H=20的单方柱流场状态，两方柱后形成共同的漩涡脱落，随着柱面的逐渐增大，，我们还发现，当间距L ≤1.4D时，方柱后的涡并没有完全分离，而是相互影响着脱落，流场比较复杂；而当L ≥1.6D时，两个涡的相互影响越来越小，逐渐呈现对称状态，当L达到一定值（L =4.0D）时，两方柱的涡完全分离，各自形成卡门涡街，涡的大小及脱落的频率接近于时的单方柱状态（见§5.2.1） 另外，我们通过观察不同间距时的流线图还发现，当L ≤1.4D，流场出现方柱绕流的奇特现象——偏流，这与实验中观察到的现象是一致的，随着L的逐渐增大，偏流现象逐渐消失，流场呈现对称状态. 下面分别对两方柱的阻力系数和升力系数进行比较，图5.2.16，图5.2.17，两方柱的阻力系数和升力系数在无量纲时间T从200步到300步的变化情况.
图5.2.16两方柱阻力系数随时间变化
注：图5.2.16至5.2.17皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 从图5.2.16、5.2.17中可以看出，当间距L ≤1.4D时，两方柱的阻力系数和升力系数并没有周期性变化，流场时而偏向方柱1，时而偏向方柱2，而当间距L ≥1.6D时，阻力系数和升力逐渐呈现周期状态.图5.2.18是特征时间T从200步到300步上、下方柱的时均阻力系数与时均升力系数随方柱间距变化曲线，可以看出随着方柱间距的逐渐变大，两方柱的阻力系数和升力系数都在逐渐变小.
图5.2.18上下方柱时均阻力系数与时均升力系数随方柱间距变化曲线 注：图5.2.18皆引用于《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英
我们考虑如图5.2.19所示的方柱绕流问题：在一个500100的渠道内放置一个边长a=10格子间距的正方形，正方形与流场上下两侧距离均为Yh=45.5 ，左端入口处来流为水平的剪切来流，速度梯度为G， 为无穷远处来流对称轴处的速度，渠道右端、上下边界及正方形的边界均取为和上文相同的边界条件.雷诺数定义为，粘性系数γ=（2-1）/6，我们确定一个较低的雷诺数Re=100，从而满足二维不可压缩流体的假设.
图5.2.19剪切来流情况下方柱绕流的流场设置 图5.2.20、图5.2.21给出当G分别等于0，0.001，0.002，0.003，0.0情况下的流线图和等涡线图.
图5.2.20 Re =100时，速度梯度G取不同值下流线图 注：图5.2.20《用格子Boltzamnn方法模拟棱柱绕流问题》，作者马金英 在均匀来流（G=0）情况下，在方柱后上下两侧形成两个大小、强度相同的涡，上方的涡顺时针转动，下方的逆时针转动，两个涡交错排列，呈现卡门涡街的特征.两个涡之间的距离也近似为常数.当G>0时，方柱后上方的涡比下方的提升得早，并且在方柱后面涡脱落的形式与均匀来流（G=0）情况不同.当G=0.001时，尾涡仍然呈周期性脱落，并在尾部依然呈现卡门涡街的特征，流线也和G=0时大致相同.随着涡的脱落，流线显示出一种周期性振动.从图5.2.20中可以看出当G=0.001时，方柱后上方和下方的涡的大小是不同的，上方的涡近似为圆形，而下方的为细长型；并且上方的涡在下游运动的速度比下方的快.这是因为我们所设的剪切流，上方的速度比下方的快.还有一点要说明的是：下方的涡的强度也随其移动越来越小，这是因为在下方形成的涡转动方向与其在流动中包含的方向相反，这样，下方的涡的耗散要比在均匀来流G=0的情况下快. 当G=0.002和G=0.003时，我们看到方柱后涡形成的区域增大，下方的涡的耗散越来越快. 流场最显著的改变出现在当G=0.004时，形成涡的区域显著增大，在方柱后形成一个细长的区域；下方逆时针旋转的涡完全消失，只剩下上方顺时针旋转的涡在方柱后形成一个单独的涡街. 5.3格子Boltzamnn方法模拟气固两相流
图5.3.1：后台阶模型草图 模拟的过程中，我们定义无量纲方程： 雷诺数，斯托克斯数，弗朗特数以及无量纲分数. 式中 ——入口处的速度； Μ——动态粘性系数； ——运动粘性系数； 和——单个粒子所对应的直径和密度； ——在体积空间中颗粒所占的体积； ——重力系数. 模拟的过程中，我们采用750×50网格，边界条件采用非平衡态外推方法来处理.
一直以来，大家都使用欧拉方法和拉格朗日方法这两种传统的方法来对带颗粒的流动进行数值模拟和研究.在文献中，Shirolkar曾经详细的对每种方法的基本特征进行了描述.特别是最近的几年，拉格朗日方法由于自身的优势，在工程实际中得到了越来越广泛的应用.从而在单相问题中使用这种运动方程可以很容易解决颗粒相的问题.在稀疏的带颗粒流动中，这个问题变得更加简单，这是因为为数众多的颗粒与颗粒的相互作用比如碰撞和合并都可以忽略不计.在本文的研究中我们假定颗粒的重量远远大于周围的流体，并且运动方程可以简化为： （5.3.1） 式中 和——对应颗粒和流体的速度； ——表征密度率，可以定义为； 系数——用来描述额外的斯托克斯拖拽力的影响：. 在Stk1的时候方程是一个刚性方程，因此我们利用一种预测-修正的指数拉格朗日轨迹跟踪方法来进行计算.可以表示为： （5.3.2） 式中 . 上标n——当前的时间步. 颗粒的速度和坐标通常可以这样表示：
（5.3.3） （2）流动相的计算 对于流动相，我们利用修正的D2G9模型来进行计算.根据前面所阐述的理论，我们引入质量力来描述颗粒相与流动相之间的相互作用，从而可以对D2G9模型进行修正.对于质量力，我们利用文献中阐述的方法来进行计算： （5.3.4） 式中 Vr——单个颗粒的体积与控制体体积的比例; ——控制体内所有的颗粒的数目; ——作用在单个颗粒上的拖拽力; F——线性分布在环绕在控制体周围的格点上. 修正后的演化方程可以表示为： （5.3.5） 式中 下标——速度分布（i=1～8）， 和t——格子间距和时间步长； ——颗粒流动的速度，. 无量纲的松弛时间τ则与动粘率有关，二者之间的关系可以用下式表示： （5.3.6） 本文的研究中，有一个稳定值0.51.则是相对应的临时平衡态分布函数，它可以用下式来定义： 参数满足如下的状态方程：，+2 =0.5；则可以用下式定义： （5.3.7） 流体的宏观压力和速度可以利用下面的式子来计算： （5.3.8） 最终，压力项可以这样计算： =0，i=1，…，8时，其他形式的也可以同理得出.同时，Navier-Stokes方程也可以用第一个方程推导得出： （5.3.9）
之间变化，斯托克斯数在1×—1×之间变化.以上的模拟都可以和前人的结论进行比较. ⑴雷诺数的影响 为了研究雷诺数的影响，我们首先模拟了后台阶模型上的单相流动.这种流动的行为特性可以简要概括如下：流体在台阶处分离，然后在下游的位置x1重新汇合起来，坐标x1几乎是呈线形增长.当雷诺数比较大的时候，将会出现一个足够大的反方向的压力梯度从而在上面出现一个再循环区域，这个也就是俗称的涡，这一现象在图5.3.1中有特别指出.上面的这个涡随着雷诺数的增大也不断的扩大，但是它的中心逐渐的下移.对于这种流动的模拟，一般是以Armaly的实验数据 以及Guj和Stella的计算结果作为比较的根据.在图5.3.2中，我们把模拟结果和实验数据进行了比较.图5.3.3中，我们比较了模拟结果与Guj和Stella的计算结果.从这两幅图可以看出，雷诺数比较低的时候，修正的D2G9模型能够得到令人满意的模拟结果，当雷诺数增大时，模拟结果出现了比较明显的误差.我们参考了文献，觉得有可能是三维的影响导致了这种偏差.
图5.3.2：与实验数据的比较
图5.3.3：与计算结果的比较 注:图5.3.2与图5.3.3皆引用于《格子Boltzamnn方法基本模型比较及应用》，作者彭伟 ⑵斯托克斯数的影响 在文献中，Barton指出，对于低斯托克斯数的颗粒具有跟随流体的趋势，并且这种粒子可以加固流体中的自由蒸汽，从而增大下面的涡；而对于高斯托克斯数的颗粒，由于受到重力的影响将会增大向下的运动趋势，从而会压缩下面的涡使得上面的涡扩大.对应于以上情况，我们在雷诺数Re=450，空隙组分的情况下，变换不同的斯托克斯数来变动工况，进行了一系列的模拟工作，来讨论斯托克斯数对于汇合和重聚点的影响，并把结果绘制在图5.3.4中.为了便于比较，我们同时绘制出单相流体的重聚和分离长度（破折线表示）.从图5.3.4中我们可以清楚看到，随着斯托克斯数的变化，带颗粒流动的x1和x2也将发生变化，不过相对单相流体却始终保持同大或者同小的性质.同时，带颗粒流动的x3将始终大于单相流体的这个长度.这些性质在文献中都能够找到近似的结论.
图5.3.4：斯托克斯数对汇合和重聚位置的影响图 注:图5.3.4引用于《格子Boltzamnn方法基本模型比较及应用》，作者彭伟 与其他模型的比较 在最近的一篇文献中，作者也利用郭照立等人提出的处理外力项的方法，修正了He-Luo模型，并且模拟了后台阶流动，也得到了令人满意的结果.有鉴于此，我们把这两个模拟的结果进行了一些比较，试图找出一个在处理外力项方面比较优越的格子Boltzamnn模型.在图5.3.5中，我们绘制出了这两种模拟的误差分析图.
图5.3.5：不同模型间的比较 注:图5.3.5引用于《格子Boltzamnn方法基本模型比较及应用》，作者彭伟 从图中我们可以看出，分别修正He-Luo模型和D2G9模型来模拟气固两相流所得到的结果非常接近，最大误差也不超过5%.从而我们认为，郭照立等人提出的处理外力项的方法也具有比较好的通用性，利用它来修正不同的格子Boltzamnn模型进行模拟的结果也比较吻合.
5.4 格子Boltzamnn方法模拟油水两相流软件设计
边界条件的处理方法在格子Boltzamnn方法中起重要的作用，对格子Boltzamnn模型的精度和稳定性都有很大的影响，因而必须对边界的选取给予足够的重视.关于格子Boltzamnn方法中边界条件的处理，学者提出过不少方法，如Zou的非平衡态反弹格式，Noble的动力学格式及其改进格式，但各种方法都有其不足，如Inamuro的反滑移速度格式通用性不好，Chen提出的外推方法的通用性较好但数值稳定性不好.受Chen的外推方法和Zou的非平衡态反弹格式方法的启发，郭照立等提出了一种新的边界处理方法—非平衡外推方法.它的数值稳定性好于根据线性外推得到的分布函数外推格式的数值稳定性，同时与分布函数外推格式一样，非平衡态外推格式也具有适用范围广，计算简单，容易实现的优点. 本文采用了郭照立的非平衡外推方法，下面以D2Q9模型为例说明非平衡态外推方法的基本原理.如图5.4.1所示，假设EOA位于边界上，BCD位于流场内，FGH位于流场外.设系统的当前时间为t.此时，流体格点C处的分布函数值都是知道的.为确定边界点O处的分布函数，将流体的分布函数分解为两部分： （5.4.1） 式中 ——非平衡态部分. 根据速度边界条件确定出平衡态部分.对速度边界条件，已知而密度未知，所以可以使用下面的虚拟平衡态分布函数来近似原来的平衡态部分： （5.4.2） 由公式5.4.1和公式5.4.2，可以推出非平衡态外推格式边界条件： （5.4.3）
图5.4.1 D2Q9速度矢量图 注:图5.4.1引用于《油水两相流流型检测方法与研究》，作者韩连福 ⑵多相流LBM模型的选择 由于格子Boltzamnn方法的微观粒子特性，它可以方便地描述不同相间的相互作用，使之在模拟多相流问题上具有常规方法所没有的优势.从而提供了研究多相流的一种有效途径.按照设计方法模拟多相流的格子Boltzamnn模型有着色模型、伪势模型、自由能模型和其他模型.Gunstensen等在1991年首先推导出LBM方法的多元流模型，Grunau等对这种模型进行了推广提出了着色模型，这个模型成功的模拟了二维相变问题，但是这个模型有一个缺点是在界面附近会产生非物理现象，并且“颜色能量"极小化过程的计算量比较大.Swift等人通过修改平衡态分布函数，提出了基于自由能函数的格子Boltzamnn模型.自由能模型解决了着色模型的缺点，也考虑了热力效应，但是不满足伽利略不变性.另一类模型是伪势模型，它可以自动追综界面的运动，相的分离过程是自动的，计算效率大大提高，可以用于复杂的多相流的系统，因而本课题采用了伪势模型，下面简单介绍伪势模型. 优势模型基于场论理论，即通过一个伪势函数反映不同相间的斥力和引力，从而导出非理想状态方程.在此方法中，每一相采用不同的分布函数表示，分布函数的演化方程为： （5.4.4） 相间的作用力为 （5.4.5） 式中 ——源于某个相互作用的伪势函数 （5.4.6） 伪势函数V与密度相关，可以为指数形式，也可以为分数形式，为方便计算本课题采用了指数形式，其伪势函数可取为： （5.4.7） 式中 G——控制相互作用的强度， 当其大于0为斥力，当其小于零为引力； ——参考密度； 本研究中取为1. ⑶模型格子的选取 实验管径为125 mm，模拟井高度为12 m.为模拟方便，在不影响结果的情况下管径为125 mm，长度为250 mm，大致如图5.4.2所示，黑点表示格点，黑线表示格点连线，我们取M=124，N=249.
图5.4.2格子选择示意图 在处理平衡态部分时，C取为离O点最近的格点， 如（i，0）（表示第i行，第0列），我们取相应的C为（i，1）. 当地声速c的计算，通过以下三式： （5.4.8） 这里L，U都无量纲化为1 （5.4.9）粘性系数v定义为 : （5.4.10） 式中 ——松弛时间， ω——松弛系数. 我们可以得到D2Q9模型的c值计算式： （5.4.11） 演化的实现 我们把连续的粒子运动离散成网格格点上的运动，运用C语言编程如： （5.4.12） 写成C程序为： id=i-e[k][0]; jd=j-e[k][1]; F[j][i][k]=f[jd][id][k]-1/τ（f[jd][id][k]-feq[jd][id][k]）;
在软件工程中，需求分析指的是在建立一个新的或改变一个现存的电脑系统时描写新系统的目的、范围和定义时所要做的所有的工作.需求分析是软件工程中的一个关键过程.在这个过程中，系统分析员和软件工程师确定顾客的需要.只有在确定了这些需要后他们才能够分析和寻求新系统的解决方法.在软件工程的历史中，很长时间里人们一直认为需求分析是整个软件工程中最简单的一个步骤，但在过去十年中越来越多的人认识到它是整个过程中最关键的一个过程.假如在需求分析时分析者们未能正确地认识到顾客的需要的话，那么最后的软件实际上不可能达到顾客的需要，或者软件无法在规定的时间里完工. 本软件的用户是实验室的学生，他们具有一定的油水两相流的基础知识，综合他们的意见定出如下功能需求和设计约束. 功能需求： （1）要能够实现数据的录入和输出功能. （2）要具有数据的存储功能，存储的文件格式为tecplot软件可以直接识别的格式. （3）要具有数据的可视化功能. （4）油井的倾斜角度为在0-90°内任意可调. 设计约束： （1）所有的单位采用油田生产中的常用工程单位，如流量用方/天表示. （2）油井的倾斜角度以水平为基准. ⑵软件的流程和界面的设计 本文采用了Windows XP下的C++6.0开发油水两相流模拟软件.应用其下的MFC AppWizad创建项目，然后进行编写代码.经软件需求分析可以知：软件系统应包含数据输入模块、数据保存模块、数据处理模块、数据可视化显示模块和LBM模拟运行模块.本次实验水流量范围为：5方/天-180方/天，步进取为0.5方/天；油流量范围为：5方/天-360方/天，步进取为2方/天.系统的程序流图如图5.4.3和图5.4.4所示.
图5.4.3系统程序流图 根据系统的需求本软件界面设置了三个区域分别为：油参数设置区，水参数设置区和实验环境记录区.油参数设置区用来输入油的密度、流速、流速步进、粘性系数和最终流速5个参数；水参数设置区用来输入水的密度、流速、流速步进、粘性系数和最终流速5个参数；环境记录区用来设置环境温度、模拟时间和油井的倾斜角度等信息，油井默认为垂直井.如果只需分析某个特定油流速和水流速下的流型，步进选项和最终流速选项可以不输入.下面介绍一下本软件的设计步骤： (1) 从C++编译器的菜单条中，选择File菜单，从列出的菜单项中选择New选项， 在出现的对话框中选择MFC AppWizad选项，输入项目名youshui选择存储位置，这样就开始建立一个项目. (2)在向导中，按照向导提示逐级选择，最后点击结束按钮就建立了一个项目. (3)在主对话框中如图5.4.5放入控件，这样一个主页面就设计完了. (4)点击各个控件创建各个控件的事件代码. (5)进行软件的编译调试工作.
图5.4.4 LBM算法程序流程图 软件模拟可视化结果 图5.4.5为水包油流型模拟图，在水包油流型中红色区域代表高密度的水，蓝色区域代表低密度的水.从图中初步分析可知：在这种流型下水包着油，并且在管径不同深度位置的油水的分布也不一样. 图5.4.6为过渡流型模拟图，在过渡流型中红色区域代表高密度的水，蓝色区域代表低密度的水.从图中初步分析可知，在这种流型下油水的分布没有什么规律，是随机的. 图5.4.7为水包油流型模拟图，在油包水流型中红色区域代表高密度的水，蓝色区域代表低密度的水.从图中初步分析可知：在这种流型下油包着水，并且在管径不同深度位置的油水的分布也不一样.
图5.4.5水包油流型模拟图
图5.4.7油包水流型模拟图 注:图5.4.5,5.4.6和5.4.7引用于《油水两相流流型检测方法与研究》，作者韩连福
式中 u——平移速度； ——旋转速度； R——颗粒的中心位置； ——边界节点. 为节点上的分布函数，而颗粒边界上的分布函数为：
式中 B——系数. 根据网格模型和质量守恒而定，当颗粒迎向流体运动时，用“+"号；当颗粒背离流体运动时，用“-"号.他模拟了有限Re下的蠕动流，其结果FDM和FEM的数值解相比是非常一致的，证明了其正确性.Aidun等在Ladd方程基础上，在上增加了源项，使得颗粒内部质量也守恒，增强了Ladd模型的理论，有效地减少了边界误差且提高了数值精度.在LBM方法中分布函数的脉动被忽略了，但在布朗运动中脉动的研究却非常重要.Dufty等在研究布朗运动时在分布函数中增加了随机项，包含了流体运动的脉动，与实验观察到的现象相吻合.
式中 =0，1， . 即局部能量守恒，Vahala等通过研究二维剪切紊流温度梯度的变化对这种模型进行了验证.Qian推导了D3 Q21和D3 Q25的LBM方法热力学（TLBM）模型，McNamara等 分析了这种TLBM模型的两大局限之处： （1）由于速度集合中的向量很少，所以温度的变化比较小； （2）由于缺少H—定理，导致这种模型存在不稳定性.Bartoloni等提出将温度作为标量，另外增加一个独立的分布函数处理温度就可部分地消除上述两个局限，然后Shan将温度作为标量，成功地求解了2D和3D的Rayleigh-Benard热传导问题，与理论分析相吻合.在以前的研究中，热传导问题和多相流问题往往是分开进行的，而LBM方法热力学模型的进一步发展将使两种问题的求解过程融为一体.这种模型的研究将开拓出新的应用领域，弥补传统计算方法的缺陷.
其中k=0代表红色液体，k=1代表蓝色液体.碰撞项由两部分组成：
总体密度其中 .用这种模型模拟了Hele-Shaw粘性指甲流的形成过程，取得了良好的模拟效果. Shan和Chen等利用宏观的相互作用，修改了与表面张力相关的碰撞操作，这种Shan-Chen模型一个重要的进步就是液体的相分离或组成的分离是自动的，Hou等通过静态气泡的演进过程，分析了表面张力、Laplace法则和交界面的厚度，认为这种Shan-Chen模型对多相流的模拟已经达到十分成熟的水平.Shan&Chen和Swift等研究了毛细波的振荡问题，所得到的离散关系与理论解是一致的.此外，Swift用自由能量法模拟了液-气两相流系统，Alexander等和Rybka等模拟了两种液体的分离过程，在复杂的交界面中充分利用LBM方法的计算特性，得到了满意的结果.He等模拟了三维Rayleigh-Taylor不稳定问题，形象而生动地描述了两种液体的交融过程.他们得到了二维模拟中不能发现的鞍面点，并且比重大的液体在下降形成尖峰，在鞍面点反向上升，形成对称的锚状；比重小的液体上升形成玫瑰花形，产生花瓣向四周扩散现象，每个断面都成规则的对称性，直到两种不同比重的流体交融.LBM方法的多相流模型还可以模拟多孔介质流的渗透系数、粘性系数和湿润度等参数，显示出诱人的应用前景. 从以上的介绍中可以看出，格子Boltzamnn方法是一种不同于传统的数值方法则很难做到这一点.另一方面，从计算的角度看，格子Boltzamnn方法属于显式时间推进方法，每个时间步的计算量O（MN）（M为离散速度数，N为计算格点数），其计算效率要高于一般的数值方法.同时，格子Boltzamnn方法（尤其是LBGK模型）的演化过程非常简单清晰，程序比较简洁；格子Boltzamnn方法中涉及的计算都是局部性的，具有天然的并行性，非常适合在大规模并行计算机上运行.正是由于具有这些优势，格子Boltzamnn方法自诞生之日起就受到包括物理学家、数学家、计算机专家和其他领域的科学家的广泛关注.LBM正处于不断的发展之中，近年来在基本理论、基本模型和应用等各方面都有所发展.
结论及展望
（a）边界条件易于处理. （b）内节点的作用原则完全是一样的，程序的实现是很简单的，并且很适合于大规模的并行运算. （c）格子气的基本原则是物理守恒定律，他的基本方法和思想一样适用于研究其他物理现象. （d）可以非常方便的实现可视化技术. 此方法不但适合于流体力学的各个方面，并且逐渐渗透到物理学的其他方向，比如相变以及晶体增长等，对于研究非线性复杂系统领域有着巨大的潜力. 格子Boltzamnn方法是最近20年来发展起来的，一种新型的介观数值模拟方法.与传统方法相比，格子Boltzamnn方法具有多种优点，而因受到国内外学者的广泛关注、并迅速在多孔介质流、多相多组分流、非牛顿流体、粒子悬浮硫、湍流、化学反应流、燃烧问题、磁流体、微/纳米尺度流等许多领域得到应用. 经过国内外众多学者的努力、格子Boltzamnn方法的理论正逐步发展，但距离最终的学科体系完善还有一个相当长的过程.目前，格子Boltzamnn方法仅用于一些较简单的物理模型.对计算流体力学领域的其他复杂问题，采用格子Boltzamnn方法模拟还有较大的困难.此外，格子Boltzamnn方法的工业应用也鲜见报道，如何将其应用于复杂系统以及更广泛的工业实际应用中，是今后可能的发展方向. 参考文献
[2] 苏铭德，黄素逸，计算流体力学基础[M]，清华大学出版社，1997． [3] 郭照立，郑楚光，李青，王能超.流体动力学的格子Boltzamnn方法[J]，湖北科学技术出版社，2002. [4] 彭伟.格子Boltzmann方法基本模型比较及应用[D].华中科技大学硕士学位论文. 2004:14—42；50—57. [5] 周光垌，严宗毅.流体力学（下册）[M]，高等教育出版社，1993． [6] 吴大猷.气体运动论及统计力学[M]，科学出版社，1983． [7] LuoY, eden JM．Flow of Non-Newtonian Fluid Through Eccentric Annul[SPE]. 16692,1987 [8] 李元.格子Boltzamnn方法的应用研究[D]. 中国科学技术大学博士学文.2009:26—37. [9] 许福，确定细胞自动机的理论和应用[D].武汉大学硕士论文.1993：22—25. [10] 康立山，全惠云，等.数值求解高维偏微分方程的分裂法[M].上海科学技术出 版社.1990:122—151. [11] 张武生,杨燕华,徐济. 格子波尔兹曼方法及其应用[J]. 现代机械,2003(3). [12] 祖迎庆. 以Boltzamnn方程为基础的流体力学数值方法[D]. 吉林大学硕士学位论文. 2005:27—46. [13] 陈明杰, 施卫平.格子Boltzamnn方法计算剪切流的方柱绕流问题[J]. 吉林大学学报. 2007(1). [14] 郭照立，郑楚光，等.格子Boltzamnn方法原理及其应用[M].北京：科学出版社，2008. [15] 何雅玲，王勇，等.格子boltzamnn方法原理及其应用[M].北京：科学出版社，2008. [16] Michael C. Sukop Daniel T. Thorne,Jr. Lattice Boltzmann Modeling[M]. Florida International University Department of Earth Sciences University .Miami,2006. 致 谢
同时也要感谢数学系的各位老师，各们老师在开题、中期检查的过程中提出了很多贵的意见，在此，向所有给予我帮助的老师表示衷心的感谢.
|
|
| ■■相 关 文 章： | 网友评论：（只显示最新10条。评论内容只代表网友观点，与本站立场无关！） |
| 没有相关毕业 |
|
||||||||||