0 引言

多孔介质中伴随化学反应的流体混合过程是渗流力学中的典型问题，普遍存在于化工生产和自然界等领域^{[1, 2, 3, 4]}，如生化制药^[2]、石油及天然气开采^[3]、二氧化碳埋存^[4]等.因此对多孔介质中反应流体混合过程的机理研究不但能促进渗流力学理论体系的发展，同时也能指导人们在实际应用中对此过程的预测和控制.

多孔介质中的流动可在孔隙(pore)尺度和表征体元(representative elementary volume,REV)尺度上进行描述^{[5, 6]}.由于多孔介质结构复杂，目前大量文献都基于REV尺度来研究流体的混合过程.对于等粘性的两流体，相关研究^{[7, 8, 9]}表明速度扰动通过产生紊乱促进其混合，分子扩散则会将该紊乱转变成流体间的有效混合.而对于粘性不同且存在粘性指进时的流体混合过程，Jha等人^{[10, 11]}提出了预测标量耗散率的数学模型，并定量分析了粘性指进与混合之间的关系.他们指出标量耗散率随时间的演化过程非单调，且粘性指进能通过增大流体接触面以及产生速度扰动促进流体混合.以上研究仅限于界面无化学反应的情况，当混溶系统中引入界面反应后，现有的研究^{[12, 13, 14, 15]}大都针对化学反应对粘性指进形态的影响.如Hornof等人^[12]研究发现粘性指进在化学反应条件下会发生明显变化，此后Sastry等人^[13]对此开展了进一步研究，发现反应能改变流体表面张力并对指进形态产生显著影响.对于互溶流体，Wit等人^{[14, 15]}在给定粘性比条件下研究了非线性化学反应对指进的影响，发现反应通过改变流体粘性使得指进的形态特征以及传播机制发生显著变化.可以看出界面反应会改变流体的物理性质(如粘性、密度、表面张力等)，进而影响指进形态特征.而指进的改变会进一步引起流体混合程度的改变，但是，目前直接研究化学反应对流体混合过程影响的报道却不多见.

应当指出的是，上述研究都是基于REV尺度的渗流模型.尽管研究成果有助于工程实际应用，但难以准确描述多孔介质内部化学反应流体的混合过程和流动特性.近几十年发展起来的格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)^{[8, 9]}，可以使用无滑移反弹格式处理复杂壁面，因此成为从孔隙尺度上研究流动及混合问题的有力工具.例如刘^[16]等人在孔隙尺度上对各向异性多孔介质中的流体混合过程进行定量分析，指出渗透率各向异性对纵向和横向标量耗散率的影响显著.在孔隙尺度上对流动过程进行数值模拟，需要使用很高的网格分辨率，导致计算量巨大.近几年提出的结合GPU的格子玻尔兹曼(Lattice Boltzmann,LB)模拟加速方法，可以大幅提高计算速度，为研究复杂渗流问题提供了可行途径.如黄等人^[17]系统分析了LBM的GPU优化性能，朱等人^[18]在GPU上实现了基于LBM的孔隙尺度多孔介质流动模拟加速.上述研究表明基于GPU的LBM可作为研究孔隙介质中流动问题的有效方法.本文结合LBM和GPU加速技术，从孔隙尺度上研究多孔介质中反应流体的粘性指进现象，并定量分析界面反应对流体混合的影响.

1 问题描述

考虑各向同性多孔介质中不同粘性反应流体间的混合过程.系统初始设置如图 1所示：在尺寸为l×l的多孔介质中，以团状散布着浓度c₀=0动力粘度为μ₀的反应溶剂a，且被浓度c₁=1粘度为μ₁＜μ₀的反应溶剂b包围.两流体各占总体积的一半，在被水平方向外力F_x驱动的过程中，接触面上会发生化学反应.此处各向同性多孔介质由四参数法(Quartet Structure Generation Set,QSGS)^[19]生成.

上述问题可由不可压Navier-Stokes (NS)方程和含源对流扩散(convection-diffusion-reaction,CDR)方程描述：

