多孔介质中含流固溶解反应的互溶驱替过程在地下水污染治理^[1]、二氧化碳埋存^[2]及温室气体强化采油^[3]等领域广泛存在, 是能源、环境等学科共同关注的典型问题. 如在地下水污染治理中, 污水会不断侵蚀和驱替地下水, 同时污水中的酸根离子会与地层中的岩石(如方解石、白云石等)发生溶解反应, 使岩石结构发生显著变化, 进而影响流体流动. 在此驱替过程中, 流动、扩散以及化学反应三者同时存在, 相互耦合. 如何有效地控制含溶解反应下的驱替是相关领域的关键技术之一.

关于互溶驱替问题的研究由来已久. 当驱替流体黏性小于被驱替流体时, 可能产生黏性指进现象. Saffman和Taylor^[4]采用Hele-Shaw装置对多孔介质内的驱替问题进行实验研究, 并提供了黏性指进指端发展的数学模型; Zimmerman和Homsy^[5]通过数值模拟研究了二维非均质多孔介质中互溶驱替过程; Wit和Homsy^[6]研究了多孔介质中流体间化学反应与黏性指进的非线性相互作用, 发现反应会使指进形态特征和传播机制发生显著的变化; 之后, Nagatsu等^[7]通过实验及理论方法研究了流体间沉淀反应对指进的影响, 并讨论了渗透率变化与指进发展之间的关系. 需要指出的是, 上述研究涉及的化学反应均是流体与流体之间的化学反应, 少有对流体与多孔介质骨架之间存在化学反应的互溶驱替现象的研究.

另一方面, 研究人员采用不同方法研究了多孔介质中含化学反应的各种问题. Békri等^[8]率先将随机模型和矩量法相结合以计算固体壁面处的反应速率; Luo等^[9]和Oltéan等^[10]通过arbitrary-Lagrangian-Eulerian (ALE)方法求解欧拉网格上的物理问题(如溶质扩散、界面反应等), 该方法需要及时更新欧拉网格; Soulaine等^[11]提出了一种基于Dancy-Brinkman-Stokes (DBS)方程的方法来模拟单相流动下固体的溶解反应, 并通过实验验证了该方法的有效性. 由于孔隙尺度下多孔介质内固体骨架形状复杂, 尤其当固体骨架随着反应发生而不断变化时, 传统数值方法处理起来难度较大. 近三十年来, 格子Boltzamnn方法(lattice Boltzmann method, LBM)作为一种介观的模拟方法, 可以更方便地处理复杂微观孔隙结构内流体与固体之间、不同流体组分之间复杂的相互作用^[12-14], 对于多孔介质内流动与反应耦合的问题, 也受到了学者们的广泛关注. 如He等^[15]将流体-壁面间的表面反应与扩散机理相耦合, 通过边界条件实现了表面反应的产生; Kang等^[16]在He模型的基础上建立了一种格子Boltzmann模型, 模拟了多孔介质中两种水合物的流动与反应过程, 并考虑了多孔介质几何形状的演化, 该模拟结果与实验观察的现象相符合; 之后又进一步研究了多组分溶液在多孔介质中的反应性传输^[17]; 张婷等^[18]详细分析了各无量纲参数对多孔介质内溶解和沉淀过程的影响; Ju等^[19]提出了一种适用于不规则表面的反应边界格式, 并成功将其拓展至三维模型; Meng等^[20]从多尺度分析了非均相反应对密度驱动流的影响, 发现界面的不稳定性由流体密度差、反应速率和孔隙率变化三者共同决定.

综上所述, 虽然目前的研究对于多孔介质互溶驱替过程和多孔介质中含化学反应的问题均有充足的认识, 但对于多孔介质内含溶解反应的互溶驱替过程研究尚不充分, 此过程中的驱替规律尚待研究. 本文采用LBM对孔隙尺度下多孔介质内含溶解反应的互溶驱替过程进行研究, 着重分析驱替过程中溶解反应对多孔介质结构和驱替效率造成的影响, 寻求控制流体驱替的方法, 为地下水污染治理, 二氧化碳混相驱油等技术提供思路.

