水沙混合物裂隙渗流特性分析

刘 玉，韩 雨，张 强，李 猛，王志飞

(1.江苏师范大学 机电工程学院，江苏 徐州 221116; 2.中国矿业大学 煤炭资源与安全开采重点实验室，江苏 徐州 221116; 3.中国矿业大学 深部岩土力学与地下工程国家重点实验室，江苏 徐州 221116)

我国煤炭资源主要集中在陕甘蒙等西部地区，其地质条件具有厚风积沙、基岩薄的特点，加之煤层埋深较浅极易发生突水溃沙事故。针对溃沙机理研究，国内外学者进行了大量工作，取得了一些成果[1-3]。李东等[4]研究了含水层沉积特征，初步查明了典型顶板水害机理，并提出疏放水和注浆改造的突水溃沙防治措施。浦海等[5]通过格子 Boltzmann 方法分析裂隙水沙两相流，分析了颗粒粒径和裂隙宽度对溃沙速度的影响，借助单裂隙的研究成果，分析溃沙的演化过程。杜锋等[6]对突水溃沙进行归类，并通过自制的试验仪器获得了孔隙度对含水层溃沙的影响，以及沙粒径对溃沙的影响。隋旺华等[7]通过试验获得土的黏性与水利坡度的关系，以及土颗粒与成分对水沙突变的影响，并得到初始水头及突水口是控制溃沙的关键因素。梁艳坤等[8]利用分形理论及离散元方法建立了垮落带破碎岩体溃沙数值模型，通过分析得知较大连通空隙形成了溃沙通道，溃沙在通道中流速分布呈纺锤形分布。范立民等[9]在GIS平台下构建影响突水溃沙的关键因素模型的评价模型，并在具体矿区进行验证，但此模型不具有通用性。杨鑫等[10]利用自制仪器进行沙粒启动试验和溃沙试验，得到了溃沙临界速度;溃沙速率与初始水力梯度呈指数关系，溃沙量与水力梯度呈线性关系。许延春等[11]建立了楔形保水压采动溃沙地质模型并通过试验进行验证，提出了煤柱留设方法以实现含水层下防沙。针对西部煤矿防止突水。孙强[12]等建立了风积沙、粉煤灰及硅酸盐水泥不同配比的二维物理模型，研究采动裂隙发育规律，并分析隔水层的稳定性。张凯等[13]通过试验获得渗流速度和沙颗粒大小呈线性分布，并得到孔隙大、渗流速度快、颗粒粒径分布大时易发生沙涌事故。刘玉和李顺才[14]从裂隙渗流的角度研究突水溃沙机理，获得了水沙在不同粒径下的渗流特性。张贵民等[15]研究不良钻孔产生突水溃沙的机理，提出关键影响因素是含水层厚度和钻孔直径，并提出预防突水的措施。张俊霞等[16]针对富含水沙的第四纪岩层进行研究，通过冻结孔内打井注浆等方式进行井筒施工，成功进行了突水溃沙控制。张金才和彭苏平[17]在现场实测的基础上，对浅埋煤层的采动煤层破坏进行研究，提出了防止溃沙煤柱的优化设计方法，以获得最大的开采效率及环境保护。范钢伟等[18]利用采前脱水和灌浆等方式，通过在沟渠周围打井的13口地面抽水井和33口地下垫层降低了地下水位。位于沟底的含水层在含水层降水方面普遍优于河岸的含水层。通过25个表面孔向埋地沟底注入化学灌浆以防止突水。张改玲[19]研究了化学注浆前后的粗沙试验的导水性试验，结果表明化学灌浆明显降低了沙土的渗透性，化学灌浆沙的导电性一般随围压的增大而降低。

笔者在前人研究的基础上，针对西部浅埋煤层垮落带突水溃沙问题，以裂隙水沙渗透特性为研究对象，采用室内试验的方法获得水沙裂隙渗流的滞后性变化特性，进一步利用数值模拟的方法研究水沙渗流场的影响因素，对掌握突水溃沙机理具有重要意义。

