单裂隙辐射流节理面破坏过程水-力耦合研究

2022-06-21杜万军柴军瑞许增光

水力发电 2022年1期

杜万军，柴军瑞，许增光，曹成

(1.杨凌职业技术学院水利工程学院，陕西杨凌 712100；2.西安理工大学省部共建西北旱区生态水利国家重点实验室，陕西西安 710048)

0 引言

在地下水利工程中，裂隙岩体在渗透力和外部压力共同作用下的稳定性和流体运动规律均比较复杂，且存在相互耦合关系。通过不同节理模型试验[1-2]分析粗糙程度和外力对裂隙水流的影响规律发现，直剪初期，裂隙渗透性发生短暂降低后继续攀升直到平缓[3- 4]，水压越高，裂隙渗透性越强[5]。剪应力在不同位移条件下出现波峰和波谷最后趋于稳定。裂隙开度的变化会引起渗流量的改变，尤其是裂隙充填物的颗粒组成和空隙率的空间分布对水流渗流通道的影响较大[6-7]。辐射流的随机性和复杂性更显著，更符合实际工况，由于裂隙主要渗流通道和剪切损伤后充填物对水流的影响显著，所以明确渗流通道和合理简化充填物的影响效应至关重要。为此，本文从宏观试验的裂隙剪切受力和渗透性出发并结合数值模拟，分析了粗糙单裂隙辐射流的剪切渗流耦合特性，可为岩体复杂节理裂隙水力特性的研究提供参考。

1 直剪渗流试验与耦合过程分析

1.1 试验设计

岩石节理直剪渗流耦合系统见图1。试验系统[8]可对试件施加垂直压力、剪切力、水压力。垂直压力和剪切力由油压系统提供，水压力由液压氮气瓶和水箱提供，试验试件采用β型高强度石膏加水拌合浇筑在固定的模具中预制而成。试件为高度80 mm、直径200 mm的圆柱体(见图2)。各项力学指标：密度为1.78 kg/m3，抗压强度为47.83 MPa，弹性模量为2.58 GPa，泊松比为0.31。

图1 岩石节理直剪渗流耦合系统

图2 粗糙裂隙试样(单位：mm)

根据不同法向压力和入口水压力设置2种试验工况，每种工况各有3种子工况，共有子工况6种。试验工况见表1。

表1 试验工况

1.2 剪切渗流耦合过程

裂隙节理剪切破坏过程见图3。从图3可知，裂隙节理的破坏主要分为2个阶段：剪切接触未破坏阶段和剪切破坏滑移阶段。发生剪切初期，节理凸起未完全接触，裂隙机械隙宽较均匀。当节理相互错动并局部发生接触时，顺剪切方向机械隙宽逐渐增大，远离端逐渐减小。当节理剪切位移继续增大，凸起尖部薄弱部位首先发生损伤破坏，节理之间继续相互错动并发生壁面滑移，大部分充填物累积在节理凹槽处，剩余部分被挤压在节理面接触位置，阻碍了上部试件节理面的法向位移。当节理面完全破坏时，顺节理方向渗流通道最大，水流几乎不沿垂直节理方向流动，故此时水-力耦合作用主要出现在损伤后的壁面接触和优势水力通道区域。

图3 裂隙节理剪切破坏过程

2 试验结果分析

2.1 裂隙剪应力

裂隙节理面受剪切作用时，剪应力出现连续波峰与波谷，首次波峰均表现出法向压力越大，剪应力越大，入口水压越大，剪应力越小的总体趋势，后期由于产生破碎颗粒充填物，导致剪应力变化规律不统一，但剪应力峰值均出现在节理咬合阶段。各工况剪应力与剪切位移的关系见图4。从图4可知：

图4 各工况剪应力与剪切位移关系

(1)节理初期咬合开始，剪应力快速爬升，最大值分别为1.65、1.98、2.93 MPa，此时节理面开始发生首次破坏并产生大量碎屑，裂隙面碾压碎屑充填物并发生滑移摩擦，导致碎屑颗粒粒径减小，出现剪涨现象[9]，细小颗粒被水流沿节理和充填物形成的渗流通道冲离裂隙面，剪应力出现降幅。破坏后节理面重新咬合开始，剪应力又迅速增高，最大值分别为3.61、3.28、4.08 MPa。节理面出现2次破坏，随后剪应力快速降低并继续在水-力共同作用下发生滑移摩擦，且第2次剪应力峰值均最高，说明初期破坏的节理面抗剪能力较弱，第2次咬合和充填物的生成持续增大了裂隙面的抗剪强度，从而出现最大峰值。随着节理面破坏程度逐步提升，剪应力产生的峰值趋于减小。

