考慮弱化效應的加高擴容尾礦壩地震可靠度分析及風險評價

李春立,李亮*,徐亮,趙民,王超,蔡鳳珍

(1.青島理工大學土木工程學院,青島 266520; 2.青島市勘察測繪研究院,青島 266033)

誘發尾礦壩潰壩的因素和機理復雜繁多,究其原因,主要與尾礦壩筑壩材料的不確定性以及地震、強降雨等復雜的因素相關[1-4]。由于中國地震頻發,給尾礦壩構成了極大威脅。地震對尾礦壩穩定性的影響主要表現為地震的慣性力使處于極限平衡狀態的尾礦壩產生變形以及地震循環剪切作用使壩體震動液化[5-10]。鑒于此,國家有關規范中明確提出[11-12],為保證尾礦壩的安全,要求對設防烈度為Ⅶ度及以上、尾礦庫等級在三級及以上的尾礦壩進行抗震性能分析。在尾礦壩地震穩定性分析方面,國內外學者開展了一系列研究,以獲得尾礦壩的安全系數、位移、液化變形及大壩的應力應變狀態。例如,楊安銀等[13]采用有限元Quake/W軟件分別研究了尾礦壩在三種地震作用下的加速度、位移、液化區域、安全系數的動力響應特性。秦曉鵬等[14]采用擬靜力法和有限差分法計算地震作用下的邊坡時程安全系數,并采用PL-Finn模型結合Flac3D研究了尾礦壩的液化變形。Naeini等[15]分別采用有限元Quake/W軟件和有限元Sigma/W軟件對大壩進行動態和應力重分布分析,并確定壩體的永久位移。朱明明[16]考慮隨機地震動和尾礦材料參數的二維空間變異性,對尾礦壩進行了動力可靠度分析和地震敏感性分析。研究成果有助于合理評估尾礦壩地震安全性能,有力保障礦山安全。

近年來,由于礦山繼續擴大再生產,尾礦庫亟需進行加高擴容,評估尾礦壩加高擴容后的地震安全性具有指導意義。尹光志等[17]采用有限元時程分析法,分別對云南大沙河尾礦庫加高擴容壩體在現狀條件下、加高中增設排滲措施和直接加高后三種工況進行了動力分析和抗震性能研究。Huaman等[18]使用FLAC2D對處于秘魯高地震活動區、采用下游法筑壩的擴容尾礦壩進行了地震災害評估。

以往的研究中重點關注尾礦壩地震安全系數或永久位移等指標,忽略了潰壩后滑動體的大小、運移過程以及堆積狀態對下游構筑物帶來的風險研究,此外,地震過程中導致尾礦壩筑壩材料的弱化,進而會加劇尾礦壩潰壩風險。因此,為合理考量上述因素對尾礦壩地震安全性能影響,現采用極限平衡擬靜力法與光滑粒子流體動力學方法(smoothed particle hydrodynamics SPH)聯合進行尾礦壩地震可靠度分析[19-27]。考慮地震對尾礦壩筑壩材料的弱化效應,利用極限平衡擬靜力法與蒙特卡羅方法尋求可能的失效樣本,進而采用自主研發的SPH程序模擬失效樣本下尾礦壩失穩過程以及最終堆積狀態,以滑動面積評估尾礦壩的風險水平。

1 尾礦壩地震可靠度分析

1.1 極限平衡擬靜力法

極限平衡擬靜力法是用于評價地震作用下巖土結構物穩定性的方法之一[28-29]。該法忽略壩體的動態響應,將地震作用簡化為一靜力作用在壩體上,結合極限平衡方法進行穩定性分析時,將地震作用等效為水平和豎向地震力分別施加于每個土條的重心。主要考慮水平地震力的作用,如圖1所示,某滑動區域被劃為n個垂直土條,其中作用于土條i的水平地震力代表值Qh,i為

Ti為土條底部法向力;Zi為土條底部剪切力;Ei和Ei+1為土條間法向應力;βi為土條底部傾角,(°),i=1,2,…,n ;r為圓弧滑動半徑;O為圓弧滑動圓心

(1)