本文研究孔隙尺度下多孔介质中含溶解反应的互溶流体驱替过程. 其中, 多孔介质固体骨架结构是通过对CT扫描得到的多孔介质灰度图进行二值化处理而获得的, 多孔介质区域宽度和高度分别为$ l_x $和$ l_y $, $ l_x = l_y $. 多孔介质孔隙率$ \varepsilon $的定义为

其中$ V_{{\rm{pore}}} $为孔隙所占体积, $ V_{{\rm{total}}} $为多孔介质的总体积. 图1所示的多孔介质中, 蓝黑色不规则部分代表的是多孔介质固体骨架B, 该多孔介质初始孔隙率$ \varepsilon_0 $为0.864, 下标“0”表示初始时刻.

系统设置如图1所示. 在初始时刻, $0.05 l_x \leqslant $$ x \leqslant l_x$区域的多孔介质中充满了黏性为$ \mu_2 $的被驱替流体(Fluid 2), 该流体处于静止状态. 剩余部分的多孔介质中充满了黏性为$ \mu_1 $的驱替流体(Fluid 1), 然后从多孔介质左侧以恒定速度$ {\boldsymbol{u}}_0 $持续注入该驱替流体, 两流体等温, 可互溶, 且均为不可压缩流体. 在驱替流体中含有溶质A, 流体无量纲浓度为$ C_1 = 1 $; 被驱替流体中不含该溶质, 流体无量纲浓度$ C_2 = 0 $. 在驱替过程中溶质A与固体骨架B发生溶解反应, A + B→P, 反应产物P为溶于流体的酸根离子.

在孔隙尺度下, 多孔介质内流体的流动可用如下方程组来刻画:

其中(2)式为不可压流体的连续性方程, (3)式为不可压Navier-Stokes方程. u为流体流过多孔介质孔隙时的速度, t为时间, p为压力, ν为流体的运动黏度.

多孔介质内两种流体间溶质的扩散采用对流扩散方程来描述. 本文考虑的是流体与固体之间的溶解反应, 不包括流体之间的化学反应. 因此, 对流扩散方程中不包含源项, 即

其中, C为流体中溶质的浓度, D为溶质间的扩散系数.

系统边界条件如下: 当$ x = 0 $时, ${\boldsymbol{u}} = (u_0, 0), $$ C = C_1$; 当$ x = l_x $时, $\dfrac{{\partial u}}{{\partial t}} + {u_0}\dfrac{{\partial u}}{{\partial x}} = 0$, ${\rm{ }}\dfrac{{\partial C}}{{\partial x}} = 0$; 上下边界的速度场满足无滑移条件, 对于上下边界的浓度场, 则有$ \nabla C = 0 $.

此外, 对于多孔介质固体骨架与溶质间的溶解反应, 也通过边界条件实现. 考虑到实际工程中溶质与固体间的溶解反应机理过于复杂, 本文对该反应进行简化处理, 采用如下的线性反应边界条件^[18]来描述:

其中, $ k_{\rm{r}} $为溶解反应的反应速率常数, n为固体壁面指向流体的法线方向, $ C_{\rm{w}} $为固体壁面处流体的浓度, $ C_{{\rm{eq}}} $代表平衡浓度, 本文中溶解反应是不可逆的, $ C_{{\rm{eq}}} $设为0.

该系统有4组无量纲参数, 分别为雷诺数(Reynolds number, $ Re $), $ Re = {u_0}{l_x}/\nu_2 $; 佩克莱数(Pèlcet number, Pe), $ Pe = {u_0}{l_x}/D $; 黏性比(viscosity ratio, M)定义为两流体的黏性之比, $ M = $$ \mu_2/\mu_1 $; 达姆科勒数(Damkohler number, Da), $ Da = k_{\rm{r}}l_y/D $. 在驱替过程中, 由于两种流体互溶, 流体黏性与浓度相关, 即

其中$ R = \ln M $. 需要注意的是, 当驱替流体黏性远小于被驱替流体时, 会产生明显的指进现象^[21,22]. 本文研究的驱替问题考虑的均是两互溶流体黏性差异较大的情况, 黏性比设定为$ M = 100 $.