1 水沙裂隙渗流方案

1.1 测试原理及系统

图1是裂隙试验系统的试验系统图，该系统由渗流仪、搅拌系统、泵送系统和采集系统构成。

图1 裂隙注水系统示意Fig.1 Injection water system of hydraulic extrusion measure

采用试验裂隙为平行粗糙裂隙，粗糙度经测定为JRC 2～4,如图2所示，长度 125 mm，宽度75 mm，面积为9 375 mm2。

图2 试样及裂隙面视图Fig.2 Specimens and fracture surface view

1.2 测试理论及步骤图

为分析裂隙沙粒径、沙质量浓度对裂隙水沙混合物渗透特性的影响规律，制定如下试验方案，方案流程如图3所示。

图3 裂隙渗透试验流程Fig.3 Test flow chart of seepage in fractured rock

具体试验步骤如下:

(1)试验系统调试:装配系统并调试，观察有无漏水现象。启动电脑及记录仪，检测数据采集情况，保证试验系统正常采集。

(2)装料:把粗糙裂隙岩样装入渗流仪，保证岩样上下对齐;加入沙粒径0.092～0.138 mm，使得水沙混合物质量分数为 20%，开动搅拌系统使之均匀混合。

(3)渗透试验:设定裂隙开度，改变电机转速，使得稳定渗透压力差在0～1.2 MPa，记录下稳态时压力差和渗流速度。重复上述试验;水沙混合物的体积质量40,60,80 kg/m3，重复上述试验;改变沙粒径0.138～0.184，0.184～0.230，0.230～0.276 mm，重复上述试验。

(4)卸载:卸料关闭螺杆输送泵，关闭阀门，一次卸下渗透仪等工具。

2 试验结果及讨论

2.1 滞环曲线的种类及转化

试验所用沙子取于河北省灵寿县，有4种粒径:0.092～0.138，0.138～0.184，0.184～0.230，0.230～0.276 mm。试验时裂隙开度设定为0.5 mm，裂隙开度与平均粒径比为2.7∶1。在沙体积质量20,40,60,80 kg/m3分别进行渗透试验。依据上述的测试原理和方法，按步骤进行渗流试验，试验发现在压力梯度上升和下降2个阶段，压力梯度与渗流速度不成一一对应关系，随着压力梯度的增加和减少，压力梯度与渗流速度在平面坐标中形成一条封闭的曲线。在压力梯度增加和下降的过程，获得4种滞环曲线，如图4中的I～IV四种类型，其中，Gp为最大滞后量，v为渗流速度。图4中完整滞回曲线分为升程段OAB和回程段BA′O两段;I型曲线升程曲线和回程曲线相交，II型曲线与磁滞回曲线类似，III型曲线升程曲线和回程曲线部分重合，IV型曲线升程曲线和回程曲线全程没有重合点。

随沙粒径和浓度增大，这种由压力梯度-渗流速度构成的滞回曲线逐渐由I型向IV型转变，见表1。

由表1可以看出，当裂隙开度b=0.5 mm时，随沙粒径和密度增大，Ⅰ型、II型、IV型、IV型滞环曲线可以由前者向后面的类型转化。原因在于，裂隙中水沙流的运动形态有多种，单相流(沙与水之间没有相对速度)、两相流(沙与水之间存在相对速度)、段塞流(沙与水之间存在一个或多个分界面)等。随沙粒径及浓度等发生变化，运动形态发生变化，导致滞环曲线发生改变。

2.2 沙粒径对滞后性参量的影响

为分析沙粒径对Ⅳ型滞环曲线滞后性参量的影响，将4组试样结果中的渗流速度-压力梯度滞环曲线最大滞后量Gp和滞环面积S汇总于表2。

