砂岩是一种主要由砂粒胶结而成的沉积岩，大量离散的颗粒相互作用形成复杂体系，具有非连续和强接触耗散等结构特性，经典的连续介质方法无法描述其颗粒体系。为了更准确地考虑颗粒结构(形状、大小、排列等)以及赋存环境等对砂岩力学特性的影响，基于离散元理论的颗粒流方法(particle flow code, PFC)应运而生，该方法克服了传统连续介质模型的宏观连续性假设，从细观角度研究砂岩的颗粒运动和力学响应。目前，线性接触模型和Hertz-Mindlin接触模型的运用最为广泛，它们是通过设置法向和切向弹簧来模拟颗粒接触力与位移之间的弹性关系，如吴良平^[1]、蒋明镜^[2]、Cui L^[3]、徐辉^[4]、贾宇峰^[5]、常在等^[6]基于上述接触模型研究了普通砂土颗粒的运动特性，并建立了相应的本构关系。

对于含有胶结物质的砂岩，胶结物限制了砂岩颗粒的转动，并对承力起着重要的作用，上述接触模型无法描述胶结物质的粘结特性。那么，要准确获得胶结砂岩的力学特性就需要运用新的数值处理技术，如Cundall^[7]、Potyondya^[7-8]、Jung-Wook Park^[9]以及Yo-Ming Hsieh等^[10-11]提出了抗转动的接触力学模型来模拟岩体的力学特性，这些新的颗粒接触模型从一定程度反映砂岩颗粒之间的胶结性。随着颗粒流数值技术的发展，模拟颗粒粘结的2种接触模型逐渐运用于实际工程中，接触粘结模型能反映颗粒之间的张拉特性，但粘结只出现于点上而不能限制颗粒的转动；具有一定接触面积的平行粘结模型不仅能限制砂岩颗粒的转动，还能传递力和力矩，能更好地模拟胶结性状对砂岩力学特性的影响。

文中以油藏储层砂岩为研究对象，研究石油开采过程中地层砂脱落的岩体结构破坏机理，其中砂岩的胶结性状对地层结构的受力特性起着重要的控制作用^[12-13]。为了更准确地描述胶结砂岩的力学特性，运用随机数学理论进行PFC3D的二次开发获得颗粒胶结物质的分布，基于平行粘结模型建立数值模型模拟胶结砂岩的宏细观力学响应和岩体结构破坏机理，详细研究砂岩的胶结性状如胶结程度(采用胶结半径比α=R/R_A表示，R为胶结平均半径，R_A为接触处颗粒中的较小半径，如图 1和图 2所示)、平行粘结刚度和胶结量(采用颗粒胶结数/模型总接触数表示)对其力学特性的影响。其研究成果能为砂岩体细观结构和本构模型的建立提供重要的理论依据，特别是对砂岩储层的出砂机理、出砂预测以及砂控策略研究具有重要的实践指导意义。

1 胶结砂岩力学响应的PFC3D模型 1.1 PFC3D中的平行粘结模型

平行粘结模型如图 1和图 2所示，该模型能描述胶着颗粒之间的本构关系，具体表现为一系列弹簧(法向、切向弹簧)来模拟颗粒之间的力学状态。不同于常规的接触模型，该粘结模型在接触处表现为具有一定面积和刚度的粘结，不仅能限制颗粒旋转，还能传递力和力矩。受力过程中，颗粒间的胶结性状会发生改变，有效刚度和接触力也会随之变化，虽然该模型不能完全模拟实际砂岩体的胶结量，但力学分析得到的胶结性状改变也从一定程度上反映了胶结砂岩的力学响应趋势。通过式(1)和式(2)获得应力值

$ {\sigma _{{\rm{max}}}} = \frac{{ - {{\bar F}^n}}}{A} + \frac{{\left| {{{\bar M}^s}_i} \right|}}{I}\bar R, $

(1)

$ {\tau _{{\rm{max}}}} = \frac{{\left| {{{\bar F}^s}_i} \right|}}{A} + \frac{{\left| {{{\bar M}^n}} \right|}}{J}\bar R, $

