蔣 魁，王 帆，蔣興中，景鵬旭

(1.中國有色桂林礦產(chǎn)地質(zhì)研究院有限公司，廣西 桂林 541004；2.桂林礦產(chǎn)地質(zhì)研究院工程有限公司，廣西 桂林 541004；3.廣西壯族自治區(qū)第一地質(zhì)隊，廣西 桂林 541004；4.中國地震應(yīng)急搜救中心，北京 100049)

當(dāng)前研究礦山邊坡穩(wěn)定性時，關(guān)于巖體非線性強度特征及其對邊坡穩(wěn)定性影響的相關(guān)研究尤為重要?，F(xiàn)場勘測結(jié)果表明：邊坡后側(cè)巖體的拉應(yīng)力區(qū)域會顯著影響礦山邊坡的穩(wěn)定性。關(guān)于抗拉強度的研究發(fā)現(xiàn)[1]，采用兩種方法可以有效地分析邊坡拉應(yīng)力區(qū)的抗拉強度，使分析結(jié)果更加趨于真實的破壞條件，其一是對巖土體的抗拉強度進(jìn)行一定程度的折減，即將線性破壞曲線優(yōu)化為非線性強度包絡(luò)線；其二則是在邊坡后緣的拉應(yīng)力區(qū)引入不同形態(tài)的裂縫。 對于后者考慮裂縫的存在，MICHALOWSKI[2]將裂縫的擴展過程與整體的滑落狀態(tài)結(jié)合，研究結(jié)果表明：對于高陡邊坡或受靜水壓力影響顯著的邊坡，裂縫的存在大大增加了失穩(wěn)的風(fēng)險。UTILI[3]分析了預(yù)設(shè)垂直裂縫邊坡的地震穩(wěn)定性的影響，研究顯示地震荷載對邊坡的穩(wěn)定性影響尤為顯著。對于前者所考慮的抗拉強度折減方法而言，當(dāng)前的處理方式，大多在變分法或極限分析法的基礎(chǔ)上進(jìn)行。BAKER等[4]對拉剪區(qū)域進(jìn)行了一定的折減，并結(jié)合變分法對強度包絡(luò)線中的張拉區(qū)域進(jìn)行了相關(guān)的穩(wěn)定性分析。PARK等[5]以極限分析上限法為理論基礎(chǔ)，選擇圓弧對經(jīng)典M-C準(zhǔn)則中的強度包線進(jìn)行優(yōu)化，討論了巖土體的抗拉強度的變化對邊坡穩(wěn)定性和破裂形態(tài)的影響。 在使用極限分析法時，僅需對邊坡達(dá)到極限平衡狀態(tài)時候進(jìn)行計算，計算過程相比有限元法更加簡潔，且計算速度快。

地震的作用往往是誘發(fā)邊坡失穩(wěn)的直接因素，TERZAGHI[6]提出了擬靜力法來分析地震邊坡的穩(wěn)定性，由于其概念清晰且計算方便，逐漸成為當(dāng)前的主流方法。 但是擬靜力法忽略了地震的動態(tài)過程，因此無法有效反應(yīng)邊坡的永久位移值。NEWMARK[7]1965年提出了計算邊坡在動力作用下的滑動位移的方法，并據(jù)此得到了邊坡的屈服加速度的推導(dǎo)和永久位移分析理論，從而有效分析地震邊坡的穩(wěn)定性。但由于Newmark法的核心為剛性滑塊法，忽視了邊坡內(nèi)部變形的影響，此后不斷有對Newmark法進(jìn)行改進(jìn)的研究[8-10]。

為研究抗拉強度折減對礦山邊坡在地震作用下穩(wěn)定性的影響，本文在前人的基礎(chǔ)上，首先通過極限分析上限法對礦山邊坡臨界的破壞形態(tài)以及功能方程進(jìn)行推導(dǎo)，計算得到不同破壞條件下的邊坡屈服加速度的變化趨勢，最后針對賀州市某礦山項目中典型礦山邊坡，選取實測地震動分析不同參數(shù)以及張拉強度折減系數(shù)對礦山邊坡穩(wěn)定性的影響。

