0 引言

采场结构的稳定性受多种因素影响，其稳定性的判据目前还没有统一的认识^[1-2]. 从单一的应力、位移角度很难综合评价采场结构的稳定性，冯夏庭等^[3]提出地下工程稳定性的综合集成智能分析与动态设计优化的新思路，并获得了成功应用. 能量法由于可以避免研究工程结构失稳破坏的中间复杂受力、应变过程，因而受到国内外一些学者的青睐，他们通过能量原理研究和探索了地下工程结构失稳破坏的相关问题，并取得了一定的成果^[4-9].

人工矿柱作为地下矿山开采中一种重要的采场结构，一方面可充分回收矿产资源，另一方面可人为改变采空区围岩应力环境^[10-11]，支撑上覆岩体和控制围岩变形.针对原生矿柱(残留矿柱)稳定性方面的研究，已经取得了一定的研究成果：唐春安等^[12]运用自主研发的岩石破裂过程分析程序(RFPA),对矿柱的变形与破坏过程及其声发射模式进行了研究；潘岳等^[13]研究了非对称开采时矿柱受偏压时的失稳破坏特征；赵奎^[14]分析了矿柱纵波速度与所受压力之间的定量关系以及矿柱破裂区分布特点，建立了基于声波测试的矿柱稳定性模糊推理系统，能够更加全面地反映矿柱稳定状态和损伤情况；周科平等^[15]对采场矿柱进行强度折减，采用矿柱塑性破坏区贯通作为相邻采场整体失稳的判据，利用矿柱强度折减的ANN-GA模型在某铜矿取得了满意效果；杨永杰等^[16]从蠕变试验的角度研究了石膏矿柱长期的稳定性问题；杨宇江和李元辉^[17]将加卸载响应比(LURR)理论引入矿柱动力稳定性研究中，对非均质矿柱模型在动力载荷循环作用下的破坏过程进行全时程动力分析. 刘洪强等^[18]、Zhou等^[19]、徐文彬等^[20]一些学者在原生矿柱稳定性方面也取得了非常有益的研究成果.

目前对人工矿柱在金属矿山地下开采支护的研究成果相对较少^[21-24]，且大多为人工矿柱的施工工艺、现场应用及定性分析方面研究. 近来，曹帅等^[25]基于弹性力学平面应变基本假设，建立阶段嗣后胶结充填体矿柱力学模型并进行理论求解；陈庆发等^[26]基于突变理论建立了金属矿山人工矿柱突发失稳的尖点突变模型，得出了满意的结果.因此，本文根据金属矿山人工矿柱尺寸和自身材料特性，从力学和能量守恒的角度研究其稳定性问题.

1 人工矿柱破坏模式及力学模型

由于金属矿山人工矿柱受多种因素的影响，其强度往往很难达到设计的要求，易发生失稳破坏. 在覆岩的载荷作用下，人工矿柱破坏模式与压杆失稳破坏模式相似，但由于覆岩、人工矿柱及载荷作用，其失稳机理比压杆更复杂. 在很多试验中，纵向劈裂破坏模式较多^[27].

由于人工矿柱与上覆岩体之间有较大的摩擦力 (图 1(a))，导致在垂直压力作用下人工矿柱顶面发生相对水平位移较小，很少产生较大的膨胀变形，而人工矿柱中部则因周围无约束力而发生"鼓形"(泊松)劈裂膨胀破坏. 继续加载时，人工矿柱中部对称位置上将产生竖向裂纹，并随载荷增大而发生贯通，使矿柱胀裂分开，并单独承受上覆岩体重力. 但此时，由于之前人工矿柱已产生横向变形而不再单纯是轴心受力，而是产生了初始偏心距的偏心受压的若干小矿柱体. 本文将人工矿柱视为压杆失稳破坏的模式，以"鼓形"(泊松)形式研究其稳定性.

取破坏后的某个小矿柱进行能量突变分析，为便于分析，假设人工矿柱长度$l$，宽度$b$，高度$h$如图 1(b)所示.由图 1可知，该单元体主要发生弯曲变形，其顶端可简化为仅在轴向有变形，在其他方向不发生变形的定向支座.在人工矿柱底端，由于底板围岩的刚度远大于人工矿柱的刚度，可视为此处无角位移和线位移，可将其简化为固支约束.因此，由人工矿柱受外载荷的力学模型可以简化，如图 2所示.

