参考文献

第5章　基于浸入边界-格子Boltzmann方法的建模

5.1　引　　言

第4章详细介绍了偏微分方程的数值解法，重点介绍了有限元方法，并采用COMSOL Multiphysics举例，采用“step by step”的方式演示了一些微流控芯片单元的建模。其中，重点介绍了流动分析、电泳进样、微混合优化、多相流液滴操控等内容，可见采用商用软件可以开展丰富的微流控芯片建模研究。尽管现有的有限元模拟软件能够胜任很多微流控芯片技术中的建模分析工作，但也存在一些功能上的短板。例如，在微流控芯片技术中，细胞分选或颗粒惯性聚焦是常见的操作模式。然而，在细胞(或颗粒)建模过程中，涉及大量复杂运动边界的模拟，这在现有的有限元软件中难以实现。

目前，格子Boltzmann方法(Lattice Boltzmann Method，LBM)在基于微流控芯片技术的细胞和颗粒操控中得到较多的应用[1-4]。但和有限元法不太一样，LBM是一种新兴的方法，市面上可用的软件非常少，其应用还主要通过编程来实现。另一方面，LBM比较灵活，易与其他方法(如浸入边界法、相场法等)相结合，形成一种具有特定功能的数值模拟框架。

编程的核心是算法。只要算法逻辑正确、步骤清晰，理论上采用哪种编程语言都可以。目前常见的算法编程语言有C/C++、Fortran、Java、Python、MATLAB等。C/C++适用面广，几乎所有操作系统都支持C/C++，其跨平台性好，计算速度比较快，但学习难度较大。Fortran语言的强项在于有长期的积累，语法简单，擅长数值计算、并行计算；弱项是其界面功能、交互功能、图形功能稍弱。Java应用范围广泛，比C/C++容易学习和理解，且提供了对Web应用开发的支持；属于解释型语言，运行效率低。Python是一种面向对象的通用开源脚本编程语言，语法简单易学，贴近人类语言，有大量第三方库，功能扩展较好；属于解释型语言，计算效率低。MATLAB是一种面向科学与工程计算的高级语言，允许用数学形式的语言编写程序，且比C语言等更加接近计算公式和人们的思维方式，它编程简单、擅长矩阵运算、易学易懂、工具箱丰富，但MATLAB也是一种解释性语言，其计算效率低。

在这些语言中，MATLAB相对简单易学，适合有理工科专业背景的人员学习和使用，尤其是在学习和理解算法方面，选择MATLAB会更高效。MATLAB的线上和线下学习资料都很丰富，因此本书不单独设置编程语言介绍模块，请读者自学。此外，MATLAB是商业软件，与之语法类似的开源软件有很多，推荐使用Octave。Octave为GNU项目下的开源软件，早期版本为命令行交互方式，4.0.0版本发布基于Qt编写的GUI交互界面。Octave语法与MATLAB语法非常接近，可以很容易地将MATLAB程序移植到Octave。

5.2　格子Bo Itzm ann方法的发展历史

格子Boltzmann方法是一种基于介观(mesoscopic)模拟尺度的模拟计算方法，主要用于计算流体力学的建模分析。这里的介观方法意为LBM在演进中涉及宏观和微观变量的相互转化，因此整体上属于介观层次方法。该方法得到的数学形式采用多尺度展开技术，可以得到与宏观的偏微分方程相同或相近的形式，因此LBM也被一些学者认为是偏微分方程的求解器。但当偏微分方程形式不同时，LBM的数学形式也不相同，所以该方法具有较高的灵活性。如果让LBM模拟流体流动，那么它的数学形式需要能够通过多尺度展开技术恢复到Navier-Stokes方程。LBM主要通过碰撞和迁移的交替演进来获得流体流动的规律，相比于其他传统的计算流体力学方法，LBM具备算法简单、易于程序实现、易于处理复杂运动边界、易于实施并行计算等优势。目前，LBM已经被广泛认为是描述流体运动与处理工程问题的有效手段，在模拟流动、流固耦合、传热、相变、对流扩散等领域获得了很好的应用。

格子Boltzmann方法最初源于格子气自动机(Lattice Gas Automata，LGA)的研究。LGA将流体视为许多只具有相同的质量而忽略体积的粒子组成，该方法在时间、空间和状态上对连续的流动介质进行离散。这些离散化的微粒根据事先给定的规则在网格点上进行碰撞，所有粒子的碰撞都遵循质量、动量和能量守恒定律。发生碰撞之后的粒子会在其所在的网格线上向相邻的节点迁移，从而描述流动的规律。

1.元胞自动机与格子气自动机

格子气自动机(LGA)实质上源自元胞自动机(Cellular Automata，CA)，前者是后者的在流体力学中的运用。

CA最初由Sanislas Ulam、John von Neumann和Konrad Zuse等学者于1950年左右提出，是一种物理系统在空间及时间上都进行离散化的数学模型。在空间上，CA将物理系统离散化成具有规则形状的格子，每个格子点都称为细胞，且格子点的状态都用一个变量来定义。在时间上，物理系统演化的时间被离散化为多个时间步长，每个格子的状态都根据时间步的推进而发生变化。下一时刻的某个格子所处的状态来源于与其相邻的格子当前时间步所处状态。

一般情况下，CA可以定义为一个四元组A=(C，S，N，f)。其中，C表示物理系统离散化后的空间点阵列，也可以称为细胞阵列，它可以是一维(或多维)的空间点阵；S是格子点可以取得的所有状态的集合，可以称为状态空间；N为与当前格子点相邻的格子点的集合；f为格子点状态的演化函数。

对于任意一个格子点P∈C，在时刻t有P(t)∈S，在t+1时刻，该细胞的值由邻域N(P)={P1，P2，…，Pn}∈S在t时刻的状态确定，即P(t+1)=f(N(P)(t))。根据此式，细胞运动状态便随着离散的时间t=0，1，2，…演化下去。

CA有几个显著的特征：能够大规模同步运行；局部相互作用；细胞结构非常简单。

CA能够满足各种不同结构计算系统的并行需求，特别是在那些分布式的计算机上。CA的计算区域可以分为多个相互独立的子区域，每个子区域内可以独立演化，而相邻子区域之间只需要交换边界上的状态数据。正因为如此，CA特别适合解决大空间规模的数学物理问题。由于CA的演化规则十分简单，却能反映复杂的物理现象，因此引起了各领域学者的广泛兴趣。

此后，很多学者对CA做了大量研究。例如，Broadwell于1964年提出了离散速度模型，并将这一模型用于对基波结构的研究，用其求解和描述一些较简单流体的Boltzmann方程。该模型首次提出将速度仅离散成几个方向，但在时间和空间上依然保持连续。

2.格子气自动机

相比于元胞自动机(CA)，格子气自动机(LGA)是一个更大胆的尝试，其不只将流体的空间离散，还将时间进行离散化处理。这样的处理使相应算法更加易于计算机实施，为计算流体力学提出了一个全新的方向。

Hardy、Pazzis和Pomeau在1973年提出了HPP模型，他们在原先速度离散模型的基础上进行了改进，提出了以正方形网格为基础的模型。通过运用此模型，他们成功地模拟了简单的流动现象。HPP模型能够满足质量和动量守恒定律，被公认为是第一个格子气自动机(LGA)模型，为格子Boltzmann方法的发展奠定了基础。但这个模型仍然存在一些问题，如由于格子缺乏各向同性，这使得模型从宏观上恢复Navier-Stokes方程时不能正确地反映非线性项和耗散效应项。

Wolfram于1983年提出新的布朗型动力系统，其在研究一维的8位元胞自动机时，发现在已经给出的初始条件下，该系统能按照事先制定的规则进行演化。在该系统中，每个细胞的运动仅与周围的细胞有关系，所以该系统具有很强的局部性，从而可以非常方便地进行并行处理。Wolfram的想法极大地促进了格子气自动机的发展。

Frisch、Hasslacher和Pomeau等学者于1986年提出了新的格子气自动机模型，该模型具有对称的二维正六边形网格，被命名为FHP模型。该模型的每个格子上都有6个方向，每个格子上仅有一个粒子，或者没有粒子。该模型中格子具有很好的对称性，很大程度上弥补了HPP模型的缺陷，基于FHP模型可以推导出Navier-Stokes方程，开创了格子Boltzmann方法研究的新局面。

格子气自动机的出现为人们研究物理问题提供了新的方法与思路，发展出的HPP模型及FHP模型，使其在模拟流动方面有了很好的应用前景。但LGA仍然存在一些严重的缺点，主要包括4个方面：

(1)由格子气自动机演化方程推导出来的动量方程不满足伽利略不变性，其主要原因是对应宏观方程的对流项与密度有关。

(2)状态方程含有速度项，与实际情况(理想气体状态方程P=ρRT中不含速度项)不符。

(3)由于采用布尔运算，因此局部量一般含有数值噪声，需要对时间和空间做平均，从而增加计算量。

(4)碰撞算子具有指数复杂性，对计算量和存储量的要求很高。

3.格子Boltzm ann方法(LBM)

为了保留格子气自动机模型的优势并消除其方法本身所存在的缺陷，学者们不断进行探索和尝试，最终促使了格子Boltzlmann方法的诞生。

1988年，McNam ra和Zanetti提出了将LGA(格子气自动机)中的布尔代数运算变成实数运算，每个格子点上的微粒数目不再用整数0或1来表示，而是采用实数f来表示整个系统所有微粒做平均后的局部分布函数，用Boltzmann方程代替LGA的演化方程，并将该模型用于流体的数值计算。这是最早的格子Boltzmann模型，从此开启了格子Boltzmann方法的大门。

1989年，Higuera和Jimenez提出了一种简化模型：引入了平衡态分布函数，将碰撞算子线性化成分布函数向平衡态演化的过程。这个模型不需要进行碰撞运算，并忽略各个粒子之间的碰撞细节，相比于之前的多粒子碰撞模型，该模型构造起来十分容易。同年，Higuera等人又进一步提出了强化碰撞算子方法，以加强模型的数值稳定性。他们提出的这两种模型统称为矩阵模型。

经历了以上两类模型的发展，格子Boltzmann方法消除了统计噪声，克服了碰撞算子指数复杂性，但此时的平衡态分布函数依然采用Fermi-Dirac平衡态函数，格子气自动机的其他缺点仍然存在。

