不同格子Boltzmann模型在毛细管流动模拟中的应用

2019-03-09马瑞程辛显康李宜强管错郭虎

石油钻采工艺 2019年5期

马瑞程辛显康李宜强管错郭虎

1. 中国石油大学(北京)提高采收率研究院；2. 中国石油大学(北京)石油工程学院；3. 中国石油大学(北京)继续教育学院

流体在多孔介质中的流动模拟是油气田开发的重要研究方向。目前，多数流动模拟方法依据连续性方程与达西定律中平均流速和压力梯度的线性关系，在宏观尺度上对求解空间和时间进行离散求解。其中，连续性方程的本质仍然是Navier-Stokes方程的特殊形式之一。实际多孔介质的孔隙空间微结构形状复杂、尺寸不一，流体在多孔介质中的实际流动情况也更加复杂。近年来，从微观层面更精细地研究渗流规律成为一个热门话题。格子Boltzmann方法(LBM)边界处理简单，在边界几何形状复杂的多孔介质中具有明显的优势［1-5］。利用LBM在孔隙尺度层面上研究流体流动过程，如测定多孔介质绝对渗透率等，能够进一步较精确地解决真实复杂的多孔介质流动问题。

1 不同格子Boltzmann模型以及边界条件

1.1 标准不可压缩格子Boltzmann模型

对不可压缩流体的流动模拟，单松弛模型(LBGK)是目前应用最广泛的格子模型。其中，文献［6］提出的DnQb模型(n代表维数，b代表速度离散方向的个数)使用最普遍，常用的模型包括D2Q9、D3Q15和D3Q19等。在这些模型中，平衡态分布函数可以表示为

对于D2Q9模型有

式中，fieq为i方向平衡态分布函数，ωi为权系数，ρ为流体密度，ci为离散速度的方向矢量，cs为格子声速，u为流体流速矢量，c为格子速度。

格子速度数值上等于离散后的空间步长和时间步长的比值。在确定了模型参数和输入参数后，即可进行流动模拟。重复以下3个步骤，即可进行不可压缩Navier-Stokes方程的模拟。

(1)碰撞步骤。

(2)迁移步骤。

(3)计算宏观量。

式中，fin为第n时间步i方向的分布函数；fieq为i方向平衡态分布函数；τ为无量纲松弛时间，其大小和流体的运动黏度有关，通常取值范围为0.5～1.5［7］；t为时间变量；δt为时间步长；x为节点坐标矢量；p为流体压力；υ为流体运动黏度。

1.2 改进的不可压缩格子Boltzmann模型

尽管LBGK模型应用广泛，然而其在模拟不可压缩流体流动时，需要在流动马赫数以及压力变化等参数合适时才能够保证足够的精度，模型的实用性大打折扣。为了克服这一问题，人们基于LBGK模型提出了一系列的改进模型，如Zou-Hou模型、He-Luo模型以及D2G9模型等［7］。重点研究二维模拟和三维模拟均适用的Zou-Hou模型以及He-Luo模型。

1.2.1 Zou-Hou模型

Zou-Hou模型的基本理论［8］：如果将D2Q9模型的平衡态分布函数表达式更改为

同时，压力和运动黏度的计算公式保持不变，而宏观密度和宏观速度的计算表达式定义为

通过Chapman-Enskog展开可以得到宏观微分方程为

可以发现，定常情况下，上述方程组与定常不可压缩N-S方程组一致。

1.2.2 He-Luo模型

该模型的基本思想是利用数学方法消除平衡态分布函数中由密度变化引起的高阶马赫数相关项［9-10］。在此模型中，新的平衡态分布函数为

引入压力函数pi

构建出以压力为平衡态的分布函数为

式中，ρ0为不可压缩流体的真实密度，pi为i方向上的压力分量，pieq为i方向平衡态压力，p0为原始状态下的基准压力。

压力表达式为

因此得到基于压力的演化方程为

基于压力分布函数模型的宏观变量计算公式为

1.3 LBM边界条件

在使用LBM进行流动模拟时，除了设定初始状况的参数，还需要设定边界条件。将多孔介质壁面设定为反弹格式边界，同时将出入口端面设定为恒速边界或定压边界。

(1)反弹格式边界。标准反弹格式是处理静止无滑移壁面的一类常用格式(常用于固液边界)。如图1所示，该格式假设粒子迁移到固体边界后，将会在下一时间步沿迁移的反方向弹回［11］。碰撞后，壁面上的分布函数为

式中，

为边界反弹方向，

为边界反弹后的分布函数，

为边界反弹前的分布函数，

为壁面格点坐标，

为流体所在格点坐标

为壁面反弹后的速度矢量，

为壁面反弹前的速度矢量。

图1 反弹边界格式Fig. 1 Sketch of rebound boundary format

(2)非平衡外推格式。除了用于壁面的反弹格式，流动入口和出口的边界条件仍需要设定。针对管流流动和渗流流动的入口和出口，定流速边界和定压力边界是常用边界。采用非平衡外推格式得到边界表达式［12］。

恒速边界为式中，ub为流体在边界的速度，ρ(xf,t)为t时邻近边界流体域格点密度，ρb为边界节点的密度。

2 LBM模型程序实现及模拟结果

2.1 程序编写及单位转换

研究中LBM模拟器采用MATLAB进行编写，其主要流程如图2所示。

图2 LBM程序执行流程Fig. 2 Execution flow of LBM program

