程 曦 ，王元戰(zhàn) ，馬殿光，肖 忠，郝林南

（1.天津大學建筑工程學院天津市港口與海洋工程重點實驗室，天津300072；2.交通運輸部天津水運工程科學研究所工程泥沙交通行業(yè)重點實驗室，天津300456）

河道土質岸坡失穩(wěn)破壞是一種幾乎存在于世界上所有江河岸坡的現象。河流上建設樞紐后，航道兩岸塌岸更為嚴重。航道岸坡失穩(wěn)引起沿岸地帶的建筑物破壞，造成財產損失；同時，岸坡失穩(wěn)將侵蝕兩岸土地資源，危害航道通航條件，也極大增加了航道的維護費用。因此，開展航道護岸穩(wěn)定的研究具有十分重要的意義。

影響航道岸坡穩(wěn)定性的因素很多，包括滲流、降雨、波浪力等，本文主要研究穩(wěn)定滲流和波浪力對岸坡的影響。在滲流方面，陳麗剛［1］基于ABAQUS對邊坡穩(wěn)定滲流進行了流固耦合分析。張曉詠［2］應用ABAQUS對壩體的滲流進行了研究，并驗證了強度折減法的可靠性。周群華、章廣成等［3-4］研究了水位變化對邊坡穩(wěn)定性的影響。在船行波方面，項菁等［5-7］對船行波波要素的特點以及現有的一些計算公式進行了總結和分析。盧無疆［8］對高速雙體客船和普通船的船行波進行了現場觀測，得到了各種波浪要素。簡文彬［9］對邊坡對循環(huán)荷載的響應進行了研究。

雖然不少學者已經對邊坡失穩(wěn)的影響因素進行了研究，但都只是針對單一影響因素研究，缺乏對多種因素共同作用下邊坡失穩(wěn)的研究，實際上邊坡失穩(wěn)是多種因素共同作用下的結果，因此多因素共同作用下邊坡的穩(wěn)定性研究也是必不可少的。本文除了分別分析穩(wěn)定滲流和波浪力這2個影響因素之外，還將其組合起來，研究這2種因素共同作用對邊坡穩(wěn)定性的影響。

本文采用大型有限元分析軟件ABAQUS，通過數值模擬對岸坡進行流固耦合分析，并采用強度折減法分析岸坡的穩(wěn)定性。在單獨考慮滲流時，先計算不同靜水位下邊坡的安全系數，然后在數值模型上施加不同的水頭差，計算穩(wěn)定滲流情況下的安全系數，最后與極限平衡方法的計算結果相對比。在單獨考慮波浪力時，參考實測的波浪要素，根據規(guī)范計算出波壓力的分布，并在模型上施加周期荷載，研究不同波高、周期的船行波對邊坡穩(wěn)定的影響。在考慮共同作用時，將之前計算的工況組合起來，研究2種因素共同作用下對岸坡穩(wěn)定性的影響，并與單因素的計算結果進行對比分析。

1 基本原理

1.1 ABAQUS流固耦合分析的基本原理

ABAQUS可對多孔介質的滲流和變形進行耦合分析，可進行非飽和土的滲流計算。ABAQUS中處理孔隙介質中的流體流動的方式是將孔隙體視為多相材料，孔隙中的流體可包含兩部分：一是液體，通常認為壓縮性相對很低；另一個則是氣體，認為是可壓縮的。土體的體積包括兩部分：土顆粒的體積和孔隙的體積，孔隙的體積等于孔隙中液體的體積與氣體體積之和。計算中有限元的網格固定在土骨架上，氣體或液體可流過網格，但需要滿足流體的連續(xù)性方程。土體的力學特性通過采用有效應力定義的本構模型來模擬，液體的滲透采用Forchheimer滲透定律模擬，常用的Darcy定律是它的簡化。

式中：δv為虛速度場；δε=sym（?δv/δx）為虛應變；σ 為真實（Cauchy）應力；t為單位面積上的表面力；為單位體積上的體積力。

滲流連續(xù)方程為

式中：vw為流體的平均流速；n為邊界S的外法線；該方程對ρw0進行了歸一化。Forchheimer滲透定律表達式為

ABAQUS能將滲流場和應力場直接藕合，無需進行滲流場和應力場的反復迭代，只要按時間過程連續(xù)求解就可得到全部結果，即通過將節(jié)點位移和孔隙水壓力作為節(jié)點自由度進行空間離散，將應力平衡方程和滲流連續(xù)方程寫成矩陣形式，并對滲流連續(xù)方程引入時間積分，得到耦合控制方程，然后在每個時間步內求解方程，并同時滿足位移邊界條件和滲流邊界條件。

