0.   引　　言

加长工作面长度能够有效提高资源利用率，降低资源损失^[1]，助推煤炭行业实现低碳升级^[2-3]。同时，加长工作面长度也是中厚煤层高产高效、实现年产千万吨的必然趋势^[4-5]。工作面长度的增加会导致来压步距^[6]、来压持续时间以及来压强度等的变化，除此之外还会产生一些异于普通工作面的新特征，具体表现为工作面中部低峰值，而中上部和中下部2个高峰值的三峰值“M”型来压，及在推进方向上工作面表现出了异步来压^[7]的现象。

王家臣等^[8-9]基于蒙特卡洛法建立了含有裂隙的三维模型，得到了裂隙在岩体中呈泊松分布的规律，并指出超长工作面更容易受到裂隙的影响。借助微震监测得到了千米深井超长工作面基本顶具有中部断裂尺寸低于两侧的破断特征，并提出了针对分区破断的支架区域化管理模式。王兆会等^[10]通过岩石力学试验，得到了超前破坏裂隙能够有效降低基本顶破断步距和破断岩块的尺寸，同时将采空区划分为了充分压实区、非充分压实区和悬露顶板支撑区3个区域。陈冬冬等^[11]构建了超长工作面薄板模型，并利用有限差分法计算了煤壁宽度、支撑系数对基本顶破断形态的、破断位置的影响，根据破断计算结果划分了非对称型的“O-X”型破断和非对称型的“C-X”型破断。肖剑儒等^[12]基于压力拱理论计算了工作面的合理长度。王新丰等^[13]基于薄板模型探究了工作面面长突变对顶板破断的特征的影响，通过理论分析与数值模拟得到了面长由小到大过程中，采场顶板“O-X”型破断的时空演化规律，为大面长工作面顶板破断的结构演化提供了支撑。徐亚军等^[14]基于弹性基础梁计算了全长工作面的顶板挠度分布情况，并通过改变计算条件得到了不同巷帮刚度、工作面长度下的挠度分布规律，并得到了与超长工作面实际情况相符的三峰值“W”型^[15-16]分布特征。ZHANG等^[17]基于弹性基础梁计算了不同工作面长度下顶板的峰值分布情况，结合现场电液控热力图，认为随工作面长度的增加，上覆岩层稳定性降低，形成了不稳定的超静定多铰拱结构，在压力分布上表现为由单峰值变为多峰值马鞍形分布。蔺星宇等^[18]基于弹性基础梁求解的挠度分布对支架工作阻力情况进行了分区，并总结了不同区域增阻特性及其对顶板控制的影响。王国法^[19]、徐亚军^[20]等基于弹性支座−连续梁模型和五弯矩方程计算了工作面面长方向的顶板挠度分布特征，将该特征与支架沿顶梁方向载荷分布相结合，提出了液压支架群组支护应力场模型，同时针对不同的工作面长度、巷帮刚度、顶板厚度等进行了分类讨论，并提出了工作面装备群组协同控制技术。

以上学者在超长工作面理论方面进行了广泛而深远的研究，但并没有详细探究支架参数对顶板挠度分布的影响。

在结构力学中通常使用力法、位移法研究超静定结构，其中位移法具有规范、便于编程计算的特点。而矩阵位移法即为位移法加上矩阵方法，矩阵的使用能够降低计算机编程计算的难度，同时矩阵的变换法则极大地提高了位移法在复杂、多样化场景中的应用。矩阵位移法的计算依托于梁的支座和边界条件、梁本身的刚度特性、梁的长度等参数，这便于从控制变量的角度分析各个参数对梁变形的影响，从而反映各个参数在实际生产中可能产生的各种影响。

本文旨在使用矩阵位移法探究不同液压支架相关参数和顶板相关参数对超长工作面顶板挠度的影响，从而为超长工作面液压支架设计提供新的思路。

1.   工作面支护模型简化与矩阵位移法应用

1.1   工作面液压支架与围岩支护模型构建