1991年，陈十一等人提出采用Maxwell-Boltzmann分布替代Fermi-Dirac分布，不但解决了非伽利略不变性问题，而且解决了压力状态方程不独立于流速的问题。1992年，钱跃竑等人引入Bhatnagar-Gross-Krook的使松弛时间控制局部分布函数逼近相应平衡态的LBGK模型，使得碰撞项变得非常简单，至今仍广泛应用。

1992年，D'humeriers提出广义LBE，即多松弛模型(MRT原型)；2000年，Lallemand和Luo论证并推广MRT模型，拓展了LBGK的雷诺数的应用范围，增强了LBM稳定性。2002年，郭照立提出D2G9模型(不可压缩模型)，有效消除了D2Q9(伪不可压缩)带来的压缩性误差，可适用于非定常问题。近十余年来，LBM算法研究则主要在于借助传统有限元、有限差分、有限体积法等的优点与LBM的融合，拓展LBM应用能力。此外，一些更加稳定和高效的LBM被陆续提出，如2008年Irina Ginzburg等人提出的双松弛因子LBM、2015年Succi等人提出的量子化LBM等。总而言之，LBM的发展方兴未艾，不但其算法本身在不断改进和提高，而且其应用领域不断向多个学科领域和工程领域渗透，是一种极有前景的数值模拟方法。

5.3　格子Bo Itzm ann方法的基本原理

LBM模型的核心概念主要有概率密度函数、格子模型、宏观量计算方法、平衡态分布函数、碰撞、迁移及边界条件设置。

5.3.1　概率密度函数

Maxwell-Boltzmann分布律的公式为

基于统计力学观点，任何(宏观)物理系统的温度都是组成该系统的分子和原子的运动的结果。这些粒子有一个不同速度的范围，而任何单个粒子的速度都因与其他粒子的碰撞而不断变化。然而，对于大量粒子来说，处于一个特定的速度范围的粒子所占的比例几乎不变，即

不发生改变。这里的f(v)就是一种“概率密度”，它的积分(从0到∞)应该等于1。

LBM中的概率密度函数在本质上和f(v)是一样的，它描述在不同方向上粒子速度的分布概率。这种概率分布的变化既反映压力的改变，也反映速度场的变化。事实上，利用概率密度函数与空间向量的更高阶矩关系组合，还可以得到其他流动演化的一些性质，如高温差、高压差下的非平衡态过程体现出的规律。LBM的概率密度函数的形式依赖于格子模型。例如，在D2Q9格子模型中，一个网格节点有9个空间方向，在这9个空间方向上都存在一个粒子速度分布概率。采用向量F来表示概率密度函数，那么在一个节点上，F将具有9个元素，这9个元素之和以1为数学期望，所以F也有“概率”的内涵，相应元素的描述函数又称“概率密度函数”。

对于处于非稳态的节点，其9个方向的概率密度函数将根据一定的规则演化(碰撞、迁移，这在后面会进行阐述)；对于达到稳态的节点，其9个方向的概率密度函数将稳定不变，此时的概率密度函数与该点的平衡态分布函数(后面会阐述平衡态分布函数的概念)相等。

概率密度函数是LBM中的“微观量”，但由这个微观量可以推导出速度、压力等宏观量。因此，概率密度函数的概念是理解LBM原理的重要基础。

5.3.2　格子模型

根据LBM离散速度设置的差异，格子模型也存在差异。在一维情况下有D1Q3、D1Q5；在二维情况下有D2Q7、D2Q9等；在三维情况下有D3Q15、D3Q19等。D指代Dimension，D后面的数值是问题的维数；Q是指钱跃竑老师姓氏的首字母，以纪念他为LBM速度离散模型做出的贡献，Q之后的数值为一个节点包含的离散速度的方向数量。格子模型包含三个主要概念，即格子向量、格子声速及格子权重。下面，以常用的D2Q9格子模型为例进行说明。

格子向量如图5.1所示，除8个箭头方向之外，第9个方向为中间点。第1～4个方向可以表示为，i=1，2，3，4；第5～8个方向可以表示为，i=5，6，7，8；第9个方向可以表示为e9=(0，0)。通过格子向量，可以对不同方向的概率密度函数进行选择性操作。

声速的定义：在0℃下，声音在一个标准大气压的空气传播速度。它既是一个场量传播速度的标准，又是一个物质传递模型的基本量，大约为334 m/s。在LBM中，声速c=Δx/Δt，即场量在场域中的固有传播速度，各向同性，一般取c=1。格子声速c s是离散Boltzmann方程恢复到Navier-Stokes方程中的一个量，来自理想气体状态方程P=ρRT，其中c s=。

格子权重ωi的导出以格子张量计算规则和质量、动量及能量守恒规则为依据，要求在具体格子向量设置下，相应的权重组合能将离散Boltzmann方程恢复至宏观的Navier-Stokes方程。不同格子模型的权重不同，D2Q9格子权重模型如图5.2所示。

图5.1　D2Q 9格子向量模型

图5.2　D2Q 9格子权重模型

下面用MATLAB代码设置LBM的格子模型，包括：格子向量、格子声速及格子权重。代码如下：

5.3.3　宏观量计算方法

LBM中，宏观量是指流场反映出的可观测的物理量，如密度、压力、速度、温度等。相比之下，概率密度函数、网格模型中的参数等都是微观量。正是由于宏观量和微观量之间的交替演化发生才使得流场在边界条件及源项(如流场受到的力源、热源作用)影响下演进。

在D2Q9格子模型中，根据格点上概率密度的宏观量计算方法如下：

式中，fi ——　在9个方向中第i方向的速度概率密度。

由动量守恒定律可得速度：

由状态方程可得压力：

有了以上公式，接下来介绍如何进行代码实现。以通道流动建模为例。假设二维的矩形流动区域为Ω，长×宽采用lx×ly个方形网格来进行建模，F是一个三维数组，存储所有网格上的概率密度函数，其中第一维对应通道长度方向的网格序列，第二维对应通道宽度方向的网格序列，第三维对应每个节点上的9个方向，因此F是一个lx×ly×9的数组。

用MATLAB代码可将式(5.3)～式(5.5)表示如下：

代码说明：

第1行，表示将F中第三维的所有元素相加，以得到密度rho，执行这个操作后，返回的rho是一个维度为lx×ly的数组。第2、3行，用于计算速度。由于是二维情形，所以速度有两个分量，一个是x方向速度分量UX，另一个是y方向速度分量UY，这两个速度分量可以合成任意方向和大小的速度。需要注意的是，在这两行代码中，我们将方向向量e x、e y的数值-1、0、1直接分配到各项，最后体现为代数加减。例如：e x=[1　1　0　-1　-1　-1　0　1　0]，表示第1、2、8个方向都为1，第4、5、6个方向都为-1，其余方向为0。与为0作乘积的代数项则不进行计算，这样可以减少计算量。此外，第2、3行还出现了“./”这个组合符号，该符号在两种情况下使用：其一，分子为一个数、分母为一个数组时，进行“./”运算的结果是返回与分母维度一样的数组，所返回数组元素为分子除以分母相应位置元素的值，如2./[5，6]=[2/5，1/3]；其二，分子、分母均为同维度数组时，“./”返回对应位置元素相除的结果，如[2，4]./[3，2]=[2/3，2]。

5.3.4　平衡态分布函数

平衡态为系统在某一时刻具有相对稳定的形态，即一个系统在没有外力的作用下随时间的变化熵不再减小时的状态。在LBM中，平衡态分布函数把节点当时的宏观量(密度、速度、压力、温度等)认为是“系统的熵不再随着时间推移减小”时的状态，即稳态。但真实情况下也许并不如此，在演化过程中，也许这里的宏观量仍然处于非稳态。由于平衡态分布函数的计算依赖于当地节点在某一时刻的具体宏观量，因此从这个意义来说，平衡态分布函数只反映一个“即时演化标准”。

1872年，Boltzmann在维也纳科学院的会议报告上发表了著名的论文Weitere Studienüber das Wärmegleichgewicht unter Gasmolekülen(《气体分子热平衡进一步研究》)。其中，他引入了量H(t)，其为

并证明了H(t)必定随时间单调递减，且在f(E，t)为Boltzmann分布时达到极小值。根据Boltzmann-H定理，对于稳态气体系统，其分子运动规律的统计平均都能得到一个确定的值，即系统处于H(t)为极小值时的状态，称为平衡态。在平衡态下，网格点上的概率密度函数fi称为平衡态分布函数，记为。

图5.3　Maxwe ll-Boltzm ann速率分布曲线

接下来，推导LBM的平衡态分布函数计算公式。图5.3表示了温度为T的一个气体系统，气体分子的速率满足Maxwell-Boltzmann分布。Maxwell-Boltzmann速率分布曲线描述了系统里气体分子的微观速度和宏观速度的关系。f(u)表示概率密度函数；阴影区域表示在一个确定的时间微元和空间微元中，速度为u的分子数占总分子数的比例(概率)。相应的Maxwell-Boltzmann分布计算公式如下：

式中，ξ——微观的分子速度；

u——宏观的分子速度；

ξ-u——描述分子的热运动。

Maxwell-Boltzmann分布是气体运动理论中最基本的物理定律，描述了气体处于平衡状态时宏观与微观速度的最概然分布，分子束技术可以足够精确地测量出分子速度的Maxwell-Boltzmann分布，证明了其正确性。

对式(5.7)进行离散，将指数项进行泰勒展开，即可得在满足Maxwell-Boltzmann分布下的概率密度函数，即平衡态分布函数(i=1，2，…，9)：

式中，N为单位体积的分子数量，在宏观下为密度ρ；ξ被简化为有限维的速度空间e i；ωi=(2πRT)，为方向权重，其值由恢复到Navier-Stokes方程时，求解参数方程组得出。

平衡态函数用在演化方程中的碰撞项中，是计算过程中多时刻流场逼近的目标。本身也是一个随时间动态更新的值，当流场达到稳态时，与当时的fi相等。

平衡态分布函数在MATLAB中的代码如下：

其中，第1行的function FEQ=FEQf(UX，UY，rho，t1，t2，t3，c_squ)表示采用一个子函数来计算平衡态分布函数的值。其中，UX、UY、rho、t1、t2、t3、c_squ是计算中所需的宏观量和参数，返回的平衡态分布函数结果存储在FEQ中。

关于上述代码有以下五点说明：

(1)MATLAB子函数可以单独存储为一个独立文件。由于该文件的文件名需要和函数名相同，因此只有存储为“FEQf.m”才能被其他程序识别和调用。

