國家天文臺500米口徑球面射電望遠鏡臺址球冠型邊坡穩定性分析

陳德茂，沈志平，姜 鵬，付君宜，劉 慧

（1.貴州正業工程技術投資有限公司，貴州 貴陽 550012；2.中國科學院國家天文臺，北京 100012）

山區工程建設不可避免地會遇到邊坡問題。邊坡穩定性評價方法目前主要包括極限平衡法、極限分析法、數值分析法三類。瑞典圓弧法作為極限平衡法的一種最早由Fellenius于1927年提出，Bishop對瑞典圓弧法進行了改進并提出了滿足力矩平衡的簡化Bishop法[1?2]；針對Bishop法不滿足力平衡條件的問題，Janbu提出了滿足力平衡條件的Janbu法[3]；鄭穎人等[4]分析了幾種極限平衡法的精度問題；麻官亮等[5]將Spencer法應用到三維邊坡模型的穩定性評價中。極限分析法由Drucker和Prager于1952年提出，該方法考慮了靜力場和速度場，根據滿足的條件不同分為上限解法（能量法）和下限解法[6?8]，當上限解和下限解相等時即為邊坡穩定性的真實解[9]。數值分析法于20世紀60年代提出，目前主要包括有限元法、離散元法、有限差分法、邊界元法。基于不同的理論，各數值分析法均有各自的優勢和局限性[10?13]。工程中最常用的數值分析方法是有限元法，能夠較好地處理小變形問題，配合強度折減法[14?16]可以獲得邊坡安全系數。

目前極限平衡法依然是工程中常用的邊坡穩定性分析方法，該類方法基于平面應變假定而未考慮邊坡外形，在評價球冠型邊坡這類有明顯外形影響的邊坡穩定性時有較大誤差，本文研究旨在消除該誤差問題。

1 FAST臺址區域概況

中國科學院國家天文臺500米口徑球面射電望遠鏡（Five-hundred-meter Aperture Spherical radio Telescope,FAST）位于貴州省平塘縣大窩凼洼地內。大窩凼洼地發育在云貴高原向廣西丘陵過渡帶的大小井地下河系統中，為“U”字型巖溶洼地，地層巖性由黏土、溶塌混合體、三疊系中統涼水井組石灰巖組成。黏土僅分布于洼地底部，厚度一般3～5 m；溶塌混合體分布于斜坡地帶及洼地底部黏土層之下，厚度一般5～50 m，平均30 m；淺表為未膠結的大塊石，中下部為半膠結、膠結的土石混合體，均勻性差；三疊系中統涼水井組白云質灰巖、含泥灰巖出露于洼地上部陡坡地段，中至厚層狀，屬較軟至較硬巖，巖體較破碎至較完整。

為契合FAST主動反射面形狀，臺址開挖完成后在反射面圈梁以下范圍內形成了大規模球冠型邊坡[17]，如圖1所示。由于淺表層未膠結的溶塌混合體極易發生破壞[18]，開挖時清除了大量該類土體，開挖完成后的球冠型邊坡主要由膠結的溶塌混合體和下伏基巖組成，局部穩定性良好，但球冠型邊坡整體穩定性評價缺少相關計算方法。

圖1 FAST臺址球冠型邊坡Fig.1 Spherical cap shape slope of FAST

考察邊坡在水平面內的形狀，可將其分為凸形、凹形和直線形，邊坡的空間形狀對其穩定性無疑會產生影響[19?21]，球冠型邊坡屬于軸對稱圓形凹坡，實際工程中通常不考慮邊坡外形，均按照直線形邊坡考慮，所得計算結果偏于保守[22?23]，需重新確定其穩定性評價方法。

目前一般認為軸對稱圓形凹形邊坡穩定性優于直線形邊坡，原因是軸對稱圓形凹坡的微單元扇形滑體側面較直線形邊坡多了一組側面力，Zhang[24]和Zhang等[25]都將該側面力按主動土壓力取值，從概念上分析軸對稱圓形凹坡的微單元扇形滑體在滑動時側面處于擠壓狀態，更接近于被動土壓力。黎莉等[26]對側面力分別按照主動土壓力、靜止土壓力、被動土壓力三種情況進行分析，結果表明側面實際壓力狀態處于靜止土壓力和被動土壓力之間。

2 軸對稱圓形凹坡的微單元扇形滑體側面力

2.1 環形條塊軸向力與外荷載的關系

根據極限平衡法劃分條塊的思路，將軸對稱圓形凹坡滑體劃分為若干個“環形條塊”（圖2），環形條塊向下滑動時會受到一個水平分力均布荷載qx，荷載方向指向對稱軸（圖3）。