如图1所示，工作面中液压支架与巷帮共同对顶板起支撑作用，液压支架的主要支撑构件是立柱^[21]，而立柱具有增阻升压、恒阻承载的工作特性，因此可以将其看作具有一定刚度的弹性构件^[22-23]，工作面两侧的巷帮则可以看作是刚度较大的弹性体。假设顶板没有产生断裂，且内部是均匀无裂隙的，即顶板是连续均质的^[24]。则可以将顶板看作是长为工作面长度、高为直接顶厚度、宽为支架控顶距的矩形长条。在本模型中，着重考虑顶板在高度方向上的变形，宽度的影响体现在截面惯性矩，即抗弯刚度中。因为顶板的跨度远远大于直接顶的厚度，且变形也符合小变形假设，因此可以将顶板看作梁模型^[25]。支架与巷帮充当弹性支座，与梁共同组成了支持在弹性支座上的连续梁模型。这样梁承担的是来自上方随动岩层破断产生的均匀分布在梁上的载荷q，而下方则是两端为高强度的弹性支座和中部依次间隔分布的普通弹性支座，如图2所示。

图  1  液压支架支护示意

Figure  1.  Hydraulic support diagram

下载: 全尺寸图片 幻灯片

图  2  支持在弹性支座上的连续梁模型示意

Figure  2.  Schematic diagram of continuous beam model on elastic supports

下载: 全尺寸图片 幻灯片

1.2   液压支架的矩阵位移法

1.2.1   对连杆进行编码

编码是为了更加规范的将单元杆整合到整体结构中。

定义长度为B的具有2个端点的杆件为单元杆，单元杆中的坐标系称为局部坐标系，单元杆左端记为原点，将杆轴向定义为$\overline x $正方向，定义$\overline y $正方向为$\overline x $正方向顺时针旋转90°。沿$\overline x $和$\overline y $方向的线位移记为$\overline \mu $和$\overline \nu $，力记为${\overline F _x}$和${\overline F _y}$，与坐标轴同向为正。角位移$\theta $和力矩$M$以$\overline x $转向$\overline y $为正方向。采用先处理法，将整个杆件细分为若干单元杆，每个杆端具有3个方向的自由度，单元杆两端约束为弹性支座，在杆端能够产生竖向位移。不考虑梁模型轴向方向的变形和力。因此采用如下编码方式：左起第1个单元杆左端记为$ \left( {1,2} \right) $，右端记为$ \left( {3,4} \right) $。第2个杆左端与第一个杆右端为同一结点，因此记为$ \left( {3,4} \right) $，右端记为$ \left( {5,6} \right) $。以此类推，第$n$个单元杆左端记为$ \left( {2n + 1,2n + 2} \right) $，右端记为$ \left( {2n + 3,2n + 4} \right) $。使用${{\boldsymbol{\lambda }}_n}$表示第$n$个单元杆的编码向量，（图3），则单元杆的编码向量为

图  3  单元杆编码

Figure  3.  Element bar coding

下载: 全尺寸图片 幻灯片

1.2.2   单元刚度矩阵

在梁端弹性支座约束处施加一刚臂^[26]作为附加约束来限制梁在此处的转动（图4），并使梁两端产生单位转角$\theta {'_A}$和$\theta {'_B}$（图5a）。称添加了附加约束的结构为基本结构，弹性支座不影响基本结构的结果。此时杆端弯矩由位移法可得：

图  4  附加约束

Figure  4.  Additional constraints

下载: 全尺寸图片 幻灯片

图  5  弯矩作用下的单元杆模型

Figure  5.  Element bar model under bending moment

下载: 全尺寸图片 幻灯片

式中：$i$为杆件的线刚度，$i$=$EI/B$。其中$EI$为抗弯刚度，E为顶板的弹性模量，Pa；$I$为顶板的截面惯性矩，m⁴；$B$为支架中心距，即单元杆的长度，m；$ {M_{A'B'}} $、$ {M_{B'A'}} $分别为$ A'B' $杆左端弯矩、右端弯矩。

令两端各产生竖向位移${\nu _1}$和${\nu _2}$，如图5b所示，此时梁端产生相应转角$\theta '{'_A}$和$\theta '{'_B}$为弦转角$\varphi $：

式中：$ {\overline \nu _1} $、$ {\overline \nu _2} $分别为局部坐标系中梁左端和右端的沿$\overline y $方向的位移。

