基于格子Boltzmann方法的分叉微通道内非牛顿流体流动特性研究*

郑馨瑶 杨 刚 张会臣

(大连海事大学船舶与海洋工程学院 辽宁大连 116026)

微通道是指水力直径在微米量级的流动通道，具有高的表面体积比、高的传热效率及较少的流体使用量等特点，广泛应用于微流体系统[1-3]。由于尺度的限制，微通道内流体的流动特性很难通过实验完全揭示清楚，并且有些实验现象的分析也亟待深入。对微通道内流体流动进行数值模拟，是对实验中难以观测记录现象和对实验现象加以分析的有效手段[4]。这些数值模拟方法主要包括计算流体动力学、格子玻尔兹曼方法以及分子动力学。相比于宏观尺度下的计算流体动力学方法与微观尺度下的分子动力学方法，格子玻尔兹曼方法是一种介于流体连续性假设与分子动力学之间的方法，节省计算资源的同时，可揭示微通道内流体流动的机制[5]。

早期一些学者借助格子Boltzmann方法模拟微通道内气体的流动。TANG等[6]采用格子玻尔兹曼方法模拟了微通道中的等温气体流动，将反弹边界和镜面反射格式组合模拟滑移边界条件[7]，推导出了松弛时间与努森数的关系。LIU和GUO[8]采用格子玻尔兹曼方法研究了长矩形微通道内压力驱动的气体流动，并与其他方法所得数据进行了对比。此外，由于格子玻尔兹曼方法易于模拟复杂边界，还被广泛应用于液滴动力学特性的模拟。梁宏等人[9]采用格子 Boltzmann 方法模拟了液滴在分叉微通道中的迁移过程,发现当毛细数足够大且壁面具有亲水性时，液滴会发生二次分裂；而毛细数足够小时，液滴在分叉处滞留并且不发生分裂，通过调节出口流量比可以使液滴实现非对称分裂。LIM、KEIVAN等[10-11]采用格子玻尔兹曼方法模拟了微通道中液滴的形成，研究表明液滴的大小以及液滴之间的距离与微通道的几何形状关系密切。

应用格子Boltzmann方法，通过局部计算剪切应力可模拟非牛顿流体的流动。DONG等[12]将流变方程引入格子Boltzmann方法中，建立了非牛顿流体两相流的格子Boltzmann模型，模拟了存在矩形粗糙单元的微通道内幂律流体非牛顿流体的位移，发现幂律指数的降低导致位移速度的增加，相对粗糙度的减小和粗糙单元间距的增大可加大位移速度，从而减少位移时间。谢驰宇等[13]将非牛顿流体模型引入格子Boltzmann方法中，研究多孔复杂结构内多相非牛顿流体的驱替过程，发现了界面张力在驱替过程中起阻碍作用，并分析了压差、润湿性对驱替过程的影响。SHI和TANG[14]采用格子Boltzmann方法模拟黏性指进现象，发现剪切增稠流体中的黏性指进增强，剪切稀化流体中的黏性指进被抑制，重力加速度导致“手指”形状不对称，高毛细数与高黏度比导致“手指”狭窄。该研究从介观的角度解释了黏性指进现象。

尽管通过实验研究与数值模拟，对微通道内流体流动的过程有了一定程度的了解，但对于流动过程中的细节仍不够清晰。本文作者采用格子Boltzmann方法，模拟不同角度的分叉微通道内剪切稀化流体的流动，通过流动过程中密度随时间的变化趋势以及稳态流动下的密度，得到微通道内压力的分布以及流动区间的压力降，着重分析溶液质量分数、入口速度与分叉角度对非牛顿流体流动特性的影响，阐明流体特性和微通道几何构型对非牛顿流体流动行为的影响机制。

1 模型构型与模拟方法

1.1 模拟构型

文中模拟的是高度和宽度均为300 μm的矩形截面多角度分叉微通道，分叉角度分别为60°、75°、90°、105°和120°。图1所示为分叉角度为120°的微通道结构示意图。微通道具有2个出口，非牛顿流体从左侧入口以固定速度流入，流经分叉点处后由2个出口流出。

图1 微通道结构示意Fig 1 Structure of microchannel

选取3种不同质量分数的聚丙烯酰胺(PAM)水溶液进行模拟，探讨流变特性对其流动特性的影响。聚丙烯酰胺水溶液为剪切稀化流体，其流变特性采用幂律模型描述，稠度系数K和幂律指数n2个参数根据流变仪测得的实验数据拟合得到，K和n及其他特性参数如表1所示。

表1 流体的特性参数

1.2 模拟方法

首先采用格子Boltzmann方法求解流体速度场，在引入BGK碰撞算子[15]后，具有二次精度的格子Boltzmann方程离散为

(1)

平衡态分布函数的具体形式一般与离散速度模型的构造有关，因此离散速度模型的构造至关重要。文中采用DnQm模型[16]，D2Q9是最具代表性的格子BGK模型，二维D2Q9模型如图2所示。

