大孤山露天礦邊坡億級(jí)自由度建模與基于RFPA3D的數(shù)值模擬

李斕堃 朱萬成 代 風(fēng) 劉溪鴿 鄧文學(xué)

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

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

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

在數(shù)值模擬計(jì)算中，數(shù)值結(jié)果依賴于網(wǎng)格尺寸，如果模型中相互接觸的網(wǎng)格尺寸比值過大，則容易造成數(shù)值模擬結(jié)果的失真。因此，諸多研究和實(shí)踐表明，數(shù)值計(jì)算模型中單元體積的變化率應(yīng)該控制在30%以內(nèi)［10］。此外，網(wǎng)格的尺度還應(yīng)該小于材料變形的最小特征尺寸，以滿足細(xì)致刻畫巖體變形和損傷演變過程的要求?；谏鲜鰞牲c(diǎn)，對(duì)于含有多種巖體的大型復(fù)雜巖石工程而言，既要滿足網(wǎng)格尺寸小于巖石變形的最小特征尺寸，又要滿足單元體積變化率控制在30%的要求，這就需要建立尺寸足夠小、數(shù)量足夠多的網(wǎng)格和單元，其自由度往往達(dá)到千萬級(jí)甚至億級(jí)。

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

本研究提出了新的結(jié)構(gòu)化六面體網(wǎng)格化方法，使用東北大學(xué)AMAX PSC-HC2X GPU服務(wù)器，建立了大孤山露天礦億級(jí)自由度三維模型，并基于RFPA3D-Parallel軟件平臺(tái)對(duì)該模型進(jìn)行了數(shù)值分析，以期為邊坡穩(wěn)定性評(píng)價(jià)及災(zāi)害防控提供參考。

1 RFPA3D-Parallel精細(xì)化數(shù)值模擬方法

1.1 RFPA3D-Parallel簡(jiǎn)介

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

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

巖石的損傷與巖石的應(yīng)力狀態(tài)密切相關(guān)。RFPA3D中采用簡(jiǎn)單的彈性損傷本構(gòu)模型，在達(dá)到破壞準(zhǔn)則之前，單元保持線彈性的力學(xué)本質(zhì)。當(dāng)巖石的應(yīng)力狀態(tài)滿足最大拉應(yīng)力準(zhǔn)則或摩爾—庫侖準(zhǔn)則時(shí)，巖石就開始發(fā)生破壞［27］，其力學(xué)模型為

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

1.2 網(wǎng)格劃分與模型構(gòu)建

為快速、自動(dòng)化實(shí)現(xiàn)復(fù)雜巖體六面體網(wǎng)格劃分，并提供RFPA3D-Parallel可直接讀取的網(wǎng)格文件，本項(xiàng)目開發(fā)了一種復(fù)雜巖體億級(jí)自由度結(jié)構(gòu)化六面體網(wǎng)格參數(shù)化建模方法?；趦|級(jí)自由度的結(jié)構(gòu)化六面體網(wǎng)格參數(shù)化建模方法包含了背景網(wǎng)格劃分、四面體網(wǎng)格數(shù)量估算、六面體網(wǎng)格更新、空屬性六面體網(wǎng)格檢索和剔除、網(wǎng)格文件導(dǎo)入等過程?；赗FPA3DParallel數(shù)值模擬軟件，針對(duì)大規(guī)模復(fù)雜巖體精細(xì)化數(shù)值計(jì)算的技術(shù)路線如圖1所示。

（1）非結(jié)構(gòu)化背景網(wǎng)格劃分。首先對(duì)大型不規(guī)則結(jié)構(gòu)的三維復(fù)雜地質(zhì)體模型進(jìn)行萬級(jí)自由度的背景網(wǎng)格劃分，并賦予編號(hào)IDt，t從1開始。該網(wǎng)格劃分過程采用非結(jié)構(gòu)化網(wǎng)格，它針對(duì)復(fù)雜地質(zhì)體的不同位置采用自適應(yīng)過渡的四面體網(wǎng)格自動(dòng)劃分方法。