(2)MATLAB子函数计算过程中的物理量在相应工作空间“workspace”中不保存，但在子函数内部是可以调用的。

(3)代码中的第2～4行为定义的过程量，这个过程量的计算结果，可以用于后续代码的计算。这样做的好处是可以简化代码和减少重复性计算。

(4)对于求取平衡态分布函数的9个方向的值，既可以如上述代码一样分别计算，也可以采用循环的方式计算。分别计算的好处是可以对不同方向平衡态分布函数的计算方法区别对待，减少不需要的计算，并且避免循环操作(MATLAB中的循环操作比较耗时)，其缺点是所需的代码多一些；而采用循环操作时，不需要区别对待具体方向都有哪些代数项参与计算，代码比较精简。

(5)应注意FEQ(：，：，1)的表示方法，其中“：”表示在某个维度上的所有序列。

5.3.5　碰撞

碰撞和迁移是LBM的两个核心步骤。

碰撞是指一个节点的概率密度函数在时间层的演化。具体可以描述为：在确定的某一时刻、某一空间点，依据当时、当地的平衡态函数，获得该点各方向速度概率分布的修正量，从而实现本点概率密度函数的演化。对于单松弛模型，碰撞项仅用一个松弛因子τ控制，公式如下：

式中，fi(x，t)——t时刻的概率密度函数；

fi(x，t+Δt)——t+Δt时刻的概率密度函数；

(x，t)——t时刻的平衡态分布函数；

τ——　松弛因子，τ=0.5+3υ，　υ为运动黏度。这里的τ的取值有一定的范围，原则上0.5＜τ＜3。这是LBM模拟不可压缩流体的要求范围，如果对此不注意，就可能导致错误的计算结果。

碰撞步骤的MATLAB代码如下：

5.3.6　迁移

迁移是流场在空间上的演化过程，如图5.4所示。迁移的作用：在确定的某一时刻和某一空间点，将该点经过碰撞的各向概率密度函数转移到相邻网格节点的相应位置，并替换该位置上的值；同时，该点也接受其他相邻网格节点转移来的值，并替换相应的本地值。在所有值迁移替换之后，就可以更新流场的密度、速度等，从而实现流场在空间上的演化。迁移公式为

图5.4　迁移示意图(书后附彩插)

(a)迁移前；(b)迁移后

注意，第9个方向(中心点处)的速度分布量在迁移中不发生变化。因此，在具体实施迁移步骤时，只考虑8个方向速度分布量的改变。相应的网格参数如图5.5所示。

图5.5　概率密度方向与概率密度分布

(a)概率密度方向；(b)概率密度分布

迁移的实施过程：假设图中位于位置“9”的网格点坐标为(i，j)，那么思考一下，f1应该迁移到什么节点呢?f2又应该迁移到什么节点呢?显然，f1应该往右迁移一个格子，所以需要将其赋值到(i+1，j)点；f2应该往右上方迁移，即需要迁移到(i+1，j+1)点。其他方向的迁移也是类似步骤。因此，从迁移的实施过程看，其本质是数组的“移位”操作。在程序实现中，可以采用循环的方式来进行移位操作，也可以用操控数组元素地址的方式来实现。对于后者，相应的MATLAB代码如下：

以第2行代码“F1(：，：，4)=F1([2：lx 1]，[ly 1：ly-1]，4)”为例，由图5.5可知，第4个方向为左上方位的方向，这个方向的值需要由点(i+1，j-1)的第4个方向迁移得到。代码等号左边表示接收迁移的节点，等号右侧则是送出迁移值的节点。这里也注意到，等号左侧、右侧均是对F1这个数组进行操作，但这不会影响修改其局部的元素值。下面来分析这一行代码是如何实现“移位”赋值的。在x方向迁移操作中的移位赋值如表5.1所示，在y方向迁移操作中的移位赋值如表5.2所示。

表5.1　在x方向迁移操作中的移位赋值

表5.2　在y方向迁移操作中的移位赋值

然而，以上移位操作也有一定的问题。例如，对于x方向的最后一个移位，以及对于y方向的第一个移位，都存在一个“首尾跨界”的问题。对于该问题，在边界处理时进行重新赋值，覆盖即可解决。

5.3.7　边界条件

研究具体的物理系统，就必须考虑对象所处的特定“环境”，而周围环境的影响常体现为系统边界上的物理状态，即边界条件。根据边界条件的设定形式不同，一般可分为三类：狄利克雷(Dirichlet)条件、诺伊曼(Neumann)条件及罗宾(Robin)条件。这在4.2.3节中已经介绍，不再赘述。

在LBM中，根据变量类型不同，边界条件分为宏观边界条件和微观边界条件。宏观边界条件主要指速度和密度(压力)。对于同一边界上的网格节点，既要设置速度，又要设置密度。原则上，速度和密度不能同时设为第一类边界，如果速度(密度)指定了，那么密度(速度)就需要采用外推方法来确定。这里的“外推”，是指将内侧一层网格的相应值赋予边界节点，这本质上是设置值为0的第二类边界条件。

例如，通道流入口边界和出口边界的宏观边界设置代码如下：

第1、2、6行为直接指定数值，属于第一类边界条件。其他行代码都采用外推方法，本质上是第二类边界条件。

微观量的设置方式比宏观量更复杂。对于入口、出口的边界条件，既可以用Zhou/He边界条件，也可以用非平衡外推方法；对于壁面边界，既可以用反弹边界(Bounce-Back)，也可以用非平衡外推方法。由于非平衡外推方法处理起来简单，适合入口、出口、壁面等方式的微观边界设置，因此应用比较广泛。

非平衡态外推格式是由郭照立等人于2002年提出的一种具有二阶数值精度的微观量外推方法[10]。该方法将概率密度函数的组成分为两部分：一部分为平衡态，另一部分为非平衡态。平衡态部分可以通过宏观边界值直接得到；非平衡态部分则以外推方法得到，计算公式为

式中，f0，f1——　边界格点和内邻格点的概率密度函数；，

——相应的平衡态分布函数。

通道流模型的四个边界上应用非平衡态外推方法的MATLAB代码如下：

其中，第1行表示左入口的非平衡外推格式，移项后即式(5.12)，也就是F(2：ly-1，1，：)-FEQ(2：ly-1，1，：)=(F(2：ly-1，2，：)-FEQ(2：ly-1，2，：))。其他三行代码分别为右出口、下壁面和上壁面三个边界上的分布函数。

非平衡外推方法主要用于处理求解域的边界条件，包括入口、出口和无滑移壁面等。在LBM中，反弹(Bounce-Back)方法也是一种常见的边界处理方法，它主要用于壁面边界的处理。

根据壁面边界的位置不同，反弹边界可以分为边界壁面(on grid)和内部壁面(mid grid)两种。在边界壁面，标准的反弹边界只有一阶精度，但经过改进后可以达到二阶精度(如修正反弹、半步反弹)，内部壁面边界则有二阶精度。在本书中，参与反弹的微观量不参与碰撞操作。若加入反弹边界条件，则需在程序流程中进行3个步骤改动，如图5.6所示。

图5.6　反弹操作在LBM计算流程中的嵌入

首先看看图5.7中所描述在边界上的“on grid”情况。

图5.7　壁面反弹边界中节点上的概率密度函数

对于边界上的节点，f6、f7、f8三个方向在迁移时是空缺的，为了在该节点上形成壁面(无滑移、无穿透)条件，在水平方向应满足无滑移条件，即f1+f2+f8=f4+f5+f6，在竖直方向应满足无穿透条件，即f6+f7+f8=f2+f3+f4，因而最简单的方法就是以f2、f3、f4分别替代f6、f7、f8，这样就使得碰撞后的微观量满足了壁面条件。此外，为了使得壁面边界的平衡态分布函数满足精度，还需要在计算宏观速度后对边界节点上的宏观速度进行修改，使得平衡态速度u eq=0。

对于“mid grid”的情况，根据宏观速度计算公式易知，反弹后微观量完全反向，因此宏观流动与上一时刻的流动趋势相反，并且由于没有增/减新的fi，可保证∑fi守恒，从而使得密度和压力的计算不引入新的误差。

“mid grid”型的Bounce-Back条件用于描述流场内部区域的固体结构，这对于构造特殊形态的通道结构(如狭窄、圆柱、漏筛等)有很大用处。

5.3.8　LBM计算步骤与数学框架

LBM模拟流动的程序框架及对应的数学公式如图5.8所示。

图5.8　LBM模拟流动的程序框架及对应的数学公式

以上是一个相对简单的LBM模拟流动的程序框架。在这个框架基础之上，可以添加外力模型、浸入边界等模型。虽然这只是一个基础框架，但它代表了LBM的基本思路，下一节将以这个框架为依据进行程序实现。

5.4　采用LBM模拟微通道流动

5.4.1　微通道流动

微通道是微流控芯片中最基本的结构。正是作为载体的流体在通道中运动，带动了溶质、颗粒、气泡等在芯片中运行。因此微通道中液体流动的模拟是微流控芯片数值建模的基本要求。在具体设计中，为了提高芯片的面积利用率，通常将通道设置成弯曲的U形通道，如图5.9所示。

流体在压力推送下，从左侧的通道入口流入，从右侧的出口流出。流体在通道中的流动符合泊肃叶流(Poiseuille flow)规律，即管道内不可压缩黏性层流。在二维情况下，泊肃叶流可以描述为一对平行平板中的流动，其特征为垂直于流动方向的速度剖面是以管道中心流速为峰值的抛物线，如图5.10所示。

图5.9　U形通道结构

图5.10　通道流在流动方向的速度剖面分布

根据泊肃叶定律，其速度抛物线分布规律与通道出入口的压力、流体的密度等有关，满足的解析关系如下：

式中，ux ——　通道在纵向位置的水平速度；

P in，P out ——　入口和出口的压力；

Lx ——　通道的长度；

υ——　流体的运动黏度；

ρ0——　流体的密度；

h ——　通道的宽度；

y ——　纵向坐标。注意：坐标原点在通道的中心线上。

以上公式是二维情况下泊肃叶流动的准确描述。通过编程实现的模拟结果需要与之保持严格一致，方能证明模拟结果的正确性。对于LBM初学者，通道流程序通常是必要的练手程序，也是学习后续模型的重要基础。

5.4.2　微通道流动的LBM程序实现

