摘要:基于不可压格子玻尔兹曼气-液两相流模型, 建立了一个新的非牛顿幂律流体气-液两相流模型, 并采用该模型研究了多孔介质内牛顿气体驱替非牛顿幂律流体液体的驱替问题, 主要探究了Ca数、动力黏度比M、固体表面润湿性θ、多孔结构几何类型及幂律指数n对驱替过程的影响. 研究发现: 不论被驱替液体为剪切变稀流体、牛顿流体还是剪切变稠流体都有随着Ca数增加, 驱替速度越快, 指进现象越明显, 驱替效率越低. 然而对于不同的幂率指数n, 驱替效率随Ca数的增加而减小的速率不同: n越大驱替效率随着Ca数增加而减小的速率越慢. 另一方面, 随着黏性比M增加, 驱替效率减小, 且幂律指数n越小, M对驱替效率的影响越大. 此外, 固体表面接触角θ对驱替过程的影响也和被驱替流体的幂率指数n相关, 虽然对于n > 1和 n < 1的情况都有随着多孔介质固体表面接触角θ增加, 驱替过程受到的阻力越小, 指进越来越不明显, 驱替完成时间和驱替效率增加, 然而当n > 1时, 随着n的增加, 接触角对驱替过程的影响越来越小. 还研究了孔隙率相同的情况下, 孔隙几何类型不同时的驱替过程, 从数值结果可以发现与多孔结构为圆形和方形障碍物相比, 当多孔结构的几何类型为三角形时, 指进现象最明显, 驱替效率最低.
关键词: 幂律两相流体/
格子Boltzmann模型/
不混溶驱替/
多孔介质

English Abstract

全文HTML

--> --> -->

剪切应力与剪切应变率之间不呈线性关系的非牛顿流体广泛存在于我们的日常生活中, 如牛奶、蛋青、血液、建筑中的泥浆、果酱等. 此外, 非牛顿流体对应的气-液两相流动在废水处理、石油开采、塑料泡沫加工以及环境保护等大量工程应用和科学研究中有着重要的作用, 例如, 在石油开采过程中^[1,2], 原油作为一种典型的非牛顿流体, 其与注入地下深处储油层的水或二氧化碳构成了非牛顿多相渗流问题. 随着能源短缺日益严峻, 开展多相非牛顿渗流相关问题的研究具有重要的科学和现实意义.
近几十年来, 众多****对微通道内的非牛顿两相流流动机理进行了一系列研究, 且大多数研究都基于幂率流体. Xu等^[3]通过实验研究了单一气泡在剪切变稀幂律流体中的上升行为, 发现剪切变稀幂律流体的流变性质对气泡上升的影响非常显著, 剪切变稀幂律流体的幂律指数与羧甲基纤维素(CMC)溶液的质量百分数有关, 其溶液质量百分数越高, 对应的幂律指数越低; 在质量百分数较低的CMC溶液中, 气泡上升过程根据其形状变化可分为三个阶段: 第一阶段, 气泡脱离喷嘴后, 形状基本为球形, 并逐渐转化成扁椭球形; 第二阶段, 气泡形状变形不连续; 最后阶段, 气泡变形连续; 随着CMC溶液质量百分数的增加, 过渡阶段逐渐消失; 在质量百分数最大的CMC溶液中, 气泡首先是平坦阶段的连续形状, 然后是形状基本不变以恒定速度沿直线上升. Du等^[4]采用实验方法研究了在流体聚焦装置中分散螺纹对剪切变稀幂律流体中液滴破裂行为的影响, 发现剪切变稀幂律流体的流变特性对液滴的破裂影响较小. Fu等^[5]通过实验研究了微通道内气泡在幂律流体中的融合过程. 在研究中观察到了两种气泡融合机理: 第一种是两个气泡融合成一个气泡, 第二种是一个气泡先分裂成两个气泡, 然后再融合成一个气泡, 该工作还分析了气泡融合的概率、气泡的位置随时间变化的关系以及气泡融合时周围的流场变化. 采用实验方法, Salehi等^[6]研究了含有聚合物的牛顿流体对牛顿液滴和非牛顿幂律流体液滴形成过程中的表面差异的影响, 发现在液滴增长的初始阶段, 牛顿流体和非牛顿幂律流体的行为相似, 然而在液滴颈缩过程中, 随着变形速率的变化, 牛顿流体和幂律流体之间存在明显的差异.
以上****运用实验方法对非牛顿流体两相流进行了研究. 随着计算机科技和数值方法的迅速发展, 采用数值方法研究复杂系统流动问题被广大****认可. 近年来, 有大量****运用数值方法研究了非牛顿流体两相流问题. Sontti等^[7]运用计算流体动力学(computational fluid dynamics, CFD)方法对T型微通道内牛顿液滴在非牛顿幂律流体气体中的形成过程进行了研究, 主要分析了幂律流体幂律指数、牛顿流体和幂律流体流量变化以及两相表面张力对液滴形成尺寸的影响, 发现了幂律流体幂律指数越大, 产生的液滴尺寸越小; 表面张力越大, 产生的液滴尺寸越小; 并发现了牛顿相和非牛顿幂律流体相的流量比对产生的液滴尺寸有显著影响. 相比于传统CFD方法, 格子Boltzmann方法(LBM)作为近几十年发展并不断完善的一种介观数值算法, 因其易于处理流体之间以及流体和固体壁面之间的相互作用、计算效率高以及不需要追踪界面等优点, 被广泛用来研究多相流问题. 目前已有众多****用LBM研究了牛顿两相流问题^[8—12]和非牛顿流体多相流问题^[13—23]. 谢驰宇等^[13]基于两相流自由能模型^[14], 构建了非牛顿流体两相流模型, 并用构建的模型研究了非牛顿流体驱替牛顿流体越过简单障碍物和多孔介质时的流动机理. Shi和Tang^[15]基于相场多相流模型^[16], 分析了光滑微通道内牛顿流体驱替非牛顿幂律流体的指进现象, 他们分析了被驱替相的幂律指数、固壁润湿性与$Bo$数对驱替指进现象的影响; 接下来, Shi和Tang^[17]又分析了在含多孔介质的微通道内牛顿液体驱替幂律流体液体的流动机理, 重点研究了幂律流体的幂律指数、毛细数、黏性比、障碍物表面润湿性以及多孔介质的排列顺序对指进现象的影响. Ba等^[18]基于颜色模型^{[19, 20]}, 构建了一个正规则格子Boltzmann (RLB)的颜色梯度模型, 并用该模型模拟两相幂律流体在两个平板之间的流动问题以及牛顿液滴在幂律流体中的变形和破裂问题. 此外, 闵琪等^[21]基于Shan-Chen (S-C)模型^[22,23]构建了幂律流体两相流格子Boltzmann模型, 模拟了牛顿液滴以及幂律流体液滴在固壁面的铺展情况.
以上研究工作提出了一些非牛顿气液两相流模型, 并用提出的模型研究了微通道内非牛顿两相流的一些运动机理. 但现有模型存在一些不足: 由于自由能模型不满足伽利略不变性^[24,25], 因此基于该模型提出的非牛顿两相流模型^[14]在数值模拟过程中会产生一些非物理效应; 基于颜色模型和基于相场模型改进的非牛顿两相流模型只能模拟密度比为1∶1的情况^[15,18]; 而在S-C模型中, 平衡态性质和网格相关^[26]; 此外, 现有的非牛顿两相流模型中流体压力的演化与流体密度直接相关, 容易产生数值不稳定. 而He等^[27]基于Enskong理论提出的不可压气-液两相流模型相对于其他两相流模型具有物理机理更清晰, 压力的演化和密度的演化相独立, 处理复杂界面变化时数值稳定性更好等优点^[25,28—30]. 鉴于此, 本文基于幂率流体, 在He等^[27]提出的针对牛顿流体的不可压多相流模型基础上, 提出了不可压非牛顿幂律流体气-液两相流模型. 此外目前对于非牛顿幂律流体两相流问题的研究大都局限在简单通道内, 而对于多孔介质内牛顿气体驱替非牛顿幂律流体液体的机理研究尚不充分. 因此, 本文又用发展的非牛顿幂律流体模型研究了多孔介质内牛顿气体驱替非牛顿幂律流体液体的问题.