1 巖土體強度包絡(luò)線

經(jīng)典的M-C準(zhǔn)則在主應(yīng)力空間可以表示為一個不規(guī)則六邊形截面的錐體，如圖1所示，當(dāng)巖體材料以化學(xué)鍵結(jié)合時，黏結(jié)性表現(xiàn)為單軸抗壓強度和單軸抗拉強度。金字塔的頂點表示三軸(各向同性)抗拉強度，而單軸抗拉強度和抗壓強度表示為主應(yīng)力軸與屈服曲面的交點。PAUL[11]對經(jīng)典的M-C屈服面進(jìn)行了修改，從而限制了抗拉強度，其屈服面如圖2所示。折減后的部分拉力截斷單軸抗拉強度ft可表示為圖3中的應(yīng)力圓C1。由圖3和圖4中的莫爾圓轉(zhuǎn)換至圖1和圖2中的主應(yīng)力空間是一個較為復(fù)雜的過程，圖2中截面ABCDEF上AB邊上的應(yīng)力狀態(tài)由圖4)中的應(yīng)力圓C2表示，而AB上各個位置由其中主應(yīng)力的大小確定(σ1≥σ2≥σ3)。例如，圖3和圖4中的應(yīng)力圓C1并且σ2=0的情況對應(yīng)于圖2中ABCDEF截面與σ3軸的交點，而圖4中點P則與圖2中角點C的應(yīng)力狀態(tài)相對應(yīng)。極限狀態(tài)時的屈服面由圖3和圖4中光滑的屈服面映射至圖1和圖2中具有棱角的屈服面。采用此拉應(yīng)力區(qū)的截斷方式適用于極限分析中的解析運動學(xué)分析方法，但因為主應(yīng)力空間中屈服面的棱角處的法向量不唯一，其在數(shù)值方法中的適用性還需深入研究。傳統(tǒng)的M-C強度包絡(luò)線如圖3所示，M-C函數(shù)為簡單的直線與圓弧部分的組合，在M-C計算中，三軸抗拉強度的大小為：f3t=c/tanφ，而單軸抗拉強度大小為ft=2ccosφ/(1+sinφ)。研究中當(dāng)考慮抗拉強度的折減時，牽引矢量從直線變化為圓弧，且張拉應(yīng)力區(qū)域的剪脹角δ從φ變化到π/2，當(dāng)δ=φ時代表線性的M-C屈服準(zhǔn)則。為了量化展示折減后抗拉強度，本文引入了抗拉強度折減系數(shù)ξ，具體見式(1)。

(1)

(2)

式中，fc=2ccosφ/(1-sinφ)，為單軸抗壓強度。

圖1 經(jīng)典M-C屈服Fig.1 The classical M-C criteria

圖2 土體受張拉折減屈服Fig.2 Tension cut-off

圖3 經(jīng)典M-C包絡(luò)線Fig.3 The classical M-C envelop

圖4 土體受張拉折減包絡(luò)線Fig.4 Tension cut-off

圖5 考慮抗拉強度折減的邊坡破壞模式Fig.5 Failure mechanism of tension cut-off slope

2 破壞機制與計算分析

基于極限分析的運動學(xué)理論，本文采用張拉剪切破壞機制的對數(shù)螺線模式[12]表示邊坡的滑動面。圖5表示在受地震荷載作用的臨界均質(zhì)礦山邊坡模型，其中，β為坡體角度，H為坡體高度；其余參數(shù)為內(nèi)摩擦角φ、黏聚力c和容重γ。邊坡的破壞機制如圖5所示，塊體BDEC以中心O為中心沿著滑裂面DEC旋轉(zhuǎn)，而整個破壞為過破趾的破壞，DEC部分為臨界滑面。研究發(fā)現(xiàn)，張拉應(yīng)力更容易出現(xiàn)在坡頂附近，而滑面上部DE部分由受拉伸區(qū)折減部分強度包絡(luò)線而引起，下方較低的部分(EC)是由線性部分的包絡(luò)線引起。DE部分剪脹角δ從點D的最大值δmax向點E的φ不斷折減，見式(3)。