2 人工矿柱能量极限状态方程

根据前述破坏模式及力学模型分析，人工矿柱的边界条件为

$x(0) = 0 ,\quad x(h) = 0 $

(1)

${x}'(0) = 0 ,\quad {x}'(h) = 0 $

(2)

由以上二式满足位移边界条件^[28]，可设矿柱的挠度曲线方程为

$ y = a_1 \Big (1 - \cos \dfrac{2\pi x}{h} \Big) + b_1 \Big (1 - \cos \dfrac{6\pi x}{h} \Big) $

(3)

式中，$a_1 $为$x$方向挠度系数，$b_1$为$y$方向挠度系数(为提高计算准确程度，上式取前两项^[28]，$a_1 $和$b_1$为自由变量，在岩土工程方面其取值范围0.001$\sim$0.004^[29])；$h$为人工矿柱单元的高度，m.

由能量原理知^[30]，人工矿柱总能量$\varPi$为

$ \varPi = - W_P - W_G + U $

(4)

式中，$W_P $为外力$P$所做的功，J；$W_G $为矿柱自重势能，J；$U$为矿柱的弯曲应变能，J.

作为小变形体的人工矿柱，由于在尺寸上，矿柱的高度大于矿柱宽度，在承受单轴压缩时，常规的破坏是"鼓形"(泊松)破坏，这里仅考虑其弯曲应变能$U$为

$U = {1 \over 2}\int_0^{{\varphi _{\max }}} M (x){\rm{ }}\varphi = {1 \over 2}\int_0^h M (x) \cdot {{M(x)} \over {EI}}ds = {{EI} \over 2}\int_0^h {{{(y'')}^2}} ds$

(5)

