榆神矿区地处我国西部干旱与半干旱地区，位于侏罗纪延安组煤层之上的全新统风积沙层和第四系萨拉乌苏组等松散孔隙含水层是地区惟一具有供水意义和重要生态价值的含水层[1]，是矿区保水采煤的目标保护层[2]，由于煤层开采形成的采动裂隙发育至基岩含水层内部甚至近地表的松散层含水层，导致地下水沿采动裂隙间接或直接渗漏失水，一般形成大而稳定的顶板失水问题[3-4]，加剧了矿区水资源短缺现状以及生态环境的负面响应[5-6]。目前相关保水采煤理论与技术研究主要以“覆岩分带理论”为研究基础，主要研究不同采煤条件下导水裂隙带发育高度与覆岩含水层空间关系[7-8]，来指导保水采煤技术实践，李文平[9]、范立民[10]、王启庆[11]等根据矿区主采煤层之上松散含水层、隔水层、基岩层以及采动导水裂隙高度的空间组合关系，对榆神府矿区“保水采煤”进行了工程地质条件、保水采煤技术、失水程度等分区，初步形成了“地下水库”[12]、“分层[13]、条带[14]、充填[15-16]”等保水采煤技术，为矿区开展具体保水采煤技术实践做出了较大贡献。而在采煤扰动下含水层地下水失水机制、失水程度等地下水响应的定量分析方面研究较少。
榆神矿区中深部煤层埋深相对较大(>200 m)，采动覆岩一般为“三带型”(垮落带、裂隙带、弯曲带)损伤变形模式，当导水裂隙带发育至基岩内部或隔水关键层底部，即位于弯曲带内松散水层与导水裂隙带之间余留一定厚度的保护层(如离石组黄土层与保德组红土层)，如图1所示,一方面由于被“导水裂隙带”揭露的含水层2地下水大量漏失，在导水裂隙内地下水被迅速疏干，含水层2以侧向直接失水为主;另一方面，由于弯曲带内含水层1与导水裂隙带内的含水层2间的水力梯度加大，加剧了含水层间的流体在压力梯度的作用下向采空区的垂向流动，出现含水层1沿保护层的垂向渗漏失水，而且弯曲带岩层变形导致保护层的渗透能力发生变异，均导致地下水流动状态发生变化。可见，榆神矿区中深煤层开采易形成典型的“侧向直接与垂向渗漏”复合失水模式。
图1 中深煤层开采含水层失水示意
Fig.1 Schematic diagram of water loss in aquifer under the influence of coal seam mining
如何构建煤层采动覆岩损伤变形与含水层失水分析模型是定量研究含水层失水规律的主要手段，笔者以COMSOL多物理场数值耦合分析为平台，提出中深煤层开采覆岩变形损伤的数值化处理方法，构建了榆神矿区曹家滩煤矿工作面尺度上煤层开采覆岩损伤变形与含水层失水耦合分析模型，定量分析了采煤导致的含水层地下水的失水规律，为我国西部干旱半干旱矿区煤-水资源协调开发提供科学依据。
1 采动覆岩变形损伤数值化处理方法
COMSOL Multiphysics为一款大型的高级数值仿真软件，可通过选择或自定义不同物理场的偏微分方程进行任意组合即可实现多物理场的直接耦合分析，且模型参数、定解条件等均支持自定义参数控制，广泛应用于多物理场耦合科学研究以及工程计算，被誉为“最专业的多物理场全耦合分析软件”[17]，笔者以COMSOL为分析平台，探索采动覆岩损伤变形与含水层失水的数值耦合分析方法。
1.1 采动弯曲变形数值化处理
由于岩体的孔隙率又与其所处的应力状态有关，笔者参考文献[18-19]，首先建立起采动应力σ与孔隙率φ变化之间数学关系:
φ=(φ0-φr)exp(-αφσv)+φr
(1)
式中,φ0为0应力状态时的孔隙率;φr为高压应力状态下的孔隙率的极限值，φr≈0;αφ为孔隙率应力敏感系数，其值参考文献[19]，可取5.0×10-8 Pa-1;σv定义为平均有效应力,Pa，压应力为正，可按式(2)计算:
σv=(σ1+σ2+σ3)/3-αp
(2)
其中，σ1，σ2，σ3为3个主应力(Pa);α(≤1)为Biot经验系数，其值取决于岩层的压缩性;p为水压力，Pa。通过式(1)和式(2)，建立了采掘扰动影响下采动应力状态与孔隙率指标关系。
根据PALMER等提出的煤岩渗透率k与孔隙率φ变化的三次指数关系[20-21]，以自定义函数形式建立起采动孔隙率φ与渗透率k变化之间数学关系:
式中，k为渗透率，m2;k0为0应力状态的渗透率，m2。
再根据渗透系数K和渗透率k的换算关系输入模型:
式中，K为渗透系数，cm/s;γ为流体的容重，N/m3;μ为流体的动力黏度系数，Pa·s，其中地下水的比重γ可以近似为γ=104 N/m3=100 Pa/cm2，水的动力黏度系数μ=10-3 Pa·s。
在COMSOL分析平台中，以上模型参数渗透系数K、渗透率k、孔隙率φ均可以自定义参数的形式参与模型的计算与输出，通过模拟计算采动应力进一步分析采动弯曲带覆岩的渗透能力变化。
对于受采动引起的覆岩弯曲变形，首先利用原剖分网格进行计算，模型通过计算出各剖分节点的位移变形量，利用COMSOL软件平台中大变形几何体自动重新剖分计算模块，重新进行网格剖分，形成变形剖分网格。
1.2 采动导水裂隙数值化处理
模型中针对导水裂隙带可采用Mohr-Coulomb岩体塑性损伤判据准则进行识别[19-20]:
I′1=σ′1+σ′2+σ′3
式中，Mohr-Coulomb剪应力破坏准则的一种变换形式，其中,σ′1,σ′2,σ′3为有效主应力，Pa;I′1为有效应力第1不变量,Pa;J′2为有效应力偏量第2不变量,Pa;α1，K1分别与岩石的内摩擦角φ和黏聚力c有关。
由于采动导水裂隙导水能力极强，分布不规则，空间变化较大，地下水流动状态不符合达西渗流特征，因此难以进行数值刻画与分析。笔者根据含水层与导水裂隙带间的地下水运动转化特征，分析可知当采动导水裂隙带发育至某一含水层下部时，如图1含水层1底部与导水裂隙之间余留有保护层，含水层2被采矿裂隙完整贯穿，矿坑涌水一部分是含水层1中地下水通过保护层垂向渗漏流失，另一部分是含水层2地下水沿采矿裂隙侧向进入采掘空间两个来源构成，形成了矿坑侧向与顶部复合涌水工况。而地下水动力学中任意一点总水头定义为
H=z+p/r
(6)
式中，H为水头值，m;z为位置水头，m;r为地下水容重,1 000 N/m3。
如图1所示,被导水裂隙完全切割的含水层2与大气联通，导水裂隙内孔隙水压力p为大气压(约等于0)，因此该接触带上地下水水头高度H近似等于该处冒裂带的位置标高[22]，未被采矿裂隙切割的含水层1通过保护层越流渗漏至矿坑中，由于弱渗透性保护层底部通过采矿裂隙与采掘空间空气联通，即在保护层底部与采矿裂隙顶部接触带上孔隙水压力p=0。
因此在COMSOL分析平台中，可根据式(6)，通过识别导水裂隙带空间范围，将导水裂隙与正常覆岩接触面自定义成地下水达西渗流的定水压边界，合理避免了导水裂隙带难以数值刻画问题，实现模拟分析采动导水裂隙引起含水层的失水规律研究。
1.3 覆岩变形损伤变形与含水层失水模型构建思路
笔者以COMSOL多物理场耦合分析软件为模拟平台，提出采动覆岩损伤变形与含水层失水模型构建思路如图2所示。
图2 采动覆岩损伤变形与含水层失水分析模型构建思路
Fig.2 Construction method for analysis of deformation damage of overburden and water loss of aquifer
(1)在岩石力学分析模块中，基于岩体弹塑性力学分析方法，通过自定义参数形式建立覆岩采动应力、孔隙率与渗透率之间函数关系(式(1)～(4))，首先进行覆岩变形损伤稳态模拟分析，模拟分析弯曲带覆岩空间变形范围，以及覆岩孔隙率与渗透率变化参数;利用Mohr-Coulomb准则(式(5))识别出采掘扰动下导水裂隙的发育高度。
(2)利用COMSOL软件平台中大变形几何体自动重剖功能，形成采动变形二次剖分网格。
(3)在达西渗流模块中，对采动导水裂隙高度进行数值边界化处理(式(6));对采动渗透率参数输入，建立含水层地下水失水分析模型，通过瞬态或稳态分析，研究采掘扰动下地下水流场的响应规律。
2 曹家滩煤矿煤层开采覆岩损伤变形与含水层失水模型构建
2.1 研究区概况
曹家滩井田位于榆神矿区一期规划区中部，毛乌素沙漠东南缘，为沙丘沙地和风沙滩地、黄土梁峁地貌，地势相对平缓，地质构造简单，矿井主采侏罗系延安组第4段顶部的2-2煤层，平均煤厚达11.20 m，埋深一般为320 m左右，采煤方法为分层大采高综合机械化采煤工艺。
井田地下水划分为新生界松散岩类孔隙潜水含水层组与基岩裂隙承压含水层组两种含水类型，松散岩类孔隙潜水含水层组，由风积沙、冲洪积层与结构松散的上更新统的萨拉乌苏组孔隙类含水层组成，一般具有统一的水力联系，地下水富水性好，水位埋深一般小于3 m，是矿区煤层开采保水目标的含水层;基岩裂隙承压含水层组主要由侏罗系中统安定组、直罗组及延安组基岩裂隙含水层构成，含水层组岩性以中细砂岩为主，富水性较差;其中在基岩顶部岩石顶面受到不同程度的风化，平均厚度26 m左右，结构杂乱，松软易碎，富水性较好，是本组的主要含水层。隔水层为第四系离石组黄土与新近系保德组红土层，平均厚度70 m，富水性差，是区内主要的相对隔水层组。
图3 覆岩变形损伤与含水层失水数值模型
Fig.3 Numerical model for analysis of deformation damage of overburden and water loss of aquifer
2.2 覆岩损伤变形与含水层失水模型构建
2.2.1 模型概化
如图3所示，根据该矿首采工作面典型覆岩结构，模型在垂向上由上之下概化为松散介质含水层(风积沙、冲洪积层与萨拉乌苏组含水层)、黏土类隔水层(离石组黄土与保德组红土层)、风化基岩含水层、安定组弱透水层、基岩裂隙含水层(直罗组与延安组砂岩层)以及2-2号煤层共6个模拟分层，各岩层主要参数见表1。为了削弱模型边界效应影响和直观分析，采煤工作面走向上的模型尺度为100 m，模型宽度设置为3 000 m，图幅显示比例调整为X∶Y∶Z=1∶10∶3，模拟分析百米进尺的覆岩损伤与失水规律。
表1 模型主要参数 Table 1 Rock mechanics and seepage parameters
岩层密度/(kg·m-3)弹性模量E/GPa黏聚力c/MPa内摩擦角φ/(°)泊松比ν初始孔隙率初始渗透系数K/(m·d-1)松散含水层2 3600.751.5350.420.261.50弱透水土层1 6801.151.2420.240.350.001风化基岩含水层2 1001.202.0360.300.251.25弱透水泥质砂质层2 7002.552.5390.290.280.01基岩含水层2 5803.343.0350.260.150.032-2煤层1 3503.253.5350.320.150.03
2.2.2 模型构建
(1)数学模型。本次以岩体弹塑性力学和地下水流体力学理论为基础建立研究问题的数学模型，包括以固体弹塑性力学为基础的位移场方程，描述采动覆岩损伤变形过程;以达西定律为基础的地下水水流方程，描述地下水动态响应机制。为了便于模型位移方程与渗流方程进行耦合，将地下水水头高度(m)均换算为水头压力(MPa)，其数学模型可简化为:
地下水渗流场方程:
岩体位移场方程:
*(u)=0
(8)
式中，K′为多孔介质渗透率,10-15m2;ρ为流体密度,g/cm3;G为切变模量，G=2E(1+ν);ν为介质的排水泊松比;E为介质弹性模量;u分别为x轴、y轴及z轴方向的位移。
(2)定解条件。模型中定义底部为“固定边界”，即限制单元在三维空间上的变形，即位移为0;模型顶部的松散层潜水位界面与地表大气联通，而井下采掘扰动形成的导水裂隙带直接与井下采掘空间的大气相通，定义潜水位界面和导水裂隙带圈闭的界面均为水压边界，其边界压力值近似为大气压(值为0.1 MPa);采宽较模型研究范围相对较小，模型两侧在水平方向变形基本可以忽略，因而定义为“辊边界”，即允许边界单位仅在垂向上有位移产生，水平方向位移为0。
3 曹家滩煤矿煤层开采覆岩损伤变形与含水层失水分析
3.1 采动覆岩损伤变形分析
如图4所示，模拟采高一次采全高工况下(采高10 m)采动应力分布、覆岩变形、孔隙率与渗透系数,以及塑性破坏区的空间变化特征。
如图4所示，基于弹塑性力学模拟得出弯曲带内各岩层最大沉陷变形位于采空区顶部，松散层、风化层、安定组最大沉陷整体变形约为3.53 m;模型中定义各含水层为各向同性，岩层的孔隙率和渗透系数初值见表1，采动后孔隙率和渗透系数模拟结果如图4(c),(d)所示，其中位于弯曲带的弱透水土层在工作面两侧由于受拉伸变形影响其渗透系数由0.001 m/d增加至0.001 2 m/d，增加了12%，而工作面顶部由于受压缩变形影响其渗透系数减少至0.000 6 m/d，总的来看弯曲带覆岩变形与覆岩渗透变异程度不大。
如图5(a)所示，模拟得出一次采全高(10 m)工况下最大发育高度为202 m，裂采比为20.2，同时模拟了上分层开采5 m工况下导水裂隙最大发育高度为128 m，裂采比为25.6(图5(b))，与文献[23]曹家滩煤矿ANSYA数值模拟成果，以及相邻井田榆树湾、金鸡滩、杭来湾煤矿的2-2煤层导水裂隙高度实测成果一致性好(表2)。
3.2 含水层失水分析
在工作面覆岩变形损伤模拟输出的基础上，对覆岩变形岩层重新剖分(图6)、孔隙率与渗透率参数自动重赋，为了便于数值处理，将识别得出的导水裂隙带范围概化成倒梯形，梯形各面定义为0水压力边界，建立采动条件下含水层地下水失水分析模型，通过瞬态模拟1 a期末研究采掘扰动下地下水流场的响应规律，在失水模型中分别对梯形各面进行达西流速积分，积分可得出来自导水裂隙侧面和顶面各含水层的地下水流失量。
如图7(a)所示，上分层开采条件下(采高5 m)导水裂隙最大发育高度为128 m，发育至直罗与延安组基岩含水层内部，直罗组合延安组基岩含水层形成明显的了以导水裂隙为中心的地下水降落区，见表3，通过积分得出各含水层向导水裂隙的失水总量为35.84 m3/h，其中来自直罗与延安组基岩含水层侧向失水量约为23.17 m3/h，由于基岩含水层与风化含水层间的水力梯度增加，来自导水裂隙顶部含水层(含风化基岩含水层、土层与松散含水层)地下水的垂向渗漏漏失量约为12.67 m3/h，地表松散含水层由于土层的保护，地下水流场基本未发生变化，其水头压力最大减少了0.015 MPa。
图4 煤层开采影响下弯曲带覆岩变形特征
Fig.4 Deformation characteristics of overburden rock affected by coal seam mining
图5 煤层开采影响下覆岩塑性损伤范围
Fig.5 Deformation and damage characteristics of overburden rock affected by coal seam mining
表2 曹家滩相邻矿井煤层开采导水裂隙高度统计
Table 2 Statistical table of water-conducting fractured zone height in adjacent coal seams of Caojiatan Coal Mine
矿井工作面测试孔号采高/m导水裂隙高度/m裂采比研究方法曹家滩[23]首采工作面—5.0123.4024.68数值模拟金鸡滩12101JSD25.5107.4919.5420104Y35.0130.5026.10实测Y55.0138.9027.78榆树湾[24]20106Y45.0137.3027.4620102—5.090.0018相似材料模拟10.0≥285≥28.5H34.5108.3224.10杭来湾30101H44.5114.3825.40实测H54.5107.8323.96H74.593.8720.86
图6 煤层开采影响下覆岩变形剖分示意(X∶Y∶Z=1∶10∶3)
Fig.6 Diagram of overburden deformation subdivision under the influence of coal seam mining
图7 煤层开采影响下覆岩含水层地下水流场特征
Fig.7 Characteristics of groundwater flow field affected by coal seam mining
表3 模型主要参数
Table 3 Rock mechanics and seepage parameter
工况导水裂隙高度/m揭露层位失水量/(m3·h-1)总失水量侧向垂向松散层水头压力变化/MPa上分层开采128基岩含水层35.8423.1712.67≤0.015一次采全高202风化基岩含水层底部130.3192.6537.66≤0.05
如图7(b)所示，一次采全高(10 m采高)条件下导水裂隙最大发育高度为202 m，直接贯通了直罗组与延安组基岩含水层，发育至风化基岩含水层底部，通过积分得出各含水层向导水裂隙的失水总量达到了130.31 m3/h，其中直罗、延安组基岩含水层与风化基岩含水层侧向失水量约为92.65 m3/h，由于风化基岩含水层与近地表松散含水层间的水力梯度增加，导水裂隙顶部的含水层(含风化基岩、土层与松散含水层)地下水的越流漏失量约为37.66 m3/h，其中松散含水层水头压力最大减少了0.05 MPa，导致地表松散含水层地下水流场发生了明显变化，在工作面顶部地下水流向逆转。
从曹家滩案例矿井的模拟分析成果可以看出，分层开采的导水裂隙不揭露风化基岩含水层时，近地表的松散含水层受采动影响极小，一次采全高的导水裂隙由于揭露风化基岩含水层，导致含水层失水总量增加明显，近地表的松散含水层由于垂向水力梯度的增加地下水流场亦发生较大变化，因此榆神矿区中深煤层通过限高、分层开采抑制导水裂隙发育高度是实现保水采煤一项重要措施。
4 结 论
(1)榆神矿区中深煤层开采易形成典型的“侧向直接与垂向渗漏”复合失水模式，笔者以COMSOL多物理场耦合分析为平台，基于岩体弹塑性力学分析方法，通过建立采动应力、孔隙率与渗透率耦合关系计算弯曲带覆岩采动渗透系数，提出了将采动导水裂隙带边界数值化处理、以及弯曲变形重新剖分的煤层开采覆岩变形损伤与含水层失水模型的构建方法。
(2)模拟分析了榆神矿区曹家滩煤矿工作面一次采全高条件下弯曲带最大沉陷量为3.5 m，弱透水土层渗透系数增加了12%，导水裂隙最大发育高度为202 m，发育至风化基岩含水层底部，失水总量达到131 m3/h，加剧了松散含水层与风化基岩含水层间的垂向流动，对松散含水层影响较大。分层开采(5 m采高)条件下导水裂隙最大发育高度为128 m，发育至直罗与延安组基岩含水层内部，含水层失水总量35.84 m3/h，对近地表松散含水层影响小。
[1] 王双明,黄庆享,范立民,等.生态脆弱区煤炭开发与生态水位保护[J].中国煤炭地质,2010,35(1):7-14.
WANG Shuangming,HUANG Qingxiang,FAN Limin,et al.Coal development and ecological water level protection in ecologically fragile areas[J].Coal Geology of China,2010,35(1):7-14.
[2] 仵拨云,彭捷,向茂西,等.榆神府矿区保水采煤受保护萨拉乌苏组含水层研究[J].采矿与安全工程学报,2018,35(5):984-990.
WO Boyun,PENG Jie,XIANG Maoxi,et al.Study on the aquifer of the Salawusu Formation protected by water-retaining mining in Yushenfu mining area[J].Journal of Mining & Safety Engineering,2018,35(5):984-990.
[3] 缪协兴,王长申,白海波.神东矿区煤矿水害类型及水文地质特征分析[J].采矿与安全工程学报,2010,27(3):285-291,298.
MIAO Xiexing,WANG Changshen,BAI haibo.Analysis of water hazard types and hydrogeological characteristics of Shendong coal mine[J].Journal of Mining and Safety Engineering,2010,27(3):285-291,298.
[4] FAN Limin,MA Xiongde.A review on investigation of water-preserved coal mining in western China[J].International Journal of Coal Science & Technology,2018,5(4):411-416.
[5] HU Zhenqi,CHEN Chao,XIAO Wu,et al.Surface movement and deformation characteristics due to high-intensive coal mining in the windy and sandy region[J].International Journal of Coal Science & Technology,2016,3(3):339-348.
[6] 顾大钊,张建民.西部矿区现代煤炭开采对地下水赋存环境的影响[J].煤炭科学技术,2012,40(12):114-117.
GU Dazhao,ZHANG Jianmin.Impact of modern coal mining on groundwater environment in western mining area[J].Coal Science and Technology,2012,40(12):114-117.
[7] 王连国,王占盛,黄继辉,等.薄基岩厚风积沙浅埋煤层导水断裂带高度预计[J].采矿与安全工程学报,2012,29(5):607-612.
WANG Lianguo,WANG Zhansheng,HUANG Jihui,et al.Prediction on the height of water-flowing fractured zone for shallow seam covered with thin bedrock and thick windblown sands[J].Journal of Mining & Safety Engineering,2012,29(5):607-612.
[8] 许家林,朱卫兵,王晓振.基于关键层位置的导水断裂带高度预计方法[J].煤炭学报,2012,37(5):762-769.
XU Jialin,ZHU Weibing,WANG Xiaozhen.New method to predict the height of fractured water-conducting zone by location of key strata[J].Journal of China Coal Society,2012,37(5):762-769.
[9] LI Wenping,WANG Qiqing,LIU Shiliang,et al.Study on the creep permeability of mining-cracked N2 laterite as the key aquifuge for preserving water resources in Northwestern China[J].International Journal of Coal Science & Technology,2018,5(3):315-327.
[10] 王双明,黄庆享,范立民,等.生态脆弱矿区含(隔)水层特征及保水开采分区研究[J].煤炭学报,2010,35(1):7-14.
WANG Shuangming,HUANG Qingheng,FAN Limin,et al.Study on the characteristics of aquifers and zoning of water-holding mining in ecologically fragile mining areas[J].Journal of China Coal Society,2010,35(1):7-14.
[11] 王启庆,李文平,李涛.陕北生态脆弱区保水采煤地质条件分区类型研究[J].工程地质学报,2014,22(3):515-521.
WANG Qiqing,LI Wenping,LI Tao.Division types of geological conditions at mining with water protection in ecological fragile area of northern Shannxi[J].Journal of Engineering Geology,2014,22(3):515-521.
[12] 顾大钊.煤矿地下水库理论框架和技术体系[J].煤炭学报,2015,40(2):239-246.
GU Dazhao.Theory framework and technological system of coal mine underground reservoir[J].Journal of China Coal Society,2015,40(2):239-246.
[13] 张杰,侯忠杰.榆树湾浅埋煤层保水开采三带发展规律研究[J].湖南科技大学学报(自然科学版),2006,21(4):10-13.
ZHANG Jie,HOU Zhongjie.Research on the development law of three zones of water-conserving mining in shallow coal seam of Yushuwan[J].Journal of Hunan University of Science and Technology(Natural Science Edition),2006,21(4):10-13.
[14] 彭小沾,崔希民,李春意,等.陕北浅煤层房柱式保水开采设计与实践[J].采矿与安全工程学报,2008,25(3):301-304.
PENG Xiaozhan,CUI Ximin,LI Chunyi,et al.Design and practice of room & pillar water-preserved mining for shallowly buried coal seam in north of Shaanxi Province[J].Journal of Mining & Safety Engineering,2008,25(3):301-304.
[15] 黄庆享,张文忠.浅埋煤层条带充填隔水岩组力学模型分析[J].煤炭学报,2015,40(5):973-978.
HUANG Qingxiang,ZHANG Wenzhong.Mechanical model of water resisting strata group in shallow seam stripfilling mining[J].Journal of China Coal Society,2015,40(5):973-978.
[16] 马立强,张东升,王烁康,等.“采充并行”式保水采煤方法[J].煤炭学报,2018,43(1):62-69.
MA Liqiang,ZHANG Dongsheng,WANG Shuokang,et al.Water-preserved mining with the method named“backfilling while mining”[J].Journal of China Coal Society,2018,43(1):62-69.
[17] PIRNIA P,DUHAIME F,ETHIER Y,et al.ICY:An interface between COMSOL multiphysics and discrete element code YADE for the modelling of porous media[J].Computers & Geosciences,2018,123:38-46.
[18] LOUIS C.Rock hydroulics[A].Rock Mechanics[C].LMuller,1974.
[19] 朱万成,魏晨慧,张福壮,等.流固耦合模型用于陷落柱突水的数值模拟研究[J].地下空间与工程学报,2009,5(5):928-933.
ZHU Wancheng,WEI Chenhui,ZHANG Fuzhuang,et al.Investigation of water inrush from karst subsidence column by using a coupled hydromechanical model[J].Chinese Journal of Underground Space & Engineering,2009,5(5):928-933.
[20] PALMER I.Permeability changes in coal analytical modeling[J].International Journal of Coal Geology,2009(77):119-126.
[21] RUTQVIST J,TSANG C F.A study of cap rock hydromechanical changes associated with CO2-injection into a brine formation[J].Environmental Geology,2002,42:296-305.
[22] 赵春虎,虎维岳,靳德武.西部干旱矿区采煤引起潜水损失量的定量评价方法[J].煤炭学报,2017,42(1):169-174.
ZHAO Chunhu,HU Weiyue,JIN Dewu.Method of quantitative evaluation on amount of groundwater loss from unconfined aquifer caused by mining disturbance in the arid area of Western China[J].Journal of China Coal Society,2017,42(1):169-174.
[23] 张晨招.曹家滩煤矿采动裂隙带与矿井涌水数值模拟研究[D].西安:西北大学,2013.
ZHANG Chenzhao.Study on Numerical Simulation of mining fracture zone and mine water inflow in Caojiatan Coal Mine[D].Xi’an:Northwest University,2013.
[24] 师本强.陕北浅埋煤层砂土基型矿区保水开采方法研究[J].采矿与安全工程学报,2011,28(4):548-552.
SHI Bengqiang.Research on water-preserved-mining in shallow seam covered with rock soil and sand in Northern Shaanxi[J].Journal of Mining and Safety Engineering,2011,28(4):548-552.
