黃 健，程 鵬，鄭小賀，馬懷發(fā)

（1．山東泰山抽水蓄能電站有限責(zé)任公司，山東 泰安 271000；2．中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；3．中國水利水電科學(xué)研究院 工程抗震研究中心，北京 100048）

1 研究背景

大型水利水電工程常建在地質(zhì)構(gòu)造復(fù)雜的高山峽谷地區(qū)，而在庫區(qū)的高邊坡穩(wěn)定與否會直接影響壩體和電站的運(yùn)行安全。一般情況下，這類邊坡的穩(wěn)定性受許多因素的影響。這些因素包括：地形地貌、邊坡巖體節(jié)里面的空間展布、巖體巖性、水文地質(zhì)條件等。地震荷載作用可導(dǎo)致邊坡整體下滑力加大，降低其穩(wěn)定性。在高烈度區(qū)，地震工況已成為高邊坡穩(wěn)定性評價的控制工況?，F(xiàn)行的水電水利工程邊坡設(shè)計規(guī)范[1]只是對巖石邊坡的靜態(tài)穩(wěn)定性，特別針對平面靜態(tài)問題給出了明確的計算方法，很少涉及地震工況下的巖石高邊坡動態(tài)穩(wěn)定性問題。而現(xiàn)行抗震設(shè)計規(guī)范要求邊坡的抗震分析和安全系數(shù)取值應(yīng)按現(xiàn)行水利水電工程邊坡設(shè)計規(guī)范的相關(guān)規(guī)定執(zhí)行。但規(guī)定對于特別重要的、地質(zhì)條件復(fù)雜的高邊坡工程，應(yīng)進(jìn)行基于動態(tài)分析的專門研究。

進(jìn)行邊坡穩(wěn)定性分析研究可采用諸多數(shù)值方法[2-3]，如剛體極限平衡方法、有限元法、紐馬克滑塊法[4]、離散元法[5]、以及不連續(xù)變形分析方法[6]等。其中有限元法已經(jīng)成為最流行的數(shù)值計算方法之一，也是巖石邊坡穩(wěn)定性分析的主要計算工具，而剛體極限平衡法在工程中的長期應(yīng)用，已形成了基于工程經(jīng)驗的安全系數(shù)概念，并建立了安全判斷準(zhǔn)則。將有限元分析與剛體極限平衡相結(jié)合，再結(jié)合接觸算法，通過選擇適當(dāng)?shù)臐撛诨瑒用?，既可回避求解?fù)雜巖體材料非線性的問題，又可有效地進(jìn)行高邊坡穩(wěn)定性分析評價。

庫區(qū)迎水巖石高邊坡的抗震安全評價需要將邊坡滑塊體、地質(zhì)結(jié)構(gòu)面、庫水作用以及近場地基作為一個結(jié)構(gòu)動力系統(tǒng)，考慮遠(yuǎn)域地基的輻射阻尼，建立三維有限元模型，在地基底部輸入地震波，計算分析邊坡的靜動綜合應(yīng)力反應(yīng)。這一結(jié)構(gòu)動力系統(tǒng)的求解自由度可能高達(dá)幾十萬或上百萬，計算時步也有幾百萬甚至千萬，迫切需要進(jìn)行大規(guī)模并行計算[7]。但目前各種通用商業(yè)軟件在求解一般常見問題時可以顯示出所謂的強(qiáng)大功能，但在考慮接觸非線性時，特別是在進(jìn)行并行計算時，效果并不明顯。因此，本文結(jié)合有限元法、極限平衡法與接觸算法，開發(fā)了具有斷層結(jié)構(gòu)面接觸問題的巖石高邊坡動力計算的專用并行程序，并將其應(yīng)用于山東泰山抽水蓄能電站上庫左岸高邊坡的靜動力穩(wěn)定計算和安全評價。

2 高邊坡接觸非線性靜動力并行計算

