陈铭，张士诚，柳明，马新仿，邹雨时，周彤，李宁，李四海

（1. 中国石油大学（北京），北京 102249；2. 中国石油长城钻探工程有限公司，北京 100101；3. 中国石化石油勘探开发研究院，北京 100083）

0 引言

水平井体积改造已成为页岩气等致密储集层油气开发的主导技术[1]。该技术常采用段塞加砂、小粒径支撑剂等工艺，缝内支撑剂多为局部支撑，且在主裂缝与分支缝交会处存在“转角支撑”现象。局部支撑、“转角支撑”导致支撑剂承受的外应力增大，使支撑剂向岩石缝面的嵌入更显著，从而引起产能递减甚至裂缝失效[2-3]。因此，体积改造背景下缝内支撑剂的嵌入行为研究及其定量分析具有重要意义。

支撑剂嵌入的研究方法包括实验研究、数值模拟和理论模型分析。许多学者针对不同岩石类型开展了支撑剂嵌入实验研究，并得到关于支撑剂嵌入的基本认识[4-8]。但实验研究通常针对某种特定类型岩石，研究结果难以推广。Alramahi等[9]采用有限元方法分析单个支撑剂颗粒与岩石表面的接触嵌入问题，发现支撑剂嵌入岩石表面主要发生于弹塑性变形阶段。Deng等[10]采用颗粒流数值模拟方法研究不同粒径支撑剂与岩石的相互作用，指出支撑剂浓度越大嵌入深度越小。虽然数值模拟方法可以求解嵌入过程，但计算量较大，不利于工程设计分析，因此目前应用较为广泛的仍是支撑剂嵌入理论模型分析方法。

理论模型经历了经验公式、半经验公式和理论解的发展历程。Huitt等[11]基于几何关系建立支撑剂嵌入深度计算公式，该公式中有两个参数需要通过实验研究来拟合确定。Volk等[12]研究支撑剂浓度、大小及分布对嵌入深度的影响，并提出基于实验结果的经验公式。Li等[13]采用双球体弹性接触模型建立了支撑剂嵌入的解析公式，该公式理论基础完善，但仅考虑了弹性变形，无法分析支撑剂弹塑性嵌入。

本文根据支撑剂嵌入岩石的力学过程建立支撑剂嵌入岩石的本构方程，该方程考虑支撑剂嵌入的弹性、弹塑性和塑性全过程。结合岩石-支撑剂体系的受力分析，给出支撑剂嵌入深度的计算方法，并将本文模型与弹性模型[13]及 Lacy等[5]、郭建春等[6]、卢聪等[7]实验结果进行对比，验证本文计算方法的可行性和准确性。最后基于本文计算方法，分析支撑剂嵌入深度的影响因素。

1 模型建立

1.1 岩石-支撑剂接触应力

取支撑剂与岩石的接触截面A（见图 1）进行分析，其截面积为At。在该截面上作用有流体压力和支撑剂接触应力[14]。根据受力平衡关系可得：

（1）式中接触应力表示与岩石接触的支撑剂承受的应力大小，为分析岩石-支撑剂体系内单个支撑剂的受力状态，需确定与裂缝面接触的支撑剂分布状态。裂缝内支撑剂分布包括 3种形式[15]：单层稀疏分布、单层紧密分布和多层分布（见图 2）。对于多层分布的支撑剂，受闭合应力影响，支撑剂在缝内处于压实状态，因此多层分布状态下与岩石接触的支撑剂可视为单层紧密分布。多层分布的支撑剂接触应力分析与单层紧密分布支撑剂接触应力分析相同。

图1 岩石-支撑剂受力分析图

图2 支撑剂分布形式示意图

采用支撑剂面密度描述支撑剂分布状态。支撑剂面密度为单位面积内与岩石接触的支撑剂数量，即：

受岩石表面粗糙度和支撑剂粒径分布影响，实际嵌入裂缝表面的支撑剂数量小于支撑剂总量[16-17]。定义支撑剂嵌入比例ε为支撑剂嵌入数量与单层支撑剂（与岩石表面接触的支撑剂）总量之比，即：

单个支撑剂平均嵌入面积为：

联立（1）式—（4）式，得到单个支撑剂平均接触应力为：