1.2 有限元強度折減法基本原理

邊坡穩(wěn)定的分析方法主要可分為兩大類：第一類可稱作傳統方法，如極限平衡法、極限分析法、滑移線法等；第二類為基于計算機技術的數值分析方法，如有限元法、離散元法、邊界元法等［10］。對于穩(wěn)定滲流的土坡穩(wěn)定分析，傳統的極限平衡條分法分成2種方法［11］：一種是根據邊界孔隙水壓力和飽和的重度考慮滲流作用的常用方法；另一種是根據滲透力和浮重力考慮滲透作用的分析方法。前者需要繪制流網求得孔隙水壓力，而后者由于計算滲透力很不方便，工作量較大，在工程實踐中較少采用。本文采用的有限元強度折減法，能夠很好地克服極限平衡法的缺點，因而已經被廣泛的應用到了有數值分析中。

抗剪強度折減系數（SSRF：Shear Strength Reduction Factor）的概念是由Zienkiewicz等1975年在《土工彈塑性有限元數值分析》一文中首次提出的，其定義為：在外載荷保持不變的情況下，邊坡內土體所能提供的最大抗剪強度與外載荷在邊坡內所產生的實際剪應力之比。在極限狀況下，外載荷所產生的實際剪應力與抵御外載荷所發(fā)揮的最低抗剪強度（即按照實際強度指標折減后所確定的、實際中得以發(fā)揮的抗剪強度）相等。當假定邊坡內所有土體抗剪強度的發(fā)揮程度相同時，這種抗剪強度折減系數相當于傳統意義上的邊坡整體穩(wěn)定安全系數Fs，又稱為強度儲備安全系數，與極限平衡法中所給出的穩(wěn)定安全系數在概念上是一致的。

折減后的抗剪強度參數表達式為

式中：c、φ為土體所能夠提供的抗剪強度；cm和φm為維持平衡所需要的或土體實際發(fā)揮的抗剪強度；Fs為強度折減系數。

判斷岸坡是否達到臨界破壞的標準主要有以下幾種［12］：（1）以數值計算收斂與否為評價標準，與有限元的算法有關；（2）以特征點的位移-折減系數曲線的拐點為評價標準；（3）以是否形成連續(xù)塑性變形貫通區(qū)為評價標準。

本文采用以是否形成連續(xù)塑性變形貫通區(qū)作為評價標準。

1.3 摩爾庫侖（Mohr-Coulomb）模型及參數的選取

為與傳統極限平衡法對邊坡巖土體進行剛塑性（理想塑性）假設相對應，本文采用的本構模型是Mohr-Coulomb準則。為避免經典Mohr-Coulomb準則在三維主應力空間和π平面上，屈服面存在棱邊和尖角（奇點），導致塑性流動方向不唯一、數值計算繁瑣和收斂緩慢的問題，ABAQUS軟件采用的是Mohr-Coulomb屈服準則的擴展。其采用Mohr-Coulomb屈服函數，包括粘聚力的各向同性硬化和軟化，但該模型的流動勢函數在子午面上的形狀為雙曲線，在π平面上沒有尖角，因此勢函數完全光滑，確保了塑性流動方向的唯一性。

數值模型中需要輸入的主要參數分別是：彈性模量E、泊松比 v、內摩擦角 φ′、粘聚力 c′、剪脹角 Ψ、飽和土的滲透系數k。除此之外還要考慮浸潤面以上的非飽和土中的滲流，計算中還需輸入非飽和土的土水特征曲線。其中為保守考慮剪脹角取值Ψ=0°，其他參數都由通過三軸固結不排水剪（測孔隙水壓力）試驗確定。土樣為飽和粘土，取自廣西壯族自治區(qū)右江上游那吉庫區(qū)段。

在考慮波浪的動荷載作用時，本文主要從粘土在循環(huán)荷載作用下抗剪強度降低這方面考慮，通過動三軸試驗，分別測出動應力與偏應力之比S1=σd1/（σ1-σ3）=0.1、S2=σd2/（σ1-σ3）=0.2時土體的抗剪強度。不同動應力與偏應力之比下土體參數如表1所示。

表1 不同動應力與偏應力之比下土體參數Tab.1 Soil parameter with different strength weakening factors