2.1 接觸模型及接觸力求解方法具有斷層或節(jié)理面的巖石高邊坡穩(wěn)定問題可轉(zhuǎn)化為接觸問題。接觸問題的特點和難點是接觸邊界和接觸力的未知性。由于接觸界面的區(qū)域大小、位置以及接觸狀態(tài)都是未知的，并且隨時間變化，因此，接觸問題表現(xiàn)出顯著的非線性特征。接觸問題的非線性決定了接觸分析過程中需要經(jīng)常插入接觸界面的搜尋判定。在采用有限元求解接觸問題時，一般將相互接觸物體的接觸界面的搜尋判定轉(zhuǎn)化為接觸體離散單元或節(jié)點的接觸判定[8]。但對于以沿接觸面法向變形為主，相對滑移較小的接觸問題，可以采用點對接觸模型用于模擬單點和另一個確定點之間的接觸，省去接觸點位置的判定。這樣只需判斷點對之間是否接觸、滑移和脫離過程，并根據(jù)接觸條件求解接觸力即可。

對于點對接觸模型，在引入了拉格朗日乘子λ表示接觸力后，含有摩擦接觸的動力學(xué)問題轉(zhuǎn)換為求如下鞍點問題[9]，其能量泛函為：

在式（1）中，g為點對間的位移約束值，U1、U2分別為主從節(jié)點的位移列陣，λi=[λn,λτ]， λl表示局部坐標(biāo)系下的乘子，其下標(biāo)n表示法向乘子，下標(biāo)τ表示切向乘子，總的拉格朗日乘子數(shù)為接觸點對個數(shù)與節(jié)點自由度的乘積。 μ為摩擦系數(shù)，N為接觸面上的初始壓力，c為凝聚力，a為節(jié)點的等效面積，當(dāng)接觸面出現(xiàn)滑動后，令ca=0。

2.2 巖石高邊坡的對等架構(gòu)并行計算程序為了實現(xiàn)大規(guī)模巖石高邊坡的靜動力計算，需要開發(fā)高性能并行計算程序。目前并行計算程序主要采用兩種結(jié)構(gòu)設(shè)計模式，即對等架構(gòu)和主從架構(gòu)[10]。主從架構(gòu)需要所有子區(qū)域都和主進(jìn)程進(jìn)行通信，而對等架構(gòu)程序只需要相鄰子區(qū)域間進(jìn)行通信，通信量相對較少。對大規(guī)模問題，對等架構(gòu)采用并行分區(qū)，可考慮前處理數(shù)據(jù)的分塊輸入，其計算結(jié)果也方便分塊存儲，并進(jìn)行分塊后處理。

在對等架構(gòu)下[11]，參與執(zhí)行的每個處理器地位基本相同，每個處理器都沒有例外地參與執(zhí)行整個計算任務(wù)中某個階段上的任何一個子任務(wù)，沒有特別指定各個子任務(wù)執(zhí)行時所在的具體處理器，即每個子任務(wù)都可能在任何一個處理器上執(zhí)行。對于稀疏線性方程組求解，不僅需要等結(jié)點相關(guān)計算任務(wù)在每個處理器上進(jìn)行均分，而且需要單元相關(guān)計算也在每個處理器上均分。因此，在具體并行計算時，不僅涉及到對結(jié)點與單元兩類計算中每一類內(nèi)部可能存在的通信問題，而且可能由于結(jié)點相關(guān)計算，在并行處理時需要按結(jié)點進(jìn)行均分；而對單元相關(guān)計算，需要按單元進(jìn)行均分，而導(dǎo)致在節(jié)點類并行計算、與單元類并行計算間進(jìn)行數(shù)據(jù)重分布所需要的通信，因此，可能對等架構(gòu)的通信結(jié)構(gòu)更為復(fù)雜，實現(xiàn)起來需要解決的問題更多，但其優(yōu)點是，由于各類計算都在所有參與計算的處理器上進(jìn)行均分，使得計算負(fù)載更為均衡，這有利于在大規(guī)模并行計算機(jī)上進(jìn)行高效可擴(kuò)展并行計算。

按照對等架構(gòu)的并行模式，本文開發(fā)出了巖石高邊坡穩(wěn)定性分析的并行計算程序，其程序模塊如圖1所示。從并行程序框架看，各個進(jìn)程子程序具有相同的副本，其功能相同，只是在計算過程中相鄰子區(qū)域的交界面進(jìn)行信息交換。

圖1 巖石邊坡穩(wěn)定性接觸非線性分析的對等架構(gòu)并行計算程序框架

