孙丽军 叶义成 胡南燕 姚 囝 李鹏程 邓兴敏

（1.中钢集团马鞍山矿山研究总院股份有限公司，安徽 马鞍山 243000；2.武汉科技大学资源与环境工程学院，湖北 武汉 430081；3.湖北省工业安全工程技术研究中心，湖北 武汉 430081）

我国中西部地区矿产资源丰富，但区域内高山、河谷地形交错分布，导致地下采空区顶板出现宏观矿压显现不对称、倾向断裂位置不明确等诸多问题［1］。因此掌握高山河谷地形下顶板变形规律对控制此类采场的矿压及围岩稳定性具有重要的现实意义［2-4］。

高山河谷地形下采空区顶板受非均布荷载作用，目前众多学者探究非均布荷载作用下层状顶板变形规律，多是采用“梁”理论进行分析，如来兴平等［5］基于梁式理论建立了急倾斜工作面顶板的非均布荷载力学模型，分析得到顶板最大沉降点位于中部偏上区域；刑望等［6］基于梁理论分析了采高、工作面长度等因素对顶板最大挠度点位置的变化关系；陈晓祥等［7］建立了非均布荷载下顶板简支梁力学模型，揭示了包括顶板倾斜角度在内的5种因素对顶板下沉的影响规律；冯强等［8］基于弹性地基梁理论获得了岩梁挠度解析表达式；赵明华、张文静［9-10］等通过简支梁理论分析了不同非均布荷载条件下顶板的变形特征。但将顶板简化为二维梁结构后不能完全反映顶板在走向与倾向上的变形特征［11］，因此也有学者指出将顶板视为弹性薄板进行分析，如黄滚、尹光志等［12，13］运用 Navier双三角级数求解了非均布荷载条件下四边简支薄板的挠度解析式；屠洪盛等［1］基于薄板理论揭示了顶板在上覆非均布荷载和充填矸石共同作用下顶板挠曲变形特征。但这些研究仅通过板理论求解了顶板的挠度方程，未系统地研究非均布荷载作用下不同因素对顶板变形特征的影响规律，而且针对高山河谷地形下采空区顶板变形的研究也鲜有报道。

以宜昌某磷矿为研究对象，通过分析高山河谷地形下采空区顶板的受力特征，构建非均布荷载顶板力学模型，基于薄板小挠度弯曲理论，采用Levy单三角级数法推导出高山河谷地形下顶板挠度方程，并据此研究不同因素对顶板最大挠度及其位置的影响规律，以期为高山河谷地形下采空区顶板支护设计提供理论依据。

1 采空区顶板力学模型构建

1.1 工程概况

宜昌某磷矿南矿块属典型的高山河谷地形下矿床，所处地形属山地地貌［14］，地表最大高差约260 m，且地表有常年性河流经过矿体正上方，见图1。整个地质体结构为单斜构造。矿体埋深在140～400 m之间，矿体走向长2.1 km，平均倾角10°，平均厚度3.09 m，属于典型浅埋缓倾斜薄矿床。矿层及围岩以碳酸盐岩为主，直接顶板为含磷白云岩，平均厚度3 m，底板为含钾页岩，平均厚度6 m，研究区可采矿体的赋存情况如图1所示。

1.2 力学模型构建

开采高山河谷地形下缓倾斜薄矿体时，一般在河流正下方设置河流保安矿柱，因此高山河谷地形下的采空区主要位于倾斜的谷坡下。设顶板所受荷载由上覆岩层重力提供，忽略矿体倾角与水平应力的影响，通过简化谷坡地形可将上覆荷载视为梯形非均布荷载。依据薄板的基本假设，由图1可知，直接顶的厚度为3 m，设空区尺寸一般在15 m以上，顶板厚度h与顶板的短边长度L满足。采空区顶板四周均留设有原岩矿柱，边界各点挠度为0，设采空区顶板的边界条件为四边简支，构建如图2所示力学模型。为方便计算，一般将梯形荷载分解为一个三角形荷载分量q（y）和均布荷载分量q0［15］。

2 顶板变形挠度方程求解

四边简支薄板需满足如下边界条件：

基于薄板小挠度弯曲理论的3个基本假设，依据弹性理论空间问题平衡、微分方程式、物理方程式及几何方程式，薄板的弹性曲面微分方程满足：

2.1 三角形荷载顶板挠度解

Navier双三角级数与Levy单三角级数是求解薄板弯曲微分方程的2种经典解法。从2种解法的解析表达式可以看出，Navier解法一般适用于均布荷载条件下的四边简支薄板，即荷载关于薄板几何中心对称，而Levy适用于条件更一般的四边简支薄板［16］。