$\begin{array}{c} \nabla \cdot u = 0,\\ \frac{{\partial u}}{{\partial t}} + u \cdot \nabla u = - \nabla p + \upsilon \left( c \right){\nabla ^2} + F,\\ \frac{{\partial c}}{{\partial t}} + u \cdot \nabla c = D{\nabla ^2}c + f\left( c \right), \end{array}$

(1)

式中u、p、υ、D、F分别为流体的速度、压力、运动粘度、扩散系数和外力.混合过程中流体运动粘度为υ=μ/ρ，而动力粘度μ与浓度c遵循指数混合法则^[20]，即μ(c)=μ₀exp(-cR)，其中流体的粘性比R=ln(μ₀/μ₁)，方程组(1)也由此耦合.由于仅在界面发生化学反应，根据文献[14]选取单浓度变量的双稳态化学反应模型，其源项表达式为f(c)=-k(c-c₀)(c-c₁)(c-c₂)，其中k为化学反应速率参数，c₂为稳态浓度.模拟时边界条件如下给定：四周均为周期边界，而多孔介质中的流固边界为u=(u,v)=(0,0).分别选取L=l、U=l/l/F_x、C=c₁为特征长度、特征速度和特征浓度，可得到表征上述问题的两个无量纲参数，即雷诺数Re=Ul/υ和贝克莱数Pe=Ul/D. 2 LBM模型及其GPU实现

采用双分布耦合LB模型求解方程组(1)，使用多松弛格子Boltzmann (MRT)模型^[6]来解流动方程.CDR方程则使用Inamuro等人^[21]提出的Lattice Kinetic Scheme(LKS)，并结合邓等人^[22]提出的含源对流扩散的单松弛格子Boltzmann(LBGK)模型来求解，下面分别予以介绍.

MRT模型在碰撞过程中使用多个松弛时间，能有效克服渗透率与粘性相关的缺点^[23]并且具有良好的数值稳定性^{[24, 25]}.其演化方程

${f_i}\left( {x + {c_i}\delta t,t + \delta t} \right) - {f_i}\left( {x,t} \right) = - {M^{ - 1}}S\left( {m\left( {x,t} \right) - {m^{{\rm{eq}}}}\left( {x,t} \right)} \right) + \delta t{F_i},i = 0,1,\cdots ,q - 1,$

(2)

式中f_i(x,t)是t时刻x处离散速度为c_i的粒子分布函数，q是离散速度个数，δt为时间步长，F_i为外力.m和m^eq分别表示分布函数在矩空间对应的矩及其平衡态，M为q×q的转换矩阵，可将速度空间分布函数映射到矩空间：m=M·f. S为对角矩阵，其元素为相应矩的松弛率.

本文采用二维九速(D2Q9)模型^[6]，离散速度为c₀=(0,0)e，c₁=(1,0)e，c₂=(0,1)e，c₃=(-1,0)e，c₄=(0,-1)e，c₅=(1,1)e，c₆=(-1,1)e，c₇=(-1,-1)e，c₈=(1,-1)e，其中e=δx/δt为格子速度，本文取值为1，δx为空间步长.该模型的变换矩阵为

$M = \left[{\begin{array}{*{20}{c}} 1&1&1&1&1&1&1&1&1\\ { - 4}&{ - 1}&{ - 1}&{ - 1}&{ - 1}&2&2&2&2\\ 4&{ - 2}&{ - 2}&{ - 2}&{ - 2}&1&1&1&1\\ 0&1&0&{ - 1}&0&1&{ - 1}&{ - 1}&1\\ 0&{ - 2}&0&2&0&1&{ - 1}&{ - 1}&1\\ 0&0&1&0&{ - 1}&1&1&{ - 1}&{ - 1}\\ 0&0&{ - 2}&0&2&1&1&{ - 1}&{ - 1}\\ 0&1&{ - 1}&1&{ - 1}&0&0&0&0\\ 0&0&0&0&0&1&{ - 1}&1&{ - 1} \end{array}} \right].$

(3)