本文采用多孔介质中流体流动及扩散的耦合格子Boltzmann (lattice Boltzmann, LB)模型^[23]来求解(2)式—(4)式组成的方程组. 该模型为双分布模型, 用两组演化方程分别描述速度场和浓度场.

速度场可用如下演化方程来描述:

其中q是离散速度的个数, $ {f_i}({{\boldsymbol{x}}}, t) $是粒子在t时刻x处离散速度为$ {\boldsymbol{c}}_i $ 的分布函数, $ \delta_t $为时间步长, $\varOmega _i$是碰撞矩阵. 流体的密度和速度可以表示为

已有研究证明多松弛(multiple-relaxtion-time, MRT)模型在模拟孔隙尺度下多孔介质内流动问题时, 相比于Bhatnagar-Gross-Krook (BGK)模型, 具有更好的数值稳定性^[24]. 故本文采用MRT模型求解流动方程.

MRT模型的碰撞算子可表示为

其中m和$ {\boldsymbol{m}}^{({\rm{eq}})} $分别是矩向量和矩空间中的平衡态, M是将分布函数映射到矩空间的$ q\times q $的变换矩阵, $ {\boldsymbol{m}} = {\boldsymbol{M}}\cdot{\boldsymbol{f}} $, ${\boldsymbol{S}} = {{\rm{diag}}}(s_0, s_1, \cdots, s_{q-1})$为松弛矩阵.

本文考虑的是二维问题, 采用二维九速(D2Q9)离散速度模型, 该模型中离散速度$ {\boldsymbol{c}}_i $定义为

其中$ c = \delta_x/\delta_t $, $ \delta_x $和$ \delta_t $分别为网格步长和时间步长, 本文取$ c = 1 $.

变换矩阵M为

其对应的矩向量为

其中T表示转置变换. (12)式中的各个矩的物理含义如下: ρ是流体密度; e和ε分别是总能和能量的平方; $ j_x $和$ j_y $分别为x和y分量上的动量, $j_x = $$ \rho u_x$, $ j_y = \rho u_y $; $ q_x $和$ q_y $则是x和y分量上的能量通量; $ p_{xx} $和$ p_{xy} $分别为应力张量的对角线和非对角线上的元素. 对于不可压流体, 流体密度近似为常数, 记作$ \rho_0 $, 而密度波动为$ \text{δ}\rho $, 从而有$ \rho = \rho_0+\text{δ}\rho $. 相应矩空间中的平衡态如下:

其中, 流体密度ρ和动量$ {\boldsymbol{j}} = \rho_0 {\boldsymbol{u}} = (j_x, j_y) $是碰撞守恒量. 而其他非守恒的平衡态是守恒量的函数, 具有如下形式:

此模型的松弛参数为

注意到流体密度和动量在碰撞过程中守恒, 这些矩对应的松弛率$ s_\rho $和$ s_j $可取任意值. 其他松弛率为$ s_e = s_\varepsilon = s_\nu = s_\nu = 1/\tau $, $ s_q = 8(2-s_\nu)/(8-s_\nu) $^[24], 其中τ与流体动力黏性ν相关,

浓度场则用另一组演化方程来描述:

其中$ g_i({\boldsymbol{x}}, t) $为浓度C的分布函数, 离散速度和上文中用于描述速度场的MRT模型中的离散速度一致; 平衡态分布函数$ g_i^{({\rm{eq}})} $为

其中$ \omega_i $为权系数, $ \omega_0 = 4/9 $, $ \omega_{1-4} = 1/9 $, $\omega_{5-8} = $$ 1/36$, 而$ c^2_s = c^2/3 $. 参数A与扩散系数D相关:

对于扩散系数较小(或Pe数较大)的情况, 可通过调整参数A来提高数值稳定性.

此模型中, 浓度定义如下:

(23)式中的浓度梯度可由该点非平衡态函数的一阶矩得到, 计算如下:

即

