基于MATLAB定量分析采场复杂力学问题的方法研究

郭勤强 王 哲 陈晓利

(河南省地矿局第一地质矿产调查院)



基于MATLAB定量分析采场复杂力学问题的方法研究

郭勤强 王 哲 陈晓利

(河南省地矿局第一地质矿产调查院)

以煤矿地下开采实测数据或数值模拟输出数据为基础，以采场顶板“梁模型假说”为背景，利用MATLAB软件在曲线拟合、多项式插值、微积分和求方程解等方面的强大功能，将实测或数值计算得出的采场顶梁载荷的离散数据拟合成曲线函数，并通过积分找出梁上剪力和弯矩的分布曲线，求出剪力值和弯矩值,同时通过求根的方法解出剪力和弯矩最大值，找出剪力和弯矩最大值的位置，实现采场顶梁受非线性分布载荷作用下的复杂力学问题的定量分析，为预测顶梁破坏提供了更为精确的方法。

MATLAB 曲线拟合 定量分析 采场顶板 复杂应力

采场周围应力分布极其复杂，由于地质条件的不确定性、岩石本身力学性质的复杂性以及采场自身不断移动、变化的特点，岩梁(板)所受载荷一般并不是均匀分布的，使得开挖空间周围岩体力学问题的定量分析异常困难。如：求矩形硐室周围应力的分布[1-3]，确定采场煤壁前方内外应力场分布范围[4-5]；计算采准巷道底板的位移量[6-7]等，在非均布载荷作用时很难计算出与工程实际相吻合的数值解，而工程中常用的做法是假定载荷均匀分布或只存在线性变化，与工程实际相比误差较大。

本研究提出在采场岩梁实测载荷和数值分析所得载荷的基础上，利用MATLAB软件将岩梁(板)上实测或数值计算得出的分布载荷离散数据拟合成曲线函数，通过积分、微分和求方程根的办法求出复杂受力问题的数值解，从而实现非线性分布载荷作用下的采场矿压定量分析。

1 基于MATLAB定量分析复杂力学问题的方法

“梁的假说”是近代矿压假说的显著成果，由于对属性梁的认识不同，有“双支梁假说”、“悬臂梁假说”、“砌体梁假说”和“传递岩梁假说”。本文以老顶岩梁受非线性分布载荷为例，介绍MATLAB实现复杂力学问题求解的方法：首先将实测的老顶分布载荷或通过数值计算得到的载荷数据，在MATLAB中进行曲线拟合，将其转化成可微积分的函数形式，然后再进行顶梁剪力、弯矩和老顶极限跨度的求解[8]。

1.1 梁上分布载荷的曲线拟合

曲线拟合是用连续曲线近似地比拟平面上离散点组所表示的坐标之间的函数关系，其思想是使拟合曲线能反映这些离散数据的变化趋势，使数据点的误差平方和最小。曲线拟合步骤：①把所给的数据画在一张图中，通过图表判断其数学形式；②确定数学形式中的待定参数；③求得数学模型后，将实际测定的数据与所用公式求出的理论值进行比较，判定其优度。

利用MATLAB软件可以直接调用regress命令:[b,bint,r,rint,stats]=regress(y,x)实现线性拟合，其中b是回归方程中的参数估计值，bint是b的置信区间；r和rint分别表示残差及残差对应的置信区间；stats包含3个数字，分别是相关系数、F统计量的值及对应的概率值[9]。

在MATLAB中也可以通过编写M文件实现线性回归分析。如参考文献[10]中编写的M文件,linregr函数可方便地实现线性回归分析：

(1)

1.2 梁剪力和弯矩的求解

以简支梁为例，运用材料力学[11]理论分析,求解梁上剪力和弯矩(图1)。

梁处于平衡状态时的主矢和对O点的主矩应为零。取梁的dx微元，由于dx很小，可视梁上载荷均匀分布，根据平衡方程ΣF=0，ΣMO=0，得：

(2)

图1 简支梁受力分析

(3)

由公式可看出，载荷q(x)是剪力Q(x)的一阶导数，是弯矩M(x)的二阶导数。在MATLAB中可用int和trapz函数求解曲线积分。其中int函数不仅能求解定积分，也可以求解不定积分,而trapz函数只能求解定积分。

int函数实现积分求解的语法是

(4)

式中，符号表达式为s关于变量x的不定积分；s为符号表达式；x为s的变量。

(5)

式中，符号表达式为s关于变量x的定积分；a、b分别为积分的上、下限。

trapz函数实现积分求解的语法是

(6)

式中，x、y分别为自变量和因变量。

在MATLAB中，还有一个函数cumtrapz可求解表达式或者函数的累计积分，实现积分求解的语法是：>>z=cumtrapz(x,y) .

(7)

利用int函数，可以快速求得积分上下限a、b上的剪力值和弯矩值。

1.3 梁剪力、弯矩最大值和最大值位置的求解

