支护力作用下双向不等压采场围岩应力和位移解析解

Analytical Solution of Surrounding Rock Stress and Displacement in Bidirectional Unequal Pressure Stope Under Supporting Stress

摘要

地下倾斜矩形采场的稳定性与顶板应力及位移的分布密切相关，而目前关于矩形开挖围岩应力及位移的解析方法应用于大宽高比倾斜矩形采场时误差较大，且并未考虑采场倾角和支护力的影响.本文基于复变函数理论，提出了大宽高比倾斜矩形采场的映射函数表达式，推导了采场在双向不等压应力和支护力作用下围岩应力和位移解析解，分析了采场倾角、宽高比和支护力对顶板应力和位移分布的影响.结果表明：解析解与FLAC模拟解的偏差小于5%；随着采场倾角的增大，顶板的应力和位移非对称分布趋势加剧，卸压程度和位移量逐渐降低；随着采场宽高比的增大，顶板卸压程度和位移量逐渐增大；采场支护体施加的支护力能够改善采场围岩的应力环境，降低顶板的下沉量.

Abstract

The stability of inclined rectangular stope is closely related to the stress and displacement distribution of stope roof in underground mines. However， the current theoretical method for stress and displacement of surrounding rock in rectangular excavation has a significant error when applied to inclined rectangular stopes with large width‑height ratio， and the influence of excavation dip angle and supporting stress is not considered. Based on the complex function theory， this paper puts forward the mapping function expression of inclined rectangular stope with large width‑height ratio， and derives the analytical solution of stress and displacement of surrounding rock in stope under bidirectional unequal pressure and supporting stress. Besides， the influence of stope dip angle， width‑height ratio and supporting stress on stress and displacement distribution of stope roof is analyzed. The results show that the deviation between analytical solution and FLAC simulation solution is less than 5%. In addition， the asymmetric distribution trend of stress and displacement around stope roof intensifies as the stope dip angle increases， and the degree of pressure relief and vertical displacement gradually decreases. With the increase of stope width‑height ratio， the pressure relief degree and the displacement of stope roof gradually increase. The supporting stress exerted by supporting bodies can improve the stress environment of surrounding rock and reduce the roof subsidence.

Graphical abstract

关键词

倾斜矩形采场 / 宽高比 / 复变函数理论 / 应力和位移 / 支护力

Key words

inclined rectangular stope / width‑height ratio / complex variable theory / stress and displacement / supporting stress

在地下采矿工程领域，为了提高矿石资源的开采效率，倾斜矩形采场结构参数逐渐向大宽高比方向发展，即沿采场倾向方向的宽度与开采高度的比值逐渐增加.然而，大宽高比采场结构顶板稳定性问题尤为突出，其更易诱发顶板冒落、垮塌等地压灾害^［1］.此外，受矿体倾角影响，倾斜矩形采场围岩应力和位移呈现偏对称分布特征，其破坏机制与水平开挖空间有较大差异^［2］.究其原因，采场的几何结构和应力环境对围岩应力位移具有较大影响，因此准确求解不同因素下采场围岩应力及位移分布特征是保证采场稳定的前提.

对于倾斜采场围岩应力和位移分布规律的研究，主要有三种研究方法.其一，数值模拟方法主要用于现场实际开挖工程问题^［3］，其能够模拟复杂地质条件下围岩的应力和位移分布特征. Ravandi等^［4］运用FLAC^3D建模计算了非静水应力场中圆形、D形和马蹄形隧道的壁面位移；Abbo等^［5］使用有限元极限分析方法研究了平面应变条件下矩形隧道的稳定性.该方法针对具体问题建立单独的计算模型，并对模型尺寸和网格密度进行敏感性分析，从而确保模拟结果的准确性，但是过程稍繁琐，且无法直观揭示各变量对结果的影响.其二，通过在现场布置监测设备能够获得开采引起的围岩应力和位移变化^［6］，但是受监测周期长、设备成本高等因素影响，限制了现场工业试验的开展.其三，理论分析方法作为理解问题本质最有效的研究方法，其能够直观地揭示不同变量对围岩应力和位移分布的影响规律^［7］.

鉴于理论分析方法的优势，许多学者采用弹塑性理论求解了地下圆形巷道开挖围岩应力和位移的弹性和弹塑性解^［8-9］.然而，对于非圆形水平孔口开挖问题，目前尚处于研究阶段.Muskhelishvili提出的复变函数方法在分析无限大弹性平面矩形孔口问题具有突出的优势^［10-11］，被广泛用于推导非圆形开挖问题中应力和位移的解析解.在矩形开挖空间围岩应力和位移求解方面，Huo等^［12］推导了宽高比为1∶1~3∶1条件下深埋水平矩形开挖空间围岩的应力和位移解；Zhao等^［13］研究了不同侧压力系数条件下方形隧道围岩的应力分布.此外，Lei等^［14］和Ng等^［15］提出了在远场倾斜单向应力作用下的水平矩形孔围岩应力位移的解析解.综上，对于地下矩形开挖空间而言，大多数研究集中在对宽高比为3∶1以下水平巷道或硐室等开挖问题的分析，其采用3项映射函数即可得到较为精确的解答，但是对于大宽高比条件下矩形采场围岩应力和位移求解时存在误差大、精度低等缺点.此外，倾角作为倾斜矩形采场的一个重要结构特征，其对采场的应力和位移分布有很大影响，目前多数研究并未考虑开挖空间自身角度的影响.

除了地下采场几何特征的影响外，为了提高大结构参数采场顶板稳定性，采场内布置的矿柱和锚杆等支护体能够对采场围岩提供一定的支护力^［16］，其对采场围岩的应力及位移分布也有较大影响.在支护体作用下巷道或采场围岩应力及位移解析解研究方面，Wang等^［17］利用复势理论计算了圆形衬砌隧道周围的应力和位移.基于弹性力学叠加原理，Wang等^［18］和Gao等^［19］提出了圆形和非圆形衬砌隧道应力和位移的解析解；Dong等^［20］推导了具有有限厚度衬砌的圆形隧道特性的解析解.然而，目前的研究主要集中于简单几何形状（圆形或椭圆形）的水平巷道或隧道，而支护力作用下的倾斜矩形采场围岩应力和位移解求解问题尚未得到解决.

鉴于目前对非圆形开挖空间围岩应力和位移解析解的研究现状，国内外学者主要研究了开挖几何结构、边界应力或支护力等单一因素的影响.然而，实际的采场开采问题是地下工程中各种因素相互作用的结果.本文旨在推导大宽高比结构的倾斜矩形采场围岩应力和位移分布的解析解，同时考虑远场双向不等压边界应力和均匀支护力的联合效应，从而用于初步预测倾斜矩形采场围岩的稳定性并对支护设计提供理论基础.

1 力学模型构建

1.1 基本假设

根据倾斜矿体赋存特征，矿体沿走向方向一般延伸几百米至上千米，远大于单个采场沿矿体倾向方向的长度，可将单个采场开挖围岩应力和位移求解问题近似视为平面应变问题.然而在矿体实际开采过程中，采场围岩的岩体条件和应力环境较为复杂，在计算之前需作出以下假设：

