陳元非，王 磊

(1. 池州學院 地理與規(guī)劃學院，安徽 池州 247000；2. 安徽理工大學 空間信息與測繪工程學院，安徽 淮南 232001)

0 引言

FLAC3D數(shù)值模擬軟件是美國明尼蘇達Itasca軟件公司開發(fā)的三維顯式有限差分程序[1]，并首先將其應(yīng)用到巖土工程力學計算中，現(xiàn)已成為巖土力學計算中重要的數(shù)值計算方法，在采礦、土木、水利工程等領(lǐng)域有著廣泛的應(yīng)用[1]。FLAC3D軟件本構(gòu)模型豐富，可以模擬各種彈塑性模型和開挖模型；具有強大的后處理功能，可以輸出位移、應(yīng)力、塑性區(qū)分布等各種信息，并且提供了內(nèi)嵌的FISH語言，方便用戶進行二次開發(fā)。但FLAC3D軟件應(yīng)用上也存在明顯的短板，一是采用命令流的方式，二是在模型構(gòu)建上的困難，一直是困擾科研工作者的難題之一[2]。

對此，國內(nèi)外也有學者展開了復雜FLAC3D模型建模和預(yù)處理的研究。高永剛[2]通過分析AutoCAD二維平面圖形與FLAC3D網(wǎng)格單元的關(guān)系，提出了基于AutoCAD快速建立FLAC3D網(wǎng)格模型的方法。李根[3]等利用ANSYS軟件的強大建模功能，驗證了ANSYS模型導入FLAC3D的具體操作可行性。羅周全[4]、林杭[5]提出基于SURPAC軟件的直觀、快速構(gòu)建FLAC3D模型的技術(shù)。荊永濱[6]等利用TetGen網(wǎng)格生成工具實現(xiàn)了FLAC3D模型的可視化構(gòu)建和交互式操作。馮艷順[7]嘗試了GIS結(jié)合EXCEL進行快速建模的方法，實現(xiàn)了礦區(qū)地層快速建模。劉秀軍[8]利用GOCAD 建模軟件，實現(xiàn)了將復雜三維地質(zhì)模型導入FLAC3D軟件中?？梢钥闯?，大量的復雜FLAC3D模型構(gòu)建的研究基本是借助于第三方軟件平臺，因此用戶需要在建模前對幾種軟件的使用較為熟悉，掌握建模方法并了解其與FLAC3D建模的區(qū)別。這無疑增加了額外的負擔和使用FLAC3D軟件的難度，打消了用戶學習和使用軟件的積極性[2]。

開采沉陷問題的FLAC3D建模具有它的特殊性，例如開采工作面一般較為規(guī)整，工作面采厚基本不變，且位于煤層中。因此可以嘗試利用簡單的六面體網(wǎng)格構(gòu)建模型，而不必要引入第三方軟件進行復雜的模型構(gòu)建。為了滿足開采沉陷問題研究中FLAC3D建模應(yīng)用方面的具體需要，同時盡可能減小FLAC3D建模難度，降低使用門檻，本文提出一種參數(shù)化構(gòu)建FLAC3D模型的方法。該方法利用大量實測鉆孔數(shù)據(jù)，通過直接參數(shù)輸入的方式，自動生成模型構(gòu)建的命令流文件，避免了命令流代碼編輯及后續(xù)建模上的困難，可以較為準確的模擬地下礦山開采實際的三維地層模型，對于傾斜煤層、變化厚度的巖層、山區(qū)地表等復雜三維模型的開采沉陷預(yù)測和規(guī)律的研究提供了方法，也對礦山開采沉陷虛擬仿真教學具有一定的參考價值。

1 FLAC3D建模方法

1.1 FLAC3D均勻格網(wǎng)建模方法

在利用FLAC3D建模進行礦區(qū)開采沉陷預(yù)計問題中，一般的會根據(jù)地表實測大量的鉆孔數(shù)據(jù)，選擇一個具有代表性的鉆孔進行巖層簡化合并，最終確定各個巖層的厚度。因此，在大量鉆孔數(shù)據(jù)被遺棄、不同巖層被簡化合并的基礎(chǔ)上，構(gòu)建出的FLAC3D地層模型往往不能準確模擬該地區(qū)實際地層情況，計算結(jié)果也因此具有較大偏差。

為了盡可能多的利用上每個鉆孔的實測數(shù)據(jù)，結(jié)合開采沉陷研究領(lǐng)域的特點，將整個模型從平面上進行單元均勻格網(wǎng)化，即將模型底面進行均勻格網(wǎng)劃分，通過插值計算得到各層格網(wǎng)節(jié)點的坐標，再進行組合建模的思路。利用單元網(wǎng)格模型采用brick六面體中8個角點定位的方式(如圖1所示)，構(gòu)建不同形狀的單元格網(wǎng)模型。再利用單元格網(wǎng)模型之間的位置邏輯關(guān)系，進行組合計算和建模，如圖2所示。