(2)

其中：A=πR²；$J = \frac{1}{2}{\rm{ \mathit{ π} }}{\bar R^4}$；$I = \frac{1}{4}{\rm{ \mathit{ π} }}{\bar R^4}$；R为接触处平行粘结接触的平均半径。

当平行粘结接触处计算的张拉应力或剪应力大于张拉强度和剪切强度时，平行粘结就宣告破坏，该粘结一旦破坏就不可恢复，这就意味着后续的数值分析中，若颗粒间还有接触，也只能运用常规的接触模型来描述砂岩体的力学特性。因此，计算中需要选择平行粘结力学参数来模拟颗粒间的胶结性状：法向和切向刚度，法向和切向强度，胶着程度等；另外，选择常规接触的力学参数进行无胶结颗粒的接触分析，如平行粘结破坏后的颗粒均需采用常规接触进行力学分析。

1.2 数值模型验证 1.2.1 数值模型设计

建立如图 3所示的颗粒流数值模型，模型高度为0.08 m，直径为0.04 m，相关参数如表 1所示。然后对该模型进行数值分析。

1) 对初始的砂岩颗粒流模型逐渐加载，一旦荷载达到10 MPa，以导入平行粘结模拟砂岩的胶结性状。

2) 基于随机数学理论设置颗粒胶结为均匀分布。采用PFC3D的Fish语言产生随机数，而每个随机数对应一个颗粒接触，对所有颗粒接触进行循环，如果产生的随机数大于实际要求的颗粒胶结量百分比，设该处颗粒接触为无胶结的接触状态，其余则为平行粘结胶结状态，然后通过PFC3D的二次开发将平行粘结的力学参数导入模拟具有一定面积、刚度和胶结量的砂岩体。

3) 开始模型剪切试验，剪切开始前的模型孔隙率为0.593。剪切试验中，设置模型顶板加载速度为0.02 mm/s，底板静止，类似于实验室的等应变加载；实际加载值与数值计算荷载之差控制于±2.5%范围，以保证整个数值模拟过程的准静止状态。另外，数值分析中通过伺服系统调节圆柱体半径值来维持侧向应力为10 MPa。

1.2.2 数值计算结果分析

由数值模型，得到剪切过程中胶结砂岩的宏观力学特性(应力比、体应变)和颗粒的细观力学响应规律(配位数、粘结破坏数)。

1) 图 4和图 5分别描述了应力比和体应变随轴向应变(ε_a)的变化趋势。应力比随着ε_a的增大而增大，应变为3.5%时，应力比达到最大峰值为0.76；峰值之后的应变软化很明显，其后软化率随着应变的增大而减小，符合岩体的应力-应变变化规律。体应变的变化趋势说明峰值应力到达之前出现了剪缩(正的体应变)和在应变软化开始后的剪胀效应。结合图 6的平行粘结破坏发展趋势，可以看出剪胀效应在平行粘结破坏开始之后出现，因为粘结破坏以后，一部分颗粒就会遵循常规接触定律而发生转动，颗粒的变形和模型整体的变形也会相应增大，因而总应变也随之增大，这一现象与实际胶结砂岩的力学行为相似。

2) 配位数是衡量颗粒接触的一个重要指标，图 7描述了配位数随轴向应变的变化趋势：初始配位数变化缓慢，大约ε_a=1.5%时配位数出现峰值，ε_a=3%之后配位数减小的速率增加，ε_a>5%后配位数减小的速率趋缓。图中可以看出配位数总的变化并不大，说明剪切过程中大多数颗粒能充分接触，因为模型中的颗粒尺寸范围较小，受力过程中颗粒“漂浮”的几率较小，因此大多数颗粒都对模型的应力状态起作用。图 6显示平行粘结破坏于ε_a=2.5%后开始，破坏的粘结数在最大峰值应力后迅速增加，说明颗粒之间的粘结破坏越多，粘结破坏的颗粒就会发生转动，运动时的约束就相应减少，则应变软化也就越明显，与图 4~图 6得出的结论一致，当ε_a>5%后粘结破坏的速率逐渐减慢，此时模型在外力作用下的受力状态逐渐趋于稳定。