图2 D2Q9模型Fig 2 D2Q9 model

对于二维D2Q9，其局部平衡分布函数的表达式为

(2)

式中：wk表示k方向的权重因子；cs为格子声速；ck为沿着格子迁移方向的速度矢量。

流体的宏观密度经平衡分布函数求和得到

(3)

压力由宏观密度得到[17]

(4)

广义牛顿流体中的幂律模型形式简单，广泛使用，且与剪切稀化流体的真实特性基本相符，所以选用幂律模型描述非牛顿流体的动力黏度

(5)

式中：K为稠度系数；n为幂律指数；γ为剪切速率。

采用根据剪切黏度调整松弛时间的直接方法，进行格子玻尔兹曼方法中非牛顿流体特性的模拟[18-20]。首先计算格子的黏度，再根据黏度调整松弛时间。其中，通过式(6)迭代计算形变速率张量，根据式(7)对剪切速率进行计算。

(6)

(7)

最后通过格子玻尔兹曼方法中松弛时间与黏度的关系调整松弛时间，实现剪切稀化特性的模拟。

(8)

在格子玻尔兹曼方法的模拟中，一般采用格子单位，与实际有关的问题模拟时需要进行物理单位与格子单位的转化。对于不可压缩流体流动问题，需要匹配宏观雷诺数和格子雷诺数。宏观雷诺数经验公式为

(9)

通常，格子速度Ulattice的数量级保持在0.1或者0.2，取值与宏观速度无关。文中Ulattice选择0.1，微通道截面的实际尺寸为300 μm，其格子单位取值为14，匹配宏观雷诺数和格子雷诺数后得到密度ρ的初始值。其中，根据微通道截面形状，a的取值为0.212 1，b的取值为0.676 6。

2 结果及分析

2.1 分叉角度对微通道内流体压力降的影响

图3所示是模拟趋于稳定后入口截面处与分叉点处的压力差值。数据整体分布在负值区域，表明从入口处流经到分叉点处时压力增大。从图中可以看出，90°分叉角度压力差值最小；当分叉角度大于90°时，压力差值随着分叉角度的增大而减小；当分叉角度小于90°时，压力差值随着分叉角度的增大而增大；压力差值最大值出现在中线偏向分叉点壁面一侧，越靠近壁面压力差值越小。从图中还可以发现，当分叉角度互补时，其压力差值结果相近，相差分别为10.72%和12.51%，且分叉角度小于90°的微通道压力差值大于分叉角度大于90°的微通道。

图3 不同分叉角度下分叉点截面处的 压力降(w=0.1%,v=2 m/s)Fig 3 Pressure drop at the cross section of the bifurcation point in different bifurcation angles(w=0.1%,v=2 m/s)

出口1截面处的压力降如图4所示，整体分布趋势呈抛物线形，压力降最大值出现在微通道中间截面处，越靠近壁面，其压力降值越小;90°分叉角度时压力降最小，压力降最大值出现在分叉角度为75°和105°处，最大值相差4.96%。从图4(a)、(b)中可以看出，入口速度为2 m/s时，压力降的最大值相差4.96%，入口速度为4 m/s时，压力降的最大值相差8.94%。从图4(a)、(c)中可以看出，溶液质量分数为0.1%时，压力降的最大值相差4.96%，溶液质量分数为0.3%时，压力降的最大值相差22.57%，说明随着溶液质量分数的增大，分叉角度对出口1处的压力降影响更大。因此，相比于入口速度，溶液质量分数使得分叉角度对出口1压力降的影响更为显著。

出口2截面处的压力降如图5所示，不论分叉角度如何，其压力降趋势一致。

2.2 入口速度对微通道内压力降的影响

图6所示是模拟趋于稳定后入口截面处与分叉点处的压力差值，可以看出，入口速度为10 m/s时压力差值最小，为2 m/s时压力差值最大，即入口速度越大，压力差值越小；而截面处压力差值最大值出现在中线偏向分叉点壁面一侧，越靠近壁面压力差值越小。

图6 不同入口速度下分叉点截面处的 压力降(w=0.1%,θ=120°)Fig 6 Pressure drop at the cross section of the bifurcation point at different inlet velocities(w=0.1%,θ=120°)

出口1截面处的压力降如图7所示。从图7(a)、(b)中可以看出，当溶液质量分数一定时，分叉角度为120°时，压力降的最大值相差2.32%，分叉角度为105°时，压力降的最大值相差5.92%。从图7(a)、(c)中可以看出，当分叉角度一定时，溶液质量分数为0.1%的情况下，压力降的最大值相差2.32%，溶液质量分数为0.3%时，压力降的最大值相差18.42%。可见，随着溶液质量分数的增大，使入口速度对出口1处的压力降影响更大。所以相比于分叉角度，溶液质量分数使入口速度对出口1压力降的影响更为显著。

图7 不同入口速度下出口1截面处的压力降Fig 7 Pressure drop at the cross section of the outlet 1 at different inlet velocities (a) w=0.1%,θ=120°； (b)w=0.1%,θ=105°;(c)w=0.3%,θ=120°

