一种采动覆岩导水通道演化的模拟方法

1.本发明属于保水采煤技术领域，具体涉及一种采动覆岩导水通道演化的 模拟方法。


背景技术：

2.晋陕蒙地区蕴含我国近60％的煤炭资源，却仅拥有我国8％的水资源。更 严重的是，高强度的煤炭地下开采会导致上覆岩层移动和破坏，从而对导水 裂隙带影响范围内的含水层产生影响，导致地下水漏失、生态环境破坏。这 为实现“黄河流域生态保护和高质量发展”目标带来巨大挑战。大同矿区燕 子山矿的工程实践表明：煤系地层包含有脆性的砂岩、灰岩以及韧性的泥岩、 高岭岩，同时这些地层中包含有大量的主、次结构面。地层复杂的力学性质、 丰富的结构面会对采动覆岩导水通道分布、高度及其演化产生重大影响，决 定了侏罗系采空区积水能否下泻至二叠系采煤工作面。采用合理的数值模拟 方法预测采动覆岩导水通道发育高度，对于优化工作面参数、保水采煤至关 重要。
3.近年来，国内外学者主要采用有限元fem和离散元dem方法开展采动覆 岩运动模拟研究。其中：
4.fem方法基于连续介质假设，采用弹性、塑性、损伤力学理论，研究采动 过程中覆岩的力场、塑性区、损伤带等的发育高度与演化规律。
5.dem方法基于离散介质假设，大多将覆岩中的结构面设置成规则的砌体结 构形式，覆岩结构面、基质多分别采用弹脆性、弹性力学本构模型，进而数 值模拟得到煤层开采过程中覆岩内水平裂隙、竖直裂隙的发育高度。但是， 覆岩中的结构面随机分布、且尺度不一，在数值计算模型中需建立多级结构 面体系及其随机生成程序；与此同时，覆岩中的泥岩、高岭岩具有显著的韧 性破坏特征。
6.上述dem方法会改变采动覆岩导水通道传播路径，并对其分布与演化规 律产生决定性作用；而fem方法会在很大程度上消耗次生应力产生的能量， 影响覆岩裂隙发育高度。


技术实现要素：

