基于Schwarz交替法求解分析浅埋近距离煤层 开采覆岩结构演化规律

李 磊，李寿君

( 1. 煤炭科学研究总院，北京 100013；2. 天地科技股份有限公司 开采设计事业部，北京 100013；3. 中煤科工开采研究院有限公司，北京 100013；4. 扎赉诺尔煤业有限责任公司，内蒙古 满洲里 021412 )

随着我国中东部地区煤炭资源的逐渐枯竭，国家煤炭资源开采重心逐渐向西部转移。同时，西部煤田一些矿井的上部煤层逐渐开采结束，矿井进入下部煤层开采阶段，比如隆德煤矿、活鸡兔矿井和布尔台煤矿等，该类型矿井煤层埋藏较浅，上下煤层间距较小，覆盖层为含水砂层。该类型矿井在开采过程中，不仅面临着第一主煤层采空区积水的威胁，而且可能面临下部煤层重复开采引起原导水裂缝带二次发育造成的矿井灾害问题，例如波及含水层造成工作面涌水溃砂等。

基于板或梁的力学分析模型，国内外专家学者已经对工作面开采覆岩结构演化规律做了大量的科学研究。1867年，德国学者舒里兹认为采空区上方岩层类似于梁结构，提出采场上覆岩层结构的“悬梁假说”[1]；苏联学者里特捷尔将采空区上方岩层类比于梁或板，提出了“压力拱假说”[2-4]；钱鸣高院士于1981年提出了“砌体梁”结构模型，研究了覆岩的“O-X”型破断规律[5-6]；1988年，宋振骐院士提出了“传递岩梁”的理论模型[7]；20世纪90年代，钱鸣高院士[8-9]等提出了岩层控制的“关键层理论”，认为岩层中存在控制岩层移动变化规律的关键岩层，并系统分析了采动损害现象和防止采动灾害；赵德深[10]提出了上覆岩层运动的“拱 板式”和“板式”两种力学平衡模型，研究了开采空间在覆岩中的传播规律；姜福兴[11]等将基本顶结构划分为“类拱式”、“拱梁式”和“梁式”；黄庆 享[12-13]提出了“拱梁”结构数学模型，并研究了 厚砂土层产生破裂的判据；窦林名[14]等分析了工作面岩层的边界差异，将覆岩运动结构划分为OX，F和T型结构；吴立新[15]等在研究大面积开采岩层 变形破坏时提出了覆岩托板控制理论；闫少宏等提出了顶板短悬臂梁-铰接岩梁结构，并研究了煤 壁前方破坏面积和支承压力规律[16]；尹希文[17]提出了浅埋煤层工作面覆岩“切落体”结构模型，并 研究了切落体结构的稳定条件及失稳类型；任艳 芳[18]提出了浅埋深工作面覆岩中“悬臂梁-铰接 岩梁”组合岩层结构，认为“铰接岩梁”结构在 煤壁处出现滑落失稳是造成浅埋长壁工作面顶 板沿煤壁切落、形成压架事故的根本原因。目前，大部分浅埋近距离煤层重复开采的研究是基于 板或梁的力学分析模型，但该类力学模型更适合分析单煤层开采情况。浅埋近距离煤层重复开采的大部分研究成果主要集中在采场矿压规律方面，而对于重复开采垮落带和导水裂缝带二次发育高度的研究并不多[18-24]。因此，本文以浅埋近距离煤层开采工作面为研究对象，提出一种更适合于近距 离煤层开采扰动的“无限平面双矩形孔口”理论分析模型，通过理论分析和相似材料模拟试验研究近距离煤层采煤工作面的覆岩结构演化规律，为 浅埋近距离煤层安全开采提供理论基础和技术保障，从而避免工作面充水溃砂等重大灾害事故的发生。

1 近距离煤层重复采动应力场理论解析

近距离煤层开采过程中，沿工作面推进方向垂直剖开形成无限个平面，每个平面的边界处受到相同的分布荷载，因此，此类问题可以看作平面应变问题。在这些平面中，上下两个采煤工作面类似于在一个半无限面内开2个矩形孔，这属于弹性力学的双连域问题—“无限平面双矩形孔口问题”。

本节采用弹性力学和复变函数理论构建近距离煤层开采扰动的“无限平面双矩形孔口”理论模型，并对其有效性进行验证。通过该理论模型分析近距离煤层开采覆岩应力场分布情况。

1.1 多连域问题的复变函数解法

