基于擬動力法考慮土體張拉截斷的邊坡抗震穩定性上限分析

梁承龍,劉 芳

(廣西交通職業技術學院 土木建筑工程學院,南寧 530012)

0 引言

近年來,邊坡失穩造成的滑坡事件屢見不鮮,其破壞性強以及發生頻率高,直接威脅著周邊人民的人身安全和正常生活,滑坡事件已成為災害治理的重點[1-3]。目前,極限平衡法、有限元法以及極限分析法是國內外學者們解決邊坡失穩問題最常見的數值和解析方法,而極限分析法相較于其他兩種常規方法有著既簡易又合理的優點,分析問題時不必考慮土體內部復雜的應力狀態,只需構建一個運動學許可的速度場便可求解極限荷載的上限解答,進而確定邊坡是否穩定。在眾多巖土體本構關系中,Mohr-Coulomb(M-C)屈服準則廣受工程師和學者們的青睞,其包含兩個抗剪強度指標:粘聚力c和內摩擦角φ。基于M-C屈服準則,CHEN[4]證明了在二維情況下邊坡最危險滑動面是一條對數螺旋線。學者們采用該對數螺旋破壞機制進行了大量邊坡穩定性分析[5-7]。戴自航等[8-9]通過分析邊坡土體內部的應力狀態得出在邊坡失穩時往往坡頂后緣下一定深度范圍內土體受到拉應力作用,導致巖土體受拉破壞,因此將整個滑動面上都視為剪切破壞所得結果偏于保守且與實際土體破壞狀態不相符,邊坡失穩時實際上是張拉-剪切復合破壞,并運用有限元法結合案例證明了考慮土體張拉破壞的必要性。因此,基于張拉-剪切破壞模型的邊坡穩定性分析不斷涌現[10-13]。例如,王偉等[14-15]采用應變軟化模型考慮邊坡失穩時張拉-剪切漸進破壞過程并計算邊坡安全系數以及搜索臨界滑動面。除此之外,DUNCAN等[16]提出兩種方法考慮邊坡后緣的受拉破壞:其一,在邊坡坡頂處引入一條垂直張拉裂縫;其二,對受拉區的巖土體抗拉強度進行折減,即張拉截斷。方法一表明滑動面由對數螺旋線和垂直張拉裂縫組成以消除邊坡后緣的受拉區。UTILI[17]綜合分析了裂縫的深度和位置對邊坡穩定性的影響;MICHALOWSKI[18]考慮了邊坡失穩時裂縫的形成過程對邊坡穩定性的影響;HE等[19]將二維情況下裂縫邊坡穩定性分析拓展到三維情況;饒平平等[20]對裂縫邊坡的三維破壞模式進行了拓展;基于擬靜力法,HE等[21]和LI等[22]對裂縫邊坡的三維抗震穩定性進行研究,繪制出了穩定性圖表;基于擬動力法,HOU等[23]基于離散運動學機構研究非均質裂縫邊坡的穩定性,分析了擬動力法參數對安全系數的影響規律,探討了擬動力法與擬靜力法的差異,表明擬動力法更符合實際,因其能考慮地震隨時空變化的動力特性和地震波在巖土層內傳播時土層振動相位差和振幅加速度放大效應。方法二是張拉截斷,即將土體的抗拉強度從M-C屈服準則中削弱或去除。DRUCKER等[24]首次將張拉截斷引入邊坡穩定性分析中;MICHALOWSKI[25]分析了土體張拉截斷對邊坡穩定性的影響規律并得出孔隙水壓力作用下完全張拉截斷會導致陡坡的穩定系數將下降至少50%;PARK等[26]分析了考慮張拉截斷的巖質邊坡三維穩定性問題;ZHANG等[27]對張拉裂縫和張拉截斷所得結果進行對比分析。基于擬靜力法,學者們對地震作用下考慮張拉截斷的邊坡臨界地震加速度進行求解[28-31],在此基礎上結合工程實例,采用Newmark法計算地震效應下邊坡的永久位移,得出考慮張拉截斷的M-C屈服準則計算所得邊坡永久位移是采用線性M-C屈服準則的計算結果2倍之多。然而,考慮張拉截斷的邊坡穩定性分析大多基于擬靜力法,未能考慮地震隨時間和空間變化的動力特征。

