露采高邊坡模型初始地應(yīng)力場快速生成方法研究*

胡 斌，方慶紅，李 京，祝 鑫，楊 為

(1.武漢科技大學(xué) 資源與環(huán)境工程學(xué)院，湖北 武漢 430081；2.冶金礦產(chǎn)資源高效利用與造塊湖北省重點(diǎn)實(shí)驗室，湖北 武漢 430081)

0 引言

在采礦工程領(lǐng)域中，初始地應(yīng)力場既是礦山工程設(shè)計施工重要因素之一，也是邊坡穩(wěn)定性模擬研究必須考慮因素[1-3]。在礦山高邊坡數(shù)值模擬試驗中，為防止初始地應(yīng)力場畸變，必須建立取至負(fù)地形的整體大模型，但大模型存在模型龐大、結(jié)構(gòu)復(fù)雜、劃分單元和節(jié)點(diǎn)數(shù)繁多、計算量巨大、計算效率低等缺點(diǎn)。

近年來，已有學(xué)者針對礦山邊坡工程初始地應(yīng)力場生成開展研究：劉俊廣等[4]提出礦山深部工程FLAC3D初始地應(yīng)力場生成速度-應(yīng)力邊界法，并驗證速度-應(yīng)力邊界法模擬初始地應(yīng)力合理性；曾范永等[5]通過對FLAC3D計算模型施加自重力，對模型邊界施加切向力和法向力發(fā)現(xiàn)，埋深是決定區(qū)域內(nèi)垂直應(yīng)力主控因素；凌影[6]基于BP神經(jīng)網(wǎng)絡(luò)，利用FLAC3D反演計算區(qū)域初始地應(yīng)力場，結(jié)果與實(shí)測地應(yīng)力基本一致；王金安等[7]充分考慮巖體自重與構(gòu)造應(yīng)力影響，對5類11種邊界條件工況進(jìn)行有限差分FLAC3D模擬加載，反演得出研究區(qū)地應(yīng)力場；裴啟濤等[8]采用遺傳神經(jīng)網(wǎng)絡(luò)法和FLAC3D對研究區(qū)域建立大尺度和小尺度模型，分別進(jìn)行初始地應(yīng)力1次反演和2次反演；祁偉強(qiáng)等[9]采用子模型技術(shù)對防滲墻進(jìn)行模擬，驗證采用子模型法模擬防滲墻具有可行性且精度較高；曹學(xué)興等[10]采用基于動力子模型法的三維非線性有限元，對長河壩心墻堆石壩壩基廊道進(jìn)行研究發(fā)現(xiàn)，壩基廊道加速度和動位移等動力反應(yīng)較小，動力條件下廊道應(yīng)力較靜力條件下有所增大。

初始地應(yīng)力場生成是數(shù)值模擬試驗必然環(huán)節(jié)，眾多學(xué)者在該領(lǐng)域獲得一定研究成果，但關(guān)于利用子模型法獲取簡化后露釆高邊坡數(shù)值模型初始地應(yīng)力場研究較少。因此，本文結(jié)合子模型法和等效應(yīng)力梯度邊界面力法，探討露采高邊坡初始地應(yīng)力場快速生成方法，并基于FLAC3D數(shù)值模擬軟件驗證數(shù)值模型簡化法可行性。研究結(jié)果可為礦山邊坡工程數(shù)值模擬分析提供理論參考。

1 數(shù)值模型簡化方法

露天采礦數(shù)值模擬過程中，因模型取至負(fù)地形僅考慮自重應(yīng)力[11]，但對于露采高邊坡數(shù)值模擬，僅考慮自重應(yīng)力將模型取至負(fù)地形，會導(dǎo)致模型龐大、結(jié)構(gòu)復(fù)雜、劃分單元和節(jié)點(diǎn)數(shù)繁多、計算量巨大、計算效率低等問題。因此，需對整體大模型進(jìn)行簡化，并快速施加初始地應(yīng)力得到簡化子模型，用等效應(yīng)力梯度邊界面力代替簡化部分水平構(gòu)造應(yīng)力，施加側(cè)邊界條件。