对于弹性力学的双连域问题，常规弹性力学解法难以求解，因此需要复变函数法进行变换求解。复变函数解法是求出双调和方程的通解，并在不同的应力边界条件下寻求特解的过程。在复变函数解法中，Airy应力函数U可用复数z的2个解析函数φ1( z )和θ1( z )来表示[25-27]。但是，直接求解函数φ1( z )和ψ1( z )仍然比较困难。因此，本节利用复变函数保角变换法[26]，把z平面上具有复杂形状边界( 如矩形、椭圆形 )的一个区域变换为另一复平面ξ 上具有简单形状的边界( 如单位圆内或单位圆外 )，如图1所示，然后把在简单边界下求得的解再变换回去，从而求得原问题的解。

图1 复平面内边界形状变换示意 Fig. 1 Schematic diagram of boundary shape transformation in complex plane

令ψ1( z )=θ1′( z )，以z为自变量的2个解析函数φ1( z )和ψ1( z )通过映射函数 z = ω （ξ ）变成以ξ 为自变量的函数。变换关系为

利用上述变换关系，应力分量的复变函数可改写为

位移分量的复变函数可表示为

则应力边界条件为

式中，E为弹性模量；μ为泊松比；Xij和 Yij分别为边界上面力的2个分量。

对于多孔问题，多连通域不可能单值地保角映射到单连通域，无法直接获得其映射函数。因此，本文将多连通域转化成一系列单连通域后，采用Schwarz交替法进行逐一求解。

1.2 近距离煤层开采扰动覆岩应力场理论解析方法

将近距离煤层开采简化为平面应变问题，并抽象为无限平面双矩形孔口问题，如图2所示。为便于分析，设上下煤层工作面尺寸相同，均为长2a，宽2b，中间岩层高为h，水平和垂直原岩应力分别为p1和p2，在上下工作面分别建立直角坐标系x1O1y1和x2O2y2。z1和z2分别表示x1O1y1和x2O2y2坐标系下的坐标。

图2 近距离煤层开采扰动平面模型 Fig. 2 Plane model of near-range coal mining disturbance

要将多连通域问题转化成一系列单连通域问题，需要用到Schwarz交替法，然后针对每个单连通域问题采用柯西积分解法求解。具体的求解流程为：① 只开采上煤层，利用柯西积分法求得只有上工作面1的应力场；② 在上工作面1的应力场作用下，下工作面2的边界必然产生多余的面力( 洞边不为零的面力 )，而洞边是要满足零面力条件的，因此，在下工作面2边界处施加与多余面力大小相等且方向相反的反面力，模拟下工作面2，此时只有下工作面2没有上工作面1，可求得只有下工作面2的应力场；③ 在上工作面1边界位置也会产生多余的面力；④ 反复进行上述模拟开挖过程，直到两洞洞边均满足零面力条件，将每次的迭代求解结果叠加，即为最终的求解结果。

1.2.1 上煤层回采后工作面复应力函数

把z平面上矩形孔洞外部区域变换成ξ平面上单位圆外部区域，根据Schwarz-Christoffel积分公式，其映射函数为

式中，n为正奇数；R为反映洞形大小的正实数，

式中，k的值与工作面的宽长比b/a有关，需进行试算，控制方程为式( 8 )。

k值应使得ε趋近于0。R的值亦与工作面尺寸有关，用式( 9 )表示。

单位圆外包括无穷远点在内的2个解析函数φ11( ξ1)和ψ11( ξ1)通式为

式 中，B =（ p1+ p2）4，B′ = （p2- p1）2，φ11( ξ1)和ψ11( ξ1)为ξ1平面单位圆外整个区域的解析函数，可表示为

在地应力p1和p2作用下只存在工作面1时，φ11( ξ1)和ψ11( ξ1)为

1.2.2 上煤层回采引起的下煤层工作面边界面力

x1O1y1坐标面简称为z1平面，x2O2y2坐标面简称为z2平面。假定φ11( ξ1)和ψ11( ξ1)在ξ2平面下的形式为φ12( ξ2)和ψ12( ξ2)，φ22( ξ2)和ψ22( ξ2)在ξ1平面下的形式为φ21( ξ1)和ψ21( ξ1)。工作面2边界点在ξ2平面下坐标为σ2，在ξ1平面下坐标为γ1。则φ11( ξ1)和ψ11( ξ1)在工作面2边界引起的多余面力f12( σ2)的表达式为

由式( 13 )求得的多余面力 f12（σ2）是下工作面2边界面力值，其分布的显式表达可利用复数级数式( 14 )进行逼近，Dk是复数逼近式中 σ2的各幂次项系数。