圖1 岸坡尺寸Fig.1 Slope dimension

2 計算實例

2.1 邊坡穩(wěn)定滲流分析

根據廣西那吉庫區(qū)地形圖可知當地為均質岸坡，坡高8 m，坡比為1：3，航道正常水位為5 m。根據三軸固結不排水試驗可知有效內摩擦角 φ′=14°，有效粘聚力 c′=14 kPa，重度r=18 kN/m3，彈性模量 E=20 MPa，泊松比 v=0.4，滲透系數k=10-8m/s，岸坡尺寸如圖1所示。

有限元模型采用孔壓/位移耦合的CPE8R四邊形八節(jié)點平面應變單元，共含有870個單元、2 757個節(jié)點。上部邊界為自由邊界，邊坡底面采用完全約束條件，左右采用水平約束條件。受力為重力和靜水壓力。網格劃分如圖2所示。

（1）算例一。計算不同靜水位下邊坡的安全系數，水位分別為4 m、5 m、6 m。

首先在不同靜水位條件下使用Geostatic分析步平衡地應力，然后在Soil分析步中選擇穩(wěn)態(tài)分析采用強度折減法計算安全系數，并與Morgenstern-Price法計算結果進行對比（表2）。

圖2 網格劃分Fig.2 Calculation model

表2 不同靜水位下的安全系數Tab.2 Safety factors of slopes under different still water levels

由表2可知，靜水位下安全系數隨著靜水位的上升而增大，有限元強度折減法計算的安全系數與Morgenstern-Price法的誤差最大為3.52%，在合理范圍內。

（2）算例二。計算不同水位穩(wěn)定滲流情況下的岸坡穩(wěn)定性，水頭差分別為1 m、2 m、3 m。

首先使用Geostatic分析步平衡不同水位時的地應力，然后在Soil分析步中選擇穩(wěn)態(tài)分析，并改變邊界的孔壓條件以達到穩(wěn)定滲流，最后采用Soil分析步中的穩(wěn)態(tài)分析進行強度折減以獲得該邊坡穩(wěn)定的安全系數，強度折減系數Fs=1.0時為邊坡初始狀態(tài)。

航道正常水位為5 m，在不同地下水位的工況下岸坡穩(wěn)定滲流的安全系數計算結果如表3所示。

水頭差與安全系數關系曲線如圖3所示。

由表3和圖3可以看出，隨著水頭差的增大邊坡安全系數呈現線性減小趨勢，有限元強度折減法計算的安全系數與Morgenstern-Price法的誤差最大為2.19%，在合理范圍內。

對比有限元強度折減法計算的安全系數與Morgenstern-Price法的計算結果，可以從表2和表3中看出最大的誤差僅為3.52%，這體現出有限元法與傳統極限平衡法的計算結果有較好的一致性，也說明將塑性區(qū)貫通作為失穩(wěn)判據確定邊坡穩(wěn)定安全系數是可行的，也是合理的。

表3 穩(wěn)定滲流下的安全系數Tab.3 Safety factors of slopes under steady seepage

圖3 水頭差與安全系數關系曲線Fig.3 Curve between head differences and safety factors

2.2 波浪對邊坡穩(wěn)定影響分析

算例三：計算不同波浪力作用下岸坡的穩(wěn)定性，此時水位為正常靜水位5 m。

根據盧無疆［8］對高速雙體客船和普通船的船行波的現場觀測結果，2種船型波各選一組波浪要素進行分析。

（1）高速雙體客船：Vma=13.11 m/s，Hmaxa=0.45 m，=2.5 s，作用時間。

（2）普通船：Vmb=6.25 m/s，Hmaxb=0.28 m，=1.8 s，作用時間。

根據海港水文規(guī)范（JTJ 213-98）8.1.4.2，波峰時，靜水面處波壓力 Psa=γHmaxa=4 500 Pa，Psb=γHmaxb=2 800 Pa，并將波浪力按照波峰時的分布簡化為正弦波。波峰時波壓力分布如圖4所示。

首先使用Geostatic分析步平衡正常水位不同土體強度下的地應力，然后采用Soil分析步中的穩(wěn)態(tài)分析進行強度折減，最后采用Soil分析步中的瞬態(tài)分析，并加上不同的波浪力循環(huán)荷載。不斷調節(jié)場變量的值進行計算，直到Soil分析步中的瞬態(tài)分析結束時岸坡剛好形成連續(xù)塑性變形貫通區(qū)，此時對應的場變量值即為邊坡的安全系數。在船行波作用下邊坡的安全系數計算結果如表4所示。