由表2可知，在裂隙开度b=0.5 mm、沙体积质量ρs=80 kg/m3时，随着沙粒径ds的增大，Ⅳ型水沙渗流速度-压力梯度滞环曲线的最大滞后量Gp有逐渐增大的趋势，并且当沙粒径ds处于0.115～0.161 mm和0.207～0.253 mm时，滞环曲线的最大滞后量Gp增加的比较缓慢;当水沙混合物中沙粒径ds处于0.161～0.207 mm时，滞环曲线的最大滞后量Gp增加的比较快，近似于线性增长。

图4 4种类型滞回曲线Fig.4 Four types of hysteretic curves

类别粒径/mmρs/(kg·m-3)204060800.092～0.138ⅠⅡⅢⅣ滞环曲线类别0.138～0.184ⅠⅢⅢⅣ0.184～0.230ⅡⅣⅣⅣ0.230～0.276ⅣⅣⅣⅣ

表2 不同沙粒径下Ⅳ型滞环曲线的滞后性指标(b=0.5 mm，ρs=80 kg/m3)Table 2 Hysteresis indicators of Ⅳ type hysteresis curves under different sediment grain size at b=0.5 mm,ρs=80 kg/m3

2.3 沙质量浓度对滞后性参量的影响

设定裂隙开度b=0.5 mm，沙粒径ds=0.230～0.276 mm，选用沙体积质量ρs为变量，通过设定20，40，60和80 kg/m3的4种不同水沙混合物,研究其裂隙渗流的滞环曲线。

试验中沙体积质量为ρs=80 kg/m3，不同沙粒径下的滞环曲线如图5所示，随沙体积质量ρs增大，Ⅳ型水沙渗流速度-压力梯度滞环曲线的最大滞后量Gp逐渐增大。在沙体积质量为20～40 kg/m3和60～80 kg/m3时，滞环曲线的最大滞后量Gp增加较快，接近于线性增长，但滞环面积S则增长缓慢。当水沙混合物中沙体积质量ρs处于40～60 kg/m3时，滞环曲线的滞环面积S增加的较快，近似于线性增长，而最大滞后量Gp增长缓慢。

图5 不同沙质量浓度下Ⅳ型滞环曲线特性(b=0.5 mm，ds=0.230～0.276 mm)Fig.5 Under different sand mass concentration type Ⅳ hys-teresis curve features at b=0.5 mm,ds=0.230～0.276 mm

3 水沙裂隙渗流场及其影响因素模拟

3.1 水沙裂隙渗流场变化

采用ANSYS Fluent软件模拟粗糙裂隙中水沙两相流动，沙粒密度ρs=2 650 kg/m3，沙粒径Dp=0.04 mm，沙粒体积分数φ=4.06%，裂隙入口速度分别为0.349，0.532，0.697和0.869 m/s。

由图6可以看出，数值模拟得到的裂隙水沙流压力梯度绝对值-渗流速度曲线与试验得到的曲线变化趋势基本一致，压力梯度绝对值与渗流速度之间呈非线性关系。数值模拟结果小于试验结果，且数值模拟结果与试验结果绝对误差接近，相对误差则随流速增大而减小。

图6 压力梯度绝对值-渗流速度曲线比较Fig.6 Pressure gradient absolute value-seepage velocity curve comparison

图7给出了裂隙入口速度为0.869 m/s时，裂隙入口处平均压力随时间变化曲线。其中，t=0～0.12 s区间为连续相流场压力曲线，t=0.12～0.20 s区间为注入沙粒后的水沙两相流场压力曲线。可以发现，裂隙入口压力在20 kPa附近剧烈波动，但是没有明显的衰减趋势，处于动态的稳定状态。

图8给出了裂隙入口速度为0.869 m/s时，t=0.17，0.22和0.27 s 三个时刻，X1=1.4 mm截面上速度、湍动能和压力分布。

