基于离散元方法的层状发育页岩中断裂扩展规律研究 卢明静， 郑雪琳， 王东英， 王康， 杨峰， 张子霖

《Processes》：Fracture Propagation Laws in Lamina-Developed Shale Based on the Discrete Element Method Mingjing Lu, Xuelin Zheng, Dongying Wang, Kang Wang, Feng Yang and Zilin Zhang

　　摘要：中国东部大陆断裂盆地中的页岩油，以济阳凹陷为代表，在产能方面取得了突破性进展。然而，深埋藏、高地层压力和原油流动性差等挑战对实现高产和稳产构成了重大障碍。为了形成复杂的裂缝网络进行刺激，需要进行水力压裂。页岩的层状结构、地质力学条件以及压裂操作参数等因素都会影响裂缝的传播。因此，本研究使用离散元软件PFC2D 6.0建立了层状页岩中裂缝传播的数值模型，在原位应力条件下进行了裂缝传播规律的模拟分析，并研究了层状结构和施工技术对裂缝复杂性的影响。结果表明，对于层状页岩，启动压力随着注射速率的增加而降低；当两个水平主应力之间的差异增大时，水力裂缝逐渐倾向于向最大主应力方向传播；在高注射压力下，形成了复杂的短裂缝网络，而在低注射压力下，主裂缝的长度会增加。高密度（9-10条/100毫米）有助于层间的穿透，促进裂缝向最大水平主应力方向延伸；低密度（4-5条/100毫米）则加强了层间的引导作用，裂缝沿层状结构靠近注射孔处传播。这项研究阐明了层状页岩中裂缝的启动和传播机制，为优化水力压裂设计提供了理论和技术支持。