各個子程序可以實現(xiàn)如下功能：

（1）模型分區(qū)程序。網(wǎng)格分區(qū)：根據(jù)網(wǎng)格分區(qū)程序parmeits的需求，讀取一部分單元信息，并進(jìn)行分區(qū)，得到整體分區(qū)信息，然后整理得到各子區(qū)域節(jié)點對應(yīng)的整體節(jié)點編號，存儲在數(shù)組nodg，并保存文件為nodg_進(jìn)程號，以供后面數(shù)據(jù)分區(qū)使用。

（2）節(jié)點及坐標(biāo)分區(qū)。首先讀取各子區(qū)域的nodg文件，根據(jù)整體網(wǎng)格圖分區(qū)信息和體單元信息，進(jìn)行重疊單元的網(wǎng)格分區(qū)，得到分區(qū)節(jié)點，并將各分區(qū)節(jié)點的坐標(biāo)信息保存文件coor0_進(jìn)程號，另外，根據(jù)重疊單元整理相鄰子區(qū)域間需要通信的節(jié)點關(guān)聯(lián)信息。

（3）物理場數(shù)據(jù)分區(qū)，將分區(qū)數(shù)據(jù)按子區(qū)域進(jìn)行存文件處理。

（4）乘子分區(qū)，將分區(qū)乘子信息按子區(qū)域（進(jìn)程號）進(jìn)行存文件處理。

（5）將數(shù)據(jù)通信修改為讀取各自子區(qū)域的數(shù)據(jù)。由Spart程序接收子區(qū)域的坐標(biāo)值；由sgetpart程序接收子區(qū)域所含節(jié)點信息得到局部單元信息、初邊值信息、節(jié)點ID信息以及子區(qū)域節(jié)點的各自由度對應(yīng)的整體方程號等；由spartlm程序接收子區(qū)域接觸點信息。

（6）由starta給出每個節(jié)點的各個自由度與將來要形成的代數(shù)方程組的變量（即方程號）的對應(yīng)關(guān)系，并得到整體剛度矩陣的存儲格式。

（7）分別進(jìn)行位移和應(yīng)力計算，通過各進(jìn)程間相互通信更新解和保存解，最后由0號進(jìn)程進(jìn)行保存結(jié)果文件操作。

并行計算程序分為兩個模塊，即計算靜力場計算程序模塊和計算地震動響應(yīng)分析的模塊，每個模塊包括位移場的計算和應(yīng)力場的計算。

2.3 并行計算效率

（1）計算模型。以泰山公司抽水蓄能電站上水庫左岸邊坡，在下文中簡稱為“邊坡”，為計算模型，考查以上并行計算程序的計算效率。該邊坡的巖性主要為黑云角閃片巖，間夾淺色花崗片麻巖，發(fā)育fs1小斷層，下盤面呈摩擦鏡面狀[12-13]，上盤巖體較完整。在邊坡開挖過程中沿此結(jié)構(gòu)面卸荷，地表裂開，對此斷層采用了預(yù)應(yīng)力錨索加固處理。在下文的穩(wěn)定計算中，巖體及結(jié)構(gòu)面的物理力學(xué)性能指標(biāo)采用文獻(xiàn)[13]的建議值，地基巖體的動態(tài)變形模量取其靜態(tài)變形模量；地基巖體的動態(tài)抗剪強(qiáng)度參數(shù)的標(biāo)準(zhǔn)值取其靜態(tài)抗剪斷參數(shù)的標(biāo)準(zhǔn)值。

有限元模型以左岸邊坡巖石邊坡滑塊體、地質(zhì)結(jié)構(gòu)面、庫水作用以及近場地基作為整體系統(tǒng)。fs1斷層作為具有摩爾庫倫特性的接觸面，預(yù)應(yīng)力錨索概化六面體直桿，計算模型取庫底平均高程372.23 m，由庫底向下延約102 m至270 m高程，面向左岸邊坡，前后跨度280 m，左右寬度190 m。如圖2所示，模型的有限元網(wǎng)格共有94521個六面體單元，975個三菱柱單元，103 465個節(jié)點，310 395個自由度。下面首先以邊坡的靜力計算考察并行計算效率。

