STUDY ON THE EFFECT OF DYNAMIC INTERFACIAL PROPERTIES OF LIQUID BRIDGES ON SPONTANEOUS LIQUID-LIQUID IMBIBITION
-
摘要: 在自然界和工业应用中, 多孔介质中的多相流动现象普遍存在, 如地下水流动、油气开采、非饱和颗粒材料等. 在这些过程中, 液桥的存在对多相流动产生显著影响. 基于改进的液−液自发渗吸理论模型及改进流−固作用力格式的两组分Shan-Chen模型, 研究了液桥的存在对毛细管内液−液自发渗吸的影响. 结果表明: 与不含液桥的渗吸过程相比, 液桥的存在使得毛细管内产生3个界面, 显著增大了系统的动态接触角, 降低自发渗吸速度. 随着润湿性的增强, 液桥的存在使得界面动态变化特性增强的幅度更大. 当液桥的黏度小于非润湿流体的黏度时, 整个系统的动态接触角随着渗吸长度的增加而增强; 若二者黏度相等, 整个系统的动态接触角为一稳定值; 而当液桥的黏度更大时, 整个系统的动态接触角随着渗吸的进行逐渐减小. 将模拟获得的实时动态接触角纳入至渗吸理论中, 可以改善由液桥界面动态变化导致的渗吸长度的理论值与模拟结果间存在的偏差, 并得到与模拟基本一致的结果. 文章还利用模拟数据, 系统地评价了Cox-Voinov模型对含有液桥以及不含液桥体系的渗吸过程中的界面动态变化特性及渗吸长度预测能力.
-
关键词:
- 液−液自发渗吸 /
- 格子Boltzmann方法 /
- 液桥 /
- 动态接触角 /
- 两相流
Abstract: The phenomenon of multiphase flow in porous media occurs widely in both natural and industrial applications, including groundwater flow, oil and gas development, and unsaturated particulate materials. Liquid bridges play a crucial role in these processes and significantly influence the dynamics of multiphase flows. In this work, we have investigated the effect of liquid bridges on spontaneous liquid-liquid imbibition in capillaries based on a modified spontaneous liquid-liquid imbibition theory model and a two-component Shan-Chen model with an improved fluid-solid interaction force format. The results indicate that the presence of liquid bridges leads to the formation of three interfaces within the capillary, which significantly increases the dynamic contact angles of the entire system and reduces the rate of spontaneous imbibition compared to processes without liquid bridges. Additionally, with the increase of wettability, the presence of liquid bridges enhances the dynamic interface change characteristics more strongly. When the viscosity of the liquid bridge is lower than that of the non-wetting fluids, an increase in the imbibition length will enhance the dynamic contact angles of the entire system; if they both have equal-sized viscosities, then the dynamic contact angles of the entire system remain a constant; and when the viscosity of the liquid bridge is larger than that of the non-wetting fluids, there is a gradual decrease in the dynamic contact angle during imbibition. By incorporating the real-time dynamic contact angle obtained from the simulation into the imbibition theory, the deviation between the theoretical and simulated imbibition lengths due to the dynamic changes of the liquid-bridge interface can be improved, and the results are consistent with the simulation. In this study, we also systematically evaluated the predictive ability of Cox-Voinov models for interface dynamic variation and imbibition length in imbibition processes with and without liquid bridges by using the data from the simulation. -
引 言
多孔介质中的液−液两相自发渗吸现象是一个由表面张力和润湿性主导的两相流动过程[1-3]. 该现象普遍发生在石油开采过程中的压裂液在页岩/致密油储层中的渗吸置换[4]、水驱/注水吞吐提高石油采收率[5-6]、地下水流动[7]以及非饱和颗粒材料[8]等领域. 在上述过程中, 由于多孔介质的复杂结构和润湿性的影响, 角流作用显著[9], 并且界面卡断现象频繁发生[10], 这会导致在孔隙内部形成液桥或者气泡[11]. 液桥的存在进一步使得液−液自发渗吸过程中产生多个界面, 致使毛细管内的力学机制也更加复杂 [12-15]. 因此, 捕捉液桥的界面动态特性并分析其对液−液自发渗吸过程的影响, 对于理解水驱提高石油采收率的微观作用机制具有重要意义.
为了深入探究界面动态特性对液−液自发渗吸过程的影响, 研究人员采用了理论、实验和数值模拟等多种方法对该过程进行了广泛的研究. 在理论研究方面, Lucas[16]和Washburn[17]最早基于泊肃叶流动定律, 提出了经典的圆管自发渗吸理论(LW方程), 该理论假设只存在一个界面且不考虑液桥的存在, 认为界面的上升高度与时间的平方根成正比. 然而, 原始的LW方程并未考虑非润湿相黏度的贡献, 因此仅适用于气−液自发渗吸的预测. 为了弥补这一缺陷, Mumley等[18]在1986年利用实验手段, 系统研究了甘油−水体系在预润湿和干燥毛细管内的渗吸过程. 在此基础上, 通过考虑非润湿流体的黏滞力的贡献, 修正了原始的LW方程, 使之能够适用于预测液−液体系的渗吸过程. 此外, 实验中还观察到两相界面呈现出动态变化特征, 接触角也随之保持动态变化. 随后, Fermigier等[19]利用实验分别测量了硅油−空气和硅油−甘油体系渗吸过程中的界面的动态接触角. 结果表明, 液−液体系的界面的动态变化特性更显著. Walls等[20]基于实验和理论分析研究了丙三醇−水体系的自发渗吸过程, 并认为忽略动态接触角影响的LW方程所预测的弯液面位置与实验数据存在较大差距. 至于毛细管内存在润湿流体液桥时的相关实验报道则相对较少. 只有André等[21]利用实验研究了水平毛细管内预先存在一个乙醇液桥时的乙醇−硅油体系的渗吸过程. 从他们的实验结果可以看出, 液桥的存在使得渗吸过程中出现了3个界面, 显著增强了动态接触角的影响. 尽管如此, 修正毛管力的LW方程仍然可以预测存在液桥时的液−液自发渗吸过程.
相比于实验, 数值模拟方法能够帮助捕捉渗吸过程中的更多细节. 其中, 格子Boltzmann方法(LBM)源于对玻尔兹曼输运方程简化形式的离散求解, 玻尔兹曼方程可以在克努森数Kn$\ll $1的条件下恢复至宏观的Navier–Stokes (NS)方程[22]. 其可以自动考虑惯性力, 入口通量(初始时刻的不完全发展流动阶段)以及捕捉界面的动态变化特性(实时动态接触角)等[23-24]. 因此, LBM十分适用于模拟两组分流体(例如油−水)的自发渗吸过程. 在4类多相流LBM模型中, Shan-Chen模型因其简便性及明确的物理背景, 被广泛应用至气−液[25]和液−液自发渗吸[26]的模拟中. 在液−液自发渗吸方面, Chibbaro等[27]基于原始的两组分Shan-Chen模型模拟了等密度比和黏度比的两相流体自发渗吸过程, 并捕捉了界面的实时动态接触角, 结果表明, 利用模拟获得的动态接触角数据可以修正LW方程, 修正后的LW方程所预测的渗吸长度与模拟结果一致. Budaraju等[28]采用原始的两组分Shan-Chen模型, 研究了具有轴向变化的毛细管内液−液渗吸过程及其动态润湿效应. 然而, 上述的模拟和理论研究大部分局限在连续流动的流动类型中. 但是有实验结果表明, 在微通道内存在液桥现象[29]. 这将导致两相界面分布更加复杂, 界面的接触角的动态变化更加显著.
为了深入研究液桥对毛细管内自发渗吸的影响, 本文首先探讨了毛细管内存在液桥时的界面特性, 并对液−液自发渗吸的控制方程进行了扩展. 接着, 本文建立了一个基于中心矩碰撞算子的两组分Shan-Chen模型, 并采用改进虚拟密度的流−固作用力格式精确表征了流体的润湿行为. 此外, 本文还通过理论和LBM模拟相结合的方法, 重点研究了液桥的存在对毛细管内液−液自发渗吸过程的影响, 并详细捕捉了液桥的界面动态特性. 同时, 本文也探讨了润湿性和黏度比对毛细管内液−液自发渗吸过程的影响. 最后, 本文通过采用模拟得出的数据, 对表征液桥的界面动态特性的经典水动力学模型进行了全面的分析和扩展. 本研究的结果可以为含水桥的多孔介质中的复杂油−水两相流动的准确表征和预测提供理论和模拟基础.
1. 模型及方法
1.1 数学模型
如图1所示, 对于二维的水平毛细管内存在一个长为l1的润湿相流体的液桥, 其渗吸过程的力学平衡关系可以表示为[30]
![]()
图 1 毛细管内存在液桥时的液−液自发渗吸示意图Figure 1. Schematic diagram of spontaneous liquid-liquid imbibition in a capillary with a liquid bridge$$ \frac{{{\text{d}}(mU)}}{{{\text{d}}t}} = {F_{{\rm{drive}}}} - {F_{{\rm{drag}}}} $$ (1) 式中, U为界面的移动速度, m为毛细管内液柱的质量, t为时间; 由于界面各处的速度并不一致, 本文中采用界面中心线处的移动速度近似代替, 即U = dl(t)/dt, l(t)指左侧润湿流体的侵入长度. 据此, 方程(1)的左边可以得到简化, 表示为
$$ \begin{split}& \frac{{{\text{d}}(mU)}}{{{\text{d}}t}} = m\frac{{{\text{d}}U}}{{{\text{d}}t}} + U\frac{{{\text{d}}m}}{{{\text{d}}t}}= \\& A\left\{ {[{\rho _w}(l + {l_1}) + {\rho _{nw}}(L - l - {l_1})]\frac{{{{\text{d}}^2}l}}{{{\text{d}}{t^2}}} + ({\rho _w} - {\rho _{nw}}){{\left(\frac{{{\text{d}}l}}{{{\text{d}}t}}\right)}^2}} \right\} \\ \end{split} $$ (2) 式中A = bH表示毛细管的侧面积, b为毛细管深度; H和L分别为微毛细管的宽度和长度, ρw和ρnw为润湿流体(蓝色)和非润湿流体(红色)的密度. Fdrive和Fdrag分别为渗吸的动力和阻力, 其中Fdrive以毛管力为主, Fdrag以黏滞阻力为主. 基于充分发展流动的假设, Fdrive和Fdrag分别表示为
$$ {F_{{\rm{drive}}}} = 2\sigma b(\cos {\theta _1} - \cos {\theta _2} + \cos {\theta _3}) $$ (3) $$ {F_{{\rm{drag}}}} = A\frac{{12[{\mu _w}(l + {l_1}) + {\mu _{nw}}(L - l - {l_1})]}}{{{H^2}}}\frac{{{\text{d}}l}}{{{\text{d}}t}}{\text{ }} $$ (4) 式中, σ表示界面张力; μw和μnw为润湿流体和非润湿流体的动力黏度; θ1, θ2和θ3分别为润湿流体左侧、中间及右侧界面的接触角, 小于90°. 在渗吸过程中左侧和右侧界面产生的毛管力与流动方向相同, 表现为动力作用;而中间界面产生毛管力与流动方向相反, 表现为阻力作用. 结合方程(1)的具体形式, 当毛细管内存在一个液桥时, 其自发渗吸的控制方程可以表示为
$$ \begin{split}& \frac{{2\sigma (\cos {\theta _1} - \cos {\theta _2} + \cos {\theta _3})}}{H} = \\&\qquad \frac{{12[{\mu _w}(l + {l_1}) + {\mu _{nw}}(L - l - {l_1})]}}{{{H^2}}}\frac{{{\text{d}}l}}{{{\text{d}}t}}+{\text{ }} \\&\qquad [{\rho _w}(l + {l_1}) + {\rho _{nw}}(L - l - {l_1})]\frac{{{{\text{d}}^2}l}}{{{\text{d}}{t^2}}} + \\ &\qquad ({\rho _w} - {\rho _{nw}}){\left(\frac{{{\text{d}}l}}{{{\text{d}}t}}\right)^2} \end{split} $$ (5) 方程(5)是典型的2阶非线性常微分方程, 本文采用4阶显式的Runge-Kutta方法对其求解, 其两个特定的边值条件为: t = 0, l = 0; t = 0, dl/dt = 0. 需要指出的是, 若不考虑动态接触角的影响, 取θ1 = θ2 = θ3 = θ0 (θ0表示润湿流体在壁面上的静态接触角).
1.2 两组分Shan-Chen模型
本文使用基于中心矩碰撞算子[31]的两组分Shan-Chen模型, 其碰撞过程可以表示为[32]
$$ \begin{split}& {f_{\alpha ,}}_i({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t + \Delta t) - {f_{\alpha ,i}}({\boldsymbol{x}},t) = \\ &\qquad {{\boldsymbol{M}}^{ - 1}}{{\boldsymbol{N}}^{ - 1}}{\boldsymbol{SNM}}[f_{\alpha ,i}^{eq}({\boldsymbol{x}},t) - {f_{\alpha ,i}}({\boldsymbol{x}},t)] + \\ &\qquad {{\boldsymbol{M}}^{ - 1}}{{\boldsymbol{N}}^{ - 1}}({\boldsymbol{I}} - {\boldsymbol{S}}/2){\boldsymbol{NM}}\left| {{R_{\alpha ,}}_i} \right\rangle \end{split} $$ (6) 式中, fα 表示流体组分α的密度分布函数, 在本文中, 下标α = 1和α = 2分别表示非润湿流体和润湿流体; ei为离散速度矢量, 对于D2Q9模型 (i = 0,1,2,···,8), 其可以表示为
$$\left. \begin{gathered} \left| {{{\boldsymbol{e}}_{ix}}} \right\rangle = {\left[ {0,1,0, - 1,0,1, - 1, - 1,1} \right]^{\text{T}}} \\ \left| {{{\boldsymbol{e}}_{iy}}} \right\rangle = {\left[ {0,0,1,0, - 1,1,1, - 1, - 1} \right]^{\text{T}}} \\ \end{gathered}\right\} $$ (7) 在方程(6)中, M和N为转换矩阵, 利用M和N可以将分布函数fα转换至对应的中心矩, 其详细的形式可以在文献[32]中找到. 同时, 利用M和N也可以将平衡分布函数$ f_{\alpha ,i}^{eq} $和外力项Rα转换至其对应的中心矩$\bar {T}_{\alpha ,i}^{eq}$和Cα, 具体表示为
$$ \left.\begin{gathered} \left| {\bar {{\boldsymbol{T}}}_{\alpha ,i}^{eq}} \right\rangle = {\boldsymbol{NM}}f_{\alpha ,i}^{eq} = \left[ {{\rho _\alpha },0,0,2{\rho _\alpha }c_s^2,0,0,0,0,{\rho _\alpha }c_s^4} \right] \\ \left| {{{\boldsymbol{C}}_{\alpha ,i}}} \right\rangle = {\boldsymbol{NM}}\left| {{R_{\alpha ,i}}} \right\rangle = \left[ {0,{F_{\alpha ,x}},{F_{\alpha ,y}},0,0,0,{F_{\alpha ,y}}c_s^2,{F_{\alpha ,x}}c_s^2,0} \right] \\ \end{gathered}\right\} $$ (8) 式中, cs为格子声速, ρα表示α流体组分的密度. 在Shan-Chen模型中, 密度ρα及组分的共同速度u 可以表示为 [33]
$$ {\rho _\alpha }{\text{ = }}\displaystyle\sum\limits_i {{f_{\alpha ,i}}} {\text{ }},\qquad {\boldsymbol{u}} = \frac{{\displaystyle\sum\limits_\alpha {\left( {\displaystyle\sum\limits_i {{{\boldsymbol{e}}_i}{f_{\alpha ,i}} + 0.5\Delta t{{\boldsymbol{F}}_\alpha }} } \right)} }}{{\displaystyle\sum\limits_\alpha {\displaystyle\sum\limits_i {{f_{\alpha ,i}}} } }} $$ (9) 在方程(6)中, S为对角矩阵, 可以表示为
$$ {\boldsymbol{S}} = {\rm{diag}}({s_0},{s_1},{s_1},{s_b},{s_2},{s_2},{s_3},{s_4},{s_4}) $$ (10) 由于质量守恒, 通常设置s0 = s1 = s4 = sb = 1.0, s3 = (16−8s2)/(8−s2); s2 = 1/τ为松弛因子, 其与流体组分的运动黏度相关(τα = 3να + 0.5). 为了使其在界面处光滑地过渡, 采用下述方法定义系统的运动黏度νs [34]
$$ \frac{1}{{{\nu _s}}} = \frac{{{\rho _1}}}{{{\rho _1} + {\rho _2}}}\frac{1}{{{\nu _1}}} + \frac{{{\rho _2}}}{{{\rho _1} + {\rho _2}}}\frac{1}{{{\nu _2}}} $$ (11) 在方程(9)中, Fα = [Fα, x, Fα, y]表示α流体组分所受的合力, 包括流体−流体作用力Ff-f, 流−固作用力Ff-s以及其他的体积力Fb. 对于两组分Shan-Chen模型, 流体−流体作用力Ff-f可以表示为 [33]
$$ {\boldsymbol{F}}_\alpha ^{f - f} = - {G_{\alpha \bar \alpha }}{\psi _\alpha }({\boldsymbol{x}})\sum\limits_i {{w_i}} {\psi _{\bar \alpha }}({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t){{\boldsymbol{e}}_i} $$ (12) 式中, ψα表示伪势, 本文取ψα = ρα; G12和G21为控制流体组分间排斥力大小的参数, 一般设置G12 = G21 > 0, 值得一提的是, G12和G21主要影响两相界面厚度及界面张力大小; wi为权重系数, 对于D2Q9格子模型(i = 0 ~ 9)具体为: |ei|2 = 0, wi = 4/9; |ei|2 = 1, wi = 1/9; |ei|2 = 2, wi = 1/36; 对于流−固作用力Ff-s, 其是表征润湿性的关键参数. 但是, 传统的基于虚拟密度的流−固作用力格式会导致壁面处存在一层非物理的虚拟密度[35]. 为此, 本文采用改进虚拟密度格式的流−固作用力抑制此类非物理现象的出现, 具体表示为[36]
$$ {\boldsymbol{F}}_\alpha ^{f - s} = - {G_{\alpha \bar \alpha }}{\psi _\alpha }({\boldsymbol{x}})\sum\limits_i {{w_i}} \phi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t){\bar \psi _{\bar \alpha }}({\boldsymbol{x}}){{\boldsymbol{e}}_i} $$ (13) 式中ϕ(x + ei)为标识函数, 对于流体节点, 其值为1; 对于固体节点, 其值为0; $ {\bar \psi _\alpha }(x) $表示固体边界层的虚拟的伪势, 是控制流−固作用强度的因素之一. 其定义为固体节点周围所有流体节点的伪势的平均值, 具体表示为[36]
$$ {\bar \psi _\alpha }({\boldsymbol{x}}) = {\chi _\alpha }\frac{{\displaystyle\sum\limits_i {{w_i}{\psi _\alpha }({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t)[1 - \phi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t)]} }}{{\displaystyle\sum\limits_i {{w_i}[1 - \phi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t)]} }} $$ (14) 式中, χα是控制润湿性的参数, 本文设置χ1 = 1 + ξ, χ2 = 1−ξ; wi为格子权重系数; 为了使得计算的平均虚拟密度更加准确, 本文采用D2Q25 (i = 0 ~ 24)的格子权重系数, 具体为: |ei|2 = 0, wi = 247/420; |ei|2 = 1, wi = 4/63; |ei|2 = 2, wi = 4/135; |ei|2 = 4, wi = 1/180; |ei|2 = 5, wi = 2/945; |ei|2 = 8, wi = 1/15120. 此外, 在Shan-Chen模型中, 当系统处于稳定状态时, 其总压力Ptotal可以表示为
$$ {P_{{\text{total}}}} = c_s^2({\rho _1} + {\rho _2}) + \frac{1}{2}c_s^2({G_{12}} + {G_{21}}){\psi _1}{\psi _2} $$ (15) 2. 模型测试
2.1 界面张力测试
本节中, 开展了LBM方法中最基本的Laplace定律测试. 其主要目的是为了获取两组分流体间的界面张力σ. 具体而言, 在Shan-Chen模型中, 通过模拟不同大小的静态圆形液滴, 当模拟达到稳定时, 液滴内部和外部的压差∆P与液滴半径R和界面张力σ满足Laplace定律, 即∆P = Pint−Pout = σ/R. 通过模拟一系列不同大小的静态圆形液滴, 可以获得一系列的压差和液滴半径的对应关系, 进而通过线性拟合确定界面张力大小[37]. 为了确定合理的体相密度ρbulk和溶解密度ρdis, 首先进行了混溶性测试, 确定了流体间作用强度G12 = G21 = 4.0 (ρbulk = 1.0, ρdis = 0.013)作为后续的模拟参数. 然后, 建立了一个101 × 101的格子区域, 四周设置为周期边界条件, 并在区域的中心初始化一个半径为R0的圆形液滴. 此外, 为了说明界面张力与黏度的无关性, 采用3种不同的黏度比M (M = μnw/μw)进行测试, 分别为 (τ1 = 1.0, τ2 = 0.55, M = 10.0; τ1 = 1.0, τ2 = 1.0, M = 1; τ1 = 0.55, τ2 = 1.0, M = 0.1). 模拟结果如图2所示. 可以看出, 对于不同半径的液滴, 其压差∆P与半径的倒数1/ R0呈现良好的线性关系, 通过拟合确定该条件下两组分的界面张力为σ = 0.088. 此外, 还可以观察得到不同黏度的液滴, 其界面张力基本一致, 说明本文模型中的界面张力与黏度无关.
2.2 静态接触角测试
上节中的Laplace定律测试主要针对流体−流体作用力的测试. 在本小节, 在此基础上考虑流−固作用力, 测试改进虚拟密度的流−固作用力格式模拟接触角的能力. 首先, 建立201 × 101的格子区域, 在区域的顶部和底部设置固体壁面, 壁面处采用反弹边界, 区域左右设置为周期边界. 模拟参数与之前保持一致. 在初始时刻, 我们在壁面上初始化一个半圆形液滴, 设置不同的接触角调节系数ξ. 具体而言, 调节系数ξ和流体−流体作用力强度的大小的乘积((1 + ξ)G21和(1 − ξ)G12), 表示壁面对于流体1和流体2的黏附力的强度. 并且, 其决定了固−液界面张力的大小(σs_1和σs_2). 根据Young氏方程可知, 接触角的余弦值cosθ可以表示为
$$ \cos \theta = \frac{{{\sigma _{{\rm{s}}\_1}} - {\sigma _{{\rm{s}}\_2}}}}{\sigma } $$ (16) 式中, σs_1和σs_2分别表示流体1和流体2与壁面的界面张力大小. 模拟稳定后, 如图3所示, 接触角θ的大小由不同的ξ控制[38]. 当ξ = 0时, 壁面对两种流体的作用力的强度相同(G12 = G21), 即, σs_1 = σs_2, 故两种流体对壁面的黏附力相同, 润湿角都应为90°; 而当ξ > 0时, (1 + ξ)G21 > (1 − ξ)G12, σs_1 > σs_2, 因此流体1润湿壁面, 润湿角小于90°, 而流体2的接触角大于90°; 相反, 当ξ < 0时, (1 + ξ)G21 < (1 − ξ)G12, σs_1 < σs_2, 壁面对流体2的黏附力较大, 因此流体2润湿壁面, 润湿角小于90°, 而流体1的接触角大于90°.
根据图3可以看出, 本文所提出的改进虚拟密度的流−固作用力方案有效地抑制了壁面上的虚拟密度层. 进一步模拟了不同黏度的液滴在不同的调节系数ξ下所对应的接触角θ. 在我们的模拟结果中, 当ξ < 0.8时,不同黏度的两种流体的接触角与黏度的大小基本无关, 但是当ξ ≥ 0.8时, 黏度比的改变导致了人工虚假速度显著增加, 致使流体在壁面上发生了局部的非物理滑动, 不同黏度比条件下的接触角呈现一定的偏差. Sedahmed等[37]在模拟中也观察到了类似的现象. 但是, 对于本文后续的模拟来说, 这种偏差是可以忽略不计的. 因为我们在模拟中对ξ的取值都是负值, 并未采用ξ ≥ 0.8时的模拟数据. 进一步, 通过拟合发现ξ-θ满足3次多项式关系, 即θ = −26.33ξ3−3.13ξ2−65.18ξ + 89.16.
此外, 由于LBM中的单位为格子单位. 因此, 需要将格子单位转化为实际的物理单位, 具体转化系数如表1所示. 需要注意的是, 为了方便描述, 本文在后续的研究中只给出了格子单位.
表 1 模型参数与实际参数的单位转换Table 1. Unit conversion of model and real parametersParameters Lattice values Physical values Scale factor density, ρw 1.0 1000 kg/m3 1000 kinematic viscosity, νw 0.1667 5 × 10−6 m2/s 3 × 10−5 lattice length, Δx 1.0 2.63 × 10−6 m 2.63 × 10−6 interface tension, σ 0.088 3 × 10−2 N/m 0.341 3. 结果与讨论
3.1 数值方法验证
在详细分析模拟结果之前, 评估本文采用的LBM方法的准确性至关重要. 然而, 目前关于毛细管内存在液桥时的液−液渗吸实验数据较为稀缺. 因此, 本文采用André等[21]开展的毛细管内液−液体系的连续渗吸实验数据作为标准数据, 验证本文的数值方法的准确性. André等 [21]在实验中以乙醇(0.839 g/cm3)和硅油(0.96 g/cm3)分别作为润湿和非润湿流体, 两种流体的黏度比约为M = 34.0. 由于实验以三维的圆管为研究对象, 而本文的模拟为二维数值模型. 为了方便对比实验和模拟结果, 采用无量纲的渗吸长度L*和无量纲的渗吸时间t*作为标度归一化实验和模拟结果. 其中, L*和t*可以表示为
$$ \left.\begin{split} & {L^*} = \frac{{{\mu _w} - {\mu _{nw}}}}{{{\mu _{nw}}}}\frac{{l(t)}}{L} \\ &{t^*} = \frac{{({\mu _w} - {\mu _{nw}})\sigma H\cos {\theta _0}}}{{\kappa \mu _{nw}^2{L^2}}}t \end{split}\right\} $$ (17) 式中, κ为转换系数, 对于三维圆管, κ = 4; 对于二维的模拟κ = 3. 为了保证模拟与实验条件基本一致, 首先建立了601 × 101的格子空间, 在150 < x < 450位置处设置长度为L = 300. 通道壁面接触角设为θ0 = 45° (ξ = −0.615), 并采用半步长反弹格式, 其余边界均采用周期边界. 在x≤150和x≥550的位置处初始化为润湿相流体, 两种流体的黏度比设置为M = 34.0 (τ1 = 2.2, τ2 = 0.55; 与实验中的黏度比保持一致). 微毛细管的长度和宽度比值为15 (与实验长宽比相同). 图4给出了模拟的渗吸过程与实验结果的对比, 需要指出的是, 由于润湿流体的黏度小于非润湿流体的黏度, 归一化后的渗吸长度的范围为: −1 ~ 0. 根据图4可以看出, 本文采用的LBM模型与实验结果所得渗吸长度的吻合性较好, 插图内t* = −0.5时的两相流动过程的快照也基本一致.
![]()
3.2 液桥长度的影响
本节通过LBM模拟研究液桥的存在对毛细管内的液−液自发渗吸过程的影响. 首先, 建立601 × 101的格子空间, 并在150 < x < 450位置处设置长度为L = 300, 宽度为H = 30的通道, 将通道壁面接触角设为θ0 = 45°, 并采用半步长反弹格式, 其余边界均采用周期边界. 在x≤150和x≥500的位置处初始化为润湿相流体. 为了突出液桥对渗吸过程的影响, 本节同时开展两组模拟. 在第1组模拟中设置液桥的长度l1 = 0 (不含液桥). 在第2组模拟中设置液桥的长度为l1 = 50 (200≤x≤250). 其余位置均设置为非润湿相流体, 两相流体的黏度比设置为M = 1. 图5对比了液桥长度分别为l1 = 0和l1 = 50在不同时刻下的渗吸过程. 可以看出, 对于l1 = 0的情况, 润湿流体(蓝色)与非润湿流体(红色)之间形成一个界面(弯液面); 相反, 对于l1 = 50的情况, 则形成3个弯液面, 其中左右两侧的弯液面在渗吸过程中形成的毛管力起到动力作用, 而中间界面形成的毛管力由于与运动方向相反, 表现为阻力作用. 此外, 相比于不含液桥的毛细管内的渗吸过程, 在相同时刻下, 存在液桥的毛细管中的左侧弯液面的推进速度明显降低, 说明液桥的存在不利于自发渗吸的进行. 这种现象在之前的气−液自发渗吸实验的研究中也有相同的报道. 例如Thamdrup等[39]通过实验研究表明, 气泡的产生会显著降低水的毛细流动速度.
![]()
图 5 毛细管内存在不同长度液桥的液−液渗吸过程Figure 5. Liquid-liquid imbibition process in capillaries with different liquid bridge lengths为了深入研究液桥的存在导致渗吸速度下降的原因, 本文在模拟中精确捕捉了两相界面及其动态接触角, 如图6所示. 图6给出了两种情况下的接触角的等效余弦值cosθ*(cosθ* = cosθ1−cosθ2 + cosθ3)随时间的演化关系(注意, 当l1 = 0时cosθ* = cosθ1). 值得一提的是, 在模拟中, 毛细管的宽度以及界面张力是不变的. 因此, cosθ*的数值大小直接决定了自发渗吸过程中的毛管力的大小. 具体而言, 对于不存在液桥的自发渗吸过程, 在t = 0时, 由于两相弯液面并未形成, 此时所有的两相界面均为平直界面. 因此在模拟中, 初始时刻的接触角表现为90° (cosθ* = 0). 随后, 弯液面快速形成并向前推进, 在此期间, 惯性力作用显著, 导致接触角存在较大波动; 伴随着渗吸距离的增加, 惯性力作用减弱, 毛管力和黏滞力占据主导作用, 渗吸速度趋于稳定, 动态接触角也趋向一个稳定值 (cosθ* = 0.6), 并且其小于静态接触角的余弦值 (cosθ0 = 0.707). 而当毛细管内存在液桥时, 在初始时刻, 其动态接触角的演化与不含液桥的情况类似. 但是, 由于液桥的存在导致其形成了3个界面, 其在稳定渗吸的过程中等效动态接触角cosθ*显著小于l1 = 0的情况 (cosθ* = 0.46). 这说明液桥的存在显著增强了系统的动态接触角效应, 导致系统渗吸时的毛管力(渗吸动力)减小, 从而抑制了渗吸速度. 当液桥的右侧界面到达出口端时, 液桥右侧弯液面趋于平直界面, 其余两个液面趋向于静态接触角, cosθ*下降至0, 渗吸停止.
![]()
图 6 渗吸过程中毛细管内存在不同长度的液桥时的动态接触角的时间演化Figure 6. Time evolution of the dynamic contact angle during imbibition for the presence of liquid bridges with different lengths in a capillary在此基础上, 我们也捕捉了左侧液体的渗吸长度随时间的演化关系, 如图7所示. 需要指出的是, 根据理论公式(5), 在不考虑液桥界面的动态变化特性(动态接触角)的情况下, 系统的静态接触角cosθ* = cosθ0. 并且由于在模拟中设置两种流体的黏度比为M = 1, 密度比也为1. 因此, 即使存在液桥, 其不考虑动态接触角的理论渗吸长度也与不含液桥的理论渗吸长度一致, 均为图7中的黑色实线. 在不考虑渗吸过程中的界面的动态变化特性时, 方程(5)预测的渗吸长度与模拟结果存在一定偏差, 尤其当液桥存在时, 偏差更为明显. 基于我们之前关于气液自发渗吸的研究思路[25], 我们尝试将图6中系统界面的实时动态接触角数据cosθ*代入方程(5)中, 并重新计算了渗吸长度. 可以看出, 在校正动态接触角后, 对于存在液桥和不存在液桥的渗吸过程, 理论预测的渗吸长度与LBM结果的偏差都得到了很大程度的改善. 这也说明, 液桥的界面动态变化特性是导致理论模型与LBM存在偏差的主要原因.
![]()
图 7 渗吸过程中存在不同液桥长度时的理论渗吸长度与LBM的比较Figure 7. Comparison of theoretical imbibition length with LBM during imbibition with different liquid bridge lengths然而, 本文建立的数学模型忽略了早期的不稳定渗吸阶段(流动发展阶段). 为了更加全面地评价校正实时动态接触角之后的方程(5)对渗吸过程预测的准确程度, 我们在图7中也给出了t ≤ 10000时间段内的LBM结果以及理论模型的局部放大图. 同时, 本文也提供了不同时刻 (t = 100, t = 1000, t = 5000, t = 10000 )的毛细管入口的速度剖面的分布(左右两个端点为壁面节点, 故速度为0), 用以判断流动是否满足充分发展的假设. 根据图7中的两个子图可以看出, 在初始时刻, LBM结果与校正动态接触角的理论值存在一定差距. 具体而言, (1) 当t = 0时, 两相界面为平直界面, 弯液面并未形成, 界面的动态接触角θd ≈ 90°, 此时渗吸的驱动力为Fdrive = 2σR(cosθ0−cosθd). 在较大的驱动力作用下, 两相界面的接触线快速向前移动, 而界面中心区域移动速度较小甚至小于零, 导致LBM模拟的渗吸长度在t ≤1000表现为负值. 而理论模型由于不能考虑弯液面的形成以及外部流体的不稳定吸入过程, 导致其预测的渗吸速度更大, 高估了渗吸长度; (2) 在1000≤ t ≤5000时间段, 伴随着渗吸时间的增加, 渗吸长度逐渐增大, 主导渗吸的惯性力作用减弱, 黏滞力的影响也逐渐增大, 入口的流动速度剖面逐渐收敛至均匀Poiseuille速度剖面. 因此, 修正动态接触角的理论模型预测的渗吸长度逐渐接近, 在t = 5000时, 二者的差距仅为2个格子长度; (3) 当5000≤ t ≤10000时, 根据图7提供的入口速度的剖面图可以确定此阶段已经完全达到流动的充分发展阶段(子图中的红色线与绿色线基本重合), 本文的理论模型的假设基本成立, 修正3个动态接触角之后的理论模型预测的渗吸长度与LBM也基本一致.
综上所述, 虽然在初始阶段流动发展不够充分, 导致理论模型与模拟结果之间存在一定差别, 但随着流动的发展 (t ≥ 5000后), 我们的理论模型在考虑界面的动态接触角之后能够较好地预测渗吸过程. 因此, 笔者认为这种差别对整个渗吸过程的理论模型的整体准确性的影响相对较小.
3.3 润湿性的影响
在之前的模拟中, 仅限于讨论壁面静态接触角为θ0 = 45°的情况. 为了更加全面地研究润湿性对于渗吸过程的影响, 本节设置了3种不同的壁面润湿性, 液桥的静态接触角分别为θ0 = 40°, θ0 = 50°以及θ0 = 60°. 其余的模拟参数均与3.2节中的保持一致. 图8给出了不同润湿性条件下, 存在液桥时的实时动态接触角随时间的演化关系以及左侧渗吸长度随时间的演化关系. 为了对比液桥的影响, 我们也给出了各自条件下对应的不含液桥的渗吸长度和动态接触角的时间演化过程. 具体来说, 图8(a) ~ 图8(c)展示的动态接触角和渗吸长度的演化过程具有一定的相似性. 即当毛细管内存在液桥时, 其渗吸长度随时间的演化速度均小于不存在液桥时的情况, 稳定渗吸时系统的动态接触角的余弦值cosθ*也减小. 不同的是, 随着壁面的润湿性增强, 在稳定渗吸时, 存在液桥时的动态接触角与不存在液桥时的动态接触角的余弦值的差值不同. 具体差值为: θ0 = 40°, ∆cosθ* = 0.145; θ0 = 50°, ∆cosθ* = 0.125; θ0 = 60°, ∆cosθ* = 0.104. 上述定量化的结果说明, 随着壁面的润湿性增强, 液桥的界面动态变化特性越强, 对渗吸速度的不利影响也越大. 而对于渗吸长度的理论预测结果与LBM模拟结果的对比也表明, 当不考虑渗吸过程中的界面动态变化特性时, 理论模型预测的渗吸长度与模拟结果存在较大的偏差, 尤其当毛细管内存在液桥时. 按照3.2节中的思路, 通过将模拟得到的实时动态接触角数据(cos θ*)直接应用至方程(5)中, 修正后的方程(5)可以得到与LBM结果基本一致的渗吸长度.
![]()
图 8 不同润湿性条件下的动态接触角的时间演化及理论渗吸长度与LBM对比Figure 8. Time evolution of dynamic contact angle and comparison of theoretical imbibition length with LBM under different wettability3.4 黏度比的影响
在上一节的模拟中, 仅限于讨论了黏度比为M = 1的情况. 为了更全面地研究黏度比对渗吸过程的影响, 本节设置了3种不同的黏度比, 分别为: M = 5.0 (τ1 = 0.75, τ2 = 0.55), M = 1.0 (τ1 = 0.75, τ2 = 0.75), M = 0.2 (τ1 = 0.75, τ2 = 1.75). 壁面的接触角设置为θ0 = 45°, 其余模拟参数与3.2节中的一致. 图9(a) ~ 图9(c)给出了不同润湿流体黏度条件下, 存在液桥时的渗吸过程中的实时动态接触角随时间以及左侧渗吸长度随时间的演化关系. 同样地, 我们也给出了各自黏度比条件下对应的不含液桥的渗吸长度和动态接触角的时间演化过程. 总体而言, 对于固定的非润湿流体黏度, 渗吸速度随着黏度比的增加而减小. 并且, 渗吸长度和界面的动态接触角的时间演化过程存在较大差别. 具体来说, (1) 当黏度比M > 1.0时 (润湿流体的黏度显著小于非润湿流体的黏度), 如图9(a)所示. 渗吸过程表现为加速动力学, 即渗吸长度随着流体吸入量的增加而增大. 其原因在于随着渗吸时间的增加, 越来越多的低黏度的流体进入毛细管, 导致整个毛细管内的黏滞阻力减小. 根据Cox-Voinov公式[40-41]可以看出, 随着渗吸速度的增加, 动态接触角与静态接触角的差距增大, 动态接触角也就越大, 系统的接触角的余弦值cosθ*表现为逐渐减小的趋势. 而当毛细管内含有液桥时(液桥为具有更小黏度的润湿流体), 其渗吸速度增加得更多, 因此其动态接触角的余弦值cosθ*下降得越快. (2) 当黏度比M = 1.0时(润湿流体的黏度等于非润湿流体的黏度), 如图9(b)所示. 可以看出, 渗吸长度与时间呈现出近似的线性关系, 说明渗吸过程中的黏滞阻力保持为一个恒定值, 故渗吸速度也表现为一个定值. 根据Cox-Voinov公式来看, 当渗吸速度基本不变时, 动态接触角基本保持在一个定值. 因此, 两种情况下的动态接触角维持在一个定值, 它们的余弦值cosθ*也基本保持不变. (3) 当黏度比M < 1.0时(润湿流体的黏度大于非润湿流体的黏度), 如图9(c)所示. 可以看出, 自发渗吸过程中渗吸长度随时间的演化速率逐渐减小, 渗吸速度也逐渐减小. 这是由于随着渗吸时间的增加, 高黏度的润湿流体占据了大部分的毛细管空间, 尤其是对于含有液桥的毛细管(液桥的黏度大于非润湿流体的黏度), 其渗吸速度随着黏滞阻力的增加而下降的更明显. 根据Cox-Voinov公式, 当渗吸速度减小时, 动态接触角会随时间减小, 其余弦值cosθ*表现为逐渐增大的趋势. 当毛细管内含有液桥时, 其速度的减小幅度更大, 导致其接触角的余弦值cosθ*增大的幅度也更大.
![]()
图 9 不同黏度比条件下的动态接触角的时间演化及理论渗吸长度与LBM对比Figure 9. Time evolution of dynamic contact angle and comparison of theoretical imbibition length with LBM at different viscosity ratios此外, 在不同黏度比M的条件下, 忽略动态接触角的理论结果均与LBM存在较大差距. 并且, 液桥的存在导致二者的差距明显增大, 这再次说明了液桥的界面动态变化对渗吸速度有着显著的抑制作用. 需要指出的是, 对于M ≠ 1的两个案例, 其存在液桥时的理论渗吸长度与不存在液桥时的不一致, 需要单独计算; 而当M = 1时, 二者基本一致(图9(a) ~ 图9(c)中的黑色实线和黑色虚线). 通过将模拟得到的实时动态接触角数据(cos θ*)直接应用至方程(5)中, 对于不同黏度比的案例, 修正由界面的动态变化特性的方程(5)可以有效减小其与LBM结果之间偏差, 得到与模拟基本一致的渗吸长度.
3.5 液桥的界面动态特性的理论分析
在前两节的模拟中, 液桥的界面动态特性的评价需要依靠在模拟中捕捉界面的实时接触角. 然而, 对于含有液桥的体系, 捕捉每个界面的实时接触角是一个十分繁琐的工作. 为此, 笔者在本节中尝试用经典的水动力学模型[40-41] (Cox-Voinov公式)对已有的6组不含液桥的渗吸过程以及考虑液桥的渗吸过程的模拟结果进行补充分析. 首先, Cox-Voinov公式可以表示为
$$\tag{18a} \theta {}_d^3(t) - 9Ca\varGamma = \theta _0^3 $$ $$ \tag{18b}\theta {}_d^3(t) + 9Ca\varGamma = \theta _0^3 $$ 式中, Ca = Uμw/σ表示毛管数, θd(t)表示界面的实时动态接触角; Γ是一个拟合参数, 一般用以拟合实验数据或者模拟数据 (分子动力学模拟或者LBM) 确定其数值大小[42-44], 其具体的形式如下
$$ \varGamma = \ln H/\varepsilon $$ (19) 式中, ε表示微观的截止长度. 需要指出的是, 方程(18 a)适用于表征不含液桥的渗吸过程中的实时接触角, 以及考虑液桥渗吸过程中的左侧和右侧界面的实时接触角的表征. 方程(18b)适用于表征考虑液桥存在时的中间界面的实时接触角. 由于本文采用了4阶显式的Runge-Kutta方法求解渗吸的控制方程, 可以很方便地在每个时间步耦合Cox-Voinov公式[44]. 为了方便比较, 笔者在拟合模拟结果时, 使用了格子单位下的微观截止长度ε. 如图10所示, 本文共采用水动力学公式拟合了6组不含液桥时的渗吸长度的模拟结果(标签顺序依照3.3和3.4节的模拟顺序). 可以看出, 在不用捕捉实时动态接触角的前提下, 仅通过调节Cox-Voinov公式中的ε的数值, 其可以得到与LBM几乎一致的渗吸长度. 但需要注意的是, 对于不同的模拟参数(例如不同的静态接触角θ0以及黏度比M), ε并不是一个固定的数值. Ahadian 等[45]的分子模拟结果以及Wu等[46]的实验结果也证实了ε并非一个固定的数值. 另外, 利用模拟结果拟合得到的ε的数值均小于毛细管宽度, 说明本文的拟合结果是合理的.
![]()
图 10 耦合Cox-Voinov公式的理论模型与LBM预测的渗吸长度的比较(不存在液桥)Figure 10. Comparison of the theoretical model coupled with the Cox-Voinov equation with the imbibition length predicted by LBM (without liquid bridges)进一步对比了水动力学模型以及模拟预测的实时接触角的时间演化过程, 如图11所示. 可以看出, 对于不同的模拟参数, Cox-Voinov公式预测的实时动态接触角与LBM的总体趋势基本一致. 但二者仍然具有一定的差距(在流动充分发展之后, 二者余弦值最大的误差约为Δcosθd(t) ≈ 0.03). 我们认为造成这种现象的原因是: Cox-Voinov公式以及理论模型的前提是两相界面处的各个位置的移动速度是相等的, 因此采用Cox-Voinov公式预测的实时动态接触角是连续的光滑的. 但是, 在LBM模拟中, 即使在渗吸稳定时, 两相界面处的各个位置的实际移动速度并非完全一致, 尤其是接触线附近的流体粒子受到壁面较强的黏附力, 造成流体黏度的变化, 增强了黏性应力导致的界面变形[24]. 另一方面, 黏度的变化也会导致速度剖面与典型的Poiseuille速度剖面形成偏差, 使得界面各处的移动速度出现轻微的偏差, 进而导致实时接触角小幅度波动. 总的来说, 对于不含液桥的自发渗吸过程来说, Cox-Voinov公式具有一定的应用价值, 一方面其可以得到与LBM近似的一致的渗吸长度, 另一方面其也可以预测充分发展流动阶段后的实时动态接触角的演化趋势以及近似的数值. 但是, 需要采用大量模拟来确定相关的ε才能与模拟结果具有较好的匹配性.
![]()
图 11 不同模型参数条件下的Cox-Voinov模型实时动态接触角与LBM的比较(不存在液桥)Figure 11. Comparison of real-time dynamic contact angle of Cox-Voinov equation with LBM for different parameters (without liquid bridges)我们也测试了当毛细管内存在液桥时, 耦合Cox-Voinov公式的方程(5)是否能够预测此类渗吸过程中的动态接触角以及渗吸长度. 由于该过程中存在着3个润湿角, 按照渗吸的方向从左到右依次为, θ1, θ2和θ3. 因此, 在拟合模拟数据时, 既要保证每个接触角的数值相对应, 也需要保证整体的渗吸长度相对应. 笔者在观察模拟得到的3个实时的动态接触角的数据之后, 决定分别采用3个截断长度ε1, ε2和ε3拟合模拟结果. 其中ε1和ε3用于拟合左侧和右侧界面的接触角, ε2用于拟合中间界面的接触角. 具体而言, 对于左侧和右侧的两个界面, 其拟合形式与方程(18a)保持一致. 此外, 根据3.3和3.4节中的模拟数据, 发现θ1和θ3的数值相对较为接近, 为了减小拟合所需的计算量, 设置ε1 = ε3. 而对于中间界面, 其拟合形式与方程(18b)保持一致. 通过调整ε1, ε2和ε3, 并通过4阶显式的Runge-Kutta方法将3个Cox-Voinov公式耦合至方程(5)中计算实时的渗吸长度及3个动态接触角的余弦值. 在经过大量尝试后, 确定了合适的微观截止长度, 其拟合的计算结果如图12所示. 可以看出, 在没有3个界面的实时接触角的前提下, 耦合Cox-Voinov公式的方程(5)同样可以适用于计算存在液桥时的毛细管内的渗吸长度. 对于不同的模型参数, 二者得出的渗吸长度基本能够保持一致.
![]()
图 12 耦合Cox-Voinov公式的理论模型与LBM预测的渗吸长度的比较(存在液桥)Figure 12. Comparison of the theoretical model coupled with the Cox-Voinov equation with the imbibition length predicted by LBM (in the presence of liquid bridges)同时, 本文还对比了在含有液桥的自发渗吸过程中, Cox-Voinov公式预测的实时接触角的余弦值cosθ* (cosθ* = cosθ1 − cosθ2 + cosθ3)和LBM模拟捕捉的实时接触角的余弦值cosθ*, 如图13所示. 总体而言, 通过调节Cox-Voinov公式中的3个截断长度ε1, ε2和ε3可以得出与LBM较为一致的结果. 但是, 当液桥的右侧界面到达毛细管出口端时, 由于在模拟中右侧界面不能通过毛细管, 致使其动态接触角逐渐减小并趋近于90°. 其余两个界面也因此逐渐趋近于静态接触角, 系统的毛管压力趋近于平衡, 导致最终的cosθ*趋近于0. 显然, Cox-Voinov公式无法预测这一过程, 从而在渗吸将要结束时, 其与LBM模拟捕捉的实时接触角存在着一定的差距. 为了更加清晰地说明这一问题, 我们也给出了Cox-Voinov公式和LBM对每一个界面的实时接触角的捕捉, 如图14所示. 同样地, 在右侧界面到达出口端之前, 二者对于3个界面的实时接触角预测趋势及数值大小基本一致. 但是当渗吸过程中液桥的右侧界面到达出口端时, 由于Cox-Voinov公式无法预测右侧界面变形过程, 导致二者存在一定偏差. 但是, 这对于整体渗吸长度的影响很小, 因此Cox-Voinov公式对毛细管内存在液桥的渗吸过程的及液桥界面动态特性的捕捉仍然有效.
![]()
图 13 不同模型参数条件下的Cox-Voinov模型实时动态接触角与LBM的比较(存在液桥)Figure 13. Comparison of real-time dynamic contact angle of Cox-Voinov equation with LBM for different parameters (in the presence of liquid bridges)![]()
图 14 不同模型参数条件下的Cox-Voinov模型预测的各个界面的实时动态接触角与LBM的比较Figure 14. Comparison of real-time dynamic contact angles predicted by the Cox-Voinov equation with LBM under different parameters值得注意的是, 为了更准确地预测液桥的界面接近出口端时其右侧接触角的变化, 可以参考我们对3个界面实时捕捉的思路, 将3个微观截止长度ε1, ε2和ε3设置为时间的变量. 然而, 这样做会失去Cox-Voinov公式的简便性, 不利于广泛应用. 未来可以尝试进行大量的LBM模拟并采用Cox-Voinov公式拟合, 从而得出不同参数条件下的截断长度的图版, 这样可以增加Cox-Voinov公式应用的简便性.
4. 结 论
本文首先针对毛细管内存在液桥时的界面特性, 对经典的渗吸方程进行了拓展. 随后, 通过在Shan-Chen模型中引入中心矩碰撞算子并纳入改进虚拟密度的流−固作用力格式, 成功建立了液−液两相流动的格子Boltzmann模拟(LBM)方法. 在此基础上, 深入研究了液桥的存在及其界面动态特性(如动态接触角)对毛细管内液−液自发渗吸过程的影响, 并得出了以下结论.
(1) 液桥的存在会导致毛细管内的液−液自发渗吸过程中形成3个界面, 相较于不含液桥的单个界面渗吸过程, 液桥的存在显著增强了液−液渗吸过程中界面的动态变化特性, 大幅降低了毛管压力, 并导致毛细管内的渗吸速度降低, 对渗吸过程产生不利影响. 此外, 若忽略渗吸过程中的界面动态变化特性, 将导致理论预测的渗吸长度明显大于模拟结果. 通过将实时捕捉的每个界面的动态接触角纳入到理论模型中, 可以有效改善由界面动态变化特性导致的理论渗吸长度与模拟结果间的偏差.
(2) 液桥的润湿性对渗吸过程中的界面动态变化特性具有一定影响. 随着毛细管内液桥润湿性的增强, 界面在渗吸过程中的动态变化特性也逐渐增强. 相较于不含液桥的渗吸过程, 随着液桥润湿性的增强, 系统的渗吸速度降低的幅度更大, 导致渗吸长度与不考虑动态接触角的理论值之间存在较大差异.
(3) 当液桥(即润湿流体)的黏度小于非润湿流体的黏度时, 随着渗吸长度的增加, 液桥的界面动态变化特性逐渐增强. 而当液桥的黏度等于非润湿流体的黏度时, 在渗吸过程中系统的动态接触角在相当长的一段时间内保持稳定. 然而, 当液桥的黏度大于非润湿流体的黏度时, 系统的界面动态变化特性随着渗吸的进行逐渐减弱. 随着液桥黏度的减小, 液桥界面的动态变化特性对渗吸长度的影响愈发显著. 为了准确预测含有液桥的毛细管内的渗吸过程, 必须捕捉界面特性并对理论模型进行修正.
(4) 无论是对于存在液桥的渗吸过程还是不存在液桥的渗吸过程, 利用Cox-Voinov公式拟合模拟数据并确定合理的拟合参数后, 都可以在不捕捉实时动态接触角的前提下, 准确预测渗吸过程中的渗吸长度和动态接触角. 但不同的是, 对于不考虑液桥的渗吸过程, 只需要拟合一个拟合参数; 而当毛细管内存在液桥时, 必须同时对3个界面的动态变化特性进行拟合, 才能与模拟得出的渗吸长度和界面的实时动态接触角保持良好的匹配性.
-
![]()
图 1 毛细管内存在液桥时的液−液自发渗吸示意图
Figure 1. Schematic diagram of spontaneous liquid-liquid imbibition in a capillary with a liquid bridge
![]()
![]()
图 5 毛细管内存在不同长度液桥的液−液渗吸过程
Figure 5. Liquid-liquid imbibition process in capillaries with different liquid bridge lengths
![]()
图 6 渗吸过程中毛细管内存在不同长度的液桥时的动态接触角的时间演化
Figure 6. Time evolution of the dynamic contact angle during imbibition for the presence of liquid bridges with different lengths in a capillary
![]()
图 7 渗吸过程中存在不同液桥长度时的理论渗吸长度与LBM的比较
Figure 7. Comparison of theoretical imbibition length with LBM during imbibition with different liquid bridge lengths
![]()
图 8 不同润湿性条件下的动态接触角的时间演化及理论渗吸长度与LBM对比
Figure 8. Time evolution of dynamic contact angle and comparison of theoretical imbibition length with LBM under different wettability
![]()
图 9 不同黏度比条件下的动态接触角的时间演化及理论渗吸长度与LBM对比
Figure 9. Time evolution of dynamic contact angle and comparison of theoretical imbibition length with LBM at different viscosity ratios
![]()
图 10 耦合Cox-Voinov公式的理论模型与LBM预测的渗吸长度的比较(不存在液桥)
Figure 10. Comparison of the theoretical model coupled with the Cox-Voinov equation with the imbibition length predicted by LBM (without liquid bridges)
![]()
图 11 不同模型参数条件下的Cox-Voinov模型实时动态接触角与LBM的比较(不存在液桥)
Figure 11. Comparison of real-time dynamic contact angle of Cox-Voinov equation with LBM for different parameters (without liquid bridges)
![]()
图 12 耦合Cox-Voinov公式的理论模型与LBM预测的渗吸长度的比较(存在液桥)
Figure 12. Comparison of the theoretical model coupled with the Cox-Voinov equation with the imbibition length predicted by LBM (in the presence of liquid bridges)
![]()
图 13 不同模型参数条件下的Cox-Voinov模型实时动态接触角与LBM的比较(存在液桥)
Figure 13. Comparison of real-time dynamic contact angle of Cox-Voinov equation with LBM for different parameters (in the presence of liquid bridges)
![]()
图 14 不同模型参数条件下的Cox-Voinov模型预测的各个界面的实时动态接触角与LBM的比较
Figure 14. Comparison of real-time dynamic contact angles predicted by the Cox-Voinov equation with LBM under different parameters
表 1 模型参数与实际参数的单位转换
Table 1 Unit conversion of model and real parameters
Parameters Lattice values Physical values Scale factor density, ρw 1.0 1000 kg/m3 1000 kinematic viscosity, νw 0.1667 5 × 10−6 m2/s 3 × 10−5 lattice length, Δx 1.0 2.63 × 10−6 m 2.63 × 10−6 interface tension, σ 0.088 3 × 10−2 N/m 0.341 
下载: 导出CSV
-
[1] 蔡建超, 郁伯铭. 多孔介质自发渗吸研究进展. 力学进展, 2012, 42(6): 735-754 (Cai Jianchao, Yu Boming. Advances in studies of spontaneous imbibition in porous media. Advances in Mechanics, 2012, 42(6): 735-754 (in Chinese) doi: 10.6052/1000-0992-11-096 Cai Jianchao, Yu Boming. Advances in studies of spontaneous imbibition in porous media. Advances in Mechanics, 2012, 42(6): 735-754 (in Chinese) doi: 10.6052/1000-0992-11-096
[2] 蔡建超. 多孔介质自发渗吸关键问题与思考. 计算物理, 2021, 38(5): 505-512 (Cai Jianchao. Some key issues and thoughts on spontaneous imbibition in porous media. Chinese Journal of Computational Physics, 2021, 38(5): 505-512 (in Chinese) doi: 10.19596/j.cnki.1001-246x.8440 Cai Jianchao. Some Key Issues and Thoughts on Spontaneous Imbibition in Porous Media. Chinese Journal of Computational Physics, 2021, 38(5): 505-512 (in Chinese) doi: 10.19596/j.cnki.1001-246x.8440
[3] Cai JC, Perfect E, Cheng CL, et al. Generalized modeling of spontaneous imbibition based on Hagen–Poiseuille flow in tortuous capillaries with variably shaped apertures. Langmuir, 2014, 30(18): 5142-5151 doi: 10.1021/la5007204
[4] Tian WB, Wu KL, Gao Y, et al. A critical review of enhanced oil recovery by imbibition: theory and practice. Energy & Fuels, 2021, 35(7): 5643-5670
[5] 朱维耀, 岳明, 刘昀枫等. 中国致密油藏开发理论研究进展. 工程科学学报, 2019, 41(9): 1103-1114 ((Zhu Weiyao, Yue Ming, Liu Yunfeng, et al. Research progress on tight oil exploration in China. Chinese Journal of Engineering, 2019, 41(9): 1103-1114 (in Chinese) Zhu Weiyao, Yue Ming, Liu Yunfeng, et al. Research progress on tight oil exploration in China. Chinese Journal of Engineering, 2019, 41(9): 1103-1114 (in Chinese
[6] 李爱芬, 何冰清, 雷启鸿等. 界面张力对低渗亲水储层自发渗吸的影响. 中国石油大学学报 (自然科学版), 2018, 42(4): 67-74 (Li Aifen, He Bingqing, Lei Qihong, et al. Influence of interfacial tension on spontaneous imbibition in low-permeability water-wet reservoirs. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(4): 67-74 (in Chinese) Li Aifen, He Bingqing, Lei Qihong, et al. Influence of interfacial tension on spontaneous imbibition in low-permeability water-wet reservoirs. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(4): 67-74 (in Chinese)
[7] Mondal PK, DasGupta D, Chakraborty S. Interfacial dynamics of two immiscible fluids in spatially periodic porous media: The role of substrate wettability. Physical Review E, 2014, 90(1): 013003 doi: 10.1103/PhysRevE.90.013003
[8] 李锡夔, 张松鸽, 楚锡华. 非饱和颗粒材料的多孔连续体有效压力与有效广义 Biot 应力. 力学学报, 2022, 55(2): 369-380 (Li Xikui, Zhang Songge, Chu Xihua. Effective pressure and generalized effective Biot stress of porous continuum in unsaturated granular materials. Chinese Journal of Theoretical and Applied Mechanics, 2022, 55(2): 369-380 (in Chinese) doi: 10.6052/0459-1879-21-483 Li Xikui, Zhang Songge, Chu Xihua. Effective pressure and generalized effective Biot stress of porous continuum in unsaturated granular materials. Chinese Journal of Theoretical and Applied Mechanics, 2022, 55(2): 369-380 (in Chinese) doi: 10.6052/0459-1879-21-483
[9] Liu Y, Berg S, Ju Y, et al. Systematic investigation of corner flow impact in forced imbibition. Water Resources Research, 2022, 58(10): e2022WR032402 doi: 10.1029/2022WR032402
[10] 张晟庭, 李靖, 陈掌星等. 气液非混相驱替过程中的卡断机理及模拟研究. 力学学报, 2022, 54(5): 1429-1442 (Zhang Shengting, Li Jing, Chen Zhangxing, et al. Study on snap-off mechanism and simulation during gas-liquid immiscible displacement. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1429-1442 (in Chinese) Zhang Shengting, Li Jing, Chen Zhangxing, et al. Study on snap-off mechanism and simulation during gas-liquid immiscible displacement. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1429-1442 (in Chinese)
[11] 杨敏, 曹炳阳. 微纳通道中牛顿流体毛细流动的研究进展. 科学通报, 2016, 61(14): 1574-1584 (Yang Min, Cao Bingyang. Advances of capillary filling of Newtonian fluids in micro- and nanochannels. Chinese Science Bulletin, 2016, 61(14): 1574-1584 (in Chinese) doi: 10.1360/N972015-00783 Yang Min, Cao Bingyang. Advances of capillary filling of Newtonian fluids in micro- and nanochannels. Chinese Science Bulletin, 2016, 61(14): 1574-1584 (in Chinese) doi: 10.1360/N972015-00783
[12] Fortes MA. Axisymmetric liquid bridges between parallel plates. Journal of Colloid and Interface Science, 1982, 88(2): 338-352 doi: 10.1016/0021-9797(82)90263-6
[13] Lian GP, Thornton C, Adams MJ. A theoretical study of the liquid bridge forces between two rigid spherical bodies. Journal of Colloid and Interface Science, 1993, 161(1): 138-147 doi: 10.1006/jcis.1993.1452
[14] Salama A, Cai JC, Kou JS, et al. Investigation of the dynamics of immiscible displacement of a ganglion in capillaries. Capillarity, 2021, 4(2): 31-44 doi: 10.46690/capi.2021.02.02
[15] Lian GP, Seville J. The capillary bridge between two spheres: new closed-form equations in a two century old problem. Advances in Colloid and Interface Science, 2016, 227: 53-62 doi: 10.1016/j.cis.2015.11.003
[16] Lucas R. Rate of capillary ascension of liquids. Kolloid Z, 1918, 23(15): 15-22
[17] Washburn EW. The dynamics of capillary flow. Phys. Rev., 1921, 17(3): 273-283 doi: 10.1103/PhysRev.17.273
[18] Mumley TE, Radke CJ, Williams MC. Kinetics of liquid/liquid capillary rise: I. Experimental observations. Journal of Colloid and Interface Science, 1986, 109(2): 398-412 doi: 10.1016/0021-9797(86)90318-8
[19] Fermigier M, Jenffer P. An experimental investigation of the dynamic contact angle in liquid-liquid systems. Journal of Colloid and Interface Science, 1991, 146(1): 226-241 doi: 10.1016/0021-9797(91)90020-9
[20] Walls PLL, Dequidt G, Bird JC. Capillary displacement of viscous liquids. Langmuir, 2016, 32(13): 3186-3190 doi: 10.1021/acs.langmuir.6b00351
[21] André J, Okumura K. Capillary Replacement in a tube prefilled with a viscous fluid. Langmuir, 2020, 36(37): 10952-10959 doi: 10.1021/acs.langmuir.0c01612
[22] Krüger T, Kusumaatmaja H, Kuzmin A, et al. The lattice Boltzmann method. Springer International Publishing, 2017, 10(978-3): 4-15
[23] Moradi B, Ghasemi S, Hosseini MA, et al. Dynamic behavior investigation of capillary rising at various dominant forces using free energy lattice Boltzmann method. Meccanica, 2021, 56(12): 2961-2977 doi: 10.1007/s11012-021-01426-z
[24] Latva-Kokko M, Rothman DH. Scaling of dynamic contact angles in a lattice-Boltzmann model. Physical Review Letters, 2007, 98(25): 254503 doi: 10.1103/PhysRevLett.98.254503
[25] 张晟庭, 李靖, 陈掌星等. 基于改进 LBM 的气液自发渗吸过程中动态润湿效应模拟. 力学学报, 2023, 55(2): 355-368 (Zhang Shengting, Li Jing, Chen Zhangxing, et al. Simulation of dynamic wetting effect during gas-liquid spontaneous imbibition based on modified LBM. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(2): 355-368 (in Chinese) Zhang Shengting, Li Jing, Chen Zhangxing, et al. Simulation of dynamic wetting effect during gas-liquid spontaneous imbibition based on modified LBM. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(2): 355-368 (in Chinese)
[26] Chen L, Kang QJ, Mu YT, et al. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer, 2014, 76: 210-236 doi: 10.1016/j.ijheatmasstransfer.2014.04.032
[27] Chibbaro S, Biferale L, Diotallevi F, et al. Capillary filling for multicomponent fluid using the pseudo-potential lattice Boltzmann method. The European Physical Journal Special Topics, 2009, 171(1): 223-228 doi: 10.1140/epjst/e2009-01032-8
[28] Budaraju A, Phirani J, Kondaraju S, et al. Capillary displacement of viscous liquids in geometries with axial variations. Langmuir, 2016, 32(41): 10513-10521 doi: 10.1021/acs.langmuir.6b02788
[29] Qian JY, Li XJ, Wu Z, et al. A comprehensive review on liquid–liquid two-phase flow in microchannel: Flow pattern and mass transfer. Microfluidics and Nanofluidics, 2019, 23: 1-30 doi: 10.1007/s10404-018-2168-8
[30] Chan WK, Yang C. Surface-tension-driven liquid–liquid displacement in a capillary. Journal of Micromechanics and Microengineering, 2005, 15(9): 1722 doi: 10.1088/0960-1317/15/9/014
[31] Hajabdollahi F, Premnath KN, Welch S WJ. Central moment lattice Boltzmann method using a pressure-based formulation for multiphase flows at high density ratios and including effects of surface tension and Marangoni stresses. Journal of Computational Physics, 2021, 425: 109893 doi: 10.1016/j.jcp.2020.109893
[32] Fei LL, Qin FF, Zhao JL, et al. Lattice Boltzmann modelling of isothermal two-component evaporation in porous media. Journal of Fluid Mechanics, 2023, 955: A18 doi: 10.1017/jfm.2022.1048
[33] Porter ML, Coon ET, Kang Q, et al. Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios. Physical Review E, 2012, 86(3): 036701 doi: 10.1103/PhysRevE.86.036701
[34] Akai T, Bijeljic B, Blunt MJ. Wetting boundary condition for the color-gradient lattice Boltzmann method: Validation with analytical and experimental data. Advances in Water Resources, 2018, 116: 56-66 doi: 10.1016/j.advwatres.2018.03.014
[35] Li Q, Yu Y, Luo KH. Implementation of contact angles in pseudopotential lattice Boltzmann simulations with curved boundaries. Physical Review E, 2019, 100(5): 053313 doi: 10.1103/PhysRevE.100.053313
[36] Coelho RCV, Moura CB, Teloda Gama MM, et al. Wetting boundary conditions for multicomponent pseudopotential lattice Boltzmann. International Journal for Numerical Methods in Fluids, 2021, 93(8): 2570-2580 doi: 10.1002/fld.4988
[37] Sedahmed M, Coelho RCV, Warda HA. An improved multicomponent pseudopotential lattice Boltzmann method for immiscible fluid displacement in porous media. Physics of Fluids, 2022, 34(2): 023102 doi: 10.1063/5.0080823
[38] Sedahmed M, Coelho RCV, Araújo NAM, et al. Study of fluid displacement in three-dimensional porous media with an improved multicomponent pseudopotential lattice Boltzmann method. Physics of Fluids, 2022, 34(10): 103303 doi: 10.1063/5.0107361
[39] Thamdrup LH, Persson F, Bruus H, et al. Experimental investigation of bubble formation during capillary filling of SiO2 nanoslits. Applied Physics Letters, 2007, 91(16): 163505 doi: 10.1063/1.2801397
[40] Cox RG. The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. Journal of Fluid Mechanics, 1986, 168: 169-194 doi: 10.1017/S0022112086000332
[41] Voinov OV. Hydrodynamics of wetting. Fluid Dynamics, 1976, 11(5): 714-721
[42] Primkulov BK, Chui JYY, Pahlavan AA, et al. Characterizing dissipation in fluid-fluid displacement using constant-rate spontaneous imbibition. Physical Review Letters, 2020, 125(17): 174503 doi: 10.1103/PhysRevLett.125.174503
[43] Zhou ZX, Yang ZB, Luo C, et al. Influence of inertia on liquid splitting at fracture intersections. Journal of Hydrology, 2023, 618: 129270 doi: 10.1016/j.jhydrol.2023.129270
[44] Hubao A, Yang ZB, Hu R, et al. Roles of energy dissipation and asymmetric wettability in spontaneous imbibition dynamics in a nanochannel. Journal of Colloid and Interface Science, 2022, 607: 1023-1035 doi: 10.1016/j.jcis.2021.09.051
[45] Ahadian S, Mizuseki H, Kawazoe Y. On the kinetics of the capillary imbibition of a simple fluid through a designed nanochannel using the molecular dynamics simulation approach. Journal of Colloid and Interface Science, 2010, 352(2): 566-572 doi: 10.1016/j.jcis.2010.09.011
[46] Wu PK, Nikolov AD, Wasan DT. Capillary rise: Validity of the dynamic contact angle models. Langmuir, 2017, 33(32): 7862-7872 doi: 10.1021/acs.langmuir.7b01762
计量
- 文章访问数: 737
- HTML全文浏览量: 219
- PDF下载量: 94