1. 引言
济阳凹陷的页岩油在产能方面取得了重大突破，极大地保障了能源供应[1,2]。然而，中国的陆上页岩油储层的特点是横向异质性强、层间胶结度低以及多层垂直叠加[3,4,5,6]。此外，低孔隙度和渗透率、高地层压力以及原油流动性差等挑战对实现高产和稳产构成了重大障碍[7,8,9]。研究表明，支撑裂缝和激活裂缝的耦合可以形成跨越纳米到厘米尺度的“双重裂缝”高渗透率流道，这有助于提高页岩油和天然气储层的开发效率[10,11]。济阳凹陷的页岩层状结构发育良好，层间的开口对于页岩发生拉伸-剪切破坏并形成复杂裂缝至关重要。
目前，关于水力压裂的研究主要依赖于理论分析、物理实验和数值模拟[12,13,14]。在理论和模型构建方面，已经建立了基于弹性理论和应力场分析的裂缝传播准则[15,16,17,18]，并提出了定量指标[19]。Hubbert等人提出了经典的裂缝传播垂直于最小主应力的理论[3]。Liu等人结合了多应力场来改进垂直裂缝准则[20]。Huang等人明确了“当最大拉应力达到抗拉强度时发生启动”的判断标准[21]。在物理实验方面，岩石力学实验[22,23]和水力压裂实验[24,25,26]与CT可视化[27,28,29]及正交实验设计[30]相结合，用来验证理论和模拟结果。由于物理实验条件和监测设备的限制，大多数研究仍仅涉及宏观现象。同时，进行实验需要耗费大量人力和时间，而数值模拟方法弥补了物理实验的不足。因此，有限元和离散元方法被用来模拟裂缝传播并分析多种因素的影响[31,32]。
有限元方法通常使用 cohesive 元素、joint 元素和网格重划分来模拟层状岩石的力学特性和裂缝传播[31,32,33]。研究揭示了页岩 bedding 平面对应力场强度和变形的影响。页岩的断裂韧性随 bedding 方向的不同而变化，垂直于 bedding 方向的断裂韧性高于平行于 bedding 方向的断裂韧性。Bedding 平面影响岩石的破坏模式（剪切破坏容易沿 bedding 平面发生），进而影响裂缝形态。然而，FEM 模拟受到裂缝尖端应力场的奇异点的影响，导致可模拟的裂缝数量有限。相比之下，离散元方法（DEM）根据粒子的接触状态来确定接触的存在与否[34,35,36]。使用 PFC 建立了流体-固体耦合模型，将粒子包围的区域视为“域”，粒子间的接触视为“管道”。该模型能够观察裂缝在介观尺度上的启动和传播动态，便于研究页岩的渗透特性和水力裂缝的扩展。尽管 DEM 已广泛应用于各种常规岩石实验，但对不同工作条件下页岩表现出切向破坏的介观特性缺乏综合分析。
优化工艺参数为裂缝传播规律提供了定量基础。在低粘度和低注射速率下，裂缝沿 bedding 均匀分布，而高位移会加速裂缝传播并增加裂缝的高度/宽度[33,34]。Feng 等人指出，注射压力的过度增加会抑制主裂缝的产生并促进次级裂缝的形成[5]。此外，岩石性质和地质环境是决定裂缝传播模式的核心因素[35]。岩石的脆性和天然裂缝对水力压裂传播和复杂裂缝的形成极为重要。Ding 等人提出了“裂缝发育系数”来量化裂缝发育的程度[36]。更强的岩石脆性和更高的天然裂缝密度有助于形成复杂裂缝，而较低的原位应力差更有利于复杂裂缝的发展[37,38,39]。地质环境对裂缝传播至关重要[40,41]。Xu 等人表明，小的水平应力差和适当的层厚可以优化压裂效果[15]。Wang 等人发现，更长的水化时间会增加裂缝复杂性[7]。Gunaydin 确认裂缝间距的均匀性会影响相互抑制[33]。页岩具有发育良好的 bedding 和显著的各向异性，这进一步增加了预测复杂裂缝传播的难度[42,43]。Bedding 平面可以改变裂缝轨迹，而 bedding 厚度、间距和应力差决定了裂缝是否沿 bedding 传播[44,45,46]。Arora 等人发现，层状页岩沿多个剪切平面发生破坏[27]，而大陆页岩中的诱导裂缝沿天然裂缝和 bedding 平面扩展[47]。
与常规页岩相比，济阳凹陷的页岩具有更高的层状发育程度，这对复杂裂缝的形成影响更大。水力压裂过程中页岩层状的不连续性以及页岩水力裂缝的微观力学特性可以通过离散元方法更好地描述。同时，可以更直观地观察水力压裂的过程和结果。因此，使用离散元方法建立了层状页岩的力学模型，并验证了其合理性。然后，建立了层状页岩的水力压裂传播模型。考虑了注射速率、注射压力、原位应力状态和 bedding 密度的影响，对层状页岩中的水力压裂传播进行了数值模拟，以阐明裂缝启动、传播和裂缝网络形成的机制，以及 bedding 平面的影响机制。本研究可以为页岩水力压裂提供理论和技术支持。

2. 层状页岩储层的裂缝传播模型建立
2.1. 页岩压裂启动和传播模型
裂缝启动遵循拉伸破坏准则和剪切破坏准则。拉伸破坏准则为：
(1)
其中 σ_t 是页岩的抗拉强度，单位为 MPa；σ_h1 和 σ_h2 是最大和最小水平主应力，单位为 MPa。
剪切破坏准则采用 Mohr-Coulomb 准则：
(2)
其中 φ 是页岩的内摩擦角，单位为度；c 是页岩的粘聚力，单位为 MPa；σ_n 是正应力，单位为 MPa；σ_s 是剪切应力，单位为 MPa。
实际情况中页岩的位置关系如图 1 所示，井筒周围岩石所受的力如图 2 所示。
图 1. 井筒与页岩储层的关系。
图 2. 井筒周围岩石的力图。
以压缩应力为正，拉伸应力为负，水平井井筒上三个主应力的公式如下：
(3)
其中 σ_h1 是井筒周围的最大主应力，单位为 MPa；σ_m 是中间主应力，单位为 MPa；σ_h2 是井筒周围的最小主应力，单位为 MPa；σ_r 是井筒的径向正应力，单位为 MPa；σ_a 是井筒的轴向正应力，单位为 MPa；σ_h 是井筒的最大水平原位应力，单位为 MPa；σ_h2 是井筒的最小水平原位应力，单位为 MPa；σ_s 是剪切应力，单位为 MPa；R 是井筒半径，单位为 m；P 是井筒内的注射压力，单位为 MPa；p 是孔隙压力，单位为 MPa。
当井筒的破坏模式为拉伸破坏时，最大拉应力为：
(4)