(3)

式中，r0和θ0分別為對數(shù)螺旋線起始點A到旋轉(zhuǎn)中心的半徑的距離AO和AO與水平參考線的夾角(如圖5所示)，其中，θh為對數(shù)螺旋線終止點到旋轉(zhuǎn)中心的半徑與水平參考線的夾角；θm為OD與水平參考線的夾角；θtc為OE與水平參考線的夾角；δ為剪脹角，其詳細(xì)的分布形式見式(4)。

(4)

極限分析上限法基于內(nèi)外能量耗散率的平衡，見式(5)。

Wr+Ws+Wqt=D

(5)

式中：Wr為土體自重狀態(tài)下的能量耗散率；Ws為擬靜力地震荷載作用下的能量耗散率；Wqt為坡頂部荷載的功率；D為滑體BDEC的內(nèi)能耗散率。外部功率Wr和Ws可表示為塊體ABC的功率減去塊體ADE的功率，見式(6)。

Wqt=qt×fq

(6)

圖5中區(qū)域OCA，AOB，BOC，OEA，AOD，DOE的能量耗散功率具體可參考文獻(xiàn)[4]，其中：fw與fs均為包含變量：β，θ0，θtc，θm，θh的函數(shù)方程。fq計算為式(7)。

(7)

(8)

式中：fd為包含變量θtc，θh，φ的函數(shù)，而fdc為包含變量θtc，θm，φ的函數(shù)。詳細(xì)的關(guān)于fw，fs，fd，fdc及fq的表達(dá)式可參考文獻(xiàn)[2]和文獻(xiàn)[3]。

將上述公式整理，即可得到屈服加速度Kc，見式(9)。

(9)

式中，Kc為一個含θ0，θtc，θm，θh以及c/γH、坡腳β和內(nèi)摩擦角φ等參數(shù)的函數(shù)，則屈服加速度Kc的獲取為一系列變量組合下的最小值。即限定條件為式(10)。

0<θ0<θm<θtc<θh<π

(10)

其中，L2/r0，見式(11)。

(11)

本文采用Matlab對上述推導(dǎo)屈服加速度式(10)進(jìn)行計算，計算時為確保效率快，精度高。采用隨機搜索法與序列二次優(yōu)化結(jié)合進(jìn)行分析。

3 屈服加速度計算

本文針對不同的礦山邊坡進(jìn)行穩(wěn)定性分析，考慮上部荷載系數(shù)qt=0.1～0.3，坡腳β=30°～90°，內(nèi)摩擦角φ=10°～30°，c/γH=0.10～0.25，以及拉應(yīng)力區(qū)折減系數(shù)ξ=0～1的破壞工況。圖6～8中上部點線部分表示過坡趾下方的破壞，其余部分線條表示過坡趾的破壞。圖6中顯示，過坡趾破壞為摩擦體的特征，而當(dāng)內(nèi)摩擦角較小時，對于坡腳較為緩和的邊坡，則更容易發(fā)生過坡趾下方的破壞。但考慮地震荷載影響時，對于內(nèi)摩擦角更大的邊坡易發(fā)生過坡趾的破壞。由圖6～8可知，考慮張拉區(qū)域折減分析，邊坡的屈服加速度折減迅速，而針對內(nèi)摩擦角較小且坡度較大的邊坡，破壞更為強烈，而頂部荷載則會加速屈服加速度的折減。

圖6 qt=0.1時礦山邊坡屈服加速度Fig.6 Yield acceleration of mine slope under qt=0.1

圖7 qt=0.2時礦山邊坡屈服加速度Fig.7 Yield accelerations of mine slope under qt=0.2