应用叠加原理联立式（2）—式（3），可得梁两端有力矩且产生了相对位移情况下的转角方程：

用线位移和角位移表示弯矩：

上式即为转角位移方程。由平衡方程求杆两端剪力$ {\overline F _{y1}} $、$ {\overline F _{y2}} $：

解得：

将梁两端弯矩表达式与剪力表达式联立为矩阵得：

与胡克定律类似，力等于刚度与形变量的乘积，因此将上式记为单元刚度方程：

式中：$ {\overline {\boldsymbol{F}} ^{\text{e}}} $为局部坐标下的杆端力向量，$ {\overline {\boldsymbol{F}} ^{\text{e}}} = {( {{{\overline F }_{y1}}}\;\;\;{{{\overline M }_1}}\;\;\;{{{\overline F }_{y2}}}\;\;\;{{{\overline M }_2}} )^{\text{T}}} $，矩阵内元素分别为梁左端的剪力和弯矩、梁右端的剪力、弯矩；$ {\overline {\boldsymbol{k}} ^{\text{e}}} $为局部坐标系下的单元刚度矩阵；$ {\overline {\boldsymbol{\varDelta }} ^{\text{e}}} $为局部坐标系下的杆端位移向量，$ {\overline {\boldsymbol{\varDelta }} ^{\text{e}}} = {(\begin{array}{*{20}{c}} {{{\overline v }_1}}&{{{\overline \theta }_1}}&{{{\overline v }_2}}&{{{\overline \theta }_2}} \end{array})^{\text{T}}} $，矩阵内元素分别为梁左端沿$\overline y $方向的位移分量、梁左端截面转角、梁右端的沿$\overline y $方向位移分量、梁右端截面转角。

单元刚度矩阵是对称矩阵，单元刚度矩阵中$\overline {{k_{cj}}} $元素为第$j$个杆端位移分量等于1（其他位移分量为0）时引起的第$c$个杆端力分量；第$c$行元素的意义是当4个杆端位移分量分别等于1时，引起的第$c$个杆端力分量的值；第$j$列元素的意义是当第$j$个杆端位移分量等于1时引起的4个杆端力的值。

研究对象为水平连续梁模型，因此单元转换矩阵为单位阵，转换后的单元刚度矩阵与原局部坐标下的单元刚度矩阵一致：

1.2.3   整体刚度矩阵

每个隔离杆件的约束条件均为两端是弹性支座约束，弹性支座的影响不体现在单元刚度矩阵而是在整体刚度矩阵中^[27]，需在整体刚度矩阵弹性支座位移所对应的主系数处增加弹簧刚度${k_{\text{s}}}$，目的在于仅当弹簧支座处产生了${\varDelta _y} = 1$的位移时，弹簧产生支反力为$k \times 1$。根据力的平衡条件求解梁时，除该处剪力外还应当加上弹性支座的支反力。使用单元集成法^[28]将单元刚度矩阵按照式(1)编码连接为整体刚度矩阵${\boldsymbol{K}}$，矩阵上侧与右侧为编码值，空白处填充0：

整体刚度矩阵为对称矩阵、带状矩阵、稀疏矩阵、非奇异矩阵。与单元刚度方程同理，整体刚度方程为

式中：$ {\boldsymbol{F}} $为节点力；$ {\boldsymbol{\varDelta }} $为节点位移。该刚度方程是根据虚设位移$ {\boldsymbol{\varDelta }} $推算节点力$ {\boldsymbol{F}} $的关系式，并不能直接反应原实际载荷。

1.2.4   等效节点载荷

当载荷单独作用于梁上而不产生位移时，基本结构中的节点约束力为$ {{\boldsymbol{F}}_{\text{P}}} $，当节点位移$ {\boldsymbol{\varDelta }} $单独作用时，基本结构中的节点约束力为$ {\boldsymbol{F}} = {\boldsymbol{K\varDelta }} $。实际载荷与位移同时作用于基本结构，且满足附加约束刚臂的约束力与约束力矩均为0时，认为基本体系转化为原结构，则有：

联立式(13)：

采用能在基本结构中产生相同节点约束力的等效节点载荷替代原载荷，则等效节点载荷$ {\boldsymbol{P}} $为