2.2. 层状页岩的裂缝传播模型
为了发生拉伸破坏，裂缝尖端的应力强度因子必须满足：
(5)
假设基质中线性裂缝传播所需的临界水压力为 φ_w，则：
(6)
其中 E 是目标页岩基质的弹性模量，单位为 GPa；l 是裂缝的半长，单位为 m；γ 是目标页岩基质的单位面积表面能，单位为 MPa·m。
裂缝扩展到层状结构时的应力状态如图 3 所示。
图 3. 裂缝扩展到层状结构。
在外场应力作用下，层状结构上的正应力和剪切应力可以表示为：
(7)
(8)
根据二维弹性理论，在内部水压力和外地应力作用下，层状结构上的正应力和剪切应力可以表示为：
(9)
(10)
在计算过程中，当层状平面处于开放状态时；当 σ_s > σ_h1 时，会发生剪切破坏，可以得到裂缝传播过程中的压力。如果设 σ_s > ρφ，则：
(11)
根据图 3 的分析，当水力裂缝扩展到层状结构时，裂缝传播的方向受层状界面强度和原位应力驱动效应的影响。根据 σ_s 和 σ_h 的关系，可以确定遇到层状结构后的裂缝方向。

3. 层状页岩水力压裂数值模型的建立和验证
页岩数值模型试样的介观参数详细列于表 1 中。
表 1. 页岩数值模型试样的介观参数。
使用 PFC2D 6.0 软件将模型直径设置为 50 mm，长度设置为 100 mm。页岩颗粒大小范围为 0.08 mm 至 0.13 mm，密度为 2650 kg/m3，接触模型设置为线性接触模型。此外，使用“contact model”命令将粒子间的接触转换为平行键合接触，建立了包含 20,353 个粒子的层状页岩试样的二维数值模型。在模拟中，壁面位移率为 0.05 m/s，需要通过 FISH 函数实现对页岩模型试样的三轴加载，包括围压和速度。图中绿色代表基质，红色代表层状结构。根据表 1 中的参数，建立了层状角度分别为 0°、45° 和 90° 的数值模型，如图 4 所示。
图 4. 不同层状倾角下的页岩数值模型。
进一步在恒定围压 50 MPa 下对数值模型进行了数值模拟，并获得了试样的应力-应变曲线，并与实验结果进行了比较，如图 5 所示。
图 5. 层状倾角为 0°、45° 和 90° 的数值模拟曲线比较。
从图 5 的比较结果可以看出，在数值模拟和实验中，层状倾角为 90° 时的最大偏应力值很低，而 0° 时的曲线显示出显著增加，最大值相对较大。数值模拟和实验显示了相同的模式，表明数值模拟与观察到的规律一致。如图6所示，建立了一个层状页岩的水力裂缝传播模型。该模型的尺寸和层参数参考了吉阳凹陷锦丰中央区块的层状页岩参数。数值样本模型的边长设置为100毫米，每层与水平方向之间的角度分别为0°、15°、45°、60°、75°和90°。为了更好地研究和模拟页岩水力裂缝的传播规律，需要对模型做出以下假设：(1) 页岩层理分布均匀且连续。实际地层中可能存在的层理缺失、局部不连续性或随机分布的微裂缝被忽略。层理的力学各向异性仅通过预设的颗粒/块体接触参数来模拟。(2) 页岩岩性是均匀的。实际地层中存在的矿物组成微观异质性对裂缝传播的影响被忽略。(3) 页岩是一种线性弹性材料。页岩在高应力或水力作用下的粘弹性/粘塑性变形被忽略。变形特性仅用颗粒/块的弹性模量和泊松比来描述。(4) 页岩的力学参数是恒定且各向同性的。实际地层中参数随深度和应力状态的动态变化被忽略。(5) 页岩的破坏遵循预设的介观裂缝准则，准则参数与宏观力学测试的结果相匹配。层间矿物颗粒之间的复杂断裂机制被忽略。(6) 破裂液是不可压缩流体。其体积变化在高压下对泵送压力和裂缝传播速率的影响被忽略。(7) 破裂液在裂缝中的流动遵循达西定律。高速流动下的惯性力和非达西效应被忽略，页岩基质中的流体渗滤/过滤也未被考虑。(8) 页岩被离散化为一系列规则颗粒，颗粒的大小和形状被设定为预设值。实际岩石的不规则介观结构被忽略。(9) 颗粒间的接触相互作用仅通过预设的本构模型传递。颗粒间的微观力（如范德瓦尔斯力和电场力）被忽略。(10) 破裂液温度变化和页岩的热膨胀/收缩对裂缝传播的影响被忽略。