求某一函数最值，数学上常用的方法是：对该函数求导，求出函数的极大值和极小值，与函数的边界值进行比较，即可得到函数最值。本文利用函数求导的办法求解梁剪力和弯矩最大值，并找出最大值的位置。

在MATLAB中提供了求解多项式根的函数roots。r=roots(c)，c为一维向量，返回值为指定多项式的所有根(包括复根)。调用格式：

>>syms x;

>>y=x^3+4*x^2;

>>c=sym2poly(y);

>>r=roots(c)

但函数roots只能求解多项式的根，对于复杂的方程只能通过编写M文件的办法求解。如参考文献[11]利用二分法求根的方法编写的M文件,bisect函数可较为精确地求解复杂方程的所有根。

利用MATLAB中的求根函数roots和bisect求解q(x)的所有根，然后将所得根代入剪力Q(x)，求得剪力Q(x)所有的极大值，与边界值进行比较，即可得出剪力最大值和最大值的位置。同理，可求出弯矩最大值及其位置。

2 工程应用分析

2.1 工程背景及数值建模

焦煤集团赵固一矿井田走向长2.0～5.5 km，倾斜宽9.5～11.0 km，井田面积约43.77 km2。主采煤层为山西组二1煤，属近水平发育的稳定型厚煤层，煤层倾角2°～6°，二1煤平均厚5.9 m。矿井采用立井单水平盘区开拓，井口标高83.8 m，水平标高-525 m，井深608.8 m。二1煤直接顶以泥岩为主，厚度1～6.5 m，平均厚度3.5 m；老顶厚度多为0.94～19.85 m,平均厚度13.5 m的粗、中、细粒砂岩；直接底为深灰色泥岩，平均厚3.4 m；老底为砂岩或砂质泥岩，平均厚14.6 m。

利用有限差分软件FLAC3D模拟赵固一矿井下开挖，记录老顶砂岩层(厚度7.1 m)的垂直应力。考虑边界效应，设计的三维计算模型长×宽×高为200 m×200 m×72 m，共划分77 400个六面体单元，生成网

格节点84 084个。X轴平行于工作面推进方向，采用Mohr-Coulomb本构模型。模型水平方向上左右两侧限制水平位移，模型底面限制水平移动和垂直移动，模型的上部施加上覆岩层的自重应力。

2.2 离散数据的处理和求解

将FLAC3D中记录的采场老顶砂岩层载荷数据导入MATLAB软件中，画出离散数据的散点图[12]，观察曲线形式(图2)。由于回采工作面沿倾斜方向的长度远大于老顶沿走向悬露跨距，可将老顶视为一端由工作面煤壁、另一端由边界煤柱支撑的固定梁。当工作面推进20 m时，判断老顶砂岩是否发生破坏。

图2 原始数据散点

由图2可以看出，x距切眼从-2～6 m、15～22 m的曲线形式，与x从6～15 m的曲线形式明显不同。因此，应采用分段拟合曲线的办法，即前9个点、后7个点、中间9个点分别进行曲线拟合，比较拟合曲线的优度。

采用MATLAB中的polyfit函数对各段进行拟合，求得老顶砂岩层载荷q(x)的分段拟合函数为

(8)

绘制出拟合曲线图(图3)，可看出拟合曲线与实际吻合。

图3 分段拟合曲线

用MATLAB软件中的trapz函数分别对q(x)在各区间上积分，求出各区间内的积分值。老顶砂岩为两端固支状态：

(9)

求出支座反力Q=57.544 kN，而砂岩的抗剪强度为37.64 MPa，意味着此时老顶砂岩在支承压力作用下已发生破断，砂岩层的受力状态由两端固支向简支转移，岩层中最大的弯曲拉应力达到其抗拉强度是岩层由弯曲沉降发展至破坏的力学条件[4]。用int函数可求出弯矩在各区间内的不定积分，然后根据边界条件求出简支梁条件下老顶砂岩的弯矩曲线函数。

图4 剪力图

由图4可看出，剪力为零的点出现在[6,15]区间内，即弯矩最大值出现在[6,15]区间内，利用多项式求根函数roots，可求解剪力Q(x)为0的点即为弯矩最大值点。代入polyval函数，求出弯矩的最大值Mmax。程序如下：

>>p=[0.0120378789,-0.72674351, 16.939958334,-179.45333767,6.061622671];

>>r=roots(p);

r=8.67564

利用MATLAB多项式求值函数polyval，求出Mmax=-156 436.587 5 kN·m，计算出此时砂岩层最大弯曲拉应力为4.76 MPa<[σt]=5.61 MPa，即此时砂岩层不会发生弯曲破断。

当工作面推进21 m时，砂岩层弯曲拉应力为5.12 MPa<[σt]，砂岩层不会发生弯曲破断。