（5）式描述了接触应力与嵌入面积、嵌入量及外载荷的关系。由于外载荷、流体压力和支撑剂分布为已知量，接触应力的求解关键在于嵌入面积的确定。

1.2 支撑剂嵌入本构方程

支撑剂嵌入岩体的过程可描述为圆球在法向载荷作用下与半无限空间体的接触力学问题。随着支撑剂嵌入深度的增加，岩体和支撑剂接触位置的局部变形经历 3种阶段——弹性阶段、弹塑性阶段和塑性阶段（见图3）。建立支撑剂嵌入过程的本构模型，包括嵌入深度与外载荷的关系、嵌入深度与接触面积的关系，是分析支撑剂嵌入深度的重要步骤。其中接触面积为支撑剂嵌入岩石的横截面积。

图3 支撑剂嵌入阶段示意图

1.2.1 弹性嵌入

支撑剂与岩石的接触面积较小时，局部变形满足弹性变形。根据 Hertz理论[18]，支撑剂嵌入岩石的深度与平均接触应力的关系为：

广义弹性模量满足如下关系式[18]：

弹性嵌入阶段，支撑剂嵌入岩石的接触面积与嵌入深度的关系为[19]：

当支撑剂与岩石的接触应力达到 1.1倍岩石屈服应力时，岩石开始发生弹塑性变形[19]。因此，根据（6）式可得弹性阶段最大嵌入深度为：

值得注意的是，根据（9）式确定的弹性阶段最大嵌入深度是理论值。

1.2.2 塑性嵌入

由于弹塑性阶段的力学特征参数值主要以数值解为主，解析公式需要根据弹性和塑性解进行插值建立，因此首先描述塑性阶段变形，再给出弹塑性嵌入的本构方程。

当接触应力达到岩石硬度时，岩石开始发生塑性变形[20]，并维持恒定的接触应力。材料硬度一般为屈服应力的3倍[21]，因此塑性阶段接触应力为：

塑性阶段支撑剂嵌入岩石的接触面积与嵌入深度的关系为[22]：

岩石发生塑性变形的临界嵌入深度的准确确定需借助纳米压痕实验。为方便理论计算，本文采用固体材料发生塑性变形的临界嵌入深度计算方法进行计算。对于固体材料，理论和实验研究表明，开始发生塑性变形的临界嵌入深度一般为弹性阶段最大嵌入深度的82.5倍[23-24]。因此，塑性变形的临界嵌入深度为：

需要注意的是，根据（12）式确定的发生塑性变形的临界嵌入深度是基于均质固体材料假设的理论结果。对于实际岩石，考虑其沉积演化等复杂过程的影响，发生塑性变形的临界嵌入深度可能会有变化，准确值需通过纳米压痕实验确定。

1.2.3 弹塑性嵌入

当dc1≤d≤dc2时，接触区域处于弹塑性变形阶段[25-27]。弹塑性变形阶段属于纯弹性和塑性变形的过渡阶段，但力学机理较为复杂，难以得到理论解。采用有限元等数值方法可以获得数值解[21]，但数值方法计算量较大，不利于工程应用。因此，本文采用插值方法建立弹塑性阶段的本构方程，方法简单，便于计算分析。

根据Francis[28]的实验研究结果，弹塑性阶段接触应力可描述为：

根据[dc1,dc2]端点的应力连续性条件可以确定系数a1、a2，由此得到：

由于弹塑性阶段为弹性和塑性的过渡阶段，因此可假设弹塑性阶段接触面积与嵌入深度的关系为：

采用三次样条插值法确定c(d)的表达式，令：

连续性条件为：

求解方程组（17）得到：

将（18）式代入（15）式可得：

为分析插值方法的可靠性，将插值法计算结果与Kogut等[27]的有限元法计算结果对比（见图 4），可以发现两者较为接近，平均误差为5.71%，说明插值方法具有较高的工程计算精度。

图4 插值结果和有限元结果对比

1.3 支撑剂嵌入深度求解流程