1） 围岩为均质、各向同性的连续介质，只产生弹性变形，且断面尺寸远小于模型整体尺寸；

2） 由于采场的开采高度远小于其埋深，因此可忽略采场附近由于重力引起的梯度效应，模型边界应力视为均匀分布；

3） 对于采场内原岩矿柱或人工矿柱等支护体，将其支护作用等效为与采场顶底板垂直的均匀支护力F.

1.2 计算原理

根据简化后倾斜矩形采场的几何特征和应力环境，建立如图1a所示力学模型.以矩形采场形心为原点建立直角坐标系xoy，其中x轴和y轴分别与采场长边和短边保持平行，采场倾角为δ.此外，采场受到竖直方向上的远场边界应力σ_v=p，水平方向上的远场边界应力σ_h=kσ_v，k为侧压力系数.

支护力作用下双向不等压倾斜矩形采场围岩应力及位移解的计算流程如图2所示.首先需要分析大宽高比矩形采场的保角变换，本文通过增加映射函数项数并修正对应系数，将物理平面上的大宽高比矩形采场变换为像平面上单位圆.随后考虑采场倾角和支护力两个因素，基于弹性力学复变函数理论求解远场应力和支护力联合作用下的应力势函数，从而计算围岩应力和位移解析分量.在求解围岩应力势函数过程中，根据弹性力学的叠加原理，支护力作用下应力势函数可以分解成两部分单独求解后叠加.首先计算远场应力作用、无支护力情况下的采场围岩应力势函数，如图1b所示；然后计算无远场原岩应力、仅有支护力情况下的采场围岩应力势函数，如图1c所示；最后将两种情况下的解析函数叠加，即可求得均匀支护力作用下采场围岩的应力复变势函数.

2 支护力作用下倾斜采场围岩应力及位移解析解推导

2.1 复变函数保角变换

在复平面z中（z=x+iy），假设深埋矩形采场宽、高分别为2a，2b，则宽高比m=a/b，根据确定的映射函数能够将物理平面z上采场外部区域变换至像平面ζ上单位圆内部区域，如图3所示.

对于地下巷道或硐室等宽高比较小的开挖结构，多数研究提出映射函数取至ζ³项即可满足精度要求.然而，地下倾斜采场沿矿体倾向方向的宽度与开采高度的比值可达9∶1，本文通过增加映射函数表达式项数至ζ⁹项，实现对大宽高比采场结构围岩应力及位移的计算分析.由于z平面倾斜矩形采场关于x和y轴对称，映射函数偶数次项系数均为0，确定映射函数基本形式为^［21］

z=\omega (\zeta )=R(\frac{\mathrm{1}}{\zeta }+c_{\mathrm{1}}\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}}+\dots ).（1）

式中：R为采场形状参数；c_{\mathrm{1}}=\frac{m-\mathrm{1}}{m+\mathrm{1}}，c_{\mathrm{3}}=\frac{-\mathrm{2}m}{\mathrm{3}{(m+\mathrm{1})}^{\mathrm{2}}}，c_{\mathrm{5}}=\frac{-\mathrm{2}m(m-\mathrm{1})}{\mathrm{5}{(m+\mathrm{1})}^{\mathrm{3}}}，c_{\mathrm{7}}=\frac{-\mathrm{2}m(m^{\mathrm{2}}-\mathrm{3}m+\mathrm{1})}{\mathrm{7}{(m+\mathrm{1})}^{\mathrm{4}}}， c_{\mathrm{9}}=\frac{-\mathrm{2}m(m^{\mathrm{3}}-\mathrm{6}{m}^{\mathrm{2}}+\mathrm{6}m-\mathrm{1})}{\mathrm{9}{(m+\mathrm{1})}^{\mathrm{5}}}.

对于ζ平面中的单位圆，取ρ为1，那么ζ=\rho {\mathrm{e}}^{\mathrm{i}\theta }=\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}\theta，则式（1）可以写为

