基于FLAC3D橫觀各向同性模型的煤礦井田初始地應(yīng)力場反演方法

余大軍，楊張杰，郭運華，楊永剛，王 波

(1.安徽省煤炭科學(xué)研究院 礦山支護中心，安徽 合肥 230001; 2.武漢理工大學(xué) 道路橋梁與結(jié)構(gòu)工程湖北省重點實驗室，湖北 武漢 430070)

巖體初始應(yīng)力場是巖石工程數(shù)值模擬計算與工程穩(wěn)定性分析的必要條件，也是巖石工程設(shè)計與施工的重要依據(jù)。地應(yīng)力是煤礦開采等地下工程圍巖變形與破壞的根本驅(qū)動力，因此應(yīng)力場研究是巖石力學(xué)與采礦工程中最基本、最重要的工作之一[1]。在地應(yīng)力研究初期，基本上將地應(yīng)力簡單地處理為求解側(cè)壓力系數(shù)，復(fù)雜地質(zhì)條件下誤差較大，根本不適合實際使用。1983年，郭懷志等[2]提出了采用實際工程地應(yīng)力實測資料，運用有限元數(shù)學(xué)模型回歸分析方法求解出巖體初始應(yīng)力場。具體作法是依據(jù)彈性力學(xué)應(yīng)力疊加原理，將復(fù)雜的邊界應(yīng)力條件分解為幾個簡單的邊界條件的線性組合，通過回歸方法求解組合系數(shù)來實現(xiàn)，奠定了多元回歸方法擬合初始地應(yīng)力方法的基礎(chǔ)。此后，諸多學(xué)者開始研究地應(yīng)力場的反分析方法，朱伯方[3]提出通過利用地下工程掘進過程中實測位移進行小范圍的地應(yīng)力反演;2000年以后，大量文獻引入遺傳算法與神經(jīng)網(wǎng)絡(luò)等優(yōu)化算法尋求地應(yīng)力反演最優(yōu)解[4-10]，在諸多巖土工程的地應(yīng)力場反演中廣泛使用這些反演方法，巖體被簡化為理想的各向同性材料，采用彈性模型分析，基本都能得出比較理想的結(jié)果[11-12]。為了進一步提高反演精度，郭運華等[13]提出了一種基于FLAC3D改進的初始地應(yīng)力場回歸方法，該方法對各應(yīng)力分量分別進行回歸，以提高局部異常應(yīng)力區(qū)域的擬合精度。

近年來，地應(yīng)力反演方法也應(yīng)用于煤礦領(lǐng)域，取得了一定的成效[14-17]。對于煤層巖體，有其特殊的層狀結(jié)構(gòu)，特別是多層構(gòu)造的煤系地層其力學(xué)性質(zhì)在垂直和平行于巖層走向的兩個方向上，差異很大，具有橫觀各向同性彈性體典型的特征。因此對于初始地應(yīng)力場的反演，可以將這種深部層狀構(gòu)造煤巖體假設(shè)為橫觀各向同性彈性體進行研究。筆者近幾年在淮南礦區(qū)，接觸到了諸多典型的層狀結(jié)構(gòu)煤層，含煤地層一般是沉積巖，除煤層外，地層中含有多層巖石夾層，其最常見的巖性有泥巖和砂巖，其厚度從幾厘米至幾米，幾十米不等，復(fù)雜結(jié)構(gòu)含煤地層中夾層數(shù)一般多達幾十層，而且?guī)r層一般具有一定傾角，傾角隨著巖層變化，巖層的力學(xué)性質(zhì)在方向上也不相同，有的甚至顯著不同，這樣的層狀結(jié)構(gòu)增加了求取力學(xué)參數(shù)和進行數(shù)值模擬的難度[18]。因此，對具有典型橫觀各向同性巖體的應(yīng)力場研究很有現(xiàn)實意義。陳宇龍和張宇寧[19]研究了橫觀各向同性巖石在不同層面傾角條件下的單軸壓縮破壞過程及聲發(fā)射特性，結(jié)果顯示互層巖體有特殊的破壞形式。張志增[20]對橫觀各向同性巖體做了系統(tǒng)的研究，提出了一種用于反演橫觀各向同性巖體地應(yīng)力的方法，該方法通過巷道內(nèi)巖體位移來反演橫觀各向同性巖體地應(yīng)力，通過測量巖體位移，根據(jù)巷道位移解析公式得到地應(yīng)力表達式。

