CONTACT ALGORITHM OF 3D MECHANICAL FINITE ELEMENT MODELING FOR SEPARATION OF OVERBURDEN LAYERS AND ITS APPLICATION
-
摘要: 针对覆岩离层这一特殊的接触问题提出了新的接触算法，并采用大规模并行计算方法实现了对离层进行真正三维的力学模拟，从而得到覆岩离层空间几何形态参数。以济宁二号煤矿西部十一采区为研究实例，验证了算法的适用性及有效性。Abstract: Based on the specialty of separation layer, this paper proposes a new contact algorithm with the adoption of massive parallel computation for the study of separation of overburden layers in coal mine production. We have done the real 3D mechanic simulation and get the potential location, shape and volume for separation of overburden layers. Meanwhile, the study takes the area of No.2 coal mine of Jining as an example, which proves the applicability and validity of this new algorithm.
-
0. 引言
在煤炭开采过程中，随着工作面推移对覆岩体自然应力平衡状态的破坏，采空区上覆岩层的次生应力场重新分布，岩层内部的应力也重新分布以达到新的平衡。在此过程中上覆岩层要发生断裂、冒落、弯曲和下沉，这种运动变化由下向上逐步发展，出现不同步的下沉，上覆岩层由于岩性、距采空区距离、厚度等方面的不同，在垂直方向出现了不同的位移差，进而形成了离层空间。离层带中经常有大量的积水，这些离层水对于下伏的岩石形成巨大压力，当压力积累到一定程度时有可能破坏断裂带，形成离层水害，是巨大的煤矿安全隐患。因此离层的量化研究对煤矿安全生产和采区人民生活的影响具有积极的意义。
科研工作者对于覆岩离层问题进行了大量的研究。从最初的开采沉陷方式[1]，到岩层运动的物理及力学模型，岩层运动过程，离层发育规律，覆岩的离层条件，以及离层所带来的危害都做了很多探讨。杨伦等[2]提出的岩层二次压缩理论，将地表沉陷与岩层的物理力学性质联系起来。宋杨等[3]曾利用有限元方法研究支承压力显现的过程，揭示其变化规律与上覆岩层运动间的关系。苏仲杰[4]等将离层岩层假设为岩板建立了岩板挠度的力学模型。梁运培[5]提出了覆岩层移动的组合岩梁理论。张健全[6]分析了覆岩离层可能出现的位置，并提出了计算最大离层的方法。然而这些研究大多停留于定性阶段，在定量研究方面却刚刚起步，对于离层进行三维模型分析的研究尚属空白，其主要困难在于覆岩薄层众多，离层量小可至毫米级大可到米级，而三维空间体积巨大，现有的三维有限元接触算法无法解决。
针对上述问题，笔者进行了两个方面的研究，其一为层状界面弱接触类三维力学有限元算法，其二为应用该算法对三维力学有限元模拟的计算机程序和网格自动生成进行研究。在研究基础上进行的大规模并行计算结果表明，上述算法可行，计算成果对煤炭安全生产具有指导性意义。
1. 模型建立与接触算法
1.1 层内模型
假定岩层内部介质变形服从弹性力学基本规律，其力学变形服从弹性力学方程。
1.1.1 平衡方程
描述弹性体内部的力学平衡性质，根据牛顿第二定律可得出：
\left\{ \begin{array}{l} \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{xx}}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{xy}}}}{{\partial y}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{xz}}}}{{\partial z}} + {f_x} = 0\\ \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{yx}}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{yy}}}}{{\partial y}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{yz}}}}{{\partial z}} + {f_y} = 0\\ \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{zx}}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{zy}}}}{{\partial y}} + \frac{{\partial {\mathit{\boldsymbol{\sigma }}_{zz}}}}{{\partial z}} + {f_z} = 0 \end{array} \right. (1) 式中，σxx，σyy，σzz，σyz，σxz，σxy为应力矢量，Pa；fx，fy，fz为内力密度（N/m3），这里为重力；{fx fy fz}={0 0 -ρg}，ρ为密度（kg/m3），g为重力加速度（N/kg）。
1.1.2 几何方程
描述应变和变形（位移）的关系。介质变形的幅度远小于空间尺度可以忽略，因此本文采用了小变形的几何方程，进行准确的模拟。几何方程如下：
\left\{ \begin{array}{l} {\varepsilon _{xx}} = \frac{{\partial u}}{{\partial x}}\\ {\varepsilon _{yy}} = \frac{{\partial u}}{{\partial y}}\\ {\varepsilon _{zz}} = \frac{{\partial w}}{{\partial z}}\\ {\varepsilon _{yz}} = \frac{{\partial w}}{{\partial x}} + \frac{{\partial v}}{{\partial z}}\\ {\varepsilon _{xz}} = \frac{{\partial w}}{{\partial x}} + \frac{{\partial u}}{{\partial z}}\\ {\varepsilon _{xy}} = \frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}} \end{array} \right. (2) 式中，εxx，εyy，εzz，εyz，εxz，εxy为应变；u，v，w分别为x, y, z方向的位移，m。
1.1.3 本构方程
描述物体应力与应变（εxx，εyy，εzz，εyz，εxz，εxy）的关系，这里假定地下介质为线弹性体，E为杨氏模量，ν为泊松比，物性参数的选择直接决定了变形的分布情况，相关物性参数由钻井数据，以及查阅资料得出。本构方程如下：
\left( \begin{array}{l} {\mathit{\sigma }_{\mathit{xx}}}\\ {\mathit{\sigma }_{\mathit{yy}}}\\ {\mathit{\sigma }_{\mathit{zz}}}\\ {\mathit{\sigma }_{\mathit{yz}}}\\ {\mathit{\sigma }_{\mathit{xz}}}\\ {\mathit{\sigma }_{\mathit{xy}}} \end{array} \right) = \frac{E}{{\left( {1 + v} \right)\left( {1 - 2v} \right)}}\left( {\begin{array}{*{20}{c}} {1 - v}&v&v&{}&{}&{}\\ v&{1 - v}&v&{}&{}&{}\\ v&v&{1 - v}&{}&{}&{}\\ {}&{}&{}&{0.5 - v}&{}&{}\\ {}&{}&{}&{}&{0.5 - v}&{}\\ {}&{}&{}&{}&{}&{0.5 - v} \end{array}} \right)\left( \begin{array}{l} {\varepsilon _{xx}}\\ {\varepsilon _{yy}}\\ {\varepsilon _{zz}}\\ {\varepsilon _{yz}}\\ {\varepsilon _{xz}}\\ {\varepsilon _{xy}} \end{array} \right) (3) 1.1.4 边界条件
第一类边界条件适用于给定边界面位移的情况：
\begin{array}{l} u = {u_0}\\ v = {v_0}\\ w = {w_0} \end{array} (4) 其中，u0，v0，w0分别为初始位移在x, y, z三个方向上的分量，m。
第二类边界条件适用于给定边界面力密度的情形：
\left\{ \begin{array}{l} {T_x} = {\sigma _{xx}}{n_x} + {\sigma _{xy}}{n_y} + {\sigma _{xz}}{n_z} = {g_x}\\ {T_y} = {\sigma _{yx}}{n_x} + {\sigma _{yy}}{n_y} + {\sigma _{yz}}{n_z} = {g_y}\\ {T_z} = {\sigma _{zx}}{n_x} + {\sigma _{zy}}{n_y} + {\sigma _{zz}}{n_z} = {g_z} \end{array} \right. (5) 其中，gx，gy，gz为表面面力密度（Pa）在x, y, z方向的分量；nx, ny, nz为表面外法向的x, y, z方向分量。
对于如上弹性问题，其变分形式为：
\begin{array}{l} \int\limits_v {{\sigma _{xx}}\delta {\varepsilon _{xx}} + } {\sigma _{yy}}\delta {\varepsilon _{yy}} + {\sigma _{zz}}\delta {\varepsilon _{zz}} + {\sigma _{yz}}\delta {\varepsilon _{yz}} + {\sigma _{xz}}\delta {\varepsilon _{xz}} + {\sigma _{xy}}\delta {\varepsilon _{xy}}\\ = \int\limits_V {{f_x}\delta u + {f_y}\delta v + } {f_z}\delta w + \int\limits_F {{T_x}\delta u + {T_y}\delta v + } {T_z}\delta w \end{array} (6) 1.2 层间模型
假定介质为n层，则第i层和第i+1层之间需满足接触边界条件。以wi，1(x，y)表示第i层下边界垂直方向位移，wi, 2(x，y)表示第i层上边界垂直方向位移，则接触边界条件表示为：
0 ＞ {w_{i,1}}\left( {x,\mathit{y}} \right) ＞ {w_{i + 1,2}}\left( {x,\mathit{y}} \right) (7) 若第n层下边界为采空区，则有：
{W_{n,2}}\left( {x,\mathit{y}} \right) ＞ - D (8) 其中，D为采空高度，m。
1.3 接触问题有限元计算方法
首先假定相邻岩层之间不会发生强制粘结，各层可以发生自由形变，因此各层可进行独立的有限元变形计算。
根据各个岩层发生的形变进行比较分析，在边界上满足接触面位移条件，即0＞wi, 1(x，y)＞wi+1, 2(x，y)。若该式成立，则两个岩层可能出现离层，继续视为两个不同表面，表面为自由表面边界条件，亦即表面力为0。若相邻两层接触面位移满足条件wi, 1(x，y)≤wi+1, 2(x，y)＜0，表明i层下表面发生的形变大于i+1层上表面，则粘合分裂的两个点，其坐标均为(x，y)，在下一次迭代过程中两个点则视为一个点。考虑最底层的下边界垂直方向位移wn，2(x，y)，如果满足wn，2(x，y)≤-D，则将该点的垂直方向位移w设为第一类边界条件，否则考虑为自由表面边界条件，如此反复迭代直至收敛，则接触边界条件得到正确处理。
2. 三维离层模拟算例
本文以济宁二号煤矿西部采区（十一采区）存在的覆岩离层问题为例进行定量化模拟研究。利用地震成像理论确定离层水平空间范围，确定计算区域和适当的边界条件（见图 1），利用有限元预处理软件建立二维模型（见图 2）。
利用生成的二维网格沿第三方向延展得到三维网格，同时记录两个相邻地层的网格对应关系，以便计算时处理接触边界。本文采用网格粘合技术（见图 3）处理接触边界条件，首先将相互连接的两个地层的同一个点分裂为2个不同的点，进行有限元计算，然后再根据计算结果接触情况进行网格粘合，图 4为研究区域的三维网格示意图。
根据采区部分钻孔的测井数据（杨氏模量、密度和泊松比等），建立岩层空间分布信息以及不同岩层的岩性，再查阅相关资料得到不同岩层的力学模拟参数（见表 1、图 5）。
表 1 采区覆岩层空间分布及力学参数Table 1. Spatial distribution and mechanical parameters of each layer
本开采工作面采掘3下煤层，采高5 m左右。利用上述三维网格以及物理参数，在大规模计算集群上进行了数值模拟。图 6列出了多次迭代收敛过程。
图 6a—图 6e代表每次迭代过程，从收敛过程可以看出，尽管第1次迭代时计算结果有上层穿插下层的情况，但是随着迭代步数增加，在第5次迭代完成后，收敛到最终的变形情况，显示出3个不同离层位置（见图 6e），分别位于累积深度-657.5 m，-707.55 m，-766.7 m处，相应的最大离层高度分别为0.2144 m、0.5261 m和0.1809 m（见表 2）。
表 2 离层数据Table 2. tatistics of separation layers
应该注意，模拟过程中应用了弹性力学方程，没有能力模拟采煤过程中冒落带形成的各种物理过程和获得其几何参数，因此离层3参数具有较大的不确定性，它可能已经位于冒落带和裂隙带范围（见图 6e）。
计算结果进行可视化，得到反映离层形态和体积的直观图，直观显示了离层形状特性以及水平位置分布。本文给出了模拟出的离层1和离层2的变形俯视图和变形三维图（见图 7、图 8），从图中可以看到离层位于采掘面后方，先变大，后闭合的过程（为了绘图的需要将x，y坐标归一化到[0, 1]）。
为了显示离层，将离层2高度放大20倍，可以清楚地看到其几何特征（见图 9）。分析表明，该形态符合离层结构力学特征。
根据唐山矿覆岩离层井下视频观测，覆岩离层的最大高度一般为几厘米到几十厘米[7]。济宁二号矿与唐山矿均位于华北地层区，开采层位均为P—C煤层，采掘深度唐山矿500 m左右，济宁二号矿700 m左右，地质和几何空间尺度等条件相似，因此上述模拟得到的覆岩离层几何参数具有可信度，在煤矿安全生产方面也具有重要参考价值。
3. 结论
在覆岩离层的接触问题中，垂直方向的位移较大，水平方向位移相对较小，利用这一特点，首先考虑相邻层只在重力作用下的变形，然后判断垂直位移的大小，如果满足接触条件，则采用网格粘合技术将其粘合视为一个点，很好地解决了这一类接触问题。通过二维网格的三维延展，实现了大规模接触问题的网格生成问题。首次使用真正三维模型进行模拟，得到了离层体积量、最大离层高度等参数。
-
表 1 采区覆岩层空间分布及力学参数
Table 1. Spatial distribution and mechanical parameters of each layer
表 2 离层数据
Table 2. tatistics of separation layers
-
[1] 何国清, 杨伦, 凌赓娣.矿山开采沉陷学[M].徐州:中国矿业大学出版社, 1991HE Guo-qing, YANG Lun, LING Geng-di. Mining subsidence[M]. Xuzhou:China University of Mining and Technology Press, 1991. [2] 杨伦, 于广明. 采矿下沉的再认识[C]//第七届国际矿测学术会议文集. 1987 http://www.oalib.com/paper/4421555YANG Lun, YU Guang-ming. Re-recognization of mining submersion[C]//Collected works of the 7th International Symposium on Mine Surveying. 1987. http://www.oalib.com/paper/4421555 [3] 宋杨, 宋振骐.采场支撑压力显现规律与上覆岩层运动的关系[J].煤炭学报, 1991, 9(1):47~56 http://cdmd.cnki.com.cn/Article/CDMD-11413-1015444972.htmSONG Yang, SONG Zhen-qi. Relationship between abutment pressure distribution laws of working faces and movement of overburden rocks[J]. Journal of China Coal Society, 1991, 9(1):47~56. http://cdmd.cnki.com.cn/Article/CDMD-11413-1015444972.htm [4] 苏仲杰, 于光明, 杨伦.覆岩离层变形力学模型及应用[J].岩土工程学报, 2002, 24(6):778~781 http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200206023.htmSU Zhong-jie, YU Guang-ming, YANG Lun. Application of mechanical model to deformation of covered rock seperation strata[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(6):778~781. http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200206023.htm [5] 梁运培.采场覆岩移动的组合岩梁理论[J].地下空间, 2001, 21(5):341~345 http://www.cnki.com.cn/Article/CJFDTOTAL-BASE2001S1000.htmLIANG Yun-pei. Study on the composed rock timber theory of overlying strata movement[J]. Underground Space, 2001, 21(5):341~345. http://www.cnki.com.cn/Article/CJFDTOTAL-BASE2001S1000.htm [6] 张健全, 廖国华.覆岩离层产生的机理及离层计算方法的探讨[J].地下空间, 2001, 21(增刊):407~417 doi: 10.3969/j.issn.1673-0836.2001.01.005ZHANG Jian-quan, LIAO Guo-hua. Investigation on formation mechanism of separated layer of rock vovering and calculation method of separated layer[J]. Underground Space, 2001, 21(Supp.):407~417. doi: 10.3969/j.issn.1673-0836.2001.01.005 [7] 贾君莹, 钟亚平, 王素华, 等.铁路下特厚煤层综放开采覆岩离层规律实测研究[J].矿山测量, 2008, (4):59~62 http://www.cnki.com.cn/Article/CJFDTOTAL-KSCL200804022.htmJIA Jun-ying, ZHONG Ya-ping, WANG Su-hua, et al. Observation research on overburden movement and bed-separation laws by fully mechanized top coal caving in extremely thick coal seam beneath railway[J]. Mine Surveying, 2008, (4):59~62. http://www.cnki.com.cn/Article/CJFDTOTAL-KSCL200804022.htm -
下载:
下载:
计量
- 文章访问数: 253
- HTML全文浏览量: 142
- PDF下载量: 7
- 被引次数: 0