\left.\begin{array}{l}x=R(\mathrm{c}\mathrm{o}\mathrm{s}\theta +c_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\theta +c_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{3}\theta +c_{\mathrm{5}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{5}\theta +c_{\mathrm{7}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{7}\theta +c_{\mathrm{9}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{9}\theta ),\\ y=R(-\mathrm{s}\mathrm{i}\mathrm{n}\theta +c_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\theta +c_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{3}\theta +c_{\mathrm{5}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{5}\theta +c_{\mathrm{7}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{7}\theta +c_{\mathrm{9}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{9}\theta ).\end{array}\right\}（2）

根据保角变换的映射关系，z平面上倾斜矩形采场右边界中点（a，0），顶板中点（0，b）和右顶角（a，b）与像平面ζ上的点（1，0），（1，-π/2）以及（1，θ）分别相互对应.结合式（2）可得

\left.\begin{array}{l}a=R(\mathrm{1}+c_{\mathrm{1}}+c_{\mathrm{3}}+c_{\mathrm{5}}+c_{\mathrm{7}}+c_{\mathrm{9}})，\\ b=R(\mathrm{1}-c_{\mathrm{1}}+c_{\mathrm{3}}-c_{\mathrm{5}}+c_{\mathrm{7}}-c_{\mathrm{9}}).\end{array}\right\}（3）

因此，根据映射函数确定的z平面中倾斜矩形采场的宽高比为

\frac{a}{b}=m\frac{\mathrm{315}{m}^{\mathrm{4}}+\mathrm{1012}{m}^{\mathrm{3}}+\mathrm{1812}{m}^{\mathrm{2}}+\mathrm{888}m+\mathrm{263}}{\mathrm{263}{m}^{\mathrm{4}}+\mathrm{888}{m}^{\mathrm{3}}+\mathrm{1812}{m}^{\mathrm{2}}+\mathrm{1012}m+\mathrm{315}}>m.（4）

根据式（4）看出，通过映射函数确定的宽高比较实际矩形孔口形状的大，会产生一定的映射偏差.为了提高计算的精度，引入修正系数ε（λ=εm）对保角映射函数进行修正，将λ替换m代入式（4），可得

\epsilon m\frac{\mathrm{315}{(\epsilon m)}^{\mathrm{4}}+\mathrm{1012}{(\epsilon m)}^{\mathrm{3}}+\mathrm{1812}{(\epsilon m)}^{\mathrm{2}}+\mathrm{888}\epsilon m+\mathrm{263}}{\mathrm{263}{(\epsilon m)}^{\mathrm{4}}+\mathrm{888}{(\epsilon m)}^{\mathrm{3}}+\mathrm{1812}{(\epsilon m)}^{\mathrm{2}}+\mathrm{1012}\epsilon m+\mathrm{315}}=m.（5）

由式（5）可计算出矩形采场在不同宽高比m条件下的修正系数ε，将λ替换m代入式（1），即可求得修正后的映射函数系数.此处，设计倾斜矩形采场开采高度为4 m，求得不同宽高比条件下映射函数系数的值，如表1所示.

由表1可知，随着矩形采场宽高比m的增大，修正系数ε逐渐减小，表明随着矩形采场宽高比的逐渐增大，未修正前映射函数计算的孔口形状与实际断面形状偏差逐渐增大.为了验证修正后映射函数的准确性，拟合了宽高比为8∶1条件下3项映射函数、未修正的6项映射函数和修正后的6项映射函数的保角映射采场形状，如图4所示.可以看出，宽高比为8∶1条件下，3项映射函数与矩形采场实际断面形状相差较大，6项映射函数降低了保角映射偏差，但其与实际形状仍然存在差异.基于本文修正后的6项映射函数获得的采场形状与预期的形状偏差可以忽略不计，其能够满足复变函数应力及位移的计算.

2.2 应力势函数推导

2.2.1 远场边界应力作用下应力势函数

在求解矩形采场围岩应力和位移解析解时，应先寻找满足边界条件的势函数.首先讨论仅在远场应力作用下倾斜矩形采场围岩应力势函数的推导.根据式（1）倾斜矩形采场保角映射函数，令ζ=σ，可得（符号上横线表示函数的共轭）

\left.\begin{array}{l}\omega (\sigma )=R(\frac{\mathrm{1}}{\sigma }+c_{\mathrm{1}}\sigma +c_{\mathrm{3}}{\sigma }^{\mathrm{3}}+c_{\mathrm{5}}{\sigma }^{\mathrm{5}}+c_{\mathrm{7}}{\sigma }^{\mathrm{7}}+c_{\mathrm{9}}{\sigma }^{\mathrm{9}}),\\ \overline{\omega (\sigma )}=R(\sigma +\frac{c_{\mathrm{1}}}{\sigma }+\frac{c_{\mathrm{3}}}{{\sigma }^{\mathrm{3}}}+\frac{c_{\mathrm{5}}}{{\sigma }^{\mathrm{5}}}+\frac{c_{\mathrm{7}}}{{\sigma }^{\mathrm{7}}}+\frac{c_{\mathrm{9}}}{{\sigma }^{\mathrm{9}}}),\\ \overline{\omega \text{'}(\sigma )}=R(-{\sigma }^{\mathrm{2}}+c_{\mathrm{1}}+\frac{\mathrm{3}{c}_{\mathrm{3}}}{{\sigma }^{\mathrm{2}}}+\frac{\mathrm{5}{c}_{\mathrm{5}}}{{\sigma }^{\mathrm{4}}}+\frac{\mathrm{7}{c}_{\mathrm{7}}}{{\sigma }^{\mathrm{6}}}+\frac{\mathrm{9}{c}_{\mathrm{9}}}{{\sigma }^{\mathrm{8}}}).\end{array}\right\}（6）

根据远场应力边界条件，有

\left.\begin{array}{l}B=\frac{\mathrm{1}}{\mathrm{4}}({\sigma }_{\mathrm{h}}+{\sigma }_{\mathrm{v}})=\frac{(k+\mathrm{1})p}{\mathrm{4}},\\ B\text{'}+\mathrm{i}C\text{'}=-\frac{\mathrm{1}}{\mathrm{2}}({\sigma }_{\mathrm{h}}-{\sigma }_{\mathrm{v}}){\mathrm{e}}^{-\mathrm{2}\mathrm{i}\delta }=-\frac{(k-\mathrm{1})p}{\mathrm{2}}{\mathrm{e}}^{-\mathrm{2}\mathrm{i}\delta }.\end{array}\right\}（7）

其中：B′，C′为应力边界函数参数.

根据复变函数理论，应力边界条件可以表述为

\begin{array}{c}{f}_{\mathrm{0}}=\mathrm{i}\int (\overline{f_x}+\overline{f_y})\mathrm{d}s-\frac{\overline{F_x}+\mathrm{i}\overline{F_y}}{\mathrm{2}\mathrm{\pi }}\mathrm{l}\mathrm{n}\sigma -\\ \frac{\mathrm{1}+\mu }{\mathrm{8}\mathrm{\pi }}(\overline{F_x}-\mathrm{i}\overline{F_y})\frac{\omega (\sigma )}{\overline{\omega \text{'}(\sigma )}}\sigma -\mathrm{2}B\omega (\sigma )-(B\text{'}-\mathrm{i}C\text{'})\overline{\omega (\sigma )}.\end{array}（8）

其中，μ为泊松比.

倾斜矩形采场内无支护力作用时，采场内边界所受面力F=0，所以\overline{f_x}=\overline{f_y}=\overline{F_x}=\overline{F_y}=\mathrm{0}（式中\overline{f_x}和\overline{f_y}分别为矩形孔口内边界在x和y方向上的面力分量）.联合式（6）~（8）并应用柯西积分公式可得

\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}\int \frac{f_{\mathrm{0}}}{\sigma -\zeta }\mathrm{d}\sigma =\frac{k-\mathrm{1}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}\delta }pR\zeta -\\ \frac{k+\mathrm{1}}{\mathrm{2}}pR(c_{\mathrm{1}}\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}}).\end{array}（9）

根据弹性力学复变函数理论，在远场边界应力作用下复变势函数可表达如下：

{\Phi }_{\mathrm{1}}(\zeta )=\frac{{\phi }_{\mathrm{1}}\text{'}(\zeta )}{\omega \text{'}(\zeta )}；{\Psi }_{\mathrm{1}}(\zeta )=\frac{{\psi }_{\mathrm{1}}\text{'}(\zeta )}{\omega \text{'}(\zeta )}.（10）

其中：

\left.\begin{array}{l}{\phi }_{\mathrm{1}}(\zeta )=\frac{\mathrm{1}+\mu }{\mathrm{8}\mathrm{\pi }}(\overline{F_x}+\mathrm{i}\overline{F_y})\mathrm{l}\mathrm{n}\zeta +B\omega (\zeta )+{\phi }_{\mathrm{0}}(\zeta )=\frac{k+\mathrm{1}}{\mathrm{4}}p\omega (\zeta )+{\phi }_{\mathrm{0}}(\zeta )，\\ {\psi }_{\mathrm{1}}(\zeta )=-\frac{\mathrm{3}-\mu }{\mathrm{8}\mathrm{\pi }}(\overline{F_x}-\mathrm{i}\overline{F_y})\mathrm{l}\mathrm{n}\zeta +(B\mathrm{\text{'}}+\mathrm{i}C\mathrm{\text{'}})\omega (\zeta )+{\psi }_{\mathrm{0}}(\zeta )=-\frac{k-\mathrm{1}}{\mathrm{2}}p{\mathrm{e}}^{-\mathrm{2}\mathrm{i}\delta }\omega (\zeta )+{\psi }_{\mathrm{0}}(\zeta ).\end{array}\right\}（11）

为了求解式（11）中应力势函数φ₁（ζ）和ψ₁（ζ），需要先确定函数φ₀（ζ）和ψ₀（ζ）：

{\phi }_{\mathrm{0}}(\zeta )+\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}{\int }_{\sigma }\frac{\omega (\sigma )}{\overline{\omega \text{'}(\sigma )}}\frac{\overline{{\phi }_{\mathrm{0}}\text{'}(\sigma )}}{\sigma -\zeta }\mathrm{d}\sigma =\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}\int \frac{f_{\mathrm{0}}}{\sigma -\zeta }\mathrm{d}\sigma ,（12a）

{\psi }_{\mathrm{0}}(\zeta )+\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}{\int }_{\sigma }\frac{\overline{\omega (\sigma )}}{\omega \text{'}(\sigma )}\frac{{\phi }_{\mathrm{0}}\text{'}(\sigma )}{\sigma -\zeta }\mathrm{d}\sigma =\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}\int \frac{\overline{f_{\mathrm{0}}}}{\sigma -\zeta }\mathrm{d}\sigma .（12b）