鑒于此,本文考慮完全張拉截斷以去除土體的抗拉強度,并基于擬動力法對邊坡進行抗震穩定性分析。在極限分析理論框架下,采用“點到點”的方法構建離散運動學破壞機構[32-33],并建立能量平衡方程,結合強度折減法對地震作用下邊坡的安全系數進行求解。通過二分法和優化算法計算所得結果與擬靜力法進行對比分析,并對比張拉截斷和張拉裂縫兩種方式下的邊坡安全系數和臨界滑動面,探討擬動力法參數和張拉截斷效應對邊坡穩定性的影響規律,以期為實際工程提供一定的參考。

1 土體屈服準則

邊坡穩定性分析中認為邊坡失穩時滑動面上任意一點土體所受到的剪應力超過其抗剪強度而發生剪切破壞。實際邊坡失穩破壞時往往在坡頂后緣一定深度范圍內會發生張拉破壞[8],而M-C屈服準則仍將這一部分張拉破壞視為剪切破壞,所得到的滑動面為剪切滑動面,未能真實地反映邊坡的失穩狀態。同時,土體的抗拉強度并不是由試驗直接測量得出的,而是根據M-C屈服準則由土體壓縮狀態向拉伸狀態外推得出的結果,遠大于土體的實際抗拉強度,不能真實地反映出土體的抗拉強度,因此在邊坡穩定性分析中需要對所得的抗拉強度進行折減[25]。如圖1所示,根據M-C屈服準則,可以得到土體三軸抗拉強度f3t=c/tanφ和單軸抗拉強度ft=2ccosφ/(1+sinφ),并引入張拉截斷折減系數ξ對土體抗拉強度進行折減[25],此時,M-C屈服準則由直線和一段圓弧組成,剪脹角δ則從φ變化到90°,折減后土體抗拉強度可表示為:

圖1 土體張拉截斷包絡線Fig. 1 Cut off envelope of soil mass by tension

(1)

式中: 0≤ξ≤1,當ξ=1時表示線性M-C屈服準則,當ξ=0時表示完全張拉截斷,即土體的抗拉強度為0,也是本文研究的重點。

以下對剪脹角δ進行說明。極限分析定理要求土體屈服時滿足相關聯流動法則,即各個點的塑性應變增量方向為該點塑性勢面的外法線方向,而張拉截斷概念的引入會造成屈服面可能有拐角或尖角,導致法線方向不是唯一確定的,但對極限分析并不造成障礙[7]。由于假設土體滿足M-C屈服準則,土體發生剪切破壞時必定伴隨著體積膨脹,根據幾何關系可得塑性應變增量方向與滑動面成角度φ,即認為土體剪脹角δ=φ,此時δ大于土體實際剪脹角,從而高估土體的體積應變,但所得結果仍是一個嚴格的上限解答。而引入張拉截斷的概念會導致δ由φ變化到90°,這是由于土體發生破壞時需滿足相關聯流動法則所造成的結果,意味著運用極限分析允許發生較大的體積應變,土體發生斷裂破壞,區別于塑性屈服,計算所得結果仍是真實極限荷載的上限解答[34-35]。PUAL[36]認為考慮張拉截斷的M-C屈服準則是描述土體發生斷裂的狀態而不是土體發生屈服的狀態。因此,在張拉截斷部分將δ理解為破裂角以區分剪脹角的概念。

2 邊坡安全系數的上限解答

極限分析上限定理可描述為對于任意一個滿足運動學許可的速度場,從能量平衡角度建立平衡方程,可以獲得極限荷載的上限解答。因此,根據所建立的運動學許可的破壞機構,計算邊坡發生失穩瞬間的外力做功功率和內能耗損率,通過建立能量平衡方程可以得到邊坡安全系數的上限解答。在本文中,邊坡土重和地震力做功為外力做功功率;內能耗散率只考慮邊坡潛在滑動面瞬時滑動所產生的能量。