考虑到梁模型中局部坐标与整体坐标正方向一致，则单元的等效节点载荷向量$ {{\boldsymbol{P}}^{\text{e}}} $为

式中：$ \overline {\boldsymbol{F}} _{\text{P}}^{\text{e}} $为单元固端约束力向量，当长度为B的梁上方施加与梁等长的均布载荷q时，梁两端的单元固端约束力向量$ \overline {\boldsymbol{F}} _{\text{P}}^{\text{e}} $的组成为

由单元集成法按照编码式(1)进行加法运算，节点载荷计算示意如图6所示，可得整体坐标系中等效节点载荷向量P：

图  6  节点等效载荷计算示意

Figure  6.  Schematic of node equivalent load calculation

下载: 全尺寸图片 幻灯片

1.2.5   解位移方程与杆端内力

由式(15)和式(16)可得基本方程：

将整体刚度矩阵${\boldsymbol{K}}$与等效结点载荷${\boldsymbol{P}}$代入式(20)可以解得全梁位移向量${\boldsymbol{\varDelta }}$。

将对应位置的竖向位移$ {\varDelta _l} $与该处弹性支座的刚度$ {k_{{\text{Rl}}}} $相乘即可求得该处弹性支座的支座反力$ {F_{{\text{Rl}}}} $：

因忽略了轴向变形和力，梁模型的各个单元杆端内力向量表达式一致：

式中：$ {{\boldsymbol{F}}^{\text{e}}} $为杆端内力向量；$ {{\boldsymbol{\varDelta }}^{\text{e}}} $为对应编码位置的位移向量。

考虑到支座为弹性支座，需要在计算杆端剪力时，在一侧的表达式中加上支座刚度，即可得到结点杆端剪力。

2.   超长工作面顶板挠度求解——以小保当二号煤矿超长工作面为例

液压支架宽度与直接顶的弹性模量均取实际数据，分别为2.05 m和35 GPa，工作面长度为450 m，液压支架中心距为2.05 m，采高平均2.55 m，岩石容重取$\gamma = 25$ kN/m³；梁的宽度取为支架控顶距长度7.024 m，梁高度为10.32 m，梁模型的截面惯性矩为643.34 m⁴；液压支架的等效刚度取0.15 GN/m，巷帮刚度取20倍的液压支架等效刚度3 GN/m。顶板承受的载荷通过关键层理论计算^[29-30]：

式中：$ {{{q}}_{t,1}} $为顶板收到来自其上随动岩层的载荷，Pa；$ {E_t} $为顶板上部第t层的弹性模量，Pa；$ {h_t} $为第t层的厚度，m；$ {\gamma _t} $为第t层的容重，kN/m³。计算其上随动岩层的载荷约为1 360 182 Pa。

将上述参数应用于整体刚度矩阵，即式(12)，可以求得位移向量$ {\boldsymbol{\varDelta }} $的结果。值得一提的是，该求解是单向的，即已知位移向量求解力是不可行的。该位移向量$ {\boldsymbol{\varDelta }} $是由挠度和转角组成的向量，提取挠度可以绘制全梁挠度曲线，如图7所示。

图  7  支护区域及对应挠度值

Figure  7.  Support area and corresponding deflection values

下载: 全尺寸图片 幻灯片

从图中可以看出，在全梁长方向上，工作面顶板挠度呈现三峰值形态，比较明显的是在40号与188号支护区域附近呈现高峰值，而在2个高峰值之间出现了位于曲线的中轴线上的极大值，该极大值相较于两侧并不明显。从整体来看，受两侧较高强度的巷帮影响，梁的两端挠度值较小，随着向中间的移动，挠度值迅速增大，且增长速度呈现递减趋势，直至遇到第1个峰值，随后挠度值开始缓慢减小至一个极小值，随后又开始非常缓慢地增加，直至增加到位于对称轴线上的低峰值。梁是对称的，因此我们将左端初始值至第1个峰值段称为第1段曲线，而在第1峰值至中轴线段称为第2段曲线。因为第2段曲线中极小值到中轴线上低峰值的增幅很小，因此将这段称为近水平段。