圖2 有限元計算模型

（2）并行計算效率。以階躍動荷載的形式施加靜力荷載，采用動力法計算，以最后的計算穩(wěn)定值作為靜力計算結(jié)果。靜力計算持時6 s，時間增量步為0.000 01 s，計算60萬步得到靜力穩(wěn)定結(jié)果，分別進(jìn)行串行計算和并行計算。串行計算在一臺服務(wù)器上進(jìn)行，其配置為Intel（R）Core（TM）i7-10700CPU@2.9GHz，耗時 45 h。并行機(jī)配置為 Intel（R）Xeon（R）CPU E5 2620 V4（2.1G/睿頻3.0G/8核16線程/20M/8GT/85W）。分別采用如表1所列的不同進(jìn)程數(shù)，完成并行計算。以2個進(jìn)程為基準(zhǔn)，并行計算效率為100%，4、6和8個進(jìn)程的并行計算效率均為91.7%，10、20、30和40個進(jìn)程計算效率依次為83.8%、54.4%、50.2%和47.1%，隨著進(jìn)程數(shù)的增加，計算耗時顯著減少，但并行計算效率逐漸降低。

表1 并行計算效率

注意到并行機(jī)一個節(jié)點為8個CPU，由于在一個節(jié)點內(nèi)集成的8個CPU間進(jìn)行通信，通信效率較高，在一個節(jié)點上的并行計算效率高達(dá)91.7%。超過8個進(jìn)程后，需要跨節(jié)點工作，節(jié)點間的信息交換導(dǎo)致計算效率降低。顯然通信效率與并行機(jī)的架構(gòu)有關(guān)，此外，還與數(shù)值模型的規(guī)模有關(guān)。隨著進(jìn)程的增加，通信界面增多，并行計算效率降低。但是，隨著進(jìn)程數(shù)的增加，總耗時一直在降低。與串行耗時45 h相比，40個進(jìn)程的計算耗時僅為70 min，顯然，并行計算大大提高了計算效率。

本文的后續(xù)工作，將利用上述并行計算程序，并依據(jù)現(xiàn)行規(guī)范，對泰山公司抽水蓄能電站上水庫左岸巖石高邊坡進(jìn)行靜力和動力穩(wěn)定性計算及安全評價。

3 邊坡穩(wěn)定性靜力計算

3.1 荷載組合工況邊坡巖體的自重、靜水壓力等為基本靜荷載，在此基礎(chǔ)上再施加地震動作用，進(jìn)行動力計算。在靜力分析完成之后，在靜力位移和應(yīng)力、應(yīng)變場的基礎(chǔ)上，考慮靜力荷載與地震荷載的共同作用，進(jìn)行動力計算。本文的程序設(shè)計將靜力荷載以階躍荷載施加于邊坡地基系統(tǒng)，在獲取穩(wěn)定的靜態(tài)位移和應(yīng)力場后，再由基巖人工邊界底部輸入地震波，對邊坡動力系統(tǒng)進(jìn)行波動反應(yīng)分析。

（1）計算工況I。認(rèn)為錨索全部失效，無錨索作用力，但考慮庫水壓力、斷層滲透水壓力，進(jìn)行左岸邊坡的靜動力穩(wěn)定分析。此工況作為邊坡靜動態(tài)穩(wěn)定評價的基本工況。

（2）計算工況Ⅱ?？紤]錨索作用力，考慮邊坡巖體自重、斷層滲透水壓力，考察錨索作用對左岸邊坡的靜動力穩(wěn)定性的影響。

對如圖2的有限元計算模型，地基周邊及其底界面采用人工黏彈性邊界[14]，迎水面施加靜水壓力，動水壓力采用附加質(zhì)量，正常蓄水位410 m，在斷層面兩側(cè)施加滲透壓力；錨索兩端施加等效壓力1500 kN，在失效狀態(tài)下取值為零。迎水面水位為410 m，背面水頭為402 m，斷層面兩側(cè)滲透壓力水頭取405 m。