2.1 離散機構的構建

在二維情況下,傳統上限法假設滑動面為對數螺旋線滑動面,進而再求解邊坡穩定系數或安全系數。不同于傳統上限法,本文不再假設滑動面形狀,通過“點到點”方法遞推出一系列離散點相連逐步構成滑動面[32-33],因此滑動面不再是光滑連續曲線,而是許多微小的直線段的組合,其優勢在于這種構成滑動面的方法能很好地處理土體內摩擦角非均質性對滑動面的影響。如圖2所示,以邊坡坡趾A為坐標原點建立直角坐標系,邊坡高度為H,坡角為β。本文假設滑動面經過A,即邊坡發生坡趾破壞,在此基礎上,建立邊坡離散運動學破壞機構,并具體給出“點到點”生成滑動面的推導過程。

圖2 離散運動學許可破壞機構[25]Fig. 2 Discrete kinematics Permit Destruction Mechanism

邊坡失穩時滑動土體視為剛體且滿足考慮土體完全張拉截斷的M-C屈服準則,并以角速度ω繞點O作旋轉破壞。根據所建立的直角坐標系,旋轉中心點O的坐標可由變量r0和θ0確定:

(2)

(3)

由于本文考慮土體抗拉強度完全折減,因此引入變量Z來描述完全張拉截斷部分滑動面的深度,深度為H-Z,如圖2所示。完全張拉截斷部分的破裂角δ采用與離散點的縱坐標呈線性關系表示[25],因此,整個滑動面的破裂角δi分布可表示為:

(4)

式中:δm為坡頂C處的破裂角,φZ為高度為Z處的土體內摩擦角。當yi≤Z,表示土體發生剪切破壞,當yi>Z時,采取線性插值方法調整yi=Z,如圖1中點D所示;以D作為起始點,土體發生斷裂,當yi>H時,同樣采取線性插值方法調整yi=H,結束滑動面的生成。

2.2 外功率

外功率包括滑動土體自重做功功率和地震力做功功率。計算簡圖如圖3所示,其中:c0和φ0分別為坡頂處粘聚力和內摩擦角,ch和φh分別為坡趾處粘聚力和內摩擦角。相較于傳統對數螺旋破壞機構,本文所構建出的“點到點”離散破壞機構的優勢在于分析內摩擦角非均質性的邊坡穩定性,例如成層邊坡或者內摩擦角隨深度變化的邊坡。因此,本文針對均質邊坡和非均質邊坡進行穩定性分析與對比。其中:均質邊坡定義為邊坡橫截面上粘聚力和內摩擦角保持不變,即c0=ch,φ0=φh;非均質邊坡定義為邊坡橫截面上粘聚力和內摩擦角隨深度線性增加。

圖3 計算簡圖Fig. 3 Calculation diagram

2.2.1 土體自重做功功率

對滑動土體ABC分割成若干個梯形塊體,即對于滑動面AC上相鄰兩個離散點Pi和Pi+1作平行線交邊坡坡面于點Qi和Qi+1從而形成梯形塊體PiPi+1Qi+1Qi。因此,滑動土體自重功率Wγ可用離散梯形塊體的自重功率之和:

Wγ=ω∑γiAi(xGi-xO)

(5)

式中:γi為第i個梯形塊體重心處的土體重度,Ai為第i個梯形塊體面積,xGi為第i個梯形塊體重心的橫坐標。

2.2.2 地震力做功功率

由于擬靜力法忽略了地震力隨時間變化及地震波在土體中傳播的特性,本文采用擬動力法來描述地震力的動力特征[23]。擬動力法假定土體的剪切模量為常數,水平和豎直地震加速度分別與地震剪切和縱波波速有關,水平和豎直地震加速度簡化為正弦函數,其大小會隨著時間呈周期變化,可以較好地體現地震波隨時間循環變化的規律,同時,擬動力法還考慮了土體對地震波的放大效應,因此擬動力法能更好地描述邊坡的動態穩定性[37-40]。因此,t時刻第i個梯形塊體的水平地震加速度系數kh和豎向地震加速度系數kv分別表示為:

(6)

(7)

式中:ah和av分別為水平和豎向地震加速度,VS和VP分別為地震剪切和縱波波速,f為土體放大系數,T為地震波周期。文獻[41]表明:計算過程中豎向地震加速度可作為水平地震加速度的一個分量表示,因此,kv=λkh,λ為豎向與水平地震加速度系數比例系數,本文λ設置為0.5。與滑動土體自重功率推導過程相同,水平和豎向地震力做功功率可表示為:

Weh=ω∑γiAi·kh·(yO-yGi)

(8)

Wev=ω∑γiAi·kv·(xGi-xO)

(9)

式中:Weh和Wev分別為水平和豎向地震力功率,yGi為第i個梯形塊體重心的縱坐標。

2.3 內能耗散率

由于滑動土體沿滑動面做剛體轉動,因此內能耗散率只發生在滑動面AC上,可表示為[25,29]:

(10)

式中:ci為點Pi處的粘聚力,Li為滑動面PiPi+1的長度,Ri為OGi的長度。

2.4 安全系數

根據能量平衡方程,令外功率等于內能耗散率,如式(11)所示,結合強度折減法計算邊坡的安全系數,通過土體的實際抗剪強度參數c和tanφ分別除以折減系數FS,可以得到土體處于極限平衡狀態下的抗剪強度參數c′和tanφ′,如式(12)所示,并定義FS為邊坡安全系數。

Wγ+Weh+Wev=D

(11)

(12)

式(11)是關于FS的非線性隱函數,變量為θ0、r0、δm、Z和t。利用MATLAB編寫迭代程序,通過二分法和優化算法可以快速精確地得出邊坡安全系數FS的值,確定滑動面的位置,具體優化過程詳見文獻[32]。

3 對比驗證

為了證明本文公式推導和計算結果的正確性和有效性,與文獻[25,33]中所得結果進行對比驗證。首先,對不考慮土體張拉截斷情況下的邊坡穩定性進行對比。圖4所示為基于擬動力法不同f,T和λ對臨界水平地震加速度khc的影響對比,khc的表達式可結合式(5)-式(11)得到:

圖4 邊坡臨界水平地震加速度對比

(13)

邊坡相關參數為:β=40°,H=10 m,γ0=14 kN/m3,γh=22 kN/m3,c0=15 kPa,ch=30 kPa,φ0=10°,φh=20°,λ=0.5,T=0.3 s,t0=0,f=1.5,VS=150 m/s,VP=280.5 m/s。其中:f的范圍為1.0～1.8,T的范圍為0.2～0.4 s,λ的范圍為-1.0～1.0,并按式(14)進行歸一化處理為無量綱形式:

(14)

式中:x′為歸一化結果,x為輸入參數,xmin為參數最小值;xmax為參數最大值。圖4表明本文計算所得khc與文獻[33]結果吻合度很高,且khc隨著f、T和λ的增大而減小,進一步說明地震周期T和土體放大系數f在邊坡抗震穩定性分析中的重要性。表1所示為在靜力狀態時考慮土體張拉截斷情況下的邊坡臨界高度Hcr的對比,文獻[24]通過構建傳統對數螺旋破壞機構求解邊坡穩定系數γH/c,為了便于計算結果比較,令c/γ=1.0,將文獻[24]結果轉化為邊坡臨界高度Hcr。邊坡相關參數為:γ=20 kN/m3,c0=ch=20 kPa,φ0=φh=20°,kh=0.0。由表1可知: 本文所得的結果與文獻[25]結果十分接近且略小于文獻[25]結果,最大誤差不超過3.1%,主要原因在于破壞機構構建的形式不同導致在土體張拉截斷區域對破裂角δ所采用的分布不同。文獻[24]認為張拉截斷區域的δ隨角度呈線性分布(文獻[24]和式(8))并表明相較于非線性分布,線性分布所得結果為最優解;而本文認為張拉截斷區域的δ與離散點高度呈線性分布。圖4和表1的對比結果說明了本文方法的正確性和有效性。