函数化的结果更利于计算的使用，若要使用该曲线，则需要对其进行数据拟合来找到该曲线的函数表达式。从曲线外形上来看，梁左端到第1个峰值段，在本例中为1～40号支护区域，挠度虽然在持续增加，但增加的幅度逐渐减小，总体变化规律与三次多项式相符，即该段可以使用$f(x) = a{x^3} + b{x^2} + cx + d$进行拟合。

经测试，使用三次多项式来进行第1段曲线部分的数据拟合能够获得非常好的拟合结果，这可能与单跨梁微分方程通解中最高次幂多为三次有关。若选择拟合目标为左端至第一峰值段，则拟合的R²为0.999 8，残差量级为10⁻⁴级，显著性水平p值均远小于0.05，这说明三次多项式对该段曲线几乎完美解释，使用曲线拟合数据点如图8a所示。而实际情况中我们需要的可能只是该曲线的一小部分。当只需要拟合曲线的一小部分，比如19～24号之间，三次多项式的拟合结果是R²为1的完美拟合，如图8b所示。利用该种情况，可以在仅知道4个点的挠度及其位置的情况下，得到第1段曲线的表达式，即可以通过监测实际生产中第1段曲线范围内的4台支架在一段时间内的平均下沉量及支架在工作面中的位置，来估算全工作面的下沉量的最大值及液压支架的理论最大工作阻力。

图  8  对“第一段曲线”进行三次多项式拟合

Figure  8.  Cubic polynomial fitting of “first segment curve”

下载: 全尺寸图片 幻灯片

3.   不同参数、工艺对顶板挠度的影响

从计算过程中看，对挠度可能产生影响的因素分别是生产系统的影响，包含设计面长；支架因素，包含支架的等效刚度、支架的数量、支架的宽度；生产环境的影响，包含顶板的弹性模量，顶板的截面惯性矩（顶板的高度与计算宽度），巷帮的刚度以及顶板之上随动岩层产生的载荷。

3.1   生产系统的影响因素

3.1.1   固定面长为450 m，动态改变支架宽度与支架数量的组合

4个组合分别为：架宽1.50 m，共300架；架宽1.75 m，共257架；架宽2.05 m，共220架；架宽2.40 m，共188架，而其他参数不变时的挠度曲线如图9a所示。第一组组合中支架数量最多，支架宽度最小，该组挠度的最值、初始值、整体水平均是最小的。第一组合的左侧高峰值位于序号50附近，约左起75 m处，挠度约为0.013 94 m，而第四组合的左侧高峰值位于序号35处，左起84 m处，挠度约为0.022 44 m。假设原矿井采用第一组合的支架数与支架宽度，工作面长度不变的情况下，改用第四组合的支架数与宽度会使最大峰值扩大到原来的1.6倍，并使最大峰值位置向工作面中部转移约9 m。这意味着更多的支架会处于高应力区域，高应力的极值也会增加。从梁模型角度分析，梁长不变，提供支反力的支座数量变少，每一个支座分摊的压力增加。这要求支架在设计时，若采用较大的支架中心距，则支架的工作阻力也应提高。

图  9  不同参数下的挠度曲线

Figure  9.  Deflection curves under various parameters

下载: 全尺寸图片 幻灯片

3.1.2   支架数量和宽度

本模型计算中，限定工作面长度时，支架数量与支架宽度成反比，且支架数量与支架宽度的乘积等于工作面长度。因此无法通过固定2个参数来观察另一个参数变化对梁整体挠度的影响，但比较不同工作面长度下，相同支架数量不同支架宽度或相同支架宽度不同支架数量依旧可以得到一些规律。

当固定支架宽度为2.05 m，而调整支架数量为200、210、220、230、240架，也可以理解为调整工作面长度为410、430.5、451、471.5、492 m时，梁的挠度曲线如图9b所示。从图中不难看出，工作面长度动态变化下，支架数量的增加对第1段曲线基本没有影响，也可以说支架数量对梁两端的挠度值没有影响，但会显著拉长第2段曲线的长度，且支架数量的增加并没有减缓第2段曲线的变化速度。支架数量的继续增加也只是延长了近水平段的影响范围，而对曲线的极值点和最值点影响不大。也就是说，支架宽度不变，支架数量与工作面长度增加时，梁上方的载荷总量虽然在增加，但对两端的支架几乎没有影响，只是增加了近水平段的影响范围，增加了处于高应力区的支架比例。