由图8(a)可以看出，在X1≥2 mm内，流速在2.7 m/s附近波动，峰值和谷值相差很小，但波峰波谷并不吻合。在X1≥10 mm段，流速陡增。流体速度重新组合，沿X1流动方向各过流断面速度分布不断变化，壁面处黏滞作用使流体减速，边界层外中间部分流体加速运动。由图8(b)可以看出，在X1≥20 mm内，湍动能k在0.03～0.07 m2/s2剧烈波动。在X1≤ 2 mm段，湍动能陡增，说明从入口处开始的流体速度重新组合，造成流场湍流强度不断增强。由图8(c)可以看出，压力虽然在局部有些起伏，但是总趋势随X1直线下降。

3.2 沙粒径对渗流场影响的影响

下面讨论沙粒密度ρp=2 650 kg/m3、沙粒体积分数φ=4.06%条件下，沙粒径对裂隙流场的影响。图9给出了压力梯度绝对值与沙粒径Dp曲线，水沙在粗糙裂隙中流动，在沙粒径为0 mm时压力梯度绝对值最大，在沙粒径为0.12 mm时压力梯度绝对值最小，在沙粒径为0.01 mm时压力梯度绝对值出现极小值。当 0.04 mm时压力梯度绝对值随沙粒径增大而单调递减。这说明水沙在粗糙裂隙流动中流动，当沙粒径较小时，压力损失随沙粒径增大而升高;当沙粒径较大时，压力损失随沙粒径增大而降低。

图9 粗糙裂隙横截面上连续相流体时均速度分布Fig.9 Uniform velocity distribution of continuous phase fluid on cross section of rough fracture

文献[17]讨论孔隙介质中湿相和非湿相流体渗透的毛细滞后效应，岩石裂隙中水沙流动行为虽然复杂，但固液接触角的变化是滞后现象的一种原因。

3.3 沙体积浓度对渗流场影响的影响

图10给出了压力梯度绝对值-沙粒体积分数曲线，在相同的沙粒密度、相同的沙粒径和相同的沙体积分数条件下，水沙在粗糙裂隙中流动的压力梯度绝对值比在光滑裂隙中的大近40倍，并且压力梯度绝对值随沙粒体积分数的变化趋势不同，说明沙粒体积分数对压力梯度的影响取决于裂隙表面形态。

图10 粗糙裂隙横截面上连续相流体时均速度分布Fig.10 Uniform velocity distribution of continuous phase fluid on cross section of rough fracture

由图10可知，水沙在粗糙裂隙中流动，裂隙横截面上连续相流体时均速度分布受沙粒体积分数影响很大。在X1=50 mm截面上，时均速度呈不对称的M形分布，有多个极值点。每个位置流体质点时均速度受沙粒体积分数影响都非常显著，并且最大时均速度位置在X1=1.2～1.5 mm变动。在X1=50.5 mm截面上，时均速度分布不再左右对称但峰值仍在中线X1=1.4 mm附近，有多个极值点;在壁面附近时均速度随沙粒体积分数剧烈变化;在其它位置，时均速度也随 变化但没有壁面附近显著。这表明，沙粒体积分数对X1=50.5 mm截面上水沙两相流动的影响主要发生在壁面附近(边界层)。

4 结 论

(1)随沙粒径和沙质量分数的增大，水沙渗流速度-压力梯度滞环曲线的最大滞后量呈类似线性增大。

(2)粗糙裂隙裂隙面的构造导致水沙渗流场具有一定随机性，原因在于水沙在粗糙裂隙中沉积，使得粗糙裂隙面、水和沙之间相互关系复杂，压力梯度-水沙渗流速度关系不惟一。

(3)而在粗糙裂隙中流动时，当沙粒径较小时，压力损失随沙粒径增大而升高;当沙粒径较大时，压力损失随沙粒径增大而降低;裂隙横截面上流体时均速度和湍动能分布受沙粒径影响很大，表现为极值点位置的偏移。

(4)而在粗糙裂隙中流动时，沙粒的存在减小了压力损失，在沙体积分数为1.02%时压力梯度绝对值最小，当沙体积分数φ≥2.07%时压力梯度绝对值变化幅度很小;裂隙横截面上流体时均速度和湍动能分布受沙粒体积分数影响很大，表现为极值点位置的偏移。