表1 完全張拉截斷邊坡穩定系數對比(c/γ=1.0)

4 結果分析

4.1 擬動力法與擬靜力法對比

由于擬靜力法概念清晰和計算簡便,被廣泛應用到實際邊坡工程的抗震穩定性分析中。但其忽略地震隨時間變化的動力特性以及土體剪切模量和剪切波波速視為無窮大會造成高估建(構)筑構的抗震穩定性。因此,首先對比均質邊坡與非均質邊坡在擬動力法和擬靜力法下的安全系數FS。對于均質邊坡,相關參數為:H=5 m,γ=20 kN/m3,c0=ch=15 kPa,φ0=φh=20°,λ=0.5,T=0.3 s,t0=0,f=1.4,VS=150 m/s,VP=280.5 m/s。對于非均質邊坡,相關參數為:H=5 m,γ=20 kN/m3,c0=10 kPa,ch=15 kPa,φ0=15°,φh=20°,λ=0.5,T=0.3 s,t0=0,f=1.4,VS=150 m/s,VP=280.5 m/s。圖5所示為具有不同坡角β邊坡安全系數FS的擬動力法和擬靜力法結果對比。從圖5可以看出: FS隨著水平地震系數kh和β的增大而減小,說明邊坡越陡且地震強度越大會大幅度降低邊坡的穩定性。同時,擬靜力法結果高于擬動力法結果,且隨著kh的增大,兩者結果的差異增大。這是由于擬靜力法將隨地震力簡化為恒定慣性力作用在邊坡上,忽略了地震波周期T,土體放大系數f等對地震加速度的影響,導致地震力做功減小,FS增大。例如對于β=60°均質邊坡,kh=0.1和kh=0.3對應的擬靜力法結果分別高出其擬動力法結果2.95%和7.73%,說明在強震作用下采用擬靜力法會大大高估邊坡的抗震穩定性。除此之外,圖5給出了折減系數ξ=0.0和ξ=1.0情況下邊坡FS的擬動力法結果對比,表明考慮土體完全張拉截斷會降低FS。由于線性M-C屈服準則將土體張拉破壞視為剪切屈服,導致高估了土體的抗拉強度;而考慮土體完全張拉截斷的M-C屈服準則使得土體張拉破壞部分能以斷裂的形式表現出來,用以區別剪切破壞,能體現出邊坡失穩瞬間的狀態。

圖5 不同坡角對邊坡安全系數的影響Fig. 5 Impact of different slope angles on slope safety factor

4.2 擬動力法參數分析

為進一步探討擬動力法對邊坡穩定性的影響,以坡角β=60°非均質邊坡為例,分別改變明地震周期T和土體放大系數f以分析FS的變化。圖6所示為不同地震波時間T對FS的變化影響。由圖可知:kh取值不同時,FS隨T的增大而減小,原因在于T增大,地震波在一個周期內持續時間增加,導致地震力做功增加,從而FS下降。除此之外,FS在T=0.1～0.2 s內變化幅度較大,并隨著kh的增大幅度變化顯著。同時,由于簡化為正弦函數的水平和豎直地震加速度在T≥0.3 s時正弦函數部分約等于1,因此水平和豎直地震加速度變化幅度很小,導致地震力做功功率基本保持不變,因此FS在T=0.3 s后基本保持不變。如圖7所示,FS隨土體放大系數f的增大呈下降趨勢,原因在于f越大,地震波在土體里傳播的加速度幅值越大,導致地震力做功功率增大,對邊坡造成不利影響加大,并且kh的增大會增加FS下降的幅度,例如。當ξ=0.0,kh=0.1時,f由1.0增大到1.8,FS下降了6.37%;當kh=0.3時,FS下降了18.36%。由此可知: 強震作用下較大的f會大幅度削弱邊坡的穩定性,造成嚴重滑坡災害。除此之外,圖6-7再次表明了土體完全張拉截斷效應會削弱邊坡穩定性,且削弱幅度隨著f和kh的增大而增大。在均質邊坡中體現出類似的規律,故不再贅述。