（2）結(jié)構(gòu)化六面體網(wǎng)格建立。根據(jù)非結(jié)構(gòu)化的四面體背景網(wǎng)格信息，確定x、y、z3個(gè)方向上的坐標(biāo)最大值和最小值xmax、ymax、zmax、xmin、ymin、zmin。從而確定包含整個(gè)復(fù)雜巖體最小長(zhǎng)方體區(qū)域的長(zhǎng)寬高X、Y和Z，根據(jù)目標(biāo)自由度的數(shù)量m，采用立方體網(wǎng)格的情況下，假設(shè)z方向有n個(gè)網(wǎng)格，則

由式（4）確定n近似值，取整后可得到x、y、z方向上網(wǎng)格數(shù)量，分別為Xnum、Ynum、Znum。故可對(duì)該長(zhǎng)方體區(qū)域進(jìn)行參數(shù)化六面體網(wǎng)分割，獲得Xnum×Ynum×Znum個(gè)結(jié)構(gòu)化六面體網(wǎng)格，并賦予編號(hào)IDs，其網(wǎng)格節(jié)點(diǎn)編號(hào)s形式為（x、y、z），按照x、y、z坐標(biāo)從小至大排序。

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

（4）空屬性六面體網(wǎng)格剔除。由于背景網(wǎng)格在巖體界面上的非連續(xù)性，以及巖體開挖運(yùn)算的存在，極易導(dǎo)致在巖體界面上出現(xiàn)零星的空屬性六面體，甚至在巖體的開挖區(qū)出現(xiàn)成團(tuán)聚集的空屬性六面體集合。此時(shí)需要人為地將上述空屬性六面體網(wǎng)格剔除或賦值，并對(duì)全部網(wǎng)格進(jìn)行重新排序，更新六面體單元的拓?fù)潢P(guān)系。

（5）億級(jí)自由度網(wǎng)格文件導(dǎo)入RFPA3D-Parallel。編寫RFPA3D-Parallel的輸出接口，導(dǎo)出節(jié)點(diǎn)坐標(biāo)、邊界約束、單元節(jié)點(diǎn)編號(hào)、節(jié)點(diǎn)力4個(gè)文本文件，將4個(gè)結(jié)果文件拷貝到RFPA文件所在目錄，新建工程并保存后即可導(dǎo)入。

2 在大孤山露天礦的應(yīng)用

2.1 工程地質(zhì)概況及數(shù)值建模

大孤山鐵礦自1916年開始開采以來逐漸成為亞洲典型的深凹露天生產(chǎn)基地。礦場(chǎng)封閉圈標(biāo)高+90 m，沿走向長(zhǎng)1 700 m，寬1 500 m，現(xiàn)已開采到-330 m標(biāo)高，開采深度已達(dá)420 m，設(shè)計(jì)開采深度為528 m。

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

為了構(gòu)建三維數(shù)值模型，進(jìn)行大孤山礦區(qū)三維數(shù)值模擬分析，利用現(xiàn)有獲取的地質(zhì)資料對(duì)礦區(qū)進(jìn)行了詳細(xì)三維地質(zhì)模型構(gòu)建。采用Rhino軟件構(gòu)建斷層、巖性分界面和開采境界的三維NUBERS曲面，建立模型區(qū)域三維長(zhǎng)方體實(shí)體模型（2 000 m×2 000 m×890 m），基于布爾運(yùn)算分割三維實(shí)體模型，最終形成圖3所示的露天礦復(fù)雜地質(zhì)體三維模型，共計(jì)19種復(fù)雜巖體結(jié)構(gòu)，27個(gè)露天境界設(shè)計(jì)臺(tái)階，圍巖中包含2條斷層，多個(gè)礦體以及包括斷層破碎帶在內(nèi)的9種不同巖性。

2.2 數(shù)值模擬