He等直接从不可压多相流的Boltzmann方程出发, 提出了一种稳定性好的适用于牛顿流体的不可压多相流模型^[27], 简称为HCZ^[31—35]模型, 下面我们将该模型推广到幂律流体.
在HCZ模型中分别采用两个分布函数f_α和g_α描述指标参数和速度/压力的演化过程, 其对应的分布函数${f_\alpha }$和${g_\alpha }$分别为

$\begin{split}&{f_\alpha }\left( {{{x}} + {{{e}}_\alpha }{\delta _t},t + {\delta _t}} \right) - {f_\alpha }\left( {{{x}},t} \right)\\ =\; & - \frac{{{f_\alpha }\left( {{{x}},t} \right) - f_\alpha ^{{\rm{eq}}}\left( {{{x}},t} \right)}}{\tau } \\& - \frac{{\left( {2\tau - 1} \right)}}{{2\tau }}\frac{{\left( {{{{e}}_\alpha } - {{u}}} \right) \cdot \nabla \psi \left( \phi \right)}}{{RT}}{\varGamma _\alpha }\left( {{u}} \right){\delta _t},\\& {g_\alpha }\left( {{{x}} + {{{e}}_\alpha }{\delta _t},t + {\delta _t}} \right) - {g_\alpha }\left( {{{x}},t} \right)\\= \; & - \frac{{{g_\alpha }\left( {{{x}},t} \right) - g_\alpha ^{{\rm{eq}}}\left( {{{x}},t} \right)}}{\tau } \\& + \frac{{2\tau - 1}}{{2\tau }}\left( {{{{e}}_\alpha } - {{u}}} \right) \cdot [{\varGamma _\alpha }\left( {{u}} \right)\left( {\kappa \rho \nabla {\nabla ^2}\rho } \right)\\& - \left( {{\varGamma _\alpha }\left( {{u}} \right) - {\varGamma _\alpha }\left( 0 \right)} \right)\nabla \psi \left( \rho \right)]{\delta _t},\end{split}$

式中$\alpha =0\;{\rm{,}}\;1\;{\rm{,}}\;2\;{\rm{,}}\; \cdots \;{\rm{,}}\;b - 1,$$b$为离散速度方向个数; ${{x}}$和$t$分别表示粒子运动的位置和时间; ${\delta _t}$代表时间步长; $\tau $为松弛时间, 其与运动黏度$\nu $相关;$\phi,\rho,{{u}}$分别代表指标函数、流体密度和流体速度; $\kappa $为决定表面张力大小的参数; $\psi \left( \rho \right) = p - \rho {c_s}^2$, 其中$p$为流体压力, 演化方程中$\psi \left( \phi \right)$和状态方程相关; 在方程(1)中函数${\varGamma _\alpha }\left( u \right)$表达式为

${\varGamma _\alpha }\left( {{u}} \right) = {\omega _\alpha }\left[ {1 + \frac{{3{{{e}}_\alpha } \cdot {{u}}}}{{{c^2}}} + \frac{{9{{\left( {{{{e}}_\alpha } \cdot {{u}}} \right)}^2}}}{{2{c^4}}} - \frac{{3{{{u}}^2}}}{{2{c^2}}}} \right],$

其中${\omega _\alpha }$代表权重系数. 演化方程中$f_\alpha ^{{\rm{eq}}}\left( {{{x}},t} \right)$和$g_\alpha ^{{\rm{eq}}}\left( {{{x}},t} \right)$是分布函数对应的平衡态, 其形式如下:

$\begin{split}& f_\alpha ^{{\rm{eq}}} = {\omega _\alpha }\phi \left[ {1 + \frac{{3{{{e}}_\alpha } \cdot {{u}}}}{{{c^2}}} + \frac{{9{{\left( {{{{e}}_\alpha } \cdot {{u}}} \right)}^2}}}{{2{c^4}}} - \frac{{3{{{u}}^2}}}{{2{c^2}}}} \right], \\ & g_\alpha ^{{\rm{eq}}} = {\omega _\alpha }\left[ {p \!+\! \rho \left( {\frac{{3{{{e}}_\alpha } \cdot {{u}}}}{{{c^2}}} \!+\! \frac{{9{{\left( {{{{e}}_\alpha } \cdot {{u}}} \right)}^2}}}{{2{c^4}}} \!-\! \frac{{3{{{u}}^2}}}{{2{c^2}}}} \right)}\!\right].\end{split}$