相应的矩为m=(ρ,e,ε,j_x,q_x,j_y,q_y,p_xx,p_yy)^T，其中T表示转置变换.ρ是流体密度，e_e和ε分别与能量和能量的平方有关，q_x、q_y以及j_x、j_y分别对应x、y方向上的能量通量和动量，p_xx和p_xy则表示应力张量的对角线和非对角线上的元素.对应的平衡态矩为

${m^{{\rm{eq}}}} = {\left( {\rho ,- 2\rho + 3{\rho _0}{u^2},\rho - 3{\rho _0}{u^2},{\rho _0}u,- {\rho _0}u,{\rho _0}v,- {\rho _0}v,{\rho _0}\left( {{u^2} - {v^2}} \right),{\rho _0}uv} \right)^{\rm{T}}},$

(4)

其中u、v分别为速度矢量u在x、y方向的分量，ρ₀是流体密度，本文取值为1，该模型的松弛参数为

$\boldsymbol{S}{\rm{ = }}\left( {{s_{\rm{\rho }}},{s_{\rm{e}}},{s_{\rm{\varepsilon }}},{s_{\rm{j}}},{s_{\rm{q}}},{s_{\rm{j}}},{s_{\rm{q}}},{s_{\rm{\upsilon }}},{s_{\rm{\upsilon }}}} \right).$

(5)

由于碰撞过程中流体密度和动量守恒，本文可给定s_ρ=s_j=0.根据文献[23]给出其余项为：s_υ=s_e=s_ε=1/τ，sq=(16-8s_υ)/(8-s_υ)，其中松弛时间τ与运动粘度υ相关

$\upsilon = c_{\rm{s}}^2\left( {\tau - \frac{1}{2}} \right)\delta t,$

(6)

式中c_s²=e²/3是格子声速，根据文献[26]对外力项定义

$F = {M^{ - 1}}\left( {I - \frac{1}{2}S} \right)M\bar F,$

(7)

其中$M\bar F$为矩空间的外力项，计算表达式为

$M\bar F = {\left( {0,6u \cdot F,- 6u \cdot F,{F_x},- {F_x},{F_y},- {F_y},2\left( {u{F_x} - v{F_y}} \right),u{F_y}v{F_x}} \right)^{\rm{T}}}.$

(8)

下面给出流体密度ρ和速度u的计算表达式

$\rho = ,{\rho _0}u = \sum\limits_{i = 0}^8 {{c_i}{f_i} + } \frac{{\partial t}}{2}F.$

(9)

此外，扩散过程通过LKS求解，并根据文献[22]引入化学反应源项，浓度分布函数g_i(x,t)的演化方程如下

${g_i}\left( {x + {c_i}\delta t,t + \delta t} \right) - {g_i}\left( {x,t} \right) = - \frac{1}{{{\tau _{\rm{D}}}}}\left( {{g_i}\left( {x,t} \right) - g_i^{{\rm{eq}}}\left( {x,t} \right)} \right) + \delta t{F_{{\rm{rei}}}}\left( {x,t} \right) + \frac{1}{2}\delta {t^2}\frac{{\partial {F_{{\rm{rei}}}}\left( {x,t} \right)}}{{\partial t}},$

(10)

方程中τ_D是与扩散相对应的松弛时间，平衡态分布函数g_i^eq以及化学反应分布函数F_rei(x,t)分别为

$g_i^{{\rm{eq}}} = c{\omega _i}\left( {1 + \frac{{{c_i} \cdot u}}{{c_{\rm{s}}^2}}} \right) + {\omega _i}A\delta t\left( {{c_i} \cdot \nabla c} \right)$

(11)

和

${F_{{\rm{rei}}}}\left( {x,t} \right) = {\omega _i}f\left( c \right)\left( {1 + \frac{{{c_i} \cdot u}}{{c_{\rm{s}}^2}}\frac{{\left( {{\tau _{\rm{D}}} - 0.5} \right)}}{{{\tau _{\rm{D}}}}}} \right),$

(12)

其中ω_i是权系数：ω₀=4/9，ω_1-4=1/9，ω_5-8=1/36.参数A和扩散系数相关，且在低Mach数条件下，对方程(10)做Chapmann-Enskog多尺度展开可得

$D = c_{\rm{s}}^2\left( {{\tau _{\rm{D}}} - \frac{1}{2} - A} \right)\delta t,$