3.2 邊坡靜態(tài)應(yīng)力計算以階躍動荷載的形式施加靜力荷載，在計算中當(dāng)接觸面出現(xiàn)滑動后，將點對的凝聚力賦值為零，以最后穩(wěn)定值作為靜力計算結(jié)果。由于動力計算是以靜力計算結(jié)果為初始條件，因此，模型的約束條件應(yīng)保持一致。

時間增量步為0.000 01 s，在并行服務(wù)器上采用40個進(jìn)程完成并行計算，分區(qū)視圖如圖3所示。在滑塊體上取一點，節(jié)點編號為40 874，坐標(biāo)為（132.40，125.00，403.80），如圖4所示，計算至6 s，60萬步，該點在三個方向的位移過程曲線已趨于靜力穩(wěn)定結(jié)果，計算耗時約70 min。

圖3 并行計算區(qū)域分解（40個分區(qū)）

圖4 滑塊上一代表點的位移時程曲線

兩種工況的靜力應(yīng)力計算結(jié)果在圖5給出。工況Ⅰ的最大拉主應(yīng)力為0.40 MPa，見圖5（a），最大靜態(tài)壓主應(yīng)力為1.42 MPa，見圖5（b）。

工況Ⅱ的最大拉主應(yīng)力為0.20 MPa，見圖5（c），圖5（d）給出工況Ⅱ的最大靜態(tài)壓主應(yīng)力為1.42 MPa。應(yīng)力水平與工況Ⅰ相比，滑塊體的最大拉主應(yīng)力降低了0.2 MPa，最大靜態(tài)壓主應(yīng)力二者基本相同。

圖5 左岸邊坡靜態(tài)應(yīng)力分布

3.3 邊坡靜態(tài)穩(wěn)定性分析評價左岸滑坡及滑塊體如圖2（c）和圖2（d）所示，滑動面與水庫水下岸坡交線高程為380.31 m，高于庫底面（庫底平均高程為372.23 m）約8.08 m?；瑒用孀詅s1斷層，由內(nèi)向外下傾，底部滑動面與水平面向下傾角成10°，滑向迎水方向?；瑒用婷娣e約為3957.21m2，滑動面的揚(yáng)壓力115.14×104kN。滑塊體處于弱風(fēng)化帶，摩擦系數(shù) f′取0.8，抗剪斷凝聚力C′取0.6 MPa。

設(shè)S（*）和N分別為滑動面的切向合力和法向合壓力，A為滑動面面積，R(*)為滑動面處巖體的抗剪斷強(qiáng)度， R(*)=f′N+C′A，安全系數(shù)K=R(*)S。根據(jù)圖5的左岸邊坡靜態(tài)應(yīng)力分布，計算得到2種工況下的靜態(tài)穩(wěn)定性安全系數(shù)如表2所列。

表2 滑動面的靜態(tài)抗滑安全系數(shù)

在無錨索作用下的工況Ⅰ，左岸邊坡靜態(tài)穩(wěn)定安全系數(shù)K為3.19，工況Ⅱ在有錨索加固作用下其靜態(tài)穩(wěn)定安全系數(shù)為3.21，穩(wěn)定性有所提高。兩種工況的靜態(tài)安全系數(shù)K均大于規(guī)范規(guī)定最小安全系數(shù)1.3，滿足規(guī)范規(guī)定的靜態(tài)抗滑穩(wěn)定要求。

4 邊坡動力計算及穩(wěn)定性分析

4.1 輸入場地震波時程泰山抽水蓄能電站工程場地的地震基本烈度應(yīng)為Ⅶ度，地基巖性為混合花崗巖，基巖場地類別應(yīng)劃為I0類場地，相應(yīng)的特征周期調(diào)整為0.25s?？紤]到上水庫左岸邊坡巖體結(jié)構(gòu)復(fù)雜性，在基本烈度基礎(chǔ)上提高1度，將設(shè)計烈度定為Ⅷ度進(jìn)行安全評價，設(shè)計水平向地震動峰值加速度代表值為0.2g，特征周期為0.25s；標(biāo)準(zhǔn)反應(yīng)譜最大值 βmax取2.5，阻尼比取5%，生成3組加速度時程[15-16]，分別如圖6、圖7和圖8所示，并由此積分得到速度和位移時程。將該人工（設(shè)計地震）地震波折半，沿近場壩基底部人工邊界的橫向和縱向輸入，再取其2/3的振幅作為豎向輸入地震荷載。