4. 水力裂缝传播规律及其在层状页岩中敏感因素的分析
4.1. 原地应力影响规律分析
设置了四组最小水平主应力，分别为50 MPa、53 MPa、56 MPa和59 MPa，最大水平主应力固定为62 MPa。对于每组，配置了0°、15°、45°、75°和90°的层理倾角，以研究启动压力与原地应力之间的关系，如图7所示。从图7可以看出，在不同原地应力下，不同层理倾角的启动压力非常相似，差异不超过6%，表明启动压力受层理倾角的影响较小。综合比较四组数据发现，启动压力随围压的增加而增加，差异超过25%。原地应力对启动压力的影响大于层理倾角。
不同原地应力下的裂缝传播模式如图8所示。可以看出，随着水平主应力差的减小，水力裂缝逐渐偏向最小水平主应力方向。在12 MPa的应力差（最小应力50 MPa）下，裂缝倾向于沿最大主应力方向传播；而在6 MPa的应力差（最小应力56 MPa）下，这种方向性减弱，裂缝向最小应力方向的延伸更为明显。对于小倾角（0°、15°），裂缝在遇到层理界面时受到显著影响，倾向于沿层理延伸，不受应力差变化的影响。对于大倾角（45°、75°、90°），裂缝更受原地应力驱动，容易穿透层理并向最大主应力方向延伸，在较大应力差下穿透比例更高。当应力差为3 MPa（最小应力59 MPa）时，主裂缝的平均长度约为41.2 ± 2.8毫米，次级裂缝数量为17 ± 2条，裂缝网络复杂度指数为0.85 ± 0.04，层理穿透率为65% ± 3%。当应力差增加到12 MPa（最小应力50 MPa）时，主裂缝的平均长度增加到57.6 ± 3.5毫米（相对增加39.8%），次级裂缝数量减少到6 ± 1条（相对减少64.7%），裂缝网络复杂度指数降至0.31 ± 0.03（相对减少63.5%），层理穿透率上升至92% ± 2%（相对增加41.5%）。小的应力差（3 MPa，最小应力59 MPa）导致裂缝形态更复杂，尤其是在75°倾角时，分支更多，方向分散。大的应力差（12 MPa，最小应力50 MPa）降低了网络复杂度，主裂缝沿最大主应力方向集中，尽管15°倾角下的裂缝仍然相对复杂。

4.2. 注射速率影响规律分析
最小主应力设定为56 MPa，最大水平主应力设定为62 MPa。裂缝液的粘度为5 mPa·s，注射速率分别设定为0.01 m/s、0.03 m/s、0.06 m/s和0.09 m/s。裂缝传播的模拟结果如图9所示。如图9所示，裂缝在注入孔周围形成裂缝网络，且裂缝受到层理的显著影响。在层理倾角为0°和15°时，裂缝遇层理时沿层理平面方向传播。在层理倾角为45°、60°、75°和90°时，裂缝主要沿层理方向延伸，但也有一些裂缝穿过层理。随着注射速率从0.01 m/s增加到0.09 m/s，裂缝网络的复杂度先增加后减少，其中0.03 m/s的注射速率下观察到最复杂的裂缝形态。