当固定支架数量为220架，而调整支架宽度为1.5、1.75、2.05、2.4 m，即固定支架数量而调整工作面长度为330、385、451、528 m时，梁的挠度曲线如图9c所示。支架宽度越大，梁上挠度的最小值和最大值越大，2个高峰值相对于低峰值越明显，近水平段的极大值与极小值差异越明显，第1段曲线的变化率越大，第22段曲线的水平段影响范围越大。从梁模型的角度分析，支架数量不变但是支架宽度增加，梁的跨度增大，梁上的载荷总量增加，分摊到每个支座上的载荷增加。

需要特别指出的时，本节所比较的参数取值较为保守，从而不易看出工作面长度变化对挠度的整体影响。当固定支架宽度为2.05 m而支架数量为100架，即工作面长度为205 m时，工作面挠度分布为单峰值，峰值出现在工作面中部。当固定支架数量为220架，而调整支架宽度为不断减少时，工作面挠度表现出了“多峰值且对称中部为极大值”到“多峰值但对称中部为极小值”再到“整体为单峰值”的转变趋势。即工作面长度的减少能够明显导致顶板挠度出现多峰值到单峰值的转变。

3.2   支架的影响因素

3.2.1   支架等效刚度

仅调整支架等效刚度为0.7×10⁸、1.0×10⁸、1.3 ×10⁸、1.6×10⁸、1.9×10⁸ N/m，而其他参数不变时的挠度曲线如图9d所示。从图中变化中可以看出，液压支架的等效刚度对挠度的影响呈递减趋势，且支架等效刚度越大，挠度的整体水平越小，挠度的变化率越小，高峰值向两侧移动，峰值间差异变小。支架等效刚度的增大更像是从上往下直接压缩了整个挠度曲线。可以理解为，支架是抵抗梁产生变形的存在，而支架等效刚度是梁“抵抗”能力的强弱，支架“抵抗”能力越强，在相同载荷情况下，梁自身的挠度变化越小，需要梁依靠自身变形抵抗的力越小，梁的挠度变化也就越小。与3.1.2中支架宽度变化对梁挠度的影响曲线对比可以发现，两者对挠度曲线的影响是相似的，只是支架宽度对挠度影响更加缓和。从模型上出发，支架宽度的变化更偏向于改变了来自梁下方的集中力分布间隔的大小，而等效刚度的变化更偏向于改变了集中力本身的大小。

3.3   生产环境的影响

3.3.1   顶板弹性模量

仅调整顶板弹性模量为0.5×10¹⁰、1.0×10¹⁰、3.5×10¹⁰、8.0×10¹⁰、13.5×10¹⁰ Pa，而其他参数不变时的挠度曲线如图9e所示。顶板弹性模量在实际应用中往往是固定的，但考虑模型应用于类似条件的矿井时，该处理可以更好的进行计算上的移植。弹性模量的增加会使梁端的挠度初始值增加，在降低高峰值的同时，使高峰值向工作面中部移动，第2段曲线的水平段缩减，高应力影响区的支架数量减少，但第2段曲线的极小值差别不大。弹性模量可以看作是梁本身抵抗变形的能力。弹性模量越大，梁抵抗突变的能力越大，最大值与最小值越向中间靠近，挠度曲线整体越“平缓”；而弹性模量越小，梁柔性越好，越容易出现局部的挠度峰值，且曲线的变化也会更加陡峭。

3.3.2   顶板截面惯性矩

仅调整顶板截面惯性矩为200、471、643、1 000、2 000 m⁴，而其他参数不变时的挠度曲线如图9f所示。截面惯性矩反映的是形状对弯曲的抵抗作用，本次计算中截面的宽度采用的是支架的控顶距，实际可能不同。该部分可以看作是直接顶的高度与宽度立方的积对挠度变化的影响。增加截面惯性矩，梁两端的挠度初始值增加，梁的挠度最大值降低并向梁内移动，第2段曲线的近水平段范围缩减，但最小值差别不大。截面惯性矩对梁的挠度的影响与顶板弹性模量对梁挠度的影响基本一致，这是因为在计算中，截面惯性矩与顶板弹性模量以乘积的形式，即抗弯刚度，共同对挠度产生影响。