出口2截面处的压力降如图8所示。从图8(a)、(b)中可以看出，不同入口速度条件下出口2截面处压力降基本相同。从图8(a)、(c)中可以看出，当分叉角度一定时，溶液质量分数为0.1%时，不同入口速度下的微通道出口2截面处的压力降趋势一致，而溶液质量分数为0.3%时，不同入口速度下压力降的最大值相差2.69%，且压力降随着入口速度的增大而减小。可见随着溶液质量分数的增大使入口速度对出口1处的压力降影响更大，所以相比于分叉角度，溶液质量分数使入口速度对出口2压力降的影响更为显著。

图8 不同入口速度下出口2截面处的压力降Fig 8 Pressure drop at the cross section of the outlet 2 at different inlet velocities (a) w=0.1%,θ=120°； (b)w=0.1%,θ=105°;(c)w=0.3%,θ=120°

2.3 溶液质量分数对微通道内压力降的影响

图9所示是模拟趋于稳定后入口截面处与分叉点处压力差值，可以看出，水的压力差值比任一质量分数剪切稀化流体的压力差值都小；溶液质量分数为0.5%时压力差值最大，随着溶液质量分数增大，压力差值增大。溶液质量分数为0.3%和0.5%时，微通道分叉点截面处的压力差值基本相同，分析其原因是在这2种质量分数下，描述其非牛顿特性的稠度系数K和幂律指数n的值相近，其流变特性的相似导致了流动过程中产生的压力差值相近。

图9 不同溶液质量分数下分叉点截面 处的压力降(v=2 m/s,θ=120°)Fig 9 Pressure drop at the cross section of the bifurcation point at different solution concentrations(v=2 m/s,θ=120°)

出口1截面处的压力降如图10所示，压力降的最大值相差3.01%;溶液质量分数为0.5%时，压力降最大，随着溶液质量分数的降低，压力降减小。由于水是非牛顿流体，没有流变特性，导致水的压力降远远小于剪切稀化流体的压力降。

图10 不同溶液质量分数下出口1截面处的 压力降(v=2 m/s,θ=120°)Fig 10 Pressure drop at the cross section of the outlet 1 at different solution concentrations(v=2 m/s,θ=120°)

出口1处压力值随时间的变化趋势如图11所示。根据图11将模拟过程分为3个阶段，在模拟初期(0～1 000)，压力值随着溶液质量分数的增大而增大，且溶液质量分数为0.3%与0.5%的微通道内压力峰值差值为27.5%，远远小于与溶液质量分数为0.1%的压力峰值差值55.21%。在这一阶段，水的压力值已经达到稳定。随着模拟的进行(1 000～3 000)，溶液质量分数大的微通道内，压力值迅速降低，但是在水的模拟过程中，模拟开始不久就达到稳定，没有压力值降低的阶段。所以剪切稀化流体在流动过程中产生了压力值大幅度降低的现象，分析其原因是由于剪切稀化流体的流变特性导致的，因为其黏度随着剪切速率的增加而减小，导致压力也随之减小。第三阶段即模拟稳定后(3 000～5 000)，压力值的降低幅度随着溶液质量分数的增大而增大，与出口截面处压力降结果一致。

图11 不同溶液质量分数下出口1处压力随时间的 变化趋势(v=2 m/s,θ=120°)Fig 11 The trend of pressure over time at the outlet 1 at different solution concentrations(v=2 m/s,θ=120°)

出口2截面处的压力降如图12所示。不同溶液质量分数条件下，压力降的最大值相差3.55%；溶液质量分数为0.1%时，压力降最大；随着溶液质量分数的增大，压力降下降。

图12 不同溶液质量分数下出口2截面处 的压力降(v=2 m/s,θ=120°)Fig 12 Pressure drop at the cross section of the outlet 2 at different solution concentrations(v=2 m/s,θ=120°)

3 结论

通过将非牛顿流体幂律模型引入牛顿流体格子Boltzmann模型的方法，数值模拟了矩形截面微通道内剪切稀化流体的流动特性，探讨了流体流变特性和微通道几何构型对非牛顿流体流动行为的影响，得到如下结论：

(1)分叉角度为90°的微通道压力降最小，当分叉角度大于90°时，压力降随着分叉角度的增大而减小，当分叉角度小于90°时，压力降随着分叉角度的增大而增大，且当微通道分叉角度在数学概念上互补时，其压力降相近。

(2)微通道内各截面处压力降呈抛物线形，最大值出现在微通道中间截面处，越靠近壁面压力降越小;在溶液质量分数与分叉角度一定时，随着入口速度的增大，压力降增大；相比于分叉角度，溶液质量分数使得入口速度对出口处压力降的影响更为显著。

(3)在入口速度与分叉角度一定时，随着溶液质量分数的增大，压力降增大，且当流体的流变特性相近时，其流动特性相近;而相比于入口速度，溶液质量分数使得分叉角度对出口处压力降的影响更为显著。