朱權潔，劉金海，王勝開，高林生，梁　娟

(1. 華北科技學院 安全工程學院，北京 東燕郊 101601；2. 馳宏鋅鍺股份有限公司，云南 曲靖 655011；3. 防災科技學院，北京 東燕郊 101601)

0　引言

礦井地應力場測量與反演是確定工程巖體力學屬性、進行圍巖穩(wěn)定性分析以及地下工程開挖設計科學化的必要前提[1-2]，也是有效防治沖擊地壓等災害、保證礦井安全有序生產(chǎn)的關鍵。目前，地應力測試都是在硬巖層中進行的，煤層應力值除可根據(jù)已知點對應巖層地應力實測結果進行估算外，目前尚無其他較好辦法，但這種方法也具有局限性和不穩(wěn)定性。要全盤掌握礦區(qū)內(nèi)不同深度、不同位置、同一煤層不同深度上的應力分布，需進行區(qū)域地應力場的反演分析，這其中地質模型的數(shù)值計算是關鍵步驟。

數(shù)值計算是當前解決科研難題、實際工程問題的有效手段[3-5]。在巖土工程領域，有限元法、離散元法以及有限差分法等是數(shù)值分析計算的重要方法之一。Itasca公司研發(fā)的FLAC3D數(shù)值模擬軟件是其中的優(yōu)秀代表之一，該軟件對于模擬不產(chǎn)生離散、連續(xù)變形模型具有良好效果，因此，在巖土工程、采礦工程等行業(yè)和領域應用廣泛。

在具有強大數(shù)值分析能力的同時，F(xiàn)LAC3D的建模尤其是復雜模型的構建能力是其最大短板。對于規(guī)則模型而言，利用fish語言編寫命令流可實現(xiàn)快速建模，但不規(guī)則、復雜的地質模型依靠命令流實現(xiàn)是非常困難的，這也是FLAC3D面臨的一大挑戰(zhàn)[6]?，F(xiàn)有研究成果中多選擇“先外部構建幾何模型、后導入FLAC3D”的方式開展[7-9]。如林杭等編制了Surpac到FLAC3D的接口程序，直接實現(xiàn)Surpac三維模型到FLAC3D的轉換；紀洪廣等[[10]提出一種利用VC++開發(fā)的復雜礦山地質形體網(wǎng)格精細劃分方法，解決了FLAC3D構建復雜地質體費時、繁瑣等缺陷；趙文等[11]提出利用MATLAB和Surfer劃分網(wǎng)格，并生成FLAC3D命令流語句，提高了FLAC3D建模速度。這些方法中前者在網(wǎng)格劃分方面存在缺陷，后者對使用者的編程能力要求較高。

此外，袁海平等[12]提出聯(lián)合AutoCAD，ANSYS和FLAC3D聯(lián)合構建三維地質模型，該方法在CAD中制作面元素，仍然受限于構建大型復雜三維模型；李翔等[13]基于鉆孔柱狀圖資料，提出利用SURPAC，CAD和ARCGIS建立復雜三維地址模型。該方法一定程度上解決上了上述問題，但在利用ANSYS的APDL語言生成模型時涉及大量命令流操作，影響模型構建的速度，如果遇到命令錯誤會導致模型無法完整生成。由此可見，如何有效、快速構建復雜三維地質模型進而開展數(shù)值計算是地應力反演的前提和重要研究課題[14-15]。

鑒于上述原因，本文提出了一種復雜三維地質模型的構建方法，該方法聯(lián)合3DMINE，GOCAD，SURFER，RHIRO，ANSYS以及FLAC3D以及MATLAB共7種軟件共同建模，減少了編程工作量、兼顧了ANSYS強大的前處理功能，彌補了大量命令流編寫、無法合理劃分網(wǎng)格等缺陷，實現(xiàn)了復雜三維地質體的快速構建。

1　建模思路

以山東鄆城煤礦為例，對整個建模思路進行介紹。為了客觀反演地應力的分布規(guī)律，綜合考慮地形、地質構造以及地層介質等因素，結合前期開展的地質勘探工作，最后確立了目標研究區(qū)域的范圍：以坐標點A(7 344.8，3 958.1)，B(9 324.8，3 958.1)，C(9 324.8，9 998.1)和D(7 344.8，9 998.1)圈定的長方形區(qū)域為計算區(qū)域，如圖1中所示。

圖1　幾何模型構建區(qū)域Fig.1　Construction area of geometric model

為了建立該的三維地質模型，使該模型更貼進實際的地形地貌，減少因數(shù)值模擬不精確帶來的地應力反演誤差，開展了礦區(qū)復雜三維地質模型精細化構建研究_建模原理基于“點-線-面-體”逐步構建的思路，在構建過程中使用了GOCAD，ANSYS，F(xiàn)LAC3D等多種軟件，具體過程及思路如圖2所示。

圖2　三維地質模型構建過程及思路Fig.2　Process of 3D numerical modeling

2　點陣數(shù)據(jù)和面的形成

常規(guī)數(shù)值模擬一般采用規(guī)則建模方法構建模型，這類建模方式由于軟件自身和建模方式不靈活從而無法完整表現(xiàn)地質體的起伏、斷層構造等信息。為了精細地描述復雜的地層結構，提出利用鉆孔信息生成巖層界面等高線，并“以點構線、以線構面、以面構體”的三級建模法構建模型，以形成更真實的復雜的地質體模型，提高數(shù)值模擬的精確度和可靠性。

2.1　等高線的生成

勘探鉆孔是構建巖層界面等值線的數(shù)據(jù)來源。鉆孔含有大量信息，通過鉆孔柱狀圖可以了解礦體、巖層界線以及斷層構造等隱藏于地下的信息。根據(jù)這些信息，可以生成相應的巖層界面等高線。利用CAD建立三維等高線操作較為麻煩，因此，選擇借助SURPAC或3DMINE等數(shù)字礦山三維軟件直接生成等高線，其操作步驟如下。

1)準備鉆孔原始數(shù)據(jù)。從柱狀圖中讀取相應的巖性信息表、鉆孔空間坐標，并導入到三維軟件之中。

