郭勤強(qiáng) 王 哲 陳曉利

(河南省地礦局第一地質(zhì)礦產(chǎn)調(diào)查院)

?

基于MATLAB定量分析采場復(fù)雜力學(xué)問題的方法研究

郭勤強(qiáng) 王 哲 陳曉利

(河南省地礦局第一地質(zhì)礦產(chǎn)調(diào)查院)

以煤礦地下開采實(shí)測數(shù)據(jù)或數(shù)值模擬輸出數(shù)據(jù)為基礎(chǔ)，以采場頂板“梁模型假說”為背景，利用MATLAB軟件在曲線擬合、多項(xiàng)式插值、微積分和求方程解等方面的強(qiáng)大功能，將實(shí)測或數(shù)值計(jì)算得出的采場頂梁載荷的離散數(shù)據(jù)擬合成曲線函數(shù)，并通過積分找出梁上剪力和彎矩的分布曲線，求出剪力值和彎矩值,同時(shí)通過求根的方法解出剪力和彎矩最大值，找出剪力和彎矩最大值的位置，實(shí)現(xiàn)采場頂梁受非線性分布載荷作用下的復(fù)雜力學(xué)問題的定量分析，為預(yù)測頂梁破壞提供了更為精確的方法。

MATLAB 曲線擬合 定量分析 采場頂板 復(fù)雜應(yīng)力

采場周圍應(yīng)力分布極其復(fù)雜，由于地質(zhì)條件的不確定性、巖石本身力學(xué)性質(zhì)的復(fù)雜性以及采場自身不斷移動(dòng)、變化的特點(diǎn)，巖梁(板)所受載荷一般并不是均勻分布的，使得開挖空間周圍巖體力學(xué)問題的定量分析異常困難。如：求矩形硐室周圍應(yīng)力的分布[1-3]，確定采場煤壁前方內(nèi)外應(yīng)力場分布范圍[4-5]；計(jì)算采準(zhǔn)巷道底板的位移量[6-7]等，在非均布載荷作用時(shí)很難計(jì)算出與工程實(shí)際相吻合的數(shù)值解，而工程中常用的做法是假定載荷均勻分布或只存在線性變化，與工程實(shí)際相比誤差較大。

本研究提出在采場巖梁實(shí)測載荷和數(shù)值分析所得載荷的基礎(chǔ)上，利用MATLAB軟件將巖梁(板)上實(shí)測或數(shù)值計(jì)算得出的分布載荷離散數(shù)據(jù)擬合成曲線函數(shù)，通過積分、微分和求方程根的辦法求出復(fù)雜受力問題的數(shù)值解，從而實(shí)現(xiàn)非線性分布載荷作用下的采場礦壓定量分析。

1 基于MATLAB定量分析復(fù)雜力學(xué)問題的方法

“梁的假說”是近代礦壓假說的顯著成果，由于對屬性梁的認(rèn)識(shí)不同，有“雙支梁假說”、“懸臂梁假說”、“砌體梁假說”和“傳遞巖梁假說”。本文以老頂巖梁受非線性分布載荷為例，介紹MATLAB實(shí)現(xiàn)復(fù)雜力學(xué)問題求解的方法：首先將實(shí)測的老頂分布載荷或通過數(shù)值計(jì)算得到的載荷數(shù)據(jù)，在MATLAB中進(jìn)行曲線擬合，將其轉(zhuǎn)化成可微積分的函數(shù)形式，然后再進(jìn)行頂梁剪力、彎矩和老頂極限跨度的求解[8]。

1.1 梁上分布載荷的曲線擬合

曲線擬合是用連續(xù)曲線近似地比擬平面上離散點(diǎn)組所表示的坐標(biāo)之間的函數(shù)關(guān)系，其思想是使擬合曲線能反映這些離散數(shù)據(jù)的變化趨勢，使數(shù)據(jù)點(diǎn)的誤差平方和最小。曲線擬合步驟：①把所給的數(shù)據(jù)畫在一張圖中，通過圖表判斷其數(shù)學(xué)形式；②確定數(shù)學(xué)形式中的待定參數(shù)；③求得數(shù)學(xué)模型后，將實(shí)際測定的數(shù)據(jù)與所用公式求出的理論值進(jìn)行比較，判定其優(yōu)度。

利用MATLAB軟件可以直接調(diào)用regress命令:[b,bint,r,rint,stats]=regress(y,x)實(shí)現(xiàn)線性擬合，其中b是回歸方程中的參數(shù)估計(jì)值，bint是b的置信區(qū)間；r和rint分別表示殘差及殘差對應(yīng)的置信區(qū)間；stats包含3個(gè)數(shù)字，分別是相關(guān)系數(shù)、F統(tǒng)計(jì)量的值及對應(yīng)的概率值[9]。