对于不同的边界条件, 均需进行边界处理. 其中左右边界均采用非平衡外推格式进行处理^[25]. 对于上下边界和多孔介质内部的固体骨架, 速度满足无滑移条件, 采用半反弹格式来处理^[26]. 标准反弹格式固体边界位于格点上, 而半反弹格式的固体边界位于格点的中间, 即$({\boldsymbol{x}}_{\rm{f}}+{\boldsymbol{x}}_{\rm{b}})/2$处, 其中$ {\boldsymbol{x}}_{\rm{b}} $是固体格点, $ {\boldsymbol{x}}_{\rm{f}} $是临近固体壁面的流体格点.

其中$ \text{δ}t $为时间步长, $ {\boldsymbol{c}}_{\overline \imath} = -{\boldsymbol{c}}_i $($ {\boldsymbol{c}}_i $指向壁面), $ {f'_i}\left( {{{\boldsymbol{x}}_{\rm{f}}}, t} \right) $为$ {f_i}\left( {{{\boldsymbol{x}}_{\rm{f}}}, t} \right) $碰撞后的值.

对于上下边界的浓度场, $ \nabla C = 0 $. 可用一种类反弹格式来实现^[27]:

其中$ {g'_i} $是$ g_i $碰撞后的分布函数, 定义见(22)式等号右边.

多孔介质固体骨架与溶质间的溶解反应, 用(5)式来表示, 采用修正半反弹格式^[28]来处理:

其中$ C_{\rm{w}} $可由有限差分格式得出:

其中, $ C_{\rm{f}} $为固体壁面附近流体节点的浓度.

溶解反应在固相上的表现为固体体积的不断变化, 同时假定固体骨架主要由易溶盐矿物组成, 能够全部溶解. 节点处固体体积随时间的变化规律满足公式^[17]:

其中, V表示多孔介质固体骨架的体积, $ V_{\rm{m}} $表示其摩尔体积, $ a_{\rm{s}} $表示化学计量数, 此处即反应表面积. 每个固体节点表示为一个控制体积, 初始时其大小为$ 1\times 1 $的单位网格单元. 在发生反应时, 固体节点的体积在每个时刻依据以下式子更新:

当固体节点的体积$ V = 0 $时, 说明该部分固体已被完全溶解, 节点被更新为流体节点, 其浓度由附近流体点浓度的平均值求得, 而速度则取为0.

本节通过模拟两平板间的互溶驱替过程, 对所使用的耦合LB模型进行验证. 如图2所示, 计算区域为$0 \leqslant x \leqslant l_x$, $0 \leqslant y \leqslant l_y$; 初始时刻, 计算区域内$0.1 l_x \leqslant x \leqslant l_x$, $0 \leqslant y \leqslant l_y$的部分充满了浓度为$ C = C_2 = 0 $, 动力黏度为$ \mu = \mu_2 $的被驱替流体, 剩余部分充满了浓度为$ C = C_1 =1$, 动力黏度为$ \mu = $$ \mu_1 $的驱替流体; 之后, 从计算区域左边界以速度$ {\boldsymbol{u}} = (0.005, 0) $持续注入该驱替流体. 本节采用的网格大小为$ 256\times32 $.

在此驱替过程中, 系统的控制方程为(2)式—(4)式组成的方程组, 本节的边界条件给定如下: 当$ x = 0 $时, ${\boldsymbol{u}} = (u_0, 0),\;\; c = c_1$; 当$ x = l_x $时, $\dfrac{{\partial {\boldsymbol{u}}}}{{\partial t}} + $$ {u_0}\dfrac{{\partial {\boldsymbol{u}}}}{{\partial x}} = 0$, ${\rm{ }}\dfrac{{\partial c}}{{\partial x}} = 0$; 上下边界为壁面, 速度满足无滑移条件, 对于浓度场, $ \nabla C = 0 $. 本节模拟了文献[21]中的4种工况: 1) $ M = 1 $, $ Pe = 5 $; 2) $M = 1$, $ Pe = 262 $; 3) $ M = 100 $, $ Pe = 5 $; 4) $ M = 100 $, $ Pe = 262 $. 模拟结果如图3所示, 图中还给出了纵截面的平均浓度分布图, 通过对比可以发现与文献中结果完全一致, 即验证了该耦合LB模型模拟互溶驱替过程的准确性.