支撑剂嵌入深度的计算方法取决于支撑剂与岩石的嵌入变形阶段，而变形阶段也是待确定的，因此采用试算法求解支撑剂嵌入深度。首先根据岩石和支撑剂基本参数计算临界嵌入深度dc1和dc2。假设支撑剂嵌入处于弹性阶段，则联立（5）式、（6）式和（8）式，得到关于d的非线性方程F1(d)，再采用牛顿法求得弹性嵌入深度，若其值小于等于dc1，则假设正确。否则假设支撑剂嵌入处于弹塑性阶段，联立（5）式、（14）式和（19）式，得到关于d的非线性方程F2(d)，再采用牛顿法求得弹塑性嵌入深度，若其值小于dc2，则假设正确。否则支撑剂嵌入处于塑性阶段，联立（5）式、（10）式和（11）式，得到关于d的线性方程F3(d)，从而求得支撑剂嵌入深度。采用Matlab软件编程实现上述求解过程。

2 模型验证

为验证本文模型准确性，将本文模型、弹性模型[13]分别与 Lacy等[5]、郭建春等[6]和卢聪等[7]实验结果对比（见图5）。

弹性模型为Li等[13]推导的支撑剂嵌入深度计算模型，该模型基于双球体弹性接触理论解建立。与实验结果对比发现，弹性模型计算值小于实验值，因此引入拟合系数修正弹性模型。然而，针对不同支撑剂和岩石，需要结合实验结果确定不同的拟合系数，且拟合系数的作用是修正弹性解，会造成对支撑剂嵌入的非弹性力学机制的忽视。而本文模型未引入拟合系数，避免了这些问题。由图 5可知，本文模型计算结果与修正后弹性模型及实验结果较为接近，验证了本文模型的准确性。因此，本文模型不但更符合支撑剂嵌入的力学机制，而且不需要借助实验引入修正系数，计算结果也比较准确，具有可行性和准确性。

图5 本文模型、弹性模型、修正弹性模型与实验结果对比

3 支撑剂嵌入深度影响因素

利用本文支撑剂嵌入深度计算方法，以昌吉油田芦草沟组致密油储集层 JHW020井参数为基础数据进行支撑剂嵌入深度影响因素分析。该井岩石和支撑剂基本参数为：岩石弹性模量为35 210 MPa，岩石泊松比为 0.22，岩石屈服应力为 120 MPa；支撑剂粒径为0.595/0.297 mm（30/50目），取调和平均值0.362 mm，支撑剂弹性模量为41 306 MPa，支撑剂泊松比为0.25，支撑剂密度为2.8 kg/m3；闭合应力为60 MPa，生产过程裂缝内流体压力为15 MPa。在对某一影响因素进行分析时，只对该影响因素取不同值，其他参数均使用基本参数值。

3.1 支撑剂铺置浓度

支撑剂铺置浓度为单位面积的支撑剂质量，根据支撑剂颗粒密度和体积可得到支撑剂铺置浓度与面密度的转化关系：

本文所用支撑剂达到最大面密度时的支撑剂铺置浓度为0.95 kg/m2。当铺置浓度大于0.95 kg/m2时，支撑剂呈现多层分布形式，但与岩石表面接触的支撑剂面密度仍为单层紧密铺置时的最大面密度，因此继续增大浓度不会增加与岩石表面接触的承压支撑剂数量。如图6所示，当支撑剂浓度大于1.00 kg/m2时，支撑剂嵌入深度和接触应力不再发生变化，与前述分析一致。图6a显示，不同缝内流体压力条件下，支撑剂浓度0.25 kg/m2时的嵌入深度是1.00 kg/m2时的嵌入深度的2～3倍，即当支撑剂为单层稀疏分布时，支撑剂浓度越大，支撑剂嵌入深度越小。

图6 不同铺置浓度的支撑剂嵌入特征

图6b显示，在本文算例条件下，不同支撑剂浓度下的接触应力均大于岩石屈服应力，并小于 3倍岩石屈服应力，表明支撑剂嵌入主要发生于弹塑性变形阶段。同时，流体压力为0～55 MPa时，接触应力变化缓慢，流体压力为55～65 MPa时，接触应力变化显著。在闭合应力60 MPa、流体压力为零的条件下，不同支撑剂浓度下的接触应力为287～315 MPa，平均值约为闭合应力的 5倍。由于支撑剂嵌入深度较小，接触面积远小于闭合应力作用面积，根据（1）式可以得到接触应力远大于闭合应力。需要注意的是，实验测试支撑剂破碎时的应力为闭合应力，并不是接触应力，因此接触应力较高并不能说明支撑剂发生破碎。