圖2 軸對稱圓形凹坡劃分環形條塊Fig.2 Dividing annulus slices of the axisymmetric circular concave slope

圖3 環形條塊水平下滑分力Fig.3 Horizontal sliding component of annulus slices

取圖3中環形條塊的一段微單元扇形滑體作為分析對象（圖4）。微單元扇形滑體的曲率半徑為環形條塊的半徑R，兩側面夾角為dφ，用r和s分別表示微單元扇形滑體的法線和切線方向，則微單元扇形滑體弧長ds=Rdφ。s和r方向的荷載集度分別為qs和qr。側面彎矩為M，側面剪力為FQ，側面軸向力為FN。

圖4 微單元扇形滑體受力圖Fig.4 Free-body diagram of the tiny element sector sliding body

由ΣFs=0得

由于dφ趨向于0，則cos（dφ/2）=1，sin（dφ/2）=dφ/2，忽略高階微量后由式（1）得

令ds=Rdφ，由式（2）整理得

由ΣFr=0同理可得

由ΣM=0得

即

由式（3）（4）（5）可知，當R趨向于無窮時，環形條塊內力公式轉變為直線形條塊。

由于環形條塊滑動時受到的水平分力qx為均布荷載，故切線荷載qs=0，法向荷載qr=qx。

設環形條塊處于無彎矩狀態，即側面彎矩M=0，側面剪力FQ=0，由式（4）可得

2.2 微單元扇形滑體側面力與水平抗滑力的關系

由式（6）可知，當環形條塊向下滑動時會在環形條塊內產生一個軸向壓力FN，FN即為微單元扇形滑體的側面力，且FN為常數。如果可以確定環形條塊軸向抗壓能力FNmax，則可以得到一個水平荷載，該水平荷載方向背離對稱軸，起到抗滑作用，屬于抗滑力。

如圖4所示，假定極限平衡狀態時，微單元扇形滑體的大主應力σ1方向為s切線方向，小主應力σ3方向為重力方向（豎直方向），則有：

考慮到邊坡坡面臨空，假設小主應力σ3=0，則有：

式中：h—微單元扇形滑體的高度；

l—微單元扇形滑體的底面寬度；

θ—微單元扇形滑體的底面與水平面夾角；

c1—微單元扇形滑體的滑體黏聚力；

φ1—微單元扇形滑體的滑體內摩擦角。

對于一般土體，式（8）中的σ1介于靜止土壓力和被動土壓力之間，符合文獻[26]的結論。

取微單元扇形滑體（圖5），由式（6）可得FNmax產生的水平抗滑力（T）的計算公式：

圖5 微單元扇形滑體抗滑力分析Fig.5 Sliding resistance force of the tiny element sector sliding body

計算大主應力σ1時由于假定σ3=0，對于坡面處的滑體部分較為合理，但對于深處的滑體部分其σ3不為0，故式（10）計算出的水平抗滑力偏小，使得最終穩定性評價結果偏于保守。

3 基于極限平衡法的軸對稱圓形凹坡穩定性評方法

3.1 改進簡化Bishop法

將2.2節中介紹的水平抗滑力加入到任意極限平衡法中，均可以使其適用于軸對稱圓形凹坡穩定性評價，這里以簡化Bishop法為例進行公式推導。

取第i個環形條塊的微單元扇形滑體i，將水平抗滑力Ti引入到簡化Bishop法中，假定水平抗滑力Ti作用在微單元扇形滑體i的重心上，微單元扇形滑體i重心所在剖面的計算簡圖如圖6所示。

圖6 改進簡化Bishop法計算簡圖Fig.6 Calculation diagram of the improved simplified Bishop method

由豎向合力ΣFz=0得

式中：γi—微單元扇形滑體i的土體容重；

Ti0、Ni—微單元扇形滑體i在滑面上的抗滑力和法向力；

ri—微單元扇形體重心到對稱軸的距離；

其余變量意義同前。

由力矩求和ΣMO=0得：

式中：Ri—微單元扇形滑體i的水平抗滑力Ti到圓弧滑面圓心O的力矩；

R—圓弧滑面半徑。

由摩爾庫倫強度準則并引入安全系數Fs可得Ti0：

式中：c2i—微單元扇形滑體i的滑面黏聚力；

φ2i—微單元扇形滑體i的滑面內摩擦角；

其余變量意義同前。

將式（11）代入式（13）整理得：

將式（14）代入式（12），約去dφ整理得：

3.2 改進簡化Bishop法驗證

3.2.1 簡化Bishop法和強度折減法對比