4.3. 注射压力影响规律分析
最小主应力设定为56 MPa，最大水平主应力设定为62 MPa。裂缝液粘度为5 mPa·s，注射压力分别设定为60 MPa、62 MPa、64 MPa和66 MPa。裂缝传播的情况如图11所示。从图11可以看出，在低注射压力（例如60 MPa）条件下，主裂缝的长度显著增加，裂缝显示从模型外侧向外扩展的趋势明显。这表明较低的注射压力更有利于主裂缝的延伸，同时抑制了次级裂缝的发展。在高注射压力（例如66 MPa）条件下，裂缝主要集中在注入孔附近，形成复杂且密集的短裂缝网络，分支裂缝数量显著增加。这是因为高注射压力可以提供更强的驱动力，促使流体在层理界面或薄弱区域诱导更多次级裂缝，从而增加裂缝网络的复杂度。无论注射压力水平如何，裂缝传播都受到层理的显著影响：当层理倾角较小时（例如0°、15°），裂缝在遇到层理后倾向于沿层理方向延伸；当层理倾角较大（例如45°、75°、90°）时，一些裂缝会穿透层理并继续延伸，而其他裂缝仍沿层理方向延伸。这反映了层理界面和原地应力共同作用下的裂缝偏转特性。注射压力通过改变流体驱动力来调节主裂缝延伸与次级裂缝发展之间的平衡。低压促进主裂缝的延伸，而高压促进形成复杂的短裂缝网络，这一过程始终受到层理结构的限制。

4.4. 层理密度影响规律分析
吉阳凹陷中的层理密度大约为每100毫米4-10条。因此，设置了两种模拟情景：低层理密度（4条/100毫米）和高层理密度（9条/100毫米）。为了分析层理对页岩裂缝传播的影响，层理密度分别设定为4条/100毫米、5条/100毫米、9条/100毫米和10条/100毫米。如图12所示，比较了高层理密度和低层理密度下不同层理倾角的启动压力。从图12可以看出，低层理密度下的启动压力高于高层理密度下的启动压力。低层理密度（4条/100毫米）下的启动压力约为59.7 MPa，而高层理密度（10条/100毫米）下的启动压力为58.5 MPa。层理密度可以降低页岩强度，改变其力学性质，从而降低裂缝启动压力。当层理密度保持不变时，增加或减少层理倾角对启动压力的影响较小，启动压力几乎不变。在高密度层理的模拟中，15°和45°层理倾角之间的启动压力差异较大。
分别设置层理密度为4条/100毫米、5条/100毫米、9条/100毫米和10条/100毫米，以研究不同层理倾角下的裂缝传播模式，如图13所示。从图13可以看出，在高层理密度（例如9条/100毫米、10条/100毫米）条件下，裂缝更有可能穿透层理界面并向最大水平主应力方向延伸。例如，当层理倾角为45°、75°和90°时，水力裂缝会穿透层理并继续延伸，受层理方向的约束较小。在低层理密度（例如4条/100毫米、5条/100毫米）条件下，裂缝围绕注入孔传播，并倾向于沿层理方向延伸，层理界面的引导作用更为明显。层理密度对裂缝路径的影响与层理倾角有协同作用。当层理倾角较小时（0°、15°），裂缝容易受层理方向的影响，但即使在高层理密度下，一些裂缝仍会穿透层理。当层理倾角较大（45°、75°、90°）时，高层理密度显著促进裂缝穿透层理；而在低层理密度下，裂缝沿层理方向延伸的趋势更为明显。在低层理密度下，裂缝沿层理的方向性传播更为突出，裂缝网络具有很强的方向性。与此同时，在高层理密度下，穿透层理的裂缝比例增加，结合原位应力的驱动力，裂缝网络倾向于更多地向最大主应力方向集中，其复杂性由层理界面和原位应力的共同作用调节。层理密度通过改变层理界面的分布密度来调节“沿层理传播”和“穿过层理传播”之间的平衡。高层理密度增强了裂缝的穿透能力，而低层理密度增强了裂缝沿层理的导向性，这一过程受到层理倾角和原位应力的协同约束。