煤系地層由于厚度及傾角變化等，通過實測得到各層巖體的力學(xué)參數(shù)具有一定的難度。筆者介紹了一種基于FLAC3D橫觀各向同性模型，通過實測地應(yīng)力參數(shù)，對具有板狀巖層構(gòu)造的礦區(qū)進行地應(yīng)力場反演的方法。

1 基本原理

1.1 多元回歸原理

多元回歸方法是建立在彈性力學(xué)應(yīng)力疊加原理基礎(chǔ)上，在巖土天然應(yīng)力狀態(tài)下，可以認為巖土體處于平衡狀態(tài)且無持續(xù)的塑性變形的變化，將地質(zhì)體所受的復(fù)雜地應(yīng)力分解為圖1中6種簡單的邊界應(yīng)力形式，通過在邊界施加單位荷載獲得內(nèi)部單元的基本初始應(yīng)力，然后將基本初始應(yīng)力當(dāng)做自變量，實測地應(yīng)力當(dāng)做因變量進行偏最小二乘回歸分析，最終求解回歸系數(shù)。

假設(shè)經(jīng)過復(fù)雜地質(zhì)運動后的地應(yīng)力由圖1中6種簡單地質(zhì)運動疊加影響形成，則地應(yīng)力可由式(1)得到。將實測結(jié)果與圖1中6種簡單加載方式求得的相同位置的地應(yīng)力結(jié)果進行回歸分析，可以求得式(1)中的系數(shù)，將實測值表示為幾種簡單加載方式計算值的組合，回歸方程為

圖1 邊界應(yīng)力分解組合示意Fig.1 Decomposition of boundary stress

σd=l1σx+l2σy+l3σm+l4τxy+l5τxz+l6τyz+e

(1)

式中，σd為初始地應(yīng)力;σm為自重應(yīng)力引起的應(yīng)力;σx，σy，τxy，τxz，τyz分別為水平x向擠壓、水平y(tǒng)向擠壓、xy平面內(nèi)剪切、xz平面內(nèi)剪切、yz平面內(nèi)剪切運動引起的應(yīng)力;e為常數(shù);l1～l6為回歸系數(shù)。

1.2 FLAC3D橫觀各向同性模型

大型三維軟件FLAC3D已經(jīng)廣泛地應(yīng)用于巖土力學(xué)工程分析中，F(xiàn)LAC3D程序自帶有橫觀各向同性彈性材料、莫爾-庫侖彈塑材料和應(yīng)變軟化和硬化塑性材料等多種本構(gòu)模型。橫觀各向同性彈性材料本構(gòu)模型主要模擬層狀結(jié)構(gòu)的彈性介質(zhì)，模型中垂直于層狀介質(zhì)的彈性模量和平行于介質(zhì)的彈性模量明顯不同。含煤地層的沉積巖正符合了橫觀各向同性彈性材料模型的特點。FLAC3D橫觀各向同性彈性材料模型包含了5個獨立的彈性常量和2個巖層物理量，其中5個獨立彈性常量為各向同性平面內(nèi)的彈性模量、垂直于各向同性平面內(nèi)彈性模量、垂直于各向同性平面內(nèi)的剪切模量以及2個不同平面內(nèi)的泊松比，2個巖層物理量為各向同性平面的傾向、各向同性平面的傾角。本文通過數(shù)值試驗得到2個方向的彈性模量、剪切模量以及泊松比;通過截面巖層的數(shù)學(xué)關(guān)系計算得到各向同性平面平均的傾向和傾角與所建坐標系的關(guān)系式。

2 應(yīng)力場加載至模型