2)三維鉆孔形態(tài)的生成。根據(jù)數(shù)據(jù)庫中存儲的鉆孔信息，生成三維形態(tài)的鉆孔立體圖，如圖3(a)所示。圖中以開孔位置為起點、終孔位置為終點連接成線，并按照鉆孔巖性表中的巖性特點分段標記。

3)三維等高線的生成。根據(jù)這些不同地層的不同標記，提取相應地層的點坐標，這些點坐標經(jīng)過分步取樣、插值即可生成等高線圖，如圖3(b)所示的煤層頂板等高線。

圖3　利用鉆孔信息生成等值線Fig.3　Contour map generation using borehole data

鉆孔數(shù)量的多少影響數(shù)值建模的精細程度，因此，可根據(jù)工程實際要求設計不同的鉆孔數(shù)量需求。如存在大量外部數(shù)據(jù)(如等高線)，可直接進行加載，導入GOCAD進行離散處理。

2.2　離散點陣的生成

由于前期鉆孔信息有限，所生成的等高線上僅存有少量離散的擬合點坐標數(shù)據(jù)，這些數(shù)據(jù)如果直接導入ANSYS或FLAC3D進行建模，一方面將造成模型偏差較大(嚴重失真)，另一方面利用少量點構建體操作性不強。因此，需要對等高線進行離散化和細分處理，生成相應的等高線擬合點數(shù)據(jù)。

經(jīng)過3DMINE處理后生成dxf文件，并導入到GOCAD中。GOCAD讀取dxf數(shù)據(jù)，提取出等高線中擬合點的數(shù)據(jù)。如圖4所示，為煤層頂板等高線的網(wǎng)格點陣圖，其中包括等高線和陣列點的空間坐標信息。

圖4　利用GOCAD生成網(wǎng)格點陣坐標(煤層頂板)Fig.4　Grid lattice coordinates using GOCAD

2.3　數(shù)據(jù)網(wǎng)格化處理

前期得到的GOCAD導出的等高線坐標點數(shù)據(jù)有限且較為離散，直接導入ANSYS構建實體模型較為復雜，因此，借助SURFER軟件對離散數(shù)據(jù)進行插值處理，使數(shù)據(jù)點均勻化，即網(wǎng)格化操作。

利用SURFER軟件“數(shù)據(jù)”選項讀取上一步獲得的離散等高線擬合點數(shù)據(jù)，選擇“克里金插值算法”對離散數(shù)據(jù)進行插值和均勻化處理，即得到等高線的點陣數(shù)據(jù)集，生成線框圖如圖5所示。