讨论
本研究阐明了原位应力、注入参数和层理特性对裂缝传播的影响规律。这些发现为页岩油储层的水力压裂工程设计提供了关键的理论基础。在实际操作中，可以根据地层中的原位应力状态和层理发育特性优化注入参数，从而有效控制裂缝传播的方向和形态，构建更理想的裂缝网络，并提高页岩油储层的采收效率[48]。这项研究在离散元模拟方面取得了一定的成果，但仍存在局限性。未来可以考虑建立三维模型，纳入更多实际因素（如岩石非均质性、天然裂缝等），并进行多因素耦合的定量研究，以进一步深入理解层理发育页岩中的裂缝传播规律，从而为页岩油储层的高效开发提供更准确的理论支持[49,50]。

层理是中国东部断陷盆地（以济阳凹陷为代表）页岩的核心结构特征，其密度范围为每100毫米4到10条。这种层状结构对水力裂缝的启动、传播和网络形成具有重要的调控作用，是决定页岩油储层刺激效果的关键因素。目前，许多研究采用了离散元方法（DEM，包括PFC 2D 6.0）来研究页岩中的水力裂缝传播，主要关注原位应力差异、注入速率和天然裂缝的影响。然而，很少有研究系统地探讨层理特性（倾角）与施工参数（注入速率、压力）和地质条件（原位应力）对裂缝行为的协同效应。这一研究空白限制了对层理发育页岩中裂缝形态的准确预测和压裂设计的优化。

为了解决这一局限性，本研究首先建立了层理发育页岩的PFC2D 6.0数值模型，并通过将模拟结果与实验应力-应变曲线进行比较来验证其合理性，然后进行了控制变量模拟以明确多种因素对裂缝传播的耦合机制。与之前仅孤立分析单一因素（例如，仅原位应力或层理角度）的离散元研究不同，本研究量化了层理特征与其他变量的交互效应。例如，在高层理密度（9-10条/100毫米）下，裂缝更容易穿透层理界面，并沿最大水平主应力方向扩展，即使在大倾角（45°、75°、90°）下也是如此。相比之下，低层理密度（4-5条/100毫米）增强了层理的导向作用，使裂缝在注水孔附近沿层理方向传播，而不论应力差异如何变化。此外，本研究揭示了注入速率与裂缝复杂性之间的非单调关系——裂缝网络在注入速率为0.03米/秒时最为复杂，而不是像之前研究中简化的那样随着速率的增加而线性增加。这些发现填补了使用离散元方法研究层理发育页岩多因素耦合的空白，并为压裂工程提供了更精准的指导。

层理是中国东部断陷盆地（以济阳凹陷为代表）页岩的核心结构特征，其密度范围为每100毫米4到10条。这种层状结构对水力裂缝的启动、传播和网络形成具有重要的调控作用，是决定页岩油储层刺激效果的关键因素。目前，许多研究采用了离散元方法（DEM，包括PFC 2D 6.0）来研究页岩中的水力裂缝传播，主要关注原位应力差异、注入速率和天然裂缝的影响。然而，很少有研究系统地探讨层理特性（倾角）与施工参数（注入速率、压力）和地质条件（原位应力）对裂缝行为的协同效应。这一研究空白限制了对层理发育页岩中裂缝形态的准确预测和压裂设计的优化。

为了克服这一限制，本研究首先建立了层理发育页岩的PFC2D 6.0数值模型，并通过将模拟结果与实验应力-应变曲线进行比较来验证其合理性，然后进行了控制变量模拟以明确多种因素对裂缝传播的耦合机制。与之前仅分析单一因素的离散元研究不同，本研究量化了层理特征与其他变量的交互效应。例如，在高层理密度（9-10条/100毫米）下，裂缝更容易穿透层理界面，并沿最大水平主应力方向扩展，即使在较大倾角（45°、75°、90°）下也是如此。相比之下，低层理密度（4-5条/100毫米）增强了层理的导向作用，使裂缝在注水孔附近沿层理方向传播，而不受应力差异的影响。此外，本研究还发现注入速率与裂缝复杂性之间存在非单调关系——裂缝网络在注入速率为0.03米/秒时最为复杂，而不是像之前研究中简化的那样随着注入速率的增加而线性增加。这些发现填补了使用离散元方法研究层理发育页岩多因素耦合的空白，并为压裂工程提供了更精准的指导。

