Safety performance of unmanned monorail cranes in mines key parameters identification research
-
摘要:
单轨吊在复杂深部矿井环境的辅助运输系统中具有不可替代的作用,由于目前无法有效精确识别单轨吊的荷载质量和轨道坡度,直接影响了运输安全性能及能量高效利用。因此,笔者提出了矿井无人驾驶单轨吊安全性能关键参数识别方法。针对单轨吊结构特性和轨道运输特点,对具有强耦合关系的荷载质量和轨道坡度建立了纵向动力学模型;基于运行数据和带有动态遗忘因子的递推最小二乘算法(DFFRLS)对纵向动力学模型参数进行实时在线识别,实现荷载质量和轨道坡度的精准解耦;并基于解耦的纵向动力学模型和识别的模型参数,动态修正当前的荷载质量识别值,以消除误差,完成荷载质量的高精度识别;由识别的纵向动力学模型参数、运行数据,应用Sage-Husa自适应扩展卡尔曼滤波算法(AEKF)对系统噪声协方差和误差协方差进行动态更新,滤除环境噪声干扰,实时调节和修正当前轨道坡度值,保证轨道坡度识别的精准性。在多种工况下,仿真与实际应用表明,基于DFFRLS-AEKF方法的荷载质量识别值与实际值的误差在3.2%以内,运行轨道坡度识别值与实际值的误差在5.3%以内。该方法可实现无人驾驶单轨吊安全性能关键参数的实时精准获取,有效减少无人驾驶单轨吊安全事故的发生,显著提升无人驾驶单轨吊的能量高效利用。
-
关键词:
- 无人驾驶单轨吊 /
- 安全性能 /
- 荷载质量 /
- 自适应扩展卡尔曼滤波(AEKF) /
- 识别误差
Abstract:Monorail cranes have an irreplaceable role in the auxiliary transport system in complex deep mine environment, which directly affects the transport safety performance and efficient energy utilization due to the current inability to effectively and accurately identify the load quality and track slope of monorail cranes. For this reason, this paper proposes a method for identifying key parameters of safety performance of unmanned monorail cranes in mines based on DFFRLS-AEKF. Based on the structural characteristics of the monorail crane and the characteristics of track transportation, a longitudinal dynamics model is established for the load mass and track slope with strong coupling relationship; the parameters of the longitudinal dynamics model are identified online in real time based on the operational data and the recursive least squares algorithm with dynamic forgetting factor to achieve the accurate decoupling of the load mass and track slope; and based on the decoupled longitudinal dynamics model and the identified model parameters, the current load mass identification is dynamically modified. Based on the decoupled longitudinal dynamics model and the identified model parameters, the current load quality recognition value is dynamically corrected to eliminate the error and complete the high-precision recognition of load quality; from the identified longitudinal dynamics model parameters and operation data, the Sage-Husa adaptive extended Kalman filter algorithm is applied to dynamically update the system noise covariance and error covariance, filter out the environmental noise interference, adjust and correct the current track slope value in real time, and ensure the accuracy of track slope recognition. Accuracy. The simulation and practical application show that the error between the load mass recognition value and the actual value is within 3.2% and the error between the running track slope recognition value and the actual value is within 5.3% under various working conditions. The method achieves real-time and accurate acquisition of key parameters of the safety performance of the unmanned monorail crane, effectively reduces the occurrence of safety accidents of the unmanned monorail crane, and significantly improves the efficient utilization of energy of the unmanned monorail crane.
-
0. 引 言
煤炭资源是我国能源构成中至关重要的部分,近些年诸多矿区资源开采逐渐由浅部转向深部[1-3],煤炭安全高效开采是保障该资源正常供应的首要问题。采煤过程中会导致原岩应力状态失去平衡[4-6],继而使得底板岩层受到一定程度的扰动,底板产生的裂隙可能贯穿至含水层,直接影响采煤工作面的安全,甚至造成严重的生命财产损失,底板破坏深度的准确判断[7-8]是矿井防治底板突水的有效手段。
目前国内外学者在煤层底板破坏深度研究方面取得的成果丰硕,理论方面以塑性滑移线场为主。鲁海峰等[9]运用瑞典条分法搜出危险滑面并获取稳定系数和滑面最大深度,为正确使用弹塑性理论求解底板最大破坏深度提供了依据;张金才等[10]根据塑性滑移线场理论构建了底板受力模型,给出了采面底板破坏形态;孙建[11] 以塑性滑移线场为依托建立了倾斜煤层的底板破坏力学模型,分析了沿煤层倾斜方向底板三区破坏形态;李昂等[12-13]针对董家河煤矿带压水上采煤问题,基于塑性滑移线场理论给出了煤层底板破坏深度解析解。上述理论计算方法中利用塑性滑移线场理论进行底板破坏深度分析时均将底板视为单一岩体,然而底板是由不同岩性岩层构成的复合结构[14-16],且塑性滑移线场理论是基于岩层平均内摩擦角和层厚等参数去求解底板最大破坏深度,若按照传统的单一岩层底板塑性滑移线场理论计算会因为忽略各岩层不同内摩擦角和岩性而产生较大误差。
随着采深的日益增加,煤矿安全开采受强烈开采扰动和高承压水的影响日趋严重化与复杂化[17-20]。虽然李昂等[21-22]对于平煤某矿深部开采建立了3层复合底板塑性滑移线场力学计算模型,并对底板破坏深度解析式进行推导,但该计算过程繁杂、工作量大,部分工况计算结果因泰勒展开式的运用只能取近似值。如果底板隔水层厚度较大、岩性较为复杂,根据柱状把岩性相近或岩层厚度较薄的岩层划分为1层,摩擦角取其平均值;岩层较厚且岩性差异较大的岩层划分为另外层,但最多只能渐划为3层,适用性较差,难以满足深部开采中多岩性结构底板破坏深度的分析计算。
对于深部煤层开采迫切需要一种计算精度高、方便快捷的多岩性结构底板最大破坏深度计算方式。基于此,在3层复合结构底板的力学模型基础上,对多层结构底板破坏深度解析解进行推导。对各计算工况进行优化,提出了多层结构底板的5种计算工况,形成了计算逻辑流程判别图,利用Matlab语言开发了多层结构底板破坏深度计算系统V1.0软件,并通过应用分析部分验证了该计算系统的严谨性与精确性,对推动深部煤岩开采面底板水害防治技术有较大作用。
1. 多层塑性滑移线场理论分析
1.1 多层结构底板超前塑性破坏长度
考虑到现场实际底板岩层的复杂性,采用n层模型计算底板最大破坏深度及采面至最大破坏深度的水平距离。而底板最大破坏深度与底板超前塑性破坏长度密切相关,采用理论值与实测值相结合的方法确定底板超前塑性破坏长度。理论值采用经典塑性滑移线力学模型计算:
$$ x_0=\frac{H_2}{\theta^*}\left\{\left[\left(K\gamma H_1+\frac{C_1}{f_1}\right)\left(\dfrac{f_1}{C_1+f_1\sigma_{\mathrm{c}}}\right)\right]^{\tfrac{\theta^*}{2k_{\mathrm{p}}f_1}}-1\right\} $$ (1) 式中:x0为底板超前塑性破坏长度;H1为煤层的埋深; H2为煤层的采高;K为支承压力峰值系数(由现场实测得到);$ \gamma $ 为煤层覆岩的平均容重;$ \theta^* $为煤层的变形角;$ C_1 $为煤层与顶底板接触面上的黏聚力;f1为煤层与顶板界面间的摩擦因数;$ \sigma_c $为单轴压缩时煤体的极限抗压强度,$ \sigma_{\mathrm{c}}=\dfrac{2C\cos\ \varphi_{ }}{1-\sin\ \varphi_{ }} $,C为煤体的黏聚力;$\varphi $为内摩擦角;$ k_{\mathrm{p}}=\dfrac{1+\sin\ \varphi_{ }}{1-\sin\ \varphi_{ }} $。
1.2 多层结构底板计算破坏深度工况划分
深部煤岩开采面底板岩性不尽相同,使得计算工况也愈发复杂。以岩体主动极限区深度H0所处岩层位置作为判别依据,对多层结构底板进行5种工况划分。n层结构底板破坏深度计算简图如图1所示。
1.2.1 工况1
n=1,即将底板看作单一岩层,学者张金才[10]认为当煤层被采出后,采空区四周岩体所承受的垂直应力增大,当垂直应力作用范围的底板 (如图2中I区所示) 下岩体所受采动应力超过其最大强度时,底板岩体结构将发生塑性破坏,由于该岩体在法向方向上受到挤压作用,将会发生水平方向的拉伸,变形后的岩体压迫过渡区域(如图2中Ⅱ区所示)内的底板岩体,并将附加应力传递到该区域。过渡区域岩体受其压迫后,将会继续压迫被动区域内的底板岩体 (如图2中Ⅲ区所示)。其中,主动和被动区各由两条直线构成,过渡区滑移线分别由对数螺旋线和自O点的放射线组成。
螺线方程为
$$ r=r_0\mathrm{e}^{\theta\tan\ \varphi_1} $$ (2) 式中:r为O到对数螺旋线BC之间的距离;r0为OB之间的距离;$ \theta $为相邻底板岩层间的夹角,此时代表OB与OE之间的夹角;$ \varphi_{1} $为底板岩层的平均内摩擦角。
根据三角函数公式能求出第1段对数螺旋线的起点半径r0:
$$ OB = {r_0} = \frac{{{x_0}}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right)}} $$ (3) 然后根据得出的起始半径r0、r和θ值,经过推导得出最终的单层结构底板最大滑移破坏深度的计算公式为
$$ h_0=\frac{x_0}{2\cos\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)}\mathrm{e}^{\left(\tfrac{\pi}{4}+\frac{\varphi_1}{2}\right)\tan\ \varphi_1}\cos\ \varphi_1 $$ (4) 由图1几何关系知开采面至最大破坏深度水平距离OG为
$$ OG=\frac{x_0}{2\cos\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)}\mathrm{e}^{\left(\tfrac{\pi}{4}+\frac{\varphi_1}{2}\right)\tan\ \varphi_1}\sin\ \varphi_1 $$ (5) 此求解方法是将底板简化为单一岩层,但实际工程中底板岩体由多种岩性的岩体组成,岩层之间内摩擦角并不相同,采用单一岩层底板计算最大破坏深度存在较大误差。
1.2.2 工况2
H0<H1,即主动极限区深度小于首层岩层厚度时,因结构底板各内摩擦角不同会导致出现多条对数螺旋线半径,采用层数n确定待求对数螺旋线半径的条数,由图1几何关系知底板主动极限区深度H0:
$$ {H_0} = \frac{{{x_0}}}{2} \cdot \tan \left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2}} \right) $$ (6) 设O点作为2段对数螺旋线的旋转中心点,且交界处螺旋线半径保持不变,根据三角函数公式能求出第1段对数螺旋线的起点半径r0:
$$ OB = {r_0} = \frac{{{x_0}}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right)}} $$ (7) 第2段对数螺旋线起始半径r1由方程组(8)联立求解:
$$ \left\{\begin{gathered}r_1=\frac{H^1}{\sin\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1\right)} \\ r_1=r_0\mathrm{e}^{\theta_1\tan\ \varphi_1} \\ \end{gathered}\right. $$ (8) 式中:$ \theta_{1}\text { 为对数螺旋线半径 } r_{0} \text { 与 } r_{1} \text { 间夹角。化简为 } $
$$ H^1=r_0\mathrm{e}^{\theta_1\tan\ \varphi_1}\cdot\sin\left(\frac{\pi}{4}+\frac{\varphi_1}{2}+\theta_1\right) $$ (9) 对于式(8)中的未知数θ1采用泰勒展开求得近似解,根据泰勒展开式:
$$ \mathrm{e}^{\theta_1\tan\ \varphi_1}\approx1+\theta_1\tan\ \varphi_1 $$ (10) $$ \sin \left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2} + {\theta _1}} \right) \approx \left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2} + {\theta _1}} \right) - \dfrac{{{{\left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1}} \right)}^3}}}{6} $$ (11) 将式(10)和式(11)代入式(9),可得
$$ H^1\approx r_0\cdot(1+\theta_1\tan\ \varphi_1)\cdot\left[\begin{array}{*{20}{c}}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1\right)- \\ \dfrac{\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1\right)^3}{6}\end{array}\right] $$ (12) 展开得
$$ \begin{array}{*{20}{c}}\left(-\dfrac{1}{6}\mathrm{tan\ }\varphi_1\right)\theta_1^4+\left[-\dfrac{1}{6}-\dfrac{1}{2}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)\mathrm{tan}\ \varphi_1\right]\theta_1^3+ \\ \left[\mathrm{\mathrm{tan}}\ \varphi_1-\dfrac{1}{2}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)-\dfrac{1}{2}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)^2\mathrm{tan}\ \varphi_1\right]\theta_1^2+ \\ \left[\begin{array}{*{20}{c}}1-\dfrac{1}{2}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)^2+\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)\mathrm{tan}\ \varphi_1 \\ \\ -\dfrac{1}{6}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)^3\mathrm{tan\ }\varphi_1\end{array}\right]\theta_1+ \\ \left[\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}-\dfrac{H_1}{r_0}-\dfrac{1}{6}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}\right)^3\right]=0\end{array} $$ (13) 在式(13)中,令
$$ \begin{gathered}a=-\frac{1}{6}\tan\ \varphi_1,b=-\frac{1}{6}-\frac{1}{2}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)\tan\ \varphi_1, \\ c=\tan\ \varphi_1-\frac{1}{2}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)-\frac{1}{2}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)^2\tan\ \varphi_1, \\ d=1-\frac{1}{2}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)^2+\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)\tan\ \varphi_1- \\ \frac{1}{6}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)^3\tan\ \varphi_1, \\ e=\frac{\pi}{4}+\frac{\varphi_1}{2}-\frac{H_1}{r_0}-\frac{1}{6}\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)^3 \\ \end{gathered} $$ 则式(13)表示为
$$ a\theta _1^2 + b\theta _1^3 + c\theta _1^4 + d{\theta _1} + e = 0 $$ (14) 由一元四次方程的费拉里解法可知,式(14)存在4个解,考虑到所求解是对数螺旋线的夹角,因此实际取值需为正数,最后经验算,实际存在的解仅1个,即
$$ \begin{array}{*{20}{c}} {{\theta _1} = \dfrac{{ - b}}{{4a}} + \dfrac{1}{2}\sqrt {\dfrac{{{b^2}}}{{4{a^2}}} - \dfrac{{2c}}{{3a}} + {\text{\Delta }}} } -\\ { \dfrac{1}{2}\sqrt {\dfrac{{{b^2}}}{{2{a^2}}} - \dfrac{{4c}}{{3a}} - {\text{\Delta }} + \dfrac{{ - \dfrac{{{b^3}}}{{{a^3}}} + \dfrac{{4bc}}{{{a^2}}} - \dfrac{{8 d}}{a}}}{{4\sqrt {\dfrac{{{b^2}}}{{4{a^2}}} - \dfrac{{2c}}{{3a}} + {\text{\Delta }}} }}} } \end{array} $$ (15) 式中
$$ \begin{gathered}{\text{\Delta }} = \frac{{\sqrt[3]{2}{{\text{\Delta }}_1}}}{{3\alpha \sqrt[3]{{{{\text{\Delta }}_2} + \sqrt[2]{{ - 4{{\text{\Delta }}_1}^3 + {{\text{\Delta }}_2}^2}}}}}} + \\ \frac{{\sqrt[3]{{{{\text{\Delta }}_2} + \sqrt[2]{{ - 4{{\text{\Delta }}_1}^3 + {{\text{\Delta }}_2}^2}}}}}}{{3\sqrt[3]{2}\alpha }}\\ {{\text{\Delta }}_1} = {c^2} - 3bd + 12ae; \\ {{\text{\Delta }}_2} = 2{c^3} - 9bcd + 27a{d^2} + 27{b^2}e - 72ace \end{gathered} $$ 式中:α为GE与对数螺旋线r之间夹角。
从而计算出第2段螺旋线起始半径r1,将所求θ1,r1代入方程组(16):
$$ \left\{\begin{gathered}r_2=\dfrac{H'+H''}{\sin\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1+\theta_2\right)} \\ r_2=r_1\mathrm{e}^{\theta_2\tan\ \varphi_2} \\ \end{gathered}\right. $$ (16) 式中:θ2 为对数螺旋线半径r1与r2间夹角。
依次计算$ \theta_I,\theta_2,\cdots\theta_{\mathrm{\mathit{n}}-2},r_1,r_2,\cdots r_{\mathit{\mathrm{\mathit{n}}}-2} $,将所求$ \theta_I,\theta_2,\cdots\theta_{\mathrm{\mathit{n}}-2},r_{\mathit{\mathrm{\mathit{n}}}-2} $代入方程组(17):
$$ \left\{\begin{gathered}r_{\mathrm{\mathit{n}}-1}=\dfrac{H^1+H^2+...+H^{n-1}}{\mathrm{sin}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1+\theta_2+\theta_3+...\theta_{\mathrm{\mathit{n}}-1}\right)} \\ r_{\mathrm{\mathit{n}}-1}=r_{\mathrm{\mathit{n}}-2}\mathrm{e}^{\theta_{\mathit{\mathrm{\mathit{n}}}-1}\mathrm{tan}\ \varphi_{\mathrm{\mathit{n}}-1}} \\ \end{gathered}\right. $$ (17) 式中:θn–1 为对数螺旋线半径rn–2与rn–1间夹角。
第2段随后的螺旋线起始半径及螺旋线夹角$ r_{\mathrm{2}},\theta_{\mathrm{2}},r_{\mathrm{3}},\theta_{\mathrm{3}},\cdots,r_{\mathrm{\mathit{n}}-1},\theta_{\mathrm{\mathit{n}}-1} $的求解方法与r1相同,最终结果采用数学计算软件Mathematic进行求解。
在△OEG中
$$ OE=r=r_{\mathrm{\mathit{n}}-1}\mathrm{e}^{\theta\tan\ \varphi_{\mathrm{\mathit{n}}}} $$ (18) $$ GH=h=r\cos\ \alpha=r_{\mathrm{\mathit{n}}-1}\mathrm{e}^{\theta t\tan\ \varphi_{\mathrm{\mathit{n}}}}\cos\left(\begin{array}{*{20}{c}}\theta+\theta_1+\theta_2+\cdots+ \\ \theta_{\mathit{\mathrm{\mathit{n}}}-1}+\dfrac{\varphi_1}{2}-\dfrac{\pi}{4}\end{array}\right) $$ (19) 式中:$ \theta $为相邻底板岩层间的夹角,此时代表对数螺旋线 r与 rn–1之间夹角。
当$ \dfrac{\mathrm{d}h}{\mathrm{d}\theta}=0 $ 时,h 达到最大破坏深度h0,解得
$$ \theta=\varphi_{\mathrm{\mathit{n}}}+\frac{\pi}{4}-\frac{\varphi_1}{2}-\theta_1-\theta_2-\cdots\theta_{\mathrm{\mathit{n}}-1} $$ (20) 将式(17)、式(18)和式(20)的解代入式(19),得出多层结构底板最大滑移破坏深度的计算公式:
$$ h_0=r_{\mathrm{\mathit{n}}-1}\mathrm{e}^{\left(\varphi_{\mathrm{\mathit{n}}}+\tfrac{\pi}{4}-\tfrac{\varphi_1}{2}-\theta_1-\theta_2-\theta_3-...-\theta_{\mathrm{\mathit{n}}-1}\right)\mathrm{tan}\ \varphi\mathit{\mathit{\mathit{_{\mathrm{\mathit{n}}}}}}}\mathrm{cos\ }\varphi_{\mathrm{\mathit{n}}} $$ (21) 开采面至最大破坏深度水平距离OG:
$$ OG=r_{\mathrm{\mathit{n}}-1}\mathrm{e}^{\left(\varphi\mathit{\mathit{_{\mathrm{\mathit{n}}}}}+\tfrac{\pi}{4}-\tfrac{\varphi_1}{2}-\theta_1-\theta_2-\theta_3-...-\theta_{\mathrm{\mathit{n}}-1}\right)\mathrm{tan}\ \varphi_{\mathrm{\mathit{n}}}}\sin\ \varphi_{\mathrm{\mathit{n}}} $$ (22) 1.2.3 工况3
H0=H1,主动极限区深度等于首层岩层厚度,即主动极限区深度H0处于第1层岩体和第2层岩体交界处。由图1几何关系知底板主动极限区深度H0:
$$ {H_0} = \frac{{{x_0}}}{2} \cdot \tan \left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2}} \right) $$ (23) 设O点作为2段对数螺旋线的旋转中心点,且交界处螺旋线半径保持不变,根据三角函数公式能求出第1段对数螺旋线的起点半径r0:
$$ OB = {r_0} = \frac{{{x_0}}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right)}} $$ (24) 各段螺旋线起始半径及螺旋线夹角$ r_1,\theta_1, r_2,\theta_2,r_3,\theta_3,\cdots,r_{\mathrm{\mathit{n}}-3},\theta_{\mathit{\mathrm{\mathit{n}}}-3} $的求解方法与工况2相同,最终结果采用数学计算软件Mathematic进行求解。
将所求$ \theta_{\mathrm{1}},\theta_2,\cdots,\theta_{\mathrm{\mathit{n}}-3},r_{\mathit{\mathrm{\mathit{n}}}-3} $代入方程组(25):
$$ \left\{\begin{gathered}r_{\mathit{\mathrm{\mathit{n}}}-2}=\dfrac{H^1+H^2+...+H^{\mathrm{\mathit{n}}-1}}{\mathrm{sin}\left(\dfrac{\pi}{4}+\dfrac{\varphi_1}{2}+\theta_1+\theta_2+\theta_3+...\theta_{\mathrm{\mathit{n}}-2}\right)} \\ r_{\mathit{\mathit{\mathrm{\mathit{n}}}}-2}=r_{\mathit{\mathit{\mathit{\mathrm{\mathit{n}}}}}-3}\mathrm{e}^{\theta_{\mathrm{\mathit{n}}-2}\mathrm{tan}\ \varphi_{\mathit{\mathrm{\mathit{n}}}-1}} \\ \end{gathered}\right. $$ (25) 式中:θn-2 为对数螺旋线半径rn-3与rn-2间夹角。
在△OEG中
$$ OE=r=r_{\mathit{\mathrm{\mathit{n}}}-2}\mathrm{e}^{\theta\tan\ \varphi_{\mathrm{\mathit{n}}}} $$ (26) $$ \begin{gathered}\qquad GH=h=r\cos\ \alpha= \\ r_{\mathit{\mathit{\mathrm{\mathit{n}}}}-2}\mathrm{e}^{\theta t\tan\ \varphi_{\mathrm{\mathit{n}}}}\cos\left(\begin{gathered}\theta+\theta_1+\theta_2+\cdots+ \\ \theta_{\mathrm{\mathit{n}}-2}+\frac{\varphi_1}{2}-\frac{\pi}{4} \\ \end{gathered}\right) \\ \end{gathered} $$ (27) 式中:$ \theta $为相邻底板岩层间的夹角,此时代表对数螺旋线 r与 rn–2之间夹角。
当$ \dfrac{\mathrm{d}h}{\mathrm{d}\theta}=0 $ 时,h 达到最大破坏深度h0,解得
$$ \theta=\varphi\mathit{\mathit{_{\mathrm{\mathit{n}}}}}+\frac{\pi}{4}-\frac{\varphi_1}{2}-\theta_1-\theta_2-\cdots\theta_{\mathrm{\mathit{n}}-2} $$ (28) 将式(25)、式(26)和式(28)的解代入式(27),得出多层结构底板最大滑移破坏深度的计算公式:
$$ h_0=r_{\mathrm{\mathit{n}}-2}\mathrm{e}^{\left(\varphi\mathit{\mathit{_{\mathrm{\mathit{n}}}}}+\tfrac{\pi}{4}-\tfrac{\varphi_1}{2}-\theta_1-\theta_2-\theta_3-...-\theta_{\mathrm{\mathit{n}}-2}\right)\mathrm{tan\ }\varphi_{\mathrm{\mathit{n}}}}\mathrm{cos}\ \varphi_{\mathrm{\mathit{n}}} $$ (29) 开采面至最大破坏深度水平距离OG:
$$ OG=r_{\mathit{\mathrm{\mathit{n}}}-2}\mathrm{e}^{\left(\varphi_{\mathrm{\mathit{n}}}+\tfrac{\pi}{4}-\tfrac{\varphi_1}{2}-\theta_1-\theta_2-\theta_3-...-\theta_{\mathrm{\mathit{n}}-2}\right)\mathrm{tan}\ \varphi_{\mathrm{\mathit{n}}}}\sin\ \varphi_{\mathrm{\mathit{n}}} $$ (30) 1.2.4 工况4
由图1几何关系
$$ \overline{N^{t} A^{\prime}}=x_{0}-\sum_{i=1}^{n-1} 2 H^{{i}} \tan \left(\frac{\pi}{4}-\frac{\varphi_{{i}}}{2}\right) $$ (31) 在△N’A’B中
$$ \overline{OA'}=\frac{\overline{N'A'}\cdot\tan\left(\dfrac{\pi}{4}+\dfrac{\varphi_{\mathit{\mathit{\mathit{\mathrm{\mathit{i}}}}}+1}}{2}\right)}{2} $$ (32) $$ H_0=\overline{O A^{\prime}}+\sum_{i=1}^{n-1} H^{{i}} $$ (33) $ \displaystyle\sum_{i=1}^{n-2}H\mathit{^{\mathrm{\mathit{i}}}} < H_0\leqslant\displaystyle\sum_{i=1}^{n-1}H^{\mathrm{\mathit{i}}} $,即主动极限区深度H0处于岩体第2层至第n–1层之间,主动极限区、被动极限区分别为2条直线,但过渡区滑移线由多段对数螺旋线组成,且各段螺旋线起始半径由主动极限区深度H0所处岩层位置确定。该工况底板层数为n层时,主动极限区深度H0所处岩层位置有i’=n–2种情况,需要对主动极限区深度H0每种情况i进行分析计算,后依次计算出相应对数螺旋线起始半径及起始半径间夹角。
第1段螺旋线起始半径r0的计算公式如下:
由图1几何关系
$$ ON = \frac{{{H^1}}}{{{\mathrm{sin}}\left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right)}} $$ (34) 在△OMN中
$$ MN = {H^1}\left[ {{\mathrm{tan}}\left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2}} \right) + {\mathrm{tan}}\left( {\frac{\pi }{4} - \frac{{{\varphi _1}}}{2}} \right)} \right] $$ (35) $$ M^{\prime} N^{\prime} = M N+ \sum_{l=2}^{i} \dfrac{H^{{{l}}}}{\tan \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{l}}}}{2}\right)} + \sum_{l=2}^{i} \dfrac{H^{{{l}}}}{\tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{{l}}}}{2}\right)} $$ (36) 在△O′M′N′中
$$ \begin{array}{c} {O^\prime }{N^\prime } = {M^\prime }{N^\prime } \cdot \cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i}} + 1}}}}{2}} \right)= \\ \left[ {\begin{array}{*{20}{c}} {{H^\prime }\left[ {\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right) + \tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _1}}}{2}} \right)} \right]}+ \\ {\left. { \displaystyle\sum\limits_{l = 2}^i {\dfrac{{{H^{{l}}}}}{{\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{l}}}}}{2}} \right)}}} + \sum\limits_{l = 2}^i {\dfrac{{{H^{{l}}}}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{l}}}}}{2}} \right)}}} } \right]} \end{array}} \right. \times \\ \cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i}} + 1}}}}{2}} \right) \\ \end{array} $$ (37) $$ \begin{array}{c} {O^\prime }P = {O^\prime }{N^\prime } \cdot \sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i}} + 1}}}}{2}} \right) = \\ \left[ {\begin{array}{*{20}{c}} {{H^\prime }\left[ {\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right) + \tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _1}}}{2}} \right)} \right]} +\\ {\left. { \displaystyle\sum\limits_{l = 2}^i {\dfrac{{{H^{{l}}}}}{{\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{l}}}}}{2}} \right)}}} + \sum\limits_{l = 2}^i {\dfrac{{{H^{{l}}}}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _l}}}{2}} \right)}}} } \right]} \end{array}} \right. \times \\ \dfrac{{\sin \left( {\dfrac{\pi }{2} + {\varphi _{{{i}} + 1}}} \right)}}{2} \\ \end{array} $$ (38) $$ N^{\prime} A^{\prime}=x_{0}-\sum_{l=1}^{{{t}}} 2 H^{{{l}}} \tan \left(\frac{\pi}{4}-\frac{\varphi_{{{l}}}}{2}\right) $$ (39) $$ N^{\prime} B=\dfrac{N^{\prime} A^{\prime}}{2 \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)}=\dfrac{x_{0}-\displaystyle\sum_{i=1}^{t} 2 H^{\prime} \tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{1}}{2}\right)}{2 \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)} $$ (40) $$ {r_0} = O'N' + N'B $$ (41) 将式(37)、式(40)代入式(41),r0为
$$ \begin{gathered} r_{0}=\left[\begin{gathered} H^{\prime}\left[\tan \left(\dfrac{\pi}{4}+\dfrac{\varphi_{1}}{2}\right)+\tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{1}}{2}\right)\right]+ \\ \displaystyle\sum_{l=2}^{i} \dfrac{H^{{{l}}}}{\tan \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{l}}}}{2}\right)}+\displaystyle\sum_{l=2}^{i} \dfrac{H^{{{l}}}}{\tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{{l}}}}{2}\right)} \end{gathered}\right] \times \\ \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)+\dfrac{x_{0}-\displaystyle\sum_{l=1}^{i} 2 H^{{{l}}} \tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{{l}}}}{2}\right)}{2 \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)} \end{gathered} $$ (42) 由图1几何关系知G′E(h′)段最大值为
$$ \begin{gathered} h_0^{\prime}=\frac{O^{\prime} P+H^{{{i}}+1}+H^{({{i}}+1)+1}+\cdots+H^{{{n}}-1}}{\sin \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}+\theta_1+\theta_2+\cdots+\theta_{{{n}}-{{i}}-1}\right)} \times \\ {\mathrm{e}}^{\left(\varphi_{{n}}+\tfrac{\pi}{4}-\tfrac{\varphi_{{{i}}+1}}{2}-\theta_1-\theta_2-\cdots-\theta_{{{n-i}}-1}\right) \tan\; {{n}}} \cos\; \varphi_{{n}} \end{gathered} $$ (43) 式中:θn−i−1为对数螺旋线半径rn−i−2与rn−i−1间夹角。
最终得出$ \displaystyle\sum_{i=1}^{n-2} H^{{{i}}}< H_{0} \leqslant \displaystyle\sum_{i=1}^{n-1} H^{\prime} $情况下的多层结构底板最大滑移破坏深度h0:
$$ \begin{gathered} {h_0} = {h^\prime }_0 + \sum\limits_{i = 1}^{n - 2} {{H^{{i}}}} - {O^\prime }P = \\ \frac{{{O^\prime }P + {H^{{{i}} + 1}} + {H^{({{i}} + 1) + 1}} + \cdots + {H^{{{n}} - 1}}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2} + {\theta _1} + {\theta _2} + \cdots + {\theta _{{{n}} -{{ i}} - 1}}} \right)}} \times \\ {{\mathrm{e}}^{ {\left({\varphi _{{n}}} + \tfrac{\pi }{4} - \tfrac{{{\varphi _{{{i}} + 1}}}}{2} { - {\theta _1} - {\theta _2} - \cdots - {\theta _{{{n - i}} - 1}}} \right)\tan \; {\varphi _{{n}}}} }}\cos\; {\varphi _{{n}}}+ \\[-6pt] \sum\limits_{i = 1}^{n - 2} {{H^{{i}}}} - {O^\prime }P \\ \end{gathered} $$ (44) 在△O′G′E中
$$ \begin{gathered} {O^\prime }{G^\prime } = r\sin \;\alpha = \\ \frac{{{O^\prime }P + {H^{{{i }}+ 1}} + {H^{({{i}} + 1) + 1}} + \cdots + {H^{{{n}} - 1}}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i}} + 1}}}}{2} + {\theta _1} + {\theta _2} + \cdots + {\theta _{{{n - i}} - 1}}} \right)}} \times \\ {{\mathrm{e}}^{ {\left({\varphi _{{n}}} + \tfrac{\pi }{4} - \tfrac{{{\varphi _{{{i }}+ 1}}}}{2} { - {\theta _1} - {\theta _2} - \cdots - {\theta _{{{{{n - i}} }}- 1}}} \right)\tan \;{\varphi _{{n}}}} }}\sin\; {\varphi _{{n}}} \\ \end{gathered} $$ (45) 在△OFN′中
$$ F N^{\prime}=\sum_{1}^{i} H^{\prime} \cdot \tan \left(\frac{\pi}{4}-\frac{\varphi_{{{i}}}}{2}\right) $$ (46) $$ PN' = O'P \cdot {\mathrm{tan}}\left( {\frac{\pi }{4} - \frac{{{\varphi _{{{i}} + 1}}}}{2}} \right) $$ (47) $$ \begin{gathered} P F=P N^{\prime}-F N^{\prime}= \\[-6pt] O^{\prime} P \cdot \tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{{i}}+1}}{2}\right)-\sum_{1}^{i} \tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{{i}}}}{2}\right) \end{gathered}$$ (48) 由图1几何关系有
$$ OG = O'G' + PF $$ (49) 将式(45)、式(48)代入式(49),知开采面至最大破坏深度水平距离OG为
$$ \begin{gathered} OG = \frac{{{O^\prime }P + {H^{{{i}} + 1}} + {H^{({{i}} + 1) + 1}} + \cdots + {H^{{{n }}- 1}}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2} + {\theta _1} + {\theta _2} + \cdots + {\theta _{{{n - i}} - 1}}} \right)}} \times \\ {{\mathrm{e}}^{ {\left({\varphi _{{n}}} + \tfrac{\pi }{4} - \tfrac{{{\varphi _{{{i}} + 1}}}}{2} { - {\theta _1} - {\theta _2} - \cdots - {\theta _{{{n - i}} - 1}}} \right)\tan\; {\varphi _{{n}}}} }}\sin \;{\varphi _{{n}}} + \\ {O^\prime }P \cdot \tan \left( {\frac{\pi }{4} - \frac{{{\varphi _{{{i }}+ 1}}}}{2}} \right) - \sum\limits_1^i {\tan } \left( {\frac{\pi }{4} - \frac{{{\varphi _{{i}}}}}{2}} \right) \\ \end{gathered} $$ (50) 1.2.5 工况5
$ \displaystyle\sum_{i=1}^{n-1} H^{{{i}}}< H_{0} $,即主动极限区深度H0处于第n层岩体内,设过渡区滑移线由1段对数螺旋线组成,主动极限区、被动极限区分别为2条直线。
第1段螺旋线起始半径r0由式(41)可知:
$$ \begin{gathered} {r_0} = \left[ \begin{gathered} {H^\prime }\left[ {\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{\mathrm{1}}}}}{2}} \right) + \tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{\mathrm{1}}}}}{2}} \right)} \right] + \\ \sum\limits_2^{n - 1} {\dfrac{{{{\mathrm{H}}^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} + \sum\limits_2^{n - 1} {\dfrac{{{{\mathrm{H}}^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} \\ \end{gathered} \right] \times \\ \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)+\dfrac{x_0-\displaystyle\sum_1^{n-1} 2 H^{{i}} \tan \left(\dfrac{\pi}{4}-\dfrac{\varphi_{{l}}}{2}\right)}{2 \cos \left(\dfrac{\pi}{4}+\dfrac{\varphi_{{{i}}+1}}{2}\right)} \end{gathered} $$ (51) 由图1几何关系知G′E(h′)段最大值为
$$ \begin{gathered} {h_0}^\prime = r\cos \;\alpha = {r_0}{{\mathrm{e}}^{\theta \tan \;{\varphi _{{n}}}}}\cos \;\left( {\theta + \dfrac{{{\varphi _{{n}}}}}{2} - \dfrac{\pi }{4}} \right) = \\ \left\{ \begin{gathered} \left[ \begin{gathered} {H^1}\left[ {\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right) + \tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _1}}}{2}} \right)} \right] + \\ \displaystyle\sum_{i = 2}^{n - 1} {\dfrac{{{H^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} + \sum\limits_{i = 2}^{n - 1} {\dfrac{{{H^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} \\ \end{gathered} \right] \\ \cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2}} \right) + \dfrac{{{x_0} - \displaystyle\sum\limits_{i = 1}^{n - 1} 2 {H^{{i}}}\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2}} \right)}} \\ \end{gathered} \right\}\times \\ {{\mathrm{e}}^{\left( {\tfrac{\pi }{4} + \tfrac{{{\varphi _{{n}}}}}{2}} \right)\tan\; {\varphi _{{n}}}}}\cos\; {\varphi _{{n}}}_{} \\ \end{gathered} $$ (52) 最终得出$ \displaystyle\sum_{i=1}^{n-1} H^{{{i}}}< H_{0} $情况下的多层结构底板最大破坏深度h0:
$$ \begin{gathered} {h_0} = {h_0}^\prime + \sum\limits_{i = 1}^{n - 1} {{H^{{i}}}} - {O^\prime }P = \\ \left\{ \begin{gathered} \left[ \begin{gathered} {H^\prime }\left[ {\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right) + \tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _1}}}{2}} \right)} \right] + \\ \sum\limits_{i = 2}^{n - 1} {\dfrac{{{H^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} + \sum\limits_{i = 2}^{n - 1} {\dfrac{{{H^{{i}}}}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}} \\ \end{gathered} \right] \times \\ \cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2}} \right) + \dfrac{{{x_0} - \displaystyle\sum\limits_{i = 1}^{n - 1} 2 {H^{{i}}}\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\varphi _{{i}}}}}{2}} \right)}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _{{{i }}+ 1}}}}{2}} \right)}} \\ \end{gathered} \right\}\times \\ {{\mathrm{e}}^{\left( {\tfrac{\pi }{4} + \tfrac{{{\varphi _{{n}}}}}{2}} \right)\tan \;{\varphi _{{n}}}}}\cos\; {\varphi _{{n}}} + \sum\limits_{i = 1}^{n - 1} {{H^{{i}}}} - {O^\prime }P \\ \end{gathered} $$ (53) 开采面至最大破坏深度水平距离OG为
$$ \begin{gathered}OG=\left\{\begin{array}{l}\left[\begin{array}{l}{H}^{1}\left[\mathrm{tan}\left(\dfrac{\pi }{4}+\dfrac{{\varphi }_{1}}{2}\right)+\mathrm{tan}\left(\dfrac{\pi }{4}-\dfrac{{\varphi }_{1}}{2}\right)\right]+\\ {\displaystyle \sum _{i=2}^{n-1}\dfrac{{H}^{{{i}}}}{\mathrm{tan}\left(\dfrac{\pi }{4}+\dfrac{{\varphi }_{{{i}}}}{2}\right)}}+{\displaystyle \sum _{i=2}^{n-1}\dfrac{{H}^{{{i}}}}{\mathrm{tan}\left(\dfrac{\pi }{4}-\dfrac{{\varphi }_{{{i}}}}{2}\right)}}\end{array}\right]\\ \mathrm{cos}\left(\dfrac{\pi }{4}+\dfrac{{\varphi }_{{{i}}+1}}{2}\right)+\dfrac{{x}_{0}-{\displaystyle \sum _{i=1}^{n-1}2}{H}^{{{i}}}\mathrm{tan}\left(\dfrac{\pi }{4}-\dfrac{{\varphi }_{{{i}}}}{2}\right)}{2\mathrm{cos}\left(\dfrac{\pi }{4}+\dfrac{{\varphi }_{{{i}}+1}}{2}\right)}\end{array}\right\}\text{×}\\ {{\mathrm{e}}}^{\left(\tfrac{\pi }{4}+\tfrac{{\varphi }_{{{n}}}}{2}\right)\mathrm{tan}\;{\varphi }_{{{n}}}}\mathrm{sin}\;{\varphi }_{{{n}}}+\\ {O}^{\prime }P\cdot \mathrm{tan}\left(\dfrac{\pi }{4}-\dfrac{{\varphi }_{{{i}}+1}}{2}\right)-{\displaystyle \sum _{1}^{i}\mathrm{tan}}\left(\dfrac{\pi }{4}-\dfrac{{\varphi }_{{{i}}}}{2}\right)\end{gathered} $$ (54) 将各工况推导公式整理后利用主动极限区深度进行逻辑流程判别,后进行底板破坏深度及开采面至底板破坏深度水平距离的计算,如图3所示。
1.3 多层结构底板破坏深度计算误差分析
因煤岩开采面底板根据柱状进行分层时各岩层内摩擦角不同,采用塑性滑移线理论分析时会产生多条对数螺旋曲线半径。对对数螺旋曲线半径的求解采用工况2加以说明,其他工况对数螺旋曲线半径采用类似的求解方式,不再赘述。
n=2,即假设底板层数为2层时,如图4所示,当滑移线主动极限区深度H0小于首层底板岩层厚度H′时,对数螺旋线将穿过上层底板岩层,此时O点将作为2段对数螺旋线的旋转中心点,r1同时作为第1段对数螺旋线的终线和第2段对数螺旋线的起点线。
对数螺旋曲线起始半径r0由图4几何关系知:
$$ {r_0} = \dfrac{{{x_0}}}{{2\cos \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2}} \right)}} $$ (55) 第2段对数螺旋线起始半径$ r_{1} $可由方程组(56)联立求解:
$$ \left\{ \begin{gathered} {r_1} = \dfrac{{{H^1}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1}} \right)}} \\ {r_1} = {r_0}{{\mathrm{e}}^{{\theta _1}\tan\;{\varphi _1}}} \\ \end{gathered} \right. $$ (56) 化简为
$$ {H^1} = {r_0}{{\mathrm{e}}^{{\theta _1}{\mathrm{tan}}\;{\varphi _1}}} \cdot \sin \left( {\frac{\pi }{4} + \frac{{{\varphi _1}}}{2} + {\theta _1}} \right) $$ (57) 对于式(57)中的未知数$ \theta_{1} $和工况2中式(9)的求解方法相同,采用泰勒展开求得近似解$ \theta_{1} $,最终可求得对数螺旋线起始半径$ r_{1} $。
n=3,即假设底板层数为3层时,如图5所示,此时O点将作为3段对数螺旋线的旋转中心点。
对数螺旋半径$ r_{1}, r_{2} $求解如下
$$ \left\{ \begin{gathered} {r_1} = \dfrac{{{H^1}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1}} \right)}} \\ {r_1} = {r_0}{{\mathrm{e}}^{{\theta _1}{\mathrm{tan}}\;{\varphi _1}}} \\ \end{gathered} \right. $$ (58) $$ \left\{ \begin{gathered} {r_2} = \dfrac{{{H^1} + {H^2}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1} + {\theta _2}} \right)}} \\ {r_2} = {r_1}{{\mathrm{e}}^{{\theta _2}{\mathrm{tan}}\;{\varphi _2}}} \\ \end{gathered} \right. $$ (59) 对$ r_{1}, r_{2} $求解方式与n=2时$ r_{1} $的求解方法相同,即通过泰勒展开式后进行一元四次方程求解,且$ r_{2} $求解需经过2次泰勒式转化分别求出$ \theta_{1} $、$ \theta_{2} $。
n=4,即假设底板层数为4层时,如图6所示,此时O点将作为4段对数螺旋线的旋转中心点。
对数螺旋半径$ r_{1}, r_{2}, r_{3} $用方程组(60)求解
$$ \left\{ \begin{gathered} {r_1} = \dfrac{{{H^1}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1}} \right)}} \\ {r_1} = {r_0}{{\mathrm{e}}^{{\theta _1}{\mathrm{tan}}\;{\varphi _1}}} \\ \end{gathered} \right. $$ (60) $$ \left\{ \begin{gathered} {r_2} = \dfrac{{{H^1} + {H^2}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1} + {\theta _2}} \right)}} \\ {r_2} = {r_1}{{\mathrm{e}}^{{\theta _2}{\mathrm{tan}}\;{\varphi _2}}} \\ \end{gathered} \right. $$ (61) $$ \left\{ \begin{gathered} {r_3} = \frac{{{H^1} + {H^2} + {H^3}}}{{\sin \left( {\dfrac{\pi }{4} + \dfrac{{{\varphi _1}}}{2} + {\theta _1} + {\theta _2} + {\theta _3}} \right)}} \\ {r_3} = {r_2}{{\mathrm{e}}^{{\theta _3}{\mathrm{tan}}\;{\varphi _3}}} \\ \end{gathered} \right. $$ (62) 则对数螺旋半径$ r_{3} $求解需经过3次泰勒式转化求出近似解$ \theta_{1} $、$ \theta_{2} $和$ \theta_{3} $,且3个近似解误差在逐次累积。可以预见若底板层数持续增加,则对数螺旋曲线半径的求解需经过更多次泰勒式转化,其产生的累积误差会随之加大,导致其计算的破坏深度误差偏大,且计算量急剧加重,不易进行人工手算。因此考虑采用Matlab内置函数编程语言进行对数螺旋线半径求解,避免泰勒展开式产生的累积误差,并且计算结果精度不受多层结构底板层数的影响。
2. 多层结构底板破坏深度计算系统软件开发
2.1 多层结构底板破坏深度计算系统V1.0开发环境
多层结构底板破坏深度计算系统V1.0 (简称计算系统)主要采用Matlab语言中的math operations模块开发而成。Matlab是用于算法开发、数据可视化、数据分析及数值计算的高级技术计算语言和交互式环境的商业数学软件,运用guide指令搭接的GUI界面主要控件见表1。
表 1 界面主要控件Table 1. Interface main controls英文名称 中文名称 控件作用 Button 按钮控件 在程序中显示按钮 Entry 输入控件 用于显示简单的文本内容 Label 标签控件 可以显示文本 Text 文本控件 用于显示多行文本 Toplevel 容器控件 用来提供一个单独对话框 MessageBox 对话框控件 用于显示应用程序的消息框 2.2 计算误差控制
使用Matlab内置函数计算多层结构底板破坏深度过程中,利用EXCEL计算部分已知参数的底板破坏深度及开采面至底板最大破坏深度水平距离的理论值,佐证Matlab编码的准确性及人为因素影响所计算出的软件值。
if strcmp(f{i}, 'MaxDegree')
v = sym(v);
elseif strcmp(f{i}, 'ReturnConditions')
f{i} = 'OutputType';
if v == true
v = evalin(symengine, '"FullMode"');
else
v = evalin(symengine, '"CompatibleMode"');
end
elseif v == true
v = evalin(symengine, 'TRUE');
else
% v is false
v = evalin(symengine, 'FALSE');
end
entries(i) = feval(symengine, '_equal', sym(f{i}), v);
end
entries(end) = evalin(symengine, 'VectorFormat = TRUE');
entries = feval(symengine, 'op', entries);
T = feval(symengine, 'table', entries);
local function warnIfParams
warn if the solution depends on parameters and conditions
function warnIfParams(parameters, conditions)
If isempty(parameters)
paramstring = char(parameters(1));
for i = 2:numel(parameters)
paramstring = [paramstring ', ' char(parameters(i))]; %#ok
end
以上为多层结构底板破坏深度计算系统软件计算过程中的部分核心代码。
从表2中可以看到对于工况1由于不涉及泰勒展开式及方程组的求解其最终计算结果与Matlab软件计算结果相同。当底板层数分别为3层和7层时,对于工况4中的理论值 (利用EXCEL对多层结构底板破坏深度推导公式进行求解)和软件值结果误差也随之加大,其原因在于求解对数螺旋线半径时多次使用了泰勒展开式,且误差随结构底板岩层层数增加误差也随之累积,符合预期设想。
表 2 多层结构底板破坏深度及开采面至破坏深度水平距离误差对比Table 2. Comparison of damage depth and horizontal distance error from mining surface to damage depth of multi-layered structure floor输入参数 输出结果 工况 h0 /m OG/m 底板层数n 底板超前塑性破坏长度X0 /m 底板层厚Hi/m 内摩擦角φi /(°) 理论值 软件值 理论值 软件值 1 7.9 25 37 16.318 16.2975 12.295 12.281 工况1 2 7.9 5.5 39 16.59 16.5868 12.528 12.5088 工况5 19.5 36 3 9.2 6.8 30 16.74 16.0903 10.494 10.8452 工况4 2.3 35 10.4 36 3 9.54 10 39 19.404 19.3943 14.103 14.0908 工况2 3.6 35 7.3 36 7 22.7 4.3 39 44.351 42.8537 30.562 29.5707 工况4 3 34 5.1 38 3.5 35 6 30 2.4 32 21 33 3 4.4 5 37 7.197 8.04518 5.026 5.0271 工况3 0.6 33 8.7 32 2.3 多层结构底板破坏深度计算系统软件主界面与功能
如图7所示,将相应参数输入到输入框,点击“计算”按钮进行最终计算。如输入底板层数“3”,页面中会显示相应层数,后在列表中依次输入相应的底板层厚和内摩擦角;然后输入底板超前塑性破坏长度x0 (实测值为空,计算取理论值;实测值有实值情况,系统默认实测值)。如输入实测值“5.6”后,点击“计算”按钮,界面会显示该3层底板最大破坏深度和开采面至最大破坏深度的水平距离分别为10.557、7.012 m。
若现场因客观原因无法测得底板超前塑性破坏长度时,需计算底板超前塑形破坏长度的理论值,则可以分别输入H1(煤层的埋深,m);H2(煤层的采高,m);K(支撑压力峰值系数);θ(煤层的变形角,(°));φc(煤体的内摩擦角,(°));φc*(煤体的残余内摩擦角,(°));C1(煤层与顶底板接触面上的黏聚力,MPa);C2(煤体的黏聚力,MPa);C*(煤体的残余黏聚力,MPa);f1(煤层与顶板界面间的摩擦力,N);r(煤层覆盖的平均容重,N/m3)。如输入图8页面中的地质参数,可计算底板超前塑形破坏长度的理论值为18.697 m。
3. 应用分析
3.1 平煤矿区软件计算结果分析
根据平煤某矿工作面底板岩体地质情况,工作面超前塑性区长度现场实测值为7.9 m。虽然下保护层层位与寒武系灰岩含水层之间有平均64 m厚的隔水层,在正常情况下隔水层可以阻隔承压水的侵入,但受开采扰动和寒武系高水头压力的影响,底板将形成扰动破坏带和承压水导升带,使得有效隔水层厚度减小,若遇大的断裂构造,承压水可通过隔水层较薄地带或构造破碎带进入工作面。参考周边矿井煤岩层工作面底板扰动破坏深度,取下保护层开采工作面底板岩层厚度H=25 m内的岩体底板进行最大破坏深度理论计算。
在《深部煤岩层复合结构底板破坏机制及应用研究》[21]中李昂等针对此工程把相近岩性或岩层厚度较薄的岩层划分为1层,摩擦角取其平均值,岩层较厚且岩性差异较大时划分为另外层,分别进行单层、双层及三层底板最大破坏深度计算,并通过现场监测数据判断底板最大破坏深度为17.8 m。
本章节将通过多层结构底板破坏深度计算系统V1.0软件对该理论值、拟合值与现场实测值进行对比分析,表3中可以看到1~3层结构底板软件计算值与理论值结果误差极小。
表 3 数值计算结果对比Table 3. Comparison of numerical calculation results底板层数n 底板超前塑性破坏长度X0/m 底板层厚Hi/m 内摩擦角φi/(°) 底板最大破坏深度h0/m 理论值 软件值 拟合值 实测值 1 7.9 25 37 16.3 16.298 18.12 17.9 2 7.9 5.5 39 16.59 16.587 19.5 36 3 7.9 5.5 39 17.08 17.077 7.2 34 12.3 38 5 7.9 4.6 36 17.33 17.761 7.17 39 6.58 35 6.19 34 2.29 39 因以往研究中以1~3层作为结构底板破坏深度分析计算,不足以应对复杂的底板岩层。因此将该开采工作面底板岩层厚度H=25 m内的岩体进行5层分析计算,依据柱状图(图9)把L8灰岩及厚度较小的煤线和砂质泥岩作为首层,厚度和为4.6 m,内摩擦角取该层位岩性平均内摩擦角36°;第2层以L7灰岩为主,层厚为7.17 m,内摩擦角39°;第3层由厚度较小的砂质泥岩、煤和泥岩构成,总厚度为6.58 m,平均内摩擦角为35°;第4层以厚度为6.19的黑灰色泥岩为主,内摩擦角34°;第5层是厚度为2.29 m的石灰岩,内摩擦角39°。将以上各参数输入到该计算系统,得到底板最大破坏深度软件值为17.761 m。从表3中可看到其软件值更接近现场实测值和拟合值,计算精度更高。虽然拟合值[23]和实测值相近,但拟合值往往取自不同矿井数据所得到的经验公式,当具体到实际工程时,可能得到较大的偏差。
3.2 渭北矿区软件计算结果分析
李昂在《带压开采下底板渗流与应力耦合破坏突水机理及其工程应用》[12]中将董家河煤矿5号煤层22507工作面作为工程研究背景,相关力学岩性指标取值为:x0=5.7 m、H=350 m、φd=35°。计算出单层底板破坏深度为10.9 m,开采面至底板破坏深度水平距离为7.6 m,现场实测底板破坏深度结果为10.8 m。
将该地质相关参数输入到计算系统软件中,如图10所示,可以看到底板破坏深度与开采面至底板破坏深度水平距离计算结果与该文中理论值相同,该软件的计算精度可达到98%。
3.3 不同矿区软件计算结果综合分析
不同钻孔揭露的岩性存在差异性,岩层厚度、层数往往存在较大差异,导致底板采动破坏不同。开展现场测试过程中,往往钻取岩芯或者借助周边较近钻孔柱状。软件分析结果与现场实测数据进行对比分析,通过软件分析所得到的结果与实测值一致。不可否认,如果换作其他矿区进行现场实测,在岩性差别较大时计算结果也会有所不同。在许多矿井没有实测值或很难开展实测的情况下,该软件依据柱状和基本力学参数便可以计算底板破坏深度的取值,且可以代替繁琐的手工计算过程,优势是明显的。
4. 结 论
1)在3层复合结构底板塑性滑移线场力学模型基础上,推导出n层结构底板破坏深度解析解。利用手算进行底板最大破坏深度计算时部分工况计算过程需多次利用泰勒展开式转化为一元四次方程式。随着层数的增加,会逐次累积误差,且计算值只能取近似解。以主动极限区深度H0所处底板位置为判别依据,对多层结构底板计算工况进行了优化,提出5种计算工况,并形成了逻辑判别流程图。
2)首次采用了Matlab内置函数运用编程语言进行多层结构底板破坏深度计算,其自带内置函数可避免手算过程中因泰勒式展开式多次利用而产生的累积误差。借助逻辑判别图开发了多层结构底板破坏深度计算系统V1.0软件,该软件通过输入相应参数可快速地计算出底板破坏深度及开采面至底板破坏深度的水平距离。
3)通过将平煤某矿及董家河煤矿的岩性力学指标参数值代入该计算系统,得到该软件计算结果与理论值相一致。并将平煤某矿底板岩层进行五层计算分析,得到底板最大破坏深度理论值与实测结果偏差0.57 m,软件值与实测值相差0.14 m,即根据岩性层数划分越细致,计算精度更高。可以应对更加复杂的底板岩性。
4)该软件具有普适性,所采用的理论是以不同实际柱状和力学参数为依据进行的分层计算,在矿井没有实测值或很难开展实测的情况下,仅仅根据柱状和基本力学参数便可以快捷准确地计算底板破坏深度的取值,且可以代替繁琐的手工计算过程,相比常规单一岩层塑性滑移线理论准确性更高。对研究我国华北型煤田煤岩层底板破坏规律具有重要的工程价值。
-
表 1 路段1期间各算法均方根误差及实时性
Table 1 Root mean square error and real-time performance of each algorithm during section 1
荷载质量−算法 RMSE 识别时间/s FFRLS 5.197 1.150 EKF 3.657 1.180 DFFRLS 3.065 1.100 坡度坡度−算法 RMSE 识别时间/s FFRLS 0.131 1.490 RLS-EKF 0.096 1.460 EKF 0.281 1.500 DFFRLS-AEKF 0.053 1.430 表 2 路段2期间各算法均方根误差及实时性
Table 2 Root mean square Error and real-time performance of each algorithm during section 2
荷载质量−算法 RMSE 识别时间/s FFRLS 5.197 1.370 EKF 4.602 1.360 DFFRLS 4.295 1.310 坡度坡度−算法 RMSE 识别时间/s FFRLS 0.149 1.600 RLS-EKF 0.057 1.550 EKF 0.058 1.610 DFFRLS-AEKF 0.031 1.520 表 3 测试车辆主要技术参数
Table 3 Parameters related to test vehicle parameters
参数 数值 参数 数值 适应坡度/(°) $ \leqslant $15 额定荷载/t 20 车型长度/m 18.1 最大转速/(r·min−1) 2 600 车型宽度/m 1.04 车型高度/m 1.45 转动惯量/(kg.m2) 2 行走轮
直径/m0.175 最大
牵引力/kN80 制动力
范围/kN120~160 空载最
大速度/(m·s−1)2.0 重载最
大速度/(m·s−1)0.5 牵引电
机功率/kW60 驱动轮
直径/mm340 表 4 实车测试各算法均方根误差及实时性
Table 4 Real vehicle tests the root mean square error and real-time performance of each algorithm
荷载质量−算法 RMSE 识别时间/s FFRLS 151.804 1.370 EKF 187.989 1.360 DFFRLS 123.234 1.320 坡度坡度−算法 RMSE 识别时间/s FFRLS 0.468 1.660 RLS-EKF 0.388 1.630 EKF 0.427 1.650 DFFRLS-AEKF 0.328 1.620 -
[1] 葛世荣,胡而已,裴文良. 煤矿机器人体系及关键技术[J]. 煤炭学报,2020,45(1):455−463. GE Shirong,HU Eryi,PEI Wenliang. System and key technology of coal mine robot[J]. Journal of China Coal Society,2020,45(1):455−463.
[2] 贺海涛, 廖志伟, 郭 卫. 煤矿井下无轨胶轮车无人驾驶技术研究与探索[J]. 煤炭科学技术, 2022, 50 (S1): 212−217. HE Haitao, LIAO Zhiwei, GUO Wei. Research and exploration on driverless technology of trackless rubber wheel vehicle in underground coal mine[J]. Coal Science and Technology, 2022, 50(S1) : 212−217.
[3] 王国法. 煤矿智能化最新技术进展与问题探讨[J]. 煤炭科学技术,2022,50(1):1−27. WANG Guofa. New technological progress of coal mine intelligence and its problems[J]. Coal Science and Technology,2022,50(1):1−27.
[4] 陈杨阳, 霍振龙, 刘智伟, 等. 我国煤矿运输机器人发展趋势及关键技术[J]. 煤炭科学技术, 2020, 48(7): 233−242. CHEN Yangyang, HUO Zhenlong, LIU Zhiwei, et al. Development trend and key technology of coal mine transportation robot in China[J]. Coal Science and Technology, 2020, 48 (7): 233−242.
[5] 鲍久圣,邹学耀,陈 超,等. 重型煤炭运输车分布式混合动力系统设计及控制策略[J]. 煤炭学报,2021,46(2):667−676. BAO Jiusheng,ZOU Xueyao,CHEN Chao,et al. Design and control strategy of distributed hybrid drive system for heavy coal trucks[J]. Journal of China Coal Society,2021,46(2):667−676.
[6] 侯 刚, 王国法, 薛忠新, 等. 煤矿辅助运输自动驾驶关键技术与装备[J]. 采矿与岩层控制工程学报, 2022, 4(3): 5−17. HOU Gang, WANG Guofa, XUE Zhongxin, et al. Key technologies and equipment for automatic driving of coal mine auxiliary transportation[J]. Journal of Mining and Strata Control Engineering, 202, 4(3): 5−17.
[7] JO K,LEE M,SUNWOO M. Road slope aided vehicle position estimation system based on sensor fusion of GPS and automotive onboard sensors[J]. IEEE Transactions on Intelligent Transportation Systems,2015,17(1):250−263.
[8] BROWN A, BRENNAN S. Model-based vehicle state estimation using previewed road geometry and noisy sensors[A]. Dynamic Systems and Control Conference [C]. New York: American Society of Mechanical Engineers, 2012: 591−600.
[9] 雷雨龙,付 尧,刘 科,等. 基于扩展卡尔曼滤波的车辆质量与道路坡度识别[J]. 农业机械学报,2014,45(11):8−13. LEI Yulong,FU Yao,LIU Ke,et al. Vehicle mass and road slope estimation based on extended Kalman filtering[J]. Transactions of the Chinese Society for Agricultural Machinery,2014,45(11):8−13.
[10] YU Z, FENG Y, XIONG L, et al. Vehicle mass estimation for four in-wheel-motor drive vehicle[A]. Electrical Engineering and Control: Selected Papers from the 2011 International Conference on Electric and Electronics[C]. Berlin: Springer Verlag, 2011: 117−125.
[11] ZHANG Y,ZHANG Y,AI Z,et al. A cross iteration estimator with base vector for estimation of electric mining haul truck's mass and road grade[J]. IEEE Transactions on Industrial Informatics,2018,14(9):4138−4148. doi: 10.1109/TII.2018.2794513
[12] 张荣立, 何国纬. 采矿技术手册[M]. 北京: 煤炭工业出版社, 2003. ZHANG Rongli, HE Guowei. Mining technical manual (1st edition)[M]. Beijing: China Coal Industry Press, 2003.
[13] LI X,MA J,ZHAO X,et al. Intelligent two-step estimation approach for vehicle mass and road grade[J]. IEEE Access,2020,8:218853−218862. doi: 10.1109/ACCESS.2020.3042656
[14] 张杰烁,刘 明,李 鑫,等. 基于递归最小二乘法的回声状态网络算法用于心电信号降噪[J]. 生物医学工程学杂志,2018,35(4):539−549. ZHANG Jieshuo,LIU Ming,LI Xin,et al. An echo state network algorithm based on recursive least square for electrocardiogram denoising[J]. Journal of Biomedical Engineering,2018,35(4):539−549.
[15] 刘 鹏,李云伍,梁新成. 基于遗忘递推最小二乘与自适应无迹卡尔曼滤波的锂电池SOC识别[J]. 汽车技术,2022,557(2):21−27. LIU Peng,LI Yunwu,LIANG Xincheng. Estimation of lithium battery SOC based on FFRLS and AUKF[J]. Automotive Technology,2022,557(2):21−27.
[16] 王 浩,郑燕萍,虞 杨. 基于动态优选遗忘因子最小二乘在线识别的磷酸铁锂电池SOC估算[J]. 汽车技术,2021,553(10):23−29. WANG Hao,ZHENG Yanping,YU Yang. Lithium iron phosphate battery SOC estimation based on the least square online identification of dynamic optimal forgetting factor[J]. Automotive Technology,2021,553(10):23−29.
[17] 孙金磊,邹 鑫,顾浩天,等. 基于FFRLS-EKF联合算法的锂离子电池荷电状态估计方法[J]. 汽车工程,2022,44(4):505−513. SUN Jinlei,ZOU Xin,GU Haotian,et al. State of charge estimation for lithium-ion battery based on FFRLS-EKF joint algorithm[J]. Automotive Engineering,2022,44(4):505−513.
[18] 雷克兵,陈自强. 基于改进多新息扩展卡尔曼滤波的电池SOC估计[J]. 浙江大学学报(工学版),2021,55(10):1978−1985, 2001. LEI Kebing,CHEN Ziqiang. Estimation of state of charge of battery based on improved multi-innovation extended Kalman filter[J]. Journal of Zhejiang University (Engineering Science),2021,55(10):1978−1985, 2001.
[19] 任志英,沈亮量,黄 伟,等. 基于AEKF的车辆质量与道路坡度实时估计[J]. 振动. 测试与诊断,2020,40(4):758−764, 827. REN Zhiying,SHEN Liangliang,HUANG Wei,et al. Real time estimation of vehicle quality and road slope based on adaptive extended kalman filter[J]. Journal of Vibration, Measurement & Diagnosis,2020,40(4):758−764, 827.
[20] 易国晶. 单轨吊机车在盘江矿井中的应用[J]. 煤炭工程,2020,52(7):88−92. YI Guojing. Application of monorail crane in panjiang mine[J]. Coal Engineering,2020,52(7):88−92.
[21] BOADA B,GARCIA-POZUELO D,BOADA M,et al. A constrained dual Kalman filter based on pdf truncation for estimation of vehicle parameters and road bank angle: analysis and experimental validation[J]. IEEE Transactions on Intelligent Transpo rtation Systems,2016,18(4):1006−1016.
[22] ZHANG Y,ZHANG Y,AI Z,et al. Estimation of electric mining haul trucks' mass and road slope using dual level reinforcement estimator[J]. IEEE Transactions on Vehicular Technology,2019,68(11):10627−10638. doi: 10.1109/TVT.2019.2943574