首先對有限元強度折減法和簡化Bishop法在計算直線形邊坡時安全系數的差異進行對比：簡化Bishop法采用Geostudio，自動搜索滑移面；有限元強度折減法采用Midas GTS NX。為控制有限元模型中滑移面為圓弧形，模型分為基巖和土體兩部分，邊坡高10 m，坡度42.5°，土體參數如表1所示，所有土體均采用摩爾庫倫強度準則，安全系數及滑移面計算結果如圖7所示。

表1 10 m高度邊坡土體參數Table 1 Soil parameters of slope with a height of 10 m

圖7 簡化Bishop法與強度折減法安全系數計算結果Fig.7 Safety factor of the simplified Bishop method and strength reduction method

由圖7可以看出兩種方法的滑移面基本一致，但有限元強度折減法的安全系數計算結果略高于簡化Bishop法[27]，這主要是因為簡化Bishop法有諸多假定，如未考慮土條間切向力、未考慮單個土條的力矩平衡等，這些假定使得計算結果偏于保守。

3.2.2 改進簡化Bishop法和強度折減法對比

建立圖7（b）模型的三維有限元模型（圖8），坡腳到對稱軸的距離Rbottom分別為5，10，20 m。計算得到安全系數和滑移面后，將滑移面提取并采用式（15）—（17）進行安全系數計算，并將有限元強度折減法計算結果和式（15）—（17）的計算結果進行對比，結果見圖9（a）。圖中H為邊坡高度，H/Rbottom=0表示直線形邊坡；同時將表1土體參數的黏聚力改為10 kPa，其余參數不變，進行安全系數計算對比，計算對比結果如圖9（b）所示。

圖8 三維有限元模型（Rbottom=20 m）Fig.8 Three-dimensional finite element model（Rbottom=20 m）

圖9 10 m高軸對稱圓形凹坡安全系數計算結果Fig.9 Safety factor of the axisymmetric circular concave slope with a height of 10 m

另建立兩組模型進行對比驗證，邊坡模型高度為30 m，邊坡坡度60°，土體參數如表2所示，安全系數計算結果如圖10（a）所示。同時將表2土體參數的黏聚力改為35 kPa其余參數不變，進行安全系數計算，計算結果如圖10（b）所示。

表2 30 m高度邊坡土體參數Table 2 Soil parameters of slope with a height of 30 m

由圖9和圖10可得，有限元強度折減法安全系數計算結果始終高于改進簡化Bishop法且兩條計算結果線段接近平行，表明兩種計算方法的安全系數差值基本保持一致，改進簡化Bishop法的計算公式較為可靠，同時可以看出軸對稱圓形凹坡的穩定性優于直線形邊坡。

圖10 30 m高軸對稱圓形凹坡安全系數計算結果Fig.10 Safety factor of the axisymmetric circular concave slope with a height of 30 m

4 軸對稱圓形凸坡的穩定性評價

軸對稱圓形凸坡的安全系數計算公式可由式（15）—（17）變換得到，當軸對稱圓形凸坡的環形條塊下滑時其軸向會產生一個拉力，一般不考慮土體的抗拉強度，則忽略式（17）得到適用于軸對稱圓形凸坡的安全系數計算公式：

對式（18）（19）進行驗證，同樣采用有限元強度折減法和改進簡化Bishop法計算得到的安全系數結果進行對比。建立有限元模型，土體單元采用表2中的土體參數，邊坡高度20 m，坡度為70°，Rbottom=20，30，40 m三種情況（圖11），兩種方法計算結果如圖12所示。

圖11 軸對稱圓形凸坡三維有限元模型Fig.11 Three-dimensional finite element model of the axis ymmetric circular convex slope

圖12 軸對稱圓形凸坡安全系數計算結果Fig.12 Safety factor of the axisymmetric circular convex slope

圖12中兩種方法安全系數計算結果比較接近，式（18）（19）適用于軸對稱圓形凸坡安全系數計算，同時軸對稱圓形凸坡的安全系數隨半徑的增大沒有發生明顯變化，與傳統概念中凸坡的穩定性差有所出入。這主要是因為上部環形條塊周長要小于下部環形條塊，即上部滑體的體積較小而下部滑體的體積大，這種情況對抗滑移是有利的，凸坡安全系數未必就比直線形邊坡小[22]。

5 結論

（1）軸對稱圓形凹坡側面力取值可以按照微單元扇形滑體小主應力為零時的大主應力進行取值。

（2）考慮側面力取值后對極限平衡法進行改進，提出改進后的簡化Bishop法安全系數計算公式，改進后的公式計算結果與數值分析結果基本一致。

（3）假定土體抗拉強度為零，提出了圓形凸坡安全系數計算方法與公式，計算結果表明圓形凸坡安全系數與長直邊坡接近。