3.2 支撑剂嵌入比例

图7 不同嵌入比例的支撑剂嵌入特征

图 7中支撑剂嵌入深度与接触应力随流体压力的变化规律与图 6类似。相同条件下，支撑剂嵌入深度和接触应力随嵌入比例减小而增大。支撑剂嵌入比例为 1时，承压支撑剂密度最大，单个颗粒的接触应力最小，因此嵌入深度最小。嵌入比例受支撑剂粒径分布影响，粒径分布分散时嵌入比例较小，而粒径分布集中的支撑剂嵌入比例会趋近于1，因此采用粒径均一或分布相对集中的支撑剂有利于减小支撑剂嵌入深度。对于段塞加砂工艺，支撑剂分布会出现局部支撑，嵌入比例会减小，从而增大了嵌入深度。因此，段塞加砂工艺要注重后期支撑剂嵌入的分析，在压裂设计阶段可以提高支撑剂浓度来减小嵌入深度。

3.3 闭合应力与缝内流体压力

当流体压力大于等于闭合应力时，支撑剂接触应力为零，不发生嵌入；随着流体返排和生产进行，裂缝内流体压力不断减小，作用在支撑剂上的接触应力不断增加，嵌入深度逐渐增大。图 8显示，流体压力为零的条件下，与闭合应力为50 MPa时相比，闭合应力为65 MPa时嵌入深度增大25.6%；在本文算例中闭合应力条件（60 MPa）下，与缝内流体压力为40 MPa时相比，缝内流体压力为零时嵌入深度增大127.5%。因此，高闭合应力情况（如深井压裂）下支撑剂嵌入会更加显著；生产过程中流体压力下降会显著增大支撑剂嵌入深度，在进行产能递减分析时不能忽视流体压力下降导致的支撑剂嵌入问题。

3.4 支撑剂粒径

图 9a显示，相同条件下，减小粒径会增大支撑剂嵌入深度。图 9b显示，不同粒径的支撑剂主要发生弹塑性嵌入，接触应力受支撑剂粒径影响，粒径越小，嵌入面积越小，接触应力则越大。因此，采用小粒径支撑剂时需要注意支撑剂嵌入引起的缝宽变化。

3.5 岩石与支撑剂弹性模量

图 10a显示，减小岩石与支撑剂弹性模量比值会增大支撑剂嵌入深度。对于煤岩、塑性页岩等，较低的弹性模量势必带来较大的嵌入深度，与实验结果[4-5,29]一致。图 10b显示，相同条件下，岩石与支撑剂弹性模量比值越大，接触应力越大。这是因为相同外载与内压作用下，支撑剂所受的接触应力与接触面积成反比，岩石弹性模量较低时支撑剂嵌入深度较大，接触面积较大，接触应力减小。

图8 不同闭合应力的支撑剂嵌入特征

图9 不同粒径的支撑剂嵌入特征

图10 不同岩石与支撑剂弹性模量比值的支撑剂嵌入特征

4 讨论

本文支撑剂嵌入深度计算模型可引入到导流能力计算公式中，从而进一步分析导流能力变化。尽管本文对嵌入过程的描述更加接近实际，但本文模型还存在以下主要问题。

①支撑剂嵌入的分析忽略了由于闭合应力过高等因素导致的支撑剂破碎问题。而支撑剂破碎会导致嵌入裂缝面的支撑剂数量和比例减小，从而加剧支撑剂嵌入[30]。

②在岩石-支撑剂体系的受力分析中，本文引入了支撑剂嵌入比例，仅分析嵌入比例对嵌入深度和接触应力的影响，未深入探讨该参数的确定方法。实际上，嵌入比例受支撑剂粒径分布和岩石表面凹凸形貌影响，确定其取值需对岩石表面形态进行准确刻画。

5 结论

基于岩石-支撑剂体系受力特征和支撑剂嵌入本构方程提出的支撑剂嵌入深度计算方法考虑了支撑剂弹性—弹塑性—塑性嵌入全过程，与弹性模型相比更加适用于支撑剂嵌入的理论分析和计算，且计算结果较准确。