在MATLAB中也可以通過編寫M文件實(shí)現(xiàn)線性回歸分析。如參考文獻(xiàn)[10]中編寫的M文件,linregr函數(shù)可方便地實(shí)現(xiàn)線性回歸分析：

(1)

1.2 梁剪力和彎矩的求解

以簡支梁為例，運(yùn)用材料力學(xué)[11]理論分析,求解梁上剪力和彎矩(圖1)。

梁處于平衡狀態(tài)時(shí)的主矢和對O點(diǎn)的主矩應(yīng)為零。取梁的dx微元，由于dx很小，可視梁上載荷均勻分布，根據(jù)平衡方程ΣF=0，ΣMO=0，得：

(2)

圖1 簡支梁受力分析

(3)

由公式可看出，載荷q(x)是剪力Q(x)的一階導(dǎo)數(shù)，是彎矩M(x)的二階導(dǎo)數(shù)。在MATLAB中可用int和trapz函數(shù)求解曲線積分。其中int函數(shù)不僅能求解定積分，也可以求解不定積分,而trapz函數(shù)只能求解定積分。

int函數(shù)實(shí)現(xiàn)積分求解的語法是

(4)

式中，符號(hào)表達(dá)式為s關(guān)于變量x的不定積分；s為符號(hào)表達(dá)式；x為s的變量。

(5)

式中，符號(hào)表達(dá)式為s關(guān)于變量x的定積分；a、b分別為積分的上、下限。

trapz函數(shù)實(shí)現(xiàn)積分求解的語法是

(6)

式中，x、y分別為自變量和因變量。

在MATLAB中，還有一個(gè)函數(shù)cumtrapz可求解表達(dá)式或者函數(shù)的累計(jì)積分，實(shí)現(xiàn)積分求解的語法是：>>z=cumtrapz(x,y) .

(7)

利用int函數(shù)，可以快速求得積分上下限a、b上的剪力值和彎矩值。

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

求某一函數(shù)最值，數(shù)學(xué)上常用的方法是：對該函數(shù)求導(dǎo)，求出函數(shù)的極大值和極小值，與函數(shù)的邊界值進(jìn)行比較，即可得到函數(shù)最值。本文利用函數(shù)求導(dǎo)的辦法求解梁剪力和彎矩最大值，并找出最大值的位置。

在MATLAB中提供了求解多項(xiàng)式根的函數(shù)roots。r=roots(c)，c為一維向量，返回值為指定多項(xiàng)式的所有根(包括復(fù)根)。調(diào)用格式：

>>syms x;

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

>>c=sym2poly(y);

>>r=roots(c)

但函數(shù)roots只能求解多項(xiàng)式的根，對于復(fù)雜的方程只能通過編寫M文件的辦法求解。如參考文獻(xiàn)[11]利用二分法求根的方法編寫的M文件,bisect函數(shù)可較為精確地求解復(fù)雜方程的所有根。

利用MATLAB中的求根函數(shù)roots和bisect求解q(x)的所有根，然后將所得根代入剪力Q(x)，求得剪力Q(x)所有的極大值，與邊界值進(jìn)行比較，即可得出剪力最大值和最大值的位置。同理，可求出彎矩最大值及其位置。

2 工程應(yīng)用分析

2.1 工程背景及數(shù)值建模

焦煤集團(tuán)趙固一礦井田走向長2.0～5.5 km，傾斜寬9.5～11.0 km，井田面積約43.77 km2。主采煤層為山西組二1煤，屬近水平發(fā)育的穩(wěn)定型厚煤層，煤層傾角2°～6°，二1煤平均厚5.9 m。礦井采用立井單水平盤區(qū)開拓，井口標(biāo)高83.8 m，水平標(biāo)高-525 m，井深608.8 m。二1煤直接頂以泥巖為主，厚度1～6.5 m，平均厚度3.5 m；老頂厚度多為0.94～19.85 m,平均厚度13.5 m的粗、中、細(xì)粒砂巖；直接底為深灰色泥巖，平均厚3.4 m；老底為砂巖或砂質(zhì)泥巖，平均厚14.6 m。

利用有限差分軟件FLAC3D模擬趙固一礦井下開挖，記錄老頂砂巖層(厚度7.1 m)的垂直應(yīng)力?？紤]邊界效應(yīng)，設(shè)計(jì)的三維計(jì)算模型長×寬×高為200 m×200 m×72 m，共劃分77 400個(gè)六面體單元，生成網(wǎng)