1.1 子模型法

子模型法又稱切割邊界位移法或特定邊界位移法，其原理是基于大型復(fù)雜整體模型結(jié)構(gòu)，截取關(guān)鍵部位建立模型，將模型網(wǎng)格進(jìn)一步精細(xì)劃分并進(jìn)行計算分析。

根據(jù)圣維南原理，若實(shí)際分布荷載被等效荷載代替，應(yīng)力和應(yīng)變只在荷載施加位置附近發(fā)生改變。應(yīng)力集中現(xiàn)象一般發(fā)生在荷載集中部位。在露采高邊坡整體大模型中，截取關(guān)鍵部位得到簡化子模型，簡化部分轉(zhuǎn)化為等效梯度應(yīng)力，作為簡化子模型側(cè)邊界，即對簡化一側(cè)施加等效梯度應(yīng)力邊界面力(將應(yīng)力轉(zhuǎn)化為節(jié)點(diǎn)力)。因地應(yīng)力隨深度呈線性增加趨勢，則面力自上而下也呈線性增加趨勢[12-14]。

子模型原理如圖1所示。由圖1可知，自整體模型結(jié)構(gòu)(1)截取子模型結(jié)構(gòu)(2)時，A部分對B部分影響轉(zhuǎn)化為施加在ab面上應(yīng)力，ab兩點(diǎn)間應(yīng)力至坐標(biāo)原點(diǎn)(0,0,0)相對距離成線性比例。

圖1 子模型原理

模型表面施加應(yīng)力視為模型最外層單元所受應(yīng)力，該應(yīng)力轉(zhuǎn)化為節(jié)點(diǎn)力作用在模型最外層單元節(jié)點(diǎn)上。模型邊界外層任一單元體如圖2所示。σ為作用在單元上的應(yīng)力，F(xiàn)為節(jié)點(diǎn)力，對該單元體任何表面有式(1)：

圖2 模型邊界外層任一單元體

(1)

式中：A為單元表面積，m2；i為節(jié)點(diǎn)個數(shù),取值范圍為1，2，…，n。

假設(shè)單元體左側(cè)外表面(面2)是模型外邊界，右側(cè)外表面(面1)與模型內(nèi)部單元接觸。初始應(yīng)力場模擬過程中，模型應(yīng)力邊界恒定，則外層單元邊界(面2)節(jié)點(diǎn)力F保持不變；表層單元體與模型內(nèi)部單元體接觸面(面1)節(jié)點(diǎn)力向模型內(nèi)部傳播，使表面應(yīng)力場分布模式向模型內(nèi)部擴(kuò)散，并最終達(dá)到平衡[15-16]。

1.2 等效應(yīng)力梯度邊界面力法

本文采用等效應(yīng)力梯度邊界面力法，通過施加等效梯度應(yīng)力邊界面力，替代簡化部分水平構(gòu)造應(yīng)力。假設(shè)圖1A部分對a點(diǎn)應(yīng)力為S0，對b點(diǎn)應(yīng)力S如式(2)所示：

S=S0+gxx+gyy+gzz

(2)

式中：S為A部分對b點(diǎn)應(yīng)力，Pa；gx,gy,gz分別為應(yīng)力在x，y，z方向變化，Pa；x，y，z為兩點(diǎn)到全局坐標(biāo)系原點(diǎn)在x，y，z方向相對距離，m。

根據(jù)A.Heim法則及側(cè)壓力系數(shù)λ得到垂直應(yīng)力，如式(3)所示：

σv=ρgZ，σh=λσv

(3)

式中:σv為垂直應(yīng)力，Pa；σh為水平應(yīng)力，Pa；ρ為巖體密度，g/m3；g為重力加速度，取10 m/s2；Z為深度，m；λ為側(cè)壓力系數(shù),反映巖體原位應(yīng)力狀態(tài)。入側(cè)壓力原則上由試驗確定，但精確度較低，不易實(shí)現(xiàn)。因此，根據(jù)彈性力學(xué)理論，側(cè)壓力系數(shù)λ如式(4)所示：

(4)

式中：u為上浮巖體泊松比。