函数φ₀（ζ）在单位圆内是ζ的连续解析函数，其解析函数可表示为

{\phi }_{\mathrm{0}}(\zeta )={\sum }_{j=\mathrm{1}}^{\infty }{\alpha }_j{\zeta }^j,j=\mathrm{1,2},\mathrm{3},\dots .（13）

将式（6），（9）和（13）代入式（12a），结合柯西积分公式和等式两边系数相等定理可得

\left.\begin{array}{l}{\alpha }_{\mathrm{1}}=pR(t_{\mathrm{10}}+\mathrm{i}{t}_{\mathrm{11}}),{\alpha }_{\mathrm{3}}=pR\left(\frac{\mathrm{2}{t}_{\mathrm{5}}{t}_{\mathrm{10}}-(k+\mathrm{1})t_{\mathrm{7}}}{\mathrm{2}{t}_{\mathrm{8}}}+\mathrm{i}\frac{t_{\mathrm{6}}{t}_{\mathrm{11}}}{t_{\mathrm{9}}}\right),\\ {\alpha }_{\mathrm{5}}=pR\left(\frac{\mathrm{2}{t}_{\mathrm{10}}(t_{\mathrm{3}}{t}_{\mathrm{8}}+\mathrm{3}{c}_{\mathrm{9}}{t}_{\mathrm{5}})-(k+\mathrm{1})(t_{\mathrm{7}}+c_{\mathrm{5}}{t}_{\mathrm{8}})}{\mathrm{2}{t}_{\mathrm{8}}}-\frac{\mathrm{i}{t}_{\mathrm{11}}(t_{\mathrm{3}}{t}_{\mathrm{9}}+\mathrm{3}{c}_{\mathrm{9}}{t}_{\mathrm{6}})}{t_{\mathrm{9}}}\right),\\ {\alpha }_{\mathrm{7}}=pR(c_{\mathrm{9}}{t}_{\mathrm{10}}-\frac{k+\mathrm{1}}{\mathrm{2}}{c}_{\mathrm{7}}-\mathrm{i}{c}_{\mathrm{9}}{t}_{\mathrm{11}})，{\alpha }_{\mathrm{9}}=-\frac{k+\mathrm{1}}{\mathrm{2}}pR{c}_{\mathrm{9}}.\end{array}\right\}（14）

其中：

\left.\begin{array}{l}{t}_{\mathrm{4}}=c_{\mathrm{1}}+\mathrm{7}{c}_{\mathrm{7}}{c}_{\mathrm{9}},t_{\mathrm{5}}=\mathrm{5}{t}_{\mathrm{3}}{c}_{\mathrm{9}}+t_{\mathrm{2}},t_{\mathrm{6}}=\mathrm{5}{t}_{\mathrm{3}}{c}_{\mathrm{9}}-t_{\mathrm{2}},\\ t_{\mathrm{7}}=c_{\mathrm{3}}+\mathrm{5}{c}_{\mathrm{5}}{c}_{\mathrm{9}},t_{\mathrm{8}}=\mathrm{1}-\mathrm{15}{c_{\mathrm{9}}}^{\mathrm{2}}-\mathrm{3}{t}_{\mathrm{3}},t_{\mathrm{9}}=\mathrm{1}-\mathrm{15}{c_{\mathrm{9}}}^{\mathrm{2}}+\mathrm{3}{t}_{\mathrm{3}},\\ t_{\mathrm{10}}=\frac{(k-\mathrm{1})\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\delta -(k+\mathrm{1})t_{\mathrm{4}}-\mathrm{3}{t}_{\mathrm{5}}{t}_{\mathrm{7}}(k+\mathrm{1})/t_{\mathrm{8}}-\mathrm{5}{c}_{\mathrm{5}}{t}_{\mathrm{3}}(k+\mathrm{1})}{\mathrm{2}(\mathrm{1}-\mathrm{7}{c_{\mathrm{9}}}^{\mathrm{2}}-\mathrm{3}{t_{\mathrm{5}}}^{\mathrm{2}}/t_{\mathrm{8}}-\mathrm{5}{t_{\mathrm{3}}}^{\mathrm{2}}-t_{\mathrm{1}})},\\ t_{\mathrm{11}}=\frac{(k-\mathrm{1})\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\delta }{\mathrm{2}(\mathrm{1}-\mathrm{7}{c_{\mathrm{9}}}^{\mathrm{2}}-\mathrm{3}{t_{\mathrm{6}}}^{\mathrm{2}}/t_{\mathrm{9}}-\mathrm{5}{t_{\mathrm{3}}}^{\mathrm{2}}+t_{\mathrm{1}})}.\end{array}\right\}（15）

因此，远场边界应力作用下函数φ₀（ζ）可以表达为

{\phi }_{\mathrm{0}}(\zeta )={\alpha }_{\mathrm{1}}\zeta +{\alpha }_{\mathrm{3}}{\zeta }^{\mathrm{3}}+{\alpha }_{\mathrm{5}}{\zeta }^{\mathrm{5}}+{\alpha }_{\mathrm{7}}{\zeta }^{\mathrm{7}}+{\alpha }_{\mathrm{9}}{\zeta }^{\mathrm{9}}.（16）

与上述求解流程类似，函数ψ₀（ζ）可以通过求解式（12b）确定.由应力边界函数表达式可得

\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}{\int }_{\sigma }\frac{\overline{f_{\mathrm{0}}}}{\sigma -\zeta }\mathrm{d}\sigma =\frac{k-\mathrm{1}}{\mathrm{2}}pR{\mathrm{e}}^{-\mathrm{2}\mathrm{i}\delta }(c_{\mathrm{1}}\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}})-\frac{k+\mathrm{1}}{\mathrm{2}}pR\zeta .（17）

将式（6），（16）和（17）代入式（12b），函数ψ₀（ζ）可以表达为如下形式：