如图3所示，由于上煤层回采引起的下煤层工作面边界多余面力的求解步骤为：① σ2通过映射变换t2=ω( σ2)得到对应z2平面下的点t2；② t2通过坐标平移T1=t2+c得到z1平面下的坐标T1；③ T1通过映射变换γ1=ω-1( T1)得到ξ1平面下的坐标γ1；④ 将γ1代入式( 13 )即可求得洞2周边的多余面力值f12( σ2)；⑤ 对多余面力f12( σ2)进行复数级数逼近。

图3 坐标点变换示意 Fig. 3 Diagram of coordinate point transformation

1.2.3 反面力作用下只存在下煤层时的复应力函数

1.2.4 应力分量迭代求解

ξ1按照式( 10 )中的坐标变换得到ξ2与ξ1的关系，ξ2=f( ξ1)。将该式代入φ22( ξ2)和ψ22( ξ2)得到φ21( ξ1)和ψ21( ξ1)。因此，Schwarz交替法的第1次迭代结果为

第1次迭代求解后，下工作面2边界满足零面力条件，而上工作面1在φ21( ξ1)和ψ21( ξ1)作用下仍可能存在多余面力。再次按照1.2.2节和1.2.3节，所述，进行反复迭代，直到2工作面的多余面力小到一定数值，从而最终将所有迭代结果叠加得到φ1( ξ1)和ψ1( ξ1)。

1.2.5 应力分量

z平面上的应力分量组合式，经变换以ξ为自变量的复变函数表达式，利用式( 3 )可得：

1.3 理论模型有效性分析

通过对比分析双矩形孔口问题的理论计算和数值模拟结果，验证该理论方法的可靠性。模型的相关参数见表1。

表1 参数设置 Table 1 Parameters setting

1.3.1 映射函数确定的合理性

利用式( 8 )，当k=0.069时，ε ≈0。此时的映射函数表达式为

根据该映射函数表达式分析可得，映射后的矩形除了在4个边角处微翘外，其余各处与原矩形的形状尺寸基本一致，说明所选映射函数是合理的。

1.3.2 零面力条件验证

无支护情况下，工作面边界应满足零面力条件，Schwarz交替法求解本身就是一个逐步消除多余面力的过程。经反复迭代后，在反面力作用下工作面2周边正应力和剪应力分布如图4所示。从图4可以看出，上煤层1回采后下煤层2边界多余面力分量虽然未全部为0，但这些面力分量均在零值附近小幅度波动，这说明Schwarz交替求解法基本满足零面力条件，该理论方法满足工程计算的要求。

图4 上煤层1回采后下煤层2边界处多余的面力 Fig. 4 Excess surface force at boundary of coal seam 2 after mining in coal seam 1

1.3.3 应力理论解与数值模拟结果比较

如图5所示，本文采用ABAQUS软件建立平面双矩形孔口计算模型，模型尺寸为400 m×400 m。上下工作面开挖尺寸均为100 m×4 m，中间岩层厚度为20 m。单元属性为弹性，单元类型为CPE4R。模型的应力边界条件为：上下边界施加3 MPa均布载荷，左右两侧施加2 MPa均布载荷。然后在上工作面的右侧和上侧选取2条应力分量监测线，从而对比分析数值模拟和理论计算2种方法获得的应力分量分布规律的差异性。

图5 有限元模型及监测点 Fig. 5 Finite element model and monitoring method

图6是2条监测线上的数值模拟和理论计算的应力分量分布规律。横坐标是距洞边的距离。对于测线1，σy的理论值和模拟值变化趋势均逐渐减小，但在工作面及其附近位置，理论值小于模拟值。σx的变化趋势是先增大后减小，模拟值的变化幅度明显。对于测线2，σx的变化趋势仍然是先增大后减小。σy的趋势则为逐渐增大，理论值普遍略大于模拟值。综上所述，理论计算和数值模拟的σx和σy的变化趋势和数值大致相同，相似程度高，说明在弹性条件下该理论方法求解双矩形孔口的应力分布问题的可靠性较好。

图6 2条测线上应力分布规律对比 Fig. 6 Comparison of stress distribution laws on two measuring lines

2 单一煤层开采扰动覆岩应力场及破坏区分析

通过理论解析方法计算近距离煤层开采覆岩应力场分布后，根据不同位置覆岩所受的应力状态分析覆岩破坏范围及其演化过程。判断标准为：垮落带为采空区上方拉张破坏区及局部拉张区，裂缝带为垮落带上部的拉裂隙区。本节将垮落带和裂缝带作为覆岩破坏范围进行分析。