为确保选取网格的计算精度和计算效率, 在数值模拟之前要进行网格无关性的验证. 计算区域的长和宽为$ l_x = l_y = 32 $, 分别选择$ 512\times512 $和$ 1024\times1024 $的网格尺寸来描述图1中的多孔介质系统. 在本次模拟中, 设定$ {\boldsymbol{u}}_0 = (0.005, 0) $, $ Re = 32 $, $ Pe = 262 $及$ Da = 0 $. 为量化不同网格间的差异, 选取$ y = 16 $处的横截面作为观测面, 图4为观测面上流体的平均浓度和平均速度随时间的演变图. 需要注意的是, 在本文的模拟中, 特征时间$ \tau = l_y/u_0 $, 无量纲时间$ t = \text{ δ} tN_{{\rm{step}}}/\tau $, 其中$ N_{{\rm{step}}} $为模拟所用时间步, 后文中t均为无量纲时间. 如图4所示, 在$ 512{\times}512 $和$ 1024{\times}1024 $ 这两种网格下演化过程符合良好, 因此选取$ 512\times 512 $的网格即可满足网格无关.

在含溶解反应的驱替过程中, Da数为表征溶解反应速率的无量纲参数, 本节将研究Da数对驱替过程的影响. 首先模拟了3种不同工况下的互溶驱替过程, 工况一为无溶解反应($ Da = 0 $)的互溶驱替过程; 工况二为含有溶解反应且反应速率较小($ Da = 4 $)的互溶驱替; 工况三中溶解反应速率较大($ Da = 32 $). 雷诺数设定为$ Re = 32 $, 佩克莱数设定为$ Pe = 262 $, 黏性比设定为$ M = 100 $. 图5给出了3种工况下不同时刻流体的浓度场图. 需要指出的是, 浓度场图中蓝黑色不规则部分代表的是多孔介质固体骨架B, 随着溶解反应的发生, 部分固体骨架被溶解, 多孔介质结构随之发生改变.

如图5(a)所示, 在驱替过程早期($ t = 0.0488 $), 驱替流体刚进入多孔介质, 与多孔介质固体骨架的接触面积较小, 溶解反应的量也较小, 因此3种工况下流体的浓度场之间无明显差异.

随着驱替过程的进行($ t = 0.2441 $, 见图5(b)), 受驱替流体和被驱替液流体间的黏性差异以及多孔介质孔隙结构两种因素的影响, 在两流体界面上产生黏性指进现象. 此外, 在工况二及工况三中, 驱替流体与多孔介质固体骨架发生流固溶解反应, 部分固体骨架逐渐被溶解, 扩大了多孔介质的孔隙率. 但由于溶解反应的发生消耗了一部分溶质, 因此Da数较大时指进发展的速度相对较慢. 总体上, 此时3种工况下流体的浓度场差异仍不大.

随后($ t = 0.4883 $, 见图5(c)), 随着驱替流体不断向前推进, 工况二和工况三中溶解反应沿指进发展路径不断发生, 从而在多孔介质内部形成裂缝, 使后续驱替流体沿裂缝前进. 与此同时, 在工况三中, 驱替流体与多孔介质固体骨架间持续发生溶解反应, 促进了裂缝的扩展, 形成两条明显的通道, 即虫洞现象^[29-31], 而工况二中由于溶解反应速率较小, 没有出现明显的虫洞. 由此可以推测, 存在一个临界$ Da_{{\rm{cr}}} $数, 当Da数达到临界值时, 互溶驱替过程中由于流固溶解反应的发生会使多孔介质中出现虫洞. 工况二中由于Da数未达到临界值, 所以驱替过程和工况一类似. 本文将在下文中对临界$ Da_{{\rm{cr}}} $数作进一步的讨论. 此外, 随着驱替的不断进行, 工况一和工况二中驱替流体仍以指状不断前进, 且随着溶质的扩散, 驱替流体与被驱替流体之间的界面逐渐增厚, 与工况三形成明显差异. 另外要注意的是, 工况二和工况三中流固溶解反应的发生同时会降低驱替流体的浓度, 减小驱替流体与被驱替流体间的黏性差异, 从而减少黏性指进的持续发展. 因此, 工况二中指进的发展较工况一略有减少.