宏观量指标参数$\phi $、压力$p$以及速度${{u}}$的统计由下面方程给出:

$\begin{split}& \phi = \sum {{f_\alpha }},\\& p = \sum {{g_\alpha }} - \frac{1}{2}{{u}} \cdot \nabla \psi \left( \rho \right){\delta _t},\\& \rho RT{{u}} = \sum {{e_\alpha }} {g_\alpha }{\rm{ + }}\frac{{RT}}{2}\left( {\kappa \rho \nabla {\nabla ^2}\rho } \right){\delta _t}.\end{split}$

流体密度$\rho \left( \phi \right)$和运动黏度$\nu \left( \phi \right)$可由指标参数$\phi $计算得到

$\rho \left( \phi \right) = {\rho _{\rm{g}}} + \frac{{\phi - {\phi _{\rm{l}}}}}{{{\phi _{\rm{h}}} - {\phi _{\rm{l}}}}}\left( {{\rho _{\rm{l}}} - {\rho _{\rm{g}}}} \right),$

$\nu \left( \phi \right) = {\nu _{\rm{g}}} + \frac{{\phi - {\phi _{\rm{l}}}}}{{{\phi _{\rm{h}}} - {\phi _{\rm{l}}}}}\left( {{\nu _{\rm{l}}} - {\nu _{\rm{g}}}} \right),$

其中ρ_g和ρ₁分别代表气相流体和液相流体密度, ${\phi _{\rm{h}}}$和${\phi _{\rm{l}}}$为指标参数的最大值和最小值, 可根据状态方程由Maxwell重构得到.
以上模型只适用于剪切应力与剪切应变率之间呈线性关系的牛顿流体的多相流动问题, 为了将该模型推广到非牛顿流体, 本文用幂率流体模型来体现流体的非牛顿特性. 对于幂率流体, 动力黏度$\eta $的表达式为

$\eta = \eta_0 {\gamma ^{n - 1}} = {\eta _0}{\left( {{S_{\alpha \beta }}{S_{\alpha \beta }}} \right)^{\left( {n - 1} \right)/2}},$

其中$\gamma $为剪切速率, ${\eta _0}$为稠度系数,

$\mathop S\nolimits_{\alpha \beta } = \frac{1}{2}\left( {\frac{{\mathop {\partial u}\nolimits_\alpha }}{{\mathop {\partial u}\nolimits_\beta }} + \frac{{\mathop {\partial u}\nolimits_\beta }}{{\mathop {\partial u}\nolimits_\alpha }}} \right)$

为应变张量, n为非牛顿流体幂律指数. 当$ n < 1$时, 流体为剪切变稀流体, 即其动力黏度$\eta $随着剪切速率$\gamma $的增大而减小; 当n > 1时, 流体为剪切变稠流体, 其动力黏度$\eta $随着剪切速率$\gamma $的增大而增大; 当n = 1时, 流体就是通常的牛顿流体, 其动力黏度$\eta $保持为一个定值, 此时$\eta = {\eta _0} = \rho v$. 因此, 根据n的取值来区分流体是牛顿流体还是非牛顿幂律流体.
通过对分布函数${g_\alpha }$进行Chapman-Enskog分析可以得到应变张量${S_{\alpha \beta }}$与分布函数之间有如下关系:

${S_{\alpha \beta }} = - \frac{{\displaystyle\sum\limits_{\alpha = 0}^8 {{e_\alpha }{e_\alpha }{g_\alpha } - p{\delta _{\alpha \beta }} - \rho {u_\alpha }{u_\beta }} }}{{\tau {\delta _t}RT\rho }}.$

从(7)式可以看出, 在该幂律流体模型中, 应变张量${S_{\alpha \beta }}$可直接用分布函数和局部点宏观量信息计算得到, 避免了差分计算, 可提高数值稳定性. 通过Chapman-Enskog分析还可以得到方程(1)对应的宏观动力学方程为

$\begin{split}& \frac{{\partial \phi }}{{\partial t}} + \nabla \cdot \left( {\phi {{u}}} \right) = - \lambda \nabla \cdot \left[ {\frac{\phi }{\rho }\nabla p\left( \rho \right) - \nabla p\left( \phi \right)} \right], \\&\quad\;\;\quad\quad\quad\quad\quad\quad \frac{1}{{\rho RT}}\frac{{\partial p}}{{\partial t}} + \nabla \cdot {{u}} = 0, \\& \rho \left[ {\frac{{\partial {{u}}}}{{\partial t}} \!+\! \left( {{{u}} \cdot \nabla } \right){{u}}} \right] \!= \!- \nabla p + \nabla \cdot { \varPi}+\! \kappa \rho \nabla {\nabla ^2}\rho,\end{split}$

其中$ \varPi $是黏性应力张量, ${ \varPi} = \rho \nu \left( {\nabla {{u}} + {{u}}\nabla } \right)$. 运动黏度系数和松弛时间$\tau $之间关系为$\nu =(\tau - $ $ 0.5){c_{\rm{s}}}^2{\delta _t}$, 这里${c_{\rm{s}}}^{\rm{2}} = {c / {\sqrt 3 }}$是与格子速度$c = {{{\rm{d}}x} / {{\rm{d}}t}}$相关的常数.
本文使用二维九数(${\rm{D2Q9}}$)模型来进行数值模拟研究, 虽然本文是开展二维研究, 但已有文献[36, 37]表明, 尽管二维得到的结果和实际三维的结果定量上可能会有所不同, 但是定性上趋势一致. 在D2Q9模型中权重系数${\omega _\alpha }$设置为: 当α = 0时, ${\omega _0} = \dfrac{4}{9}$; 当α = 1—4, ${\omega _\alpha } = \dfrac{1}{9};$当α = 5—8, ${\omega _\alpha } = \dfrac{1}{{36}}.$${{{e}}_\alpha }$为离散速度, 其表达式为^[38]