在程序进行初始化时，要将输入的国际单位制物理参数转化为格子单位。lu代表格子长度，ts代表时间步长，mu代表格子质量，依据量纲导出格子系统中其他物理量的基本单位。

不可压缩流体在毛细管轴向剖面上的流动是一个二维问题，因此采用D2Q9模型进行求解。如图3所示，假设毛细管长度为L，半径为R，流动流体黏度为μ，在压差Δp下为黏滞性渗流，渗流入口速度为Uin。流体润湿毛细管壁的情况下，管壁出的液体流速为0，在管中轴线处流速最大。

图3 流动模拟情境Fig. 3 Flow simulation scenario

考虑到黏滞力和驱替压力的平衡关系，得出Poseuille流动的中的流速表达式为

在沿径向的剖面上，流动平均速度可以表达为

式中，r为毛细管中某一点距毛细管中轴线的径向距离，m；为距离毛细管中轴线的径向距离为r处的流速，m/s；pin为入口压力，Pa；pout为出口压力，Pa；R为毛细管半径，m；μ为流体黏度，Pa· s ；L为毛细管长度，m;为毛细管中流体平均流速，m/s。

通过上式可以看出，入口速度剖面应当为一条抛物线，而平均流速和驱替压差呈线性关系。

将输入参数进行转化后即可进行Poseuille流动模拟。其中，模拟需要的松弛时间和驱替压差折算等效密度为

式中，μD为流体黏度，ρD为无量纲流体密度，pDin为入口压力，ρDin为入口等效密度，pDout为出口压力，ρDout为出口等效密度。

2.2 程序模拟结果及误差分析

编写模拟器并运行后，可得到二维Poseuille流动时速度分布图像，从图4可以看出，在求解定常问题时，随着时间步数的增加，速度场的分布逐渐均匀，最终趋于收敛。此外，径向上的速度剖面为一抛物线，且毛细管壁处的速度为0，这与Poseuille解析解的基本假设一致。为了确认模拟结果是否正确，设毛管半径为3 μm，毛细管长度为7.5 μm，流体密度为1 000 kg/m3，流体黏度为1 mPa · s，两端的驱替压差为700 Pa。经过计算，LBM的平均模拟结果为0.034 4 m/s，而解析解计算结果为0.035 m/s，二者相对误差为1.71%。求解收敛时的速度场分布图和毛细管中各个位置的速度剖面情况如图5和图6所示。计算结果表明，LBM在模拟液体在毛细管中的低雷诺数流动上具有足够的准确性，加之LBM在编程上较为简洁且易于并行计算，因此，其在处理微观层面上的多孔介质流动问题具有明显优势。

图4 不同时间步数无量纲速度场分布Fig. 4 Distribution of dimensionless velocity field of different time steps

图5 真实算例无量纲速度场分布Fig. 5 Distribution of dimensionless velocity field of a real case

图6 真实不同位置出口速度剖面Fig. 6 Exit velocity section of different real positions

为了进一步对3种LBM模型的性能进行评价，选取相同的Poseuille流动算例对3种不同模型的计算结果与解析解结果进行比对。在计算过程中，为便于评价和分析，算例中的输入参数全部为已经预先处理好的无量纲参数，模拟输入参数及模拟结果如图7所示，可以看出，随着驱替压差折算后的当量密度不断增加，3种模型计算出的无量纲平均速度尽管也随之增加，但都与解析解结果产生了误差。其中，Zou-Hou模型与解析解契合最好，He-Luo模型次之，LBGK模型最差。

为了对其精确性和稳定性进行分析，作出不同模型在不同马赫数下计算结果的相对误差曲线，结果如图8所示，可以看出，3种模型的相对误差均随着马赫数的增加而增加，误差增幅从小到大排序依次为：Zou-Hou模型＜He-Luo模型＜LBGK模型。

图7 不同模型及解析解计算结果对比Fig. 7 Comparison between the calculation results of different models and the analytical solutions

图8 不同模型相对误差随马赫数变化曲线Fig. 8 Variation of the relative error of different models with the Mach number

在标准LBM模型中，动量方程与不可压缩NS模型不完全一致［13］。如果流动马赫数极小时，则LBGK模型对应的宏观方程与标准的不可压缩Navier-Stokes方程近似。在流体力学中，常用马赫数作为衡量流体压缩程度的无量纲准数，其定义为流体流速和声速的比值。流体力学原理表明，当马赫数大于0.3时，流体的压缩性将不可忽略。而LBGK方法是求解不可压缩N-S方程的一种人工压缩解法，其相对误差会随着马赫数的增加而逐渐增加，因此LBGK模型适用的马赫数范围较小。Zou-Hou模型和He-Luo模型的误差尽管也会随着马赫数的增加而增加，但二者均可以通过数学方法消除与马赫数有关的误差项，相对误差增幅相对较低。

在LBM中，时间步长往往和格子长度满足一定的相关关系［14］。在多孔介质流动模拟中，如果直接对物理问题进行折算，时间步长通常在10-12s到10-10s数量级之间，模拟效率很低。如图9所示，为了提高运行效率，需要保证雷诺数不变，提高马赫数对物理问题进行放缩从而增大时间步长［15］。此时，可以放缩的最大倍数受到模型精确性的制约。3种模型的运算效率从高到低依次为：Zou-Hou模型＞He-Luo模型＞LBGK模型。因此，在编写二维和三维通用模拟器时，如果求解问题为定常问题，推荐使用Zou-Hou模型进行求解；而如果求解问题为非定常问题，则推荐使用He-Luo模型进行求解。