圖6 不同地震波周期對非均質邊坡安全系數的影響Fig. 6 Impact of different seismic wave periods on the safety factor of heterogeneous slopes

圖7 不同土體放大系數對非均質邊坡安全系數的影響Fig. 7 Effect of different soil magnification factors on the safety factor of heterogeneous slope

4.3 不同折減方式對比

圖8所示為消除土體抗拉強度的兩種方法的FS對比:其一為本文方法,即完全張拉截斷;其二為在坡頂處引入一條垂直張拉裂縫。由圖8可知:兩種方法所得FS非常接近,相差不超過4%,在工程允許范圍之內。由于邊坡引入張拉裂縫建模更為簡便,因此在工程實踐中可采取此方法快速獲取裂縫邊坡的安全系數。同時,邊坡穩定性分析另一個核心內容為求解邊坡失穩瞬間的臨界滑動面。圖9給出了kh=0.0和kh=0.2情況下兩種方法所得的臨界滑動面。圖9表明雖然兩種方法所得FS基本相同,但其臨界滑動面差異很大。當kh=0.0時,均質邊坡和非均質邊坡的臨界滑動面基本一致;而兩種方法所得的臨界滑動面差異明顯。完全張拉截斷區域深度小于張拉裂縫的深度約20%,造成的滑動土體體積高出張拉裂縫滑動土體體積約14%。當kh=0.2時,均質邊坡臨界滑動面比非均質邊坡臨界滑動面更淺,說明強震作用會放大土體的非均質性對臨界滑動面的影響。但完全張拉截斷區域深度大于張拉裂縫的深度約60%,造成滑動土體體積高出張拉裂縫滑動土體體積約10%。由此可知:引入張拉裂縫的臨界滑動面比完全張拉截斷的臨界滑動面更淺,邊坡失穩時滑動土體體積較少,邊坡的破壞規模較小。因此,引入張拉裂縫方法會低估邊坡的破壞規模,低估了對邊坡失穩后所造成的災害,導致防災措施不充分,造成大量人員傷亡和財產損失。因此,對邊坡破壞規模的評估應采用完全張拉截斷方法。

圖8 不同折減方式對邊坡安全系數的影響Fig. 8 Impact of different reduction methods on slope safety factor

圖9 不同折減方式對邊坡臨界滑動面的影響Fig. 9 Effect of different reduction methods on the critical sliding surface of slope

5 結論

本文通過引入完全張拉截斷以去除Mohr-Coulomb屈服準則中高估的土體抗拉強度。基于極限分析上限定理,以“點到點”的方式構建離散運動學破壞機構解決土體內摩擦角非均質性對邊坡穩定性影響的問題,運用擬動力法對邊坡的抗震穩定性進行研究。采用二分法和優化算法并結合強度折減法求解邊坡安全系數及臨界滑動面。得出主要結論如下:

1) 由于忽略地震的動力特性和土體放大系數,擬靜力法所得結果高于擬動力法所得結果,且在強震作用會加大兩者計算結果的差異,說明強震作用下采用擬靜力法會明顯高估邊坡穩定性。同時,邊坡安全系數隨土體放大系數f的增大而增大,并在T≥0.3 s基本保持不變。

2) 相較于線性M-C屈服準則,考慮土體完全張拉截斷可以將滑動面張拉破壞區域以斷裂的形式表示,區別塑性剪切屈服區域,并且土體完全張拉截斷會降低邊坡安全系數。

3) 考慮土體完全張拉截斷與引入垂直張拉裂縫兩者所得邊坡安全系數基本一致,但臨界滑動面區別明顯。在靜力狀態下完全張拉截斷區域深度小于張拉裂縫深度;反之,在動力狀態下大于張拉裂縫深度。考慮土體完全張拉截斷的臨界滑動面更深,滑動土體體積更多,表明邊坡失穩時破壞規模更大。引入垂直張拉裂縫建模簡單,可快速獲得裂縫邊坡安全系數;考慮土體完全張拉截斷可用來評估滑坡規模。