大孤山露天礦邊坡億級自由度建模與基于RFPA3D的數值模擬

李斕堃 朱萬成 代 風 劉溪鴿 鄧文學

（東北大學資源與土木工程學院，遼寧 沈陽 110819）

近年來，數值模擬被廣泛應用于采礦工程中的巖體穩定性分析。大型巖土工程計算區域常達千米級，并且地質構造復雜、巖石類型多樣以及涉及到多場耦合等問題。Park D等［1］建立了長達1 000 m的隧道模型，并得到隧道沿軸線方向的變形響應特征；尚振華等［2］采用FLAC3D軟件對尺寸為1 000 m×1 600 m×660 m的某鎢礦模型進行模擬，將采空區穩定性問題量化；劉春等［3］建立了尺寸為2 800 m×1 500 m×1 600 m的茂縣新磨村滑坡模型，模擬了新磨村滑坡啟動、高速下滑和堆積的全過程；張成良等［4］采用三維有限元模擬了3種開采方案對地表公路的影響，得到當開采高度超過1 300 m水平時，公路上部圍巖將處于不穩定狀態；漆祖芳等［5］建立630 m×1 850 m×900 m的有限元模型，模擬了大崗山水電站壩肩邊坡開挖支護過程；馬克等［6］用RFPA3D對大崗山水電站壩址區右岸實際微震監測范圍（690 m×400 m×1 600 m）等比例的三維模型進行模擬，對抗剪洞加固前、后的邊坡穩定性進行分析。

考慮到計算成本及技術方面的問題，前人對計算模型進行了大量簡化，且受到計算模型網格數量受限，致使模型網格尺寸較大，導致了數值模型難以反映實際工程的現場真實情況［7］。因此，對于大規模、復雜巖石工程的數值模擬，大多數情況下會將關鍵位置網格和非關鍵位置網格差異化對待，即關鍵位置網格剖分需要細化，而非關鍵位置可采用過渡單元，將單元尺寸從核心研究區域向邊界區域逐漸粗化［8-9］。網格單元尺寸的差異化將客觀地對巖石工程內部應力的傳播、損傷演化及裂紋擴展等路徑造成不可忽視的影響。因此，在不考慮計算機性能和計算成本的前提下，數值模擬計算的期望是將計算單元或網格合理地細化，以期更為準確地表征圍巖特征。

在數值模擬計算中，數值結果依賴于網格尺寸，如果模型中相互接觸的網格尺寸比值過大，則容易造成數值模擬結果的失真。因此，諸多研究和實踐表明，數值計算模型中單元體積的變化率應該控制在30%以內［10］。此外，網格的尺度還應該小于材料變形的最小特征尺寸，以滿足細致刻畫巖體變形和損傷演變過程的要求。基于上述兩點，對于含有多種巖體的大型復雜巖石工程而言，既要滿足網格尺寸小于巖石變形的最小特征尺寸，又要滿足單元體積變化率控制在30%的要求，這就需要建立尺寸足夠小、數量足夠多的網格和單元，其自由度往往達到千萬級甚至億級。

隨著大型工程對精細數值分析要求的不斷增加，超過千萬自由度的數值分析不斷增多［11］，然而，受到計算機性能，如CPU浮點計算能力、GPU圖像處理能力以及內存開銷等硬件因素，采用常規數值模擬軟件計算一些大的模型需要幾天甚至幾個月的時間［12］，這樣就極大提高了工程分析的時間成本，也使得分析人員不得不減小計算規模。為此，一些研究人員利用高性能并行計算機對復雜巖土工程問題進行模擬。范宣華等［13］構建了大規模模態分析并行計算體系，并驗證3種算法構建的并行求解體系均可在1 h內求解千萬級自由度量的大規模模態分析問題。張友良等［14］采用網格加密方法生成5 000萬自由度的有限元模型，利用對偶原始有限元撕裂內聯法（FETI-DP）進行求解。K.Garatani等［15］采用GeoFEM系統對簡單形狀進行線彈性分析，該計算自由度高達一億。劉耀儒等［16-17］考慮到耦合分析時步多、計算量大的問題，采用基于element-by-element策略的有限元并行計算方法進行數值模擬。針對巖石破裂及大型巖體失穩的求解問題，張永彬等［18］結合現代有限元方法和數值計算方法，在消息傳遞并行環境下，利用區域分解和主從編程模式，建立了巖石破裂過程分析RFPA3D-Parallel并行分析系統。

本研究提出了新的結構化六面體網格化方法，使用東北大學AMAX PSC-HC2X GPU服務器，建立了大孤山露天礦億級自由度三維模型，并基于RFPA3D-Parallel軟件平臺對該模型進行了數值分析，以期為邊坡穩定性評價及災害防控提供參考。

1 RFPA3D-Parallel精細化數值模擬方法

1.1 RFPA3D-Parallel簡介