(2)首次峰值最大剪应力依次为3.14、1.98、0.80 MPa，第2峰值依然为整个剪切过程最大值，但波谷值变化不大，主要出现在摩擦滑移节理咬合前阶段，说明碎屑充填物的产生在一定程度上降低了节理咬合深度，减小了剪切破坏范围，水流与多孔介质破损颗粒充填物共同组成了缓冲区域，导致上下裂隙面之间发生滑移，使得壁面之间的抗摩擦强度较稳定。

2.2 裂隙渗透性

渗流量和等效渗透系数可以反映剪切破坏过程中的裂隙宏观渗透特性。各工况裂隙渗流量与剪切位移的关系见图5，各工况等效渗透系数与剪切位移的关系见图6。

图5 各工况裂隙渗流量与剪切位移关系

图6 各工况等效渗透系数与剪切位移关系

(1)从图5a和图6a可知，法向压力越大，裂隙渗流量和等效渗透系数越小，并随剪切位移的增大出现多个峰值。渗流量最大值分别为57.5、48.5、37 cm3/s，等效渗透系数最大值分别为0.26、0.22、0.17 m/s，首次峰值出现在剪切位移10mm左右处，恰好与剪应力最小值位置对应，说明此时裂隙渗流能力较强。随着剪切位移的增大，充填物颗粒继续受力挤压、密度增大，节理开始2次咬合，并把充填物挤压在节理咬合处，导致渗流量和等效渗透系数出现第2次峰值，且与首次峰值大小基本一致。说明剪切位移为20 mm左右(剪应力最小值)时，继续出现机械隙宽的突增，导致裂隙渗透性快速升高。

(2)从图5b和图6b可知，不同入口水压下，等效渗透系数的连续波动性不大。水压力越大，渗流量和等效渗透系数越大，渗流量最大值依次为65.5、48.5、44.5 cm3/s，等效渗透系数最大值依次为0.37、0.32、0.23 m/s。由于剪切力相同条件下，节理面破坏程度趋势差异不大，裂隙最大渗透性出现在12～15 mm剪切位移处，说明随着水压力增大，水流带走了粒径相对大一点的破碎颗粒，裂隙机械隙宽的增大导致裂隙渗透能力出现较大提升。当剪切位移为25～35 mm时，渗流量依次保持在32.3、37.5、39.2 cm3/s，等效渗透系数依次保持在0.063、0.099、0.115 m/s左右，两者均相对稳定。主要由于裂隙面保留的大粒径颗粒依然充满节理面且形成了由充填物组成的水流通道，并与增大的水压力共同缓解了节理之间相互咬合和破坏作用，此时上下节理面主要处在轻度咬合和剪切滑移阶段，导致裂隙渗透性较稳定。

3 数值模拟分析与讨论

3.1 数值模拟理论依据

基于平均雷诺RANS方程，由于流体的惯性力比黏性力大很多，故计算模型选择RANSk-ε方程[9]，即

ρ(u·∇)u=∇·[-ρI+(μ+μT)(∇u+(∇u)T)]+F

(1)

ρ∇·(u)=0

(2)

(3)

(4)

(5)

Pk=μT[∇u:(∇u+(∇u)T)]

(6)

式中，ρ为流体密度；u为流速张量；I为单位张量；μ为粘性系数；μT为紊流粘性系数；F为体力张量；k为紊流动能；σk为湍流动能的普朗特数；Pk为平均速度梯度产生的湍流能量；ε为吸入率；σε为耗散率的湍流普朗特数；Cε1、Cε2和Cμ为计算出的系数，Cε1=1.44、Cε2=1.92、Cμ=0.09。

3.2 裂隙优势水力通道数值模拟

根据试件尺寸和节理几何特征建立三维数值计算模型，见图7。图7中，1、2、3、4为特征截线位置。数值计算采用进口和出口恒定压力边界条件，入口边界设置为恒定水压，大小根据试验工况计算得到，出口边界条件为充分自由发展流动，外部压力设为0。裂隙壁面设置为滑移边界条件，水流从下部试件中心注入，自辐射中心由内向外自由流动，裂隙面处在XZ平面，下部注水方向沿Y轴。模型裂隙宽度为0.25 mm，入口水压分别为表1工况2中的1.5×104、2.3×104、3.6 ×104Pa。由于裂隙充填物的随机性和复杂程度无法直接建模计算，而试验阶段已对充填物做了宏观分析，故数值模拟过程不考虑充填物力学特性的影响，主要以不同水压下节理破坏前后裂隙形态和开度变化对水力特性的影响为主，适应于节理破坏程度不严重的裂隙。

图7 数值模拟模型