3.3.3   巷帮刚度

仅调整巷帮刚度为0.25×10⁹、0.5×10⁹、1.0×10⁹、2.5×10⁹、5.0×10⁹ N/m，而其他参数不变时的挠度曲线如图9g所示。巷帮可以看作是梁模型两侧的约束条件，将其视作不同于支架的等效弹簧，可以查看巷帮等效刚度对梁挠度的影响。当巷帮刚度与支架等效刚度差别不大时，梁整体挠度曲线变化很小，三峰值差距不大。随着巷帮刚度的提高，梁两端挠度初始值明显增大，三峰值的差异逐渐明显，而近水平段基本没有变化，高峰值的位置基本不变。从图中可以明显看到，巷帮刚度对两侧的挠度影响明显，从梁模型角度分析，巷帮刚度越大，巷帮的变形越小，巷帮侧的梁模型表现越趋向于悬臂梁。而巷帮的刚度越小，巷帮侧就越趋向于简支梁，挠度的变化就越平缓。

3.3.4   均布载荷

仅调整直接顶上方载荷为0.8×10⁶、1.0×10⁶、1.36×10⁶、1.8×10⁶、2.5×10⁶ Pa，而其他参数不变时的挠度曲线如图9h所示。载荷的大小对梁挠度的影响与支架等效刚度对梁挠度影响基本一致，原因是2个都是施加于梁的力，只不过方向和分布不同。以梁为研究对象，顶部随动岩层施加于梁向下的均布载荷，而支座给予梁规律间隔的向上的集中力，两者在影响上是相同的，此处不再赘述。

4.   超长工作面支护数值模拟

4.1   超长工作面模型构建与支护单元布置

采用3DEC进行超长工作面的数值模拟，以x轴正方向为工作面倾向方向，y轴正方向为推进方向，z轴为模型的高度。模型总长度为550 m，其中左右两侧各50 m的煤柱充当巷帮，中间为450 m工作面；模型宽度为200 m，模拟开挖100 m，分10次进行；模型高度为100 m，其中底板高5 m，煤层高3 m，上部岩层按照采高30倍设计为90 m，实际使用了92 m，这样模型整体高度为100 m整。模型的4个侧边限制横向移动，即限制节点x和y方向速度为0，模型的底边限制纵向移动，即限制节点z方向速度为0，模型整体施加沿z轴负方向的重力加速度，设置顶面自由并施加向下的垂直载荷模拟上覆岩层重力。模型开挖后的整体情况如图10a所示。

图  10  3DEC中的超长工作面布置

Figure  10.  Layout of ultra-long working face in 3DEC

下载: 全尺寸图片 幻灯片

模型所用岩层参数以实际岩石试验报告为准，在模型中插入了随机DFN（Discrete Fracture Network）离散裂隙网格用以切割模型产生随机裂隙，保证节理具有一定的随机性，图10a中的红色圆片即为切割模型的随机裂隙平面fracture。采用结构单元桩模型pile充当液压支架，该模型在本次应用中只考虑了沿轴向的受力情况，不考虑桩本身的倾斜、弯曲，因此桩模型可以看作弹簧而具有与支架等效的作用。

桩模型的模拟原则是从50 m处开始，实际工作面所用液压支架最大控顶距为7.024 m，顶梁长5.530 m，柱间距0.970 m，支架中心距为2.05 m。基于上述参数，在支架顶梁的4个角处布置四根桩结构单元充当支架来模拟对顶板的支护，桩模型的布置如图11所示。以该种方式在450 m工作面布置桩结构单元，共计安放的桩结构组数为220组，其中一组为单独的2个桩，则共有桩模型数量为878个。设置桩结构单元的刚度与第二章理论计算中所采用支架刚度一致，为0.15 GN/m，桩结构单元的布置如图10b、图10c所示。

图  11  桩单元布置俯视图

Figure  11.  Top view of pile element arrangement

下载: 全尺寸图片 幻灯片

4.2   数值模拟结果分析