圖6 第一組場地（設(shè)計地震）地震波

圖7 第二組場地（設(shè)計地震）地震波

圖8 第三組場地（設(shè)計地震）地震波

4.2 邊坡的動態(tài)應(yīng)力計算將工況Ⅰ的靜力計算結(jié)果作為動力計算的初始條件，分別輸入圖6、圖7和圖8三組場地地震波，用附加質(zhì)量考慮迎水面邊坡的動水壓力，進(jìn)行左岸邊坡在設(shè)計地震作用下動力響應(yīng)計算。其地震動持時為30 s，時間增量步為0.000 01 s，計算300萬時間步，用40個進(jìn)程完成動力并行計算，計算耗時約9小時40分。計算得到如圖9、圖10和圖11所給出的三組地震荷載作用下的應(yīng)力分布。由此可以看出，三組人工地震波作用的最大主拉應(yīng)力和最大主壓應(yīng)力均存在局部應(yīng)力集中區(qū)。

圖9 在第一組場地地震波作用下的計算結(jié)果

圖10 在第二組場地地震波作用下的計算結(jié)果

圖11 在第三組場地地震波作用下的計算結(jié)果

考慮到邊坡基巖體為斑紋狀混合巖或混合花崗巖，在強(qiáng)風(fēng)化區(qū)的抗壓強(qiáng)度也高達(dá)90 MPa，抗拉強(qiáng)度保守估值（抗壓強(qiáng)度的8%）也在7 MPa以上。為了清楚地顯示滑塊主體的應(yīng)力分布梯度，應(yīng)力場的最大值控制在7.00 MPa以下，可以清楚地看出高應(yīng)力區(qū)僅占滑塊體邊角的很小區(qū)域。

4.3 邊坡的動態(tài)穩(wěn)定分析評價現(xiàn)行水工建筑物抗震設(shè)計標(biāo)準(zhǔn)規(guī)定，邊坡的抗震分析和安全系數(shù)取值應(yīng)按《水電水利工程邊坡設(shè)計規(guī)范》DL/T 5353的相關(guān)規(guī)定執(zhí)行。邊坡穩(wěn)定分析評價采用傳統(tǒng)的確定性方法，即將作用荷載和結(jié)構(gòu)抗力看作為確定性值，計算動態(tài)穩(wěn)性安全系數(shù)K，與靜態(tài)穩(wěn)定性分析不同的是，這里計算得到的穩(wěn)定安全為滑動面的安全系數(shù)時程。

在進(jìn)行巖石高邊坡動力時程分析時，通過有限元計算得到各時刻滑動面上的應(yīng)力分布，將其積分得到作用在塊體上的滑動力與抗滑力，按剛體極限平衡法給出滑動面的安全系數(shù)時程，再基于工程經(jīng)驗的安全系數(shù)概念以最小安全系數(shù)評價巖石邊坡的穩(wěn)定性。

由圖9、圖10和圖11三組場地地震波作用下的動應(yīng)力計算結(jié)果，可以相應(yīng)地計算得到如圖12、圖13和圖14場地地震波作用下滑動面的穩(wěn)定安全系數(shù)時程曲線。

圖13 在第二組場地地震波作用下滑動面的安全系數(shù)時程

圖14 在第三組場地地震波作用下滑動面的安全系數(shù)時程

在三組地震波作用下的滑動面的動態(tài)抗滑安全系數(shù)最小值分別為2.65、2.70和2.75，動態(tài)抗滑安全系數(shù)均遠(yuǎn)大于規(guī)范規(guī)定的安全系數(shù)1.05，邊坡動態(tài)抗滑穩(wěn)定滿足規(guī)范的要求。

4.4 邊坡錨索對邊坡穩(wěn)定性作用評價工況Ⅱ的動力計算主要考察錨索作用力對邊坡動力響應(yīng)的影響。以工況Ⅱ的靜力計算結(jié)果作為動力計算的初始條件，輸入第一組如圖6的場地地震波，用附加質(zhì)量考慮迎水面邊坡的動水壓力，進(jìn)行動力響應(yīng)計算，以考察錨索作用的效果。如圖15所示，類似于工況Ⅰ的動應(yīng)力分布，在位于滑塊體棱角的小范圍內(nèi)出現(xiàn)高于7.00 MPa的應(yīng)力集中現(xiàn)象，其他大部分區(qū)域的應(yīng)力都在7.00 MPa以下。