動應力與偏應力之比與安全系數關系曲線如圖5所示。

圖4 波壓力分布示意圖Fig.4 Sketch of wave pressure distribution

表4 船行波作用下的安全系數Tab.4 Safety factors of slopes under wave pressure

從表4和圖5可以看出，高速雙體客船對岸坡安全系數的影響要比普通船大0.015左右。波浪力在土體強度未折減時對邊坡的穩(wěn)定性影響不大，但是粘土在循環(huán)動荷載作用下發(fā)生強度弱化對邊坡的穩(wěn)定性有較大影響，隨著弱化程度的加劇，邊坡的安全系數基本呈線性減小的趨勢。

2.3 穩(wěn)定滲流和波浪共同作用下邊坡穩(wěn)定影響分析

算例四：選取穩(wěn)定滲流水頭差為3 m的工況與船行波的各種工況組合進行分析。

首先使用Geostatic分析步平衡正常水位不同土體強度下的地應力，然后采用Soil分析步中的穩(wěn)態(tài)分析進行強度折減，再采用Soil分析步中選擇穩(wěn)態(tài)分析，并改變邊界的孔壓條件以達到穩(wěn)定滲流，最后采用Soil分析步中的瞬態(tài)分析，并加上不同的波浪力循環(huán)荷載。不斷調節(jié)場變量的值進行計算，直到Soil分析步中的瞬態(tài)分析結束時岸坡剛好形成連續(xù)塑性變形貫通區(qū)，此時對應的場變量值即為邊坡的安全系數。計算結果如表5所示。

動應力與偏應力之比與安全系數關系曲線如圖6所示。

從表5和圖6可以看出，高速雙體客船和普通船船行波作用下的岸坡安全系數之間的差別較小，滲流和波浪力的共同作用對岸坡的穩(wěn)定性已有較大影響，在土體強度未弱化時，邊坡的安全系數的降幅達到0.12以上。粘土在循環(huán)動荷載作用下發(fā)生強度弱化時，隨著弱化程度的加劇，邊坡的安全系數也是基本呈線性減小的趨勢。

有無穩(wěn)定滲流下邊坡的安全系數對比結果如圖7所示。

由圖7可以看出，土體強度的弱化對岸坡穩(wěn)定會產生較大影響，最大差值達到0.223。無論是在有無滲流的作用下，粘土在循環(huán)動荷載作用下發(fā)生強度弱化時，隨著弱化程度的加劇，邊坡的安全系數減小的趨勢是相同的，基本為線性分布，且相互平行。由此可見考慮土體強度弱化時，穩(wěn)定滲流和波浪力這2個影響因素是相互獨立的。

圖5 動應力與偏應力之比與安全系數關系曲線Fig.5 Curve between reduction factors and safety factors

表5 共同作用下的安全系數Tab.5 Safety factors of slopes under interaction

圖6 動應力與偏應力之比與安全系數關系曲線Fig.6 Curve between strength weakening factors and safety factors

圖7 動應力與偏應力之比與安全系數關系曲線Fig.7 Curve between strength weakening factors and safety factors

算例五：不考慮土體強度弱化，計算不同水頭差穩(wěn)定滲流和船行波的各種工況組合進行分析。

首先使用Geostatic分析步平衡正常水位下的地應力，然后采用Soil分析步中的穩(wěn)態(tài)分析進行強度折減，再采用Soil分析步中選擇穩(wěn)態(tài)分析，并改變邊界的孔壓條件以達到不同水頭差下的穩(wěn)定滲流，最后采用Soil分析步中的瞬態(tài)分析，并加上不同的波浪力循環(huán)荷載。不斷調節(jié)場變量的值進行計算，直到Soil分析步中的瞬態(tài)分析結束時岸坡剛好形成連續(xù)塑性變形貫通區(qū)，此時對應的場變量值即為邊坡的安全系數。計算結果如表6所示。

共同作用時不同水頭差下邊坡的安全系數對比結果如圖8所示。

由圖8可以看出，不考慮土體強度弱化時，船行波對岸坡穩(wěn)定性的影響最大只有0.034，在不同水頭差的穩(wěn)定滲流船行波對安全系數的影響基本為一定值，由此可見在不考慮土體強度弱化時，穩(wěn)定滲流和波浪力這2個影響因素也可以看作是相互獨立的。