RFPA3D-Parallel并行分析系統采用分布存儲稀疏線性迭代并行求解方法，在Linux機群上實現應力分析模塊中有限元計算的并行處理，通過Windows和Linux協調處理策略，有效地把原有的前后處理功能和機群系統強大的計算能力結合起來，建立巖石破裂過程分析平臺［19］。經過后期實驗驗證，RFPA3DParallel在確保高效、準確、強大的計算能力的同時，實現了巖石介質破壞過程的準確還原［20-21］。

為了在數值模擬中表達天然巖石的非均質特征，假設巖石由許多細觀單元組成，其力學性質服從Weibull分布［22］。利用Weibull分布確定巖石參數的方法已被大量應用于數值模擬，并取得了很好的結果［23-26］。

巖石的損傷與巖石的應力狀態密切相關。RFPA3D中采用簡單的彈性損傷本構模型，在達到破壞準則之前，單元保持線彈性的力學本質。當巖石的應力狀態滿足最大拉應力準則或摩爾—庫侖準則時，巖石就開始發生破壞［27］，其力學模型為

式（1）為最大拉應力準則；式（2）為經典的摩爾庫倫準則；式（3）中ω為損傷變量，E、E0分別為損傷單元和未損傷單元的彈性模量。

1.2 網格劃分與模型構建

為快速、自動化實現復雜巖體六面體網格劃分，并提供RFPA3D-Parallel可直接讀取的網格文件，本項目開發了一種復雜巖體億級自由度結構化六面體網格參數化建模方法。基于億級自由度的結構化六面體網格參數化建模方法包含了背景網格劃分、四面體網格數量估算、六面體網格更新、空屬性六面體網格檢索和剔除、網格文件導入等過程。基于RFPA3DParallel數值模擬軟件，針對大規模復雜巖體精細化數值計算的技術路線如圖1所示。

（1）非結構化背景網格劃分。首先對大型不規則結構的三維復雜地質體模型進行萬級自由度的背景網格劃分，并賦予編號IDt，t從1開始。該網格劃分過程采用非結構化網格，它針對復雜地質體的不同位置采用自適應過渡的四面體網格自動劃分方法。

（2）結構化六面體網格建立。根據非結構化的四面體背景網格信息，確定x、y、z3個方向上的坐標最大值和最小值xmax、ymax、zmax、xmin、ymin、zmin。從而確定包含整個復雜巖體最小長方體區域的長寬高X、Y和Z，根據目標自由度的數量m，采用立方體網格的情況下，假設z方向有n個網格，則

由式（4）確定n近似值，取整后可得到x、y、z方向上網格數量，分別為Xnum、Ynum、Znum。故可對該長方體區域進行參數化六面體網分割，獲得Xnum×Ynum×Znum個結構化六面體網格，并賦予編號IDs，其網格節點編號s形式為（x、y、z），按照x、y、z坐標從小至大排序。

（3）空屬性六面體網格檢索。每個六面體網格質心坐標與背景網格中對應的四面體（如圖2），通過其在背景網格中的網格編號IDt，檢索其所從屬的復雜巖體網格編號IDs。判斷空間中一點是否位于四面體內，可計算空間點與四面體任意3個頂點組成的4個小四面體體積之和，并比較其與該四面體體積的大小以作判斷：當兩者相等時，空間點位于四面體內或者四面體邊界上；當兩者不等時，空間點位于四面體外。將位于所有四面體外的六面體定義為空屬性六面體網格，對于上億級自由度的網格，該過程可采用并行計算加速六面體網格的復雜巖體屬性檢索過程。

（4）空屬性六面體網格剔除。由于背景網格在巖體界面上的非連續性，以及巖體開挖運算的存在，極易導致在巖體界面上出現零星的空屬性六面體，甚至在巖體的開挖區出現成團聚集的空屬性六面體集合。此時需要人為地將上述空屬性六面體網格剔除或賦值，并對全部網格進行重新排序，更新六面體單元的拓撲關系。

（5）億級自由度網格文件導入RFPA3D-Parallel。編寫RFPA3D-Parallel的輸出接口，導出節點坐標、邊界約束、單元節點編號、節點力4個文本文件，將4個結果文件拷貝到RFPA文件所在目錄，新建工程并保存后即可導入。

2 在大孤山露天礦的應用

2.1 工程地質概況及數值建模

大孤山鐵礦自1916年開始開采以來逐漸成為亞洲典型的深凹露天生產基地。礦場封閉圈標高+90 m，沿走向長1 700 m，寬1 500 m，現已開采到-330 m標高，開采深度已達420 m，設計開采深度為528 m。

大孤山露天礦地層主要由混合巖巖體、太古代花崗巖巖體、花崗斑巖巖體、綠泥石英片巖巖體、千枚巖巖體、玢巖巖體、小孤山礦體、大孤山礦體和低品位礦石條帶礦體等組成，各地層有明顯的巖性分界面。