最后($ t = 0.7324 $, 见图5(d)), 工况一和工况二中驱替流体继续沿着指进的路径不断驱离多孔介质区域, 且工况二中指进数量明显少于工况一; 而工况三中驱替流体则沿着虫洞通道驱离. 在3种工况中, 随着溶质间的扩散仍不断进行, 指进和虫洞均被不断拓宽.

为研究临界$Da_{{\rm{cr}}}$数的存在, 图6给出了$ t = $$ 0.7324 $ 时刻$ Da = 4 $和8的浓度场图. 可以看出, $ Da = 4 $时由于反应速率较小, 驱替过程中无明显的虫洞现象; 通过实验发现, 当$ Da = 8 $时, 驱替过程中可见虫洞的形成. 由此可推出, 存在临界达姆科勒数$ Da_{{\rm{cr}}} $, 当$Da\geqslant Da_{{\rm{cr}}}$时, 驱替过程中会有虫洞现象生成. 由于驱替现象会随着多孔介质结构的改变而呈现出不同的细节, 不同多孔介质下$ Da_{{\rm{cr}}} $的具体数值会有所差异. 本节各工况中的多孔介质结构见图1, 对应的$ Da_{{\rm{cr}}} \approx 8 $.

下文以$ Da_{{\rm{cr}}} $为界, 讨论不同Da数下, 多孔介质孔隙率$ \varepsilon $和驱替效率η随时间的变化情况, 其余无量纲数与4.1节保持一致. 图7(a)为$ Da = 0, 2, 4, 8 $时多孔介质孔隙率随时间的变化曲线. 可以看出, 当$ Da = 0 $, 即无溶解反应时, 多孔介质孔隙率保持不变; 当$ Da = 2, 4, 8 $, 即溶解反应较小时, 由于反应的发生, 孔隙率随时间不断增大. 且随着Da数的增大, 同时刻多孔介质的孔隙率越大. 图7(b)为$ Da = 16, 32, 160, 320 $时多孔介质孔隙率随时间的变化曲线, 此时溶解反应会生成明显的虫洞. 从图7(b)可以看出, 与Da数较小时不同的是, 孔隙率随时间的变化趋势为s型曲线. 初始阶段, 驱替流体与多孔介质固体骨架接触面积小, 溶解多孔介质骨架的量小, 多孔介质孔隙率缓慢增长; 伴随驱替流体的逐渐进入, 驱替流体与多孔介质接触面积变大, 溶解多孔介质骨架的量逐渐变多, 多孔介质孔隙率越来越大, 且孔隙率的增长也逐渐加快, 在$ t = 0.35 $左右孔隙率的变化速率达到顶峰, 此时为驱替流体初次到达流场出口的时间; 之后, 由于大量固体骨架已经被溶解, 驱替流体只能持续性拓宽虫洞周围, 因此孔隙率变化较为缓慢, 这在图5(d)中也可体现. 此外, 随着Da数的增大, 同时刻多孔介质的孔隙率越大, 这与Da数较小时情况相同. 这是因为当$ Re $数, Pe数和黏性比M保持不变时, 增大Da数等价于增大反应速率常数$ k_{\rm{r }}$. 反应速率越大, 相同时间内溶解的多孔介质骨架越多, 因此孔隙率越大. 具体地, 当$ t = 0.7324 $, $ Da = 4 $时孔隙率为0.9122, 较初始孔隙率提高了5.49%; 当$ t = $$ 0.7324 $, $ Da = 32 $时孔隙率为0.9421, 较初始时刻孔隙率提高了9.03%; 当$ t = 0.7324 $, $ Da = 320 $时孔隙率为0.9524, 较初始时刻孔隙率提高10.22%, 说明此时Da数的影响已经减弱.

为了定量描述驱替效率, 定义驱替效率