式(1)中:ξ為地震作用效應折減系數,依據《水工建筑物抗震設計標準》(GB 51247—2018)規定[30],除采用動力法計算鋼筋混凝土結構外,ξ取0.25;αi為

土條i地震慣性力的動態分布系數,需根據圖2(a)所示的壩體地震慣性力動態分布系數示意圖確定;ah為水平地震加速度;kh為水平地震力系數,可通過查閱表1、表2確定;g為重力加速度;Wi為土條i的重力。鑒于Slope/W軟件中的極限平衡擬靜力法未考慮地震慣性力的動態分布,為便于計算,本文中將圖2(a)所示地震慣性力動態分布圖簡化圖2(b)所示的矩形分布圖。

表1 工程場地地震烈度與地震動峰值加速度(PGA)對照表[31]

表2 工程場地地震烈度與水平地震力系數(kh)對照表[31]

H為壩高,αi為某壩高Hi時地震慣性力的動態分布系數;am為壩體頂部的分布系數

以Bishop法為例,采用極限平衡擬靜力方法計算壩體安全系數FS,計算公式為

(2)

式(2)中:φi為土條底部的內摩擦角,(°);λi為土條長度,m;c為土條黏聚力,kN。在給定的尾礦壩筑壩材料參數下,可以根據式(1)和式(2)計算地震作用下壩體的安全系數,為考慮筑壩材料參數的隨機性,需要結合蒙特卡羅抽樣方法生成一系列樣本,分別采用式(1)和式(2)計算每一樣本下的壩體安全系數值,進而計算壩體的失效概率和潰壩風險。為提高計算效率,基于Win-BatchTM平臺,編寫批處理程序調用Slope/W軟件進行壩體失效概率計算。

1.2 基于Win-BatchTM語言的失效概率計算

采用1.1節極限平衡擬靜力法計算的安全系數確定壩體的極限狀態函數,進而結合蒙特卡羅方法計算尾礦壩的失效概率,即

G(X)=FS(X)-1=0

(3)

式(3)中:G(X)為極限狀態函數,X=(x1,x2,…,xm)為尾礦壩地震可靠度分析時考慮的巖土體參數變量,譬如黏聚力c和內摩擦角φ等,其概率分布特征通常假定為對數正態分布。依據土體強度參數的平均值、標準差及分布特征進行N次抽樣,產生N組隨機樣本,X1,X2,…,XN。依次將每組隨機樣本輸入Slope/W中計算尾礦壩的安全系數,若Gi(X)<0,則隨機樣本Xi視為失效樣本,統計失效樣本個數M,尾礦壩的失效概率pf≈M/N。pf精度計算式為

(4)

式(4)中:pf為尾礦壩的失效概率;Cpf為失效概率pf的變異系數;N為總抽樣次數。

基于Win-BatchTM平臺編寫批處理程序,調用Slope/W軟件批量計算尾礦壩的安全系數FS,并結合蒙特卡羅方法計算失效概率的基本流程如圖3所示,主要步驟如下。

圖3 基于Win-BatchTM語言失效概率計算流程圖

(1)確定尾礦壩的幾何模型,如壩高、土層分層,土體強度參數和工程場地地震烈度,在Slope/W軟件中建立尾礦壩數值分析模型并確定kh,將尾礦壩分析模型保存為.xml源文件。

(2)根據尾礦壩土體強度參數平均值和標準差,抽取N組隨機樣本,X1,X2,…,XN,將N組隨機樣本保存為N個.xml源文件,并批量導入N個文件夾中;

(3)通過批處理程序調用Slope/W軟件依次計算N組隨機樣本的安全系數FS1,FS2,…,FSN;

(4)基于尾礦壩的極限狀態函數,統計失效樣本個數M,估算失效概率pf=M/N。

2 考慮弱化效應的尾礦壩失穩風險評價

2.1 土體強度的地震弱化效應

研究表明,地震時循環往復的剪切作用會導致土體強度發生弱化,即“土體軟化”[32-33]。為合理考慮土體強度的弱化效應對尾礦壩失穩風險的影響,定義土體強度的弱化系數為