推进100 m后，工作面处的中部桩结构单元顶部受力情况如图10 d所示，将桩的顶端位移数据使用fish导出，以桩的编号为横轴，桩顶端的位移为纵轴作图12。因左前柱与右前柱图形基本一致，选取左前柱作为支架前柱的代表，后柱同理。在模型中存在DFN随机节理的情况下，距煤壁1.5 m处表现出了两侧挠度最小，而在中上部和中下部存在挠度极大值，在中部出现了明显的极小值。图12a中右侧的高峰值相较于左侧不是十分明显，考虑到DFN切割的随机性、划分块体的不均衡以及pile桩布置的不对称问题可能是导致右侧高峰值低于左侧高峰值的原因。图中下沉量并不明显，这可能与模型块体的划分较大有关，较大的块体更容易在采空后形成铰接平衡结构^[31]，从而弱化底部支护体的效果。同时值得提出的一点是，在同一个支架的4个pile柱模型中，两根后柱顶端下沉量基本一致，两根前柱顶端下沉量基本一致，在中上部峰值位置，后柱顶端下沉量水平约为前柱的1.67倍。

图  12  桩结构单元顶端挠度分布

Figure  12.  Deflection distribution at top of pile structural element

下载: 全尺寸图片 幻灯片

从结果上来看，在左前柱位置，即距煤壁约1.50 m处，顶板的下沉表现出明显的两侧最小值、中上部和中下部存在极大值，而在中部出现了极小值的三峰值“M”型特征。距离煤壁7.024 m处的后柱处依旧表现出了两侧最小值，中部存在极小值的特征，但相较于前柱，中上部和中下部的极值表现不再明显，顶板挠度整体增大，且表现出向中部积聚的现象。

5.   实际应用分析

以2023年7月26日小保当2号井132202工作面电液控数据为为准，每隔3个支架记录1次，以每个支架20:00所在循环的循环末阻力作为该支架的数据值，每间隔1 d绘制1条全工作面的电液控数据拟合曲线图。使用R语言大数据包^[32]绘制结果如图13所示。从图中能够非常明显的看出，在以7月26日—7月29日期间的工作面工作阻力呈现出了“M”型来压特征，具体表现为工作面两端的工作阻力较小，而工作面中部产生了2个极大值，在125号支架附近出现了一个极小值。

图  13  连续4 d 20:00整电液控数据统计

Figure  13.  Electro-hydraulic control data statistics

下载: 全尺寸图片 幻灯片

同时值得一提的是，在28日20:00循环末阻力曲线中，160号支架附近出现了一个相对较弱的峰值，经与源数据比对发现，在155号、159号、161号支架的工作阻力均大于40 MPa，甚至159号支架出现了远大于左侧峰值的50.9 MPa的情况，但由于该段处于高峰值的支架数量较少且间隔取值容易漏掉极值情况，使得曲线拟合时弱化了该极值而呈现了一个较弱的峰值。从整体结果上看，现场数据与上述分析基本一致。

6.   结　　论

1)基于矩阵位移法的结构力学计算方法能够得到超长工作面的顶板挠度两端小，中间大的“M”型特征，同时可以计算梁的杆端内力。

2)三次多项式可以精确的拟合第1段曲线，利用拟合函数能够预测最大峰值前支架的工作阻力情况，同时为细化支架顶梁载荷计算提供理论支持。

3)在矩阵位移法中的“常数量”均会对顶板挠度产生不同的影响，其中支架宽度、支架等效刚度和载荷强度的影响体现在挠度的整体大小上，而支架数量、顶板的弹性模量、顶板的截面惯性矩、巷帮刚度的影响多体现在对挠度曲线形状的影响，比如两最大值之间间距、高应力区的影响范围、两端最小值的大小、第1段曲线的变化幅度、第2段曲线的极值情况。

4)当生产系统要求的工作面设计长度固定时，提高支架的宽度、降低液压支架的布置数量会导致顶板整体挠度水平显著提高，梁端的初始挠度提高，同时使两峰值间距减小。

在本次研究中，理论计算的结果与数值模拟结果、实际生产中的电液控数据三者均具有相同的特征，以上3个手段共同作用可以作为超长工作面三峰值“M”型来压特征的有力证明。