其中$ V(t) $表示在t时刻驱替流体所占的体积, $ V_{{\rm{total}}} $为多孔介质的总体积. 需要说明的是, 驱替效率应为驱替流体所占体积与流体所占总体积之比. 本文中溶解反应的发生会改变多孔介质结构, 从而改变流体所占的总体积. 为使驱替效率中所对比的量保持一致, 本文设定驱替效率为驱替流体所占的体积与多孔介质总体积之比. 此外, 由于溶质间会发生扩散, 此处定义溶质浓度$ C \geqslant 0.5 $时为驱替流体.

图7(c)为$ Da = 0, 2, 4, 8 $时驱替效率随时间的变化曲线. 该过程可以分前后两个部分进行讨论. 在$ t = 0—0.35 $时, 驱替流体并未流出多孔介质区域, 驱替效率增长较快. 由于此时Da数较小, 不同Da数下驱替效率相差不大. 当$ t > 0.35 $ 时, 部分驱替流体已流出多孔介质区域, 驱替效率在微降后缓慢增长. 此阶段, 随着Da数的增大, 驱替效率越小. 这是因为Da数较大时, 指进波及范围减少, 虫洞现象逐渐明显(如图6所示), 驱替效率越低.

图7(d)为$ Da = 16, 32, 160, 320 $时驱替效率随时间的变化曲线. 当$t = 0—0.35$时, 驱替效率同样快速增长. Da数越大, 驱替效率越小. 当$ t > 0.35 $时驱替效率会有一个小幅下降的过程. 驱替效率出现短时下降一方面是因为一部分驱替流体已经离开了流场, 另一方面是因为此时驱替流体波及范围较广, 溶解反应发生部位较多, 反应消耗的驱替流体较多, 导致流入多孔介质驱替流体的量少于流出的量加上反应的量, 驱替效率下降. 随着时间的发展, 如图5所示, 驱替流体波及范围逐渐增大, 驱替效率会逐渐变大. 另外, 在驱替过程后期, 随着Da数的增大, 驱替效率也越大(如图7(d)中插图所示). 这是因为Da数越大, 生成的虫洞相对越宽, 驱替流体波及范围相对越广, 驱替效率越大. 具体地, 在$ t = 0.7324 $时刻, $ Da = 32 $时驱替效率为0.3478, 而$ Da = 320 $的驱替效率为0.3710. 可见虽然Da数越大驱替效率相应越大, 但驱替效率增加幅度不明显.

综上所述, 当流场内未生成明显虫洞时, Da数越大, 驱替效率越小; 当流场内生成明显虫洞时, 随着时间发展, Da数越大, 驱替效率会越大, 但是仍远小于无反应时的驱替效率.

Pe数表征驱替过程中对流速率与扩散速率的相对大小, 具有重要研究意义. 模拟中设定$ Da = $$ 32 $, 通过改变Pe数来观察其对驱替过程的影响. 图8给出了$Pe {=} 5$, 65, 524, $t {=} 0.2441$时的浓度场瞬时图. 可以看出, 在Pe数很小时($ Pe = 5 $, 图8(a)), 流体的扩散能力强, 使得溶质可充分扩散, 在两流体间的界面上无法产生指进现象. 此时驱替流体以平界面形式沿驱替方向溶解固体骨架, 图中左半侧的多孔介质骨架已被溶解. 随着Pe数的逐步增大($ Pe = 65 $, 图8(b)), 流体间溶质的扩散减弱, 对流作用增强, 开始产生指进现象, 但此时指进个数较少, 多孔介质骨架被溶解的量也随之减少. 而当Pe数较大时($ Pe = 524 $, 图8(c)), 出现了明显的指进现象, 一些指进被多孔介质骨架分裂开来, 形成一系列更小的指进, 即发生了指进分裂现象. 其他研究人员也对此类现象进行了研究^[32]. 此外, 注意到此工况下多孔介质骨架被溶解的量最少.