圖2 相鄰格網(wǎng)間的連接

在圖1中，直接由8個角點確定一個單元brick六面體模型，它表示某一巖層的單元矩形格網(wǎng)模型。P0-P1P2P3方向為X、Y、Z軸方向。單元格網(wǎng)模型的尺寸可以設(shè)定，根據(jù)設(shè)置的整體模型XYZ方向的總尺寸和SIZE參數(shù)，確定各單元格網(wǎng)的尺寸。對于圖1所示單元格網(wǎng)模型，設(shè)置得該巖層單元格網(wǎng)模型SIZE=1,1,3。

多個單元格網(wǎng)模型組合形成整體模型，各單元之間有其緊密連接的條件。如圖2所示，單從X方向上看，單元1和單元2有4個點共用的約束條件，分別是[1-P1；2-P0]；[1-P4；2-P2]；[1-P6；2-P3]；[1-P7；2-P5]。因此Y方向上各個相鄰單元之間的約束條件分別是[1-P2；3-P0]；[1-P4；3-P1]；[1-P5；3-P3]；[1-P7；3-P6]。從Z方向上，下一層巖層的頂面又是上一層巖層的底面，共用下一巖層的P3、P5、P6、P7四個點。因此，采用均勻格網(wǎng)的方法建模時，各節(jié)點XY坐標的計算較為簡單，可以直接通過均勻格網(wǎng)的劃分得到。Z坐標則可以采用多鉆孔數(shù)據(jù)中的巖層厚度信息經(jīng)過空間插值的方法計算。這樣就可以構(gòu)建由眾多單元格網(wǎng)模型組合構(gòu)成的FLAC3D整體模型。

1.2 網(wǎng)格節(jié)點巖層厚度數(shù)據(jù)插值

多鉆孔數(shù)據(jù)的建模方法的思想是將整體模型格網(wǎng)化，分解成一個個便于建立的單元模型，再將單元模型組合成一個整體模型的思路。其中最為重要的工作是準確插值計算每個單元模型中各個節(jié)點的Z值。對于插值方法國內(nèi)外已經(jīng)有過大量研究。其中比較典型的插值方法有：克里金(Kriging)插值法、反距離加權(quán)插值法(IDW)、樣條插值法(Spline)等。其中克里金插值法在地質(zhì)探礦、地理空間插值、地學統(tǒng)計與建模中應(yīng)用十分廣泛[9,10]。因此，本研究采用克里金插值作為鉆孔巖層數(shù)據(jù)的插值方法。

多鉆孔數(shù)據(jù)建模工作主要是計算單元格網(wǎng)模型的8個角點的三維坐標，其難點在于插值各巖層上界面4點的Z坐標。通過單元格網(wǎng)的底面Z坐標，利用克里金插值方法插值出該節(jié)點處的各個巖層厚度，并從下往上相加即可。底面4個節(jié)點則可以根據(jù)相鄰下一巖層單元模型的頂面推算，以此推算出所有節(jié)點的Z值。采用克里金插值方法插值的流程如圖3所示。

圖3 克里金插值法對任一單元節(jié)點P的層厚計算流程

1.3 參數(shù)化建模方法

均勻格網(wǎng)的建模方法可以實現(xiàn)基于多鉆孔數(shù)據(jù)的巖層模型建模，但是整個建模過程十分繁瑣，需要給每一巖層的每一個單元格網(wǎng)模型進行建模、組合、建立接觸面、參數(shù)賦值、開挖、求解等過程，并且涉及到網(wǎng)格節(jié)點的巖層厚度插值計算等內(nèi)容。因此僅僅利用手工計算編輯命令流的方式是不可能實現(xiàn)的。為了降低使用難度，結(jié)合EXCEL VBA程序編譯功能，對程序功能進行封裝，用戶可以輸入建模需要主要參數(shù)，從而直接生成FLAC3D建模命令流文件，實現(xiàn)建模計算功能的參數(shù)化。在利用程序生成命令流文件時，需要輸入的主要參數(shù)如表1所示。

表1 程序設(shè)計輸入?yún)?shù)內(nèi)容