圖5　點陣數(shù)據(jù)構建的線框圖Fig.5　Wireframe generation using grid lattice coordinates

需要注意的是，由于地質模型在導入ANSYS中時會考慮模型大小問題，因此，需要對鉆孔的絕對坐標進行偏移，以保證模型的順利生成。坐標系偏移為：X偏移-20479 680，Y偏移-3911 100(這是后文進行反演時的重要參數(shù)，以此對比CAD圖紙上任一點與幾何模型坐標的關聯(lián)性)。

2.4　不規(guī)則面的生成

SURFER提供的點坐標數(shù)據(jù)無法直接導入ANSYS模型，借助命令流過程復雜。因此，為減少編程帶來的工作量，將SURFER生成的dat文件導入RHINO三維軟件，并利用RHINO提供的srfptgrid命令將數(shù)據(jù)點連線成面，如圖6所示，最后將該面文件轉換為ANSYS可讀的iegs格式文件(該格式可供ANSYS直接讀取)。如圖6所示，利用RHINO軟件生成的三維面模型，分別展示了該面模型的俯視圖、前視圖、右視圖和透視圖。

圖6　模型的生成與格式轉換Fig.6　Face model generation and data format conversion

3　模型的建立與數(shù)值計算

FLAC3D是采礦、巖土工程等領域的專業(yè)數(shù)值模擬軟件，在數(shù)值計算和分析方面有先天優(yōu)勢，但其建模功能受制于命令流式的手動建模方式，對于構建復雜模型過程復雜。ANSYS有限元程序具有強大的前處理功能，因此，借助ANSYS軟件建立三維實體模型、劃分網(wǎng)格單元，然后導入FLAC3D進行數(shù)值計算，將大大簡化建模過程，提高數(shù)值模擬的效率和精度。

3.1　幾何模型構建

ANSYS提供了2種建模方式，一是創(chuàng)建或導入實體模型，然后進行網(wǎng)格劃分，并生成有限元模型；另一種是利用單元和節(jié)點編號，依靠ANSYS提供的命令流進行生成操作。本文選擇前者，其建模過程可簡述如下。

1)將前面處理好的關鍵巖層分界面iges文件導入WORKBENCH的Geometry之中，形成如圖7(a)線框圖，(b)面域圖所示的面文件。

2)將相鄰兩面選取，點擊Skin命令蒙皮成體，“面面構體”生成體單元，如圖7(c)所示，完成后組成新的Part部分。

3)轉入到“Model”模式，選擇合理參數(shù)，對模型進行網(wǎng)格劃分。

4)將Model中劃分網(wǎng)格后的模型導入到Mechanical APDL模式。

5)然后利用河海大學鄭文棠博士編寫的ANSYS_TO_FLAC3D接口程序將ANSYS模型信息導出，另存為節(jié)點(Node)和單元(Element)2個文件，并經(jīng)過轉換文件轉換為*.flac3d文件。FLAC3D通過Import Grid即可直接讀取該文件，得到如圖8所示的三維地質模型。

圖7　ANSYS中由面構體形成的三維地質模型Fig.7　3D model generation using the face data

3.2　網(wǎng)格劃分

礦山工程巖體應力具有三維特征，然而三維數(shù)值計算及成圖復雜，對計算軟件要求高，再加上反演區(qū)域范圍大，劃分單元多，給計算帶來難度。為此，模型大小及網(wǎng)格劃分尤為關鍵。

ANSYS提供了強大的網(wǎng)格劃分功能，相比FLAC3D更為智能和便利。綜合考慮有限元數(shù)值模擬的速度和構建效果，網(wǎng)格劃分時選定的網(wǎng)格大小為80 m，最小網(wǎng)格為0.5 m，采用六面體網(wǎng)格劃分法將計算模型進行劃分，共劃分得到114 240個單元，488 179個節(jié)點，網(wǎng)格劃分完畢后的效果如圖8所示。

圖8　網(wǎng)格劃分結果Fig.8　Mesh partition

3.3　FLAC3D中模型建立與數(shù)值計算

在ANSYS WORKBENCH中建立好三維地質模型后，對實體模型進行單元劃分，將模型導入到ANSYS APDL中，然后利用河海大學鄭文棠博士編寫的ANSYS_TO_FLAC3D接口程序將ANSYS模型信息導出，另存為節(jié)點(Node)和單元(Element)兩個文件，并經(jīng)過轉換文件轉換為*.flac3d文件。FLAC3D通過Import Grid即可直接讀取該文件，得到如圖9所示的三維地質模型。