分析还确定了两种典型层理密度情景（低：4-5条/100毫米；高：9-10条/100毫米）在裂缝路径中的不同调控作用。低层理密度增强了“层理导向效应”：富粘土的层理界面作为薄弱平面，使裂缝沿层理偏转并形成定向网络。高层理密度增强了“应力驱动的穿透效应”：密集的层理界面增加了裂缝与层理相互作用的可能性，在原位应力的驱动下，裂缝穿透更多层理并集中在最大主应力方向。这种区分补充了之前仅基于“薄弱平面理论”的层理页岩研究，并为裂缝传播模式提供了更全面的机理解释。

由于时间和模型的限制，本研究仍有一些方面需要进一步改进。首先，当前的PFC2D模型将页岩简化为二维系统，无法完全反映层理的三维空间分布以及显著影响水平井实际压裂效果的各种裂缝因素的垂直传播行为。未来的工作将构建三维离散元模型，以模拟裂缝的垂直扩展范围及其与多层理的相互作用。其次，本研究主要关注裂缝传播规律的定性和半定量分析，但缺乏一个综合考虑层理密度、倾角和原位应力的定量本构模型来预测裂缝启动压力和网络复杂性。建立这样的模型对于将研究成果转化为实际压裂设计参数至关重要。最后，模型未考虑济阳凹陷页岩中广泛存在的天然裂缝。将天然裂缝纳入数值模型并探索其与层理的耦合效应将进一步提高裂缝传播预测的准确性。

由于研究范围和模型设置的限制，一些工作仍需进一步推进。本研究阐明了层理发育页岩中水力裂缝传播的变化规律及其影响因素（特别是层理密度和层理倾角的耦合效应），并提供了这些变化规律的理论解释。然而，尚缺乏宏观物理实验的验证，仍需进一步探索数值模拟结果与实际储层裂缝传播特征之间的一致性。未来将在层理发育页岩岩心中进行实验室水力压裂物理实验，并结合CT扫描和声发射监测技术，以验证数值模拟模型的可靠性。尽管本研究对单因素和多因素耦合作用下的裂缝传播规律进行了定性和定量分析，但尚未建立考虑层理特性的裂缝传播预测模型。下一步的研究方向将是建立一个综合考虑层理密度、层理倾角、原位力和注入参数的层理发育页岩水力压裂传播模型，这对于提高济阳凹陷页岩油储层的水力压裂设计精度至关重要。

结论
本文主要利用离散元方法分析了页岩在载荷下的力学性质和裂缝传播规律。主要结论如下：
(1) 随着最小水平主应力的增加，启动压力显著上升。水平主应力差异主导了裂缝方向，较大的应力差异增强了向最大主应力的偏转。小倾角（0°、15°）使裂缝沿层理方向传播，而大倾角（45°、75°、90°）则受应力差异驱动，较大的应力差异增加了穿透比例。
(2) 注入速率与启动压力呈负相关。当注入速率从0.01米/秒增加到0.09米/秒时，启动压力从59.34兆帕降低到58.88兆帕，尽管减少幅度在后一阶段变缓。注入压力对裂缝形态和传播模式有显著影响：低注入压力（例如60兆帕）有利于主裂缝的扩展，而高注入压力（例如66兆帕）促进了注水孔附近形成复杂密集的短裂缝网络。
(3) 层理密度和倾角协同控制裂缝路径。高层理密度（9-10条/100毫米）增强了层理穿透能力，有利于向最大水平主应力方向扩展；低层理密度（4-5条/100毫米）增强了层理导向作用，使裂缝在注水孔附近沿层理传播。小倾角增强了层理导向作用；大倾角下，高密度促进穿透，低密度则增强沿层理的扩展。