(13)

以及浓度梯度的计算表达式

$\nabla c = \frac{{\sum {{c_i}{g_i} - \sum {{c_i}g_i^{{\rm{eq}}}} } - {\tau _{\rm{D}}}\delta t\sum {{c_i}{F_{{\rm{rei}}}}} }}{{c_{\rm{s}}^2\delta t\left( {A - {\tau _{\rm{D}}}} \right)}}.$

(14)

在标准LB模型中，当τ_D趋近于0.5时，会出现数值不稳定现象.而此处可通过调整参数A保证模型在D较小(即Pe较大)时的数值稳定性，本文取τ_D=1，该模型的浓度定义如下

$c\left( {x,t} \right) = \sum\limits_{i = 0}^8 {{g_i}\left( {x,t} \right).} $

(15)

边界处理是LBM 的一个关键问题，对于本文多孔介质中固体的速度边界，采用半反弹格式^[27-28]来实现的无滑移边界条件

${g_{\bar i}}\left( {{x_{\rm{f}}},t + \delta t} \right) = {f'_i}\left( {{x_{\rm{f}}},t} \right) - 2{\omega _i}\rho \frac{{{c_i} \cdot {u_{\rm{w}}}}}{{c_{\rm{s}}^2}},$

(16)

其中x_f是临近固体的流体格点，c_i=-c_i(c_i指向壁面)，f'_i(x_f,t)为f_i(x_f,t)碰撞后的值，u_w为壁面速度.当给定壁面浓度c_w时，类似给出固体浓度边界

${g_{\bar i}}\left( {{x_{\rm{f}}},t + \delta t} \right) = - {{g'}_i}\left( {{x_{\rm{f}}},t} \right) - 2{\omega _i}{c_{\rm{w}}}.$

(17)

在LBM计算过程中，每个格点有多个分布函数，但其碰撞及迁移过程的计算却相对简单且具有良好的局部性，因此影响程序性能的瓶颈主要在于数据的访存.基于这一特点，本文在GPU上对LBM进行优化的主要目标为利用GPU的带宽优势，提高访存效率.CUDA是由NVIDIA推出的通用并行计算构架，通过Grid-Block-Thread的三层次结构，将计算任务映射到GPU的大量处理单元上执行并行运算.以C语言为例，基于CUDA的程序计算基本步骤如下：首先，分别在GPU和CPU上分配内存，并初始化相关数据.其次，设置Grid和Block的大小，在GPU上进行并行计算.最后，当计算结果满足终止准则时，将数据拷贝到CPU上，得到最终结果，并释放GPU和CPU的内存.

在程序中，为实现连续、合并访问GPU，应保证相邻格点同一离散速度的分布函数在GPU中的位置连续.因此，我们采用形式为f(i×NX×NY+y×NX+x)、F(i×NX×NY+y×NX+x)的两套一维数组来存储坐标为(x,y)的格点沿i方向的分布函数，其中NX×NY表示网格大小.在执行碰撞迁移过程时，对于奇数时间步，每个格点从f中加载分布函数，经碰撞流动后存储在F中；对于偶数时间步则反过来执行.下面将对模型的可行性以及优化性能进行测试.

3 模型验证

为验证模型的可行性，我们在R=0、Pe=1及Re=0.32(大粘性流体)时，对Taylor-Aris弥散问题^{[29, 30]}进行计算.该问题中两流体的粘性和密度均相同，且界面无化学反应，可由NS方程和CD方程描述，即方程组(1)中f(c)=0.此处选取网格数为32×256以及U=0.005，ρ₀=1，l_y=32(以上物理量均为格子单位).在Pe足够小时，截面平均浓度c_A(x,t)有解析解

${c_{\rm{A}}}\left( {x,t} \right) = \frac{1}{2}\left[{1 + {\rm{erf}}\left( {\frac{{{x_0} - x}}{{2\sqrt {Dt} }}} \right)} \right],$

(18)

其中${\rm{erf}}\left( z \right) = \int_0^z {{{\rm{e}}^{ - {\eta ^2}}}{\rm{d}}\eta } $是误差函数，x₀为c(x,t)=0.5界面的平均位置.模拟时，截面平均浓度c_A(x,t)的计算如下

${c_{\rm{A}}}\left( {x,t} \right) = \left( {\int_0^H {c\left( {x,y,t} \right){\rm{d}}y} } \right)/{l_y}.$

(19)

图 2给出了计算时间为t=100、400、800(格子单位)时，计算结果与解析解的对比.可以发现二者吻合良好，验证了上述耦合LB模型的正确性.

我们也测试了GPU加速算法的性能，计算环境为Intel xeon E5-2670 CPU，GPU为Nvida Tesla k20.CUDA toolkit版本为 6.0，操作系统为red hat 6.2，程序使用CUDA扩展C语言.我们模拟了多孔介质中粘性不同的互溶流体间的混合过程.该问题可由NS方程和对流扩散方程来描述，即方程组(1)中f(c)=0.系统初始位置如图 1所示，流体在尺寸为1×1的多孔介质内做周期运动.其中多孔介质由QSGS方法生成，给定固体基质各向生长概率为0.005.根据文献[16]选取网格数为800×800、R=3、Pe=1×10⁵、F_x=2.5×10^-3、Re=120(大粘性流体)及孔隙率P_o=0.8.模拟时根据特征速度U=F_xl可得特征时间T=l/U以及无量纲时间t=δtN_step/T，此处取计算步数N_step=1×10⁶.大粘性流体的运动粘度υ₀和扩散系数D由给定的Re和Pe确定，即υ₀=Ul/Re及D=Ul/Pe.

为了定量表征计算过程中流体的混合，我们定义混合度为

$\chi = 1 - {\sigma ^2}/\sigma _{\max }^2,$

(20)

其中σ²=〈c²〉-〈c〉²为浓度方差，符号〈·〉表示对计算区域的空间平均.最大方差σ_max²则对应两流体完全独立的状态，当〈c〉=0.5时σ_max²=0.25.流体混合程度越高，浓度场的整体方差σ²越小，混合度χ越大；当流体完全混合时σ²=0，χ=1.我们还定义了标量耗散率ε来表征流体的混合率

$\varepsilon = D\left\langle {{{\left| {\nabla c} \right|}^2}} \right\rangle .$

(21)

由上式可以看出，ε与流体的浓度梯度$\nabla c$相关：当D固定时，较小的ε值意味着整个空间内流体的浓度分布比较均匀.

图 3给出了GPU和CPU程序计算结果的对比，可以看出二者关于ε和χ的结果都吻合良好.表明本文所采用的基于GPU的耦合LBM对模拟多孔介质中流体的混合过程是可靠的.此外我们对相同条件下二者的计算效率进行比较，计算速度衡量标准为每秒百万格点更新速率(MFLUPS)^[17]，其计算表达式为

$V = {N_{{\rm{step}}}}{N_{\rm{f}}}/\left( {t \times {{10}^6}} \right),$

(22)

其中N_f为流体格点，计算结果如表 1所示.

可见相对于CPU平台，使用GPU加速算法的加速比可达34，很大程度地提高了计算效率.

4 计算分析

基于上述LBM-GPU方法，我们对多孔介质中包含界面化学反应的流体混合过程进行数值模拟和分析.计算模型及边界条件与前文一致，并选取相关物理量：R=3、Pe=1×10⁵、F_x=2.5×10^-3、ρ=1、l=1及Re=120(大粘性流体).本文考虑对称且速率不同的化学反应对混合过程的影响，即给定稳态浓度c₂=0.5，并在一定范围内调节反应速率参数k.

4.1 网格无关性验证

为确定能保证精度和计算效率的最佳网格，首先进行网格无关性验证.此处分别使用400×400，800×800和1 600×1 600的网格来离散图 1中的计算区域.分别计算了k=0.02、0.04及孔隙率P_o=0.8时的混合过程，并测试了计算区域ε和χ随时间的演化，结果如图 4所示.

由图 4可以观察到，网格数为400×400的计算结果与两个较细网格的结果差异明显，而800×800与1 600×1 600这两种网格下的结果彼此间变化趋势一致且吻合良好，表明800×800的网格能够保证计算精度且可以节省计算量，因此以下的数值模拟均采用该网格.

4.2 孔隙率的变化

由于工程应用中不同多孔介质的孔隙率变化范围广，而孔隙率的变化会影响流体在多孔介质中的流通.因此在上述选取的网格下，还需确定本文计算采用的孔隙率.我们模拟计算了孔隙率P_o=0.4、0.6、0.8及k=0.00时的流体混合过程，并测试计算区域ε和χ随时间的演化，结果如图 5所示.

由图 5可以观察到，不同孔隙率下ε和χ的演化过程均可分为非单调的四个阶段，与文献[16]报道的结果一致.该结果表明孔隙率的改变对流体混合的具体过程影响显著，但流体混合整体趋势的变化却不大.此外我们还发现孔隙率P_o=0.8时各阶段区分明显，且孔隙率越大流体越容易流通，因此本文以下的数值模拟选取P_o=0.8进行计算.

4.3 反应流体的混合过程

现在研究反应速率k对反应流体粘性指进和混合过程的影响.图 6为工况k=0、0.02、0.04下t/T=0.5、5、30、60时的浓度场分布，反映了不同反应速率时多孔介质中流体混合的具体情况.

如图 6(a)所示，k=0(即无界面反应)时，流体主要经历扩散过渡、指进形成和快速混合的过程，且最终流体将均匀分布在多孔介质中，与文献[16]报道的结果一致.与之对比，当存在化学反应时，由图 6(b)和6(c)可以发现，在扩散过渡和指进形成阶段(即t/T=0.5、5)，流体界面较薄，化学反应作用极小，引入反应前后的浓度分布基本不变.随着指进的重叠以及扩散的继续，流体混合区域(0＜c＜1)迅速增加.由源项表达式f(c)=kc(c-1)(c-0.5)可知，反应使得混合区域中大于稳态浓度(即c=c₂=0.5)的区域内浓度增加，而小于稳态浓度的区域内浓度减小.从而导致流场中两流体界面变薄，混合区域减少，抑制流体混合.同时还发现，反应使得浓度c=1的小粘性流体所占比例增加，且随着反应速率的增加，作用越明显.为进一步定量描述化学反应的影响，下面我们将分别分析平均浓度、标量耗散率、混合度等参数的变化规律.

1) 平均浓度〈c〉

在连续的计算时间内对二维的浓度分布c(x,y,t)求平均，可得到计算区域平均浓度〈c〉关于时间的演变过程，计算表达式为

$\left\langle c \right\rangle \left( t \right) = \frac{1}{{l \cdot l}}\int_0^l {\int_0^l {c\left( {x,y,t} \right){\rm{d}}x{\rm{d}}y,} } $

(23)

其大小表征了两流体所占比例，〈c〉越大表示c=1的小粘性流体所占比例越大.其变化趋势反应了化学反应对两流体比例的影响.

计算结果如图 7所示，可以看出当k=0时，曲线维持〈c〉=0.5不变，表明在周期边界的条件下，流体混合过程中两流体比例保持不变.引入化学反应后，在扩散过渡和指进发展阶段，反应作用不明显，两流体比例不发生改变.如图中所示曲线在初期都能维持〈c〉=0.5，这与图 6中t/T=0.5、5时的浓度分布保持一致.随后指进引起的紊乱使流体快速混合，界面开始变厚，化学反应将显著影响混合区域两流体所占比例.根据文献[15]可知，在本文所考虑的对称反应下，团状散布的流体最终将被包裹它的流体消耗，即此处化学反应会使c=1的小粘性流体所占比例增加.如图中k>0时的曲线所示，〈c〉随着时间的增加而增加且最终维持在最大值〈c〉_max处，且随反应速率越大，〈c〉增大越明显且最终〈c〉_max也越大.同样由图 6中t/T=30、60时的结果也可以看出相同的浓度变化趋势.

2) 标量耗散率ε

为量化多孔介质中流体的混合率，我们测量了不同反应速率条件下表征流体混合率的标量耗散率ε，如图 8所示.由公式(20)可知ε反映了计算区域浓度分布的均匀程度，浓度不均匀区域Λ以及其中的浓度梯度$\nabla c$越小，ε的数值越小，流体混合程度越高.

由图(8)可以观察到在扩散过渡和指进形成阶段，不同工况下ε随时间的变化规律一致：先减小到局部极小值，然后增大到局部极大值.这是因为此时界面较薄，化学反应作用不明显.根据图 6中t/T=0.5时的浓度分布可以看到，在扩散过渡阶段，Λ面积变化不大.而$\nabla c$随着扩散而明显减小，导致ε下降.随着粘性指进的发展，如图 6中t/T=5时的结果所示，Λ迅速增加.即使扩散作用使$\nabla c$变小，但其程度远小于Λ的增加幅度.因此ε会从局部极小值开始上升，第一次呈现“反混合”现象.

当指进发展到一定程度，混合带开始重叠，如图 6中t/T=30时，Λ增加幅度变小.此前指进引起的扰动，促进大面积的扩散使得流体不断均值化，$\nabla c$逐渐减小，系统进入快速混合的前期.在指进和扩散的作用下，ε将会从局部最大值开始下降，如图中k=0.00的曲线所示.比较图 6可以发现，此时反应将显著作用于流体混合区域，导致界面明显变薄且c=1的小粘性流体增加.从而使得Λ明显减少而$\nabla c$显著增加，抑制ε的下降程度，且反应速率越大越明显，如图中k=0.02、0.04的曲线所示.

在快速混合的后期，如图 6中t/T=60时，当c=1的小粘性流体增加到一定程度后，c=0的大粘性流体已经很少.Λ减少程度将会下降，而$\nabla c$仍会增加，这使得ε出现小幅度的回升，再次出现“反混合”现象.在短暂的“反混合”过程后，化学反应会使得流体达到分离的极限，即$\nabla c$增加不明显，而Λ持续减少导致ε继续下降，如图中k=0.02所示.而反应速率很大时，反应使得c=1的流体迅速增加，几乎布满计算区域，不会出现“反混合”现象，ε会迅速下降，如图中k=0.04所示.

3) 混合度χ