結(jié)合式(3)～(4)得到式(5)：

(5)

由式(1)與式(5)可得式(6)：

(6)

2 工程概況

礦區(qū)主要開采礦層為二疊系下統(tǒng)棲霞組(P1q)和茅口組(P1m)水泥用灰?guī)r。在棲霞組底部為灰黑色薄層狀含粉砂質(zhì)泥質(zhì)灰?guī)r，局部夾中厚層狀泥粉晶灰?guī)r，粉晶泥質(zhì)結(jié)構(gòu)，間夾波狀泥質(zhì)條帶。含粉砂質(zhì)泥質(zhì)灰?guī)r層理發(fā)育，構(gòu)成軟弱結(jié)構(gòu)面，巖層呈微弧形單斜層狀產(chǎn)出。受黃山斷層F1及魚洞口斷層F2構(gòu)造影響，于東部羅溝近峨眉山玄武巖一帶略有倒轉(zhuǎn)，向東倒轉(zhuǎn)逐漸加劇，產(chǎn)狀直立乃至倒轉(zhuǎn)。黃山某石灰石礦山邊坡巖體位于礦區(qū)中部羅溝高程720 m以上，開挖邊坡最高高程約895 m，殘坡積物(Q4)主要分布在高程895 m，底部主要為灰?guī)rP15?；轮饕泼鏋槎B系炭質(zhì)泥頁巖形成的軟弱夾層面，該軟弱夾層力學(xué)強(qiáng)度低、水理性質(zhì)差。礦區(qū)夏季多雨，為防止礦區(qū)發(fā)生季節(jié)性泥石流，在坡底640平臺預(yù)設(shè)防泥石流內(nèi)排凹槽。同時，該地層為單斜地層，南高北低，南部存在陡崖。地質(zhì)地形及剖面如圖3所示，剖面取至負(fù)地形的3個剖面跨度明顯。

圖3 地質(zhì)地形及剖面

3 模型計算分析

3.1 構(gòu)建FLAC3D模型

1)FLAC3D數(shù)值模擬流程

利用數(shù)值模擬軟件建立地質(zhì)力學(xué)模型，根據(jù)試驗數(shù)據(jù)及地質(zhì)地形資料，對參數(shù)進(jìn)行合理賦值，并進(jìn)行應(yīng)力應(yīng)變計算分析。FLAC3D數(shù)值模擬流程如圖4所示。

圖4 FLAC3D數(shù)值模擬流程

2)地質(zhì)模型建立

本文選取典型剖面3-3′構(gòu)建計算模型，其工程地質(zhì)剖面如圖5所示。

圖5 典型剖面3-3′工程地質(zhì)剖面

由于預(yù)開挖防泥石流內(nèi)排凹槽所處位置特殊，計算剖面取至負(fù)地形跨度大，導(dǎo)致整體計算模型龐大，且很大一部分防泥石流內(nèi)排凹槽不在研究范圍，對礦山高邊坡影響作用無法凸顯，開挖造成地應(yīng)力重新分布。因此，截取邊坡已開挖防泥石流內(nèi)排凹槽最高處為頂點(diǎn)的子模型，使研究部分更易被觀察。邊坡滑坡主要滑移面為二疊系炭質(zhì)泥頁巖形成的軟弱夾層面，其中3條夾層面埋深較深，作用不明顯，因此，只考慮最上面1條軟弱夾層。截取計算模型突出部位，如圖6所示。

圖6 截取模型突出部位

3)計算模型及網(wǎng)格剖分

為消除邊界效應(yīng)影響，整體大模型剖面3-3′計算區(qū)域取x軸0～710 m，z軸豎直向上取0～286 m；簡化子模型剖面3-3′計算區(qū)域取x軸計算范圍0～241 m，z軸豎直向上取0～310 m。