3) 图 8为Y-Z平面对应不同ε_a时的力接触网络(X∈(-0.001 m, 0.001 m))。剪切开始之前，其颗粒之间的接触为均匀分布，随着轴向应变的增大，颗粒之间的接触力随之增大。图 8(b)显示ε_a=2%时的接触力增大较多，与图 7对应的配位数发展趋势相吻合，说明此时的颗粒接触形成的力链结构充分发育，形成复杂的荷载传递路径；达到峰值应力之后，应变软化，部分力链方向发生变化和减弱，如图 8(c)所示ε_a=5%时，平行粘结破坏始于顶部和底部2个不同的位置，形成了明显的剪切带，而剪切带由于平行粘结的破坏导致颗粒转动约束随之消失，颗粒间的接触力也随之发生变化；随着应变逐渐增大到ε_a=8%、ε_a=10%和ε_a=12%，如图 8(d)、图 8(e)和图 8(f)所示，剪切带的宽度相应扩大，颗粒之间的粘结相应减少，力链传递荷载的能力也越来越弱，但相对于图 8(b)~图 8(d)的剪切带变化，此时的剪切带宽度变化较小而逐渐趋于稳定，直至模型最终趋于破坏。上述结果说明力链的改变能清楚地阐述图 4~图 7显示的胶结砂岩力学变化规律，因此力链结构体系在颗粒模型力学性质的研究中起着重要的作用。

因此，数值分析得到的胶结砂岩力学变化规律说明试样受剪时的非均匀响应能通过颗粒流数值模拟得到充分的体现。由此说明文中采用随机数学方法模拟砂岩颗粒的胶结分布合理，考虑平行粘结接触模型进行胶结砂岩力学分析可行，为后续不同胶结性状的储层砂岩研究提供了重要的研究手段。

2 胶结性状对砂岩力学特性的影响

储层砂岩具有不同的胶结性状，为了真实反映砂岩的结构破坏机理，现场试验或室内试验一般会针对实际的油藏岩心，但由于试验条件的限制和较高的造价导致岩样较少，模型单一化，得到成果通常适用于特定环境中的储层砂岩。而表征胶结砂岩特征的胶结性状对岩体结构受力起着重要的作用，那么不同胶结性状的研究能为描述石油开采过程中砂岩结构破坏机理和出砂机理奠定基础。

2.1 胶结程度对砂岩力学特性的影响

选用上述的颗粒流数值模型，保持其他参数不变，模拟胶结半径比变化时的砂岩力学特性如图 9~图 12所示。当胶结半径比大于0.5时，图 9表现更大的峰值应力比，对应的应变软化也更明显一些，说明胶结半径比越大，颗粒之间的约束越强，形成的力链网络越强，因此反映在模型上的应力就越大一些；当小于0.2时，由于平行粘结约束较弱，其力学特性接近于松砂，无明显的峰值应力比和应变软化形态，说明当胶结半径比越小时，模型的力学响应主要由摩擦接触控制，类似于无粘结的松砂模型。因此，胶结性状控制着胶结砂岩的力学特性，胶结程度越高，模型的初始刚度和峰值应力比也越大，说明实际胶结砂岩颗粒之间的约束和强度也越大，但是胶结程度半径比越大时，胶结越强，砂岩体的脆性也会越大。

2.2 胶结刚度对砂岩力学特性的影响

砂岩力学特性也会受到胶着物质刚度的影响，如图 13~图 16所示。图 13显示峰值应力比随着轴向应变的增加而增加，平行粘结刚度越大，最大峰值应力比对应的应变值越小，说明粘结刚度越大，颗粒刚性越强，形成的力链网络也越强，颗粒传递力的能力就越强，相同荷载作用下，粘结刚度越大的试样就会在更短的时间内达到较大的应力值，颗粒变形相对较小，因而对应的轴向应变值也较小。另外，平行粘结刚度越大，峰值应力之后其应变软化很明显，其软化率随着应变的增加而减小，特别是对于刚度较小的试样，基本上不会出现应变软化的情况。图 14描述的体应变随轴向应变的变化趋势说明粘结刚度越大，体积剪缩值越小，且对应的轴向应变也越小，对于平行粘结刚度越小的模型，颗粒一直处于剪缩状态，主要在于颗粒刚度越小颗粒传递外力的能力较弱，颗粒之间的粘结越容易破坏，此时颗粒压缩变形也就越大。