由于弹性力学理论求解近距离煤层开采扰动覆岩应力场问题的复杂性，求解不同力学性质和垂直各向异性的层状煤岩层力学模型变得非常困难，甚至无法得到理论解析解。因此，将实际的层状煤岩层地质模型简化为均匀各向同性的弹性体模型。单煤层开采理论计算的相关力学参数见表2，分析其应力场分布及其破坏范围。

基于Mohr-Coulomb准则的推导式( 19 )，建立覆岩塑性极限平衡区计算方法，并以此作为判断覆岩的破坏区。当 0f ＜ 时，围岩应力状态处于抗剪强度包线下方，此时围岩未破坏；当 =0f 时，围岩应力状态在包线上，处于极限平衡状态；当 0f ＞ 时，围岩应力状态超出包线，围岩发生破坏。

表2 理论模型的参数设置 Table 2 Parameter setting of mechanical model

2.1 单一煤层开采覆岩应力场分布规律

图7为单一煤层开采覆岩应力分量σx和σy的分布云图。在图7( a )中，覆岩σx整体分布大致呈“X”状，顶底板处于拉应力状态，呈拱形分布，随着远离顶底板，σx由拉应力变为压应力。左右两侧应力呈椭圆状分布，存在压应力集中区。在图7( b )中，工作面附近σy应力场分为3个区域：顶底板钟型分布区、左右两侧耳状分布区、外部原岩应力分布区。顶底板、左右两侧均为压应力，左右两侧耳状分布区压应力值要普遍大于顶底板钟型分布区，越靠近工作面的左右两侧，压应力值越大。

2.2 单一煤层开采覆岩破坏区分布

由图7可知，随着采煤工作面开挖，工作面前方的覆岩破坏范围也随着扩大。采空区上方岩层处于拉应力状态，极值拉应力为0.56 MPa，上覆岩层破坏范围达到32 m。

3 近距离煤层开采扰动覆岩应力场及破坏区分析

基于上述理论力学模型，分析在不同中间岩层厚度( h=12，20，36 m )条件下覆岩应力场、破坏区，以 及其与中间岩层厚度和力学参数的关系。除中间岩层厚度不同外，近距离煤层开采扰动分析的相关参数见表2。

图7 单一煤层开采覆岩应力场云图 Fig. 7 Overburden stress field nephogram in single coal seam mining

3.1 覆岩应力场随中间岩层厚度的演化规律

中间岩层h=20 m时，重复开采覆岩内水平和垂直应力场如图8( c )，( d )所示。σx整体呈现“X”型分布，左右两侧呈椭圆状分布，上下顶板附近呈橄榄状分布。与单一煤层不同的是，左右两侧存在的4个集中区分布于远离中间岩层一侧角点位置。σy呈中心对称分布，顶底板附近呈现钟型分布，左右两侧呈耳状分布。工作面两侧存在应力集中区，两煤层集中区相互连接。

图8( a )，( c )，( e )为不同中间岩层厚度下σx云图。随着中间岩层厚度的增大，“X”型分布区域有所变化，上煤层区域的水平应力σx的拉应力范围缩小；当h＞20 m时，X型区域出现断连，椭圆状分布区域和橄榄状分布区域略有增大。中间岩层左右两侧的压应力逐渐减小，逐渐趋于中间岩层内部应力值。图8( b )，( d ) ，( f )为不同中间岩层厚度下σy云图。随着中间岩层厚度的增加，顶底板附近压应力区没有明显变化。

图8 不同中间岩层厚度下应力演化规律 Fig. 8 Evolution law of stress field with different intermediate rock thickness

3.2 覆岩破坏区随中间岩层厚度的演化规律

由图8可知不同中间岩层厚度( 12，20，36 m )下破坏区演化规律，当中间岩层为12～20 m时，上下2工作面破坏区向外发育并交汇在一起，形成大面积贯通性破坏区，此时中间岩层垮落，并引起上煤层工作面原有破坏区的增大。当中间岩层为36 m时，破坏区有一定发育，但并未贯穿中间岩层，说明随着中间岩层厚度增大，2工作面相互作用导致破坏效果受到限制，中间岩层的支撑作用增强。

中间岩层h=20 m时，随着采煤工作面开挖，工作面前方的覆岩破坏区也随着扩大。采空区上方岩层处于拉应力状态，极值拉应力为0.82 MPa，上覆岩层破坏范围达到58 m。与上煤层开采引起的32 m裂缝带破坏高度相比，破坏区高度明显增大。说明下工作面采动影响范围到达上工作面的煤层采空区，对已破坏覆岩产生二次影响，上工作面原破坏区发生二次发育。

4 近距离煤层重复开采覆岩裂隙演化物理模拟