(5)

式(5)中:su為地震前尾礦壩各土層黏聚力和內摩擦角,可依據實際地勘報告或工程經驗確定;s′u為土體弱化后的黏聚力和內摩擦角,土體強度弱化系數η可由現場試驗或實驗室試驗進行測定。本文中假定不同的土體強度弱化系數η分別為0.75、0.5和0.15,結合工程實例探討其對尾礦壩失穩風險的影響。

2.2 光滑粒子流體動力學方法確定滑動面積

尾礦壩失穩后的運動過程需通過大變形方法來模擬。光滑粒子流體動力學方法作為一種無網格方法,克服了傳統有限元方法網格畸變對計算結果的不利影響,通過一系列攜帶質量、位移、速度等變量的粒子模擬計算區域,各粒子之間的相互作用通過核函數來模擬,基于質量守恒和動量守恒,引入不同宏觀本構,該方法可模擬壩體滑動過程。限于篇幅,關于SPH方法的基本原理在此不再贅述,讀者可參閱文獻[34-37]。本團隊自主研發了SPH程序并成功地進行了香港翡翠道滑坡模擬,并與現場監測結果對比驗證了所研發程序。除此之外,應用所研發的SPH程序進行了一系列的滑坡過程模擬。

采用SPH方法計算壩體滑動面積的大體思路如下:針對M組失效樣本,分別將每組樣本值輸入自主研發的SPH程序,保存每次迭代計算的粒子滑動位移uj,對粒子滑動位移進行降序排列,u1,u2,…,uq,q為粒子總數。設置粒子的臨界位移值θ,當某粒子滑動位移uj超過臨界位移值θ時,該粒子滑動面積計入尾礦壩的失穩滑動面積Si。表達式為

(6)

式(6)中:Si為尾礦壩的失穩滑動面積;dj為第j個粒子面積;q為粒子總數;uj為粒子的滑動位移;θ為粒子的臨界位移值,可由工程經驗給出或者取尾礦壩壩高H的某一分數值,本文中暫取θ=5%H、10%H、15%H。

2.3 尾礦壩失穩風險計算

尾礦壩失穩風險R可通過其失效概率pf和失穩后果C的乘積來表征。采用蒙特卡羅方法計算尾礦壩M種失效模式的失效概率pf,i,由第1節可知,尾礦壩失效概率pf=M/N,其中每種失效模式的失效概率為pf,i=1/N(i=1,2,…,M)。用式(6)計算所得Si量化其失穩后果Ci(i=1,2,…,M),尾礦壩的失穩風險R為

(7)

本文所提尾礦壩失穩風險的計算流程,具體如圖4所示。

圖4 尾礦壩失穩風險的計算流程圖

3 案例分析

3.1 工程概況

云南省斑毛溝擴容尾礦庫位于哀牢山脈西側,屬山谷型尾礦庫。在原有舊尾礦壩(簡稱“舊壩”)的下游選址建設高達198 m的新尾礦壩(簡稱“新壩”),使得原有尾礦庫由三級庫擴容為新二級尾礦庫,總壩高H=230 m,總擴容達1 087.3萬m3。根據地勘報告顯示,尾礦壩土層分別為壩基、初期壩、子壩、尾粉土和尾粉質黏土,場地地震動峰值加速度PGA為0.1g,地震動反應譜特征周期為0.45 s。擴容尾礦壩施工期間主要分為三個階段,分別為新庫區擴容至原庫區初期壩壩頂、新庫區擴容至原庫區子壩中點、新庫區擴容至原庫區壩頂。本文分別對擴容尾礦壩的三個施工階段進行穩定性分析。

3.2 擴容尾礦壩擬靜力穩定性分析