图 15描述的配位数变化规律说明初始配位数的增加缓慢，平行粘结刚度增大，其配位数达到峰值后的减小趋势也大，说明当平行粘结刚度越大，颗粒之间的刚性越强，脆性也越强，一旦颗粒发生平行粘结破坏，颗粒之间的转动约束就随之消失，整体的应力下降也就越多。图 16说明平行粘结刚度越小，颗粒刚性越低，其破坏的粘结数会越多，在相同的荷载作用下，承受张拉和剪切力的能力降低，所以破坏的粘结数就会增加更多，与图 13~图 15的结论一致。

2.3 胶结量对砂岩力学特性的影响

胶结量描述试样的颗粒胶结接触数占颗粒总接触数的比例，其比例的大小对砂岩力学特性的影响很大，如图 17~图 19所示。3种工况的应力比随着轴向应变的增大而增大，峰后应变软化率随应变的增大而减小。当胶结含量越大时，应力比的斜率越大，说明初始加载对应的弹性模量也愈大，模型刚性越强，得到的峰值应力也越大，但峰值后的应变软化更明显，也说明胶结含量越大，其砂岩体承载能力越强，但其脆性也越强。体应变随着轴向应变的趋势曲线说明胶结含量愈大，峰前剪缩越大，峰后软化也越强，也说明了胶结量越大岩体脆性也越强，与上述的结论一致。配位数的变化趋势同样也反映了胶结量越大，试样强度越大但是脆性越强的力学特性。

3 结论

文中建立了基于平行粘结接触的三维颗粒流数值模型研究胶结砂岩的力学特性，从岩体的宏细观力学响应验证了该数值模型的可行性，并分析了不同的胶结性状如胶结程度、平行粘结刚度和胶结量对砂岩力学特性的影响。

1) 基于随机数学理论和PFC3D的二次开发实现砂岩体的胶结性状分布，建立了PFC3D数值模型，通过描述应力比、体应变、配位数和粘结破坏数与轴向应变的关系阐述剪切过程胶结砂岩的力学特性，满足胶结砂岩体受力剪切破坏的一般规律；并从细观角度分析颗粒运动和颗粒接触形成的接触网络，其力接触网络的演化清晰地描述了试样剪切破坏过程，说明了力链在传力过程中所起的重要作用。由此也说明文中建立的数值模型是合理可行的，为建立储层砂岩剪切破坏机理提供了新的研究思路。

2) 比较了岩体胶结性状如胶结半径比、平行粘结刚度以及胶结量对胶结砂岩力学特性的影响，其成果显示上述因素的影响均较大。在胶结量一定的情况下，胶结半径比和平行粘结刚度越大，岩体的承载力越强，脆性也越强；前两因素一定的情况下，胶结量越大说明被约束的颗粒越多，岩体的承力能力越强，同样也会表现很强的脆性。由此可见，胶结性状对储层砂岩结构受力起着重要的控制作用，其研究成果为砂岩结构破坏机理的研究奠定了基础。那么，针对不同胶结性状的实际砂岩储层开采，需要建立相应的数值模型对砂岩力学响应规律进行验证和标定，才能确保数值分析模型在实际工程运用中的准确性和真实性。

文中提出的PFC3D数值模型试验可真实模拟储层砂岩宏细观力学响应和结构破坏机理。特别是对于规模和难度较大的现场或室内岩心试验，能突破常规模型在仪器设备能力和试验条件上的局限性，为石油开采过程中储层出砂机理、出砂预测及砂控策略的研究提供新的研究思路和手段。