论文部分内容阅读
摘 要：采用流体力学格子Boltzmann方法分别建立酸化中流動过程模型及传质过程模型。通过速度场对两个过程进行耦合，实现对酸化过程的微观尺度数值模拟。结果表明：①格子Boltzmann方法能实现对酸化过程的微观模拟，为研究碳酸盐岩油气藏酸化提供新手段;②孔隙尺寸的不均匀分布是造成碳酸盐岩非均匀溶蚀的重要原因。
关键词：酸化;格子Boltzmann;非均匀溶蚀;流动;传质
酸化是解除碳酸盐岩储层近井筒污染，提高油气井产量的重要措施[1-6]。大量学者通过实验手段及数模手段对宏观尺度上蚓孔扩展规律进行研究[7-10]，认为非均匀溶蚀是形成蚓孔的重要因素。本文利用格子Boltzmann方法，研究微观尺度上盐酸与碳酸盐岩非均匀溶蚀问题，证实格子Boltzmann方法对酸化过程模拟的可行性，为酸化数值模拟研究向微观方向发展提供新思路。
1 流动过程
1.1 流动过程的D2Q9模型（图1）
格子Boltzmann方程包含两部分：迁移和碰撞。式（1）为LB方程，表征流体流动过程，即分布函数的演化过程：
[fαx+eαΔt，t+Δt-fαx，t=-1τf αx，t-feqαx，t]…（1）
式中：[fα]为[α]方向的粒子速度分布函数;[Δt]为时间增量;[τ]为弛豫时间;[feqα]为相应的平衡态分布函数。
上式可分解为两个过程，即碰撞过程和迁移过程。
碰撞过程：
[fαx，t+Δt=fαx，t1-1τ+1τfeqαx，t]
迁移过程：
[fαX+Δx，t+Δt=fαx，t+Δt]
对于D2Q9模型，松弛时间与运动粘度关系为：
[v=Δx2c2sΔtτ-12]
为保证运动粘度为正值，弛豫时间的选取大于0.5，[c2s]为格子声速的平方。
平衡态分布函数计算式：
[feqαX=wαρX1+eα?uc2s+eα?u22c4s-u22c2s]…（2）
式中：[ρX]为[X]处宏观密度;[u]为宏观速度。
1.2 MRT模型
上述介绍的基于单松弛时间的LBGK模型（SRT）是最简单的格子玻尔兹曼模型，由于简单，LBGK模型为最流行的格子玻尔兹曼模型。
SRT模型的两个不足：一是数值稳定性差，二是会造成粘性相关的边界。一旦模型用于模拟多孔介质中流动时，可能导致计算得到的渗透率与粘性相关，这有悖于渗透率的物理意义。为克服BGK模型的上述不足，引入MRT模型模拟多孔介质内的流体流动[11]。
MRT模型在碰撞中使用多个松弛参数，在矩空间执行碰撞过程，在速度空间执行迁移过程。变换矩阵M[12]：
[ M=111111111-4-1-1-1-122224-2-2-2-21111010-101-1-110-20201-1-110010-111-1-100-20211-1-101-11-1000000000111-1]
[m]为速度矩矢量，[f]为分布函数，[meq]为平衡态速度矩矢量。关系：[m=Mf]。
[m=ρ，e，ε，jx，qx，jy，qy，pxx，pxyT]
[meq=ρ，-2ρ+3u2x+u2y，ρ-3u2x+u2y，]
[ux，-ux，uy，-uy，u2x-u2y，uxuyT]
各矩的物理意义，其中，[e]为能量;[ε]为能量平方;[jx，jy]分别为[x]和[y]方向的质量流量;[qx]，[qy]分别为[x]和[y]方向的能量流量;[pxx，pxy]分别为黏度应力张量对角和非对角线部分;[ux，uy]分别为[x]和[y]方向的宏观速度。
松弛矩阵[S]为对角矩阵，
[S=diags0，s1，s2，s3，s4，s5，s6，s7，s8];其中，[sv=1/τ]，用于确定剪切黏度系数。不同的矩可采用不同松弛时间，这是符合物理规律的。当[si=1/τi=0，1，2，L，8]时，MRT模型退回为LBGK模型。
多松弛模型中不同松弛参数组合对结果有较大影响，通过比较给出模拟多孔介质时最佳参数组合，即[s0，s3，s5]取0，[s1，s2，s7，s8]取[1/τ]，[s4，s6]取[82τ-18τ-1]，本文松弛矩阵的选取就采用此种组合。
1.3 边界条件处理
为将格子Boltzmann方法应用于实际问题，还要加上一定的边界条件，本文使用的边界条件有反弹边界、Zou-He边界、非平衡外推边界[13-15]。
1.3.1 反弹边界
反弹边界是LBM方法中最简单的边界处理法，处理复杂几何边界时，如多孔介质中流动时有着广泛应用。它将边界看成具反弹边界条件的边界，分布函数看成运动的刚性粒子，粒子于边界发生碰撞后，以原速度反方向返回，其值大小不变。表达式为：
[fα，Xf，t+Δt=fαxb，t]
[α，]表示[α]的反方向，[Xf]表示临近边界的流体格点，[xb]表示边界格点。
1.3.2 Zou-He速度压力边界
速度边界 实际应用中，最常见的是边界给定某一速度。图2是边界函数分布图。
以North边界为例，北部边界，未知分布函数为[f4，7，，8]，密度[ρ]，由宏观密度方程和动量方程可得：
[ρ=f0+f1+f3+2f2+f5+f61+v0]…（3）
[f4=f2-23ρv0]…（4） [f7=f5+12f1-f3-16ρv0]…（5）
[f8=f6+12f1-f3-16ρv0]…（6）
由此可得未知分布函数[f4，7，，8]及密度[ρ]，其他3个方向的边界条件求解过程同上。
1.3.3 非平衡外推边界
边界处分布函数分为平衡与非平衡两部分，其中非平衡部分与相邻节点处的非平衡部分近似，可用于处理速度与压力边界，基于此种格式处理的速度和压力边界有二阶精度。
假设边界格子节点为[O]，与之相邻的内部节点为[G]，[O]点处速度[u]已知，密度[ρ]未知，[fi]为未知方向分布函数。节点[O]处粒子分布函数由平衡态和非平衡态两部分组成，即：
[fiO，t=feqiO，t+fneqiO，t]…（7）
确定节点[O]处的平衡态粒子分布函数所需的密度可用格子节点[G]处的流体密度代替，即：
[feqiO，t=feqiρG，t，t]…（8）
对非平衡态部分[fneqiO，t]，则用邻近流体格子节点[G]处的非平衡态分布函数近似得到，即：
[fneqiO，t=fneqiG，t=fiG，t-feqiG，t]…（9）
该处理方法的的特点是边界格子节点上所有粒子分布函数值全重新赋值，数值稳定性较好（图3）。
2 传质过程
2.1 传质过程的D2Q4 模型
与流动过程类似，LBM也被用于模拟传质过程。由于控制方程中只有质量守恒方程，因此模拟中需要的速度方向更少。与D2Q9模型相比，在保证精度和稳定性前提下，D2Q4模型具计算速度快、消耗内存少等优点。本文采用D2Q4模型对传质过程进行模拟。
碰撞迁移过程：
[gαx+dx，t+dt-gαx，t=-1τgαx，t-geqαx，t]…（10）
式中：[t]为时间，[dt]时间增量，[dx]位置增量;[eα]为[α]方向发格子速度，矢量;[gα]为[α]方向的分布函数;[geqαx，t]为位置x处、时间t时的對应的平衡态分布函数。
格子速度：
[eα=cosα-1π2，sinα-1π2]c，
[α=1，2，3，4，] [c=dx/dt]
平衡态分布函数的计算：
[geqα=0.25C+0.5Ceα?u/c2]…（11）
式中[C]为浓度。
扩散系数与松弛时间关系：
[D=τ-0.5dx22dt]…（12）
2.2 边界条件处理
边界条件的一般表达式[16]：
[aαCαn+bC=d]…（13）
式中[a，b]系数，[C]浓度，[n]为边界的法方向。据[a，b]的不同取值，边界条件可分为3类：
①当[a]=0，[b]≠0时，为第一类边界条件或狄里克莱（Dirichlet）条件，即边界处浓度为定值;
②当[a]≠0，[b]=0时，为第二类边界条件或诺依曼（Neumann）条件，即给出边界处的浓度梯度;
③当[a]≠0，[b]≠0时，为第三类边界条件或混合边界条件。
本文模拟过程中，需处理的边界有定浓度边界，即第①种当[a]=0，[b]≠0时;零通量边界即第②种当[a]≠0，[b]=0时;酸岩接触面的反应边界即第③种当[a]≠0，[b]≠0时。下面分别介绍上述遇到的边界处理方法。
迁移过程之后要对未知方向分布函数进行处理。边界条件的处理是根据边界条件求出边界处的浓度值，再由浓度值求出未知方向的分布函数。
2.2.1 定浓度边界
定浓度边界即边界处浓度为定值，即[C=C0]。假设边界处[gα]未知，Huber等根据公式浓度与分布函数关系[17]，在边界处可求未知方向的[gα]，即
[gα=C-i≠αngα]…（14）
2.2.2 零通量边界
当模拟区域边界为固体边界或自由出口时，满足零通量边界条件，即[αCαn=0]。
同流动过程一样，非平衡外推边界也可用于计算浓度边界。本次模拟采用非平衡外推边界方法计算浓度边界[18]。由[αCαn=0]可得到边界处浓度，如上边界处：
[Cx，y=L，t=Cx，y=L-Δx，t]…（15）
将分布函数分为平衡态与非平衡之和，即[g4=geq4+gneq4]，非平衡态部分[gneq4]由相邻节点处的非平衡态得到。
2.2.3 酸岩反应边界
假设酸岩反应为一级反应，酸岩接触面满足：
[D?C=krC-Ceq]…（16）
式中，[D]酸液的扩散系数，[kr]酸岩反应常数，[Ceq]平衡浓度。
本次模拟采用张婷提出的半反弹边界处理格式[20]，边界处：
[gi，xf，t+Δt=-g+ixf，t+2wiCw]…（17）
式中：[g+i]为碰撞之后的分布函数，[gi，]为[gi]的相反方向，[wi]为对应的权值，[Cw]为边界处的浓度。
3 酸岩非均质反应
3.1 非均匀反应的理论原因
不均衡因素，如渗透率大小不同、孔隙大小不同、裂缝的存在、矿物成分不均等，会造成酸的分布和反应不均衡。流动阻力小的地方促使流体流量增大，流量增大有利于易溶物质的溶解，而物质溶解又使得介质孔隙度渗透率不断增大，这样的反馈结果是不均一溶蚀。下面给出理论推导和数值验证。 流体在孔道中流动符合Poiseuille流动规律，流量公式为：
[Q=πr4?p8μl]…（18）
式中：[Q]为体积流量;[r]为孔道半径;[μ]为流体粘度;[?p/l]压力梯度。
假设在同一流动区域两半径分别为[r1]和[r2]的长度相等的孔道，由于并联的压差是相等的，同一种流体在两孔道中的流量之比为：
[Q1Q2=πr41?p8μlπr4?p8μl=r41r42]…（19）
一定体积的酸液代表一定的溶蚀能力。假设通过通道的酸液完全消耗，则通过通道的酸液体积与该通道中溶蚀的岩石体积成正比。两孔道中溶蚀的岩石体积[V1]，[V2]的关系为[V1/V2]=[r1/r24]，即通道半径之比的4次方。假设[r1/r2=2]，则[V1/V2=16]，即通道1溶蚀掉的岩石是通道2的16倍，所以通道1半径扩展远远大于通道2。
3.2 模型建立
模拟区域为2维矩形结构，左侧为入口端，右侧为出口端，上下侧为固体边界。假设酸液传质过程对流动没有影响，可将酸蚀蚓孔形成的整个过程分为酸液在多孔介质中的流动、由于浓度及速度引起的酸液的传质、酸液与岩石的反应、岩石颗粒的减小及孔隙空间的扩展等过程。用第1节介绍的模拟流动过程的D2Q9格子的MRT-LBM模型模拟注入酸液的速度场，模拟区域的入口用Zou-He边界处理，出口用非平衡外推方法处理，上下边界及内部流固边界用反弹边界方法处理。用第2节介绍的模拟传质过程的D2Q4格子的LBM模型模拟酸液的浓度场，模拟区域入口用定浓度边界，出口及上下边界用零通量边界处理，流固表面用酸岩反应边界处理。同一时间内，传质过程将采用流动过程计算出的速度场，从而完成两个过程的耦合。
特征长度0.001 mm，粘度为1 mpa·s，流体格子浓度为1.0，松弛时间为1.0，盐酸浓度为1 mol/l，盐酸密度为1 015 kg/m3，反应速率常数为0.002 m/s，扩散系数为3.6×10-9 m2/s，注入速度为1 mm/s。
3.3 模型验证
模拟区域大小为114×50，图4表示的是孔隙介質结构，红色部分为固体基质，蓝色部分为孔隙空间。
图4分别为酸岩反应前、后的孔隙结构分布。从图中可看出，入口处颗粒变小，孔隙空间增大，远处由于酸液浓度很低，颗粒基本保持原样。酸液浓度随距离入口端长度增加逐渐降低。模拟结果与实验及前人研究成果相符，所以该模型能准确模拟酸蚀过程。
3.4 非均匀反应过程模拟
建立大小为100×70的模拟区域，分别模拟等径与不等径情况下的反应结果。
当孔径大小为1∶1时，图5所示为通酸前后的孔道结构，红色部分为固体颗粒，蓝色部分为孔隙空间。模拟时间为[t=1.33 s]。
当孔径大小为2：1时，图6所示为通酸前后的孔道结构，红色部分为固体颗粒，蓝色部分为孔隙空间。模拟时间为[t=1.33 s]。
从上述模拟可看出，当两通道大小为1∶1时，两通道的速度、浓度分布及最终的孔道结构没差别;当两通道大小为2∶1时，大直径通道内的速度大约是小直径通道内速度的4倍，大直径通道中的浓度高，小通道中除入口处外浓度基本为0。孔隙空间扩展结果可看出，大直径通道扩展明显，小直径通道基本没扩展。因此认为，通道大小的差异导致了酸蚀结果差异。
4 结论
（1） LBM方法能实现酸化过程中流动及传质过程的模拟及耦合，从而实现对酸化过程的微观尺度模拟。
（2） 孔隙尺寸的不均匀分布是造成碳酸盐岩非均匀溶蚀的重要原因。
参考文献
[1] G.Daccord.Chemical dissolution of a porous media by a reactive fluid[J].Physical Review Letters，1987，58：479-482.
[2] M.L.Hoefner，H.S.Fogler.Porous evolution and channel formation during flow and reaction in porous media[J].American Institute of Chemical Engineers Journal，1988，34（1）：45-54.
[3] C.N.Fredd，H.S.Fogler.Optimum conditions for wormhole formation in carbonate porous media：influence of transport and reaction[J].SPE，1999，4（3）：196-205.
[4] G.Daccord，R.Lenormand.Fractal patterns from chemical dissolution[J]. Natrual， 1987， 325（1）：41-43.
[5] k.M.Hung，A.D. Hill，K.Sepehrnoorl.A mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing[J].Journal of Petroleum Technology，1989：59-66.
[6] 柳明，张士诚，牟建业. 酸蚀蚓孔的分形性和酸液类型对蚓孔的影响[J].石油勘探与开发，2012，39：591-596 [7] Daccord， G.， Lenormand， R.， and Lietard， O.， 1993a， “Chemical Dissolution ofa Porous Medium by a Reactive Fluid-1. Model for the ‘Wormholing’ Phenomenon，”Chem. Eng. Sci.， 48（1）， pp. 169-178.
[8] Bazin， B.， 2001， “From Matrix Acidizing to Acid Fracturing： A LaboratoryEvaluation of Acid/Rock Interactions，” SPE Prod. Facil.， 16（1）， pp. 22-29.
[9] Ortoleva， P. J.， Chadam， J.， Merino， E.， and Sen， A.， 1987， “Geochemical Self-Organization-Part II： The Reactive-Infiltration Instability，” Am. J. Sci.，287（10）， pp. 1008-1040.
[10] Wei， C.， and Ortoleva， P.， 1990， “Reaction Front Fingering in Carbonate-Cemented Sandstones，” Earth-Sci. Rev.， 29（1-4）， pp. 183-198.
[11] Pan C.Luo L S.Miller C T.An evaluation of Lattice Boltzmann scheme for porous medium flow simulation[J].Compute&Fluids，2006，35：898-909.
[12] 楊帆，施徐明，郭雪岩，等. 多弛豫时间格子玻尔兹曼方法的分块算法[J]. 排灌机械工程学报，2013，31（1）：56-60.
[13] A.A.Mohamad.Lattice Boltzmann Method：Fundamentals and Engineering Applications with Computer Codes[M].2011
[14] 何雅玲，王勇，李庆，等.格子Boltzmann方法的理论及应用[M].北京：科学出版社，2008.
[15] QisuZou，XiaoyiHe.On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[J].Phys.Fluids，1997，9（6）：1591-1598.
[16] 张婷.多孔介质内多组分非均相反应流的格子Boltzmann方法研究[D].湖北：华中科技大学，2012.
[17] Christian Huber，BabakShafei，AndreaParmigiani.A new pore-scale model for linear and non-linear heterogeneous dissolution and precipitation[J].Geochimica et Cosmochimica Acta，2014，124：109-130.
[18] Hou SL， Zou Q， Chen SY， et al. Simulation of Cavity Flow by the Lattice BoltzmannMethod[J].Journal of Computational Physics， 1995， 118（2）：329-347.
New Application of Lattice Boltzmann Method in Simulation of Oil & Gas Revervoir Acidizing
Gao Qingchun1，Wang Zhiming1，Zeng Quanshu1，WeiWei2，ZhaiYongfan2，Mi liping2
（1.State Key Laboratory of Petroleum Resources and Prospecting， China University of Petroleum - Beijing，Beijing，102249， China;2. Laojunmiao Production Plant，Yumen Oilfield Compamy，Jiuquan，Gansu，735019，China）
Abstract：The flow model and transport model in the acidizing process is established by lattice Boltzmann method， respectively. And then， the two models are coupled by the velocity field， thereby， the numerical simulation of acidizing process at microscopic scale is obtained.The simulation results suggest that： The lattice Boltzmann method is able to simulate the acidizing process at microscopic scale accurately， and be a new tool to study acidizing in carbonate formations; The heterogeneity in pore size is the key factor resulting in acid-rock non-uniform dissolution.
Key words：Acidizing;Lattice Boltzmann;Non-uniform dissolution;Flow;Transport
关键词：酸化;格子Boltzmann;非均匀溶蚀;流动;传质
酸化是解除碳酸盐岩储层近井筒污染，提高油气井产量的重要措施[1-6]。大量学者通过实验手段及数模手段对宏观尺度上蚓孔扩展规律进行研究[7-10]，认为非均匀溶蚀是形成蚓孔的重要因素。本文利用格子Boltzmann方法，研究微观尺度上盐酸与碳酸盐岩非均匀溶蚀问题，证实格子Boltzmann方法对酸化过程模拟的可行性，为酸化数值模拟研究向微观方向发展提供新思路。
1 流动过程
1.1 流动过程的D2Q9模型（图1）
格子Boltzmann方程包含两部分：迁移和碰撞。式（1）为LB方程，表征流体流动过程，即分布函数的演化过程：
[fαx+eαΔt，t+Δt-fαx，t=-1τf αx，t-feqαx，t]…（1）
式中：[fα]为[α]方向的粒子速度分布函数;[Δt]为时间增量;[τ]为弛豫时间;[feqα]为相应的平衡态分布函数。
上式可分解为两个过程，即碰撞过程和迁移过程。
碰撞过程：
[fαx，t+Δt=fαx，t1-1τ+1τfeqαx，t]
迁移过程：
[fαX+Δx，t+Δt=fαx，t+Δt]
对于D2Q9模型，松弛时间与运动粘度关系为：
[v=Δx2c2sΔtτ-12]
为保证运动粘度为正值，弛豫时间的选取大于0.5，[c2s]为格子声速的平方。
平衡态分布函数计算式：
[feqαX=wαρX1+eα?uc2s+eα?u22c4s-u22c2s]…（2）
式中：[ρX]为[X]处宏观密度;[u]为宏观速度。
1.2 MRT模型
上述介绍的基于单松弛时间的LBGK模型（SRT）是最简单的格子玻尔兹曼模型，由于简单，LBGK模型为最流行的格子玻尔兹曼模型。
SRT模型的两个不足：一是数值稳定性差，二是会造成粘性相关的边界。一旦模型用于模拟多孔介质中流动时，可能导致计算得到的渗透率与粘性相关，这有悖于渗透率的物理意义。为克服BGK模型的上述不足，引入MRT模型模拟多孔介质内的流体流动[11]。
MRT模型在碰撞中使用多个松弛参数，在矩空间执行碰撞过程，在速度空间执行迁移过程。变换矩阵M[12]：
[ M=111111111-4-1-1-1-122224-2-2-2-21111010-101-1-110-20201-1-110010-111-1-100-20211-1-101-11-1000000000111-1]
[m]为速度矩矢量，[f]为分布函数，[meq]为平衡态速度矩矢量。关系：[m=Mf]。
[m=ρ，e，ε，jx，qx，jy，qy，pxx，pxyT]
[meq=ρ，-2ρ+3u2x+u2y，ρ-3u2x+u2y，]
[ux，-ux，uy，-uy，u2x-u2y，uxuyT]
各矩的物理意义，其中，[e]为能量;[ε]为能量平方;[jx，jy]分别为[x]和[y]方向的质量流量;[qx]，[qy]分别为[x]和[y]方向的能量流量;[pxx，pxy]分别为黏度应力张量对角和非对角线部分;[ux，uy]分别为[x]和[y]方向的宏观速度。
松弛矩阵[S]为对角矩阵，
[S=diags0，s1，s2，s3，s4，s5，s6，s7，s8];其中，[sv=1/τ]，用于确定剪切黏度系数。不同的矩可采用不同松弛时间，这是符合物理规律的。当[si=1/τi=0，1，2，L，8]时，MRT模型退回为LBGK模型。
多松弛模型中不同松弛参数组合对结果有较大影响，通过比较给出模拟多孔介质时最佳参数组合，即[s0，s3，s5]取0，[s1，s2，s7，s8]取[1/τ]，[s4，s6]取[82τ-18τ-1]，本文松弛矩阵的选取就采用此种组合。
1.3 边界条件处理
为将格子Boltzmann方法应用于实际问题，还要加上一定的边界条件，本文使用的边界条件有反弹边界、Zou-He边界、非平衡外推边界[13-15]。
1.3.1 反弹边界
反弹边界是LBM方法中最简单的边界处理法，处理复杂几何边界时，如多孔介质中流动时有着广泛应用。它将边界看成具反弹边界条件的边界，分布函数看成运动的刚性粒子，粒子于边界发生碰撞后，以原速度反方向返回，其值大小不变。表达式为：
[fα，Xf，t+Δt=fαxb，t]
[α，]表示[α]的反方向，[Xf]表示临近边界的流体格点，[xb]表示边界格点。
1.3.2 Zou-He速度压力边界
速度边界 实际应用中，最常见的是边界给定某一速度。图2是边界函数分布图。
以North边界为例，北部边界，未知分布函数为[f4，7，，8]，密度[ρ]，由宏观密度方程和动量方程可得：
[ρ=f0+f1+f3+2f2+f5+f61+v0]…（3）
[f4=f2-23ρv0]…（4） [f7=f5+12f1-f3-16ρv0]…（5）
[f8=f6+12f1-f3-16ρv0]…（6）
由此可得未知分布函数[f4，7，，8]及密度[ρ]，其他3个方向的边界条件求解过程同上。
1.3.3 非平衡外推边界
边界处分布函数分为平衡与非平衡两部分，其中非平衡部分与相邻节点处的非平衡部分近似，可用于处理速度与压力边界，基于此种格式处理的速度和压力边界有二阶精度。
假设边界格子节点为[O]，与之相邻的内部节点为[G]，[O]点处速度[u]已知，密度[ρ]未知，[fi]为未知方向分布函数。节点[O]处粒子分布函数由平衡态和非平衡态两部分组成，即：
[fiO，t=feqiO，t+fneqiO，t]…（7）
确定节点[O]处的平衡态粒子分布函数所需的密度可用格子节点[G]处的流体密度代替，即：
[feqiO，t=feqiρG，t，t]…（8）
对非平衡态部分[fneqiO，t]，则用邻近流体格子节点[G]处的非平衡态分布函数近似得到，即：
[fneqiO，t=fneqiG，t=fiG，t-feqiG，t]…（9）
该处理方法的的特点是边界格子节点上所有粒子分布函数值全重新赋值，数值稳定性较好（图3）。
2 传质过程
2.1 传质过程的D2Q4 模型
与流动过程类似，LBM也被用于模拟传质过程。由于控制方程中只有质量守恒方程，因此模拟中需要的速度方向更少。与D2Q9模型相比，在保证精度和稳定性前提下，D2Q4模型具计算速度快、消耗内存少等优点。本文采用D2Q4模型对传质过程进行模拟。
碰撞迁移过程：
[gαx+dx，t+dt-gαx，t=-1τgαx，t-geqαx，t]…（10）
式中：[t]为时间，[dt]时间增量，[dx]位置增量;[eα]为[α]方向发格子速度，矢量;[gα]为[α]方向的分布函数;[geqαx，t]为位置x处、时间t时的對应的平衡态分布函数。
格子速度：
[eα=cosα-1π2，sinα-1π2]c，
[α=1，2，3，4，] [c=dx/dt]
平衡态分布函数的计算：
[geqα=0.25C+0.5Ceα?u/c2]…（11）
式中[C]为浓度。
扩散系数与松弛时间关系：
[D=τ-0.5dx22dt]…（12）
2.2 边界条件处理
边界条件的一般表达式[16]：
[aαCαn+bC=d]…（13）
式中[a，b]系数，[C]浓度，[n]为边界的法方向。据[a，b]的不同取值，边界条件可分为3类：
①当[a]=0，[b]≠0时，为第一类边界条件或狄里克莱（Dirichlet）条件，即边界处浓度为定值;
②当[a]≠0，[b]=0时，为第二类边界条件或诺依曼（Neumann）条件，即给出边界处的浓度梯度;
③当[a]≠0，[b]≠0时，为第三类边界条件或混合边界条件。
本文模拟过程中，需处理的边界有定浓度边界，即第①种当[a]=0，[b]≠0时;零通量边界即第②种当[a]≠0，[b]=0时;酸岩接触面的反应边界即第③种当[a]≠0，[b]≠0时。下面分别介绍上述遇到的边界处理方法。
迁移过程之后要对未知方向分布函数进行处理。边界条件的处理是根据边界条件求出边界处的浓度值，再由浓度值求出未知方向的分布函数。
2.2.1 定浓度边界
定浓度边界即边界处浓度为定值，即[C=C0]。假设边界处[gα]未知，Huber等根据公式浓度与分布函数关系[17]，在边界处可求未知方向的[gα]，即
[gα=C-i≠αngα]…（14）
2.2.2 零通量边界
当模拟区域边界为固体边界或自由出口时，满足零通量边界条件，即[αCαn=0]。
同流动过程一样，非平衡外推边界也可用于计算浓度边界。本次模拟采用非平衡外推边界方法计算浓度边界[18]。由[αCαn=0]可得到边界处浓度，如上边界处：
[Cx，y=L，t=Cx，y=L-Δx，t]…（15）
将分布函数分为平衡态与非平衡之和，即[g4=geq4+gneq4]，非平衡态部分[gneq4]由相邻节点处的非平衡态得到。
2.2.3 酸岩反应边界
假设酸岩反应为一级反应，酸岩接触面满足：
[D?C=krC-Ceq]…（16）
式中，[D]酸液的扩散系数，[kr]酸岩反应常数，[Ceq]平衡浓度。
本次模拟采用张婷提出的半反弹边界处理格式[20]，边界处：
[gi，xf，t+Δt=-g+ixf，t+2wiCw]…（17）
式中：[g+i]为碰撞之后的分布函数，[gi，]为[gi]的相反方向，[wi]为对应的权值，[Cw]为边界处的浓度。
3 酸岩非均质反应
3.1 非均匀反应的理论原因
不均衡因素，如渗透率大小不同、孔隙大小不同、裂缝的存在、矿物成分不均等，会造成酸的分布和反应不均衡。流动阻力小的地方促使流体流量增大，流量增大有利于易溶物质的溶解，而物质溶解又使得介质孔隙度渗透率不断增大，这样的反馈结果是不均一溶蚀。下面给出理论推导和数值验证。 流体在孔道中流动符合Poiseuille流动规律，流量公式为：
[Q=πr4?p8μl]…（18）
式中：[Q]为体积流量;[r]为孔道半径;[μ]为流体粘度;[?p/l]压力梯度。
假设在同一流动区域两半径分别为[r1]和[r2]的长度相等的孔道，由于并联的压差是相等的，同一种流体在两孔道中的流量之比为：
[Q1Q2=πr41?p8μlπr4?p8μl=r41r42]…（19）
一定体积的酸液代表一定的溶蚀能力。假设通过通道的酸液完全消耗，则通过通道的酸液体积与该通道中溶蚀的岩石体积成正比。两孔道中溶蚀的岩石体积[V1]，[V2]的关系为[V1/V2]=[r1/r24]，即通道半径之比的4次方。假设[r1/r2=2]，则[V1/V2=16]，即通道1溶蚀掉的岩石是通道2的16倍，所以通道1半径扩展远远大于通道2。
3.2 模型建立
模拟区域为2维矩形结构，左侧为入口端，右侧为出口端，上下侧为固体边界。假设酸液传质过程对流动没有影响，可将酸蚀蚓孔形成的整个过程分为酸液在多孔介质中的流动、由于浓度及速度引起的酸液的传质、酸液与岩石的反应、岩石颗粒的减小及孔隙空间的扩展等过程。用第1节介绍的模拟流动过程的D2Q9格子的MRT-LBM模型模拟注入酸液的速度场，模拟区域的入口用Zou-He边界处理，出口用非平衡外推方法处理，上下边界及内部流固边界用反弹边界方法处理。用第2节介绍的模拟传质过程的D2Q4格子的LBM模型模拟酸液的浓度场，模拟区域入口用定浓度边界，出口及上下边界用零通量边界处理，流固表面用酸岩反应边界处理。同一时间内，传质过程将采用流动过程计算出的速度场，从而完成两个过程的耦合。
特征长度0.001 mm，粘度为1 mpa·s，流体格子浓度为1.0，松弛时间为1.0，盐酸浓度为1 mol/l，盐酸密度为1 015 kg/m3，反应速率常数为0.002 m/s，扩散系数为3.6×10-9 m2/s，注入速度为1 mm/s。
3.3 模型验证
模拟区域大小为114×50，图4表示的是孔隙介質结构，红色部分为固体基质，蓝色部分为孔隙空间。
图4分别为酸岩反应前、后的孔隙结构分布。从图中可看出，入口处颗粒变小，孔隙空间增大，远处由于酸液浓度很低，颗粒基本保持原样。酸液浓度随距离入口端长度增加逐渐降低。模拟结果与实验及前人研究成果相符，所以该模型能准确模拟酸蚀过程。
3.4 非均匀反应过程模拟
建立大小为100×70的模拟区域，分别模拟等径与不等径情况下的反应结果。
当孔径大小为1∶1时，图5所示为通酸前后的孔道结构，红色部分为固体颗粒，蓝色部分为孔隙空间。模拟时间为[t=1.33 s]。
当孔径大小为2：1时，图6所示为通酸前后的孔道结构，红色部分为固体颗粒，蓝色部分为孔隙空间。模拟时间为[t=1.33 s]。
从上述模拟可看出，当两通道大小为1∶1时，两通道的速度、浓度分布及最终的孔道结构没差别;当两通道大小为2∶1时，大直径通道内的速度大约是小直径通道内速度的4倍，大直径通道中的浓度高，小通道中除入口处外浓度基本为0。孔隙空间扩展结果可看出，大直径通道扩展明显，小直径通道基本没扩展。因此认为，通道大小的差异导致了酸蚀结果差异。
4 结论
（1） LBM方法能实现酸化过程中流动及传质过程的模拟及耦合，从而实现对酸化过程的微观尺度模拟。
（2） 孔隙尺寸的不均匀分布是造成碳酸盐岩非均匀溶蚀的重要原因。
参考文献
[1] G.Daccord.Chemical dissolution of a porous media by a reactive fluid[J].Physical Review Letters，1987，58：479-482.
[2] M.L.Hoefner，H.S.Fogler.Porous evolution and channel formation during flow and reaction in porous media[J].American Institute of Chemical Engineers Journal，1988，34（1）：45-54.
[3] C.N.Fredd，H.S.Fogler.Optimum conditions for wormhole formation in carbonate porous media：influence of transport and reaction[J].SPE，1999，4（3）：196-205.
[4] G.Daccord，R.Lenormand.Fractal patterns from chemical dissolution[J]. Natrual， 1987， 325（1）：41-43.
[5] k.M.Hung，A.D. Hill，K.Sepehrnoorl.A mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing[J].Journal of Petroleum Technology，1989：59-66.
[6] 柳明，张士诚，牟建业. 酸蚀蚓孔的分形性和酸液类型对蚓孔的影响[J].石油勘探与开发，2012，39：591-596 [7] Daccord， G.， Lenormand， R.， and Lietard， O.， 1993a， “Chemical Dissolution ofa Porous Medium by a Reactive Fluid-1. Model for the ‘Wormholing’ Phenomenon，”Chem. Eng. Sci.， 48（1）， pp. 169-178.
[8] Bazin， B.， 2001， “From Matrix Acidizing to Acid Fracturing： A LaboratoryEvaluation of Acid/Rock Interactions，” SPE Prod. Facil.， 16（1）， pp. 22-29.
[9] Ortoleva， P. J.， Chadam， J.， Merino， E.， and Sen， A.， 1987， “Geochemical Self-Organization-Part II： The Reactive-Infiltration Instability，” Am. J. Sci.，287（10）， pp. 1008-1040.
[10] Wei， C.， and Ortoleva， P.， 1990， “Reaction Front Fingering in Carbonate-Cemented Sandstones，” Earth-Sci. Rev.， 29（1-4）， pp. 183-198.
[11] Pan C.Luo L S.Miller C T.An evaluation of Lattice Boltzmann scheme for porous medium flow simulation[J].Compute&Fluids，2006，35：898-909.
[12] 楊帆，施徐明，郭雪岩，等. 多弛豫时间格子玻尔兹曼方法的分块算法[J]. 排灌机械工程学报，2013，31（1）：56-60.
[13] A.A.Mohamad.Lattice Boltzmann Method：Fundamentals and Engineering Applications with Computer Codes[M].2011
[14] 何雅玲，王勇，李庆，等.格子Boltzmann方法的理论及应用[M].北京：科学出版社，2008.
[15] QisuZou，XiaoyiHe.On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[J].Phys.Fluids，1997，9（6）：1591-1598.
[16] 张婷.多孔介质内多组分非均相反应流的格子Boltzmann方法研究[D].湖北：华中科技大学，2012.
[17] Christian Huber，BabakShafei，AndreaParmigiani.A new pore-scale model for linear and non-linear heterogeneous dissolution and precipitation[J].Geochimica et Cosmochimica Acta，2014，124：109-130.
[18] Hou SL， Zou Q， Chen SY， et al. Simulation of Cavity Flow by the Lattice BoltzmannMethod[J].Journal of Computational Physics， 1995， 118（2）：329-347.
New Application of Lattice Boltzmann Method in Simulation of Oil & Gas Revervoir Acidizing
Gao Qingchun1，Wang Zhiming1，Zeng Quanshu1，WeiWei2，ZhaiYongfan2，Mi liping2
（1.State Key Laboratory of Petroleum Resources and Prospecting， China University of Petroleum - Beijing，Beijing，102249， China;2. Laojunmiao Production Plant，Yumen Oilfield Compamy，Jiuquan，Gansu，735019，China）
Abstract：The flow model and transport model in the acidizing process is established by lattice Boltzmann method， respectively. And then， the two models are coupled by the velocity field， thereby， the numerical simulation of acidizing process at microscopic scale is obtained.The simulation results suggest that： The lattice Boltzmann method is able to simulate the acidizing process at microscopic scale accurately， and be a new tool to study acidizing in carbonate formations; The heterogeneity in pore size is the key factor resulting in acid-rock non-uniform dissolution.
Key words：Acidizing;Lattice Boltzmann;Non-uniform dissolution;Flow;Transport