接下来，给出一个采用单松弛因子模型(LBGK)的MATLAB通道流程序。该程序包括一个主程序(Poiseuille.m)和两个子程序(Stream.m和FEQf.m)，子程序分别负责迁移计算和平衡态分布函数的计算。

1.主程序Poiseuille.m

2.子程序Stream.m

3.子程序FEQ f.m

5.4.3　计算结果分析

在τ=0.509，0.56，1.95和2.6四种情况下的数值解和解析解如图5.11所示。

图5.11　通道流速度剖面验证(书后附彩插)

(a)τ=0.509；(b)τ=0.56

图5.11　通道流速度剖面验证(续)(书后附彩插)

(c)τ=0.95；(d)τ=2.6

由于，说明运动黏度υ与松弛因子τ正相关，从图5.11中4幅子图的纵坐标可以看出，随着υ的增大(运动黏度的增加)，在相同的驱动压力差下，流体的稳态流速逐渐减小。

通过调节τ，获得＜10-8时(系统达到稳态时)通道中心线上数值解和解析解的相对误差如图5.12所示。结果表明，当τ＜0.536和τ＞2.9时的相对误差超过了1%，这表明在采用单松驰因子LBM(LBGK)时，τ的范围最好在0.54～2.6之间选择。此外，如图5.13所示，图中统计了系统达到稳态所需的迭代次数。相应数据说明随着τ的增大，系统达到稳态的速度就越快。

图5.12　误差随τ值的改变规律

图5.13　系统达到稳态所需的迭代次数

这里，我们设置的网格数量是60×31，那么，锁定长宽比之后，就可以设置为120×62吗?或者设置为180×93?这当然是可以的，LBM中采用的是无量纲计算，理论上只要离散精度足够，那么模拟结果都是一样的，这称为网格无关性，但网格规模更大时，会导致更大的计算量。然而，如果网格数量太小，虽然计算会很快，但会导致不可接受的数值误差或出现数值稳定性的问题。因此，在很多情况下，为了证明采用的网格规模大小适宜，需要进行“网格无关性验证”，也就是表明模拟结果与网格的划分方式基本无关。第4章的图4.6所描述的结果就可以证明所采用的网格划分方式均满足模拟要求。在满足精度要求的情况下，当然是采用规模较小的网格划分方案能使计算更“经济”。

5.4.4　量纲转化问题

在LBM建模过程中，采用无量纲参数比较方便，但在进行具体的物理量分析时，通常需要采用更加直观的物理量纲，因此量纲转化是建模中的一个基础步骤。由于LBM的量纲转化涉及一些公式描述的对应关系，这些对应关系通常都是一些比值关系，不便理解，因此量纲转化一直是一个经常困扰初学者的问题。接下来，举例说明如何进行量纲转化。

5.4.2节的程序所模拟的通道流，其通道长度为59 μm，宽度为30 μm，对于采用的60×31网格规模，一个格子的宽度就是1 μm，这就是空间尺度的对应关系。假设流动的流体为常温水，则其密度为1 000 kg/m3，其运动黏度υ=1×10-6 m2/s，并假设通道两端压力差为4.8 Pa，那么根据泊肃叶定理(式(5.13))，相应的最大速度u′max为

在LBM下，设υ=1/6，则相应的无量纲速度计算公式为

物理量纲与LBM无量纲之间的对应关系如表5.3所示。

表5.3　物理量纲与LBM无量纲之间的对应关系

续表

基于以上转换关系，从LBM中得到速度、压力差、受力等无量纲结果后，需要通过量纲转化，方可得到相应物理量纲尺度下的数值。

5.5　采用LBM模拟不规则通道中的流动

5.5.1　不规则通道流动模型

在微流控芯片中，直通道是最简单的通道结构。根据设计需要，有时还要在通道中构造特殊的结构，如狭窄结构、柱状阻挡结构、漏筛阻挡结构等。这些结构的共同特征就是在通道的指定位置设置固体区域，以改变流动方向。在LBM中，通常采用反弹(Bounce-Back)方法来完成通道内特殊结构的构建。

在5.4节，我们举了一个采用LBM模拟微通道流动的例子，对这种方法进行了大致介绍。本节在涵盖5.4节内容的基础上，将相对详细地剖析一个LBM程序，以便大家可以更加深入地理解LBM，并掌握LBM模拟不规则通道流动的方法。

本节的建模任务：在通道中建立两个固体圆柱，并通过改变雷诺数来观察流动状态的改变。通道结构如图5.14所示。图中，黑色区域为固体区域，包括上下壁面和通道内部的两个圆柱。观察两个圆柱区域可知，其边缘有锯齿，并不平滑。这是因为，在构造圆柱时，我们只能保证在一个圆形范围内的网格节点为“固体”节点，这些节点位于背景网格的交叉点上，在空间是均匀和离散分布的。图5.15所示为圆柱半径R=2，6，10和15的情况，对比可知，半径越大，圆柱边缘就越贴近正圆形。

本节模型与于5.4节通道流模型的另一个区别在于宏观边界设置(微观边界均采用非平衡外推方法)。本节模型中通道的入口边界条件为速度边界(5.4节中设置的是压力边界)；通道的出口边界为自由出口边界(5.4节中设置的是压力边界)。

图5.14　圆柱绕流模型的通道结构

图5.15　不同半径约束下的圆柱构型

(a)R=2；(b)R=6；(c)R=10；(d)R=15；

5.5.2　雷诺数与流动状态

在这个模型中，流动的状态主要和雷诺数Re有关。1883年，英国物理学家雷诺观察了圆管内的流动状态，首先提出圆管中的流态由层流向湍流的过渡取决于以下比值：

式中，L ——　特征尺度，指代体现影响流动状态的关键尺寸(单位为m)，如通道流中的通道宽度、圆柱绕流模型中的圆柱尺寸等；

U—— 特征速度(单位为m/s)，如通道入口的最大速度；

υ——　流体的运动黏度(单位为m2/s)，，μ是流体 的动力黏 度(单位为Pa·s)。

关于通道流的流态，主要分为以下三种情况：

(1)当流速很小时，流体分层流动，互不混合，称为层流，又称稳流或片流。

(2)逐渐增加流速，流体的流线开始出现波浪状的摆动，摆动的频率及振幅随流速的增加而增加，这种流动称为过渡流。

(3)当流速增加到很大时，流线不再清楚可辨，流场中有许多小漩涡，层流被破坏，相邻流层间不但有相对滑动，还有混合。这时的流体作不规则运动，有垂直于流管轴线方向的分速度产生，这种流动称为湍流，又称乱流、扰流或紊流。

具体而言，当Re＜2 300时，管道中的流动总是层流；Re＞4 000时，流动一般为湍流；其间为过渡区，流动既可能是层流，也可能是湍流，取决于外界条件。当有圆柱在通道中时，情况则有所不同。

对于圆柱绕流，流动的模式就更复杂一些。在不同雷诺数下，对应的流线示意如图5.16所示。

图5.16　不同雷诺数下的圆柱绕流的流线示意图

(a)无涡尾流(Re≤1，大致范围)；(b)对涡尾流(4＜Re＜40)；(c)规则脱涡尾流(60＜Re＜300)；(d)不规则脱涡尾流(Re＞300)；(e)湍流尾流(Re＞105)

如图5.16(a)所示，在低雷诺数情况下(Re≤1，大致范围)，流场中的惯性力与黏性力相比居次要地位，圆柱上下游的流线前后对称，阻力系数近似与Re成反比。

如图5.16(b)所示，随着Re增大，圆柱上下游的流线逐渐失去对称性。当Re＞4时，沿圆柱表面流动的流体在到达圆柱顶点(90°)附近就离开了壁面，分离后的流体在圆柱下游形成一对固定不动的对称漩涡(附着涡)，涡内流体自成封闭回路而成为“死水区”。随着Re增大，死水区逐渐拉长，圆柱前后流场的非对称性逐渐明显。

如图5.16(c)所示，Re增大到40以后，附着涡瓦解，圆柱下游流场不再是定常的，圆柱后缘上下两侧有涡周期性地轮流脱落，形成规则排列的涡阵，这种涡阵称为卡门涡街。此时的阻力系数在1～2之间。

如图5.16(d)所示，Re增大到300以后，圆柱后的“涡街”逐渐失去规则性和周期性，但分离点(约82°)前圆柱壁面附近仍为层流边界层，分离点后为层流尾流。

如图5.16(e)所示，随着Re继续大幅增大，如当Re＞105时，层流边界层随时有可能转涙为湍流，分离点后移至100°以后，湍流时绕流尾迹宽度减小，阻力系数发生骤减。

5.5.3　运行结果与分析

对于本节的模型，考虑Re=100时的流动，流场的涡量图如图5.17所示。

图5.17　涡量图(书后附彩插)

涡量是描述漩涡运动最重要的物理量之一，其定义为流体速度矢量的旋度，涡量的单位是秒分之一(s-1)。涡旋通常用涡量来度量其强度和方向。在流体中，只要有“涡量源”，就会产生尺度大小不一的涡旋。二维情况下的涡量计算公式为

相应的MATLAB代码如下：

上述代码多次用到的“'”是计算过程中需要用到的“矩阵转置”计算，diff()是求一阶差分函数。

二维情况下的速率场如图5.18所示，其中入口平均流速U in=0.1。

图5.18　速率场(书后附彩插)

速率计算公式为

相应的MATLAB代码如下：

注意到以上流场中最大的无量纲速率超过了0.2，虽然程序能稳定运行，但是实际分析时，需要考虑到当流动超过了0.3Ma后，流动会体现出可压缩性。因此，如果考虑不可压缩流动，那么最大的无量纲速率一般不超过0.17，即0.3×c s≈0.17。在这种情况下，为了控制流场最大速度不超过0.17，在设置入口流速时，要调低当前U in=0.1的设置，以U in=0.05为宜。

5.5.4　代码分析

参考5.4节的代码，这里重点对有变动的代码进行注释说明。

第一部分，参数初始化：

这部分代码主要设置通道的尺寸，定义雷诺数Re，根据雷诺数和特征速度，以10为特征长度(小圆柱直径)，可以计算得到运动黏度nu，进一步获得松弛因子τ的倒数omega。第7～10行定义格子参数，包括声速cs、方向权重t，方向向量ex和ey。

第二部分，变量初始化：