${e_\alpha } = \left\{ \begin{aligned}& 0, \quad \quad \quad \quad \quad \quad\quad \quad\quad\quad\quad\;\; \alpha = 0,\\& \left( {\cos \left[ {\left( {\alpha - 1} \right){\text{π}}/2} \right],\sin \left[ {\left( {\alpha - 1} \right){\text{π}}/2} \right]} \right)c, \\&\quad\quad\quad\quad \quad \quad \quad \quad\quad\quad\quad\quad \alpha = 1{\text{—}}4,\\& \sqrt 2(\cos \left[ {\left( {\alpha - 5} \right){\text{π}}/2 + {\text{π}}/4} \right],\\ & \sin \left[ {\left( {\alpha - 5} \right){\text{π}}/2 + {\text{π}}/4} \right])c, \quad \alpha = 5{\text{—}}8.\end{aligned} \right.$

演化方程中梯度和拉普拉斯算子的离散方法均采用二阶中心各向同性方法(ICS)^[39].

$\begin{split} & \nabla \chi \left( x \right) \approx {\nabla _{\rm{c}}}\chi \left( x \right) = \sum\limits_{i \ne 0} {\frac{{{\omega _i}{c_i}\chi \left( {x + {c_i}{\delta _t}} \right)}}{{c_{\rm{s}}^{\rm{2}}{\delta _t}}}{\rm{, }}} \\& {\nabla ^2}\chi \left( x \right) \approx {\nabla ^2}_{\rm{c}}\chi \left( x \right) = \sum\limits_{i \ne 0} {\frac{{2{\omega _i}\left[ {\chi \left( {x + {c_i}{\delta _t}} \right) - \chi \left( x \right)} \right]}}{{c_s^2{\delta _t}^2}}} .\end{split} $

演化方程中$\psi \left( \phi \right)$由状态方程决定, 在本文中采用Carnahan-Starling状态方程^[27], 其对应的$\psi \left( \phi \right)$可写为如下形式:

$\psi \left( \phi \right)={\phi ^2}RT\frac{{4 - 2\phi }}{{{{\left( {1 - \phi } \right)}^3}}} - a{\phi ^2},$

其中$a$决定分子间相互吸引力强度, $R$为气体常数, $T$为流体温度.

2

-->

3.1.拉普拉斯定律验证

本小节采用Laplace定律来验证模型的正确性. 初始时在网格数为$128 \times 128$区域中心内放置半径$ r $, 密度${\rho _{\rm{l}}}{\rm{ = }}\;0.5$的静止幂律流体圆形液滴, 其余区域充满着密度为${\rho _{\rm{g}}}{\rm{ = }}\;0.1$的牛顿气体. 计算域四周的边界条件均为周期性边界条件. 根据Laplace定律可知, 当系统达到稳定时表面张力$\sigma $恒定, 且液滴内外的压力差${P_{\rm{i}}} - {P_{\rm{o}}}$与半径的倒数$ {1 / r} $呈线性关系, 关系式如下:

$ {P_{\rm{i}}} - {P_{\rm{o}}} = \frac{\sigma }{r}. $

为了验证Laplace定律, 在数值模拟中分别考虑了$ r= $20, 22, 24, 26, 28和30六种情况. 为了保证数值模拟结果具有一般性, 对于每一种情况都研究了n = 0.7, n = 1.0与n = 1.3三种情形. 图1给出了在不同的幂律指数下得到的液滴内外的压力差${P_{\rm{i}}} - {P_{\rm{o}}}$与半径的倒数$ {1 / r} $的关系, 可知, 对于所有的n(流体分别为剪切变稀流体、牛顿流体和剪切变稠流体), 模拟结果与Laplace定律均一致.
图 1 液滴内外压力差${P_{\rm{i}}} - {P_{\rm{o}}}$和半径倒数1/r之间的关系
Figure1. Relationship between pressure jump across the droplet interface${P_{\rm{i}}} - {P_{\rm{o}}}$ and inverse of droplet radius 1/r.

2 -->

3.2.润湿性固壁面静态接触角验证

润湿现象在自然界中广泛存在, 如水滴在玻璃板上将迅速铺展开, 而水银滴在玻璃板上将呈现为球滴状, 这就是润湿性不同程度的结果. 润湿性反映流体和固体之间相互作用力的强度. 在复杂微通道内, 其是影响气-液-固或液-液-固三相动态行为的重要指标参数. 在格子Boltzmann方法(LBM)中, 通过润湿性边界条件来描述壁面的润湿性. 本文考虑复杂微通道内固体表面润湿性对牛顿气体驱替液相非牛顿幂律流体的驱替动态影响, 润湿性边界条件采用Davies等^[40]提出的方法, 在该方法中壁面的润湿强度采用表面亲和性${\alpha _{\rm{s}}}$来刻画, 并把表面亲和性与指标参数$\phi $联系起来, 其关系式为

${\alpha _{\rm{s}}}=\frac{{\phi - {\phi _{\Sigma} }}}{{{\phi _l} - {\phi _{\Sigma}}}},$

其中${\phi _{\Sigma} } = \dfrac{1}{2}\left( {{\phi _{\rm{h}}} + {\phi _{\rm{l}}}} \right)$. 则${\alpha _{\rm{s}}}$取值范围位于–1到1之间, ${\alpha _{\rm{s}}}= - 1$对应完全亲气表面, ${\alpha _{\rm{s}}} = 1$对应完全亲水表面. 静态接触角${\theta _{{\rm{eq}}}}$与${\alpha _{\rm{s}}}$关系式为

$\cos \left( {{\theta _{{\rm{eq}}}}} \right) = \frac{{{\sigma _{{\rm{s}}2}} - {\sigma _{{\rm{s}}1}}}}{{{\sigma _{12}}}} = \frac{{{\alpha _{\rm{s}}}}}{2}\left( {3 - \alpha _{\rm{s}}^{\rm{2}}} \right),$

其中${\sigma _{12}}$为气-液表面张力; ${\sigma _{{\rm{s}}1}}$, ${\sigma _{{\rm{s2}}}}$分别代表固-液表面张力与固-气表面张力.
下面对幂律流体静态液滴的静态接触角进行验证. 数值模拟中网格数为128 × 65, 在计算区域下边界中心处放置半径$ r = 20 $, 密度${\rho _{\rm{l}}}=0.5$以及幂律指数n = 0.7的非牛顿幂律流体半圆液滴, 液滴周围充满了${\rho _{\rm{g}}}{\rm{ = }}\;0.1$的牛顿气体, 初始两相的运动黏度均为${\nu _{\rm{g}}} = {\nu _{\rm{l}}} = 1/6$. 计算域的边界条件设置为: 上下壁面是无滑移边界条件, 左右是周期边界条件.
图2展示了当数值模拟达到稳定时, 壁面静态接触角${\theta _{{\rm{eq}}}}$分别为$60^\circ $, $90^\circ $与$120^\circ $时, 模拟所得到的稳态接触角$\theta $, 其结果分别为57.6°, 87.3°与117.7°. 模拟结果和理论值之间的相对误差分别为4.0%, 3.0%与1.9%. 图3展示了数值模拟得到的壁面稳态接触角$\theta $与固壁面上的指标参数${\phi _{{\rm{wall}}}}$的关系, 结果与(13)和(14)式符合得较好.
图 2 不同初始静态接触角${\theta _{{\rm{eq}}}}$时得到的稳态接触角$\theta $　(a) ${\theta _{{\rm{eq}}}}{\rm{ = }}{60^{\rm{o}}}$; (b) ${\theta _{{\rm{eq}}}}{\rm{ = }}{90^{\rm{o}}};$ (c) ${\theta _{{\rm{eq}}}}={120^{\rm{o}}}$
Figure2. Steady state contact angles $\theta $ obtained with the different values of static contact angles ${\theta _{{\rm{eq}}}}$: (a) ${\theta _{{\rm{eq}}}}{\rm{ = }}{60^{\rm{o}}}$; (b) ${\theta _{{\rm{eq}}}}{\rm{ = }}{90^{\rm{o}}};$ (c) ${\theta _{{\rm{eq}}}}{\rm{ = }}{120^{\rm{o}}}$.

图 3 稳态接触角$\theta $与指标参数${\phi _{{\rm{wall}}}}$的线性关系
Figure3. Linear relationship between steady state contact angle $\theta $ and the order parameter of a solid wall ${\phi _{{\rm{wall}}}}$.

2 -->

3.3.T型通道内液滴分离

本小节验证在不同$Ca$数下, T型微通道内形成液滴的大小. T型微管道结构如图4所示, 其中W₀ = W₁ = 30, L = 520, Y₁ = 75, Y₂ = 120. 初始分散相的子管道充满牛顿液相流体, 连续相主通道充满幂律流体, 其幂律指数n = 0.4. 边界条件设置为: 连续相与分散相的进口是速度进口边界, 连续相出口是对流边界条件^[30], 管径的固壁面都采用无滑移边界条件. 图5给出了在不同的$Ca$数下, 分离的牛顿流体液滴(黄色). 图6给出了在不同的$Ca$数下, 分离的牛顿流体液滴尺寸, 并与数值结果^[41]进行对比, 得到了一致的结果.
图 4 T型通道问题物理模型
Figure4. Physical model for the case of T shape channel.

图 5 不同Ca数对应的液滴形态　(a) Ca = 0.06370; (b) Ca = 0.06835; (c) Ca = 0.07300; (d) Ca = 0.07750; (e) Ca = 0.0820; (f) Ca = 0.08650; (g) Ca = 0.0910
Figure5. Droplet morphology obtained under various values of Ca: (a) Ca = 0.06370; (b) Ca = 0.06835; (c) Ca = 0.07300; (d) Ca = 0.07750; (e) Ca = 0.0820; (f) Ca = 0.08650; (g) Ca = 0.0910.

图 6 在剪切变稀幂律流体中, 不同的$Ca$数下形成液滴的无量纲直径(其中D是形成的液滴的直径, H是管径的直径)
Figure6. Droplet dimensionless diameters at different values of $Ca$ in shear thinning power-law fluid. D is diameters of the droplet and H is width of the main channel.

本节研究多孔介质内幂律流体气液两相流驱替问题, 所得到的结论对石油开采、二氧化碳地质埋存等相关问题均有一定的指导意义, 由于本文主要关注牛顿流体驱替非牛顿幂律流体的机理研究, 如最近文献[42—46]中的处理方法一样, 将实际地质结构简化为一种理想的多孔结构, 其示意图见图7. 计算域的网格数为M × N, 多孔障碍物的尺寸分别为A与B, 两障碍物中心线之间的水平距离为X, 垂直距离为Y, 最靠近进口的障碍物中心离进口的水平距离为S₁. 初始时刻在$x < S$处充满密度${\rho _{\rm{g}}} = 0.1$, 运动黏度${\nu _{\rm{g}}} = 1/6$, 动力黏度${\eta _{\rm{g}}} = $ 0.01667的牛顿流体(蓝色), 在$x > S$处充满被驱替液相(黄色)非牛顿幂律流体, 其密度为${\rho _{\rm{l}}} = 0.5$. 边界条件设置为: 左边入口处采用速度边界条件, 右边出口处采用出口边界条件^[47], 多孔表面与上下固壁面均采用无滑移边界条件.
图 7 多孔介质驱替模型
Figure7. The model for porous media displacement problem.

需要特别指出的是下文中需要用到如下两个变量: 无量纲毛细数$Ca$, 其定义为$Ca = u{\eta _{\rm{g}}}/\sigma $, 其中u是驱替相的速度, ${\eta _{\rm{g}}}$是驱替相动力黏度, $\sigma $是气液两相表面张力; 驱替效率$De$, 其定义为当驱替流体到达出口时通道内剩余的被驱替相体积与初始被驱替相体积的比值.
2

-->

4.1.Ca数的影响

本小节研究Ca数对不混溶幂律流体驱替过程的影响. 在数值模型中M × N = 240 × 600, A = 30, B = 20, S = 6, S₁ = 40, X = 60, Y = 60, 孔隙率$\xi = 0.8507$. 初始动力黏度比M = 5.0, 即被驱替液初始运动黏度${\nu _{\rm{l}}} = 1/6$, 动力黏度${\eta _{\,\rm{l}}} = 0.08333$, $\sigma = 0.0056$. 固体表面均为中性润湿($\theta = {90^{\rm{o}}}$).
图8给出了在不同Ca数下被驱替液为剪切变稀、牛顿以及剪切变稠三种流体时得到的驱替完成时指进形态图. 从图8可以看出, 不论被驱替液是剪切变稀(图8(a)—(c) n = 0.7)、牛顿(图8(d)—(f) n = 1.0)还是剪切变稠(图8(g)—(i) n = 1.3)流体, 都有Ca数越大, 指进现象越明显^[17], 驱替完成时花费的时间越少. 具体而言, 当n = 0.7时(图8(a)—(c)), Ca = 0.0298对应的驱替完成时间t = 36.2, 当Ca数增加到0.0877时, 驱替时间减少到了t = 11.1, 减少了69.3%; 而当n = 1.0时(图8(d)—(f))和n = 1.3 (图8(g)—(i))时, Ca数从0.0298增加到0.0877, 驱替时间分别减少了68.8%和67.7%. 发生以上现象的原因是由于Ca数是表征黏性力与表面张力比值的无量纲参数, Ca数越大说明表面黏性力越大而表面张力越小, 而随着黏性力的增加, 驱替过程受到的阻力增加, 随着表面张力减小, 气液界面更容易发生变形, 因此Ca数越大指进现象越明显. Shiri等^[48]以及Liu等^[49]在研究多孔介质内的驱替问题时也得到了相似的结论. 另一方面, 从图8还可以发现当Ca数相同时, 当被驱替相的幂率指数n越大时, 指进现象越明显, 驱替完成所需时间越短^[17]. 例如当Ca = 0.0298时(图8(a)、图8(d)、图8(g)), n = 0.7, 1.0, 1.3对应的驱替完成时间分别为36.2, 32.4, 29.1; Ca = 0.0595时(图8(b)、图8(e)、图8(h)), n = 0.7, 1.0, 1.3对应的驱替完成时间分别为16.9, 15.3, 14.2; 而Ca = 0.0877时(图8(c)、图8(f)、图8(i)), n = 0.7, 1.0, 1.3对应的驱替完成时间分别为11.1, 10.1, 9.4. 因此对应这三种Ca数的情况, 随着幂率指数n的增加驱替时间分别减小了19.6%, 16.0%和15.3%. 也就是说, 幂率指数n越大, Ca数的增加导致的驱替时间减少的速率越来越慢. 导致该现象的原因是对于剪切变稠流体随着驱替过程的进行其动力黏度会大于初始时刻的动力黏度, 剪切变稀流体的动力黏度会小于初始时刻的动力黏度, 而牛顿流体的动力黏度会保持不变. 为了说明这一现象, 图9给出了Ca = 0.0298时对应驱替完成时刻的两相流体的动力黏度分布. 从图9可以看出, 对于研究的所有n的值, 当驱替过程完成时, 驱替气相的动力黏度一直保持在初始值0.01667附近, 而对应幂律指数n = 0.7 (图9(a))、1.0 (图9(b))以及1.3 (图9(c))的被驱替液得到的动力黏度分别为0.04409, 0.08333与0.1579. 即剪切变稀、牛顿以及剪切变稠流体被驱替过程中, 其分别对应的两相动力黏度比M小于5、等于5以及大于5. 而两相动力黏度比越大, 两相流体间黏性力影响变大, 即所受黏性阻力越大, 因此指进现象越明显, 驱替越快^[15,17].
图 8 不同的$Ca$数下, 被驱替液为剪切变稀、牛顿与剪切变稠流体时得到的指进形态图　(a)?(c) n = 0.7; (d)?(f) n =1.0; (g)?(i) n = 1.3
Figure8. Final finger patterns obtained under different values of Ca for shear thinning, Newtonian and shear thickening fluids: (a)? (c) n = 0.7; (d)?(f) n =1.0; (g)?(i) n = 1.3.

图 9 驱替完成时, 不同幂律指数情况下得到的气液两相动力黏度示意图　$({\rm{a}})\;n = 0.7$; $({\rm{b}})\;n = 1.0$; $({\rm{c}})\;n = 1.3$
Figure9. Schematic diagram of gas-liquid two phase dynamics viscosity obtained under different values of power-law exponent: $({\rm{a}})\;n = 0.7$; $({\rm{b}})\;n = 1.0$; $({\rm{c}})\;n = 1.3$.

为了定量分析Ca数以及幂律指数n对驱替过程的影响, 图10给出了不同Ca数以及幂律指数n下得到的驱替效率. 从图10可以看出, 不论被驱替液是剪切变稀、牛顿还是剪切变稠流体, 都有Ca数越大, 驱替效率De越低. 具体而言, n = 0.7时, Ca数从0.0298到0.0877时驱替效率De的值从0.744减小到0.688, 减小了7.52%; n = 1.0时, Ca数从0.0298到0.0877时驱替效率De的值从0.663减小到0.617, 减小了6.93%; n = 1.3时, Ca数从0.0298到0.0877时驱替效率De的值从0.590减小到0.550, 减小了6.78%. 从以上分析以及图10中曲线变化可以发现, 不论被驱替相是牛顿流体还是幂律流体, 驱替效率随着Ca数的增加而减小, 然而驱替效率减小的速率随着n的增加而减小. 另外, Ca数一定时, 幂律指数n越大, 驱替效率De越低.
图 10 $Ca$ 数和幂律指数n对幂律流体驱替效率的影响
Figure10. Effects of $Ca$and power-law exponent n on power-law fluid displacement efficiency.

2 -->

4.2.动力黏度$M$的影响

本小节研究动力黏度比M对不混溶幂律流体驱替过程的影响, 这里黏性比的增加是通过改变被驱替液的黏性实现的. 在本小节Ca数设置为0.0446, 其余参数与4.1节相同.
图11给出了被驱替液为剪切变稀、牛顿以及剪切变稠三种流体在不同初始动力黏度比M的情况下, 驱替完成时指进形态图. 从图11可以发现, 对于被驱替液是剪切变稀流体的情况(图11(a)—(c) n = 0.7), 黏度比从2.5增加到12.5, 驱替时间从26.5减少到19.7, 减少了25.7%. 对于被驱替液为牛顿流体(图11(d)—(f) n = 1.0)和剪切变稠(图11(g)—(i) n = 1.3)流体的情况, 动力黏度比从2.5增加到12.5时对应的驱替时间分别减少了19.6%和9.3%. 因此对于所有n的情况都有初始动力黏度比M越大, 指进现象越明显, 驱替完成时所花费的时间越少^[15,17]. 这是因为黏性比越大说明被驱替液黏性越大, 而一个轻流体(黏性较小)驱替一个重流体(黏性较大)是非常困难的, 因此指进现象会更明显. 其他****^[48—51]也发现了类似的现象. 从以上数据以及指进形态图还可以发现驱替相的幂律指数n越大, 黏性比的增加对驱替过程的影响越小. 另一方面从图11可以观察到当初始动力黏度比M相同时, 被驱替相的幂律指数n越大, 指进现象越明显, 对应的驱替完成所花费的时间越少. 这与4.1节流体越黏稠, 流体越难被驱替结论一致^[15,17]. 如M = 2.5 (图11(a)、图11(d)、图11(g))时, 对应n = 0.7, 1.0和1.3的驱替完成时间分别为26.5, 23.5和20.3, 此时随着幂率指数n的增加, 驱替完成时间减小了23.4%.
图 11 不同的动力黏性比M下, 被驱替液为剪切变稀、牛顿与剪切变稠流体时得到的指进形态图　(a)?(c) n = 0.7; (d)?(f) n = 1.0; (g)?(i) n = 1.3
Figure11. Final finger patterns obtained under different values of viscosity ratios M for shear thinning, Newtonian and shear thickening fluids: (a)?(c) n = 0.7; (d)?(f) n = 1.0; (g)?(i) n = 1.3.

M = 7.5 (图11(b)、图11(e)、图11(h))时, 对应以上三个幂率指数的驱替完成时间分别为21.4, 19.6和18.5. M = 12.5 (图11(c)、图11(f)、图11(i))时对应的驱替完成时间分别为19.7, 18.9和18.4. 因此对应这两种黏性比的情况, 随着幂率指数增加驱替时间分别减小了13.6%和6.6%. 从以上数据分析可以发现随着黏性比的增加, 驱替速率减少, 且当黏性比较小时, 幂律指数n越大, 流体越难被驱替, 而随着黏性比的增加, 被驱替相是牛顿流体和幂律流体的驱替结果之间的差异越来越小.
为了进一步分析初始动力黏度比M以及幂律指数n对驱替过程的影响, 图12给出了不同初始动力黏度比M以及幂律指数n下得到的驱替效率. 从图12可以看出, 当不论被驱替液是剪切变稀、牛顿还是剪切变稠流体, 都有动力黏度比M越大, 驱替效率De越低, 而且随着n的增加得到的驱替效率曲线的斜率越来越平缓, 说明黏性比M对于驱替效率的影响随着n的增加而减小. 具体而言, n = 0.7时, M从2.5到12.5时驱替效率De的值从0.827减小到0.596, 减小了27.93%; n = 1.0时, M从2.5到12.5时驱替效率De的值从0.730减小到0.562, 减小了23.01%; n = 1.3时, M从2.5到12.5时驱替效率De的值从0.626减小到0.535, 减小了14.54%. 因此, 不论被驱替相是牛顿流体还是幂律流体, 驱替效率随着M的增加而减小, 然而驱替效率减小的速率随着n的增加而减小. 另外, M一定时, 幂律指数n越大, 驱替效率De越低.
图 12 动力黏度比M和幂律指数n对幂律流体驱替效率的影响
Figure12. Effects of viscosity ratio M and power-law exponent n on power-law fluid displacement efficiency.

2 -->

4.3.润湿性的影响

本小节研究固体表面润湿性对幂律流体驱替效率的影响. 在数值模拟中动力黏度比$M = 5.0$, 毛细数$Ca = 0.0446$, 其余参数与4.2节相同.
图13给出了被驱替液为剪切变稀、牛顿以及剪切变稠三种流体在四种不同润湿性($\theta = {45^{\rm{o}}}, $ ${60^{\rm{o}}},{120^{\rm{o}}}$与135°)情况下, 驱替完成时得到的指进形态图. 从图13可以发现, 对于所有的n都有随着固体表面接触角$\theta $增加, 指进现象越不明显. 这是因为当固体面的接触角越大, 固体表面的疏水性越强, 则被驱替液受到固壁的黏附力越小, 所以被驱替液更容易驱替, 则指进现象越不明显. 同时由于接触角越大被驱替流体越容易被驱替, 表现为指进长度变短, 所以被驱替液到出口的时间, 也即驱替完成的时间增加. 因此驱替流体到达出口的时间越长. 其他****^[49,52]也发现了类似的现象. 具体而言, 对于所研究的四种接触角的情况($\theta = $45°, 60°, 120°与135°), 当n = 0.7 (图13(a)、图13(d)、图13(g)、图13(j))对应的驱替完成时间分别为23.1, 23.5, 25.5和26.8; n = 1.0 (图13(b)、图13(e)、图13(h)、图13(k))对应的驱替完成时间分别为20.7, 20.9, 22.2, 22.6; n = 1.3 (图13(c)、图13(f)、图13(i)、图13(l))对应的驱替时间分别为19.0, 19.1, 19.9和20.3. 通过计算可以得到接触角$\theta $从${45^{\rm{o}}}$增加到${135^{\rm{o}}}$时, n = 0.7, 1.0和1.3对应的驱替时间分别增加了16.0%, 9.2%和6.8%, 说明随着幂率指数n的增加, 固体表面润湿性对驱替完成时间段的影响逐渐减小. 另一方面, 从图13可以发现当固体表面润湿性$\theta $相同时, 被驱替相的幂律指数n越大, 指进越明显, 又一次说明了被驱替相流体越黏稠越不容易被驱替.
图 13 不同的润湿性角度$\theta $下, 被驱替液为剪切变稀、牛顿与剪切变稠流体时得到的指进形态图　(a)?(c) $\theta = {45^{\circ}}$; (d)? (f) $\theta = {60^{\circ}}$; (g)?(i) $\theta = {120^{\circ}}$; (j)?(l) $\theta = {135^{\circ}}$
Figure13. Final finger patterns obtained under different values of contact angles $\theta $ for shear thinning, Newtonian and shear thickening fluids: (a)?(c) $\theta = {45^{\circ}}$; (d)－(f) $\theta = {60^{\circ}}$; (g)?(i) $\theta = {120^{\circ}}$; (j)?(l) $\theta = {135^{\circ}}$.

为了定量分析固体表面接触角$\theta $以及幂律指数n对驱替过程的影响, 图14给出了不同固体表面润湿性以及幂律指数n下得到的驱替效率. 从图14可以看出无论剪切变稀、牛顿以及剪切变稠流体, 都有固体表面接触角越小, 驱替效率De越低. 具体而言, n = 0.7时, 固体表面接触角$\theta $从45°到135°时驱替效率De的取值为0.611到0.838, 增加了37.15%; n = 1.0时, 固体表面接触角$\theta $从45°到135°时驱替效率De的取值为0.530到0.703, 增加了32.64%; n = 1.3时, 固体表面接触角$\theta $从45°到135°时驱替效率De的取值为0.471到0.620, 增加了31.63%. 同时观察到, 固体表面润湿性固定时, 幂律指数n越大, 驱替效率De越低.
图 14 润湿性$\theta $和幂律指数n对幂律流体驱替效率的影响
Figure14. Effects of contact angles $\theta $ and power-law exponent n on power-law fluid displacement efficiency.

2 -->

4.4.多孔介质几何类型的影响

本小节研究多孔介质几何类型对不混溶幂律流体驱替过程的影响. 在数值模拟中M × N = 240 × 600, S = 6, S₁ = 40, X = 60, Y = 60, 初始动力黏度比M为5.0, Ca = 0.0446, 固体表面都是中性润湿($\theta = {90^{\rm{o}}}$). 障碍物的几何类型分别为正方形、圆形与等边三角形三种情况, 当障碍物的几何类型为正方形时, A = B = 30; 当障碍物的几何类型为圆形时, 障碍物都是半径为16.93的圆形; 当障碍物为等边三角形时, 其边长为49.96. 保证对于所研究的所有的情况孔隙率都为$\xi = 0.78125$.
图15给出了多孔介质几何类型和被驱替液为剪切变稀、牛顿以及剪切变稠三种流体驱替完成时指进形态图. 从图15可以发现, 对于所有的幂率指数n的情况, 都有当多孔介质类型为三角形时指进现象最明显, 因此驱替完成时所需的时间最短. 当n = 0.4时(图15(a)—(c)), 对应多孔介质结构为方形、圆形和三角形的驱替完成时间分别为25.0, 24.5和20.9; 当n = 0.7时(图15(d)—(f)), 对应上述多孔介质结构的驱替完成时间分别为22.1, 22.0和19.0; n = 1.0时(图15(g)—(i)), 对应上述多孔介质结构的驱替完成时间分别为20.7, 20.3和17.7; n = 1.3时(图15(j)—(l)), 对应上述多孔介质结构的驱替完成时间分别为19.8, 19.1和17.1; n = 1.6时(图15(m)—(o)), 对应上述多孔介质结构的驱替完成时间分别为19.3, 18.5和16.9. 从上述数据可以发现相对于其他两种情况, 多孔介质结构为三角形时驱替完成时间明显减少, 这是因为在相同的孔隙率情况下, 相对于其他两种多孔结构, 流体通过三角形结构的多孔介质结构时所对应的通道最小, 因此驱替过程受到的阻力最大. 从图15还可以发现对于相同的多孔介质结构被驱替相的幂律指数n越大, 指进现象越明显, 这是因为被驱替相的幂率指数n越大, 驱替过程的阻力就越大, 因此指进越明显, 对应的驱替完成时间越短.
图 15 不同的障碍物几何类型, 被驱替液为剪切变稀、牛顿与剪切变稠流体时驱得到的指进形态图　(a)?(c) n = 0.4; (d)? (f) n = 0.7; (g)?(i) n = 1.0; (j)?(l) n = 1.3, (m)?(o) n = 1.6
Figure15. Final finger patterns obtained under different geometric type for shear thinning, Newtonian and shear thickening fluids: (a)? (c) n = 0.4; (d)?(f) n = 0.7; (g)?(i) n = 1.0; (j)?(l) n = 1.3; (m)?(o) n = 1.6.

为进一步定性分析障碍物几何类型以及幂律指数n对驱替过程的影响, 图16给出了在多孔介质几何类型不一样时不同幂律指数n下得到的驱替效率. 从图16可以看出对于所研究的多孔介质结构, 驱替效率都随着n的增加而减小, 该结论与图15的指进现象得到对应: n越大指进现象越明显. 这是因为随着幂率指数的增加, 被驱替流体更黏稠, 不利于驱替过程的进行. 从图16还可以看出, 对于所有的幂率指数n的情况, 圆形多孔结构和方形多孔结构得到的驱替效率比较接近(与图15中的驱替时间对应), 而三角形多孔介质结构得到的驱替效率明显减小.
图 16 障碍物几何类型和幂律指数n对幂律流体驱替效率的影响
Figure16. Effects of geometric type and power-law exponent n on power-law fluid displacement efficiency.

基于HCZ两相流模型, 建立了一个不可压非牛顿幂律流体两相流模型, 并用该模型研究了在含多孔介质的微通道内牛顿气体驱替非牛顿幂律流体液体的流动机理, 主要分析了非牛顿流体的幂律指数n、Ca数、初始动力黏度比M、多孔介质表面润湿性$\theta $以及多孔介质障碍物几何类型对驱替过程的影响, 得出以下结论:
1) 幂律指数n越大, 驱替过程的指进现象越明显, 驱替效率越低, 且幂率指数n越大, 其驱替过程受Ca数的影响, 黏性比的影响越小;
2)对于剪切变稀, 牛顿与剪切变稠流体, 都有随Ca数和黏性比的增加指进现象越明显, 驱替效率越低;
3)固体表面接触角越小, 其表面残留的液体就越多, 驱替效率越低;
4)在孔隙率$\xi $相同的情况下, 相对于多孔介质内障碍物为圆形和方形的情况, 障碍物为等边三角形时, 指进现象最明显, 驱替效率最低.