当工作面推进22 m时，砂岩层弯曲拉应力为5.93 MPa>[σt]，砂岩层发生弯曲破断。可以看出砂岩层应在工作面推进[21,22]区间内发生弯曲破断，由于工作面推进距离较大，假定老顶砂岩层发生弯曲破断时的载荷曲线与工作面推进21 m时的载荷曲线一致，解出此时砂岩层发生破断的推进距离为L21=21.36 m；同理，假定老顶砂岩层发生弯曲破断时的载荷曲线与工作面推进22 m时的载荷曲线一致，解出此时砂岩层发生破断的推进距离为L22=21.78 m；二者求均值，得工作面老顶砂岩层弯曲破断的推进距离即老顶初次来压步距：L=21.57 m，距切眼r=9.876 m处发生破断。

2.3 对比线性载荷作用下梁的剪力和弯矩

用regress函数拟合出线性载荷q(x)曲线(图5)。

图5 线性分布的载荷曲线

用上述方法求出在线性载荷q(x)作用下工作面推进20 m时,赵固一矿老顶砂岩层最大剪力值Q线max=59.687 kN,最大弯矩值M线max=164.598 kN·m,最大弯矩值位置r线max=9.687 m。与实际载荷q(x)作用下工作面推进20 m时的最大剪力值Qmax、最大弯矩值Mmax及最大弯矩值的位置rmax相比，误差为

线性载荷作用下的老顶砂岩层初次来压步距L线=19.26 m，而非线性载荷作用下的老顶砂岩层初次来压步距L非线=21.37 m，误差为9.87%。

3 结 论

(1)用MATLAB绘出基础数据点集拟合的函数曲线，观测其变化趋势，同样，将求出的应力、弯矩等曲线绘出，直观具体地表现地下工程复杂的力学特征。

(2)以梁模型为例，将离散的非线性数据点集拟合成可微积分的曲线函数，通过MATLAB中的积分运算求解最大剪力和弯矩值，进而可以对梁是以何种方式发生破坏(剪切破坏还是拉伸破坏)做出明确判断，还可以通过求方程根的方法找出最大剪力值和最大弯矩值的位置，判断梁在何处发生破坏。

(3)数据处理和计算都在MATLAB软件中完成，不仅能够保证计算工程的精度，而且能实现快速求复杂力学问题数值解的目的。

[1] 赵 凯,刘长武,张国良.用弹性力学的复变函数法求解矩形硐室周边应力[J].采矿与安全工程学报,2007,24(3):361-365.

[2] 袁 林，高召宁，孟祥瑞.基于复变函数法的矩形巷道应力集中系数黏弹性分析[J].煤矿安全，2013,44(2):196-200.

[3] 王振武，牛铮铮，冯秀苓.地下矩形硐室应力分布的复变函数解[J].北华航天工业学院学报,2010,20(4):1-6.

[4] 宋振骐.实用矿山压力控制[M].徐州：中国矿业大学出版社，1988.

[5] 宋振骐，刘义学，陈孟伯，等.岩梁断裂前后的支承压力显现及应用的探讨[J].山东矿业学院学报，1984(1):27-39.

[6] 王卫军,黄成光,侯朝炯.综放沿空掘巷底鼓的受力变形分析[J].煤炭学报,2002,27(1):26-30.

[7] 孟祥瑞,徐铖辉,高召宁，等.采场底板应力分布及破坏机理[J].煤炭学报,2010,35(11):1832-1836.

[8] 乔立山，王玉兰，曾锦光.实验数据处理中曲线拟合方法探讨[J].成都理工大学学报：自然科学版,2004,31(1)：91-95.

[9] 魏国祥.基于MATLAB的桥梁健康监测数据处理与可靠度分析[D].兰州：兰州交通大学，2012.

[10] Steven C，Chapra.工程与科学数值方法的MATLAB实现[M].唐玲艳，田尊华译.2版.北京：清华大学出版社，2009.

[11] 刘鸿文.材料力学[M].5版.北京：高等教育出版社，2011.

Study on the Analysis the Quantitatively Method of the Stope Complex Mechanical Problem Based on MATLAB

Guo Qinqiang Wang Zhe Chen Xiaoli

(No.1 Institute of Geological & Mineral Resources Survey of Henan Province)

Based on the measured data or numerical simulate results of underground mining of coal mine,taking the "beam model hypothesis" of stope roof as the research background,using the functions of curve fitting,polynomial interpolation,differential and integral calculus and solve equation of MATLAB software,the curve function of the discrete data of stope roof beam load obtained by actual measurement or numerical calculation are obtained,the distributive curve of shear and bending moment on the beam is obtained by integral operation,the shear value and bending moment value are acquired,besides that,the maximum value of shear and bending moment are calculated by using the roof value calculation method,and the location of the maximum value of shear and bending moment are also determined.The goal of the quantitative analysis of complex mechanical problems of stope roof under the nonlinear distributed loadings is achieved,which provide a more accurate method to predict roof beam damage.

MATLAB,Curve fitting,Quantitative analysis,Stope roof,Complex stress

2016-05-10)

郭勤强(1990—)，男，硕士，471000 河南省洛阳市洛龙区龙门大道573号。