圖9　有限元數(shù)值模型的建立與導入Fig.9　Establishment and loading of finite element model

根據(jù)實際情況對不同巖層分別賦予材料屬性，并進行分組處理；然后依據(jù)前期分析確立相應的邊界條件和計算參數(shù)，最后導入模型進行數(shù)值計算。

4　地應力反演結果展示與分析

以山東鄆城煤礦為例，根據(jù)該礦前期測試的4個測點的地應力數(shù)據(jù)，通過建立三維數(shù)值計算模型，利用先進的數(shù)值模擬軟件(FLAC3D)和應力反分析技術，反演和重構鄆城煤礦一采區(qū)初始地應力場。最終獲得的鄆城煤礦初始地應力場回歸模型，回歸模型的相關性系數(shù)為88.26%，表明擬合效果較好。

圖10所示為利用回歸系數(shù)重構的初始地應力場SZZ的分布情況，取其中3組巖層上覆巖層、煤層和底部巖層展示。

圖10　煤巖層SZZ應力分布規(guī)律Fig.10　SZZ stress distribution of coal seam and rock strata

圖11　-800 m最大主應力SIG1分布云圖Fig.11　Nephogram of maximum principal stress on -800 m level

為了更清晰展示鄆城煤礦一采區(qū)區(qū)域地應力概況，將反演重構的煤層地應力場投影到CAD圖上，得到如圖11所示的最大主應力SIG1分布云圖?？梢钥闯?，最大主應力與埋深和構造密切相關。圖中深色區(qū)域為主應力較大區(qū)域，這些區(qū)域正好位于褶曲構造位置，由此可以推斷，褶曲是引起區(qū)域內(nèi)巖體應力高度集中的主因；同時，區(qū)域內(nèi)的水平應力大于垂直應力。

利用上述建模方法，實現(xiàn)了對復雜礦山地質模型的快速構建，這為準確計算和展示地應力場分布規(guī)律提供了依據(jù)，為區(qū)域地應力場反演奠定了基礎，同時也為礦井深部開采和采掘規(guī)劃提供了依據(jù)。

5　結論

1)利用3DMINE，SUFER，ANSYS以及FLAC3D等軟件聯(lián)合建立大型精細三維有限元模型的方法，通過“點-線-面-體”的建模思路，實現(xiàn)了復雜三維地質模型的構建。

2)以山東鄆城煤礦為例，通過構建了該礦的三維地質模型，結合現(xiàn)場實測測點地應力值，反演出該礦的地應力場分布，并以應力云圖的方式進行直觀展示，效果良好。

3)實際工程應用表明，該方法可實現(xiàn)大尺度、復雜三維地質模型的構建，有效解決了復雜模型構建效率低、網(wǎng)格劃分難度大等問題，通過數(shù)據(jù)的逐步對接，減少了中間過程的人為干擾和誤差，一定程度縮短了建模時間、提高了數(shù)值計算精度。

[1]康紅普，林健，張曉. 深部礦井地應力測量方法研究與應用[J]. 巖石力學與工程學報，2007，26(5)：929-933.