未破坏节理不同水压下水流速度场云图见图8。从图8可知，裂隙未破坏前节理面完整，水流沿距离辐射中心近处的节理平行方向凹槽辐射流出，形成优势水力通道。优势水力通道范围内，越靠近辐射中心位置，流速分布趋势主要表现为中部流速高，依次径向减小，且初期对称性不强。优势水力通道范围之外的裂隙水流速度较小甚至无流动，垂直节理方向的远端几乎无水流通过。随着入口水压的增大，优势水力通道范围略有扩大，呈基本中心对称的辐射状“双U”形分布，裂隙开度越大，形成的优势水力通道的范围相对越小，越集中于顺节理方向。主要由于机械隙宽的增大提升了渗流通道，并束窄了优势水力通道扩散范围。随着节理面破坏程度的加剧，凸起的节理相互错动并产生滑移摩擦，导致水流“跨越”凸起“障碍”后沿顺节理凹槽流出，优势水力通道的范围开始沿垂直节理径向小幅扩散，但范围有限。当优势水力通道积累的充填物增多且裂隙挤压导致其密实时，水流绕过固体颗粒继续寻找新通道，从而出现局部优势水力通道范围扩大的现象。

图8 未破坏节理不同水压下水流速度场云图

破坏节理不同水压下水流速度场云图见图9。从图9可知，节理破坏后，由于壁面之间机械隙宽突增，挤压性更强，导致优势水力通道的范围略有减小，流速梯度有所提高。说明节理破坏虽然引起了水力通道横向辐射范围减小，但流速整体有所升高。

图9 破坏节理不同水压下水流速度场云图

为了解节理垂直方向轴线位置水流分布，选择图7b中3个特征位置，按照隙宽0.25mm和入口水压2.3×104Pa计算，得到特征位置水流速度分布，见图10。从图10可知，流速分布呈现典型的漏斗状，辐射中心位置由于水流法向高速流入时受到上壁面的正面阻挡，流速反而较低(水压最大)。靠近辐射中心的位置1，在节理破坏前和破坏后的最小流速分别为0.5、0.52 m/s，差异不大。但未破坏节理的最大流速出现在靠近辐射中心附近，破坏后的节理最大流速出现在末端，说明节理破坏在一定程度上释放了对水流的约束作用，更容易越过近处节理发生流动，优势水力通道范围虽小，但流速整体较高。由于节理破坏引起优势水力通道范围束窄，导致位置2、3处最大流速较小，而节理未破坏时却能达到1～1.6 m/s。

图10 特征位置水流速度分布

分别计算裂隙面破坏前后图7b中截线4位置距离辐射中心20mm内的流速分布，结果见图11。从图11可知，入口水压从1.5×104Pa增至3.6×104Pa，裂隙流速均表现为递增趋势，但节理破坏后裂隙水流辐射范围比破坏前的5mm提高至15mm左右，且整体流速较高，最大值为1.49 m/s，是破坏前最大流速0.75 m/s的近2倍，说明节理的破坏在某种程度上提升了水流垂直节理方向的辐射范围。主要原因在于节理凸起高度降低，导致垂直节理轴线方向的裂隙开度增大，水流在该方向的辐射力度增强所致。

图11 不同水压截线4位置流速分布

3.3 节理破坏与优势水力通道关系

工况1、2下部试件节理面破坏分别见图12、13。从图12、13可知：

图12 工况1各子工况下部试件节理面破坏

图13 工况2各子工况下部试件节理面破坏

(1)工况1节理面整体破坏程度要高于工况2，说明法向压力越大，节理破坏程度越严重，生成的裂隙充填物也越多。剪切滑移产生的碎屑充填物主要集中于节理面的顺剪切方向的两侧，尤其是颗粒较大的碎屑分布在已破损的节理凹槽处，个别部位在水-力耦合作用下密实度升高，降低了水流的辐射流出强度。靠近中部位置节理相对较完整，成为水流优势路径的首选，充填物沿着主要渗流通道流失，裂隙面在水流浮托力作用下抑制了节理的深度破坏，进而形成了优势水力通道，与数值模拟结果相符。当法向压力增大到一定程度时，节理抗剪强度和水压力无法维持抵御作用，进而导致节理损伤程度进一步加剧，此时水流主要从充填物形成的多孔弹性介质中辐射流出。

(2)入口水压越大，附着在节理面的充填物数量显著减少，且主要分布在优势水力通道范围外，说明水流从优势水力通道流出的过程中即降低了节理的损伤程度又冲走了部分充填物颗粒，而充填物的流失又可能造成节理损伤程度进一步加深。优势水力通道范围之外几乎无水流流动，此区域充填物几乎不发生流失，对节理面起到一定的保护作用。滞留在靠近顺节理中部位置的充填物也会引起优势水力通道的局部扩散，当节理破坏程度极度严重时，优势水力通道扩散的随机性增大，数值模拟结果则很难适用。