由高山河谷地形下顶板的受力特征可知，三角形荷载关于x=a/2轴对称，顶板沿y轴呈线性函数递增，更适用于Levy单三角级数解法。因此将薄板的挠度函数展开成单三角级数形式，得到三角形荷载下顶板挠度表达式：

其中，Ym是关于y的任意函数；m为正整数。显然级数式（3）能满足式（1）中x=0及x=a两边的边界条件。

因此，只需选择函数Ym，使式（3）能同时满足式（2）和式（1）中y=0及y=b两边的边界条件。现将式（3）代入式（2），得

将式（4）右边的q(x,y)/D展为的级数，由傅里叶级数展开公式得

与式（4）左侧对比，可见有

解得非齐次方程（6）的通解为

式中，fm(y)是非齐次微分方程（6）的一个特解；前4项是相应齐次方程的通解，系数Am，Bm，Cm，Dm是待定常数，决定于y=0及y=b两边的边界条件。将式（7）代入式（3），即得挠度w的表达式：

现只考虑沿y轴方向的非均布荷载下顶板挠度变形，即q(x,y)是关于y的函数，与x无关，则

于是微分方程（6）的特解可以取为

根据式（1）中与y相关的4个边界条件构建4个方程，可以确定出4个未知数。

分类讨论，当m=1，3，5，…时

因此，当m=1，3，5，…时，求出的系数为有效解，将其代入式（8），得出三角形荷载下顶板挠度w1：

2.2 均布荷载顶板挠度解

文献［16］基于Levy单三角级数求解了坐标原点位于（0，0.5b）处的均布荷载顶板挠度解析式，转换成与2.1节相同的坐标系，即将y坐标向下平移0.5b，得

2.3 非均布荷载顶板挠度解

根据顶板挠度可叠加原理，得到高山河谷地形下采空区非均布荷载顶板挠度方程：

对式（13）顶板挠度方程求偏导，即可得到顶板的最大挠度及最大挠度点坐标，且在计算过程中，一般地取m=1即可满足工程实际的精度要求［12］。

3 顶板变形影响因素分析

从式（11）和式（12）可以看出，影响高山河谷地形下顶板变形的主要因素有顶板自身特性（其中包括弹性模量E、厚度h、泊松比υ等）、顶板所受非均布荷载特征和采空区尺寸。现有研究顶板自身特性对顶板变形的影响规律方面的成果［17-19］已比较成熟，不再赘述。非均布荷载特征对顶板变形的影响主要体现在三角形荷载的倾角θ和均布荷载分量q0值，结合高山河谷地形的具体情况分析，其谷坡倾角范围一般为15°～75°，浅埋缓倾斜矿体的均布荷载分量取值范围为0～4 MPa。此外，分析式（12）与式（13）可知，走向长度a与倾向长度b的比值对计算结果影响较大，令倾向长度为60 m，通过改变走向长度实现a/b值的变化，其变化范围为0.25～4。参考宜昌某磷矿顶板力学参数进行计算，见表1。

3.1 谷坡倾角θ

为了探究高山河谷地形谷坡倾角对顶板变形的影响规律，控制顶板走向与倾向长度均为60 m，均布荷载分量q0=2 MPa，分别设计θ=15°、30°、45°、60°、75°5个水平；其中非均布荷载分量q(y)与谷坡倾角的关系为

将式（14）代入式（11）即可得出顶板挠度函数。图3为倾角θ=75°时的顶板变形特征。

由图3可以看出，求解出的顶板挠曲面符合顶板的边界条件及连续条件，且顶板最大挠度为478.92 mm，位于（30，31.82），即在走向上位于跨中，但倾向上逐渐向高山侧偏移。

图4为不同谷坡倾角下顶板变形特征，当谷坡倾角由15°增大到45°时，顶板最大挠度由216.62 mm缓慢增加到271.86 mm。当倾角继续增大后，顶板最大挠度增速越来越明显，且当倾角为75°时，顶板最大挠度达到478.92 mm，分析认为出现这种变化趋势是由于采空区顶板上覆累积荷载量的变化导致的，此处定义为倾角荷载累积效应。

图 5 中S1、S2、S3、S4、S5为倾角每增加 15°后的荷载增加量，由图可知，当倾角在15°～45°之间时，荷载增加量S1、S2、S3缓慢增加，当倾角大于 45°后，荷载增加量S3、S4、S5变化越来越明显。因此会出现随着倾角的增大，顶板最大挠度先缓慢增加再急剧增加的特征。随着倾角从15°增加到75°，最大挠度点位置逐渐向高山侧偏移，偏移量范围在0.505b～0.53b之间，可认为谷坡倾角对顶板最大挠度点位置影响较小。