第15行，meshgrid()是MATLAB中用于生成网格坐标序号的函数，函数格式为[X，Y]=meshgrid(xgv，ygv)。meshgrid()函数生成的是X、Y大小相等的矩阵，xgv、ygv是两个行向量。X的值通过将xgv复制length(ygv)行(严格意义上是length(ygv)-1行)得到；Y的值通过先对ygv进行转置得到ygv'，再将ygv'复制(length(xgv)-1)次得到。

第17行，对数组中所有元素进行并行操作。代码“bnds=bnds+(((x'-Cen_r(i，1)).^2+(y'-Cen_r(i，2)).^2)＜=Cen_r(i，3)^2)；”是 一 个“比较”操作，对于满足条件的元素返回结果1，若不满足条件则返回结果0。举一个例子，令“A=[1 2 4 5 7 8]”，那么“M=(A＞=5)”返回的结果即“M=[0 0 0 1 1 1]”。本行代码的内涵是计算所有点距圆心的距离，将距离小于等于半径的点标识为1，否则标识为0。由于本例中设置了两个圆柱(两个圆柱无公共点)，因此可以用累加方式“bnds=bnds+…”获得两个圆柱所有节点标识的结果。

第19行，直接将通道的上下壁面处节点均标识为1，这样就设置了固体壁面。

第20行，把二维向量转换为一维向量，并获得固体节点在一维向量中的位置序号。举个例子，令

若要返回4的整倍数的元素位置给变量b，那么首先需获得标识矩阵，其MATLAB代码为“b=mod(a，4)==0”，即

之后，采用“[bx，by]=find(b)”获得值为1的元素位置坐标，可得bx=2，3，4，by=1，2，3。其中bx为元素1的行号、by对相应列号。但如果返回值只有一个，即采用代码“obs=find(b)”，则返回的结果obs为2、8、14。这个结果如何来的呢?看看图5.19中箭头指示方向，沿着指示的方向就可以发现，在第2、8、14个位置的值为1。因此在MATLAB中，采用的是“列优先”的方式进行二维数组元素的排序。

同时应注意到，此处的obs和bx、by之间有换算关系“obs=(by-1)*5+bx”。有了以上的基础，就可以更好地理解第23～30行中关于m和mb的操作。

第21行，给定反弹序列，导致在反弹操作前后，节点的概率密度函数方向发生如图5.20所示的变化。

图5.19　二维数组一维化的列优先原则示意

图5.20　反弹操作前后概率密度函数的分布

(a)反弹前；(b)反弹后

因此，对于流场内部的固体节点，反弹是在其节点上给了一个阻止流动的趋势，从而达到模拟无滑移、无穿透效果的目的。需要注意的是，对于本书使用的LBM程序架构，用于反弹的概率密度函数无须经过碰撞操作，因此在进行碰撞操作之前，应先存储需要反弹的概率密度函数。在所有节点进行碰撞操作之后，替换相应固体节点上的概率密度函数值。

接下来，分析第23～30行代码。这段代码在整体上是一个双重循环，外面的一层循环变量为i，内部一层的循环变量为j。所以，位于内层循环的表达式“n=n+1”返回的是两层循环的累计循环次数结果。

首先来看内层的循环的“m(n)=obs(j)+(i-1)*lx*ly；”。注意：obs是从lx×ly的数组中获得的元素序号序列，因此obs中的最大值不超过lx×ly。那么，在后面为何要加“(i-1)*lx*ly”呢?其原因是概率密度函数F()是一个三维数组，从第24行代码“F(：，：，i)=rho*t(i)”可以看出，F()实际上包括9个维度为lx×ly的二维数组。如果将每个二维数组当成一个二维平面，那么F()的结构就可以表示成如图5.21所示的形式。

【思考】已知点p在F(：，：，5)中的obs=1088(px=180，py=8)，求其在全局中的值m?

由于这里lx=300，ly=60，因此答案为：m=1 088+(5-1)×300×60=19 088。得到这个m值后，在MATLAB中通过在代码中用“F(m)”就可直接定位到F(180，8，5)指代的位置。因此，采用一个数值，即可定位到F()的三维空间指定坐标。而在“mb(n)=obs(j)+(oppo(i)-1)*lx*ly”中，mb存储了反弹方向的位置，因此，令F(m)=F0(mb)，就可以精准地将F0中存储的反弹方向的概率密度函数赋予F中对应位置，以实现反弹操作(参见第三部分的第56行代码)。这种做法的优点是：计算过程不涉及无须参与反弹的流体节点，因此可以避免冗余计算。

图5.21　概率密度函数数组结构图

第三部分，LBM主循环体：

第三部分是程序的核心部分，这里大部分代码在介绍LBM方法时就已经介绍了，此处不再赘述。本部分程序与5.4.2节的主程序Poiseuille.m比较，有两个主要变化。变化一，边界条件设置。首先，入口被设为速度边界，出口被设为自由出口边界；其次，处理上下壁面的微观边界条件时，直接用反弹边界代替非平衡外推方法，因此所用的代码只有两行(第59、60行)，而在Poiseuille.m中所用的代码有4行(第52～55行)。变化二，增加了反弹条件设置，用于构造固体壁面或固体区域，主要包括三处：①第41、42行，对固体节点的速度进行矫正；②第54行，在碰撞之前备份反弹量；③第56行，执行反弹操作。

第四部分，计算结果可视化：

其中，第64～66行显示速率场，第67～70行显示涡量场。若需显示涡量图，则取消语句前的注释符并注释屏蔽第64～66行代码即可。

对比5.4节和5.5节的程序内容，可以发现大部分代码重复，为了获得所需的通道结构，只需对基础的通道流程序进行局部修改。至于是否改得正确，则通常需要进行定量验证，只有定量验证通过了，才可确认代码的修改是正确的。比如，在5.4节中，将通道流速度的数值解与解析解进行对比(图5.11)，就是一种定量化的验证方法。

5.6　应用IB-LBM模拟不规则通道中的细胞运动

5.6.1　流体力学与计算流体力学概念

流体是有固定质量而无固定形状的物体，或承受(任意大小)剪应力作用后会产生连续不断的变形的物体(如气体和液体)。根据受力后变形规律的差异，流体可以分为牛顿流体和非牛顿流体。牛顿流体是剪切力与剪切应变率之间是线性关系的流体，如空气、水等。非牛顿流体是剪应力与剪切应变率之间不是线性关系的流体，如血液、泥浆等。研究流体(液体和气体)的力学运动规律及其应用的学科称为流体力学。20世纪50年代以来，随着计算机的发展，产生了一个介于数学、流体力学和计算机之间的交叉学科，其主要研究内容是通过计算机和数值方法来求解流体力学的控制方程，对流体力学问题进行模拟和分析，该学科称为计算流体力学。

5.6.2　流固耦合

流固耦合是研究变形固体在流场作用下的各种运动规律以及固体位形对流场影响的一门科学，其重要特征是流体-固体两相介质之间的相互作用，变形固体在流体载荷作用下会产生变形或运动；变形或运动又反过来影响流体运动，从而改变流体载荷的分布和大小。由于流固耦合问题的交叉性，其在学科方面涉及流体力学、固体力学、界面动力学等学科知识，在技术方面与不同工程领域(如土木、航空航天、船舶、动力、海洋、石化、机械、微生物工程等)均有关系。其研究问题涉及面甚广，因此在许多领域中都有非常重要的应用。

在本章中，流固耦合主要指流体和具有弹性的细胞界面之间的相互作用。这里的细胞，是一个“胶囊”状的弹性结构，其胶囊壳由“……节点-弹簧-节点-弹簧……”的结构交替形成，节点之间的距离依靠弹簧来保持，相邻弹簧之间的角度相对固定。因此，在细胞受到流体挤压时会发生形变；当挤压力取消后，细胞会回复其初始形状。这里的“弹性结构”由浸入边界(Immersed Boundary，IB)法来模拟实现，流体则由LBM进行计算。通过构造IB和LBM之间的参数传递机制(如力、位移)，可实现流固耦合计算，这种流固耦合计算框架称为IB-LBM流固耦合框架。

5.6.3　欧拉网格与拉格朗日网格

要更好地理解流固耦合的实现机理，须先理解两种类型的网格。一种是网格节点不运动，但节点上的参数发生变化，这种类型的网格称为欧拉网格；另一种是网格节点会在流动影响下发生运动，这种网格称为拉格朗日网格。

欧拉网格侧重于表述“场”，把流体的性质定义为空间位置与时间的函数。通俗而言，就是把空间分成一个个屋子，屋子的位置不变，流体可以自由进出这些屋子；在计算时，给每个屋子一个“门牌号”；在演化后，我们关注的是每个门牌号的屋子里到底在发生什么：想知道某个屋子里的情况，按门牌号对应“敲门”即可。

拉格朗日网格侧重于表述“质点”(又称“流体微元”)，把流体的性质按质点/流体微元来逐个定义；定义的方法通常是把这些性质写成初始坐标的函数，即用质点的坐标改变来描述其运动。通俗而言，就是把流体划分成一个个包裹，包裹的位置一直在移动，大小也会变化，但流体不会穿透包裹袋。首先，给包裹们编号；在演化后，按包裹的编号找包裹，再拆开包裹看内容。

这两种类型的网格有各自不同的适应场景。例如，欧拉网格适合描述一个区域内固定点上的参数改变，如流动中的密度、速度、温度、熵、焓，甚至单位流体中的磁通量，等等；拉格朗日网格适合描述雾霾粒子、河流中的泥沙、水中的气泡、风洞中的示迹的烟，乃至风中的树叶、大洋中的鱼群等的运动规律。这两套网格所表述的方法在一定条件下是可以互相转换的。

综上，结合我们关注的问题，LBM采用欧拉网格，IB则采用拉格朗日网格。

5.6.4　浸入边界法

Peskin在1972年提出浸入边界(Immersed Boundary，IB)法，并将IB应用于血液在可收缩的心脏瓣膜中流动的数值模拟。IB的产生为解决复杂边界情况下的流固耦合提供了新的思路和方法。在流体运动方程中，IB添加了体积力项来描述物体边界和流场之间的相互作用，并在数值计算中使用两套网格体系。一方面，使用固定不变的欧拉网格计算流场；另一方面，使用一系列拉格朗日点来表示物体的边界。然后通过两套网格的信息交换来演化固体边界和流场之间的相互作用。当前，IB在模拟血液流动、湍流的直接模拟以及多相流等方面取得了巨大的成功，在流体力学领域有着很好的发展前景。

图5.22　浸入边界及节点之间的关系[11]