\begin{array}{c}{\psi }_{\mathrm{0}}(\zeta )=\frac{k-\mathrm{1}}{\mathrm{2}}pR{\mathrm{e}}^{-\mathrm{2}\mathrm{i}\delta }(c_{\mathrm{1}}\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}})-\frac{k+\mathrm{1}}{\mathrm{2}}pR\zeta -[{\alpha }_{\mathrm{1}}{f}_{\mathrm{2}}(\zeta )-(\mathrm{3}{t}_{\mathrm{1}}{\alpha }_{\mathrm{3}}+\mathrm{5}{t}_{\mathrm{2}}{\alpha }_{\mathrm{5}}+\mathrm{7}{t}_{\mathrm{3}}{\alpha }_{\mathrm{7}}+\mathrm{9}{c}_{\mathrm{9}}{\alpha }_{\mathrm{9}})\zeta +\\ \mathrm{3}{\alpha }_{\mathrm{3}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{2}}-(\mathrm{5}{t}_{\mathrm{1}}{\alpha }_{\mathrm{5}}+\mathrm{7}{t}_{\mathrm{2}}{\alpha }_{\mathrm{7}}+\mathrm{9}{t}_{\mathrm{3}}{\alpha }_{\mathrm{9}}){\zeta }^{\mathrm{3}}+\mathrm{5}{\alpha }_{\mathrm{5}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{4}}-(\mathrm{7}{t}_{\mathrm{1}}{\alpha }_{\mathrm{7}}+\mathrm{9}{t}_{\mathrm{2}}{\alpha }_{\mathrm{9}}){\zeta }^{\mathrm{5}}+\mathrm{7}{\alpha }_{\mathrm{7}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{6}}-\mathrm{9}{t}_{\mathrm{1}}{\alpha }_{\mathrm{9}}{\zeta }^{\mathrm{7}}+\mathrm{9}{\alpha }_{\mathrm{9}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{8}}].\end{array}（18）

因此，将式（16）和（18）代入式（11），即可得到应力势函数φ₁（ζ）和ψ₁（ζ），根据式（10）即可求得远场边界应力作用下复变应力势函数.

2.2.2 仅均匀支护力作用下应力势函数

在没有远场边界应力作用下，根据式（7）可知B=B′=C′=0.此时，采场顶底板围岩只受到均匀支护力F的作用，即矩形孔口内边界受到沿y轴方向上大小相等的均布面力F，如图1c所示.由于矩形采场内边界上所有面力在x轴和y轴方向上的主矢量为0，则\overline{F_x}=\overline{F_y}=\mathrm{0}.此外，

\left.\begin{array}{l}\overline{f_x}=F\mathrm{c}\mathrm{o}\mathrm{s}(N,x)=\mathrm{0},\\ \overline{f_y}=F\mathrm{c}\mathrm{o}\mathrm{s}(N,y)=-F\frac{\mathrm{d}x}{\mathrm{d}s}.\end{array}\right\}（19）

根据式（8）可得应力边界函数为

\begin{array}{c}{f}_{\mathrm{0}}=\mathrm{i}\int (\overline{f_x}+\mathrm{i}\overline{f_y})\mathrm{d}s=\int F\mathrm{d}x=\frac{\mathrm{1}}{\mathrm{2}}F(z+\overline{z})=\\ \frac{\mathrm{1}}{\mathrm{2}}FR\left((c_{\mathrm{1}}+\mathrm{1})\sigma +\frac{c_{\mathrm{1}}+\mathrm{1}}{\sigma }+c_{\mathrm{3}}{\sigma }^{\mathrm{3}}+\frac{c_{\mathrm{3}}}{{\sigma }^{\mathrm{3}}}+c_{\mathrm{5}}{\sigma }^{\mathrm{5}}+\right.\\ \left.\frac{c_{\mathrm{5}}}{{\sigma }^{\mathrm{5}}}+c_{\mathrm{7}}{\sigma }^{\mathrm{7}}+\frac{c_{\mathrm{7}}}{{\sigma }^{\mathrm{7}}}+c_{\mathrm{9}}{\sigma }^{\mathrm{9}}+\frac{c_{\mathrm{9}}}{{\sigma }^{\mathrm{9}}}\right).\end{array}（20）

由柯西积分公式，可得

\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}{\int }_{\sigma }\frac{f_{\mathrm{0}}}{\sigma -\zeta }\mathrm{d}\sigma =\frac{\mathrm{1}}{\mathrm{2}}FR\left[\left(c_{\mathrm{1}}+\mathrm{1}\right)\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}})\right].（21）

与2.2.1节复函势函数求解思路类似，需要根据式（12a）和（12b）确定函数φ₀（ζ）和ψ₀（ζ），将式（6），（13）和（21）代入式（12a），此处将符号α替换为β以作区分，可得

\left.\begin{array}{l}{\beta }_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}}FR{t}_{\mathrm{12}},{\beta }_{\mathrm{3}}=\frac{\mathrm{1}}{\mathrm{2}}FR\frac{t_{\mathrm{7}}+t_{\mathrm{5}}{t}_{\mathrm{12}}}{t_{\mathrm{8}}},\\ {\beta }_{\mathrm{5}}=\frac{\mathrm{1}}{\mathrm{2}}FR\left[c_{\mathrm{5}}+\frac{\mathrm{3}{c}_{\mathrm{9}}{t}_{\mathrm{7}}+t_{\mathrm{12}}(\mathrm{3}{c}_{\mathrm{9}}{t}_{\mathrm{5}}+t_{\mathrm{3}}{t}_{\mathrm{8}})}{t_{\mathrm{8}}}\right],\\ {\beta }_{\mathrm{7}}=\frac{\mathrm{1}}{\mathrm{2}}FR(c_{\mathrm{7}}+c_{\mathrm{9}}{t}_{\mathrm{12}}),{\beta }_{\mathrm{9}}=\frac{\mathrm{1}}{\mathrm{2}}FR{c}_{\mathrm{9}}.\end{array}\right\}（22）

其中t_{\mathrm{12}}=\frac{t_{\mathrm{4}}+\mathrm{3}{t}_{\mathrm{5}}{t}_{\mathrm{7}}/t_{\mathrm{8}}+\mathrm{5}{c}_{\mathrm{5}}{t}_{\mathrm{3}}+\mathrm{1}}{\mathrm{1}-\mathrm{7}{c_{\mathrm{9}}}^{\mathrm{2}}-\mathrm{3}{t_{\mathrm{5}}}^{\mathrm{2}}/t_{\mathrm{8}}-\mathrm{5}{t_{\mathrm{3}}}^{\mathrm{2}}-t_{\mathrm{1}}}.

因此，仅考虑均匀支护力作用下函数φ₀（ζ）表达式为

{\phi }_{\mathrm{0}}(\zeta )={\beta }_{\mathrm{1}}\zeta +{\beta }_{\mathrm{3}}{\zeta }^{\mathrm{3}}+{\beta }_{\mathrm{5}}{\zeta }^{\mathrm{5}}+{\beta }_{\mathrm{7}}{\zeta }^{\mathrm{7}}+{\beta }_{\mathrm{9}}{\zeta }^{\mathrm{9}}.（23）

此外，根据柯西积分定理可得

\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\mathrm{i}}{\int }_{\sigma }\frac{\overline{f_{\mathrm{0}}}}{\sigma -\zeta }\mathrm{d}\sigma =\frac{\mathrm{1}}{\mathrm{2}}FR[(c_{\mathrm{1}}+\mathrm{1})\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+\\ c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}}].\end{array}（24）

将式（6），（23）和（24）代入式（12b），并将符号α替换为β，可得

{\psi }_{\mathrm{0}}(\zeta )=FR[(c_{\mathrm{1}}+\mathrm{1})\zeta +c_{\mathrm{3}}{\zeta }^{\mathrm{3}}+c_{\mathrm{5}}{\zeta }^{\mathrm{5}}+c_{\mathrm{7}}{\zeta }^{\mathrm{7}}+c_{\mathrm{9}}{\zeta }^{\mathrm{9}}]/\mathrm{2}-