通過多元回歸擬合，得到地應(yīng)力回歸方程。地應(yīng)力加載即將滿足回歸方程的應(yīng)力分量載入模型單元及求解邊界節(jié)點力實現(xiàn)模型力平衡，用于下一步分析設(shè)計、數(shù)值計算的過程。在已有的地應(yīng)力場回歸反演文獻中，基本上都只介紹了求解回歸方程的過程，而對地應(yīng)力場如何加載至模型很少提及。郭運華等[13]提出如下一套加載方案:通過上述提到的6種簡單的加載方式，得到每種加載方式下模型單元的應(yīng)力分量，并將其作為回歸方程的一個自變量，實測值作為因變量進行偏最小二乘回歸，求解得到回歸方程系數(shù)。按回歸方程計算各單元地應(yīng)力分量，并賦值給模型單元應(yīng)力分量。此時模型單元只有內(nèi)力，節(jié)點力均為0，并不能達到平衡。此時需要計算模型節(jié)點不平衡力，將求得的節(jié)點荷載(與不平衡力大小相等，方向相反)加載至模型中，這樣便滿足了模型的邊界應(yīng)力條件，同時應(yīng)力場也滿足實測地應(yīng)力回歸估計值。

3 潘一礦東區(qū)地應(yīng)力反演分析

3.1 工程概況

潘一礦東區(qū)位于安徽省淮南市北部潘集區(qū)。南以淮河與淮南老礦區(qū)相隔，西與張謝礦區(qū)相鄰。潘一礦井東起第0勘探線與潘二礦毗鄰，西至第Ⅸ勘探線及技術(shù)邊界線與潘三礦相接，北以F2、潘集背斜軸、F4、F5、F5-1斷層與潘二、潘北礦為界，南至13-1煤層-800 m底板等高線地面投影線。礦井東西走向長14.6 km，傾斜寬約4 km，面積約58.4 km2，如圖2所示。

圖2 潘一東礦區(qū)范圍示意Fig.2 Map of Panyi east mining area

潘一礦東區(qū)位于潘集背斜及東部傾伏轉(zhuǎn)折端東南翼。總體形態(tài)為褶曲構(gòu)造，發(fā)育縱貫全區(qū)的軸向東西轉(zhuǎn)南東向的潘集背斜，其間次級褶曲較為發(fā)育;區(qū)內(nèi)共發(fā)育斷層142條，其中正斷層113條，逆斷層29條，構(gòu)造復(fù)雜程度中等，斷層的走向以北西西向為主，其次為北東向等;斷層傾向以南南西向為主，其次為北西向等，筆者主要考慮了規(guī)模較大的F32，F(xiàn)S2，F(xiàn)S8，F(xiàn)S9四條斷層，其走向基本為EW向，傾向為S或N向，傾角40°～70°，如圖2所示。已有研究成果[4]顯示淮南礦區(qū)應(yīng)力場以水平應(yīng)力為主，構(gòu)造應(yīng)力占絕對優(yōu)勢。最大主應(yīng)力為NEE—EW向，部分區(qū)為NW—SE向。

根據(jù)潘一礦東區(qū)煤層及巷道的布置以及初始地應(yīng)力測點的位置，本次模擬選取包含了大部分采區(qū)區(qū)域，礦區(qū)寬度3 500 m，長度5 000 m，深度1 300 m，坐標原點取在整個礦區(qū)右下角，y軸沿正北方向，x軸沿正東向，z軸沿深度方向，向上為正。

3.2 地質(zhì)模型概化

勘查報告顯示潘一東區(qū)地層自老至新有石炭系、二疊系、第四系，主要是由30多層層厚從1 m至100多米的砂質(zhì)泥巖、粉砂巖、中砂巖、泥巖、細砂巖、煤層組成，地層綜合巖性見表1。數(shù)值模型中很難詳盡地反映地層情況，有些巖層厚度0.5～2.0 m，而且?guī)r性交互變化，因此必須找到一種有效可行的簡化建模方法。筆者將力學(xué)性質(zhì)相近的地層合并，對于層厚較小、層數(shù)較多且層面平行的多種巖性交互的地層，則采用等效模型、等效參數(shù)代替。具體做法是:首先將力學(xué)性質(zhì)相近、層厚較小的巖層合并;從地面至深部，可按照上覆層(包括強風(fēng)化巖層)、巖性及交互規(guī)律一致的巖層、新巖體等概化地層;同時模型也主要考慮了規(guī)模較大的F32，F(xiàn)S2，F(xiàn)S8，F(xiàn)S9四條斷層。通過數(shù)值模擬計算得到上述交互巖層等效的各向同性和橫觀各向同性材料參數(shù)。根據(jù)通過上述方法，可得潘一礦東區(qū)地層概化模型如圖3所示。