為便于數值模擬,簡化其壩基部分進行建模(圖5)。考慮到尾礦庫內水位對尾礦壩穩定性的影響,首先利用有限元Seep/W軟件對各施工階段進行滲流穩定性分析,之后再耦合Slope/W軟件進行安全系數的計算。采用極限平衡擬靜力法模擬擴容尾礦壩的地震性能,場地地震動峰值加速度為0.1g,根據中國地震動參數區劃圖(GB 18306—2015)[31]的地震動峰值加速度(peak ground acceleration,PGA)與工程場地地震烈度對照表(表1)確定尾礦壩的地震烈度為Ⅶ度,由《水工建筑物抗震設計標準》(GB 51247—2018)[30]可知,壩頂動態分布系數am=3.0。再依據表2確定水平地震力系數kh為0.1,地震作用效應折減系數ξ=0.25。本算例壩高H=230 m,0.6H高程的動態分布系數取1.0+(am-1)/3,壩底動態分布系數取1.0,等效后的地震慣性力動態分布系數為1.73,為安全起見,最終確定地震慣性力的動態分布系數為2.0。

圖5 擴容尾礦壩簡化模型

以擴容尾礦壩的三個施工階段為分析工況,工況一為新庫區擴容至原庫區初期壩壩頂、工況二為新庫區擴容至原庫區子壩中點、工況三為新庫區擴容至原庫區壩頂。根據地勘報告的各土層參數(表3),分別計算三種工況條件下的最小安全系數FS。

表3 擴容尾礦壩土層物理力學參數

圖6匯總了三種工況條件下擴容尾礦壩的最小安全系數及其滑動面。其中,工況一下最小安全系數為1.56,其對應的滑動面(稱之為臨界滑動面)位于舊壩;工況二下最小安全系數為1.60,其臨界滑動面位于新壩;工況三下最小安全系數為1.22,其臨界滑動面位于新壩。結果表明:隨尾礦壩擴容高度的上升,舊壩的穩定性基本保持不變,而新壩的穩定性持續降低,自2.03降低至1.22。該尾礦壩總庫容1 087.30萬m3,壩高230 m,由構筑物抗震設計規范(GB 50191—2012)[12]尾礦壩抗震等級可知(表4),其抗震等級為二級,由表5可知,規范要求的最小安全系數為1.15。三種工況條件壩體最小安全系數均大于1.15,滿足規范要求。但需要注意的是,工況三下,其最小安全系數與規范要求值接近,需要重點關注其穩定性,下節重點針對工況三進行可靠度分析與風險評價。

表4 尾礦壩抗震等級[12]

表5 尾礦壩地震穩定性最小安全系數值[12]

圖6 擴容尾礦壩三種工況下尾礦壩地震穩定安全系數

3.3 擴容尾礦壩工況三失效概率計算

考慮尾礦壩筑壩材料的不確定性,采用蒙特卡羅方法對工況三下擴容尾礦壩進行可靠度分析。假定子壩、尾粉土、尾粉質黏土的黏聚力c、內摩擦角φ均為對數正態分布,以表2中的各土層參數為平均值,變異系數均假定為0.1,隨機生成1 000組樣本。采用圖3所示流程計算1 000組隨機樣本下壩體的最小安全系數并計算其失效概率。圖7以最小安全系數為橫軸,以該安全系數的頻率為縱軸,繪制最小安全系數直方圖。由圖7可見,最小安全系數的平均值為1.218,標準差為0.104,失效樣本個數為10個,其失效概率為1%,由式(4)可知,計算所得失效概率變異系數Cpf≈0.3。

圖7 工況三下安全系數的直方圖

3.4 擴容尾礦壩失穩風險評價

3.4.1 尾礦壩失穩過程模擬

尾礦壩一旦潰壩,將會導致尾礦庫內大量的有害物質流向下游,對下游的居民、環境、水源等造成不可估量的損失。基于此,以第10個失效樣本為例,取弱化系數η=0.15,圖8分別給出了t=1,3,5,10 s時的滑動狀態。t=1 s時,尾礦壩的最大滑動位移為2.17 m,小于尾礦壩的臨界位移值θ,此時整體處于穩定狀態;t=3 s時,最大滑動位移為18.99 m,主要位于新壩處的子壩和尾粉質黏土部分,舊壩部分的滑動位移為4～6 m;t=5 s時,最大滑動位移為37.79 m,舊壩部分的滑動位移約為8～10 m,此時部分壩體已經越過初期壩,有向下游滑動的趨勢;t=10 s時,最大滑動位移為241.56 m,部分壩體已經滑落至初期壩的邊緣,擴容尾礦壩整體開始變形。

