考慮坡體后緣張拉應力折減的超載邊坡穩定性分析

王漢敬

(中國土木工程集團有限公司 北京 100038)

1 引言

地震是自然界最具破壞性的地質災害之一，不斷有學者研究地震對土工構筑物的影響[1－4]。其中，對于當前研究邊坡穩定性，通常采用線性的Mohr-Coulomb屈服準則(簡稱M-C準則)來評估其安全性，鮮有關于巖土體非線性強度特征及其對邊坡穩定性影響的相關研究。但土工試驗結果表明:采用線性的強度包絡線，即將土體的黏聚力和內摩擦角用以計算抗剪強度，進而外推得出土體的抗拉強度，但此強度并非是土體真實的抗拉強度[5]。另一方面，通過抗拉強度的研究發現[6]，采用兩種方法可以有效地分析邊坡拉應力區的抗拉強度，使分析結果更加趨于真實的破壞條件，其一是對巖土體的抗拉強度進行一定程度的折減，即將線性破壞曲線優化為非線性強度包絡線；其二則是在邊坡后緣的拉應力區引入不同形態的裂縫。對于后者考慮裂縫的存在，Michalowski等[7－8]將裂縫的擴展過程與整體的滑落狀態結合，研究結果表明:對于高陡邊坡或受靜水壓力影響顯著的邊坡，裂縫的存在大大增加了失穩的風險。Zhao等[9]分析了線性準則下不同大小的抗剪強度參數對帶裂縫邊坡的地震穩定性的影響，研究顯示地震荷載對邊坡的穩定性影響尤為顯著。Baker等[10－11]對拉剪區域進行了一定的折減，并結合變分法對強度包絡線中的張拉區域進行相關的穩定性分析。康孝森等[12]基于非飽和土力學和最大伸長線應變理論推導出入滲影響下黃土邊坡坡頂張拉裂縫極限深度控制方程。

實際上，超載邊坡是廣泛存在的，且目前關于超載邊坡尚未有考慮土體拉應力區強度折減后的超載邊坡動態全過程分析的相關研究。為研究抗拉強度折減對超載邊坡在地震作用下穩定性的影響，本文在一定研究成果的基礎上，首先通過極限分析上限法對超載邊坡臨界的破壞形態以及功能方程進行推導；然后對計算求解算法進行深度優化，結合隨機搜索法(Monte Carlo Method)與二次序列規劃程序(SQP)進行求解；最后針對真實地震中邊坡破壞案例，選取實測地震動時程，分析不同超載系數以及張拉強度折減系數對邊坡穩定性的影響。

2 土體強度包絡線

對邊坡穩定性進行分析時，土體一般采用理想均質材料，而M-C準則并不能精確反映土體的抗拉強度。大量的試驗結果表明，非線性的強度破壞包絡線能更真實地反映巖土體的抗拉強度。

經典的M-C準則在主應力空間可以表示為一個不規則六邊形截面的錐體，如圖1a所示。后來學者對經典的M-C屈服面進行多次修改，從而限制了抗拉強度，其屈服面如圖1b所示。折減后的部分拉力截斷單軸抗拉強度ft可表示為圖2a中的應力圓C1。由圖2中的莫爾圓轉換至圖1中的主應力空間是一個較為復雜的過程，圖1b中截面ABCDEF上AB邊上的應力狀態由圖2b中的應力圓C2表示，而AB上各個位置由其中主應力的大小確定(σ1≥σ2≥σ3)。例如，圖2中的應力圓C1在σ2＝0的情況下對應于圖1b中ABCDEF截面與σ3軸的交點。極限狀態時的屈服面由圖2中光滑的屈服面映射至圖1中具有棱角的屈服面。采用此拉應力區的截斷方式適用于極限分析中的解析運動學分析方法，但因為主應力空間中屈服面棱角處的法向量不唯一，其在數值方法中的適用性還需深入研究。傳統的M-C強度包絡線如圖2a所示，M-C函數為簡單的直線與圓弧部分的組合。在M-C計算中，三軸抗拉強度的大小為f3t＝c/tanφ，而單軸抗拉強度大小為ft＝2ccosφ/(1＋sinφ)。研究中當考慮抗拉強度的折減時，牽引矢量從直線變化為圓弧，且張拉應力區域的剪脹角δ從φ變化到π/2，當δ＝φ時代表線性的M-C屈服準則。為了量化展示折減后的抗拉強度，本文引入了抗拉強度折減系數ξ:

其中，ξ＝0～1。當ξ＝1時為線性M-C準則下的抗拉強度；當ξ＝0時表示完全不考慮抗拉強度時的破壞狀態，如圖2b所示。采用此拉伸截斷方式令極限分析的解析運動學分析過程大大簡化。

圖1 屈服準則模型

圖2 不同包絡線示意

3 破壞機制與計算分析

3.1 破壞機制分析

圖3為受地震荷載作用的臨界均質超載邊坡模型。塊體BDEC以中心O為中心沿著滑裂面DEC旋轉，而整個破壞為過坡趾的破壞，DEC部分為臨界滑面。研究發現，張拉應力更容易出現在坡頂附近，而滑面上部DE部分由受拉伸區折減部分強度包絡線而引起，下方較低的部分(EC)是由線性部分的包絡線引起。DE部分剪脹角δ從D點的最大值δmax向E點的φ不斷折減。圖4則展示了DE面上的一個微元，可得:

沿DE積分可以得到:

圖3 受地震荷載作用的臨界均質超載邊坡模型

圖4 邊坡上滑面計算單元

圖3～圖4中，r0和θ0分別為對數螺旋線起始點A到旋轉中心半徑的距離AO和AO與水平參考線的夾角(見圖2)；θh為對數螺旋線終止點到旋轉中心的半徑與水平參考線的夾角；θm為OD與水平參考線的夾角；θtc為OE與水平參考線的夾角；δ為剪脹角，其詳細分布形式可表示為:

極限分析上限法基于塊體內外能量耗散率的平衡，功能方程為:

式中，Wr為土體自重狀態下的能量耗散率；Ws為擬靜力地震荷載作用下的能量耗散率；Wqt為超載作用的功率；D為滑體BDEC的內能耗散率。外部功率Wr和Wr為塊體ABC的功率減去塊體ADE的功率，分別表示為:

公式中帶有下標的參數分別表示圖3中區域OCA、AOB、BOC、OEA、AOD、DOE 的能量耗散功率；fw與fs均為包含變量β、θ0、θtc、θm、θh的函數方程。分析超載部分功率時，引入超載系數qt，fq為超載函數。根據主應變率，由張拉區與折減產生的內部能量耗散為，則整個滑裂面的做功可表示為:

式中，fd為包含變量θtc、θh、φ的函數；fdc為包含變量θtc、θm、φ的函數。將上述公式整理，即可得到屈服加速度Kc的表達式:

Kc為一個含θ0、θtc、θm、θh以及c/γH、坡腳β和內摩擦角φ等參數的函數，則屈服加速度Kc的獲取為一系列變量組合下的最小值。

3.2 計算優化

為搜索Kc在θ0、θtc、θm、θh、c/γH、β、φ等參數下的臨界值，需對Kc進行求導分析，具體限定以及計算式為:0＜θ0＜θm＜θtc＜θh＜π

在限定的基礎上，進行屈服加速度的求解。但采取傳統的計算方式存在一些缺陷，如遍歷法計算速度慢、隨機搜索法計算精度低以及優化分析易于陷入局部最優解等。本文采用隨機搜索法和非序列二次優化法(SQP)結合進行分析求解，隨機搜索法精度雖然較低但是可以很好地給出初值的區間范圍，而得到較為合理的初值后可采用序列二次優化進行計算臨界值搜索。部分驗證結果見表1～表2。其中表1為參數設計，表2為與經典解的對照。可以看出計算結果與經典解非常吻合。

表1 邊坡參數設計

表2 本文計算結果與經典解對比結果

4 屈服加速度計算及地震破壞案例分析

4.1 屈服加速度計算

本文針對不同程度的超載邊坡進行穩定性分析，考慮超載系數qt＝0 ～0.3、坡腳β＝30°～90°、內摩擦角φ＝30°、c/γH＝0.1～0.25，以及拉應力區折減系數ξ＝0～1的破壞工況。圖3中上部點線部分表示過坡趾下方的破壞(Below-toe failure)；其余部分線條表示過坡趾的破壞(Toe failure)；圖中原點為過渡點。由圖3可知，過坡趾破壞為摩擦土體的特征，而當內摩擦角較小時，對于坡腳較為緩和的邊坡，則更容易發生過坡趾下方的破壞。但考慮地震荷載影響時，對于內摩擦角更大的邊坡易發生過坡趾的破壞。而由圖5可知，考慮張拉區域折減分析，邊坡的屈服加速度折減迅速，而針對內摩擦角較小且坡度較大的邊坡，破壞更為強烈。超載作用則會加速屈服加速度的折減。圖中ξ＝1表示線性M-C屈服分析結果。

圖5 不同超載系數下的邊坡屈服加速度

4.2 地震破壞案例分析

選取2018年5月12日的里氏8級汶川地震為研究對象，記錄臺站為512MZQ，加速度放大系數為0.868，采樣間隔為0.005 s。為進行全面分析，本文選取破壞工況為c/γH＝0.15、β＝55°、φ＝30°。圖6所示為超載邊坡在不同程度張拉折減系數下的永久位移。結果顯示，考慮張拉折減的影響，將對邊坡永久位移增加2倍甚至2倍以上，如qt＝0.3，邊坡在M-C準則與考慮ξ＝0時產生的位移分別為123 cm與195 cm，考慮強度折減后永久位移增加了59%。而圖7顯示，當非線性系數ξ取0.4～1.0時，非線性準則對于邊坡永久位移影響并不明顯；而當ξ逼近0時，邊坡永久位移迅速增加，影響顯著。而結果曲率顯示，考慮超載作用對于邊坡永久位移的折減程度顯著大于非超載邊坡。

圖6 汶川地震作用下考慮張拉折減邊坡永久時程位移分布

圖7 汶川地震波作用下不同張拉折減系數邊坡永久位移對比

5 結論

(1)采用Monte-Carlo與SQP綜合優化法，顯著提高運算效率，且本文計算結果與經典解進行對比驗證相吻合。

(2)張拉應力區域折減系數ξ對于臨界加速度Kc的影響顯著，隨著超載作用不斷增加，影響在不斷放大，而坡度較大、內摩擦角較小的脆性邊坡更容易受到破壞。

(3)考慮張拉折減與超載共同作用時，邊坡產生的永久位移將超過傳統M-C準則下兩倍，且當抗拉強度折減程度不斷增加時，邊坡位移積累增加迅速，而超載作用則會加速位移的上升。