7.鉴于此，为解决上述背景技术中所提出的问题，本发明的目的在于提供 一种采动覆岩导水通道演化的模拟方法。
8.为实现上述目的，本发明提供如下技术方案：一种采动覆岩导水通道演 化的模拟方法，包括：
9.s1.确定覆岩中结构面的几何参数，并建立关于覆岩中结构面的三维网络 几何模型；所述结构面包括主要结构面及次要结构面；
10.s2.建立包含结构面的采场覆岩导水通道演化的数值计算模型；
11.s3.确定零厚度裂缝单元的尺寸，分别在所述数值计算模型的覆岩结构面 处及基质内实体单元边界处嵌入所述零厚度裂缝单元，以此模拟覆岩中主要 结构面及次要结构面的随机分布；
12.s4.分别获取岩层基质的弹塑性特征、结构面的脆性/韧性断裂力学特征、 完全断裂后岩块的剪切摩擦特征，将各特征赋予所述数值计算模型的结构面 和裂缝单元中，并采用fem与dem耦合计算；
13.s5.获取并输入所述数值计算模型的材料参数、边界条件和初始条件；
14.s6.编制煤层位置处单元随时间和空间位置删除的vusdfld程序，嵌入所 述数值计算模型，以模拟采煤工作面回采过程；
15.s7.求解所述数值计算模型的模拟结果，并输出；
16.s8.分析所述模拟结果，得到采动覆岩导水裂隙带空间分布及发育高度随 回采参数的演化规律，将所述演化规律与地表、地下含水层分布状态相结合， 优化煤炭开采方案。
17.优选的，在所述步骤s1的确定覆岩中结构面的几何参数中包括：
18.根据电阻率测井法确定各岩层中主要结构面的间距；
19.根据钻井取芯过程中各岩层岩芯的平均长度确定次要结构面的间距；
20.计算结构面的平均间距、间距均方差、以及在预设倾角范围内结构面的 出现频率。
21.优选的，在所述步骤s1的建立关于覆岩中结构面的三维网络几何模型中 包括：
22.根据所需求建立的数值计算模型的大小，确定在x、y、z三个正交方向 上voronoi多边形的个数；
23.根据结构面的间距均方差，在python中生成一组满足统计分布的空间离 散点；
24.将所有的空间离散点连接成delaunay三角形；
25.依据delaunay空间外接球准则找出三角形三点之外的第四点，并构成一 个delaunay四面体；
26.由四面体相邻面的中垂面围成voronoi多面体，由此生成覆岩中结构面 的三维网络几何模型。
27.优选的，在所述步骤s3的确定零厚度裂缝单元的尺寸中包括：
28.断裂过程区的长度l：其中，π为圆周率常数，μ为材料 泊松比，d
0，c
为弹性模量，g
n
为断裂能，为材料的抗拉强度；
29.钻井取芯过程中各岩层岩芯的平均长度t；
30.在大尺度结构面中，裂缝单元的尺寸小于等于l；
31.在小尺度结构面中，裂缝单元的尺寸为平均长度t与长度l中的小者。
32.优选的，在所述步骤s3的嵌入所述零厚度裂缝单元中包括：
33.(1)更新实体单元节点编号：
34.对所述步骤s2中建立的数值计算模型进行全局网格划分，划分后，在 job模块下的write input中生成.inp脚本文件，在脚本文件中的关键字 *nset和*elset下方数字中分别找到实体单元的最大节点编号、最大单元 编号，并分别记为n
max
，e
max
；
35.将实体单元编号为e
j
的第n
i
(n
i
≠1)个节点复制a个，且a为共用实 体单元的个数；
36.在节点坐标相同的前提下，保持最小实体单元编号的节点编号不变， 剩余的a
‑
1个节点编号修改为：n
i
＝(a
‑
1)
×
10
ten(nmax)
+n
i
；其中，ten(n
max
)为n
max
的十进制位数；具体的，
经过上述变换，第n
i
个节点被复制了a次，编号 依次为n1,n2,
…
,n
a
；而共用第n
i
个节点的单元也有a个，单元编号依次为 e1,e2,
…
,e
a
，将新的节点编号复制在对应单元编号位置处，即可形成一一 对应的两个有序数组elem＝(e1,e2,
…
,e
b
)和node＝(n1,n2,
…
,n
b
)，将更新后 的有序数组复制在.inp脚本文件中，实现实体单元节点号的更新；
37.(2)创建裂缝单元节点编号、裂缝单元编号
38.在初始.inp脚本文件中，复制实体单元编号相同的节点，并按照下式 确定的顺序组成新的数组：
39.group＝(e
j
,e
k
,n
shared01
,n
shared02
,n
shared03
,n
shared04
,node
01
,node
02
,node
03
,node
04
)；其中，
40.数组group中的e
j
和e
k
为相邻实体单元编号；n
shared01
‑
n
shared04
为初始.inp文 件中的相邻实体单元之间的相同编号的节点；在node
01
‑
node
04
中，仅保留 由n
shared01
‑
n
shared04
更新的实体单元节点编号，然后逆序排列，得到单元e
j
和 e
k
之间的裂缝单元节点编号；
41.裂缝单元编号根据实体单元最大编号确定：e
coh,min
＝k+e
i
；其中，k为 10的整数幂次方倍，且k大于实体单元最大编号；例如，实体单元最大节 点编号为756，k值可取1000，也可取10000；
42.将裂缝单元节点编号与裂缝单元编号组合形成一个裂缝单元编号集 合；具体的，在.inp文件中添加关键字“*element type＝coh2d4”(以四 边形单元为例，若为六面体单元，则修改为“*element type＝coh3d8”)， 同时在关键字“*part”与“*assembly”之间添加“*elset＝cohelemall”， 而后按照先单元编号、后节点号的顺序写入相应数字，从而将所有裂缝单 元组成一个名称为“cohelemall”的集合；
43.(3)在所述数值计算模型的主要结构面处嵌入裂缝单元
44.复制更新后的目标裂缝单元编号e
j
以及对应节点坐标n
j
(x
j
，y
j
，z
j
)， 并采用内积法求得该目标裂缝单元的方向向量d
coh
＝(d1,d2,d3)；
45.判断向量d
coh
是否与所述数值计算模型中主要结构面的方向向量d
ip
共 线、判断向量d
coh
中的任一节点是否位于主要结构面上，是则判定所述目标 裂缝单元位于主要结构面内。
46.具体的，对于满足上述条件的目标裂缝单元编号、节点编号复制在.inp 文件中，在关键字“*part”与“*assembly”之间添加“*elset＝cohelemdis”， 从而将所有结构面上的单元组成一个名称为“cohelemdis”的集合。
47.优选的，在所述步骤s4的获取岩层基质的弹塑性特征中包括：
48.确定岩层基质在采动过程中发生的弹塑性变形特征
[0049][0050][0051]
其中：
[0052]
σ
s
为应力矢量、为塑性应变、e0为弹性刚度矩阵、ε
s
为总应变；
[0053]
a＝(σ
b0
/σ
c0
‑
1)/(2σ
b0
/σ
c0
‑
1)；b＝σ
s
/σ
t
·
(1
‑
a)
‑
(1+a)；c＝3(1
‑
k
c
)/(2k
c
‑
1)；
[0054]
p
eff
＝trace(σ
s,eff
)/3,s
eff
＝σ
‑
p
eff
i；
[0055]
(i＝1,2,3)是x、y、z三个正交方向上的最大主应力；
[0056]
σ
b0
/σ
c0
是三轴与单轴压缩条件下屈服应力之比；
[0057]
σ
t
和σ
s
是等效拉压应力；k
c
取值为0.667。
[0058]
优选的，在所述步骤s4的获取结构面的脆性/韧性断裂力学特征、完全断 裂后岩块的剪切摩擦特征中包括：
[0059]
(1)确定覆岩中的非贯通结构面出现的弹性变形特征
[0060]
σ
c
＝d
0,c
ε
c
；ε
c
＝s/t0；
[0061]
其中、d
0，c
为弹性刚度矩阵、ε
c
为应变矢量、t0为裂缝单元的本构厚度， 且t0取值为实体单元特征尺寸的0.01倍，s为裂缝单元的位移矢量；
[0062]
(2)确定覆岩中的非贯通结构面出现的脆性断裂与韧性断裂特征
[0063]
脆性断裂特征：ε
c
＝s/t0；其中，σ
c,l
和分别为牵引力及牵引力峰值，l表示为法向和两个切向方向；
[0064]
韧性断裂特征：
[0065][0066]
其中，s
n
，s
s
和s
t
均为变量，并分别表示法向和两个切向方向的位移；γ
n
,γ
s
为 断裂能常数，且表达式为：
[0067][0068]
g
s
＝g
s
+g
t
；χ
n
＝s
n,p
/s
n
,χ
s
＝s
s,p
/s
s
；
[0069]
式中，g
n
与g
s
分别为法向与总切向断裂能；g
s
和g
t
分别为两个切向方向 的断裂能；β与γ分别为纯i型和纯ii型断裂实验的应力
‑
应变曲线形状参 数；s
n,p
和s
s,p
分别为纯i型和纯ii型断裂实验的应力
‑
应变曲线峰值位移； s
n
和s
s
分别为纯i型和纯ii型断裂实验的应力
‑
应变曲线断裂位移；
[0070]
(3)确定完全断裂后岩块的剪切摩擦特征
[0071]
σ
s
＝μ1σ
nn
；其中，σ
s
为剪切应力；μ1为摩擦系数；σ
nn
为法向应力。
[0072]
优选的，在所述步骤s4的采用fem与dem耦合计算中包括：
[0073]
确定边界条件，并通过fem求解所述数值计算模型中覆岩基质内实体单 元的节点力；
[0074]
通过实体单元与裂缝单元的共享节点传递力场数据，并采用显式求解法 计算裂缝单元的位移分量；
[0075]
在s＞s
f
时，判定裂缝单元断裂，且dem响应；其中s
f
为混 合型断裂模式下的断裂位移；
[0076]
在dem中，通过newton
‑
raphson迭代方法求解，且迭代求解的表达式为： jacobi(u
j
)du
j+1
＝
‑
r
j
；其中，jacobi(u
j
)为第j迭代步时的雅可比矩阵，r
j
为 第j迭代步时的残差；迭代收敛条件为：
[0077]
||r||2＜tolerance并且||u||2＜tolerance；其中||||2为向量的二型范数，且容差 tolerance取值为10
‑8时，迭代2
‑
5次后收敛。
[0078]
优选的，在所述步骤s6的模拟采煤工作面回采过程中包括：
[0079]
在materials模块中定义状态变量参数depvar＝1；
[0080]
在step模块中设置输出变量中包含完全断裂的裂缝单元状态参量 status；
[0081]
在vusdfld程序中输入总时间totaltime、分析步时间steptime、单元 号nblock、单元积分点坐标coordmp，且当煤层中的裂缝单元满足设定条件 时，指定statenew(nblock,1)＝0，以将该完全断裂的裂缝单元删除；
[0082]
所述设定条件为：
[0083]
f(totaltime,coordmp)＝steptime
×
v+coordmp(nblock,1)
[0084]
其中v表示采煤工作面推进速度；coordmp(nblock,1)中的“1”表示所 述数值计算模型中的煤层开挖方向。
[0085]
本发明与现有技术相比，具有以下有益效果：
[0086]
(1)本发明为地下煤层开采引起覆岩破断、导水裂隙带分布与演化提供了 有效的模拟工具。具体，采用本发明所提供的fem+dem方法，可以在覆岩的 所有主要、次要结构面上嵌入裂缝单元，实现采动过程中裂缝在覆岩中的任 意扩展。另外，fem+dem方法是对传统有限元、离散元、拉格朗日元方法的重 大改进，并结合覆岩基质的弹塑性特征，以及反映覆岩结构面及微裂隙的脆 性、韧性断裂力学特征，从计算方法和力学理论上，实现采动过程中覆岩从 完整
→
破裂
→
沿破裂面剪切滑移的全过程数值模拟；综上，有效获得采动过 程中覆岩垮落带、导裂带的分布及演化，为开采参数的优化提供模拟方法。
[0087]
(2)本发明方法针对实体单元采用隐式时间积分求解算法，针对裂缝单元 采用基于显式时间积分求解算法，可以大大提高数值解的收敛性，同时加快 计算效率。
[0088]
(3)本发明方法以测井资料与岩芯长度为依据，分别建立主要结构面与次 要结构面模型，在两个位置处嵌入不同材料参数的裂缝单元，从而有效模拟 大型结构面的剪切滑移、结构面之间完整岩体的随机碎裂过程，由此实现开 采过程中，煤层覆岩中宏细观裂缝的随机扩展、随机分布。
[0089]
(4)本发明可用于分析不同开采参数、不同地质参数影响下采动覆岩导水 裂隙分布及演化规律，从而优化开采参数、实现保水采煤。
附图说明
[0090]
图1为本发明模拟方法的流程图；
[0091]
图2为本发明覆岩中结构面的三维网络几何模型示意图；
[0092]
图3为本发明全局嵌入裂缝单元的方法原理图；
[0093]
图4a为本发明主要结构面处嵌入的裂缝单元示意图；
[0094]
图4b为本发明主要结构面处嵌入的裂缝单元示意图；
[0095]
图4c为本发明基质内实体单元示意图；
[0096]
图5为采动覆岩裂隙分布数值计算结果示意图；
[0097]
图6为主要导水裂隙带分布位置示意图；
[0098]
图7及图8为覆岩mises应力云图及相应的采动裂隙分布示意图；
具体实施方式
[0099]
下面结合附图对本发明做进一步说明，但本发明的保护范围不局限于以 下所述。
[0100]
s1.以大同矿区燕子山矿为背景；
[0101]
根据电阻率测井法确定各岩层中主要结构面的间距；
[0102]
根据钻井取芯过程中各岩层岩芯的平均长度确定次要结构面的间距；
[0103]
计算结构面的平均间距、间距均方差、以及在预设倾角范围内结构面的 出现频率；
[0104]
根据所需求建立的数值计算模型的大小，确定在x、y、z三个正交方向 上voronoi多边形的个数；具体，并结合燕子山矿工程地质情况，x＝870m， 高y＝451m，z＝1m；
[0105]
根据结构面的间距均方差，在python中生成一组满足统计分布的空间离 散点；
[0106]
将所有的空间离散点连接成delaunay三角形；
[0107]
依据delaunay空间外接球准则找出三角形三点之外的第四点，并构成一 个delaunay四面体；
[0108]
由四面体相邻面的中垂面围成voronoi多面体，由此生成如图2所示的 覆岩中结构面的三维网络几何模型。
[0109]
s2.建立包含结构面的采场覆岩导水通道演化的数值计算模型。
[0110]
s3.确定零厚度裂缝单元的尺寸，分别在数值计算模型的覆岩结构面处及 基质内实体单元边界处嵌入零厚度裂缝单元，以此模拟覆岩中主要结构面及 次要结构面的随机分布；其中
[0111]
确定零厚度裂缝单元的尺寸时：
[0112]
断裂过程区的长度l：其中，π为圆周率常数，μ为材料 泊松比，d
0，c
为弹性模量，g
n
为断裂能，为材料的抗拉强度；其中弹性模 量、断裂能和抗拉强度均通过三点弯曲实验获得；
[0113]
钻井取芯过程中各岩层岩芯的平均长度t；
[0114]
上述，在大尺度结构面中，裂缝单元的尺寸小于等于l；在小尺度结构面 中，裂缝单元的尺寸为平均长度t与长度l中的小者。
[0115]
嵌入零厚度裂缝单元时：
[0116]
在mesh模块下采用六面体网格对数值模型整体划分网格；而后通过图3 所示的全局嵌入裂缝单元的方法，分别在数值模型的主要结构面处(命名为 cohelemdis)、次要结构面处(命名为cohelemall)嵌入三维8节点的裂缝 单元(如图4a和图4b)，其他剩余部分为8个节点的六面体实体单元(命名 为baseelem，如图4c)。
[0117]
s4.在数值计算模型，根据岩层实际分布情况，通过减选方法将单元集合 baseelem、cohelemall、cohelemdis分别分解为baseelem
‑
sy(砂岩实体单 元集合)、
cohelemall
‑
sy(砂岩次要结构面裂缝单元集合)、cohelemdis
‑
sy (砂岩主要结构面裂缝单元集合)，以及相应的泥岩(ny)、煤(m)、高岭 岩(gl)等的裂缝单元、实体单元集合；分别输入各岩层基质(表1)、主要 结构面(表2)、次要结构面(表3)的材料参数，并赋予相应的实体单元、 裂缝单元。
[0118]
s5.在assembly模块下对数值计算模型除上表面以外5个边界的法向位 移进行约束，上表面为自由边界；而后根据采煤工作面的推进速度和模型尺 寸，在step模块下设定70个分析步(表示70个开挖步)，每个分析步均设 置为几何非线性，时间长度为1(表示1d)，时间增量设置为0.001；设置输 出变量为应力、应变、摩擦位移、摩擦应力、状态变量status等；通过自平 衡方式确定模型的初始应力。
[0119]
s6.在materials模块的variable number controlling element deletion 选项中定义状态变量参数depvar＝1。在vusdfld子程序中传入主程序计算得 到的总时间totaltime、分析步时间steptime、单元号nblock、单元积分点 坐标coordmp。
[0120]
按照采煤工作面推进速度为10m算，设定子程序中的关键语句为
[0121]
do k＝1,nblock
[0122]
if(f(totaltime,coordmp).eq.(steptime
×
10m+coordmp(nblock,1)))
[0123]
statenew(nblock,1)＝0
[0124]
end if
[0125]
end do
[0126]
由上，即可实现任意分析步下指定位置处单元删除，从而实现采煤工作 面以10m/d的速度向前推进的目的。
[0127]
s7.在job模块下新建一个任务，在user subroutine file选项卡中输 入编制好的vusdfld子程序；在edit job选项卡中设置并行计算核心数量为 48核，提交计算，求解数值计算模型的模拟结果，并输出。
[0128]
s8.分析模拟结果，得到采动覆岩导水裂隙带空间分布及发育高度随回采 参数的演化规律，将演化规律与地表、地下含水层分布状态相结合，优化煤 炭开采方案。
[0129]
具体的，在步骤s7中计算输出模拟结果包括：
[0130]
如图5所示的采动覆岩裂缝分布与演化规律；从计算结果及图示可以看 出，模型中上四层煤(侏罗系煤层)开采会引起其上覆岩层完全碎裂，导水 裂隙带贯通含水层、直达地表；而模型最下方煤层(石炭系煤层)开采后， 也会产生相应的覆岩裂隙，垮落带高度约为26m，导水裂隙带延伸至侏罗系采 空区，可能引发采空区积水下泻。但由于高岭石显著的塑性变形，导致高岭 石岩层区域的采动裂隙不甚发育；同时考虑到高岭石本身渗透率很低，而且 吸水膨胀后会导致裂隙开度降低，因此，判断采空区积水会下泻，但涌水量 有限。
[0131]
如图6所示的主要导水裂隙带分布位置；由图可知，采空区范围内垮落 的覆岩被其上部岩层压实紧密，导致渗透率较低；而在采空区边缘，往往形 成宽度较大的导水裂隙带，因此该位置为防治水重点治理区域。
[0132]
如图7及图8所示的覆岩mises应力云图及相应的采动裂隙分布；具体， 图7为第52此迭代分析所形成的覆岩mises应力云图及相应的采动裂隙分布 示意图，图8为第63此迭代分析所形成的覆岩mises应力云图及相应的采动 裂隙分布示意图。由图可知，当采煤工作面接近直接顶的主要结构面时，会 导致顶板中mises应力快速增加；而当采煤工作面、直
接顶的主要结构面、 基本顶的主要结构面三者相接近时，mises应力值达到最大、且范围最广，导 致覆岩导水裂隙高度增加。
[0133]
综上所述，本发明与现有覆岩导水裂隙模型相比，一方面考虑了主要、 次要结构面分布对导裂带分布演化的影响，更加符合岩石几何非连续的特点； 另一方面，同时考虑了覆岩基质的弹塑性力学特征、覆岩结构面的脆性/韧性 断裂、以及完全断裂后的摩擦效应，结合fem
‑
dem数值模拟方法，成功模拟 了覆岩从连续介质转化为离散介质的过程，从而为采动覆岩导水裂隙演化模 拟提供了有力的工具。
[0134]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而 言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行 多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限 定。
[0135]
表1：覆岩基质力学参数
[0136] e/gpaμσ
b0
/σ
c0
a
s
/mpaσ
t
/mpak
c
δψ/
°
土0.020.312.22.50.120.6670.130粉砂岩10.70.211.5623.00.6670.140细砂岩13.80.221.5854.80.6670.142中砂岩20.60.221.41075.70.6670.142高岭岩7.40.301.8352.60.6670.130煤1.60.351.9141.10.6670.134泥岩8.00.241.6434.60.6670.135
[0137]
表2：主要结构面力学参数
[0138]
[0139]
表3：次要结构面力学参数
[0140]