X—浸入边界点的位置；x—流场节点的位置；h—欧拉点的节点间距

接下来，详细介绍IB的基本原理及程序实现。浸入边界与流场节点之间的关系如图5.22所示。图中的空心点称为欧拉点(在LBM中的网格节点)，该点本身不发生移动，只描述某个固定的位置，在这些位置上的流体速度和压力由Navier-Stokes方程描述；图中的实心点称为拉格朗日点，该点可在流体作用下发生移动，并且相邻点之间存在内力(拉伸、压缩、弯曲、剪切等)约束，其运动规律由流体和各相邻点之间的内力共同控制，由这些点组成的边界即浸入边界。

浸入边界法的核心公式为

式中，δ(·)——Dirac delta函 数，简 称Delta函数。

Delta函数的形式不止一种，本质上是以一种以空间距离为依据的加权函数。一种常用的格式为

式中，r——浸入边界点与欧拉点在坐标轴方向的距离。

浸入边界法有两个核心步骤。

步骤1，将浸入边界上的内力(体积力)分配到周围，即实施公式f(x，t)=∫F(s，t)δ(x-X)d s的描述，如图5.23所示。

图5.23　拉格朗日点上的力分配至附近的欧拉点示意(书后附彩插)

对图5.23左图中间的灰黑圆点A，假如其位置坐标是(5.62，7.55)，该点上x方向的内力分量是0.002、y方向的内力分量是0.001。要将这两个力分量体现在对流场的影响上，就需要分配和转移其到流场节点x(欧拉点)上。图5.23右图是转移结果示意。其转移规则用Delta函数控制。式(5.22)中的r是节点X与周围欧拉点x在x方向或y方向的距离，的最大值是2；若超过2，则δ(r)=0。因此，点A实际上只与周围的16个欧拉点存在关系，即图5.23中红色虚线所圈定的范围。

浸入边界上的节点距离也是受到限制的，一般情况下，相邻点距离d s在0.4～0.8个网格距离。如果同时考虑浸入边界上所有点上内力对流场节点的影响，那么分配到流体上的力则分布在浸入边界两侧，如图5.24所示。图中，红色的“*”代表浸入边界上的点，蓝色的点则是浸入边界节点上的力所影响的周围流体节点范围。

下面列出IB步骤1的核心代码，其功能是将浸入边界节点的内力转移分配到流场节点上。

图5.24　浸入边界(拉格朗日点)与附近流体节点(欧拉点)之间的关系(书后附彩插)

其中，第2、3行中的floor()函数表示向下取整，如floor(2.3)=2、floor(5.89)=5。因此，对于本例点A的坐标为(5.62，7.55)的情况，“ii=floor(IBx(1，i))-1：floor(IBx(1，i))+2”返回的结果为4、5、6、7，而“jj=floor(IBy(1，i))-1：floor(IBy(1，i))+2”返回的结果为6、7、8、9。这分别是点A周围16个欧拉点的x坐标和y坐标。第6～9行，计算Delta函数，也就是式(5.22)的代码形式，这里有δ(xi)=1。第10、11行，将分配得到的拉格朗日力累加到欧拉点上，是公式f(x，t)=∫F(s，t)δ(x-X)d s的程序表述。

步骤2，将点A附近16个欧拉点上的速度，采用δ(·)函数加权后，累加到点A上，从而使得点A也具有一个运动速度——式(5.21)中的U(s，t)。这一过程如图5.25所示。

图5.25　欧拉点的速度加权合成拉格朗日点的速度示意(书后附彩插)

从图5.25中可以看到，速度的传递方向正好和步骤1中浸入边界点内力的传递方向相反。此外，如果将这个运动速度与时间步长相乘，就可以更新点A的位置，对所有浸入边界节点进行位置更新，就可以实现运动边界的模拟，即公式=U(s，t)的表述。

对照步骤1的代码，步骤2的MATLAB代码如下：

容易看出，该代码与步骤1代码的差异在于第10、11行，这里的10、11行是公式U(s，t)=∫u(x，t)δ(x-X)d x的代码表述。

5.6.5　细胞力学模型

血液中的红细胞是一种高弹性体，在模拟中常用“胶囊”模型来描述其结构和力学特性。“节点-弹簧”结构是用来实现胶囊结构模型的基本方法，其中，细胞膜界面力学结构可简化描述为多个等距节点与弹簧的交替连接。采用“节点-弹簧”模型构造的细胞模型如图5.26所示。

图5.26　采用“节点-弹簧”模型构造的细胞模型

在浸入边界法的步骤1代码中，第10、11行涉及Fx、Fy两个力分量，这两个力分量来自细胞膜发生变形时，节点上产生的使膜复原的内力。该内力是一个合力，由切向抗拉伸压缩、法向抗剪切和法向抗弯曲三个力分量共同作用，可表示为

式中，s ——　细胞边界的拉格朗日节点坐标；

t ——　时间；

F s(s，t)——　抗拉压力；

F sh(s，t)——　抗剪切力；

Fb(s，t)——　抗弯曲力。

Fs(s，t)控制细胞膜两个相邻节点间的距离，依据新胡克定理，Fs(s，t)表示为式中，El ——　拉压模量；

ε——　伸缩率，ε=ds/ds0，ds0是相邻节点的平衡态距离，ds是细胞在运动过程中相应相邻节点之间的动态距离。

相应的MATLAB代码如下：

其中，circshift()函数为移位函数，可以在MATLAB命令窗口中通过“help circshift”了解其用法和格式。

抗剪力切Fsh、抗弯曲力Fb与细胞膜节点的曲率κ相关，在获得节点的坐标(x，y)后，可以通过下式求得节点的曲率：在得到细胞膜节点的坐标分量cellx和celly后，通过Calcur()函数可以计算各节点的实时曲率。相应的MATLAB代码如下：

抗剪切力F sh是代表节点在法向的抗剪切力学作用，公式为[11]

式中，E s ——　细胞膜结构的弯曲模量；

κ0——　细胞没有发生形变时的节点曲率；

n ——　节点的单位法向量。

相关的MATLAB代码如下：

抗弯曲力F b(s，t)用于表示细胞膜结构抵抗弯曲效应的力学作用，其具体公式为[11]

续接F sh计算所用过程量，相应的MATLAB代码如下：

最终，将这几个力的分量相加，就得到了节点上的合力F。相应的MATLAB代码如下：

由于本书中只涉及二维建模的情况，因此速度、受力只有x、y两个方向的分量。在三维模型中，会有x、y、z三个方向的分量。

5.6.6　IB-LBM流固耦合框架

前面的章节对于LBM和IB两种方法已进行了介绍。本节重点介绍这两种方法是如何结合的。由于前面介绍的LBM并不能处理外力对流场的作用，因此要想实现LBM和IB两种方法的结合，还需要对原有框架进行改进。

LBM中有多种考虑外力作用的方法。在此，选择郭照立在2002年提出的LBM外力作用格式，该外力格式易于理解、稳定性好，是目前常用的方法[10]。在前面介绍的LBM模型基础上，主要有两处改动：

第1处：碰撞步骤变为fi(x，t+Δt)=fi(x，t)-Δt·G i，比式(5.10)增加了Δt·G i项，通常设置Δt=1。其中，

式中，g i ——　流体节点上受到的外力。第2处：计算速度的公式变为

与式(5.4)相 比，式(5.29)在 分 式 的 分 子 上 多 了 一 项0.5g iΔt。式(5.28)、式(5.29)的程序转化将在5.7.3节中进行。

通过分析可知，采用浸入边界法的步骤1计算得到流体节点上的外力，就是为了提供上述公式中的g i。因此，LBM和IB的结合方式可概括为两方面：一方面，浸入边界将其节点的拉格朗日力转移分配到流体节点上，以影响流体的流动状态；另一方面，在浸入边界附近的流场速度加权累计到相应的节点上，赋予节点一个速度，从而使浸入边界具备运动属性。

5.7　一种基于IB-LBM的细胞微流控分离方案设计

5.7.1　模型介绍

5.2～5.6节介绍了关于IB-LBM的数学模型和程序原理，本节将采用IB-LBM编写程序，建立一个完整的细胞操控方案。这里主要利用层流分流的原理，将体积不一样的细胞进行分离。图5.27描述了所建立的通道结构。分析可知，这是由两个相同的操作单元的串联，因此建模时可以只考虑一个单元。在设置时，将其流动边界设置为循环边界(循环边界指流入和流出边界共享同一边界)，以便减少计算量。流体采用力场驱动，这样出口和入口的密度能保持连续性，方便细胞在出口与入口处跨界。

图5.27　细胞分离通道结构与分离效果图(书后附彩插)

通道中的白色部分为栅栏区域和椭圆阻挡区域，这些区域为固体结构，可以采用反弹条件建立。本通道结构设计的原理如下：在栅栏结构中施加一个侧向的力场，以驱动栅栏内的流体从下往上流动，这样流经栅栏结构的细胞将由于分流原理而最终贴向上端栅栏结构的底缘，从而达到对齐的目的。在尺寸不同的细胞沿着栅栏底缘对齐后，细胞中心则位于不同的流层，这些流层在椭圆结构处发生分离，从而将细胞分离。

这个例子的代码，除了用到前面介绍的知识外，还要用到一些其他知识，如循环边界的概念、采用cell型数据来设置细胞存储空间、细胞流动过程的跨界操作等。请读者通过分析5.7.3节中的程序来自学理解。

5.7.2　程序框架

本实例的程序由1个主程序和7个子程序组成：

(1)主程序Main.m。

(2)子程序Ibf.m：浸入边界将力从拉格朗日节点转移分配到欧拉节点。

(3)子程序Ibu.m：浸入边界将速度从欧拉节点加权累计到拉格朗日节点。

(4)子程序Calcur.m：计算节点曲率。

(5)子程序M force.m：计算细胞的弹性力。

(6)子程序FigureOut.：在计算过程中输出图像。

(7)子程序Stream.m：迁移步骤(程序同5.4.2节的Stream.m)。

(8)子程序FEQf.m：计算平衡态分布函数(程序同5.4.2节的FEQf.m)。

程序流程图如图5.28所示。

图5.28　程序流程图

5.7.3　程序代码

1.主程序M ain.m

2.子程序Ib f.m

3.子程序Ibu.m

4.子程序Calcur.m

5.子程序M force.m

6.子程序FigureOut.m

以上是本模型全部代码，其中Stream.m和FEQf.m的代码与5.4.2节的代码相同，因此不再列出。