參數(shù)化建模過程直接利用命令流文件中網(wǎng)格坐標數(shù)據(jù)之間的邏輯聯(lián)系，通過程序設(shè)計手段，生成完整的命令流文件。由于整個模型構(gòu)建是由大量單元格網(wǎng)模型的組合形成，因此文件中會為每一單元格網(wǎng)模型編寫一行由brick六面體構(gòu)建的相應(yīng)的命令流。最終完成整體模型命令流文件的生成。以平面5×5格網(wǎng)建模為例，其建模過程如圖4所示。

圖4 建模流程示意圖

(1) 首先建立最底層初始單元格網(wǎng)模型a，單元網(wǎng)格模型8個角點坐標中，平面坐標由模型X、Y方向坐標與X、Y單元格數(shù)SIZE劃分確定。模型Z坐標由鉆孔數(shù)據(jù)中底層巖層厚度數(shù)據(jù)進行插值計算得到；

(2) 依次沿X、Y軸方向建立其它底層單元格網(wǎng)模型b～f；

(3) 繼續(xù)依次沿Z軸方向建立其它巖層單元格網(wǎng)模型g。該層單元格網(wǎng)模型底部Z值與相鄰下部巖層模型該點頂部Z值相等；頂部Z值由鉆孔數(shù)據(jù)中各層巖層厚度數(shù)據(jù)進行插值，并從底部向上累計計算得到；

(4) 最終完成整體模型構(gòu)建h。

2 工程案例

現(xiàn)模擬礦區(qū)計劃開采一個走向長2000 m，傾向長300 m的工作面。在該區(qū)域范圍內(nèi)，共有巖層鉆孔數(shù)據(jù)123個。該區(qū)域巖層從下至上分別為粉砂巖(底板)、煤層、中砂巖 、細砂巖 、表土層共5層，煤層厚度4 m。各巖層的力學參數(shù)可以參考以往歷史數(shù)據(jù)?，F(xiàn)可以根據(jù)需求構(gòu)建一個平面大小為4000 m×4000 m區(qū)域的模型。模型開挖煤層中范圍為1000 m

利用Excel vba編譯的程序自動生成FLAC3D命令流文件，并用FLAC3D軟件調(diào)用命令流文件，生成礦區(qū)該地區(qū)三維地層模型，如圖5所示。可以看出，該地區(qū)地表有一定的起伏狀態(tài)，地表點距離模型底部高度在170～500 m距離不等，坡度在1/10左右。為分析數(shù)據(jù)插值效果，對比了直接利用鉆孔數(shù)據(jù)的Surfer8繪制的地表圖，如圖6所示。開挖工作面后最終計算的下沉和水平移動云圖如圖7所示。該模型在Core i7、32G內(nèi)存的計算機上運行，代碼自動生成和建模過程用約3 min，模型解算過程約用1 h。可以看出參數(shù)化建模方法建模過程所消耗的時間相對于整體解算過程消耗時間可以忽略，驗證了該方法的效率。

圖5 礦區(qū)FLAC3D三維模型

圖6 程序插值地表(a)與Surfer8軟件繪制的地表云圖(b)對比

圖7 模型開挖之后地表下沉(a)及水平移動云圖(b、c)效果

從模型建立和開挖之后的變形結(jié)果來看，模型地表的三維結(jié)果與Surfer8基于鉆孔數(shù)據(jù)的插值結(jié)果一致。開采變形運行結(jié)果沉降最大值約為3.5 m，最大水平移動值約為0.8 m，下沉系數(shù)和水平移動系數(shù)分別為0.88和0.22，符合一般礦區(qū)的經(jīng)驗，驗證了該方法可靠性和有效性。

3 結(jié)論

(1) 基于開采沉陷預(yù)測領(lǐng)域FLAC3D建模問題進行分析，提出利用均勻格網(wǎng)分割、插值計算、組合為一體的參數(shù)化建模方法，實現(xiàn)了多鉆孔實測數(shù)據(jù)的建模應(yīng)用?？梢暂^為準確的模擬地下礦山開采實際的三維地層模型，對于傾斜煤層、變化厚度的巖層、山區(qū)地表等復雜三維模型的開采沉陷預(yù)測和規(guī)律的研究提供了方法，也對礦山開采沉陷虛擬仿真教學具有一定的參考價值。

(2) 針對網(wǎng)格節(jié)點坐標計算和組合問題采用克里金插值方法插值巖層厚度數(shù)據(jù)，并通過編程實現(xiàn)了克里金插值、參數(shù)化建模功能，避免了人工處理大量鉆孔數(shù)據(jù)編輯命令流文件的困難。

(3) 利用礦區(qū)實測鉆孔數(shù)據(jù)進行參數(shù)化建模應(yīng)用，通過分析三維模型和變形模擬結(jié)果，驗證了該方法的可靠性和有效性。