為了構建三維數值模型，進行大孤山礦區三維數值模擬分析，利用現有獲取的地質資料對礦區進行了詳細三維地質模型構建。采用Rhino軟件構建斷層、巖性分界面和開采境界的三維NUBERS曲面，建立模型區域三維長方體實體模型（2 000 m×2 000 m×890 m），基于布爾運算分割三維實體模型，最終形成圖3所示的露天礦復雜地質體三維模型，共計19種復雜巖體結構，27個露天境界設計臺階，圍巖中包含2條斷層，多個礦體以及包括斷層破碎帶在內的9種不同巖性。

2.2 數值模擬

對大孤山三維地質體模型進行萬級自由度的四面體網格劃分，因其不直接參與有限元計算，無需保證各地質體交界面上網格的嚴格連續性。因此，全過程可自動化實現，且較少的網格即可準確地表達出復雜地質體的幾何形態，如圖4所示。經統計，網格數量為119 077個，節點數量為32 515個。

根據1億自由度的網格目標，采用式（4）可確定x、y和z方向上各有420、420和188個網格，共計33 498 549個節點，此為達到1億自由度所需的最少節點，建立結構化六面體網格后，采用并行計算加速六面體網格的復雜巖體屬性檢索過程，對于空屬性的六面體集合，根據有限元的計算需要剔除或者賦值開挖屬性，空屬性網格的空間分布如圖5所示。

將剔除空屬性的六面體群的結構化網格后，對節點重新編號，并生成網格計算文件，導入RFPA3DParallel軟件中，如圖6所示，單元數量為33 163 200個，節點33 498 549個，自由度達100 495 647。

大孤山鐵礦礦體及其圍巖涉及的全部巖石材料的物理力學性質見表1。

2.3 數值模擬結果分析

通過數值求解，得到大孤山鐵礦損傷場和彈性模量場如圖7所示。首先進行重力平衡計算，隨后采用離心加載法［28］逐步增加重力加速度，離心加載系數取值為0.4。計算結果顯示，F14斷層帶與F15斷層帶之間的楔形體礦石條帶出現了滑移現象，如圖7所示，數值模擬結果與大孤山鐵礦觀測到的小范圍滑動破壞的位置、規模基本相符。

圖7表明大孤山鐵礦的西北幫為滑坡災害發生的高風險區域，并顯示了該區域楔形體滑移形成過程中損傷演化和裂紋擴展等過程，據此可將整個楔形體滑移的孕育和發生過程分為多個階段。

Step1～Step6，模型在自重應力下，整體逐漸發生沉降和壓實，損傷場和彈性模量場幾乎沒有任何變化。

Step7～Step13，因模型各部分巖體強度存在差異，以及各部分巖體之間結構面的存在。模型各部分巖體在自重應力下相互擠壓、緩慢變形和錯動，這造成了坡體斷層面處首先出現應力集中，繼而出現縱向剪切損傷和裂紋，但由于其坡腳下部臺階較厚大，抑制了邊坡的整體變形。

Step14～Step20，楔形體上部已有裂隙在巖體自重應力作用下持續發育、擴展和貫通，部分裂隙貫穿巖體，坡體表層出現滑塌后，坡體失穩加快。礦石條帶沿著F14斷層帶與F15斷層帶組成的楔形滑面發生整體錯動，局部巖塊滾落，并且在滾落過程中迅速解體，滑至坡腳時已碎裂成塊。

大孤山鐵礦邊坡的西北幫的變形及局部滑塌主要是由F14斷層帶與F15斷層帶組成的楔形造成的。巖體揭露后，坡體暫時沒有出現明顯的損傷和裂紋擴展行為；在自重應力的持續作用下，礦石條帶與F14斷層、F15斷層相互擠壓、錯動，導致應力集中，繼而在坡體上部表面出現損傷區和裂紋；隨著損傷區和裂紋的繼續發育、擴展和貫通，楔形體在自重應力作用下發生了整體滑移。上述現象表明，大孤山露天礦西北幫F14斷層與F15斷層之間的楔形礦石條帶存在較大滑坡風險。

3 結論

針對礦山三維復雜巖體大尺度建模與求解計算難點，本研究在進行大孤山露天礦億級自由度建模的基礎上，利用RFPA3D-Parallel軟件平臺研究了大孤山露天礦西北幫滑坡災害形成過程，得到如下主要結論：

（1）實現了露天礦工程尺度復雜巖體的結構化六面體網格精細構建，且基于并行計算平臺實現了大孤山露天礦邊坡穩定性的億級自由度的計算分析。

（2）數值模擬結果表明，大孤山露天礦西北幫的F14斷層帶與F15斷層帶之間楔形礦石條帶出現局部損傷區，該結果與現場實際觀測情況完全吻合，即該楔形體有可能引起滑坡災害，應采取加固、削坡等方式進行治理。

（3）大孤山露天礦邊坡楔形體的損傷具有自上而下的發展演化特征，即損傷首先發生在上部臺階，隨著裂隙向下逐漸擴展、貫通，將最終導致楔形體沿斷層發生錯動、滑移，這一認識對現場滑坡監測預警及防控治理具有重要的指導意義。