5.7.4　结果分析

在运行过程中，列出了4个时间点的细胞运动图像，如图5.29所示。这些图像是由FigureOut.m自动生成的图像序列，其中T为程序迭代次数。

图5.29　不同时刻的细胞分布情况(书后附彩插)

(a)T=1 500；(b)T=15 000；(c)T=40 000；(d)T=55 000

图5.29(a)所示为初始化后的情况，大小细胞混合，此时细胞还没有开始分离。如图5.29(b)～(d)所示，可观察到各细胞倾向于贴着上栅栏结构底缘运动，这说明栅栏内通道侧向流动发生了作用，对主通道流动进行了分流。如图5.29(c)(d)所示，可见小细胞倾向流向椭圆区域的上侧通道，而大细胞倾向于流向椭圆区域的下侧通道，这样就实现了不同尺寸细胞的分离。

本模型作为一种不同尺寸细胞分离的概念模型，其分离结果与主、侧通道的流动速度、椭圆区域在通道中的垂直距离以及细胞混合物的疏密程度有关，也和串联的分离单元数量有一定关系。一方面，这说明本模型的工作机理复杂，影响因素多，很难用解析方法来计算分离效果；另一方面，这说明建模分析对于开发此类模型的重要性——如果采用实验方法来进行方案设计和参数调节，不论是时间成本还是物资成本都将很高。

5.8　应用IB-LBM进行细胞分选的案例

在微流控芯片中，颗粒、细胞操控是常见操作。由于有限元法、有限体积法等方法很难处理复杂运动边界，所以直接采用流固耦合来研究弹性颗粒运动的报道并不多。因此，IB-LBM模拟流动中颗粒运动的优势就显得尤为重要和突出。由于IB-LBM是近年来出现的一种解决复杂运动边界流固耦合问题的新兴方法，且这种方法有较好的应用价值，因此本书以大篇幅在第5章介绍其数学原理、学习方法、程序实现等。这里要强调的是，LBM不仅是一种方法，而且是一套强大的数值模拟体系，IB-LBM只是其中的一个应用领域，基于LBM可以研究多相流、多孔介质、传热、对流扩散、沸腾、激波传递、非平衡态过程等复杂物理现象；同时，LBM具有比较开放的框架，不但可以和浸入边界法结合模拟流固耦合，还可以和phase-filed方法结合模拟晶体生长；而模拟油水界面分离和气泡模型的双分布LBM模型可以看成LBM与自身相结合。由于LBM变化多样，加之成熟时间并不长(或者说尚不够成熟)，所以目前还没有很好的软件来支持这种方法的应用，直接编程是学习和应用这种方法的基本要求。

作为IB-LBM方法在细胞操控中的应用，笔者课题组也做了一些研究工作，并获得了一定的成果，接下来进行简单介绍。

5.8.1　错流法分离不同体积的细胞

Huang等[12]提出一种称为Deterministic Lateral Disp lacement(DLD)的新型分离策略，在交错排列的圆柱矩阵中，通过调整微柱之间的间隙和相邻两行的微柱偏移量来达到将不同尺寸的微粒分离的目的。对于一组阵列排布的微圆柱，其微柱的间隙和偏移量可决定一个临界尺寸，称为临界直径D c。当颗粒的尺寸大于该临界直径时，它将以一种被称为Disp lacement mode的方式在微阵列中运动；当颗粒的尺寸小于该临界直径时，它将以一种被称为Ziazag mode的方式在微阵列中运动，如图5.30所示。混合着不同尺寸颗粒的流体从分离装置的左侧进入，不同尺寸的颗粒将采取不同的方式流经装置内的微柱阵列，它们从右侧流出的位置不同，进而达到分离的目的。

文献[3]基于DLD对连续产生的多个血细胞的分离过程进行了模拟，通过调节细胞的变形能力、圆柱排列的倾斜角度、圆柱之间的孔隙率及细胞释放频率等操作，对细胞的分离效果产生影响。研究结果表明，通过参数调节可得到较好的分离结果，如图5.31所示。

图5.30　DLD方法分离不同尺寸细胞示意(书后附彩插)

图5.31　DLD方法分离细胞结果(书后附彩插)绿色-白细胞；红色-红细胞

以上研究为实现DLD进行细胞的高效分离提供了优化设计方案，也对血液中基于体积差异的癌细胞筛选提供了一定的参考。

5.8.2　夹流法分离不同体积的细胞

PFF(Pinched Flow Fractionation，夹流分离法)自从提出后就以其操作简单、适用范围广等优点引起人们的关注，很多学者用实验、理论研究以及数值模拟的方法对其进行研究。我们也采用IB-LBM开展了相应研究，芯片设计方案如图5.32所示。本模拟结构的入口采用速度边界，出口采用压力出口。芯片的总长度为x0、总宽度为y0，入口Inlet 1和Inlet 2的流量分别为Q in 1和Q in 2，夹流区宽度为w0，出口Outlet 1和Outlet 4通道为直通道，且宽度为w b。为了方便边界条件的设置，我们将出口Outlet 2和Outlet 3设为弯曲通道(分为两段)，前一段通道的宽度为w b，后一段通道变窄，宽度为w e。为了使4个通道的流阻相同，我们在模拟中改变4个通道的长度，使得4个出口压力相同时的流量相同。这样，改变出口Outlet 4的压力时，出口Outlet 1、Outlet 2和Outlet 3的流量相同。当Inlet 1与Inlet 2的流量之比(Q in 1∶Q in 2)大于特定值时，夹流区处的流体会迫使细胞紧贴下壁面流动，以此可以使细胞在夹流区的位置相对固定，因此可以调节Q in 1∶Q in 2的值，使特定大小的细胞紧贴下壁面。芯片的结构尺寸及网格参数设置如下：

结构尺寸：x0=459 μm，y0=398 μm，w0=30 μm，w b=26 μm，w e=23 μm，w i=70.11 μm。

网格 参 数：x0=459Δx，y0=398Δx，w0=30Δx，w b=26Δx，w e=23Δx，w i=70.71Δx，Δx=1。

图5.32　夹流法分离不同尺寸细胞的通道结构设计

这里，改变Outlet 4的压力，就可以改变其流量，并能控制细胞的出口。分离效果如图5.33所示，通过调节Outlet 4的流出流量占比β来分离不同尺寸的细胞。

关于本研究的详细情况，可参阅文献[2]。

5.8.3　U形筛捕获细胞系统优化设计

这项工作中，我们采用IB-LBM方法研究了在微流控装置中实现细胞的捕获，研究方案如图5.34所示。该装置在一个圆形腔室(半径为r4=117.6 μm)中有7个U形筛子，一个筛子放置在腔室中心，其他筛子对称放置于围绕室中心的圆上(半径为r3=76.8 μm)。每个筛子的半径为r2，由6个半径为r1的圆柱组成。设置150个细胞，在入口随机位置一个接一个释放，细胞在经过筛子时，会被筛子捕获。模拟中，设置影响捕获效率的参数有入口宽度h，筛子的半径为r2及圆柱半径为r1。

图5.33　通过调节Outlet 4的流出流量占比β来分离不同尺寸的细胞(书后附彩插)

(a)β=0.1，分离8 μm和20 μm细胞；(b)β=0.2，分离16 μm和24 μm细胞；(c)β=0.4，分离8 μm、16 μm和20 μm细胞；(d)β=0.5，分离8 μm、16 μm和24 μm细胞；(e)β=0.6，分离8 μm、16 μm、20 μm和24 μm细胞

细胞捕获和培养的三个步骤如图5.34(b)所示。其中，b1为装载步骤。为了培养细胞，使其在数量上增殖，每个筛子需要几个细胞作为种子。这是细胞培养的第一步，也是必要的步骤。种子细胞的分布是否均匀，对后续细胞的培养、生长和增殖有显著影响。b2为细胞培养步骤。细胞从周围液体中获得营养以生长并逐渐繁殖成更多细胞，然后新细胞结合在一起或附着在相邻的实体壁上。b3为胰蛋白酶化步骤，在该步骤中，附着在其他表面的细胞可以通过在设备中加入一定量的酶试剂被释放到流体中，呈悬浮状态。最后实施冲洗，其中大部分悬浮细胞被水流冲走和提取。在冲洗过程中，少数细胞可能再被筛子拦截，返回b1，被拦截的细胞成为下一轮培养周期的种子细胞。在这个过程中，我们主要关注b1，即种子细胞的加载步骤。

图5.34　细胞捕获结构设计(书后附彩插)

(a)细胞捕获结构设计；(b)细胞培养流程示意

设圆形腔室的直径为D，在其他条件不改变的情况下，只改变入口、出口的宽度h，我们可以观察到各腔室和U形筛在150个细胞释放完毕后的装载情况，如图5.35所示。显然，更窄的通道入口、出口更有利于细胞的装载。

图5.35　改变通道入口和出口尺寸后的细胞捕获结果对照(书后附彩插)

(a)D/h=1.5；(b)D/h=2；(c)D/h=2.5；(d)D/h=3；(e)D/h=3.5；(f)D/h=4

关于其他参数变化对细胞装载效率的影响，可参见文献[1]。

5.9　学习LBM的策略

总体上，LBM比传统的计算流体力学思路更加简单，在MATLAB下模拟二维通道流动，仅需几十行代码即可实现。但是，LBM不是一种模式固定算法的代称，它针对不同的模拟对象和要求，在参数处理和设置上有很多技巧，此外，LBM还通常需要与其他算法相结合，以解决特定问题。因此，学习LBM需要一定的策略，方能快速准确地熟悉和使用这种方法。根据笔者的经历和感受，给出一套学习LBM的策略，供大家参考。

1.获取可供学习的文献资源

1)介绍算法机理的书籍

通常，介绍算法的书籍能较全面和系统地阐述算法的机理，这有助于我们从算法本源的角度来把握算法的本质。但是，书籍的描述(尤其是数学模型的描述)往往是精炼的、抽象的，很多学习者(尤其是非本专业领域的学习者)发现，仅通过阅读书籍，有很多知识点难以理解，以至于看到后面就越来越看不懂。因此，通过看书的方法来学习算法，通常只能掌握算法的概貌，可定位为“导读”或“背景阅读”。对于学习LBM，华中科技大学郭照立[5]和西安交通大学何雅玲[6]所著的书籍，以及意大利著名教授Succi关于LBM的书籍[7]是比较受推崇的。此外，近年来也出现了一些偏算法程序实现的相关书籍，如A.A.Mohamad所著的书籍[8]，这类书籍不要求读者有厚实的数理理论知识，非数学物理学科专业的研究生更易于理解。