圖15 工況Ⅱ的動力計算結(jié)果

由圖15的動應(yīng)力計算結(jié)果，計算得到如圖16在場地地震波作用下滑動面的穩(wěn)定安全系數(shù)時程曲線。由此得到滑動面的動態(tài)抗滑安全系數(shù)最小值為2.75。

圖16 有錨索作用時的滑動面抗滑穩(wěn)定安全系數(shù)的時程

在表3中分別列出了在有無錨索作用時的靜動態(tài)穩(wěn)定安全系數(shù)。有錨索和無有錨索作用時的靜態(tài)穩(wěn)定的安全系數(shù)分別為3.21和3.19，有錨索時提高了0.02；動態(tài)穩(wěn)定的安全系數(shù)分別為2.75和2.65，有錨索時提高了0.10。

表3 錨索作用對滑動面抗滑穩(wěn)定性的影響

5 結(jié)論

本文將巖石邊坡的斷層作為接觸面，采用點對接觸模型，開發(fā)出了巖石高邊坡靜動態(tài)穩(wěn)定性的并行計算程序。程序包括靜力場計算模塊和地震動響應(yīng)分析模塊，每個模塊都具有位移場和應(yīng)力場的計算功能。由于采用了對等架構(gòu)的并行模式，在并行計算時各個并行分區(qū)僅需相鄰子區(qū)域的交界面進(jìn)行信息交換。計算實例顯示，該程序具有很高的并行計算效率，最高效率可達(dá)91.7%，同時發(fā)現(xiàn)并行計算效率不僅與并行機(jī)的硬件架構(gòu)有關(guān)，還與數(shù)值模型的規(guī)模有關(guān)。隨著進(jìn)程的增加，通信界面增多，會降低并行計算效率，因此，需要合理選擇與計算規(guī)模相匹配的并行計算進(jìn)程數(shù)。

作為工程應(yīng)用實例，本文詳細(xì)地介紹了山東泰山抽水蓄能電站上庫左岸邊坡穩(wěn)定性計算和安全評價的全過程。該邊坡有限元模型由邊坡潛在滑動體、斷層面、近場地基及人工黏彈性邊界構(gòu)成。靜力計算以階躍動荷載的形式施加靜力荷載，以穩(wěn)定值為靜力計算結(jié)果。動力計算在底部人工邊界橫向和縱向輸入地震波。由有限元計算得到滑塊界面在每一時間步長內(nèi)的應(yīng)力分布，再將其應(yīng)力積分得到作用在塊體上的滑動力與抗滑力，按剛體極限平衡法計算滑動面的安全系數(shù)時程，以安全系數(shù)最小值，基于現(xiàn)行規(guī)范進(jìn)行邊坡穩(wěn)定性評價。

對泰山抽水蓄能電站上庫左岸邊坡穩(wěn)定安全評價得到如下結(jié)論：（1）邊坡深層靜態(tài)抗滑安全系數(shù)為3.19，動態(tài)抗滑安全系數(shù)最小值為2.75，邊坡滿足規(guī)范規(guī)定的靜態(tài)和動態(tài)抗滑穩(wěn)定要求；（2）錨索作用使靜態(tài)和動態(tài)穩(wěn)定安全系數(shù)計算值分別提高了0.02和0.10，局部預(yù)應(yīng)力錨索的加固措施，對邊坡深層抗滑穩(wěn)定性的貢獻(xiàn)不太顯著。

不同工程的計算結(jié)果可能不同，但本文的穩(wěn)定性計算過程和評價方法可為類似工程提供參考。另外，需要指出的是，盡管本文有限元數(shù)值計算僅考慮了材料線彈性，但結(jié)合極限平衡法和接觸算法，通過選擇合適的滑動面，基于現(xiàn)行規(guī)范，可以應(yīng)用于大多數(shù)巖石高邊坡的穩(wěn)定性評價。