此外我们还测试了不同工况下的混合度χ，以此来量化流体的浓度方差.与ε的演化对应，不同反应速率条件下，χ在扩散过渡和指进发展阶段也能良好吻合，如图 9所示.同时可以看出此过程中χ持续增加，即指进发展阶段的“反混合”过程实质上是促进混合的.而由图 6中t/T=0.5、5时的浓度分布也可以看出，此过程中两流体持续混合.

此后，系统开始快速混合并使得χ持续增加，而由图 6(b)和6(c)可知，此时化学反应将抑制流体混合，χ的增加幅度减小.如图 9中k=0.02、0.04的曲线所示，反应速率越大抑制作用越强，χ的增加也幅度越小.在快速混合后期，k=0.02时曲线会出现小幅度的下降，表明此时化学反应抑制流体混合，出现了本质上的反混合过程.而k=0.04时曲线迅速上升，由图 6(c)中t/T=60时的浓度分布可以看出，c=1的小粘性流体已经几乎布满整个区域，所以不会再出现反混合现象.

5 结论

采用格子Boltzmann方法并结合GPU计算技术，在孔隙尺度上数值模拟了多孔介质中反应流体的粘性指进.定量分析了不同反应速率界面反应对混合过程的影响，结论如下：

1) 在扩散过渡和指进发展阶段，流体间界面较薄，化学反应作用极小，不同反应的混合过程差别不大，统计的平均浓度〈c〉、标量耗散率ε和混合度χ随时间变化的趋势一致.

2) 随着指进的重叠，在其引起的紊乱下流体开始快速扩散，流体愈加混合.此时化学反应开始显著作用于混合过程：首先，反应导致不同粘性流体的比例发生变化，以团状散布的小粘性流体逐渐被消耗，使得〈c〉增加.其次，通过测量ε以及χ随时间的演化，可以发现反应能抑制混合过程，甚至会出现反混合现象.并且随着反应速率的增加，上述影响越剧烈.