[{\beta }_{\mathrm{1}}{f}_{\mathrm{2}}(\zeta )-(\mathrm{3}{t}_{\mathrm{1}}{\beta }_{\mathrm{3}}+\mathrm{5}{t}_{\mathrm{2}}{\beta }_{\mathrm{5}}+\mathrm{7}{t}_{\mathrm{3}}{\beta }_{\mathrm{7}}+\mathrm{9}{c}_{\mathrm{9}}{\beta }_{\mathrm{9}})\zeta +

\mathrm{3}{\beta }_{\mathrm{3}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{2}}-(\mathrm{5}{t}_{\mathrm{1}}{\beta }_{\mathrm{5}}+\mathrm{7}{t}_{\mathrm{2}}{\beta }_{\mathrm{7}}+\mathrm{9}{t}_{\mathrm{3}}{\beta }_{\mathrm{9}}){\zeta }^{\mathrm{3}}+

\mathrm{5}{\beta }_{\mathrm{5}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{4}}-(\mathrm{7}{t}_{\mathrm{1}}{\beta }_{\mathrm{7}}+\mathrm{9}{t}_{\mathrm{2}}{\beta }_{\mathrm{9}}){\zeta }^{\mathrm{5}}+

\mathrm{7}{\beta }_{\mathrm{7}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{6}}-\mathrm{9}{t}_{\mathrm{1}}{\beta }_{\mathrm{9}}{\zeta }^{\mathrm{7}}+\mathrm{9}{\beta }_{\mathrm{9}}{f}_{\mathrm{2}}(\zeta ){\zeta }^{\mathrm{8}}].（25）

由于B=B′=C′=0，根据式（11）可求得仅在均匀支护力作用下应力双调和函数为

{\phi }_{\mathrm{2}}(\zeta )={\phi }_{\mathrm{0}}(\zeta );{\psi }_{\mathrm{2}}(\zeta )={\psi }_{\mathrm{0}}(\zeta ).（26）

因此，在不考虑远场边界应力、仅有采场顶底板边界支护力作用下，倾斜矩形开挖空间围岩复变势函数为

{\Phi }_{\mathrm{2}}(\zeta )=\frac{{\phi }_{\mathrm{2}}\mathrm{\text{'}}(\zeta )}{\omega \mathrm{\text{'}}(\zeta )};{\Psi }_{\mathrm{2}}(\zeta )=\frac{{\psi }_{\mathrm{2}}\mathrm{\text{'}}(\zeta )}{\omega \mathrm{\text{'}}(\zeta )}.（27）

2.3 采场围岩应力和位移解析分量

基于弹性力学应力叠加原理，可以求得倾斜矩形采场在远场边界应力和顶底板均匀支护力联合作用下的应力势函数，即

\left.\begin{array}{l}\Phi (\zeta )={\Phi }_{\mathrm{1}}(\zeta )+{\Phi }_{\mathrm{2}}(\zeta ).\\ \Psi (\zeta )={\Psi }_{\mathrm{1}}(\zeta )+{\Psi }_{\mathrm{2}}(\zeta ).\end{array}\right\}（28）

利用推导出的复变应力势函数，可以计算支护力作用下任意倾角矩形采场围岩的应力分布.用复变量表示的采场围岩应力分量可以推导如下：

\left.\begin{array}{l}{\sigma }_y+{\sigma }_x=\mathrm{4}\mathrm{R}\mathrm{e}[\Phi (\zeta )],\\ {\sigma }_y-{\sigma }_x+\mathrm{2}\mathrm{i}{\tau }_{xy}=\frac{\mathrm{2}}{\omega \text{'}(\zeta )}[\overline{\omega (\zeta )}\Phi \text{'}(\zeta )+\omega \text{'}(\zeta )\Psi (\zeta )].\end{array}\right\}（29）

其中：σ_x 和σ_y 是xoy坐标系下的应力分量；Re表示函数的实部.此外，xoy坐标系下矩形采场围岩的位移分量可以表达为

u_x+\mathrm{i}{u}_y=\frac{\mathrm{1}}{\mathrm{2}G}[\kappa {\phi }_{\mathrm{0}}(\zeta )-\frac{\omega (\zeta )}{\overline{\omega \text{'}(\zeta )}}\overline{{\phi }_{\mathrm{0}}\text{'}(\zeta )}-\overline{{\psi }_{\mathrm{0}}(\zeta )}].（30）

其中：κ=3-4 μ；G=E/2（1+μ），E和μ分别为围岩岩体的弹性模量和泊松比.由式（30）可以计算得到仅远场应力作用下倾斜矩形采场围岩的位移分量u_x1和u_y1，以及仅均匀支护力作用下位移分量u_x2和u_y2.基于弹性力学叠加原理，可以得到倾斜矩形采场支护力作用下的围岩位移为

\left.\begin{array}{l}{u}_x=u_{x\mathrm{1}}+u_{x\mathrm{2}},\\ u_y=u_{y\mathrm{1}}+u_{y\mathrm{2}}.\end{array}\right\}（31）

需要指出的是，应力分量σ_x，σ_y 和τ_xy 以及位移分量u_x 和u_y 是在图5所示的xoy坐标系中定义的，其与正常坐标系存在一定的夹角.根据弹性力学坐标转换方法，将应力和位移分量转换成图5所示x'oy'坐标系的相应分量，具体为

\left.\begin{array}{l}{\sigma }_{x\mathrm{\text{'}}}={\sigma }_x\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\delta +{\sigma }_y\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\delta +\mathrm{2}{\tau }_{xy}\mathrm{s}\mathrm{i}\mathrm{n}\delta \mathrm{c}\mathrm{o}\mathrm{s}\delta ,\\ {\sigma }_{y\mathrm{\text{'}}}={\sigma }_x\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\delta +{\sigma }_y\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\delta -\mathrm{2}{\tau }_{xy}\mathrm{s}\mathrm{i}\mathrm{n}\delta \mathrm{c}\mathrm{o}\mathrm{s}\delta ,\\ {\tau }_{x\mathrm{\text{'}}y\mathrm{\text{'}}}={\tau }_{xy}(\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\delta -\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\delta )+({\sigma }_y-{\sigma }_x)\mathrm{s}\mathrm{i}\mathrm{n}\delta \mathrm{c}\mathrm{o}\mathrm{s}\delta .\end{array}\right\}（32）

同时，直角坐标系x'oy'中围岩的水平位移和垂直位移可以计算为

\left.\begin{array}{l}{u}_{x\text{'}}=u_x\mathrm{c}\mathrm{o}\mathrm{s}\delta +u_y\mathrm{s}\mathrm{i}\mathrm{n}\delta ,\\ u_{y\text{'}}=u_y\mathrm{c}\mathrm{o}\mathrm{s}\delta -u_x\mathrm{s}\mathrm{i}\mathrm{n}\delta .\end{array}\right\}（33）

3 采场围岩应力及位移解析解验证