圖8 qt=0.3時礦山邊坡屈服加速度Fig.8 Yield accelerations of mine slope under qt=0.3

圖9 考慮張拉折減下的礦山邊坡位移Fig.9 Displacement of mine slope under tension cut-off

4 工程案例分析

在設(shè)計中通常采用擬靜力來考慮邊坡上的地震荷載。為分析地震過程(加速度歷史)，對邊坡對影響，并深入了解不同頂部荷載系數(shù)、張拉區(qū)折減系數(shù)等因素的影響，本節(jié)將針對不同的地震記錄得出邊坡永久位移的計算結(jié)果。

在極限分析上限定理的框架下，CHANG等[13]提出了一種分析具有旋轉(zhuǎn)破壞面的邊坡永久位移對方法。本文采用該模型進(jìn)行分析，永久位移的計算過程簡要描述如下：當(dāng)?shù)卣鹣禂?shù)K超過屈服加速度Kc時，塊體開始滑落。通過平衡功率方程，可獲得旋轉(zhuǎn)角加速度，見式(12)。

(12)

式中：G為滑塊的重力(即圖5中的BDEC部分)；l為滑塊旋轉(zhuǎn)中心O到滑塊重心的距離。坡腳水平方向永久位移由式(12)的二重積分計算，見式(13)。

(13)

項目礦山位于廣西賀州市八步區(qū)里松鎮(zhèn)造紙廠東側(cè)約2 km一帶的山坡上，隸屬賀州市八步區(qū)里松鎮(zhèn)管轄。礦區(qū)距賀州市城區(qū)直距約24 km。賀州市城區(qū)至里松鎮(zhèn)有水泥公路相通，從里松鎮(zhèn)到礦區(qū)有簡易公路相通，路況較差，交通條件一般。

根據(jù)現(xiàn)場勘測結(jié)果，本文選取代表性的破壞工況為c/γH=0.15，坡腳β=55°，內(nèi)摩擦角φ=30°。對邊坡模型進(jìn)行穩(wěn)定性分析。2016年7月31日17時18分，廣西梧州市蒼梧縣發(fā)生了5.4級地震，震源深度達(dá)到了10 km，本文選取類似地震數(shù)據(jù)進(jìn)行輸入并深入討論其在地震力作用下產(chǎn)生的動態(tài)永久位移。

圖10 考慮張拉區(qū)域折減的礦山邊坡位移分布Fig.10 Distribution of permanent displacement ofmine slope considering tension cut-off

圖9展示了礦山邊坡在不同程度張拉折減系數(shù)下的永久位移。結(jié)果顯示，考慮張拉折減的影響，將對邊坡永久位移增加將近60%，如針對qt=0.3，邊坡在M-C準(zhǔn)則與考慮ξ=0時產(chǎn)生的位移分別為123 cm與195 cm，考慮強度折減后永久位移增加了58%。而圖10顯示，當(dāng)非線性系數(shù)ξ取0.4～1時，非線性準(zhǔn)則對于邊坡永久位移影響并不明顯;而當(dāng)ξ逼近0時，邊坡永久位移產(chǎn)生迅速增加，影響顯著。而結(jié)果曲率顯示，考慮頂部荷載作用對于邊坡永久位移的折減程度顯著大于無荷載作用邊坡。

5 結(jié) 論

1) 張拉應(yīng)力區(qū)域折減系數(shù)ξ對于屈服加速度Kc的影響顯著，考慮坡頂荷載作用則會增加該影響，因此在實踐中需要著重考慮邊坡后緣拉應(yīng)力區(qū)的影響，并進(jìn)行適度的加固或減少超載。

2) 永久位移的計算結(jié)果顯示，考慮張拉折減時，邊坡產(chǎn)生的永久位移將達(dá)到M-C準(zhǔn)則下近兩倍，且頂部系數(shù)的存在會顯著提高張拉強度折減帶來的不利影響。