2)学位论文及期刊论文

硕士和博士学位论文展示了研究生对算法运用的具体结果，论文对于算法的学习也很有参考价值。虽然论文中关于算法的描述部分与有些书籍中的算法描述会有重复，但由于需要进行算法的具体实现，所以可能会有一些关于算法实现的细节描述。这些被强调的细节对于算法的实现有时具有很重要的参考价值。在LBM方法中，LBGK(Lattice Bhatnagar-Gross-Krook)模型是最常见和基础的模型，其中的松弛因子τ是非常重要的一个系数。然而，关于这个系数的取值范围，在经典的书籍中往往没有描述，但是在某些学位(期刊)论文中能找到参考。期刊论文在算法描述方面往往比较精简，就算法实现来说，其算法细节不会很翔实。然而，在需要对算法的性能进行探讨或对算法本身进行创新时，一些期刊论文能为改进算法提供指导思路。例如，Luo等[9]对比了LBM中几种常见的模型，对比了各种模型的特定适用范围，这对于深入了解各模型的性能特征有很好的指导作用。

2.获得可供参考的样本程序

学习算法，既可以按“基础理论-数学形式-计算机程序”的顺序进行，也可以按“计算机程序-数学形式-基础理论”的倒序进行。通常，计算机程序是算法的最终体现形式，因此只要能结合对算法的理解读懂程序，也就容易理解一个算法的来龙去脉了。一些新近的算法一般在网络上都有基础的、开源的、可参考的样本程序，这些程序往往是我们高效学习算法的重要指导。在LBM的学习过程中，有些重要的网站不仅提供样本程序，还提供开源的软件包。这些网站提供的程序资源有助于我们快速学习LBM。在palabos网站的W IKI中提供了LBM的圆柱绕流、多相流、热对流等的MATLAB语言样本程序，这些样本程序对于理解LBM算法的碰撞、迁移机理、平衡态分布函数、宏(微)观边界设置方法提供了范例。还有一些网站开发了开源的LBM程序软件包，提供了基础的算法平台，如palabos、OpenLB等，这些平台为大家提供了方便和快捷的程序开发工具，但在使用这些开发平台之前，仍需对LBM的基础机理有充分理解。此外，还有一些专门的程序网站(如程序员、GitHub等平台)上也有很多宝贵的程序资源，涉及与LBM相关的各类诸如LBM-MRT模型、3D的LBM模型等可借鉴的程序代码，可为我们学习LBM提供极大的帮助。从笔者学习LBM的经历来看，一边看文献，一边读程序，两者相互对照可显著提高学习效率。

3.获得“在线”学习支持

在算法的学习过程中，通常会遇到一些难以逾越的“关卡”，这些“关卡”既可能是算法的难点部分或新拓展部分，也可能是数学公式向程序转化过程中的细节处理技巧。当这些“关卡”通过文献(或其他方式)难以克服时，学习者则希望能直接请教已经逾越这些关卡的前辈，以得到指点。如果在学习者身边有可方便询问的前辈，问题的解决可能就会变得简单，但如果身边没有这样的前辈，那么获得“在线”支持就十分重要。这里所说的在线支持主要包括三种方式：

(1)QQ群、微信群。LBM方法涉及广泛的理工科学科领域，随着近年来的应用扩展，出现了一批活跃的QQ群、微信群。这些群里有资深专家、青年教师，还有一些研究生和初学者。对于算法和程序中的关键问题，在群里提问一般都会得到合适的解答。例如，LBM的无量纲转化问题令很多初学者头疼；LBM的状态方程来源于稀薄气体动力学，为什么LBM又可以用于液体的流动模拟；等等。

(2)专业的网络论坛。有时，论坛也是很好的在线帮手，在论坛中，问题如果能得到广泛关注，往往能被讨论得更加充分。此外，我们也可以通过论坛的搜索功能，搜索在学习过程中可能遇到的共性问题，通过了解相关的问题以及回复，学习者可能会发现自己所关注问题的解决方案。在LBM学习中，流体中文网的LBM模块、Palabos的论坛等，有很多可供学习的问题和解答。

(3)E-mail问询。有时，我们在科技论文中发现有关某种新算法的报道，而这种算法很有希望能解决我们面临的科学研究问题，然而由于是创新的算法，在其他渠道找不到可供学习的参考资料。在这种情况下，我们可以直接联系论文的通讯作者，通过E-mail的方式问询算法中的相关问题。在LBM的学习中，笔者曾采用E-mail问询的方式询问过论文作者的算法问题，也解答过来自其他学者的提问。

4.积极参加学术活动，获得现场学习机会

对于算法学习，现场研讨也是非常重要的方式。现场算法研讨主要针对最新的算法模型及应用，这对于科学创新是非常重要的。现场研讨常见的形式包括：

(1)资深专家的专场报告。专场报告通常聚焦在新算法的机理及创新应用方面，不会太专注算法的细节，比较宏观。专场报告有利于开阔我们的研究视野，同时学习者也可以抓住机会提问，结合报告内容和自身需求向专家征求解答，解决自身在学习和研究中遇到的问题。

(2)学术会议。对于主题明确的学术会议，通常聚焦领域内的重要问题和前沿问题，有较多的同行齐聚和报告，建议学习者应尽量参加。在参加学术会议的过程中，一些报告内容可以激发我们的研究灵感，而且还有认识同行的机会，大家通过互相了解，可在一起询问和探讨算法和模型的细节问题。在LBM领域中，比较重要的学术会议有两个：介观方法国际会议ICMMES(International Conference for Mesoscopic Methods in Engineering and Science)、ICDSFD(International Conference on Discrete Simulation of Fluid Dynamics)。

(3)学术论坛。在LBM领域，一年一度的“LBM与微流动学术研讨会”自2013年组织以来，已经举办九届(截至2021年)，笔者是该论坛的主要组织者之一。该论坛以LBM算法及其应用为基础，聚集了国内一大批来自理工科学科领域的中青年学者和研究生，可为他们提供一个学习和使用LBM方法的重要平台。

5.总结、拓展、传承学习成果

算法学习的目标是理解、掌握和运用。对算法的充分了解，应该同时体现在三方面：掌握算法的基本理论；数学形式；程序实现。对于LBM的学习而言，其核心是碰撞和迁移步骤，同时，平衡态分布函数、边界设置等也是重要环节。但是需要说明的是，LBM的算法相对于其他确定性的算法(如C-均值聚类、遗传算法、小波分析、支持向量机等算法)不同的是，LBM本身有着更丰富的内容。除了空间维度不同(包括一维、二维、三维)，还包括不同的格子形式(如D2Q7、D2Q9、D3Q15、D3Q19等)，不同的松弛因子处理形式(如LBGK、MRT)，不同的边界设置(如反弹、外推、非平衡外推、Zhou/He边界等形式)，不同的算法升级模式(如LBM多相流模型、浸入边界-LBM模型、LBM动理学模型等)，以及不同的程序架构(如串行、MPI并行、GPU并行等)。这些内容之间相互组合，使得LBM可作为一类算法。还需要说明的是，我们在设计程序时，程序语言的选择也需要技巧，根据笔者的经验，在初学LBM时，MATLAB是较好的语言选择，因为MATLAB程序语言精练，便于把握算法的步骤和全局结构，输出、调试等都十分方便。但在程序的应用过程中，MATLAB的运行速度是一个大障碍，尤其是三维程序的计算非常耗时。这时选择C++或Fortran就是更好的选择。基于LBM的系统性和复杂性，在学习LBM的过程中，就需要对LBM进行有条理的总结，最好对常用的具体方法亲自编写程序。实际上，由于一些步骤和代码是通用的，当从一种具体算法(如LBGK)改变为另一种算法(如MRT)时，只需要对程序进行局部调整即可。在对感兴趣的LBM算法已经进行程序实现并进行验证后，就可以相对自由地在原有程序基础上进行适当拓展，以满足各自领域的研究需求。同时，我们对算法和程序本身进行有条理的总结，对于一个大课题组来说，这是良好的传承。这种传承，将是整个课题组的资源积累，后面的同学可以在这个积累的基础上少走弯路，提高学习效率，从而有更多的时间去发展和应用原有的算法，做出更高水平的创新。

参考文献

[1]MA J T，ZU W H，TANG X Y，et al.An IB-LBM design of a microfluidicsbased cell capture system[J].ARCHIVE Proceedings of the Institution of Mechanical Engineers Part C Journal of Mechanical Engineering Science 1989-1996(vols 203-210)，2020：095440622094606.

[2]MA J T，XU Y Q，TANG X Y.A numerical simulation of cell separation by simplified asymmetric pinched flow fractionation[J].Computational and Mathematical Methods in Medicine，2016：25645844.

[3]WEI Q，XU Y Q，TANG X Y，et al.An IB-LBM study of continuous cell sorting in deterministic lateral displacement arrays[J].Acta Mechanica Sinica，2016，32(6)：1023-1030.

[4]MA J T，XU Y Q，TIAN F B，et al.IB-LBM study on cell sorting by pinched flow fractionation[J].Bio-Medical Materials and Engineering，2014，24(6)：2547-2554.

[5]郭照立，郑楚光.格子Boltzmann方法的原理与应用[M].北京：科学出版社，2009.

[6]何雅玲，王勇，李庆.格子Boltzmann方法的理论及应用[M].北京：科学出版社，2009.

[7]SUCCI S.The lattice Boltzmann equation for fluid dynamics and beyond[M].New York：Oxford University Press，2001.

[8]MOHAMAD A A.格子玻尔兹曼方法——基础与工程应用[M].杨大勇，译.北京：电子工业出版社，2015.

[9]LUO L S，LIAO W，CHEN X W，et al.Numerics of the lattice Boltzmann method：effects of collision models on the lattice Boltzmann simulations[J].Physical Review E，2011，83：056710.

[10]GUO Z L，ZHENG C G，SHI B C.Nonequilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method[J].Chinese Phys.，2002，11：366-374.

[11]ZHANG J，JOHNSON P C，POPEL A S.An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows[J].Physical Biology，2007，4(4)：285-295.

[12]HUANG L R，COX E C，AUSTIN R H，et al.Continuous particle separation through deterministic lateral displacement[J].Science，2004，304(5673)：987-990.