根據(jù)現(xiàn)場實(shí)測地形地質(zhì)資料，建立整體計算模型。計算域剖分58 950個4邊形單元，共計71 928個節(jié)點(diǎn)，對凹槽、裂縫產(chǎn)生區(qū)域及軟弱夾層采用較密集單元，模型其余部分采用網(wǎng)格劃分技術(shù)進(jìn)行過渡。同理，簡化子模型計算域剖分28 385個4邊形單元，共計34 728個節(jié)點(diǎn)。整體計算模型與簡化子模型結(jié)構(gòu)對比見表1。選取最易發(fā)生變化的7個坡頂點(diǎn)(編號1～7)作為監(jiān)測點(diǎn)，對其水平及豎直方向應(yīng)力實(shí)施監(jiān)測。整體計算模型與簡化子模型計算數(shù)值網(wǎng)絡(luò)及監(jiān)測點(diǎn)布設(shè)如圖7所示。

表1 整體計算模型與簡化子模型結(jié)構(gòu)對比

圖7 整體計算模型和簡化子模型計算數(shù)值網(wǎng)格及監(jiān)測點(diǎn)布設(shè)

由表1可知，簡化子模型較整體計算模型剖分單元及節(jié)點(diǎn)分別簡化51.85%，51.72%，可有效節(jié)約計算時間，提高計算效率。

4)計算參數(shù)確定

根據(jù)工程地質(zhì)條件、巖石力學(xué)試驗結(jié)果及現(xiàn)場勘察報告，某石灰石礦巖體物理力學(xué)參數(shù)見表2[17]。

5)邊界條件

簡化數(shù)值網(wǎng)格模型計算域底部采用固定約束，右側(cè)采用法向約束，左側(cè)施加梯度應(yīng)力邊界面力。由于力的施加，在模型平衡過程中產(chǎn)生極大位移，可以在初始應(yīng)力場形成后，對模型所有節(jié)點(diǎn)位移清零再進(jìn)行后續(xù)計算[18-19]。由表2可知，石灰?guī)rμ1=0.24，ρ1=2.68 g/cm3，軟弱夾層μ2=0.28，ρ2=1.92 g/cm3，g取10，簡化部分取平均高度Z=4.2 m，將數(shù)值帶入式(6)可得S0=2.4×105Pa，gz=5.68，左側(cè)面力梯度利用式(2)編入FLAC3D命令。

表2 某石灰石礦巖體物理力學(xué)參數(shù)

3.2 計算結(jié)果分析

采用FLAC3D數(shù)值模擬軟件，輸出水平及豎直方向應(yīng)力云圖如圖8～11所示。由圖8～11可知，整體大模型和簡化子模型開挖后，水平及豎直方向應(yīng)力云圖相似，簡化子模型在施加面力影響下，邊界附近水平方向應(yīng)力相對較大，符合圣維南原理，誤差隨截取子模型邊界距負(fù)地形距離減小而減小。

圖8 整體大模型水平方向應(yīng)力

圖9 整體大模型豎直方向應(yīng)力

圖10 簡化子模型x水平方向應(yīng)力

圖11 簡化子模型z豎直方向應(yīng)力

整體大模型與簡化子模型監(jiān)測點(diǎn)水平及豎直方向應(yīng)力值見表3，2種模型對比折線如圖12～13所示。定義相對誤差，如式(7)所示：

圖12 水平方向應(yīng)力值對比

表3 整體大模型及子模型開挖監(jiān)測結(jié)果表

(7)

整體大模型與簡化子模型監(jiān)測點(diǎn)開挖后水平及豎直方向應(yīng)力值整體相差較小，由式(7)計算相對誤差可知，整體誤差均在工程誤差允許范圍內(nèi)。

圖13 豎直方向應(yīng)力值對比

4 結(jié)論

1)簡化子模型可對礦山高邊坡整體大模型進(jìn)行有效簡化，快速獲得初始地應(yīng)力場。

2)通過對比整體大模型與簡化子模型初始地應(yīng)力場模擬結(jié)果可知，整體大模型與簡化子模型水平及豎直方向應(yīng)力相對誤差較小，驗證本文模型簡化法合理有效。

3)整體大模型存在體積(面積)過大、節(jié)點(diǎn)單元過多、計算效率低下等問題；相較整體大模型，簡化子模型單元數(shù)和節(jié)點(diǎn)數(shù)減少約52%，有效縮短計算時間，提高計算效率。