基于分區有限元與塊體界面元混合算法的重力壩深層抗滑穩定分析

王 凱，李同春，程 井

(1. 河海大學水利水電學院，南京 210098；2. 河海大學農業工程學院，南京 210098； 3. 河海大學水資源高效利用與工程安全國家工程研究中心，南京 210098)

0 引 言

近幾十年來，水電事業得到了蓬勃發展，重力壩作為重要壩型得到了廣泛應用。隨著資源不斷開發利用，不可避免地需要在地質條件復雜的巖基上修建重力壩。由于巖體中普遍存在節理、夾層、斷層等軟弱結構面，特別是緩傾角的結構面在重力和水推力的作用下，有可能形成滑移通道，導致大壩失穩破壞，因此重力壩深層抗滑研究具有實際的工程意義。

巖體中存在的軟弱結構面將巖體分割成不連續體，節理、夾層、斷層等軟弱結構面構成了不連續結構面。壩體失穩往往伴隨著巖體中這些不連續結構面地張開、閉合和滑移等現象，屬于典型的局部位移非線性的接觸問題，因此深層抗滑穩定分析的關鍵在于如何有效地模擬巖石裂隙的力學行為[1,2]。近年來，主要采用剛體極限平衡法和基于連續介質理論的有限單元法[3-5]進行重力壩深層抗滑穩定分析。剛體極限平衡法由于其概念清晰，計算簡便，被大量的工程采用，積累了豐富的工程經驗，同時具備一套成熟的設計準則，是一種十分常見的穩定分析方法。但剛體極限平衡法為了求解方便，在力學分析上作了較大的人為假定，與實際情況偏差較大。隨著計算機技術的發展，以有限單元法為代表的數值計算方法得到廣泛運用。然而有限元法屬于連續介質方法，必須滿足變形相容性條件，無法真實模擬裂縫地張開滑移。為了解決此類局部位移非連續問題，李同春等[6]提出了分區有限元與塊體界面元混合算法，可以較好地模擬不連續結構面的閉合、滑移和張開，同時僅將非線性迭代收縮到接觸界面上進行，大大提高了計算效率。本文在此方法基礎上，針對某重力壩深層抗滑穩定問題，結合強度折減法，利用接觸點對的狀態來模擬深層滑動面的力學行為，以失效點對首次貫通作為失穩判據，求解得到深層抗滑穩定安全系數。該方法將重力壩深層抗滑問題考慮為接觸問題，為今后的深層抗滑穩定研究提供了一種新的求解思路和方法。

1 深層抗滑穩定分析方法

1.1 剛體極限平衡方法

剛體極限平衡法的求解思路是根據巖基中存在的軟弱結構面假設可能的深層滑動面，一般的假設滑面形式包括單滑面、雙滑面以及多滑面。然后將滑裂體看作不變形的剛體，不考慮力矩平衡，只通過力的平衡求出滑面上的滑動力和抗滑力，將抗滑力與滑動力比值作為抗滑穩定安全系數。最常見的滑面形式為雙滑面，根據設計規范采用等安全系數法，用試算法或迭代法求得整體抗滑穩定安全系數K及抗力Q值。該方法應用十分廣泛，但也存在很多缺點，比如計算時需首先假定滑動面，通過試算確定最危險滑面，但對于地質條件復雜的工程，試算工作量巨大；并且將巖體考慮為剛體，不能考慮其受力變形所帶來的影響，極限狀態與允許的工作狀態有較大的出入[3]；不能確定滑面上的應力分布，不能探索破壞的機理及其變化發展過程；等等。但它是經過大量的工程實踐驗證過的，具有很高的可信度，被規范[7]所采用。

1.2 有限單元法

有限單元法可以考慮巖體的非線性彈塑性本構關系，以及變形對應力的影響，結合強度折減法，可以模擬巖體漸進破壞到整體剪切破壞過程。通過有限元計算得到壩基的塑性區圖，基于最大最小值理論及有限元的最小勢能原理，可以近似認為塑性區圖上塑性應變值最大的點的連線(對平面問題)即為臨界滑動面[8]，即可得到重力壩的失穩模式。為了獲得抗滑穩定安全系數，一般需要結合失穩判據確定。常用的失穩判據包括特征點位移突變和塑性區貫通[5]，但這兩種失穩判據均存在較大的人為性，比如位移突變拐點常常不明確，得到的抗滑穩定安全系數為一個區間，以及塑性值達到多大才算貫通一直沒有統一標準，判據還存在較大爭議[9]。同時，有限元法基于連續介質模型，需滿足位移協調條件，無法模擬巖石裂隙地張開滑移，不能真實反映滑動面上的力學特性，而且彈塑性有限元法還存在計算不容易收斂、接近臨界失穩狀態無法求解等問題，因此有必要尋找一種既能反映巖體的連續-非連續位移特征，又兼具較高計算效率的分析方法。

