2. 中国石油塔里木油田分公司勘探开发研究院, 新疆 库尔勒 841000;
3. 中国石油长庆油田分公司勘探开发研究院, 西安 710018;
4. 长江大学石油工程学院, 武汉 430100
2. Research Institute of Exploration and Development, Tarim Oilfield Company, PetroChina, Korla 841000, Xinjiang, China;
3. Exploration & Development Research Institute, Changqing Oilfield Company, Petrochina, Xi'an 710018, China;
4. College of Petroleum Engineering, Yangtze University, Wuhan 430100, China
0 前言
岩石物性应力敏感指的是岩石孔隙度、渗透率等随有效应力的变化[1]，以渗透率应力敏感研究最为热门，研究重点以实验测试与数据处理方法为主[2]。以往学者考虑典型孔隙形状，如圆形[3]、椭圆形[4]、圆锥形[4]、双曲边三边形[5]、双曲边四边形 (即星形)[6]，建立了孔隙渗透率应力敏感理论模型；但这些模型仅考虑了单一形状孔隙模型，不足以从机理上解释渗透率应力敏感。
岩石孔隙结构决定了其物理性质及其相互之间的关系，更控制了岩石物性随有效应力的变化规律，通过建立岩石孔隙网络模型可以构建其孔隙结构。孔隙网络模型有基于孔隙空间成像的网络模型[7]和规则孔隙网络模型[8]。基于孔隙空间成像的网络模型实验成本昂贵[7, 9]，而规则孔隙网络通过与逾渗理论相结合，采用一定的分布函数设定孔隙喉道大小及分布，即可建立与真实岩石同样复杂的孔隙网络模型[10]。本次研究基于QT Creator平台、采用C++编程生成随机孔隙网络模型，通过各种孔隙形状不同比例的组合，从微观孔隙网络模型角度着手，研究无因次半径随应力的变化；并通过建立渗透率与无因次半径的函数关系以实现渗透率应力敏感机理解释。
1 理论基础 1.1 逾渗理论逾渗理论主要用于描述随机的结构与连通性，其实质是用统计方法来研究流体在随机无序介质，如随机孔隙网络模型中的分布和流动规律。本次研究的随机孔隙网络模型有如下假设：模型中连线代表具有一定体积和流动阻力的喉道，它具有一定的分布特征，如均匀分布；而连线的交点则视为没有体积和流动阻力，只起连接的作用。所以，逾渗参数主要针对连线的特性进行计算，而不考虑节点。逾渗理论服从概率论，统计的样品必须达到一定要求才可靠，所以基于逾渗理论的网络模拟要求三维网络模型的节点数至少应达到20×20×20[8, 10]，且模型越大，越能反应宏观属性的一般性。
1.2 水电相似原理网络模拟的基础为水电相似原理。流体在多孔介质中流动，可以采用类似于电路中的网络结构进行模拟分析。在研究网状电路时，假定电路中所有电阻器由充满传导流体的圆形截面管束组成，由Ohm定律得
(1) 式中：I为电流；ΔE为电压；r为管束半径；ρ为管束电阻率；l为电阻器或管束长度。
假设流体流动为层流，以最简单的圆柱形管束为例，其Poiseuille方程为
(2) 式中：q为流量；μ为流体黏度；Δp为管束两端的压差。
电流电路网络与流体孔隙网络均由管束构成，对比式 (1) 与式 (2)：Ohm定律和Poiseuille方程都表述了电流或流体的流动与管束两端的压差及管束阻力的关系。这种流动关系表现为水电相似性，因此流体流动的Poiseuille方程可以简化为类似于电流中的Ohm定律的形式：
(3) 其中，
由上述水电相似原理可知，流体在多孔介质网络中的流动与电流在电路网络中的流动存在相似性；因此在研究孔隙网络模型时，可以将电路网络中的电流动借鉴到流体流动模拟中，并且可以将电流的分析方法用于分析、模拟流体流动。
1.3 Kirchoff定律电流电路网络中的Kirchoff定律有以下两种形式：①Kirchoff电流定律，对于任意一时刻，流入电路中任一节点的电流与流出的电流在数值上相等；②Kirchoff电压定律，对于任意一个电路回路，按照电流流动的方向绕行一周，各电路之路中的电压在任意一个瞬间恒为0。对于孔隙网络模型，一般采用Kirchoff电流定律建立模型方程组。
2 随机孔隙网络模型 2.1 部分饱和岩石有效应力方程推导本次研究的是准静态两相流随机孔隙网络模型，模型中流体流动完全由毛管压力控制，并且假设流体不可压缩，忽略黏滞力的影响。这种假设在物理模型中表现为，如果不发生驱替行为，模型中的油水或气水界面保持相对静止，与大多数多孔介质的低速渗流情况一致。在网络模拟中，只要所加的非湿相流体压力大于或等于某一喉道的入口毛管压力，驱替行为就开始，并且非湿相很快将湿相流体驱替出这一喉道，未被驱替到的喉道仍充满湿相流体，并保持压力值 (0) 不变。在这种情况下，被驱替喉道中的毛管压力为
(4) 式中：pc为毛管压力；pn为非湿相流体驱替压力。
被驱替的喉道中只存在一种流动流体，处于角隅处的湿相流体被认为是静止的、不可流动的，其压力值保持为0；管束中的流动可以看成是单相流体流动，孔隙流体压力即为非湿相驱替压力。由于未被驱替喉道中的流体为湿相水，压力为0，没有非湿相与湿相的压力差，不存在毛管压力，故本文不对其进行研究。因此，本次研究中有效应力表达式为
(5) 式中：peff为有效应力；pov为围压；pf为孔隙流体压力；α为有效应力系数。
这一理论可通过文献[11]得以验证。对于部分饱和多孔介质，Bishop[11]通过实验定义了广义的有效应力：
(6) 式中：Sw表示湿相饱和度；pw表示湿相压力。
当Sw为两个端点极值 (0和1) 时，Sw(pn-pw) 为0。当其取一个较低值时，比如0.1或者更小，Sw(pn-pw) 出现峰值；然而，它仍然是很小的，以至于与有效应力比起来在绝大多数情况下都可以忽略[11]。因此，式 (6) 化简为式 (5)。
式 (5) 为两相流体的有效应力方程，其系数α一直是研究的热点与难点[12-13]：如果α=0，有效应力就指围压；如果α=1，此时的有效应力为净应力[12-13]。为了简化研究，这里取α=0。在这种情况下，认为孔隙流体压力对有效应力不产生影响，或者说是孔隙流体压力相比于围压，其对有效应力的影响作用微弱。这样，本文的岩石微观孔隙结构随有效应力变化研究就停留在不同围压下的岩石孔隙结构变化上，这与实验方法也能很好地吻合。
2.2 程序核心方程整理本程序的核心方程包括流体在喉道中的水力传导率方程和无因次半径方程，各种形状喉道的核心方程如表 1所示。
构建运算表达式是网络模型模拟的核心部分，其建立的基本思想是Kirchoff定律、节点以及连线的传导性。Noble[14]给出了n个节点网络的方程：
(7) 其中：q=[q1, q2，..., qn]′；V=[V1, V2，..., Vn]′。式中：q为流量矩阵；qn为节点n的流量；A为关联矩阵；K为对角矩阵；V为电压矩阵；Vn为节点n的电压。
由统计原理，网络模型中的节点、连线数需要达到足够大才能保证模拟的可靠性，运用Cholesky分解系数矩阵，精度会受到舍入误差的影响，并且程序编写代码冗赘[8]。为此，采用迭代法求解矩阵。下面以最简单的二维正方形网络模型为例，对其求解过程进行介绍。
1) 给节点赋初始值。由于最终解与给定的初始值无关，因此初始值可以任意设定，但迭代中收敛的速度取决于初始值估计的准确程度。假定在正方形网络中，流动方向从左到右，这样最左边边界节点赋值电压V1，最右边边界节点赋值电压V2，其余节点赋值0。
2) 建立方程。以图 1中的0节点为例，根据Kirchoff电流定律有：
(8)
|
| 图 1 正方形网络模型的一个节点单元 Figure 1 One node cell of square network model |
|
|
式中：Vi, j为节点 (i, j) 的电压；Gi, j为节点 (i, j) 的电导。
对于网络中的每一个节点，都可以写出类似方程，这样就得到了方程数与节点数相等的线性方程组，变形为
(9) 3) 迭代求解。进行每一次迭代求解后，计算出网络中输入/输出电流大小。迭代结束的要求是输入/输出电流相等或其差值满足设定的误差精度。当迭代结束后，求出输入/输出电流，由Ohm定律计算出电阻大小，进而可以求取网络模型中的其他参数。
2.4 程序实现根据网络模型理论基础及核心编程点，基于QT平台采用C++语言编写随机网络模型程序。网格大小 (x、y、z方向大小)、最大/最小喉道半径、最大/最小纵横比、地层水黏度、网格单位长度 (管束长度)、连通概率，以及各形状孔隙所占比例、纵横比及其比例由程序界面输入，如图 2所示。
|
| 图 2 程序输入参数界面 Figure 2 Input parameters interface |
|
|
程序生成的三维规则立方体随机孔隙网络模型如图 3所示，模型节点数为80×80×80。
|
| 图 3 三维孔隙网络模型 Figure 3 3-D pore network model |
|
|
为了对比研究纵横比对计算结果的影响，程序设定为单一形状喉道，通过设置不同纵横比，进而研究不同纵横比条件下无因次半径随有效应力的变化情况，结果如图 4所示。
|
| a.椭圆形喉道；b.锥形喉道；c.双曲边三边形喉道；d.星形喉道。 图 4 不同纵横比条件下无因次半径随有效应力的变化 Figure 4 Dimensionless radius response to the effective stress under different aspect ratio |
|
|
由图 4可知，纵横比越小，无因次半径受到有效应力影响越明显。以星形喉道为例进行解释：由星形喉道无因次半径方程知，纵横比越小，无因次半径随应力变化效果越显著 (表 1)；星形喉道在受到应力变化时发生的形变也可以说明，纵横比越小，星形喉道越接近裂缝形状 (图 5)，显然裂缝受到应力作用反映明显。
|
| a. ε=0.200；b. ε=0.005。 图 5 不同纵横比的星形喉道形状 Figure 5 Star-shaped throat at different aspect rati |
|
|
下面讨论不同纵横比按不同比例组合时无因次半径与有效应力的关系。以星形喉道为例，考虑单一形状喉道，纵横比取值1.000、0.100、0.050、0.001，所占比例分别为20%、20%、30%、30%的组合方案，运行结果如图 6所示。图 6表明，不同纵横比按不同比例组合时，无因次半径随有效应力的变化更加明显。
|
| 组合方案：纵横比取值1.000、0.100、0.050、0.001，所占比例分别为20%、20%、30%、30%。 图 6 星形喉道无因次半径与有效应力的关系 Figure 6 Dimensionless radius response to the effective Stress by star-shaped throat |
|
|
由前面分析可知，在有效应力的作用下，岩石的孔隙结构发生变化，进而引起其物性变化；因此，渗透率应力敏感机理研究的本质为孔隙半径对有效应力的响应，即
(11) 式中：r为孔隙半径。
由表 1可知，无因次半径是有效应力的函数。而模型程序中孔隙半径与无因次半径存在如下关系：
(12) 式中：r0为初始孔隙半径。
这样便建立了孔隙半径与有效应力的函数。水力传导率是孔隙半径的函数 (表 1)，程序中输入/输出流量又是水力传导率的函数，而基于达西定律模拟计算的渗透率是输入/输出流量的函数，从而建立起了渗透率与有效应力的函数。
组合40%椭圆形喉道和60%星形喉道：椭圆形喉道纵横比取值0.010 0、0.050 0、0.001 2，所占比例分别为60%、30%、10%；星形喉道纵横比取值0.100 0、0.050 0、0.001 2，所占比例分别为20%、30%、50%。将程序运行结果与岩心无因次渗透率与有效应力实验数据[2]进行对比，拟合度高 (图 7)，因此从网络模型机理解释渗透率应力敏感。
|
| 组合方案：40%椭圆形喉道、60%星形喉道：椭圆形喉道纵横比取值0.010 0、0.050 0、0.001 2，所占比例分别为60%、30%、10%；星形喉道纵横比取值0.100 0、0.050 0、0.001 2，所占比例分别为20%、30%、50%。 图 7 渗透率应力敏感机理解释 Figure 7 Mechanism interpretation of permeability stress sensitivity |
|
|
1) 推导了部分饱和的岩石有效应力方程，并应用迭代法求取网络模型的核心运算表达式。
2) 纵横比越小，无因次半径受到应力的影响越明显。
3) 基于网络模型程序可以联立渗透率与无因次半径的关系，进而建立渗透率与有效应力的关系，并通过多种形状孔隙组合有效地解释了渗透率应力敏感，可以推广为基于本文网路模拟的岩石物性应力敏感机理解释。
| [1] | Tiab D, Donaldson E C. Petrophysics[M]. New York: Elsevier, 2012. |
| [2] | Li M, Bernabe Y, Xiao W L, et al. Effective Pressure Law for Permeability of E-Bei Sandstones[J]. Journal of Geophysical Research:Solid Earth, 2009, 114(B7): 223-223. |
| [3] | Jaeger J C, Cook N G W, Zimmerman R W. Funda-mentals of Rock Mechanics[M]. London: Blackwell Publishing, 2007. |
| [4] | Seeburger D A, Nur A. A Pore Space Model for Rock Permeability and Bulk Modulus[J]. Journal of Geophysical Research, 1984, 89(B1): 527-536. DOI:10.1029/JB089iB01p00527 |
| [5] | Yale D P.Network Modeling of Flow, Storage, and De-formation in Porous Rock[D].Stanford Calif:Stanford University, 1984. |
| [6] | Sigal R F. The Pressure Dependence of Permeability[J]. Petrophysics, 2002, 43(2): 92-102. |
| [7] | Fredrich J T. 3D Imaging of Porous Media Using Laser Scanning Confocal Microscopy with Application to Microscale Transport Processes[J]. International Journal of Rock Mechanics Sciences & Geomechanics Abstracts, 1993, 24(99): 551-561. |
| [8] | 陶正武. 毛管压力曲线数值实验研究[D]. 成都: 西南石油大学, 2013. Tao Zhengwu.Numerical Experiments of Capillary Pressure Curve[D].Chengdu:Southwest Petroleum University, 2013. |
| [9] | 张婷, 徐守余, 王子敏. 储层微观孔喉网络图形识别方法[J]. 吉林大学学报 (地球科学版), 2011, 41(5): 1646-1650. Zhang Ting, Xu Shouyu, Wang Zimin. Identificaiton Method of Reservoir Microscopic Pore-Thort Network[J]. Journal of Jilin University (English Science Edition), 2011, 41(5): 1646-1650. |
| [10] | Tsakiroglou C D, Fleury M. Pore Network Analysis of Resistivity Index for Water-Wet Porous Media[J]. Transport in Porous Media, 1999, 35(35): 89-128. |
| [11] | Bishop A W. The Principle of Effective Strength[J]. Teknisk Ukeblad, 1959, 106(39): 859-863. |
| [12] | Gray W G, Schrefler B A. Thermodynamic Approach to Effective Stress in Partially Saturated Porous Media[J]. European Journal of Mechanics:A:Solids, 2001, 20(4): 521-538. DOI:10.1016/S0997-7538(01)01158-5 |
| [13] | Nuth M, Laloui L. Effective Stress Concept in Unsa-turated Soils:Clarification and Validation of a Unified Framework[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(7): 771-801. DOI:10.1002/(ISSN)1096-9853 |
| [14] | Noble B. Applied Linear Algebra[M]. New Jersey: Englewood Cliffs, 1969. |