圖3 潘一東礦區(qū)地質(zhì)概化分區(qū)Fig.3 Geological overview division of Panyi east mining area

淮南礦業(yè)集團及有關(guān)單位對潘一礦東區(qū)地應(yīng)力場分析做了大量實測工作，其地應(yīng)力特征實測數(shù)據(jù)轉(zhuǎn)換為模型坐標系統(tǒng)見表2。

表2 地應(yīng)力實測值轉(zhuǎn)換計算結(jié)果Table 2 Results of the conversion of measured stress values

3.3 材料等效參數(shù)計算

自然界中大部分巖石都是各向異性的，特別是沉積巖具有顯著的層狀構(gòu)造特征。目前在理論上可將層狀巖體處理成橫觀各向同性介質(zhì)，其力學(xué)關(guān)系包含5個彈性常數(shù)[20]。

潘一東礦區(qū)地層主要是由30多層層厚從1 m至100多米的砂質(zhì)泥巖、粉砂巖、中砂巖、泥巖、細砂巖、煤層組成。地層綜合巖性見表1。為了探索解決這種多層巖體參數(shù)的問題，本文采用數(shù)值模擬的方法，計算得到巖層的等效參數(shù)[21]。

表1 地層綜合巖性Table 1 Comprehensive lithology

首先根據(jù)各巖層的力學(xué)性質(zhì)、巖性互層規(guī)律等進行概化分區(qū)，將礦區(qū)的巖層自上而下簡化分為6層，見表3。根據(jù)本文等效假設(shè)，再對每一層進行數(shù)值加載試驗，得出各巖層的等效參數(shù)。

表3 巖層分區(qū)簡化Table 3 Simplified rock partition

以第3層為例，該層層厚約150 m，由粉砂巖與砂質(zhì)泥巖互層組成，傾角隨著位置而變化。巖層累計10層，其每層的厚度4～10 m，建立數(shù)值計算模型，如圖4所示。

圖4 層狀巖體等效模型Fig.4 Equivalent model of layered rock mass

模型中各巖層參數(shù)分別由實驗獲得，分別通過對模型軸向加載試驗、橫向加載試驗等數(shù)值試驗(圖5)，得出模型的等效參數(shù)，見表4,5。

圖5 加載試驗應(yīng)力應(yīng)變曲線Fig.5 Stress and strain curves of loading test

表4 各向同性材料力學(xué)參數(shù)Table 4 Mechanical parameters of isotropic materials

表5 橫觀各向同性材料力學(xué)參數(shù)Table 5 Mechanical parameters of transverse isotropic materials

另外通過截面巖層的數(shù)學(xué)關(guān)系計算得到各向同性平面平均的傾向和傾角與所建坐標系的關(guān)系式。

FLAC3D橫觀各向同性材料模型中包括2個物理常數(shù)，分別為各向同性巖層面的傾向和傾角。研究范圍內(nèi)，巖層的傾向SE—SW，傾角為5°～23°。為了簡化計算，巖層的傾向取優(yōu)勢方向0°，傾角根據(jù)yz平面內(nèi)巖層線的幾何關(guān)系求得，如圖6所示。

圖6 巖層yz平面剖面Fig.6 Section of the rock formation in the yz plane

第3層巖層線在yz平面內(nèi)擬合的方程式為:z=-8×105y2+0.297y+65.53，則巖層的傾角為:θ=-16×10-5y+0.297，rad，其中y為模型中單元的y軸坐標值。通過傾角的關(guān)系式即可編寫接口文件進行計算。

斷層的寬度根據(jù)經(jīng)驗及斷層附近巷道掘進過程中揭露的斷層影響范圍，取為40 m，其力學(xué)參數(shù)為:彈性模量為5 GPa，泊松比為0.35，密度為2 000 kg/m3。至此，全部參數(shù)都已求得，可進行地應(yīng)力反演計算。

3.4 地應(yīng)力反演