支撑剂主要发生弹塑性嵌入。支撑剂单层分布时，铺置浓度越大，接触应力越小，嵌入深度越小；多层分布时，与岩石接触的承压支撑剂紧密分布，接触应力和嵌入深度不再随铺置浓度增大而变化。采用粒径均一或分布相对集中的支撑剂有利于减小支撑剂嵌入深度。相同条件下，支撑剂嵌入深度随闭合应力增大而增大，随缝内流体压力下降而增大。支撑剂粒径及岩石与支撑剂弹性模量比值越小，嵌入深度越大。

符号注释：

a1，a2，c1，c2，c3，c4——待定系数；A——支撑剂与岩石的接触截面；Ap——支撑剂与岩石的接触面积，m2；Ape——单个支撑剂平均嵌入面积，m2；Ap,e，Ap,p，Ap,ep——弹性阶段、塑性阶段和弹塑性阶段支撑剂与岩石的接触面积，m2；′——弹性阶段、塑性阶段和弹塑性阶段接触面积对嵌入深度的一阶导数，m；At——截面A总面积，m2；c(d)——待定函数，1≤c≤2；cp——支撑剂铺置浓度，kg/m2；d——支撑剂嵌入岩石的深度，m；dc1——弹性阶段最大嵌入深度，m；dc2——塑性变形的临界嵌入深度，m；E′——广义弹性模量，MPa；Ep，Er——支撑剂、岩石的弹性模量，MPa；F1(d)，F2(d)，F3(d)——关于d的方程；n——单位面积内支撑剂嵌入数量，个/m2；nt——支撑剂面密度，个/m2；N——嵌入岩石的支撑剂总量；Nt——与岩石接触的支撑剂总数量；pf——流体压力，MPa；R——支撑剂颗粒半径，m；ε——支撑剂嵌入比例，无因次；ρ——支撑剂颗粒密度，kg/m3；σc——闭合应力，MPa；σp——支撑剂接触应力，MPa；σp,e，σp,p，σp,ep——弹性阶段、塑性阶段和弹塑性阶段的接触应力，MPa；σy——岩石屈服应力，MPa；υp，υr——支撑剂、岩石的泊松比。

[1]胥云, 陈铭, 吴奇, 等. 水平井体积改造应力干扰计算模型及其应用[J]. 石油勘探与开发, 2016, 43(5): 780-786, 798.XU Yun, CHEN Ming, WU Qi, et al. Stress interference calculation model and its application in volume stimulation of horizontal wells[J]. Petroleum Exploration and Development, 2016, 43(5):780-786, 798.

[2]CHEN D, YE Z H, PAN Z J, et al. A permeability model for the hydraulic fracture filled with proppant packs under combined effect of compaction and embedment[J]. Journal of Petroleum Science and Engineering, 2016, 149: 428-435.

[3]CUI A, GLOVER K, WUST R. Elastic and plastic mechanical properties of liquids-rich unconventional shales and their implications for hydraulic fracturing and proppant embedment: A case study of the Nordberg Member in Alberta Canada[R].Minneapolis: 48th US Rock Mechanics/Geomechanics Symposium,2014.

[4]LACY L, RICKARDS A, ALI S. Embedment and fracture conductivity in soft formations associated with HEC, borate and water-based fracture designs[R]. SPE 38590, 1997.

[5]LACY L, RICKARDS A, BILDEN D. Fracture width and embedment testing in soft reservoir sandstone[J]. SPE Drilling & Completion,1998, 13(1): 25-29.

[6]郭建春, 卢聪, 赵金洲, 等. 支撑剂嵌入程度的实验研究[J]. 煤炭学报, 2008, 33(6): 661-664.GUO Jianchun, LU Cong, ZHAO Jinzhou, et al. Experimental research on proppant embedment[J]. Journal of China Coal Society,2008, 33(6): 661-664.

[7]卢聪, 郭建春, 王文耀, 等. 支撑剂嵌入及对裂缝导流能力损害的实验[J]. 天然气工业, 2008, 28(2): 99-101.LU Cong, GUO Jianchun, WANG Wenyao, et al. Experimental research on proppant embedment and its damage to fractures conductivity[J]. Natural Gas Industry, 2008, 28(2): 99-101.

