王 凱，李同春，程 井

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

0 引 言

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

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

1 深層抗滑穩(wěn)定分析方法

1.1 剛體極限平衡方法

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

1.2 有限單元法

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

1.3 分區(qū)有限元與塊體界面元混合算法

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

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

針對某典型重力壩工程的深層滑動問題，采用分區(qū)有限元與塊體界面元混合算法進行分析?；谠摶旌纤惴ǖ乃悸?，首先需要將求解系統(tǒng)分解為若干塊體和非連續(xù)界面，需預先確定可能的深層滑動面位置，故求解該重力壩深層滑動問題的主要步驟如下：①首先根據(jù)巖基的地質(zhì)構(gòu)造建立有限元模型，采用彈塑性有限元強度折減法確定最可能的滑動面。②修改初始有限元模型，在滑動面位置建立接觸點對。③結(jié)合強度折減法，采用分區(qū)有限元與塊體界面元混合算法對修改后的模型進行計算，以接觸點對全部失效作為失穩(wěn)判據(jù)，求出重力壩抗滑穩(wěn)定安全系數(shù)。

2 典型工程應用

2.1 工程概況

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

2.2 材料參數(shù)及計算工況

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

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

2.3 數(shù)值計算

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

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

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

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

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

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

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

2.3.2 分區(qū)有限元與塊體界面元混合算法

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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

□