一种煤层开采覆岩沉陷三维半解析分析方法.pdf

本发明公开了一种煤层开采覆岩沉陷三维半解析分析方法，利用半解析数值分析方法来解决层状岩体矿山开采覆岩沉陷的三维分析。具体的，该分析方法包括如下步骤：(1)首先建立层状覆岩沉陷三维半解析模型；(2)然后利用一些层、柱单元的标准解析解，并结合有限元的离散化分析，完成条单元、柱单元和层单元之间的转化，实现对模型各部分的位移、单元刚度矩阵、单元载荷矩阵和应力的计算。本发明方法不仅克服了纯解析的理论分析在数学上的困难及应用上的局限性，又大大降低了基于全离散原理的纯数值方法的计算工作量。

权利要求书
1.  一种煤层开采覆岩沉陷三维半解析分析方法，其特征在于，包括如下步骤：
s1、建立层状覆岩沉陷三维半解析模型，其中：
模型上部进行单元层分割，表土为水平层，视为横观各向同性材料，为水平条单元，采 用二维解析的层单元，部分的局部坐标与整体坐标相重合；
模型中部进行三角块分割，呈三角状的倾斜块，离散为三棱柱单元；
模型下部进行单元层分割，包括开挖在内的呈倾斜状的各岩层，为倾斜条单元，采用二 维解析的层单元，部分的局部坐标与整体坐标有一倾角，在局部坐标系下二者的计算公式有 相同的形式；
s2、利用一些层、柱单元的标准解析解，并结合有限元的离散化分析，完成条单元、柱 单元和层单元之间的转化，实现对层状覆岩沉陷三维半解析模型各部分的位移、单元刚度矩 阵、单元载荷矩阵和应力的计算。

2.  根据权利要求1所述的煤层开采覆岩沉陷三维半解析分析方法，其特征在于，所述步 骤s2中，位移、单元刚度矩阵、单元载荷矩阵和应力的计算公式如下：
(1)使用模型层单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p Σ n = 1 q [ N ] m | N yz { δ } mn | - - - ( 1 ) ]]>
其中，
[ N ] mn = cos mπx s 0 0 cos nπy b 0 0 0 sin mπx a 0 0 sin nπy b 0 0 0 sin mπx a 0 0 sin nπy b , ]]>
[ N yz ] = z 2 - z h [ N ] n z - z 1 h [ N ] n , ]]>
[ N ] n = cos nπy b 0 0 0 sin nπy b 0 0 0 sin nπy b , ]]>
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 w 2 T , ]]>
a、b、h分别为层单元的长、宽、厚，u、v、w为边界位移；z1，z2为边界点坐标值；p， q为所划分层单元的数量；
棱柱单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p [ N ] m N yz { δ } mn - - - ( 2 ) ]]>
其中，
[Nyz]＝[L1[I] L2[I] L3[I]]，
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 u 3 v 3 w 3 T , ]]>
L1、L2、L3为三角形内的面积坐标，[I]为三阶单位矩阵；
(2)使用模型层单元刚度矩阵计算公式，如下：
{ ϵ } = Σ m = 1 p Σ n = 1 q [ B ] mn { δ } mn - - - ( 3 ) ]]>
其中，
[ B ] = ∂ u ∂ x ∂ v ∂ y ∂ w ∂ z ∂ u ∂ x + ∂ v ∂ y ∂ v ∂ y + ∂ w ∂ z ∂ w ∂ z + ∂ u ∂ x ; [ B ] mn = [ B m B n ] ; ]]>
棱柱单元刚度矩阵计算计算公式，如下：
{ ϵ } = Σ m = 1 p [ B ] mn { δ } mn [ K ] m , k = ∫ ∫ ∫ v [ B ] m T [ D ] [ B ] k dxdydz - - - ( 4 ) ]]>
其中， [ D ] = v 1 E 1 E 2 0 0 0 E 1 v 1 E 2 0 0 0 E 2 E 2 v 2 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 2 ]]>为横观各向同性或各向同性材料的弹性矩阵，E 为岩石弹性模量，G为剪切弹性模量；
(3)使用模型层单元载荷列阵计算公式，如下：
{ F } mn = ∫ z 1 z 2 Σ j = 1 jy ∫ y qj y zj Σ i = 1 x zi ∫ x qi x zi [ N ] yz T [ N ] m T 0 0 γ dxdydt - - - ( 5 ) ]]>
棱柱单元载荷列阵计算公式，如下：
{ F } = ∫ ∫ ∫ v [ N yz ] T [ N ] m T 0 0 γ dxdydz - - - ( 6 ) ]]>
(4)使用模型层单元应力计算公式，如下：
{σ}＝ΣΣ[D][B]mn[T]T{δe}mn    (7)
棱柱单元应力计算公式，如下：
{σ}＝ΣΣ[D][B]m{δe}mn(8)
其中，[T]T＝{δ}T为转置矩阵，6阶矩阵[D]是横观各向同性或各向同性材料的弹性矩阵。