對(duì)大孤山三維地質(zhì)體模型進(jìn)行萬級(jí)自由度的四面體網(wǎng)格劃分，因其不直接參與有限元計(jì)算，無需保證各地質(zhì)體交界面上網(wǎng)格的嚴(yán)格連續(xù)性。因此，全過程可自動(dòng)化實(shí)現(xiàn)，且較少的網(wǎng)格即可準(zhǔn)確地表達(dá)出復(fù)雜地質(zhì)體的幾何形態(tài)，如圖4所示。經(jīng)統(tǒng)計(jì)，網(wǎng)格數(shù)量為119 077個(gè)，節(jié)點(diǎn)數(shù)量為32 515個(gè)。

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

將剔除空屬性的六面體群的結(jié)構(gòu)化網(wǎng)格后，對(duì)節(jié)點(diǎn)重新編號(hào)，并生成網(wǎng)格計(jì)算文件，導(dǎo)入RFPA3DParallel軟件中，如圖6所示，單元數(shù)量為33 163 200個(gè)，節(jié)點(diǎn)33 498 549個(gè)，自由度達(dá)100 495 647。

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

2.3 數(shù)值模擬結(jié)果分析

通過數(shù)值求解，得到大孤山鐵礦損傷場(chǎng)和彈性模量場(chǎng)如圖7所示。首先進(jìn)行重力平衡計(jì)算，隨后采用離心加載法［28］逐步增加重力加速度，離心加載系數(shù)取值為0.4。計(jì)算結(jié)果顯示，F(xiàn)14斷層帶與F15斷層帶之間的楔形體礦石條帶出現(xiàn)了滑移現(xiàn)象，如圖7所示，數(shù)值模擬結(jié)果與大孤山鐵礦觀測(cè)到的小范圍滑動(dòng)破壞的位置、規(guī)?；鞠喾?/p>

圖7表明大孤山鐵礦的西北幫為滑坡災(zāi)害發(fā)生的高風(fēng)險(xiǎn)區(qū)域，并顯示了該區(qū)域楔形體滑移形成過程中損傷演化和裂紋擴(kuò)展等過程，據(jù)此可將整個(gè)楔形體滑移的孕育和發(fā)生過程分為多個(gè)階段。

Step1～Step6，模型在自重應(yīng)力下，整體逐漸發(fā)生沉降和壓實(shí)，損傷場(chǎng)和彈性模量場(chǎng)幾乎沒有任何變化。

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

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

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

3 結(jié)論

針對(duì)礦山三維復(fù)雜巖體大尺度建模與求解計(jì)算難點(diǎn)，本研究在進(jìn)行大孤山露天礦億級(jí)自由度建模的基礎(chǔ)上，利用RFPA3D-Parallel軟件平臺(tái)研究了大孤山露天礦西北幫滑坡災(zāi)害形成過程，得到如下主要結(jié)論：

（1）實(shí)現(xiàn)了露天礦工程尺度復(fù)雜巖體的結(jié)構(gòu)化六面體網(wǎng)格精細(xì)構(gòu)建，且基于并行計(jì)算平臺(tái)實(shí)現(xiàn)了大孤山露天礦邊坡穩(wěn)定性的億級(jí)自由度的計(jì)算分析。

（2）數(shù)值模擬結(jié)果表明，大孤山露天礦西北幫的F14斷層帶與F15斷層帶之間楔形礦石條帶出現(xiàn)局部損傷區(qū)，該結(jié)果與現(xiàn)場(chǎng)實(shí)際觀測(cè)情況完全吻合，即該楔形體有可能引起滑坡災(zāi)害，應(yīng)采取加固、削坡等方式進(jìn)行治理。

（3）大孤山露天礦邊坡楔形體的損傷具有自上而下的發(fā)展演化特征，即損傷首先發(fā)生在上部臺(tái)階，隨著裂隙向下逐漸擴(kuò)展、貫通，將最終導(dǎo)致楔形體沿?cái)鄬影l(fā)生錯(cuò)動(dòng)、滑移，這一認(rèn)識(shí)對(duì)現(xiàn)場(chǎng)滑坡監(jiān)測(cè)預(yù)警及防控治理具有重要的指導(dǎo)意義。