表6 共同作用時不同水頭差下的安全系數Tab.6 Safety factors of slopes under head differences

圖8 水頭差與安全系數關系曲線Fig.8 Curve between head differences and safety factors

3 結論

通過以上算例的計算與分析可以得到以下結論：

（1）本文采用的ABAQUS有限元強度折減法計算的安全系數與傳統極限平衡Morgenstern-Price法的計算結果無論是在靜水位還是穩(wěn)定滲流情況下都只產生了較小誤差，驗證了本文采用方法的可靠性，同時也說明以是否形成連續(xù)塑性變形貫通區(qū)作為邊坡破壞的評價標準是可行的。

（2）岸坡的安全系數隨著水位的增長而增大。在岸坡穩(wěn)定滲流時，水差越大安全系數越小。高速雙體客船船行波比普通船船行波對岸坡的影響要大。粘土在循環(huán)動荷載作用下發(fā)生強度弱化會對邊坡的穩(wěn)定性產生較大影響。

（3）當穩(wěn)定滲流和波浪力共同作用時，邊坡安全系數最小，其變化趨勢與無滲流的情況相同。當穩(wěn)定滲流和波浪力共同作用時，2個因素之間不互相影響，可以看作是相互獨立的。

本文建立的數值模型與采用的土體參數均來自廣西壯族自治區(qū)右江上游那吉庫區(qū)，這對當地岸坡的穩(wěn)定性研究具有重要的意義，同時也為研究邊坡在多因素共同作用下的穩(wěn)定性提供參考。

［1］陳麗剛.基于ABAQUS滲流與應力耦合作用的邊坡穩(wěn)定性分析［D］.鄭州：鄭州大學，2010.

［2］張曉詠，戴自航.應用 ABAQU 程序進行滲流作用下邊坡穩(wěn)定性分析［J］.巖石力學與工程學報，2010，29（1）：2 927-2 934.ZHANG X Y，DAI Z H.Analysis of Slope Stability under Seepage by Using ABAQUS Program［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（1）：2 927-2 934.

［3］周群華，閆澍旺，鄧衛(wèi)東，等.水庫水位變化對某岸區(qū)公路邊坡穩(wěn)定性影響分析［J］.公路交通科技，2011，28（6）：32-39.ZHOU Q H，YAN S W，DENG W D，et al.Analysis of Influence of Reservoir Water Level Fluctuation on Stability of a Bank Slope［J］.Journal of Highway and Transportation Research and Development，2011，28（6）：32-39.

［4］章廣成.水位變化對滑坡穩(wěn)定性的影響研究［D］.武漢：中國地質大學，2005.

［5］項菁，石根娣.天然航道船行波波高計算方法［J］.河海大學學報，1994，22（2）：47-48.XIANG J，SHI G D.Calculation of Ship Wave Heights in Navigable Channels［J］.Journal of Hohai University，1994，22（2）：47-48.

［6］張璠，張緒進，尹崇清.船行波與運河岸坡的研究綜述［J］.中國水運：學術版，2006，6（5）：45-50.ZHANG F，ZHANG X J，YIN C Q.Research Summary of Ship Wave and Canal Slope［J］.China Water Transport：Academic Version，2006，6（5）：45-50.

［7］潘寶雄，蔣宗燕.船行波問題的研究綜述［J］.河港工程，2001（1）：14-19.PAN B X，JIANG Z Y.Research Summary of Ship Wave［J］.River Port Engineering，2001（1）：14-19.

［8］盧無疆.高速雙體客船船行波特性現場觀測［J］.水利水運科學研究，1995（2）：191-196.LU W J.On Prototype Observation of Ship Wave Characteristics of High-speed Double Volume Ship［J］.Hydro-Science and Engineering，1995（2）：191-196.

［9］簡文彬，胡忠志，樊秀峰，等.邊坡對循環(huán)荷載的響應研究［J］.巖石力學與工程學報，2008，27（12）：2 563-2 567.JIAN W B，HU Z Z，FAN X F，et al.Research on Responses of Slope Under Cyclic Load［J］.Chinese Journal of Rock Mechanics and Engineering，2008，27（12）：2 563-2 567.

［10］Abramson L W，Lee T S，Sharma S，et al.Slope stability and stabilization methods［M］.New York：Wiley，2001.

［11］楊進良.土力學［M］.北京：中國水利水電出版社，2000.

［12］費康，張建偉.ABAQUS在巖土工程中的應用［M］.北京：中國水利水電出版社，2009.