通過以上得到的橫觀各向同性模型中7個參數(shù)，編制FLAC3D計算文件，進行地應(yīng)力反演計算。整個反演過程通過編寫的接口文件讀寫操作完成。得到橫觀各向同性條件下對應(yīng)的地應(yīng)力回歸方程為

σd=2.334σx+2.5σy+1.357σm+1.703τxy+

6.576τxz-8.318τyz+7 611 039.518

反演結(jié)果表明初始地應(yīng)力90%以上的變動都可以被該模型所解釋，擬和度較高。由此表明顯著性檢驗是合理的。經(jīng)應(yīng)力場平衡計算，得到計算坐標系下測點的地應(yīng)力實測值與回歸值對比見表6，7。

表6 測點實測值與回歸值對比Table 6 Comparison of measured and regression values of points MPa

表7 測點實測值與回歸值對比Table 7 Comparison of measured and regression values of points

由表6,7可知，采用橫觀各向同性材料進行反演計算的擬合結(jié)果誤差較小，大部分回歸值與實測值之間的絕對誤差在30%以內(nèi)，這些誤差產(chǎn)生的原因主要是簡化計算模型與實際地質(zhì)原型之間的差異，需要通過細化計算模型、詳盡地質(zhì)勘察來減少誤差。結(jié)果還表明采用橫觀各向同性模型進行計算更能體現(xiàn)具有層狀構(gòu)造的煤系地層力學(xué)性質(zhì)的方向差異性，能提高局部區(qū)域的反演精度。

按回歸方程計算得到各單元應(yīng)力值并加載至模型，得到精確滿足實測地應(yīng)力分布的地應(yīng)力場。通過編寫接口文件，計算得到-850 m水平標高平面主應(yīng)力等值線圖(圖7)、-850 m水平標高平面最大主應(yīng)力傾角云圖(圖8)和-850 m水平標高平面最大主應(yīng)力方位角云圖(圖9)。

圖7 -850 m水平標高平面最大主應(yīng)力等值線Fig.7 -850 m contour map of maximum principal stress in depth

圖8 -850 m水平標高平面最大主應(yīng)力傾角云圖Fig.8 -850 m maximum principal stress inclination cloud indepth

圖9 -850 m水平標高平面最大主應(yīng)力方位角云圖Fig.9 -850 m maximum principal stress azimuth in depth

由圖7可知，-850 m水平標高平面巖層的最大主應(yīng)力值為20～45 MPa，大部分區(qū)域最大主應(yīng)力在35～40 MPa，屬高應(yīng)力水平;總體來看東北區(qū)域和西南區(qū)域較西北、東南區(qū)域最大主應(yīng)力要稍大。分析原因，所研究的4條斷層所在區(qū)域與應(yīng)力降低的區(qū)域基本一致，進一步說明斷層對地應(yīng)力的影響。

將最大主應(yīng)力傾角定義為最大主應(yīng)力矢量與xy平面的夾角(向上為正)，方位角定義為沿正北向(y軸正向)順時針旋轉(zhuǎn)至最大主應(yīng)力矢量在xy平面上的投影的角度(0°～360°)。由圖9可知，最大主應(yīng)力方位角優(yōu)勢方向為100°～150°和300°～350°，即主應(yīng)力的優(yōu)勢方向為NW-SE向。圖8表明最大主應(yīng)力傾角為-24°～32°，大部分傾角為15°～30°，向下傾，近似水平。因此可知最大水平主應(yīng)力即最大主應(yīng)力，在布置巷道及煤層開采走向設(shè)計時應(yīng)予以考慮。

4 結(jié) 論

(1)選擇了分區(qū)域橫觀各向同性模型與各向同性模型相結(jié)合的計算方法，反演精度進一步得到提升。

(2)通過煤層分布的幾何關(guān)系，對煤層傾角作了動態(tài)計算，有利于提高數(shù)值計算精度。

(3)通過對復(fù)雜的地層進行地質(zhì)概化，建立與實際巖層一致的計算模型，進行數(shù)值試驗，得到了多層巖體等效力學(xué)參數(shù)，該方法優(yōu)化了建模過程，提高了數(shù)值計算效率。