1.3 分區有限元與塊體界面元混合算法

分區有限元與塊體界面元混合算法將求解系統分解為若干塊體與非連續界面，每個塊體作為系統有限元的一個分區，將不連續結構面視為接觸界面，在接觸界面上生成接觸點對，通過接觸點對的閉合、滑移、張開模擬界面的破壞情況。以塊體內的位移作為分區有限元求解的基本變量，以塊體形心的剛體位移作為剛體運動變量，通過界面上力與變形之間的關系和剛體平衡方程，形成以界面上的接觸內力和塊體形心剛體位移為混合變量與界面結點總位移進行迭代求解的混合方程[10,11]。具體的方程推導見文獻[6]。

分區有限元與塊體界面元混合算法計算結果只能得到接觸點對的狀態、結點的應力應變和塊體的剛體位移。因此將混合算法與強度折減法相結合，不斷折減滑動面上的抗剪斷強度參數，以接觸點對首次全部失效作為失穩判據，將此時的折減系數作為抗滑穩定安全系數。

針對某典型重力壩工程的深層滑動問題，采用分區有限元與塊體界面元混合算法進行分析。基于該混合算法的思路，首先需要將求解系統分解為若干塊體和非連續界面，需預先確定可能的深層滑動面位置，故求解該重力壩深層滑動問題的主要步驟如下：①首先根據巖基的地質構造建立有限元模型，采用彈塑性有限元強度折減法確定最可能的滑動面。②修改初始有限元模型，在滑動面位置建立接觸點對。③結合強度折減法，采用分區有限元與塊體界面元混合算法對修改后的模型進行計算，以接觸點對全部失效作為失穩判據，求出重力壩抗滑穩定安全系數。

2 典型工程應用

2.1 工程概況

圖1 擋水壩段剖面圖(單位：m)Fig.1 Retaining dam section profile

2.2 材料參數及計算工況

混凝土壩體分三期施工，一期、二期和三期分別采用C30、C25、C20混凝土進行澆筑。巖基主要為大理巖夾云母石英片巖。在壩踵附近發現有一條傾向下游的夾層，編號f511，寬度約為15cm。壩體、地基及軟弱夾層的物理力學參數如表1所示。

選取正常蓄水位工況進行分析計算。重力壩整體結構所受的荷載主要包括壩體自重、泥沙壓力、上下游水壓力以及地基受到的滲透壓力。初始地應力僅考慮為地基自重所引起的地應力場。上下游正常蓄水位高程分別為2 702.0和2 544.32 m，水重度取9.81 kN/m3。淤沙高程2 595.35 m，浮容重7.8 kN/m3，內摩擦角12°。對于滲透壓力，考慮壩體為不透水材料，僅考慮地基受滲流力作用，先根據各分區滲透系數等計算參數，進行穩定滲流場的計算，得到滲流體積力，然后將計算得到的滲透壓力作為結點荷載施加到地基中[12]。

2.3 數值計算

2.3.1 彈塑性有限元強度折減法

建立重力壩初始有限元模型，如圖 2所示。壩基范圍分別向上下游、底部延伸2倍壩高。整體采用4結點單元進行離散，共計2488個節點，2368個單元。壩體材料采用線彈性本構模型，軟弱夾層采用MC屈服準則，地基巖石采用DP屈服準則。壩基兩側采用法向約束，底部采用固定約束。

(1)可能滑動面搜索。首先進行壩基滲流場計算，考慮防滲帷幕的作用，上下游河谷分別施加相應水頭，其他作為不透水邊界，得到滲透壓力分布如圖 3所示。將滲透壓力作為節點荷載施加到地基中，并施加其他荷載，不斷折減地基材料的抗剪斷強度參數，塑性區不斷擴展并最終貫通。根據之前的分析，將塑性區圖上塑性應變值最大的點連接起來即為最危險滑動面，求得的可能滑動面標記如圖 4所示。可以看到軟弱夾層自然構成深層滑移的主滑面，滑移剪出面角度約20°左右，從壩趾處附近剪出。滑面為典型的雙滑面形式。