选取辽宁省丹东市某金矿浅部20 m中段采场为研究对象，采场埋深为320 m，采场沿矿体倾向宽度和开采高度分别为24 m和4 m，即宽高比为6∶1，倾角为20°.根据我国地应力场分布特征^［22］，缓倾斜采场受到的竖直方向的远场边界应力为8.7 MPa，水平构造应力为7.4 MPa，侧压力系数为0.85.由于采场内矿柱仅对顶底板施加支护力，在验证实例中将支柱的支护力等效为均匀支护力垂直施加于采场顶底板，支护力F设为1.0 MPa.

为了验证倾斜矩形开挖围岩应力和位移解析解的准确度，建立FLAC数值计算模型，如图5所示.数值模型尺寸为x×y=200 m×200 m，该模型内均采用四边形单元划分.模型采用弹性本构模型，岩体弹性模量为17.5 GPa，泊松比为0.23.为了提高开挖采场围岩应力和位移的计算精度，在采场附近布置精细网格，距离采场较远的区域布置相对粗糙的网格.在数值模拟中，模型上下左右边界均施加应力边界条件，采场顶底板位置施加应力边界模拟支护力，模拟过程中先对整个模型施加初始应力状态，然后进行开挖.

对于金属矿山地下采场开采，采场顶板的稳定性对于矿体的安全开采具有重要的意义.将数值模拟与不同映射函数确定的顶板位置应力和位移解析解进行对比，绘制如图6所示曲线图，其中横轴对应于图5中A点至B点的位置坐标.规定压应力为负，拉应力为正，位移沿x′和y′轴的正方向为正.可以看出，基于修正的6项映射函数计算的围岩应力和位移分布特征和量值与FLAC结果偏差最小.

为了进一步验证解析解的正确性，给出解析解与数值模拟解的平均相对偏差系数s：

s=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n\left|X_{\mathrm{t}\mathrm{h}}-X_{\mathrm{s}\mathrm{i}}\right|/\left|X_{\mathrm{s}\mathrm{i}}\right|.（34）

式中：s为平均相对偏差系数；n为顶板测点个数；X_th和X_si分别代表各测点的解析解和模拟解.

以数值模拟结果为对比基准，求得顶板区域应力和位移的平均相对偏差系数，结果如表2所示. 3项映射函数计算的顶板围岩应力及位移均存在较大的偏差，其中水平应力相对偏差系数达到130%以上，垂直位移约为17%.随着映射函数项数的增加，6项映射函数计算结果的平均相对偏差有了大幅度降低，说明映射函数项数的增加能够提高解析解的精度.采用修正后的6项映射函数得到的解析解偏差最小，应力和位移偏差能够控制在5%之内.

4 不同几何形态采场应力和位移分析

4.1 不同矿体倾角条件下应力和位移分析

为了研究不同矿体倾角条件下采场顶板应力及位移分布特征，现分别取采场倾角δ为0°，10°，20°和30°，其他参数保持不变，且采场内边界支护力设为0.不同倾角条件下采场顶板应力及位移解析解分布曲线如图7所示，同时添加了FLAC模拟结果作为对比.

如图7所示，不同矿体倾角条件下，采场顶板应力和位移的FLAC数值解和复变解析解在量值和变化规律方面基本吻合.在应力分布方面，对于水平矿体开采，采场顶板垂直应力和位移均沿采场中央呈现对称分布，其中垂直应力在0 MPa附近波动，容易发生拉应力诱发的拉伸破坏.随着采场倾角的增大，一方面，采场顶板垂直应力在量值上呈现逐渐增加的趋势，改善了顶板的应力环境；另一方面，垂直应力由水平采场条件下的对称分布逐渐向大倾角条件下的非对称分布特征演变，且这种现象随着倾角的增加而不断加剧.具体表现为顶板应力释放程度较大的区域向采场中上部转移，而在采场顶板下部区域应力急剧增加.在位移分布方面，水平采场顶板垂直位移较大，峰值位移位于采场顶板中央位置，约为12 mm.随着采场倾角的逐渐增大，采场顶板垂直位移逐渐减小，顶板垂直位移峰值逐渐向下部区域转移，呈现出非对称分布形态.综上所述，矿体倾角对顶板应力和位移的影响较大，近水平矿体开采顶板卸压程度和下沉量较大.随着矿体倾角的增大，采场顶板卸压程度和下沉量逐渐降低.

4.2 不同宽高比条件下应力和位移分析

为了分析不同宽高比条件下倾斜采场顶板应力及位移分布特征，此处对比采场的宽高比分别为2∶1，3∶1，4∶1，5∶1，6∶1，7∶1和8∶1.保持采场的结构特征和岩体参数不变，采场倾角为20°，且采场内无支护应力.根据解析解计算公式，得到不同宽高比采场顶板区域应力及位移的分布特征，如图8所示.

不同采场宽高比条件下，倾斜采场顶板垂直应力均呈现出上部小、下部大的偏对称分布特征，即采场顶板上半部分区域卸压程度较大，这是由采场倾角引起的.随着倾斜矩形采场宽高比的逐渐增大，采场斜长逐渐增加，顶板垂直应力逐渐减小，围岩应力降低程度逐渐增加.同时，采场顶板两侧肩角区域应力值随着采场斜长的增加而逐渐增大，容易诱发高应力剪切破坏.对于采场顶板位移分布特征，不同采场宽高比条件下顶板垂直位移分布规律并未发生变化，均呈现近似U形分布，如图8b所示.随着采场宽高比的增加，顶板垂直位移在量值上逐渐增加，这是由于采场顶板卸压程度增大引起的.因此，当缓倾斜采场开采长度较大时，顶板容易发生因大变形诱发的破坏.

根据解析计算和数值模拟的应力和位移结果，采场宽高比小于3∶1条件下，解析公式求解的顶板垂直应力存在上下波动的现象；当宽高比大于3∶1时，顶板应力解析解与数值模拟解具有较好的吻合度.与应力偏差不同之处在于，顶板位移解析解在不同宽高比条件下均与数值模拟结果具有较好的吻合度.因此，本文的研究结果在宽高比小于3∶1的采场围岩应力分析上稍有不足.

4.3 支护力对顶板应力和位移分布影响

根据上述分析，倾斜采场在大宽高比条件下，顶板容易发生拉应力诱发的破坏.研究表明，通过在采场内架设人工矿柱等支护体能够改善采场围岩的应力环境，限制围岩的位移，从而提高采场顶板的稳定性.为了分析支护力大小对采场顶板应力和位移分布的影响，设置采场顶底板均匀支护力分别为0~4 MPa，间隔为1 MPa.此外，矩形采场的宽、高尺寸为24 m×4 m，远场边界应力和岩体力学参数与前文保持一致，汇总数值模拟和解析计算结果如图9所示.可以看出，不同支护力作用下，采场顶板应力和位移在量值和变化规律方面与FLAC结果基本一致，偏差较小.

由图9可以看出，在不同支护力作用下，采场顶板应力和位移分布特征并未发生变化，但是在量值上变化较大.随着对采场顶底板支护力的增大，顶板垂直应力逐渐增大，其中垂直应力最小值由-0.07 MPa增大至-3.62 MPa.因此，对采场围岩施加支护力改善了顶板围岩的应力状态，降低了产生拉应力破坏的风险.此外，不同支护力条件下顶板垂直位移近似呈对称分布形态.随着采场顶底板支护应力从0 MPa逐渐增大至4 MPa，顶板围岩垂直变形量由11.0 mm逐渐降低至5.8 mm.采场顶板位移随着支护力的增大逐渐减小，避免了围岩因过大变形而诱发垮塌等地压灾害.