KANG Hongpu, LIN Jian, ZHANG Xiao. Research and application of in-situ stress measurement in deep mines[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(5): 929-933.

[2]楊德林. 巖土工程問題的反演理論與工程實踐[M]. 北京：科學出版社，1996.

[3]龔曉南．對巖土工程數(shù)值分析的幾點思考[J]．巖土力學，2011，32(2)：321-325.

GONG Xiaonan. Reflections on numerical analysis of geotechnical engineering[J]. Rock and Soil Mechanics, 2011, 32(2): 321-325.

[4]朱志潔，張宏偉，韓軍，等. 不同應力條件下巷道穩(wěn)定性研究[J]. 中國安全生產(chǎn)科學技術，2015，11(11)：11-16.

ZHU Zhijie, ZHANG Hongwei, HAN Jun, et al. Research on stability of roadway under different conditions of stress field[J]. Journal of Safety Science and Technology, 2015, 11(11): 11-16.

[5]劉寧，朱維申，辛小麗. 雙江口水電站初始地應力場反演回歸分析[J]. 山東大學學報(工學版)，2008，38(6)：121-126.

LIU Ning, ZHU Weishen, XIN Xiaoli. Back regression analysis on the initial geostress field of the Shuangjiangkou hydropower station[J]. Journal of Shandong university: Engineering Science, 2008, 38(6): 121-126.

[6]鄭文棠，徐衛(wèi)亞，童富果，等．復雜邊坡三維地質可視化和數(shù)值模型構建[J]．巖石力學與工程學報，2007，26(6)：1633-1644．

ZHENG Wentang, XU Weiya, TONG Fuguo, et al. 3D geological visualization and numerical modeling of complicated slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1633-1644.

[7]裴啟濤，李海波，劉亞群，等. 復雜地質條件下壩區(qū)初始地應力場二次反演分析[J]. 巖石力學與工程學報，2014，33(S1)：2779-2785.

PEI Qitao, LI Haibo, LIU Yaqun, et al. Two-stage back analysis of initial geostress field of dam areas under complex geological conditions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S1): 2779-2785.

[8]崔芳鵬，胡瑞林，劉照連，等．基于Surfer平臺的FLAC3D復雜三維地質建模研究[J]．工程地質學報，2008(5)：699-702.

CUI Fangpeng, HU Ruilin, LIU Zhaolian, et al. Surfer software platform based on complex three dimensional geological digital models for pre-processing of FLAC3D[J]. Journal of Engineering Geology, 2008(5): 699-702.

[9]廖秋林，曾錢幫，劉彤，等．基于ANSYS平臺復雜地質體FLAC3D模型的自動生成[J]．巖石力學與工程學報，2005，24(6)：1010-1013.

LIAO Qiulin, ZENG Qianbang, LIU Tong, et al. Automatic model generation of complex geologic body with FLAC3D based on ansys platform[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(6): 1010-1013.

[10]向鵬，紀洪廣，鄒靜，等. 復雜礦山地質形體的三維數(shù)值網(wǎng)格劃分自動化平臺[J]. 巖土力學，2015，36(4)：1211-1216.

XIANG Peng, JI Hongguang, ZOU Jing, et al. An automation platform for generating the 3D numerical grids of mining fields with complex geomorphology[J]. Rock and Soil Mechanics, 2015, 36(4): 1211-1216.

[11]趙文，段魁，王彥東，等. 基于MATLAB、SURFER平臺的FLAC3D復雜三維地質建模研究[J]. 工程地質學報，2016，24(s1)：208-213.

ZHAO Wen, DUAN Kui, WANG Yandong, et al. Study on complex three-dimensional geological modeling of Flac3D based on Matlab and Surfer platform[J]. Journal of Engineering Geology, 2016, 24(s1): 208-213.

[12]袁海平，王金安，蔡美峰. 復雜地應力場幾何跨尺度反演與重構方法研究[J]. 采礦與安全工程學報，2011，28(4)：589-595.

YUAN Haiping, WANG Jinan, CAI Meifeng. Geometric cross-scale back-analysis and reconstruction menthod for complicated geo-stress field[J]. Journal of Mining & Safety Engineering, 2011, 28(4): 589-595.

[13]李翔，王金安，張少杰. 復雜地質體三維數(shù)值建模方法研究[J]. 西安科技大學學報，2012，32(6)：676-681.

LI Xiang, WANG Jinan, ZHANG Shaojie. 3D modeling method of complicated geological body[J]. Journal of Xi’an University of Science and Technology, 2012, 32(6): 676-681.

[14]王金安，李飛. 復雜地應力場反演優(yōu)化算法及研究新進展[J]. 中國礦業(yè)大學學報，2015，44(2)：189-205.

WANG Jinan, LI Fei. Review of inverse optimal algorithm of in-situ stress and new achievement[J]. Journal of China University of Mining & Technology, 2015, 44(2): 189-205.

[15]趙自強，王學潮，隨裕紅. 南水北調西線工程深埋隧洞巖爆與地應力研究[J]. 華北水利水電學院學報，2002，23(1)：48-50.

ZHAO Ziqiang, WANG Xuechao, SUI Yuhong. The deep tunnel rock blasting and crustal stress research of west route project of water transfer from south to north[J]. Journal of North China Institute of Water Conservancy and Hyadroelectric Power, 2002, 23(1): 48-50.