格節(jié)點(diǎn)84 084個(gè)。X軸平行于工作面推進(jìn)方向，采用Mohr-Coulomb本構(gòu)模型。模型水平方向上左右兩側(cè)限制水平位移，模型底面限制水平移動(dòng)和垂直移動(dòng)，模型的上部施加上覆巖層的自重應(yīng)力。

2.2 離散數(shù)據(jù)的處理和求解

將FLAC3D中記錄的采場老頂砂巖層載荷數(shù)據(jù)導(dǎo)入MATLAB軟件中，畫出離散數(shù)據(jù)的散點(diǎn)圖[12]，觀察曲線形式(圖2)。由于回采工作面沿傾斜方向的長度遠(yuǎn)大于老頂沿走向懸露跨距，可將老頂視為一端由工作面煤壁、另一端由邊界煤柱支撐的固定梁。當(dāng)工作面推進(jìn)20 m時(shí)，判斷老頂砂巖是否發(fā)生破壞。

圖2 原始數(shù)據(jù)散點(diǎn)

由圖2可以看出，x距切眼從-2～6 m、15～22 m的曲線形式，與x從6～15 m的曲線形式明顯不同。因此，應(yīng)采用分段擬合曲線的辦法，即前9個(gè)點(diǎn)、后7個(gè)點(diǎn)、中間9個(gè)點(diǎn)分別進(jìn)行曲線擬合，比較擬合曲線的優(yōu)度。

采用MATLAB中的polyfit函數(shù)對各段進(jìn)行擬合，求得老頂砂巖層載荷q(x)的分段擬合函數(shù)為

(8)

繪制出擬合曲線圖(圖3)，可看出擬合曲線與實(shí)際吻合。

圖3 分段擬合曲線

用MATLAB軟件中的trapz函數(shù)分別對q(x)在各區(qū)間上積分，求出各區(qū)間內(nèi)的積分值。老頂砂巖為兩端固支狀態(tài)：

(9)

求出支座反力Q=57.544 kN，而砂巖的抗剪強(qiáng)度為37.64 MPa，意味著此時(shí)老頂砂巖在支承壓力作用下已發(fā)生破斷，砂巖層的受力狀態(tài)由兩端固支向簡支轉(zhuǎn)移，巖層中最大的彎曲拉應(yīng)力達(dá)到其抗拉強(qiáng)度是巖層由彎曲沉降發(fā)展至破壞的力學(xué)條件[4]。用int函數(shù)可求出彎矩在各區(qū)間內(nèi)的不定積分，然后根據(jù)邊界條件求出簡支梁條件下老頂砂巖的彎矩曲線函數(shù)。

圖4 剪力圖

由圖4可看出，剪力為零的點(diǎn)出現(xiàn)在[6,15]區(qū)間內(nèi)，即彎矩最大值出現(xiàn)在[6,15]區(qū)間內(nèi)，利用多項(xiàng)式求根函數(shù)roots，可求解剪力Q(x)為0的點(diǎn)即為彎矩最大值點(diǎn)。代入polyval函數(shù)，求出彎矩的最大值Mmax。程序如下：

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

>>r=roots(p);

r=8.67564

利用MATLAB多項(xiàng)式求值函數(shù)polyval，求出Mmax=-156 436.587 5 kN·m，計(jì)算出此時(shí)砂巖層最大彎曲拉應(yīng)力為4.76 MPa<[σt]=5.61 MPa，即此時(shí)砂巖層不會(huì)發(fā)生彎曲破斷。

當(dāng)工作面推進(jìn)21 m時(shí)，砂巖層彎曲拉應(yīng)力為5.12 MPa<[σt]，砂巖層不會(huì)發(fā)生彎曲破斷。

當(dāng)工作面推進(jìn)22 m時(shí)，砂巖層彎曲拉應(yīng)力為5.93 MPa>[σt]，砂巖層發(fā)生彎曲破斷?？梢钥闯錾皫r層應(yīng)在工作面推進(jìn)[21,22]區(qū)間內(nèi)發(fā)生彎曲破斷，由于工作面推進(jìn)距離較大，假定老頂砂巖層發(fā)生彎曲破斷時(shí)的載荷曲線與工作面推進(jìn)21 m時(shí)的載荷曲線一致，解出此時(shí)砂巖層發(fā)生破斷的推進(jìn)距離為L21=21.36 m；同理，假定老頂砂巖層發(fā)生彎曲破斷時(shí)的載荷曲線與工作面推進(jìn)22 m時(shí)的載荷曲線一致，解出此時(shí)砂巖層發(fā)生破斷的推進(jìn)距離為L22=21.78 m；二者求均值，得工作面老頂砂巖層彎曲破斷的推進(jìn)距離即老頂初次來壓步距：L=21.57 m，距切眼r=9.876 m處發(fā)生破斷。

