毕业设计（论文）-格子Boltzmann方法及其在非牛顿流体力学的应用.docx

格子Boltzmann方法及其在非牛顿流体力学的应用大连理工大学本科毕业设计(论文)格子Boltzmann方法及其在非牛顿流体力学的应用The lattice Boltzmann method and its application in non-Newtonian fluid dynamics学院(系)：能源与动力学院专业：制冷与低温工程学生姓名：学号： 201172169 指导教师：评阅教师：完成日期： 2015年7月12日大连理工大学Dalian University of Technology摘要作为一种计算流体动力学技术，格子Boltzmann方法(LBM)已经成为研究气固两相流的有力工具。与传统方法相比，该方法有很多优点：程序简单、计算高效、天然并行性使其成为处理复杂流动现象的有力工具。格子Boltzmann方法能够直接计算局部剪切速率并可以达到二次精度，因此在非牛顿流动数值模拟中展现出一定优势。同时，非牛顿流体广泛地存在于自然界中，90%以上的流体都具有非牛顿特性，非牛顿流体在工业生产和日常生活中占有重要地位，因此，深入研究非牛顿流体有很大的实用价值。本文利用格子Boltzmann方法模拟了幂律流体(power law fluid)在直通道中的流动，并重点研究了宾汉姆流体(Bingham fluid)流经突扩通道的流动特性，应用了修正的宾汉姆模型(Papanastasiou模型)，通过增加一个足够大的应力增强指数来克服宾汉姆流体模型的不连续性，根据修正的宾汉姆模型模拟宾汉姆流体的流动特性。主要工作包括：研究LBGK模型及其边界处理和作用力模型以及非牛顿流体的格子Boltzmann模型，然后用C+编程对直通道中幂律流体的流动进行模拟，并且与解析解进行了比较。结果显示数值解与解析解吻合良好。随后对宾汉姆流体在2:1突扩通道中的流动过程进行了模拟，固定宾汉姆数(Bn)，选取不同的雷诺数(Re)数来考察屈服区域和非屈服区域的分布情况以及流场情况，最后通过无量纲漩涡长度定量分析了模拟结果的正确性。关键词: 格子Boltzmann方法; 非牛顿流体; 二维直通道; 突扩通道. The lattice Boltzmann method and its application in non-Newtonian fluid dynamicsAbstractThe lattice Boltzmann method (LBM),as an novel technology of the Computational Fluid Dynamics (CFD),has become a useful tool for studying gas-solid two-phase flow. Compared with traditional methods, this method has a lot of advantages, such as the simplicity of program, high efficiency of computation, natural parallelism which is a strong tool to simulate the complex flows. LBM can be directly used in calculations of the local shear rate to the second-order accuracy, so it has some advantages in the simulation of the non-Newtonian behavior of flowsMeanwhile, non-Newtonian fluid widely exists in nature. More than 90% of the fluid has non Newtonian characteristics. Since non-Newtonian fluids have played an important role in industry and our daily life, further study of non-Newtonian fluid has great practical value.In this paper, the flow of power law fluid in a straight channel and the flow characteristics of Bingham fluid, flow through sudden-expansion channel is simulated via the lattice Boltzmann method. In this study, the modified Bingham model (Papanastasiou model) is employed, in which a stress enhancement index is introduced to overcome the discontinuity of Bingham fluid model, therefore flow characteristics of Bingham fluid is simulated based on modified Bingham model. The main work includes the research on LBGK model and its boundary processing and force model,and the study on the Boltzmann model of non-Newtonian. The code is written on the basis of the C+ language programming to study the flow characteristics of power law fluid in a straight channel. Comparison of numerical solution with analytical solution verifies the accuracy of the simulation. Then the flow of Bingham fluid, flow through a 2:1 sudden-expansion channel is simulated. The distribution of yield area and non-yield area, and the flow field are studied, under the conditions of fixed Bingham number and different Reynolds numbers. Finally, the dimensionless vortex length is used to quantitative analysis the accuracy of simulation results.Key Words: Lattice Boltzmann method; Non-Newtonian fluid; Two-dimensional straight channel; Sudden-expansion channel.目录摘要IAbstractII目录III1 文献综述11.1 引言11.2 格子Boltzmann方法介绍31.2.1 格子Boltzmann方法的发展历史和研究现状31.2.2 格子Boltzmann方法的基本结构41.3 本文的主要工作52 LBGK模型及边界处理72.1 LBGK模型介绍72.1.1 格子Boltzmann-LBGK方程72.1.2 DnQb离散模型82.2 初始条件及边界条件处理112.2.1 初始条件112.2.2 边界条件处理122.3 格子Boltzmann方法的几种作用力模型152.3.1 平衡态分布的压力校正方法152.3.2 平衡态分布的速度校正方法162.3.2 在演化方程中增加作用力项162.4 格子Boltzmann方法的计算步骤173 非牛顿流体的格子Boltzmann模型193.1 非牛顿流体概述193.1.1 非牛顿流体的基本特征193.1.2 非牛顿流体的粘性193.1.3 非牛顿流体的弹塑性193.2 幂律流体与宾汉姆流体本构方程203.3 非牛顿流体LBM模型的计算过程224 二维直通道内幂律流体流动的LBM模拟234.1 模型描述234.2 数据分析234.2.1 流速分析244.2.2 剪切速率分析264.2.3 粘度分析285 Bingham流体在2:1突扩通道内的LBM模拟305.1 程序验证305.2 问题描述315.3 数据分析325.3.1 剪切应力分布325.3.2 流线分布346 总结36参考文献37致谢38IV1 文献综述1.1 引言自然界中存在着大量非牛顿流体，例如油脂、油漆、牛奶、牙膏、动物血液和泥浆等。非牛顿流体力学在化学纤维工业、塑料工业、石油工业、化学工业、轻工业和食品工业等许多部门有广泛的应用1。由此可见，我们无时无刻都离不开非牛顿流体，准确预测非牛顿流体的流动行为有着不可估量的工程应用价值，所以对非牛顿流体的研究具有重大意义。非牛顿流体力学的研究对象主要是流体，它要研究的是流体的流动与变形。因此，非牛顿流体力学就是研究流体流变学的科学，也可称为流体流变学。按照剪切应力与变形率之间的关系，可将流体分为牛顿流体和非牛顿流体2。对于牛顿流体，剪切应力与流速梯度dudy成线性关系，即=dudy，其中的是在任意给定温度、压强条件下牛顿流体流动的特征性比例系数，该比例常数称为流体粘度。而非牛流体的物质构成远比牛顿流体复杂，反映在粘性上，同一种非牛顿流体的粘性不仅与温度、压力等环境条件有关，而且还与其流动状态有关，一般很难用单一的系数来表征其实际的粘性3。依照本构方程的不同可以将非牛顿流体做以下分类：广义牛顿流体、有时效的非牛顿流体和粘弹性流体。非牛顿流体中的与应力历史无关的流体就称为广义牛顿流体。我们生活周围的大多数非牛顿流体均为广义牛顿流体。人们根据粘性的不同特性将广义牛顿流体分为以下几种情况4：(1)假塑性流体(Pseudoplastic)：粘性系数随着剪切率的增大而减小，也称剪切变稀流，包含大多数高分子聚合物的溶液或溶质。(2)膨胀流体(Dilatant)：黏性系数随着剪切率的增大而增大，也称剪切增稠流，玉米面糊就属于膨胀流体。但这类流体在我们生活中并不常见。(3)宾汉姆流体(Bingham Plastics)：当切应力低于某一屈服力临界值时，流体呈固状；超过这屈服力后，流体由静止表现出流动现象。常见的塑性体有：烂泥、油墨和血液等。其中假塑性体和膨胀流体归属于幂律型流体，所以广义牛顿流体包括两大类型：幂律型流体和宾汉姆流体。管道突扩流动在化工、生物工程、医药、食品、物料混合、热交换、燃烧和核反应等工程实际中被广泛应用。它作为流体力学的一类基础问题，一直受到人们的重视。管道突扩流动在几何上看似简单，但由于其流与管道内受限射流的相互作用，使得问题变得非常复杂。在气-固或液-固两相突扩湍流中，颗粒相的加入会影响流体相的流场结构和湍动特性等，使问题变得更为复杂。管道运输中，突扩管道中流体的流动是一种常见的实际状况。由于流线自身所拥有的特殊性质，导致流体流经突扩管道时将会产生流动的分离和再附现象，这是学者们一直在关注的问题。因为产生的分离和再附会引发涡旋，从而导致压力的降低和能量的损失。本文将会对突扩通道内的宾汉姆流体流动进行模拟计算。图1.1 突扩管的物理模型从方法论的角度讲，对于流体的流动可以分别从宏观、介观和微观三个层次进行。在宏观层次上，可以将流体假设为连续的介质。流体的运动满足质量守恒、动量守恒和能量守恒，并由Euler方程组、Navier-Stokes方程组等描述。现有的大多数场模拟方法都属于宏观方法，如有限差分法(finite-difference method)、有限元法(finite-element method)、有限容积法(finite-volume)、有限分析法(finite-analytic method)、边界元法(boundary element method)和谱方法(spectral method)等。这些方法的发展已经较为成熟，既可以进行物理问题的机理研究，还可以发现一些新的物理现象。一些基于宏观方法的通用大型商业软件也应运而生，如FLUENT、STAR-CD、CFX和PHOENICS等，也被应用于解决工程实际问题。在微观层次上，流体是由大量的离散分子组成，分子运动特性是由分子间相互作用力以及外加作用力的影响。因而，一个很直接的想法就是通过模拟每一个分子的运动，再基于不同的法则进行统计平均，以获得流体流动的宏观规律。1957年，Alerder和Wainwright首先通过硬球模型，采用分子动力学研究液体的状态方程5，从而开创了分子动力学模拟(molecular dynamics simulation)的先例。在介观层次上，流体被离散成一系列的流体粒子。这些粒子的级别比分子要大。但从宏观上看，又无限小，其质量比有限容积法中的控制容积质量要小得多。单个分子的运动细节并不影响流体的宏观特性，通过构造符合一定物理规律的演化机制，让流体粒子进行演化计算，并获得与物理规律相符合的数值结果。格子Boltzmann方法6源于格子气自动机，将流体离散成流体粒子，并且物理区域也被离散成一系列的格子，时间也被离散成一系列的时步。格子Boltzmann方法具有物理意义清晰，程序易于实施，边界条件处理简单，模型健壮性高，并行性能好等特点。所以，它被认为是最具有前途的数值模拟方法之一7。下面将对格子Boltzmann方法(LBM方法)进行详细的介绍。1.2 格子Boltzmann方法介绍1.2.1 格子Boltzmann方法的发展历史和研究现状格子Boltzmann方法(lattice Boltzmann method，LBM)自诞生以来，在理论、模型和应用等各方面都得到了迅速的发展，业已成为模拟各种流动现象的有力工具8。与有限差分法或有限元法等传统数值模拟方法不同，格子Boltzmann方法属于介观模拟方法。流体不再被假设为连续介质，而是抽象为离散的粒子，通过数值求解分布函数获得宏观的流动信息9。由于LBM具有演化过程清晰，并行性好，适合大规模数值计算等优势，迅速在微纳米流体、多孔介质和非牛顿流体等领域得到广泛应用10。20 世纪 50 年代开始，元胞自动机的概念开始萌芽，并有学者将其应用于流体力学的研究中。不同于一般的动力学模型，元胞自动机不是由严格定义的物理方程或函数确定，而是用一系列模型构造的规则构成，是一种时间、空间和状态均离散的动力系统。随后不久，Broadwell等学者提出了离散速度模型，用离散化思想将流动速度空间朝着某些特定方向进行分解，但从宏观上来说，时间和空间的连续性依然存在。20 世纪 70 年代，为了研究流体的输运特性，Hardy、Pomeau和Pazzis提出了第一个完全离散模型，即HPP 模型。该模型同时考虑了流体、时间和空间的离散化，并将其离散到二维正方形格子。但是由于正方形格子的非完全对称性，使得模型中的应力张量无法满足各项同性，对应的宏观方程不能够表示准确的非线性以及耗散效应。1986 年 Frisch、Hasslacher和Pomeau提出了具有足够对称性的二维正六边形格子气自动机模型，即FHP 模型。同年dHUmieres、Lallemand和Frisch提出了四维面心立方模型，即FCHC模型。这两种模型都成功克服了HPP模型对称性不足的缺点，在一定条件下可以得到二维和三维的不可压缩Navier-Stokes 方程。自格子气元胞自动机模型(Lattice Gas Automata，LGA)出现以来，得到了一定的发展和进步。LGA模型将流体视为驻留在一个规则格子上的大量假想粒子，这些例子按照一定的规格在格子上进行碰撞和迁移，通过统计学的知识得到宏观平均值。虽然LAG方法物理清晰直观，计算较为稳定，但是仍然存在一些问题需要解决。例如结果中会含有统计噪声，需要对时间和空间进行数值平均，从而在很大程度上增加了计算量。同时推导出的动量方程不满足Galiean恒定性，并且算子具有指数复杂性，对存储空间也有较为严格的要求。为了解决这些问题，LAG方法也在不断被修正和改进。1988年，McNamara 和 Zanetti提出在LGA中使用布尔变量的统计平均量粒子分布函数进行演化以消除统计噪声。但是，LGA方法的指数复杂性仍未得到解决。随后为了进一步对LGA方法进行优化，Higuere和Jinenez提出了线性化碰撞算子模型，近似代替计算量很大的指数型算子。在这种模型中，引入分布函数距离其平衡态分布函数，即。其中是Fermi-Dirac分布的展开形式，为非平衡部分。在低速条件下，可以进一步将按照速度项进行展开：(1.1)其中为包含速度uk的项。将碰撞算子在全距平衡态处展开， (1.2)略去高阶小量可得简便形式如下： (1.3)取，同时略去高阶小量，可得 (1.4)于是得到了一个近似线性化的碰撞算子，其中K为线性化碰撞矩阵。由于采用了线性化手段，HJ模型大大提高了计算效率，为LBE模型的发展奠定了重要基础。1991-1992年，一些专家学者又相继提出了单松弛模型(SRT模型)和BGK模型。在格子BGK模型(即LBGK模型)中，碰撞过程用区域某一平衡态的松弛过程代替，而碰撞矩阵由松弛时间决定。这一模型完全克服了格子气自动机的一系列缺点，因此也是目前最基本且应用最广泛的计算模型11。1.2.2 格子Boltzmann方法的基本结构格子Boltzmann方程模型包含三个要素：(1)流体粒子的离散速度集合，(2)格子结构，(3)演化方程。它描述了具有离散速度的流体粒子分布函数在一个固定格子上的演化运动过程。 (1.5)其中x是格子上的一个格点，ci (i=1，2，3，4，b)是流体粒子的离散速度集合，t为离散时间步长，t为当前时间步；fi是以速度ci运动的速度分布函数，i是碰撞算子，表示分子间的碰撞对速度分布函数的影响。流体的密度、速度u和内能e可由离散分布函数的速度矩得到。(1.6)其中D为空间维数。演化方程又称为格子Boltzmann方程，其中的碰撞算子反映了微观流体的相互作用。因此碰撞算子对模型能否准确刻画出流体系统的物理规律起到了决定性的作用。目前LBE方法一般采用线性化碰撞算子，即 (1.7)其中K是一个bb的碰撞矩阵，是依赖于宏观物理量的平衡态分布函数。因此，碰撞矩阵和平衡态分布函数完全决定了碰撞算子的特性，也就相应决定了一个LBE模型刻画流动问题的精度。从上述分析中可以发现，格子Boltzmann方法是一种不同于传统数值方法的流体计算和建模方法。从离散的网格来说，LBE具有Euler方法的属性，而从离散的粒子观点来说，LBE又具有Lagrange方法的特点。这种特性使得LBE方法具有以下几个显著特点：(1)虽然LBE模型主要用于模拟宏观连续流动，但它不是基于连续模型的N-S方程，而是基于介观模型，本身没有连续介质条件的假设。因此，只要设计的参数合理，可以用于描述微观尺度和稀薄流等非连续流动问题。(2)格子Boltzmann方法可以较为直观地处理流体内部以及流体和周围环境的相互作用，从而给多组分、多相态系统、界面动力学、渗流等复杂流动现象的研究提供了方便；(3)从计算的角度来看，LBE方法计算简单，易于用编程实现，具有较好的并行性和可扩展性，对于大规模流动问题的计算具有很大优势。1.3 本文的主要工作首先是通过对LBM理论方法的学习，掌握LBM方法的整个计算流程，以及整个方法计算过程中各部分之间的联系，包括粒子的碰撞迁移过程，和如何选取边界条件，不同的边界格式的适用条件对于程序计算过程稳定性的影响。其次是非牛顿流体LBM模型的构建，包括幂律流体和Bingham流体。在LBGK模型的基础上，结合非牛顿流体的本构方程，建立了本文的非牛顿流体计算模型。由于Bingham流体本构方程本身的不连续性，对模拟造成了一定的困难，所以模拟中选取了Papanastasiou模型，通过模拟直通道中幂律流体的流动，并且与解析解进行了比较，数值解与解析解吻合良好。最后对Bingham流体在2:1突扩通道中的流动过程进行了模拟，固定Bn数，选取不同的Re数来考察屈服区域和非屈服区域的分布情况以及流场情况，最后通过无量纲漩涡长度定量分析了模拟结果的正确性。2 LBGK模型及边界处理2.1 LBGK模型介绍现在已经发展完善的格子Boltzmann模型中，一般包含了三个部分11：离散速度型，平衡态分布函数，分布函数的演化方程。这三大部分是彼此相关相互联系的，构造格子Boltzmann模型的关键是选择合适的平衡态分布函数，而平衡态分布函数的具体形式又与构造离散速度模型有关，离散速度的对称性决定了相应的格子Boltzmann模型能否还原到所要求解的宏观方程，因此，选择建立合理模型对准确还原到宏观方程至关重要。单松弛模型或LBGK模型是至今应用最为广泛的LBM模型，由Qian，Chen等人于1992年分别独立地提出的，是用一种简单的松弛方法来简化复杂的碰撞过程。其中最具代表性的是由他们提出DnQb模型12(n代表空间维数，b代表离散速度数)，被视为格子Boltzmann方法的基本模型，其他常用模型都是在此模型基础上发展演化而形成的，下面详细介绍一下格子Boltzmann方法常用的模型。2.1.1格子Boltzmann-LBGK方程根据微观气体动力学原理和统计学方法，可以描述微观运动的格子方程，流体特性由粒子分布函数描述，该函数给出了虚拟流体质点在离散速度上，t时刻位于网格点x处的粒子密度分布函数。格点上碰撞和迁移所引起的平衡态分布函数的演化过程是由离散的Boltzmann方程所决定的，如下所示 (i=0，1，.，8) (2.1)式中，是时间步长，这里引入了一个碰撞算子是关于的函数，表示的是碰撞造成粒子分布函数的变化所占比率，定量描述了分子间的微观相互作用。在LBGK模型中利用了一个简单的碰撞算子代替LBE中原有碰撞项：(i=0，1，.，8) (2.2)式中，为速度平衡态分布方程，为无量纲的弛豫时间，只与流体的物性参数(如：粘性系数、热传导系数、流体的质量扩散系数等)有关。将LBE方程进行离散，粒子的速度在相空间离散为有限维的速度空间，由此可以得到含有外力项的格子Boltzmann方程的有限差分格式： (2.3)此方程极大简化了原来的LBE方程，而且适用于不同平衡态分布函数和多物理场问题，方程具有Lagrange特性，一阶精度的LBGK方程实际上具有二阶精度，令该模型的计算精度基本达到二阶精度，使很多很复杂的边界条件都能够很容易得到解决，因此该模型具有较强的边界问题处理能力。2.1.2 DnQb离散模型在DnQb系列模型中，n代表空间维数，b代表离散速度数，根据n和 b的不同，可以分为D2Q5、D2Q9、D3Q15、D3Q19模型等，这些模型中所采用的建模控制方程和平衡态分布函数的求解方程形式均相同，区别在于离散速度的具体配置情况。本文所采用的是D2Q9模型。在此模型中，我们将不可压缩流体流动和一个九速度模型表示在一个二维格子(D2Q9)上，网格采用标准的正方形格子，根据速度大小和方向，其离散速度配置如图2.1所示。其中D代表空间维度，而Q则代表在一个节点的不同速度数量。每个格点只与相邻的8个格点接触，LBGK方法的思想是，用粒子分布函数来代表各个格点的流体，在流体流动过程中，在该区域内的流体离散分向9个方向，当i=0时，fi停留在原地，而当i=18时，则粒子将沿如图2.1所示的8个方向以速度ei运动，经过时间步长到达相应格点。经历了一个后，在内部相互的作用下，各格点所代表的粒子又汇聚成了新的流体群，上述流动过程将会随着时间继续不断发生。图2.1 D2Q9离散速度模型速度配置如下：(2.4)其中，为格子速度，和分别为网格步长和时间步长，在通常情况下方向和方向的网格步长相同，即=，在模拟过程中两者的取值均为1。在DnQb系列模型中，平衡态分布函数有一个统一的表达形式： (2.5)其中，是代表粒子速度函数权系数，与声速有关(或称为格子声速)，这两个参数是决定LBGK模型的关键参数。依赖于所使用的格子类型，对于D2Q9模型相关参数设置如下： (2.6) (2.7)此模型中的宏观压力与密度之间关系如下：(2.8)模型中宏观密度和和速度可以分别由定义：(2.9)(2.10)使用Chapman-Enskog方法展开便可推导出该模型所对应的宏观方程组： (2.11)其中运动粘度系数的关系如下 (2.12)可以看出，当流体密度为常数且满足低马赫数流动时，式(2.12)可以转化为标准不可压缩流体的N-S方程组。事实上，无论是导出的LBGK方程的泰勒展开还是以上的推导过程，都有一个共同的假设，那就是马赫数即Ma充分小，上述LBGK模型的结果只限于流体的低马赫数流动过程，在流体力学中一般讲Ma小于0.3的流动视为不可压缩的，因此除了压力项外，密度为常数。换言之，LBGK方法实际上是求解不可压缩流体N-S方程的一种人工压缩法。以下列出常见的DnQb模型离散速度、格子声速(假定c = 1)和权系数的参数值。(1) D1Q3模型： (2.13)(2) D1Q5模型： (2.14)(3) D2Q7模型： (2.15)(4) D2Q9模型： (2.16)(5) D3Q15模型： (2.17)(6) D3Q19模型： (2.18)2.2 初始条件及边界条件处理初始条件和边界条件是流体力学研究中的重要课题，一般说来，初始条件和边界条件总是按着宏观量(压力、速度、温度)给出，但在LBM中基本变量是分布函数，虽然根据分布函数可以简单的确定宏观流动变量，但是根据宏观量确定分布函数却不简单，因此如何根据给定的初始和边界条件构造离散分布函数的相应条件，是LBM用于实际计算的前提，大量实验研究表明，实现初始和边界条件的不同方法对LBM的计算精度、数值稳定性和计算效率都有很大的影响。可见，初始和边界条件的处理是LBM的基本问题，本节将介绍如何确定初始和边界条件2.2.1 初始条件对稳态和准稳态流动，初始条件对最终结果影响不大，此时一般可以将分布函数设平衡态，即，其中的参数是初始时刻的宏观变量，但对非稳态和对初始条件敏感的强非线性流动(如湍流、多相流等)，如何准确实现初始条件非常关键，目前LBE方法的初始化有两类方法，一种方法是非稳态校正法，最初由Skordos P针对等温的单松弛BGK模型提出13，该方法利用Chapman-Enskog技术展开，通过对分布函数求解获得高阶项近似，进而可以得到近似的计算式。另一种方法是迭代计算方法14，其主旨是借助格子Boltzmann方程求解初始压力对应的Poisson方程，并未相应的分布函数赋初值，由此获得与速度场相匹配的初始分布函数。2.2.2 边界条件处理与初始条件类似，在应用模型时也需要给出分布函数的边界条件，对于流动问题和传热问题，边界条件起着非常重要的作用，直接影响并决定整体数值的计算精度、稳定性和计算效率15。以稳态问题为例，当系统达到稳定后，流场和温度场与初始条件无关，而主要取决于边界条件。在格子Boltzmann方法中，边界条件同样起着重要作用，给定一个初始条件和边界条件，就可以给出离散粒子分布函数的一些对应条件。边界条件处理的方式对计算的效率、计算数值的稳定性和计算结果的精度都有很大影响。应用LBM进行数值模拟，通过每个时间步长的碰撞和迁移过程可以得到内部节点的粒子分布函数，而边界节点上的粒子分布函数还并未全部已知，只有确定了边界节点上的粒子分布函数，才能进行下一阶段的计算。因此，需要从已知的宏观上的边界条件来试图确定边界节点上的粒子分布函数，在这个过程所应用的方法称为格子Boltzmann方法的边界处理方法，所设计出的计算格式称为边界处理格式。根据格式处理方法可以分为启发式格式，动力学格式和外推格式，根据边界条件的类型可以分为速度边界和压力边界，其中速度边界包括了平直边界和曲面边界。本文主要应用的是平直边界格式，下面介绍几种常用的边界格式。 (1)启发式格式启发式格式主要根据边界上的宏观特性，如对称性，周期性，充分发展等，通过粒子的运动规则来确定边界节点的粒子分布函数，与其他的边界处理格式相比，他不需要复杂的数学推导和公式求解，启发式格式主要包括周期性边界处理格式、对称边界处理格式、充分发展边界处理格式，以及适用用于固体表面的反弹格式、镜面反射格式、反弹与镜面反射混合格式等。周期性边界处理格式：数值模拟时，如果流场在空间内呈现周期性变化或在某一方向无限大时，则可以将周期性单元取出作为模拟区域，并在相应的那些边界上采用周期性边界。能严格保证整个模拟区的质量和动量守恒，是一种最容易实现的边界处理格式，其基本思想是假设粒子从一端离开流场后，在下一个时间步长以由其相反的另一端对称位置进入流场另一边界，整个计算区域保持封闭再循环。周期性边界处理格式相比于其他常用的边界处理格式，在数值计算的稳定收敛性方面保持最好。图2.2 周期边界示意图对称边界处理格式：对于对称性问题，处于节约资源的考虑，可以取其物理模型的一半作为模拟区域，并在对称轴上采用对称边界处理格式。充分发展边界处理格式：针对的是流体在管内流动达到充分发展段后，此时流体的速度和密度在主流方向上都不再发生改变，因而可以采用这种数值稳定性较好的边界处理格式。反弹边界处理格式：也成为无滑移壁面边界处理格式，是常用来处理简单边界(如固体壁面)的格式，也是一种典型的处理静止流体与固体接触面边界处理格式。包括标准反弹格式、半步长反弹格式及修正反弹格式。基本思想是假设壁面处的流体粒子宏观速度为零，该格式假设粒子以一个速度到达边界，并在下一时间步长以相同的速度原路返回，此过程中没有能量损失，且边界节点不发生碰撞过程。该格式易于操作计算简单，且能满足系统的各物理量守恒，但是精度不高，只能达到一阶精度，且对运动边界等复杂边界很难处理，为了改进这一缺点，人们提出了修正反弹格式和半步长反弹格式：将物理壁面移至格点中间，在边界与格点之间应用修正反弹格式和半步长反弹格式。这两种改善格式能够达到二阶计算精度，可以有效弥补标准反弹格式精度的不足。反弹格式在处理复杂的不规则界面具有很大优势。图2.3 反弹边界示意图(2)动力学格式动力学格式主要包含Nobel格式、非平衡反弹格式等，其是直接求解边界格点粒子分布函数通过宏观物理量的定义。Nobel处理格式：通过构造压力约束来求解边界上未知粒子的分布函数，其基本思想是：在边界节点速度已知的情况下，根据边界上分布函数与宏观量之间的关系，确定边界节点的分布函数，不管是速度边界还是压力边界，实际上就是对流体施加一个外力作用，引起动量变化，非平衡反弹格式：1997年Zou和He对二维和三维的LBGK模型进行了研究，提出来这种新的边界处理格式，在实验研究过程中通过与其他处理格式比较，发现这种格式具有良好的数值稳定性。(3)外推式格式外推式格式是将边界扩展，假设边界上的点也与流场内的所有点一起参与碰撞，然后通过线性外推来获得边界上的粒子分布函数，该方法虽然可以达到二阶精度，但是其数值稳定性较差。(4)Chen格式Chen格式在时间和空间上都具有二阶精度，很好的保持了格子Boltzmann方法的整体精度，无需对边界的流体分布函数作任何假设，因此适用非常广泛，但也有其缺点，就是数值稳定性较差。(5)非平衡态外推格式非平衡态外推格式是基于外推格式和非平衡反弹格式，郭照立教授等人于2002年提出了非平衡态外推格式，这种格式将边界节点上的分布函数分解为两部分，平衡态和非平衡态，即和，平衡态分布函数，可以通过相应具体给定的边界条件定义近似获得；而非稳态分布函数，可以直接通过外推方法来确定。在D2Q9模型中详细介绍该方法，如图2.1，假设301是边界，625位于流场内，748位于流场外，t是当前时刻，已知2点粒子分布函数，求0点粒子分布函数，该方法将0点粒子分布函数分为两部分： (2.19)节点2粒子分布函数，宏观速度及密度都已知，因此2点非平衡态粒子分布函数可以计算为 (2.20)0点平衡态部分近似为 (2.21)通过非平衡外推可得0点粒子分布函数 (2.22)非稳态外推格式在时间和空间上都能够达到二阶精度，同时具有很好的数值稳定性，并且兼顾了其他外推格式的诸多优点，且其计算不复杂，较容易实现，因此适用范围十分广泛。2.3 格子Boltzmann方法的几种作用力模型在格子Boltzmann方法中，粒子的速度被离散化，而粒子速度的改变则体现了外力或内部力对粒子运动的影响。因此，如何在离散速度的框架下描述作用力是一个十分重要的问题。LBGK的作用力模型一般有三种，即平衡态分布的压力校正方法、平衡态分布的速度校正方法和在演化方程中增加作用力项。2.3.1平衡态分布的压力校正方法应用平衡态分布的压力校正方法必须满足的条件上，作用力可以表示成势函数的形式，同时作用力导致的密度变化可以忽略不计。LBGK模型中的平衡态粒子分布函数是密度和速度的函数，可统一表示为 (2.23)其中，为权系数，与声速有关。标准的LBGK方程对应的宏观流动方程如下：(2.24) (2.25)其中 (2.26)将作用力表示为势函数的形式，即。假设密度空间的变化可以忽略不计，则可得。因此，可通过适当修改平衡态粒子分布函数将状态方程修改为(2.27)变换过程如下所示： (2.28)这样就得到了包含作用力F的LBGK模型。2.3.2平衡态分布的速度校正方法平衡态分布的速度校正方法的主要思想是希望改变平衡态粒子分布函数中的宏观速度来研究作用力的影响。假设发生一次碰撞的平均时间是t，则由作用力引起的液体动量变化为Ft。为考虑这种影响，Shan-Chen模型采用如下平衡态粒子分布函数： (2.29)其中称为“平衡态速度”(2.30)为了得到更精确的结果，流体速度v应为碰撞前后速度及=的平均值，即 (2.31)2.3.2在演化方程中增加作用力项这种作用力模型可用下式表示：(2.32)其中，F表示作用在流体粒子上的总作用力F的离散形式，它不能直接对连续的作用力项进行离散处理。这类方法的代表有基于LGA的模型(包括LGA模型和修正的LGA模型)、Buick-Greated模型、He-Shan-Doolen模型、二阶矩模型及修正的二阶矩模型等。2.4 格子Boltzmann方法的计算步骤应用LB方法模拟具体流场的算法程序流程图如图所示。格子Boltzmann方法常用的离散模型有D2Q9(二维正方形网格)，D2Q7(二维正六边形网格)，D3Q15(三维正方体网格15点模型)，D3Q18(三维正方体网格18点模型)等。在计算开始之前，需要选择适合的网格类型。图2.4 格子Boltzmann算法程序流程图在确立了DnQb网格模型后，对粒子密度分布函数赋初值并计算对应的密度和速度，从而确定平衡态粒子分布函数，进行碰撞和迁移过程。和只进行迁移过程的边界节点不同，内部节点同时进行碰撞和迁移过程。然后选用合适的边界条件进行边界处理。常用的边界条件有周期边界、反弹边界、速度边界、外推边界等。在完成了边界节点的处理后，由得到的粒子分布函数来计算速度、密度和压力等宏观物理量。至此主要计算完成。若模拟的物理过程为非定常过程，则需要判断是否达到循环达到要求次数来确定输出结果还是继续循环。若流场属于定常问题，则需要判断相邻时间步长的各物理量是否达到收敛条件，若满足收敛判据，则输出结果，否则，继续进行循环。3 非牛顿流体的格子Boltzmann模型3.1 非牛顿流体概述3.1.1 非牛顿流体的基本特征流动是一种连续变形的运动，通常把不符合牛顿内摩擦定律的流体称为非牛顿流体。物体受外力作用，当外力除去后能恢复原状的变形，则称为弹性变形，不能恢复原状的变形称为塑性变形。有的非牛顿流体既有粘性又有弹性，称为粘弹性流体。非牛顿流体的种类很多，如血液、蛋白液、油漆、泥浆以及高分子聚合物的溶液等。还应该指出的是，非牛顿流体一般均为液体，也就是不可压缩流体。3.1.2 非牛顿流体的粘性流体的分子结构决定流体的粘性。牛顿流体的额分子结构比较简单，所以用牛顿内摩擦定律表征的粘性系数是温度与压力的函数，与流动状态无关。非牛顿流体的物质构成要比牛顿流体的要复杂得多。在粘性的反映上，非牛顿流体的粘性不仅与温度、压力有关，而且还与流动状态有关，一般来讲很难用单一的系数表征其实际的粘性，通常情况下，以牛顿内摩擦定律为模式来定义非牛顿流体的粘度，用来表示：(3.1)为剪切应力，为剪切速率，称为非牛顿流体的表观粘度。这里不同于牛顿流体，其中不是常数，而是剪切速率的函数。根据表观粘度随剪切速率而变化的关系可将非牛顿流体分为：假塑性流体、胀塑性流体、粘塑性流体、触变性流体和粘弹性流体。3.1.3 非牛顿流体的弹塑性非牛顿流体一般为液体，属于不可压缩流体。以牛顿流体的观点而言，不可压缩流体是没有弹性的，流体的流动变形是不能恢复的，但是，许多非牛顿流体不仅具有粘性，而且还有弹性，这就是说，非牛顿流体的变形，既有不可恢复的变形部分，也有可恢复的变形部分，在流动变形过程中，具有弹性恢复效应，或出现应力松弛现象，这种非牛顿流体称为粘弹性流体。在外力作用下，固体一般先产生弹性变形，当应力超过屈服极限时就产生塑性变形，不是有任意微小剪切力作用，就会产生流动变形，它有一个屈服限，称为屈服应力。当剪切力小于屈服应力时，此类流体是不会流动的，只有当剪切应力超过屈服应力值时候，此类流体才会流动，像牛顿流体一样流动，这种流体称为宾汉姆(Bingham)流体16。在LBM中，对于非牛顿流体，由于其运动粘度与剪切速率有关，所以每个格子节点的粘度都要根据与剪切速率关系进行计算，计算出粘度之后，才能得出各个节点的松弛时间。以下，给出两种非牛顿流体的本构方程。3.2 幂律流体与宾汉姆流体本构方程许多种非牛顿流体的实验表明，剪切应力虽然与剪切变形速率不呈线性函数关系，但是可以表示为简单的幂函数关系。在简单的剪切流动中，本构方程只需要剪切应力与剪切变形速率k来表示即可。这类流体的本构方程是(3.2)式中，为剪切应力；为稠度系数，或称幂律系数，是粘度的度量，其随着粘度的增大而增大；n代表流性指数，又称幂律指数，反映偏离牛顿流体的程度。适用于这类本构方程的流体称为幂律流体。工程应用中常将(3.2)式改写成类似牛顿流体的本构方程(3.3)式中(3.4)称为幂律流体的表观粘度。表观粘度是剪切变形速率的函数。当n=1时，流体即为牛顿流体。当n值越低或越高，非牛顿性越强，切应力与剪切速率曲线也越弯曲。当n1时，表观粘度随剪切变形速率的增大而增大，称为膨胀型性流体，也称为剪切增稠流。宾汉流体的表观粘度为：(3.5)可以看出，宾汉流体的表观粘度是随流速梯度而变化的。宾汉姆塑性体是一种简单的也是比较常见的粘塑性流体，根据塑性流体的流变曲线，可以写出如下关系式:(3.6)式中：为极限动切应力，称为结构粘度(或称塑性粘度)。上式称为宾汉方程，符合宾汉方程的流体称为宾汉流体，塑性流体也称为宾汉流体。对于宾汉流体而言，当其所受到的切应力小于屈服应力时，流体表现出固体性质，不会流动。当切应力大于屈服应力之后，流体开始流动，其剪应力与剪切速率的关系与牛顿流体类似。可将Bingham流体的本构方程写成 (3.7)式中，0-屈服应力；p-塑性黏度。其中是D1的不变量，通常定义为(3.8)实际上不能精确的确定，它的值可能依赖于用于确定它的仪器。将Papanastasiou(修正的宾汉姆)模型引入LBM模拟中，为了克服不连续性这一缺点，Papanastasiou