图9为相应情形下多孔介质内水平方向上的平均浓度图. 可以看出, $ Pe = 5 $时, 平均浓度的变化十分平滑. 随着Pe数的增大($ Pe = 65 $), 浓度曲线有了细微的波动, 说明此时已有指进生成. 当Pe数为524时, 浓度曲线上出现了明显的波动, 说明高Pe数下驱替过程中产生大量指进, 流体界面不稳定性增加.

接下来选取指进现象较为明显的情况作进一步讨论. 图10为$ Pe = 131 $, 262, 524时, 多孔介质的孔隙率和驱替效率随时间的变化情况. 如图10(a)所示, Pe数越大, 同一时刻下孔隙率越小. 这是因为Pe数越大, 扩散的作用越小, 指进现象越明显, 指进的增长速度也越快, 溶解的多孔介质骨架越少, 孔隙率越小. 此外, 不同Pe数下孔隙率的变化趋势相同, 均为s型曲线. 如图10(b)所示, 当Pe数越大时, 初始阶段驱替效率越大. 这是因为Pe数越大, 同时刻下指进数量越多, 驱替流体波及范围越广(如图10(b)中插图①和②所示), 驱替效率相应越大. 但随着时间的发展, Pe数越小, 虫洞越宽, 且入口和虫洞附近溶解的固体骨架越多(如图10(b)中插图③和④所示), 驱替流体占据更多的空间, 驱替效率越大.

在上文的模拟中, 所选取的多孔介质均如图1所示. 为验证上述结果具有普适性, 本节选取同一多孔介质CT扫描灰质图中的其他部分进行二值化处理, 所获取的多孔介质初始孔隙率分别为$ \varepsilon_0 = 0.8 $和0.75, 并分别模拟了两种多孔介质中含溶解反应的驱替过程. 图11为两种多孔介质下孔隙率和驱替效率随时间的变化图. 可以看出, 孔隙率的变化虽然会显著影响孔隙率和驱替效率的数值大小, 但并不会改变整体的变化趋势, 即孔隙率变化依然呈s型, 驱替效率随时间先增大然后减小最后增大. 这表明上文研究的驱替过程的基本规律在不同的孔隙率下仍然适用, 研究结果不失一般性.

此外, 如图11(b)所示, 随着多孔介质初始孔隙率的增加, 驱替效率也相应增加. 这是由于同一多孔介质中, 初始孔隙率越大的部分, 流体越容易通过多孔介质, 为驱替的进行创造更好的条件, 从而为溶解反应在流固界面的进行也提供了更多的机会, 使得驱替效率随之增加. 因此不同多孔介质结构对含溶解反应的驱替有着显著影响, 初始孔隙率较大的多孔介质中驱替效率也越高.

本文采用LBM模拟了孔隙尺度下多孔介质内含有溶解反应的互溶驱替过程, 观察了含流固溶解反应时多孔介质内流体浓度场随时间的演化情况, 并与不含溶解反应的工况进行对比, 研究溶解反应对驱替过程的影响; 之后通过改变Da数和Pe数研究孔隙率和驱替效率的变化情况, 得到如下结论.

1)在两流体黏性差异较大的前提下, 在驱替过程早期, 溶解反应会消耗驱替流体的浓度, 抑制指进的继续发展; 随着时间的增加, 溶解反应会不断溶解多孔介质骨架, 在多孔介质内部生成虫洞, 驱替流体会直接沿着该虫洞离开多孔介质, 减小了驱替效率.

2)随着Da数的增大, 反应速率变大, 孔隙率变化越大, 所生成的虫洞越宽, 驱替效率逐渐增大.

3)随着Pe数的增大, 指进增长的速率越快, 反应的量越小, 孔隙率变化越小, 驱替效率也随之变小.

4)不同的多孔介质结构下获得的孔隙率和驱替效率随时间的变化趋势一致, 说明以上结论具有普适性. 此外, 初始孔隙率较大的多孔介质中驱替效率也越高.

在互溶驱替过程中, 若两种流体黏性差异较大且无法避免流固溶解反应发生的情况下, 为获取较高的驱替效率, 应选择与多孔介质发生溶解反应速率高的驱替流体, 但无需一味追求反应速率. 此外, 系统的Pe数越小驱替效率越高. 为减缓驱替, 可采取相反的策略.