2.3 對比線性載荷作用下梁的剪力和彎矩

用regress函數(shù)擬合出線性載荷q(x)曲線(圖5)。

圖5 線性分布的載荷曲線

用上述方法求出在線性載荷q(x)作用下工作面推進(jìn)20 m時(shí),趙固一礦老頂砂巖層最大剪力值Q線max=59.687 kN,最大彎矩值M線max=164.598 kN·m,最大彎矩值位置r線max=9.687 m。與實(shí)際載荷q(x)作用下工作面推進(jìn)20 m時(shí)的最大剪力值Qmax、最大彎矩值Mmax及最大彎矩值的位置rmax相比，誤差為

線性載荷作用下的老頂砂巖層初次來壓步距L線=19.26 m，而非線性載荷作用下的老頂砂巖層初次來壓步距L非線=21.37 m，誤差為9.87%。

3 結(jié) 論

(1)用MATLAB繪出基礎(chǔ)數(shù)據(jù)點(diǎn)集擬合的函數(shù)曲線，觀測其變化趨勢，同樣，將求出的應(yīng)力、彎矩等曲線繪出，直觀具體地表現(xiàn)地下工程復(fù)雜的力學(xué)特征。

(2)以梁模型為例，將離散的非線性數(shù)據(jù)點(diǎn)集擬合成可微積分的曲線函數(shù)，通過MATLAB中的積分運(yùn)算求解最大剪力和彎矩值，進(jìn)而可以對梁是以何種方式發(fā)生破壞(剪切破壞還是拉伸破壞)做出明確判斷，還可以通過求方程根的方法找出最大剪力值和最大彎矩值的位置，判斷梁在何處發(fā)生破壞。

(3)數(shù)據(jù)處理和計(jì)算都在MATLAB軟件中完成，不僅能夠保證計(jì)算工程的精度，而且能實(shí)現(xiàn)快速求復(fù)雜力學(xué)問題數(shù)值解的目的。

[1] 趙 凱,劉長武,張國良.用彈性力學(xué)的復(fù)變函數(shù)法求解矩形硐室周邊應(yīng)力[J].采礦與安全工程學(xué)報(bào),2007,24(3):361-365.

[2] 袁 林，高召寧，孟祥瑞.基于復(fù)變函數(shù)法的矩形巷道應(yīng)力集中系數(shù)黏彈性分析[J].煤礦安全，2013,44(2):196-200.

[3] 王振武，牛錚錚，馮秀苓.地下矩形硐室應(yīng)力分布的復(fù)變函數(shù)解[J].北華航天工業(yè)學(xué)院學(xué)報(bào),2010,20(4):1-6.

[4] 宋振騏.實(shí)用礦山壓力控制[M].徐州：中國礦業(yè)大學(xué)出版社，1988.

[5] 宋振騏，劉義學(xué)，陳孟伯，等.巖梁斷裂前后的支承壓力顯現(xiàn)及應(yīng)用的探討[J].山東礦業(yè)學(xué)院學(xué)報(bào)，1984(1):27-39.

[6] 王衛(wèi)軍,黃成光,侯朝炯.綜放沿空掘巷底鼓的受力變形分析[J].煤炭學(xué)報(bào),2002,27(1):26-30.

[7] 孟祥瑞,徐鋮輝,高召寧，等.采場底板應(yīng)力分布及破壞機(jī)理[J].煤炭學(xué)報(bào),2010,35(11):1832-1836.

[8] 喬立山，王玉蘭，曾錦光.實(shí)驗(yàn)數(shù)據(jù)處理中曲線擬合方法探討[J].成都理工大學(xué)學(xué)報(bào)：自然科學(xué)版,2004,31(1)：91-95.

[9] 魏國祥.基于MATLAB的橋梁健康監(jiān)測數(shù)據(jù)處理與可靠度分析[D].蘭州：蘭州交通大學(xué)，2012.

[10] Steven C，Chapra.工程與科學(xué)數(shù)值方法的MATLAB實(shí)現(xiàn)[M].唐玲艷，田尊華譯.2版.北京：清華大學(xué)出版社，2009.

[11] 劉鴻文.材料力學(xué)[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)

郭勤強(qiáng)(1990—)，男，碩士，471000 河南省洛陽市洛龍區(qū)龍門大道573號(hào)。