说明书一种煤层开采覆岩沉陷三维半解析分析方法
技术领域
本发明属于煤矿开采覆岩运动地表变形预计领域，特别涉及一种适用于计算覆岩沉陷三 维演变规律的煤层开采覆岩沉陷三维半解析分析方法。
背景技术
开采覆岩沉陷研究发展100多年来，各种理论分析方法复杂多样。其中，波兰学者 Litwinszyn基于开采引起的岩层和地表移动规律与作为随机介质的非连续颗粒体介质模型所 描述的规律在宏观上相似原理，于1954年提出的随机介质理论方法，由于其数理概念明确， 实用性强，得到了广泛的应用，我国在50年代末由刘宝琛、廖国华等学者引入此方法，并逐 渐发展成为目前在国内外应用最广泛，理论与实践上最为成熟的概率积分法。但由于该方法 未能给出材料介质位移和应力响应方程，使得这种主要用于地表移动变形分析的方法对于开 采覆岩沉陷问题的岩体内部应力、应变计算产生了困难。
半解析数值方法是指在数值分析中采用与引入部分解析解或解函数，而计算得到的仍是 一系列离散化数值结果的分析方法。它通过应用解析手段来帮助减少数值计算的维数，即由 人完成诸如一维解函数、基本解、通解等较为规则、简单情况下的理论分析，而同时又通过 数值离散手段用计算机来处理降维后的低维离散化方程的求解，弥补人为所取解析函数的不 足，以逼近真实解。基于上述原因，本发明所要解决的技术问题是，如何利用半解析数值分 析方法来解决层状岩体矿山开采覆岩沉陷的三维分析。
发明内容
本发明的目的在于提出一种煤层开采覆岩沉陷三维半解析分析方法，该分析方法不仅能 够克服纯解析的理论分析在数学上的困难及应用上的局限性，又能够大大降低基于全离散原 理的纯数值方法的计算工作量。
为了实现上述目的，本发明采用如下技术方案：
一种煤层开采覆岩沉陷三维半解析分析方法，包括如下步骤：
s1、建立层状覆岩沉陷三维半解析模型，其中：
模型上部进行单元层分割，表土为水平层，视为横观各向同性材料，为水平条单元，采 用二维解析的层单元，部分的局部坐标与整体坐标相重合；
模型中部进行三角块分割，呈三角状的倾斜块，离散为三棱柱单元；
模型下部进行单元层分割，包括开挖在内的呈倾斜状的各岩层，为倾斜条单元，采用二 维解析的层单元，部分的局部坐标与整体坐标有一倾角，在局部坐标系下二者的计算公式有 相同的形式；
s2、利用一些层、柱单元的标准解析解，并结合有限元的离散化分析，完成条单元、柱 单元和层单元之间的转化，实现对层状覆岩沉陷三维半解析模型各部分的位移、单元刚度矩 阵、单元载荷矩阵和应力的计算。
在上述步骤s2中，位移、单元刚度矩阵、单元载荷矩阵和应力的计算公式如下：
(1)使用模型层单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p Σ n = 1 q [ N ] m | N yz { δ } mn | - - - ( 1 ) ]]>
其中，
[ N ] mn = cos mπx a 0 0 cos nπy b 0 0 0 sin mπx a 0 0 sin nπy b 0 0 0 sin mπx a 0 0 sin nπy b , ]]>
[ N yz ] = z 2 - z h [ N ] n z - z 1 h [ N ] n , ]]>
[ N ] n = cos nπy b 0 0 0 sin nπy b 0 0 0 sin nπy b , ]]>
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 w 2 T , ]]>
a、b、h分别为层单元的长、宽、厚，u、v、w为边界位移；z1，z2为边界点坐标值；p， q为所划分层单元的数量；
棱柱单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p [ N ] m [ N yz ] { δ } mn - - - ( 2 ) ]]>
其中，
[Nyz]＝[L1[I] L2[I] L3[I]]，
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 u 3 v 3 w 3 T , ]]>
L1、L2、L3为三角形内的面积坐标，[I]为三阶单位矩阵；
(2)使用模型层单元刚度矩阵计算公式，如下：
{ ϵ } = Σ m = 1 p Σ n = 1 q [ B ] mn { δ } mn - - - ( 3 ) ]]>
其中，
[ B ] = ∂ u ∂ x ∂ v ∂ y ∂ w ∂ z ∂ u ∂ x + ∂ v ∂ y ∂ v ∂ y + ∂ w ∂ z ∂ w ∂ z + ∂ u ∂ x [ B ] mn = [ B m B n ] ; ]]>
棱柱单元刚度矩阵计算计算公式，如下：
{ ϵ } = Σ m = 1 p [ B ] mn { δ } mn [ K ] m , k = ∫ ∫ ∫ v [ B ] m T [ D ] [ B ] k dxdydz - - - ( 4 ) ]]>
其中， [ D ] = v 1 E 1 E 2 0 0 0 E 1 v 1 E 2 0 0 0 E 2 E 2 v 2 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 2 ]]>为横观各向同性或各向同性材料的弹性矩阵，E 为岩石弹性模量，G为剪切弹性模量；
(3)使用模型层单元载荷列阵计算公式，如下：
{ F } mn = ∫ z 1 z 2 Σ j = 1 jy ∫ y qj y zj Σ i = 1 ix ∫ x qi x zi [ N ] yz T [ N ] m T 0 0 γ dxdydz - - - ( 5 ) ]]>
棱柱单元载荷列阵计算公式，如下：
{ F } = ∫ ∫ ∫ v [ N yz ] T [ N ] m T 0 0 γ dxdydz - - - ( 6 ) ]]>
(4)使用模型层单元应力计算公式，如下：
{σ}＝∑∑[D][B]mn[T]T{δe}mn    (7)
棱柱单元应力计算公式，如下：
{σ}＝∑∑[D][B]m{δe}mn     (8)
其中，[T]T＝{δ}T为转置矩阵，6阶矩阵[D]是横观各向同性或各向同性材料的弹性矩阵， 可以从有关书中查得。
与现有技术相比，本发明具有如下优点：
本发明方法，首先建立层状覆岩沉陷三维半解析模型；然后利用一些层、柱单元的标准 解析解，并结合有限元的离散化分析，完成条单元、柱单元和层单元之间的转化，实现对层 状覆岩沉陷三维半解析模型各部分的位移、单元刚度矩阵、单元载荷矩阵和应力的计算。本 发明方法不仅克服了纯解析的理论分析在数学上的困难及应用上的局限性，又大大降低了基 于全离散原理的纯数值方法的计算工作量，具有广阔的推广应用性。
附图说明
图1为本发明中层状覆岩沉陷三维半解析模型图。
图2为半解析分析法中典型的层单元图。
图3为本发明具体实例中层状覆岩沉陷三维半解析模型图。
图4为本发明具体实例中沿x轴方向计算下沉曲线图。
图5为本发明具体实例中地表下沉等值线图。
图6为本发明具体实例中地表下不同深度水平应力曲线图。
图7为本发明具体实例中走向主剖面覆岩下沉等值线图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
结合图1所示，一种煤层开采覆岩沉陷三维半解析分析方法，包括如下步骤：
s1、建立层状覆岩沉陷三维半解析模型，其中：
模型上部进行单元层分割，表土为水平层，视为横观各向同性材料，为水平条单元，采 用二维解析的层单元，部分的局部坐标与整体坐标相重合；
模型中部进行三角块分割，呈三角状的倾斜块，离散为三棱柱单元；
模型下部进行单元层分割，包括开挖在内的呈倾斜状的各岩层，为倾斜条单元，采用二 维解析的层单元，部分的局部坐标与整体坐标有一倾角，在局部坐标系下二者的计算公式有 相同的形式。
s2、利用一些层、柱单元的标准解析解，并结合有限元的离散化分析，完成条单元、柱 单元和层单元之间的转化，实现对层状覆岩沉陷三维半解析模型各部分的位移、单元刚度矩 阵、单元载荷矩阵和应力的计算。
具体的，在步骤s2中，位移、单元刚度矩阵、单元载荷矩阵和应力的计算公式如下：
(1)使用模型层单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p Σ n = 1 q [ N ] m | N yz { δ } mn | - - - ( 1 ) ]]>
其中，
[ N ] mn = cos mπx a 0 0 cos nπy b 0 0 0 sin mπx a 0 0 sin nπy b 0 0 0 sin mπx a 0 0 sin nπy b , ]]>
[ N yz ] = z 2 - z h [ N ] n z - z 1 h [ N ] n , ]]>
[ N ] n = cos nπy b 0 0 0 sin nπy b 0 0 0 sin nπy b , ]]>
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 w 2 T , ]]>
a、b、h分别为层单元的长、宽、厚，u、v、w为边界位移；z1，z2为边界点坐标值；p， q为所划分层单元的数量。
棱柱单元位移计算公式，如下：
{ f } = u v w = Σ m = 1 p [ N ] m [ N yz ] { δ } mn - - - ( 2 ) ]]>
其中，
[Nyz]＝[L1[I] L2[I] L3[I]]，
{ δ } mn T = u 1 v 1 w 1 u 2 v 2 u 3 v 3 w 3 T , ]]>
L1、L2、L3为三角形内的面积坐标，[I]为三阶单位矩阵；
(2)使用模型层单元刚度矩阵计算公式，如下：
{ ϵ } = Σ m = 1 p Σ n = 1 q [ B ] mn { δ } mn - - - ( 3 ) ]]>
其中，
[ B ] = ∂ u ∂ x ∂ v ∂ y ∂ w ∂ z ∂ u ∂ x + ∂ v ∂ y ∂ v ∂ y + ∂ w ∂ z ∂ w ∂ z + ∂ u ∂ x [ B ] mn = [ B m B n ] ; ]]>
棱柱单元刚度矩阵计算计算公式，如下：
{ ϵ } = Σ m = 1 p [ B ] mn { δ } mn [ K ] m , k = ∫ ∫ ∫ v [ B ] m T [ D ] [ B ] k dxdydz - - - ( 4 ) ]]>
其中， [ D ] = v 1 E 1 E 2 0 0 0 E 1 v 1 E 2 0 0 0 E 2 E 2 v 2 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 1 0 0 0 0 0 0 G 2 ]]>为横观各向同性或各向同性材料的弹性矩阵，E 为岩石弹性模量，G为剪切弹性模量。
(3)使用模型层单元载荷列阵计算公式，如下：
{ F } mn = ∫ z 1 z 2 Σ j = 1 jy ∫ y qj y zj Σ i = 1 ix ∫ x qi x zi [ N ] yz T [ N ] m T 0 0 γ dxdydz - - - ( 5 ) ]]>
棱柱单元载荷列阵计算公式，如下：
{ F } = ∫ ∫ ∫ v [ N yz ] T [ N ] m T 0 0 γ dxdydz - - - ( 6 ) ]]>
(4)使用模型层单元应力计算公式，如下：
{σ}＝ΣΣ[D][B]mn[T]T{δe}mn    (7)
棱柱单元应力计算公式，如下：
{σ}＝ΣΣ[D][B]m{δe}mn     (8)
其中，[T]T＝{δ}T为转置矩阵，为6阶矩阵[D]是横观各向同性或各向同性材料的弹性矩 阵，可以从有关书中查得。
为进一步理解本发明，下面给出一个具体实例进行说明：
在具体实例中，煤矿采区走向长度为600m，斜长为320m，煤层倾角26°，煤厚2m，利 用完全垮落法处理采空区。
计算范围内共有7种不同的岩层组成，共划分为12层，建立层状覆岩沉陷三维半解析模 型，如图3所示，地表两层由属于横观各向同行的第一种材料组成，其余各岩层的物理参数 如表1所示，各层的厚度如表2所示。
表1 覆岩层物理力学参数
岩层 弹性模量/(104×MPa) 密度/(kg/m3) 泊松比u 上部砾岩 4.48 2680 0.2 下部砾岩 4.92 2680 0.2 红层 3.21 2650 0.22 泥岩 1.21 2650 0.22 砂岩 3.88 2660 0.23 煤 1.1 1400 0.27 破碎带 0.098 2000 0.1   K1＝10.8MPa，K2＝10.8MPa    
表2 计算层厚
层号 1 2 3 4 5 6 7 8 9 10 11 12 厚度/m 50 50 - - 40 30 20 10 10 10 20 30
计算所得沿x轴方向计算下沉曲线如图4所示，地表下沉等值线图如图5所示。为分析 地面斑裂机理及止裂深度，计算走向剖面地表以下不同深度水平应力，曲线如图6所示，可 得到距采空区外侧边界约185m处地表下20m左右岩体附加应力已经超过了砾岩的抗压强度 (1.85MPa)，砾岩产生拉断裂缝，形成地面斑裂，这为分析地面斑裂机理提供了依据。实例 中走向主剖面覆岩下沉等值线如图7所示。
本发明方法不仅克服了纯解析的理论分析在数学上的困难及应用上的局限性，又大大降 低了基于全离散原理的纯数值方法的计算工作量，具有广阔的推广应用性。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说 明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变 形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。