3.2 顶板两向长度比a/b

为了探究采空区顶板两向长度比对顶板变形的影响规律，分别设计a/b=0.25、0.5、1、2、4 5个水平，控制河谷倾角为60°，均布荷载分量q0=2 MPa，顶板倾向长度b=60 m，则采空区走向长度分别为15 m、30 m、60 m、120 m、240 m，代入式（13）得出顶板挠度函数。

图6为顶板两向长度比a/b=0.25时的顶板变形量，由图可知，顶板最大挠度点位于（7.5，40.98），已严重偏离倾向跨中位置，但由于走向上采空区跨度较小，使得顶板最大挠度仅4.34 mm。

从图7中可以看出，顶板最大挠度值随a/b值的增加呈“S”形趋势递增，当0.25＜a/b＜0.5时，顶板最大挠度受采空区走向尺寸的影响，此时采空区走向尺寸在15～30 m之间，顶板最大挠度值变化缓慢；当0.5＜a/b＜1时，顶板最大挠度值增加速度变快；当a/b＞1时，采空区倾向尺寸为薄板的短轴，顶板最大挠度值主要受倾向尺寸的影响，由于倾向尺寸保持为60 m不变，因此顶板最大挠度值开始缓慢减小，最终增加速度趋于平稳。

从顶板最大挠度点位置变化趋势来看，顶板最大挠度点位置随着a/b值的增加逐渐向跨中偏移，且偏移速度先快后慢，且变化范围较大，其位置在0.508b～0.683b之间，说明采空区顶板最大挠度点位置受采空区顶板尺寸的影响较大。

3.3 均布荷载分量q0

分别设计q0=0、1 MPa、2 MPa、3 MPa、4 MPa等5个水平，并控制顶板走向与倾向长度均为60 m，河谷倾角为60°，将参数代入式（13）得出顶板挠度函数，图3（a）为q0=0时的顶板挠曲面。

特别地，当q0=0时，即顶板上覆荷载为三角形荷载，顶板最大挠度值为131.72 mm，最大挠度点位置为（30，33.02），即沿倾向y=0.55b处，这与前人研究结果［12］相一致。

由图9可以看出，保持非均布荷载分量q（y）不变时，顶板最大挠度值随q0呈线性递增关系，此外，随着q0值的增加，受非均荷载分量影响逐渐减小，顶板最大挠度点逐渐向跨中偏移，其变化范围在0.513b～0.55b之间，说明非均布荷载分量q0对最大挠度点位置的影响较小。

3.4 实例验证

宜昌某磷矿矿段南矿块地下采空区沿走向长86 m、倾向长68 m，每沿走向隔11 m、倾向隔8 m设置一个4 m×4 m的点柱；地表河谷倾角为68°，均布荷载量q0=3.78 MPa，根据式（13）得出顶板最大挠度位置为y=37.87 m。

现场调研结果显示，位于采空区正中心的3#矿柱出现非对称剥落破坏，靠近高山侧剥落程度较河谷侧更为严重，如图10所示，由此反映出高山河谷地形下采空区上覆顶板变形不对称，使得矿柱受力不均匀；此外，在y=36.72 m处产生了垂直于走向的裂缝，其长度4.3 m、宽5 cm，分析认为是顶板最大弯矩力超过其极限抗拉强度引起的拉伸破坏［15］。对于四边简支薄板来讲，顶板最大挠度与顶板最大弯矩两者所在位置相同，则推断顶板实际最大挠度点位于沿倾向y=36.72 m处，这与理论计算结果基本相符。因此研究得到的顶板挠度计算结果可为高山河谷地形下采空区顶板支护设计提供理论依据。

4 结 论

（1）依据高山河谷地形下采空区顶板受力特征，构建了非均布荷载顶板力学模型，基于Levy单三角级数推导得出了非均布荷载条件下顶板的挠度方程。

（2）分析了河谷倾角θ、采空区两向长度比a/b、均布荷载分量q0对顶板最大挠度及其位置的影响规律。高山河谷地形下采空区顶板最大挠度w与谷坡倾角θ呈先缓慢递增再急剧递增的关系，与采空区两向长度比a/b呈“S”形递增关系，与均布荷载分量q0呈线性递增关系。顶板最大挠度点位置均向高山侧偏移，其偏移距离受采空区两向长度比a/b影响较大，受地表谷坡倾角θ及均布荷载分量q0影响较小。

（3）以调查分析为基础，根据宜昌某磷矿地下采空区顶板及矿柱的破坏特征，分析了顶板最大挠度所在位置为y=36.72 m处，与理论计算结果37.87 m基本相符，验证了研究结果的可靠性。