式中，$I$为单个人工矿柱的惯性矩，$I= bh^3/12$，m$^{4}$；$E$为人工矿柱材料的弹性模量，MPa；$M(x)$为$x$高度处的截面弯矩，N$\cdot $m；$d s$为微弧段长，$d s = \sqrt {1 + ({y}')^2} d x$，rad.

将$d s = \sqrt {1 + ({y}')^2} d x$按泰勒公式展开

$ \sqrt {1 + ({y}')^2} \approx 1 + \dfrac{1}{2}({y}')^2 $

则

$ d s = \left[{1 + \dfrac{1}{2}({y}')^2} \right] d x $

(6)

将式(6)代入式(5)得

$ U = \dfrac{EI}{2}\int_0^h ({y}")^2\left[{1 + \dfrac{1}{2}({y}')^2} \right] d x $

(7)

由于实际问题中，$\theta $ ($\theta = {y}' = \tan \theta $)较小，并$\theta^2$ ($\theta ^2 = {y}'^2$)与1相比为高阶微量，可略去不计，简化后

$ U = \dfrac{EI}{2}\int_0^h ({y}")^2 d x $

(8)

将式(3)分别进行一阶求导与二阶求导得

$ \left. {y}' = \dfrac{2\pi a_1 }{h}\sin \dfrac{2\pi x}{h} + \dfrac{6\pi b_1 }{h}\sin \dfrac{6\pi x}{h} \\ {y}" = \dfrac{4\pi ^2a_1 }{h^2}\cos \dfrac{2\pi x}{h} + \dfrac{36\pi ^2b_1 }{h^2}\cos \dfrac{6\pi x}{h} \!\!\right\} $

(9)

将式(9)代入式(8)得

$\eqalign{ & U = {{EI} \over 2}\int_0^h ( {{4{\pi ^2}{a_1}} \over {{h^2}}}\cos {{2\pi x} \over h} + {{36{\pi ^2}{b_1}} \over {{h^2}}}\cos {{6\pi x} \over h}{)^2}dx = \cr & {{4EI{\pi ^4}} \over {{h^3}}}(a_1^2 + 81b_1^2) \cr} $

(10)

令$\lambda = \dfrac{1}{2}\int_0^h ({y}')^2 d x$，并将式(9)代入得

$\eqalign{ & \lambda = {1 \over 2}\int_0^h ( {{2\pi {a_1}} \over h}\sin {{2\pi x} \over h} + {{6\pi {b_1}} \over h}\sin {{6\pi x} \over h}{)^2}dx = \cr & {{{\pi ^2}} \over h}(a_1^2 + 9b_1^2) \cr} $

式中，$\lambda$ 为与$P$相应的广义位移.

则外载荷$P$所做的功

$ W_P = P\lambda = \dfrac{P\pi ^2}{h}(a_1^2 + 9b_1^2 ) $

(11)

矿柱自重势能为

$ W_G = \dfrac{1}{2}mgh = 0.5\rho gsh^2 $

(12)

式中，$\rho $为人工矿柱的密度，kg/m$^{3}$；$s$为人工矿柱的水平截面积，m$^{2}$；$g$为重力加速度，m/s$^{2}$.

将式(10) $\sim$式(12)代入式(4)得

$ \varPi = - W_P - W_G + U=\\ \dfrac{4EI\pi ^4}{h^3}(a^2_1 + 81b_1^2 ) - \\ \dfrac{P\pi ^2}{h}(a_1^2 + 9b_1^2 ) - 0.5\rho gsh^2 $

(13)

此式即为结构功能(势能)函数的数学表达式. 由突变理论的相关知识可知^[31-32]，用势函数$\varPi (x)$的空间曲面图形来表示如图 3.

图 3(a)所示的是一个具有光滑折痕的曲面，其上一点代表所研究系统的一种平衡状态. 在该曲面上叶和下叶上，满足${\varPi}"(x) >0$，即系统势能取极小值，则平衡状态是稳定的；在曲面的中叶上${\varPi}"(x) ＜0$，即系统势能取极大值，则平衡状态是非稳定的.

在曲面上叶和下叶与中叶的交界，即图 3(b)所示的光滑折痕$OA,OB$上，${\varPi}"(x)= 0$，即为临界状态.

如前所述，当${\varPi}"(x) ＞ 0$时，平衡状态是稳定的；当${\varPi }"(x) ＜0$ $时，系统状态是非稳定的；当$ {\varPi }"(x) = 0$时，系统处于临界状态，此时的外力$P$为临界载荷，记为$P_{\rm cr} $.

令$\varPi = - P\lambda - W_G + U = 0$得

$ P = \dfrac{4EI\pi ^2}{h^2}\cdot \dfrac{a_1^2 + 81b_1^2 }{a_1^2 + 9b_1^2 } - \dfrac{\rho gsh^3}{2\pi ^2}\cdot \\ \dfrac{1}{(a_1^2 + 9b_1^2 )} $

(14)

从上式中可看出，$\dfrac{a_1^2 + 81b_1^2 }{a_1^2 + 9b_1^2 } ≥1$，当且仅当$b_1 = 0$时，取最小值1.

外载荷$P$的极大值，即临界载荷

$ P_{\rm cr} = \dfrac{4EI\pi ^2}{h^2} - \dfrac{\rho gsh^3}{2\pi ^2a_1^2 } $

(15)

式中，$h$为人工矿柱高度，m；$b$为人工矿柱宽度，m；$g$为重力加速度，N/kg；$I$为单个人工矿柱的惯性矩，$I =bh^3/12$，N$\cdot $m；$a_1$为$x$方向挠度系数，查阅《建筑结构静力计算手册》中的有关内力系数和挠度系数，0.001 527 5.

相应的，临界应力

$ \sigma _{\rm cr} = \dfrac{p_{\rm cr} }{A} = \dfrac{ 4EI\pi ^2}{h^2bl} - \dfrac{\rho gh^3}{2\pi ^2a_1 ^2} $

(16)

3 人工矿柱设计稳定性分析

根据非线性科学的理论^[33]，该系统的定态方程为

$ \left.\begin{array}{l} {\dfrac{\partial \varPi }{\partial a_1 } = 0} \\ {\dfrac{\partial \varPi }{\partial b_1 } = 0} \end{array} \right \} $

(17)

将式(13)代入得

$ \left.\begin{array}{l} {a_1 \Big(\dfrac{8EI\pi ^4}{h^3} - \dfrac{2\pi ^2}{h}P \Big) = 0} \\ {b_1 \Big(\dfrac{36EI\pi ^4}{h^3} - \dfrac{\pi ^2}{h}P \Big) = 0} \end{array} \!\! \right \} $

(18)

解此方程组可得

$①\left\{ {\matrix{ {{a_1} = 0} \cr {{b_1} = 0} \cr } } \right.②\left\{ {\matrix{ {{a_1} = 0} \cr {P = {{36EI{\pi ^2}} \over {{h^2}}}} \cr } } \right.③\left\{ {\matrix{ {P = {{4EI{\pi ^2}} \over {{h^2}}}} \cr {{b_1} = 0} \cr } } \right.$

当第①组解代入式(3)有$y =0$，与人工矿柱在力$P$的作用实际发生了数值不为零的位移相矛盾；通过比较②和③组解，易知第③组解中的$P$值小于②组解中的$P$值，故可根据载荷由小到大逐次增大知，临界载荷$P_{\rm cr} = \dfrac{4EI\pi ^2}{h^2} - \dfrac{\rho gsh^3}{2\pi ^2a_1^2}$，此时系统处于临界状态.

再由式(14)讨论外载荷$P$变化过程对结构系统

稳定性的影响，即$P = \dfrac{4EI\pi ^2}{h^2}\cdot \dfrac{a_1^2 + 81b_1^2 }{a_1^2 + 9b_1^2 } -\dfrac{\rho gsh^3}{2\pi ^2}\cdot \dfrac{1}{ a_1^2 + 9b_1^2 }$，可知外载荷$P$是关于$a_1 $和$b_1 $的函数，故

可用其偏导数来研究外载荷$P$的增减趋势.

$ \left. \dfrac{\partial P}{\partial a_1 } = - \dfrac{a_1 (576\pi ^4EIb_1^2 - \rho gsh^5)}{\pi ^2h^2(a_1^2 + 9b_1^2 )^2} \\ \dfrac{\partial P}{\partial b_1 } = \dfrac{9b_1 (64EI\pi ^4a_1^2 + \rho gsh^5)}{h^2\pi ^2(a_1^2 + 9b_1^2 )^2} \!\!\right\} $

(19)

根据式(19)可得载荷力(压力)$P$的变化规律：

当$a_1 ＞ 0$时，$\dfrac{\partial P}{\partial a_1 } <0$，即压力$P$随$a_1 $增大而减小；当$a_1 ＜0$时，$\dfrac{\partial P}{\partial a_1 } ＞ 0$，即压力$P$随$a_1 $增大而增大.

同理，当$b_1 ＞ 0$时，压力$P$随$b_1 $增大而增大；当$b_1 <0$时，压力$P$随$b_1 $增大而减小.

由材料力学可知，矿柱产生突变时,必须满足 $P ＞ P_{\rm cr} $.

考虑式(13)及式(18) 可知：

当$P < P_{\rm cr} = \dfrac{4EI\pi ^2}{h^2} - \dfrac{\rho gsh^3}{2\pi ^2a_1^2 }$时，$\varPi$随$\left| {a_1 } \right|$的增大而增大，即系统处于稳定平衡状态，势函数$\varPi $取极小值；当$P ＞ P_{\rm cr}$时，势函数$\varPi$随$\left| {a_1 } \right|$的增大而减小，即$\varPi $取极大值，系统此时处于非稳定的平衡状态.

在外载荷(压力)$P$从0增大至$P_{\rm cr} $(极小值)的过程中，系统$\varPi $取极小值，始终处于稳定平衡状态；$P >P_{\rm cr} $时系统势函数$\varPi $开始取极大值，这之间发生了稳定状态的跳跃，即平衡状态的突变.

综上可知，欲使人工矿柱系统处于稳定平衡状态，则作用在人工矿柱上的外载荷小于临界载荷，即$P < P_{\rm cr} = \dfrac{4EI\pi ^2}{h^2} - \dfrac{\rho gsh^3}{2\pi ^2a_1^2 }$，也即是要确保人工矿柱承受的应力值小于其临界应力$\sigma _{\rm cr} =\dfrac{4EI\pi ^2}{h^2bl} - \dfrac{\rho gh^3}{2\pi ^2a_1^2}$. 在外载荷不变的情况下要使人工矿柱应力减小，可采取增大人工矿柱尺寸面积，或增大人工矿柱强度，如采取高强度等级水泥、加筋等措施.

4 工程实例

根据上述针对人工矿柱承受外载荷的稳定性分析，为了能将该理论研究应用到矿山工程实践中，作者所在课题组对焦冲金矿^[22]采用人工混凝土矿柱代替原生矿柱的尺寸参数进行了优化设计，并对支护效果进行了研究.

该矿$-430$ m中段矿体走向长度85 m，经某设计单位设计人工矿柱4个，每个矿柱宽8 m，矿房5个，矿房宽度为10.6 m.为了降低充填成本，同时确保采场结构的稳定，课题组根据前述构建的人工矿柱稳定性能量极限状态方程，结合人工矿柱弹性模量及高宽比等因素的临界应力函数关系式对该矿山人工矿柱进行了优化设计.经理论研究和室内相关实验，最终确定充填体弹性模量8.5 GPa，容重22 kN/m$^{3}$，设计该中段矿柱为4个，每个矿柱宽6.5m，矿房个数为5个，矿房宽度为11.8 m.根据课题组所设计的开采方案，为了监测矿山覆岩对人工矿柱的压力，以便实时观测人工矿柱是否能有效支撑覆岩载荷，在充填人工矿柱之前在分条上山底部布设混凝土压力计(图 4)，经过数月的监测其最大垂直应力值为$P_{\hbox{\tiny测应力}} =0.78$ MPa (图 5).从压力计监测曲线图来看，所设计的人工矿柱在矿房回采初期的两个月内，人工矿柱承受上覆岩层的压力持续增大, 且上升速度较快；当矿房回采完毕后，压力值变化曲线逐渐趋于平缓，表明人工矿柱在开采扰动结束后，上覆岩层的载荷分布也趋于稳定，人工矿柱承受的应力也逐渐平稳.

根据矿柱和矿房的布置，每个矿柱承载左右相邻每个矿房一半面积的上覆岩层载荷，一个矿柱应该承载的顶板面积

$ s_1 = \dfrac{b}{2} \left( {b + w} \right) $

(20)

式中，$w$为矿房宽度，m；即每个矿柱实际承载载荷 $P = P_{\hbox{\tiny测应力}} s_1 = 4.64\times 10^4$ kN.

将上述参数代入式(15)，可得人工矿柱的临界载荷

$ P_{\rm cr} = 3.90 \times 10^5 {\rm kN} $

(21)

从上述现场实测数据可知，课题组优化设计的人工矿柱所承受覆岩载荷为464 MN，根据前述所建立的人工矿柱临界应力函数关系式所计算出的人工矿柱临界载荷为390 MN，实测人工矿柱承受载荷约为理论设计临界载荷的十分之一，远小于其临界载荷的支撑力，课题组所设计的人工矿柱是稳定的，符合工程安全规范要求，能够满足矿山安全生产的需要.为了进一步优化人工矿柱尺寸，减少充填成本，将该中段其他剩余人工矿柱的宽度又适当减小后进行施工.截止到2014年12月底，通过5年多的矿山生产实践，$-430$ m中段人工矿柱较稳定，未发现大的变形和剥落；此中段覆岩也未出现大的沉降及位移.

通过现场实验和多年的矿山生产实践，表明上述所建立的能量极限状态方程可以作为金属矿山人工矿柱设计的力学依据.

5 结论

(1) 将人工矿柱受力环境进行分析简化，根据能量守恒原理得出考虑外载荷做功、自重势能和应变能条件下的人工矿柱总能量方程，推导出了结构功能(势能)函数的表达式. 根据突变理论的相关知识，建立了在考虑弹性模量、高宽比等参数条件下的能量极限状态方程，根据势函数的稳定判据得出了外载荷的临界应力函数关系式.

(2) 根据临界载荷关系，从力学角度定量评判了人工矿柱在外载荷变化过程对其稳定性的影响. 当外力$P$从0增大至$P_{cr} $ (极小值)的过程中，系统$\varPi $取极小值，始终处于稳定平衡状态；$P ＞ P_{\rm cr} $时系统$\varPi$开始取极大值，期间发生了稳定状态的跳跃，即平衡状态的突变.欲使人工矿柱始终处于稳定状态，应采取一定措施，确保人工矿柱所承受的外载荷小于其临界载荷.

(3) 根据理论研究结果，对优化设计后的焦冲金矿人工矿柱承受载荷现场实测，其实测结果约为理论设计临界载荷的十分之一，远小于其临界载荷的支撑力. 通过5年多的矿山生产实践，人工矿柱及覆岩都较稳定，能够满足矿山安全生产的需要. 表明所建立的人工矿柱能量极限状态方程，可为地下金属矿山人工矿柱设计提供一定的力学参考依据和理论指导.