5 结 论

1） 通过增加保角映射函数的项数并引入修正系数，提出了大宽高比条件下矩形采场映射函数表达式.基于弹性力学复变函数理论推导了大宽高比倾斜矩形采场在双向不等压远场边界应力和均匀支护力作用下围岩应力和位移的解析解，其与FLAC数值模拟结果相对偏差小于5%.

2） 不同矿体倾角、采场宽高比和支护力条件下，顶板应力位移解析解与数值模拟结果在量值和演化规律方面均具有较好的一致性.采场倾角对围岩应力及位移分布影响较大，随着倾角的逐渐增大，倾斜矩形采场顶板应力和位移呈现出越来越剧烈的非对称分布特征.

3） 在矿体倾角较小、宽高比较大的条件下，采场顶板卸压程度和垂直位移量较大，容易发生拉应力诱发的垮塌灾害.随着矿体倾角的增加或宽高比的减小，采场顶板应力逐渐增大，垂直位移逐渐减小.支护体对围岩施加支护力能够显著改善围岩的应力环境，限制围岩的变形.

参考文献

[1] Yang K Y， Chen C X， Xia K Z，et al.Stability analysis and failure mechanism of the hanging wall in a sublevel caving mine based on microseismic monitoring［J］.Bulletin of Engineering Geology and the Environment，2024，83：259.

[2] Xie P S， Luo Y， Wu Y P， Gao X C，et al.Roof deformation associated with mining of two panels in steeply dipping coal seam using subsurface subsidence prediction model and physical simulation experiment［J］.Mining，Metallurgy & Exploration，2020，37（2）：581-591.

[3] 施耀斌.缓倾斜页岩矿床上向条带回采采场围岩变形特征及破坏判据［D］.武汉：武汉科技大学，2022.

[4] Shi Yao‑bin.Deformation characteristics and failure criterion of stope surrounding rock in upward strip mining of gently inclined shale deposit［D］.Wuhan：Wuhan University of Science and Technology，2022.

[5] Ravandi E G， Rahmannejad R.Wall displacement prediction of circular，D shaped and modified horseshoe tunnels in non‑hydrostatic stress fields［J］.Tunnelling and Underground Space Technology，2013，34：54-60.

[6] Abbo A J， Wilson D W， Sloan S W，et al.Undrained stability of wide rectangular tunnels［J］.Computers and Geotechnics，2013，53：46-59.

[7] Zheng M Z， Li S J， Feng Z J，et al.Three‑dimensional stress variation characteristics in deep hard rock of CJPL-II project based on in‑situ monitoring［J］.International Journal of Mining Science and Technology，2024，34（2）：179-195.

[8] Wang H N， Zeng G S， Jiang M J.Analytical stress and displacement around non‑circular tunnels in semi‑infinite ground［J］.Applied Mathematical Modelling，2018，63：303-328.

[9] 肖建清，冯夏庭，张腊春，等.均匀地应力场下圆形隧道静态弹塑性解析方法［J］.岩石力学与工程学报，2013，32（sup2）：3466-3477.

[10] Xiao Jian‑qing， Feng Xia‑ting， Zhang La‑chun，et al.Static elastoplastic analytical method of circular tunnel under uniform geostress field［J］.Chinese Journal of Rock Mechanics and Engineering，2013，32（sup2）：3466-3477.

[11] 覃爱民，骆汉宾.基于Mohr-Coulomb准则寒区隧道围岩应力弹塑性解析［J］.地下空间与工程学报，2018，14（2）：395-402.

[12] Qin Ai‑min， Luo Han‑bin.Analytical elasto‑plastic solution for stress of surrounding rock based on Mohr-Coulomb criterion［J］.Chinese Journal of Underground Space and Engineering，2018，14（2）：395-402.

[13] Muskhelishvili N I.Some basic problems of the mathematical theory of elasticity［M］.Groningen：Noordhoof Ltd，1954.

[14] 李元辉，刘明飞，李坤蒙，等.矩形采场支护作用下围岩周边应力复变函数解析解［J］.金属矿山，2021（4）：20-26.

[15] Li Yuan‑hui， Liu Ming‑fei， Li Kun‑meng，et al.Analytical solution of complex function of surrounding rock stress under supporting action of rectangular stope［J］.Metal Mine，2021（4）：20-26.

[16] Huo H， Bobet A， Fernández G，et al.Analytical solution for deep rectangular structures subjected to far‑field shear stresses［J］.Tunnelling and Underground Space Technology，2006，21（6）：613-625.

[17] Zhao G P， Yang S L.Analytical solutions for rock stress around square tunnels using complex variable theory［J］.International Journal of Rock Mechanics and Mining Sciences，2015，80：302-307.

[18] Lei G H， Ng C W W， Rigby D B.Stress and displacement around an elastic artificial rectangular hole［J］.Journal of Engineering Mechanics，2001，127（9）：880-890.

[19] Ng C W W， Lei G H.An explicit analytical solution for calculating horizontal stress changes and displacements around an excavated diaphragm wall panel［J］.Canadian Geotechnical Journal，2003，40：780-792.

[20] 李坤蒙，李元辉，王者超，等.硬岩采场新型预应力膨胀支柱的研发及应用［J］.东北大学学报（自然科学版），2021，42（2）：213-219.

[21] Li Kun‑meng， Li Yuan‑hui， Wang Zhe‑chao，et al.Development and application of a new pre‑stressed expandable pillar in hard rock stope［J］.Journal of Northeastern University （Natural Science），2021，42（2）：213-219.

[22] Wang M B， Li S C.A complex variable solution for stress and displacement field around a lined circular tunnel at great depth［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2009，33（7）：939-951.

[23] Wang H N， Jiang M J， Zhao T，et al.Viscoelastic solutions for stresses and displacements around non‑circular tunnels sequentially excavated at great depths［J］.Acta Geotechnica，2019，14（1）：111-139.

[24] Gao X， Wang H N， Jiang M J.Analytical solutions for the displacement and stress of lined circular tunnel subjected to surcharge loadings in semi‑infinite ground［J］.Applied Mathematical Modelling，2021，89：771-791.

[25] Dong X J， Karrech A， Qi C C，et al.The critical behaviour of finite thickness lining systems in tunnels［J］.European Journal of Environmental and Civil Engineering，2022，26（8）：3313-3330.

[26] Savin G N.Stress concentration around holes［M］.New York：Pergamon Press，1961.

[27] Yang S X， Huang L Y， Xie F R，et al.Quantitative analysis of the shallow crustal tectonic stress field in China mainland based on in situ stress data［J］.Journal of Asian Earth Sciences，2014，85：154-162.

基金资助

国家重点研发计划项目(2022YFC2903801)

国家重点研发计划项目(2022YFC2903804)

国家青年科学基金资助项目(52004054)

国家自然科学基金资助项目(52274115)