圖3 巖基滲透壓力云圖(單位：kPa)Fig.3 Pore-pressure nephogram of rock bed

圖4 巖基塑性區通道示意圖(折減系數Fs=2.4)Fig.4 Plastic zone of rock bed (reduction factor Fs=2.4)

(2)基于彈塑性有限元強度折減法的抗滑穩定安全系數。選取壩頂和壩趾為特征點，圖 5為特征點水平位移隨折減系數的過程線。從曲線圖可以看出，兩特征點位移拐點皆不明顯，無法直接準確獲得抗滑穩定安全系數，塑性值的選取也沒有統一的標準，該方法的失穩判據在確定抗滑穩定安全系數數值時具有不可避免的人為因素。對于本工程，結合位移拐點和塑性區貫通判據，綜合判斷抗滑穩定安全系數約在2.3～2.5之間。

圖5 特征點位移過程線Fig.5 Curve of x displacement-reduction factor for feature points

2.3.2 分區有限元與塊體界面元混合算法

基于分區有限元與塊體界面元混合算法，在搜索得到的可能滑動面上建立接觸點對，修改后的有限元模型見圖6。模型一共包含2 931個節點，2 804個單元，滑動面上一共生成53對接觸點對。壩體和部分地基作為滑動體，滑動面以下的地基作為基礎。計算時僅考慮滑動面上的強度特性，滑動塊體和基礎塊體按彈性體考慮。采用分區有限元與塊體界面元混合算法對系統進行計算，得到接觸點對的狀態。接觸點對所處狀態定義有3種：0為張開，1為閉合，2為滑移，計算時認為接觸點對處于張開或滑移狀態判定為失效。折減系數從1.00開始逐漸遞增，初始每次遞增0.10，在接觸點對幾乎全部失效，即將貫通時，變為每次遞增0.01，統計不同折減系數下的接觸點對的狀態，結果如表2和圖7所示。可以看出隨著折減系數不斷增大，接觸點對的失效比例逐漸上升。從圖8失效點對的位置可以發現，軟弱夾層抗剪斷強度參數低，其上的接觸點對首先失效。折減系數不斷增大時，剪出面上的接觸點對開始失效，并最終全部失效貫通，此時判定大壩失穩。

圖6 重力壩修改后的有限元模型(含滑動面)Fig.6 Modified FEM meshes of gravity dam (sliding surface included)

序號折減系數接觸點對失效組數/接觸點對總組數11.0024/5321.3025/5331.4030/5341.5037/5351.6043/5361.8046/5371.9049/5382.0150/5392.1051/53102.1353/53

圖7 接觸點對失效比例隨折減系數變化過程線Fig.7 Curve of failure ratio of nodal pairs-reduction factor

圖8 不同折減系數下接觸點對狀態圖(Δ為失效)Fig.8 States of nodal pairs of different reduction factor (Δ represents failure)

根據統計結果，當折減系數為2.13時，接觸點對首次全部失效，認為重力壩失穩，求得抗滑穩定安全系數為2.13。結果與彈塑性有限元強度折減法求得結果相比，數值上十分接近，證明了此方法的正確性。但前者比后者結果要小，原因在于分區有限元與塊體界面元混合算法很好地考慮了滑動面上的力學行為，模擬了裂縫的滑移與張開。而有限單元法基于連續介質模型，必須滿足位移協調性，滑動面無法張開或滑移，無法及時進行應力重分配，導致法向壓應力相對增加和切向合力減小[1]，滑動面的安全裕度較大。

3 結 論

本文將重力壩深層抗滑問題考慮為局部位移非連續的接觸問題，采用分區有限元與塊體界面元混合算法對某重力壩進行深層抗滑穩定分析，并與傳統的彈塑性有限元強度折減法進行對比，得到以下結論。

(1)分區有限元與塊體界面元混合算法可以較好地模擬不連續結構面的力學行為以及破壞過程。

(2)前者僅將非線性迭代收縮到接觸界面上進行，具有很好的收斂性和計算效率。

(3)判據明確，避免了彈塑性有限元強度折減法在抗滑穩定安全系數取值上的人為性。

(4)由于該方法需預先知道可能的滑動面，以在接觸區域設置接觸點對，故適用于已知滑動面的深層滑動問題。對于未知滑動面問題，可采用彈塑性有限元強度折減法對滑動面進行搜索，或者預設可能的滑動面進行求解。

□