基于分區(qū)有限元與塊體界面元混合算法的重力壩深層抗滑穩(wěn)定分析

王 凱，李同春，程 井

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

0 引 言

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

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

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

1.1 剛體極限平衡方法

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

1.2 有限單元法

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

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

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

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

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

2 典型工程應(yīng)用

2.1 工程概況

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

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

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

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

2.3 數(shù)值計(jì)算

2.3.1 彈塑性有限元強(qiáng)度折減法

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

(1)可能滑動(dòng)面搜索。首先進(jìn)行壩基滲流場(chǎng)計(jì)算，考慮防滲帷幕的作用，上下游河谷分別施加相應(yīng)水頭，其他作為不透水邊界，得到滲透壓力分布如圖 3所示。將滲透壓力作為節(jié)點(diǎn)荷載施加到地基中，并施加其他荷載，不斷折減地基材料的抗剪斷強(qiáng)度參數(shù)，塑性區(qū)不斷擴(kuò)展并最終貫通。根據(jù)之前的分析，將塑性區(qū)圖上塑性應(yīng)變值最大的點(diǎn)連接起來(lái)即為最危險(xiǎn)滑動(dòng)面，求得的可能滑動(dòng)面標(biāo)記如圖 4所示。可以看到軟弱夾層自然構(gòu)成深層滑移的主滑面，滑移剪出面角度約20°左右，從壩趾處附近剪出?；鏋榈湫偷碾p滑面形式。

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

圖5 特征點(diǎn)位移過(guò)程線Fig.5 Curve of x displacement-reduction factor for feature points

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

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

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

序號(hào)折減系數(shù)接觸點(diǎn)對(duì)失效組數(shù)/接觸點(diǎn)對(duì)總組數(shù)11.0024/5321.3025/5331.4030/5341.5037/5351.6043/5361.8046/5371.9049/5382.0150/5392.1051/53102.1353/53

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

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

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

3 結(jié) 論

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

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

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

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

(4)由于該方法需預(yù)先知道可能的滑動(dòng)面，以在接觸區(qū)域設(shè)置接觸點(diǎn)對(duì)，故適用于已知滑動(dòng)面的深層滑動(dòng)問(wèn)題。對(duì)于未知滑動(dòng)面問(wèn)題，可采用彈塑性有限元強(qiáng)度折減法對(duì)滑動(dòng)面進(jìn)行搜索，或者預(yù)設(shè)可能的滑動(dòng)面進(jìn)行求解。

□