[8]TAN Y, PAN Z, LIU J, et al. Experimental study of permeability and its anisotropy for shale fracture supported with proppant[J]. Journal of Natural Gas Science & Engineering, 2017, 44: 250-264.

[9]ALRAMAHI B, SUNDBERG M. Proppant embedment and conductivity of hydraulic fractures in shales[R]. ARMA 12-291,2012.

[10]DENG S, LI H, MA G, et al. Simulation of shale-proppant interaction in hydraulic fracturing by the discrete element method[J].International Journal of Rock Mechanics & Mining Sciences, 2014,70(9): 219-228.

[11]HUITT J, MCGLOTHLIN B. The propping of fractures in formations susceptible to propping-sand embedment[R]. API 58-115, 1958.

[12]VOLK L, RAIBLE C, CARROLL H, et al. Embedment of high strength proppant into low-permeability reservoir rock[R]. SPE 9867,1981.

[13]LI K W, GAO Y P, LU Y C, et al. New mathematical models for calculating proppant embedment and fracture conductivity[J]. SPE Journal, 2015, 20(3): 496-507.

[14]ZHANG J J, KAMAMENOV A, ZHU D, et al. Laboratory measurement of hydraulic fracture conductivities in the Barnett shale[R]. SPE 163839-PA, 2014.

[15]CHEN C, MARTYSEVICH V, OCONNELL P, et al. Temporal evolution of the geometrical and transport properties of a fracture/proppant system under increasing effective stress[J]. SPE Journal, 2015, 20(3): 527-535.

[16]MUELLER M, AMRO M. Indentation hardness for improved proppant embedment prediction in shale formations[R]. SPE 174227,2015.

[17]AWOLEKE O, ZHU D, HILL A D. New propped-fracture-conductivity models for tight gas sands[J]. SPE Journal, 2016, 21(5): 1-10.

[18]徐芝纶. 弹性力学简明教程[M]. 4版. 北京: 高等教育出版社,2013: 228-230.XU Zhilun. Short course of elasticity[M]. 4th ed. Beijing: Higher Education Press, 2013: 228-230.

[19]JOHNSON K L. Contact mechanics[M]. Cambridge: Cambridge University Press, 1985.

[20]BIG-ALABO A, HARRISON P, CARTMELL M P. Contact model for elastoplastic analysis of half-space indentation by a spherical impactor[J]. Computers & Structures, 2015, 151: 20-29.

[21]SONG Z, KOMVOPOULOS K. Elastic-plastic spherical indentation:Deformation regimes, evolution of plasticity, and hardening effect[J].Mechanics of Materials, 2013, 61(8): 91-100.

[22]BARTIER O, HERNOT X, MAUVOISIN G. Theoretical and experimental analysis of contact radius for spherical indentation[J].Mechanics of Materials, 2010, 42(6): 640-656.

[23]STRONGE W J. Impact mechanics[M]. Cambridge: Cambridge University Press, 2000.

[24]STRONGE W J. Contact problems for elasto-plastic impact in multi-body systems[J]. Lecture Notes in Physics, 2000, 551:189-234.

[25]ZHAO Y W, MAIETTA D M, CHANG L. An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow[J]. Journal of Tribology, 2000, 122(1): 86-93.

[26]BRAKE M R. An analytical elastic-perfectly plastic contact model[J].International Journal of Solids & Structures, 2012, 49(22):3129-3141.

[27]KOGUT L, ETSION I. Elastic-plastic contact analysis of a sphere and a rigid flat[J]. Journal of Applied Mechanics, 2002, 69(5):657-662.

[28]FRANCIS H A. Phenomenological analysis of plastic spherical indentation[J]. Journal of Engineering Materials & Technology Transactions of the ASME, 1976, 98(3): 272-281.

[29]ZHANG J J, OUYANG L, HILL A D, et al. Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in shale reservoirs[R]. SPE 170775, 2014.

[30]HAN J H, WANG J Y, PURI V. A fully coupled geomechanics and fluid flow model for proppant pack failure and fracture conductivity damage analysis[R]. SPE 168617, 2014.