为了验证理论解析方法对分析近距离煤层上覆岩层破坏区域范围问题的可靠性，对灵露煤矿 Ⅱ2-1和Ⅱ3近距离煤层开采进行相似材料模拟试验，并对比分析重复采动覆岩应力场解析方法的破坏范围与相似材料模拟试验结果，从而证明理论解析法的有效性。

4.1 模型建立

以综合钻孔柱状图作为参考依据，模拟分析 Ⅱ2-1，Ⅱ3两层煤开采条件下上覆岩层的变形破坏情况( 包括覆岩垮落带和导水裂缝带高度 )。其中，上煤层Ⅱ2-1的模拟开采厚度为4 m，下煤层Ⅱ3模拟开采厚度为4 m，煤层间距20 m。铺设模型高度为139 m，用等效荷载代替74.5 m的岩层( 第四纪和泥岩 )厚度。Ⅱ2-1和Ⅱ3近距离煤层的地质条件和煤层厚度与中间岩层h=20 m的理论分析模型相近。

根据灵露煤矿Ⅱ2-1、Ⅱ3煤层采煤工作面地质条件与开采技术条件，相似模型的几何相似比设置为CL=180∶1，容重比为Cγ=1.6∶1，则其应力比为Cσ=288∶1；煤层工作面开采时间与相似模型开挖时间的时间相似常数为：。灵露煤矿真实的煤岩层岩性和厚度及力学性质见表3。

表3 灵露煤矿煤岩层的力学性质 Table 3 Mechanical properties of coal and rock strata

4.2 覆岩垮落带和裂缝带特征

上煤层开采破坏影响区域：当上煤层开挖空间达到42 m时，上覆岩层基本顶岩层发生大范围垮落，垮落高度达到7 m，同时，工作面上方的岩层断裂角很明显，认定42 m为工作面的初次来压。当开挖空间达到63 m时，工作面上覆顶板再次垮落，判定为第2次顶板的周期来压步距( 13 m )。如图9( a )所示，垮落高度达到13 m，裂隙富集区主要集中在图中的椭圆形内，裂隙发育高度32 m。当上煤层开挖空间达到99 m时，工作面上方的岩层垮落角为50°，如图9( b )所示，垮落高度达到12 m，裂隙发育高度35 m。

下煤层开采破坏影响区域：当下煤层工作面推进到45 m时，如图10( a )所示，上覆岩层离层顶板垮落，工作面上方的岩层垮落角为55°，垮落高度发育至16 m，中间岩层未完全破坏。当下煤层工作面推进至终采线位置，如图10( b )所示，上覆岩层垮落高度发育至20 m，即整个中间岩层发生破坏。下煤层开采引起的裂隙富集区主要集中在图中的椭圆形内，与上煤层开采引起的35 m裂缝带高度相比，导水裂缝带发育高度明显增大，达到61 m。Ⅱ3煤层顶板平均周期来压步距为12 m。说明下煤层采动影响范围达到Ⅱ2-1煤层采空区，对已破坏覆岩产生二次影响，上覆岩层裂缝带发生二次发育。重复采动覆岩应力场解析方法的破坏区求解结果与相似材料模拟试验结果相差不超过10%，满足工程应用要求。

图9 上部煤层开采覆岩破坏情况 Fig. 9 Fractured zone of overlying strata of upper coal seam

图10 下部煤层开采覆岩破坏情况 Fig. 10 Fractured zone of overlying strata of lower coal seam

5 结 论

( 1 ) 针对我国西部浅埋近距离煤层重复开采问题，提出一种新的“无限平面双矩形孔口”理论解析模型，并采用Schwarz交替求解法推导了近距离煤层开采扰动覆岩应力场理论计算公式。该理论模型的建立和求解为求解分析近距离煤?层开采扰动的覆岩结构演化规律提供了理论基础。

( 2 ) 近距离煤层重复开采覆岩应力场分布特征：① 近距离煤层重复开采后，煤岩层内的水平应力σx呈现“X”型分布，顶底板均处于拉应力状态，采空区两侧岩体应力呈椭圆状分布；② 随着厚度的增大，中间岩层内应力集中范围和程度降低，下工作面开采引起的覆岩破坏区域二次发育情况也不同。

( 3 ) 近距离煤层重复开采覆岩破坏范围分布规律：当中间岩层厚度达到一定程度，下工作面开采并不会引起原覆岩破坏高度的明显增大。当中间岩层厚度h=20 m时，下工作面开采后上覆岩层破坏区贯穿整个中间岩层，采空区上方岩层处于拉应力状态，下工作面开采引起的原裂缝带发育高度明显增大，下煤层工作面采动引起了上工作面原破坏区域 的二次发育。