圖8 η=0.15的擴容尾礦壩失穩滑動過程

根據上述尾礦壩的失穩滑動過程,對尾礦壩最大滑動速度進行量化。其中,1 s內壩體尚未滑動,其滑動速度為0,3～5 s內的最大滑動速度為9.40 m/s,5～10 s的最大滑動速度為40.95 m/s。在5 s內,尚未滑至初期壩,滑動速度相對較小,整體處于可控狀態。5～10 s內部分尾礦壩逐漸滑落至初期壩的底部邊緣,其滑動速度迅速增加,導致滑動位移急劇增大。10 s后滑動速度將持續增大,壩體及庫內有害物質將迅速殃及下游部分區域。

3.4.2η=1時尾礦壩失穩風險評價

圖9 不同θ下擴容尾礦壩失效樣本對應的風險Ri

在本例中,取θ=15%H=34.5 m時,只有當失穩后尾礦壩的滑動位移超過34.5 m時,才會產生風險,如圖10所示,10個失效樣本下的尾礦壩的最大位移僅為25.69 m,未超過34.5 m,因此風險為0。圖10還表明,隨著失效樣本安全系數的降低,尾礦壩滑動位移值呈非線性增長趨勢,當考慮地震弱化效應后,失效樣本安全系數會降低,因此需要重點關注考慮弱化效應后的尾礦壩失穩風險量化。

圖10 擴容尾礦壩失效樣本的最大滑動位移

3.4.3η對尾礦壩失穩風險的影響

為研究η對計算結果的影響,取η=0.75、0.5、0.15,分別采用SPH模擬不同弱化系數下,尾礦壩失效樣本的失穩滑動面積,并計算在不同臨界位移值下的風險。以土體強度弱化系數為橫軸,尾礦壩的風險為縱軸,繪制不同臨界位移值下,其風險隨弱化系數的變化曲線,如圖11所示。臨界位移值θ=5%H時,η=1.0、0.75、0.5、0.15時,尾礦壩的失穩風險R分別為309.17、380.31、437.72、842.95 m2;臨界位移值θ=10%H時,η=1.0、0.75、0.5、0.15時,尾礦壩的失穩風險R分別為36.14、300.97、418.23、824.74 m2;臨界位移值θ為15%H時,η=1.0、0.75、0.5、0.15時,尾礦壩的失穩風險R分別為0、205.12、398.84、807.16 m2。

圖11 不同θ下擴容尾礦壩風險隨η的變化曲線圖

圖12 不同η下θ取值對擴容尾礦壩風險的影響

4 結論

結合某擴容尾礦壩工程實例,考慮地震對尾礦壩筑壩材料的弱化效應,采用極限平衡擬靜力法與光滑粒子流體動力學方法結合進行尾礦壩地震可靠度分析及風險評價。利用極限平衡擬靜力法與蒙特卡羅方法尋求可能的失效樣本,進而采用自主研發的SPH程序模擬失效樣本下尾礦壩失穩過程以及最終堆積狀態,以滑動面積評估尾礦壩的風險水平,并得出以下結論。

(1)隨尾礦壩擴容高度的上升,舊壩的穩定性基本保持不變,而新壩的穩定性持續降低,擬靜力法所得安全系數自2.03降低至1.22。

(2)在工況三下,即擴容高度完成后,尾礦壩最小安全系數與規范要求值接近,需要重點關注其穩定性,經可靠度分析,當尾礦壩筑壩材料變異系數為0.1時,其相應的失效概率為1%,該失效概率值的變異系數為0.3左右。

(3)不考慮地震作用對土體的弱化效應,擴容尾礦壩的失穩風險R=309.17 m2,尾礦壩的失穩風險R隨臨界位移值θ的增大而減小。

(4)考慮弱化效應后,隨弱化系數減小,擴容尾礦壩的失穩風險增大;弱化系數大于0.5時,臨界位移值θ對擴容尾礦壩失穩風險影響顯著。