錦屏一級水電站猴子崖邊坡三維非連續(xù)變形數(shù)值模擬

秦永濤 宋慧 游湘 楊 濤

（1 中國水電顧問集團(tuán)成都勘測設(shè)計研究院，四川 成都 610072；2 西南交通大學(xué)土木工程學(xué)院，四川 成都 610031）

1 引 言

邊坡穩(wěn)定性是巖土工程學(xué)中的一大類問題，國內(nèi)外眾多學(xué)者對此都做過研究，應(yīng)運(yùn)而生了許多有價值的研究成果，對此很多人做過論述[1]～[6]，為邊坡工程的設(shè)計、施工起到了很好的指導(dǎo)作用。近幾十年來,隨著大型水利水電工程、公路工程、鐵路遂道工程和采礦工程等的迅速發(fā)展, 巖土力學(xué)作為巖土工程基礎(chǔ)學(xué)科得到了較快的發(fā)展。特大型的工程實(shí)踐為巖土力學(xué)的發(fā)展賦予了巨大的推動力,同時也產(chǎn)生了許多復(fù)雜的巖土力學(xué)新課題。特別是西部大開發(fā)以來，西南地區(qū)水利水電工程、公路鐵路基礎(chǔ)設(shè)施建設(shè)中遇到的大量復(fù)雜地質(zhì)情況、高地應(yīng)力、高陡巖質(zhì)邊坡和地下洞室群問題，使得對巖土工程的認(rèn)識有了更深的要求，出現(xiàn)了更多的工程技術(shù)難題，更多的現(xiàn)有理論需要進(jìn)一步完善。

解決巖土力學(xué)問題的方法主要有實(shí)驗(yàn)方法、理論分析方法和數(shù)值模擬方法三大類。這三類方法相輔相成,互為補(bǔ)充。與其它兩種方法相比，數(shù)值模擬是解決巖土工程問題較有效的手段,已被學(xué)術(shù)界和工程界廣泛接受,作為一種力學(xué)狀態(tài)的分析工具,它越來越多的應(yīng)用于巖土體的穩(wěn)定性分析、巖土工程設(shè)計和巖土工程基本問題分析中。但是前面兩種方法是數(shù)值分析的前提和基礎(chǔ)，數(shù)值模擬是對實(shí)驗(yàn)結(jié)果和理論分析成果的佐證。

幾十年以來,巖土力學(xué)數(shù)值分析方法得到了迅速發(fā)展, 出現(xiàn)了有限單元法、離散元法、有限差分法、邊界元法、無限元法、拉格朗日元法、非連續(xù)變形分析方法、無單元法、流形元法及其耦合的數(shù)值計算方法和以數(shù)值模擬為主的漸進(jìn)破壞模型等[7]～[14]數(shù)值分析技術(shù)。其中基于連續(xù)介質(zhì)大變形分析的拉格朗日元法在實(shí)際工程中也得到了較好的應(yīng)用。拉格朗日元法避免了有限元法進(jìn)行大型矩陣的復(fù)雜計算，可以同時考慮巖土體的材料非線性和幾何非線性,并能跟蹤物體變形的全過程，適于分析巖土力學(xué)中的大變形問題。本文運(yùn)用連續(xù)介質(zhì)大變形分析理論計算了錦屏一級水電站右岸下游河道霧化區(qū)猴子崖邊坡(簡稱“猴子坡”)在天然狀況和泄洪霧化雨作用下的應(yīng)力應(yīng)變變化及位移變化趨勢，評價了邊坡的穩(wěn)定性。

2 連續(xù)介質(zhì)大變形理論

三維快速拉格朗日（3D-FLAC）分析是一種基于顯式有限差分法的數(shù)值分析方法,是非連續(xù)變形分析方法中的一種。該法將計算區(qū)域劃分為若干六面體單元，單元網(wǎng)格可以隨著材料的變形而變形，即所謂的拉格朗日算法。這種算法基本原理類同于離散單元法，但它卻能像有限元那樣適用于多種材料模式與邊界條件的非規(guī)則區(qū)域的連續(xù)問題求解。還能針對不同的材料特性，使用相應(yīng)的本構(gòu)方程來比較真實(shí)地反映實(shí)際材料的動態(tài)行為。可以準(zhǔn)確地模擬材料的屈服、塑性流動、軟化直至大變形,尤其在材料的彈塑性分析、大變形分析以及模擬施工過程等領(lǐng)域有其獨(dú)到的優(yōu)點(diǎn)[15]。

2.1 空間導(dǎo)數(shù)的有限差分近似

快速拉格朗日分析采用混合離散方法，將區(qū)域離散為常應(yīng)變六面體單元的集合體，又將每個六面體看作以六面體角點(diǎn)為角點(diǎn)的常應(yīng)變四面體的集合體，應(yīng)力、應(yīng)變、節(jié)點(diǎn)不平衡力等變量均在四面體上進(jìn)行計算，六面體單元的應(yīng)力、應(yīng)變?nèi)≈禐槠鋬?nèi)四面體的體積加權(quán)平均。如一四面體，節(jié)點(diǎn)編號為1到4，第n面表示與節(jié)點(diǎn)n相對的面，設(shè)其內(nèi)一點(diǎn)的速率分量為vi，由高斯公式得

其中V為四面體的體積，S為四面體的外表面，nj為外表面的單位法向量分量。

對于常應(yīng)變單元，vi為線性分布，nj在每個面上為常量，由式（1）可得

式中上標(biāo)l表示節(jié)點(diǎn)l的變量，（l）表示面l的變量。

2.2 運(yùn)動方程

快速拉格朗日分析以節(jié)點(diǎn)為計算對象，在時域內(nèi)求解。節(jié)點(diǎn)運(yùn)動方程為

式中，F(xiàn)li(t)為在t時刻l節(jié)點(diǎn)在i方向的不平衡力分量，可由虛功原理導(dǎo)出。ml為l節(jié)點(diǎn)的集中質(zhì)量，對于靜態(tài)問題，采用虛擬質(zhì)量以保證數(shù)值穩(wěn)定，而對于動態(tài)問題則采用實(shí)際的集中質(zhì)量。

將（3）式左端用中心差分來近似，則可得

2.3 應(yīng)變、應(yīng)力及節(jié)點(diǎn)不平衡力

快速拉格朗日分析由速率來求某一時步的單元應(yīng)變增量，即

有了應(yīng)變增量，即可由本構(gòu)方程求出應(yīng)力增量，進(jìn)而得到總應(yīng)力。

2.4 阻尼力

對于靜態(tài)問題，在式（3）的不平衡力中加入了非粘性阻尼，以使系統(tǒng)的振動逐漸衰減直至達(dá)到平衡狀態(tài)（即不平衡力接近零）。此時式（3）變?yōu)?/p>

阻尼力為

式中 α為阻尼系數(shù)，其默認(rèn)值為0.8，而

3 工程概況

錦屏一級水電站位于四川省涼山彝族自治州鹽源縣和木里縣境內(nèi)，是雅礱江干流水能資源最富集的中、下游河段五個梯級水電開發(fā)的第一級。壩址位于普斯羅溝與手爬溝之間1.5km長的河段上，河流流向N25°E，河道順直而狹窄。電站以發(fā)電為主，兼有防洪、攔沙等作用。壩址控制流域面積10256km2，多年平均流量1200m3/s。水庫正常蓄水位1880m，攔河大壩為混凝土雙曲拱壩，壩高305m，水庫庫容77.6億m3，為年調(diào)節(jié)水庫，電站裝機(jī)容量3600MW，年發(fā)電量166.2億kW·h。

猴子坡為壩區(qū)右岸Ⅳ線下游與Ⅲ線之間一三面臨空的凸出山梁。地貌上該斜坡呈SEE向延伸的山脊?fàn)?，脊?fàn)钚逼碌腟W側(cè)和NE側(cè)均呈陡坎，下部較窄。山脊走向與河流近于垂直，臨江坡度50°～80°。坡體由雜谷腦組第二段第6、7層大理巖組成，第6層為薄～中厚層狀條帶狀大理巖，局部夾綠片巖夾層，第7層為厚層條紋狀大理巖；巖層產(chǎn)狀N40°～55°E，NW∠35°～40°。如圖1所示。

圖1 猴子坡數(shù)值模擬區(qū)域

該邊坡距大壩約600～780m，根據(jù)泄洪霧化模型試驗(yàn)報告，猴子崖邊坡處于右岸泄洪霧化區(qū)邊坡范圍內(nèi)，加上有多種不利結(jié)構(gòu)面組合，需高度重視在泄洪霧化降雨作用下的邊坡穩(wěn)定性分析。

4 三維數(shù)值模擬分析

4.1 三維數(shù)字模型的建立

根據(jù)右岸猴子巖山梁區(qū)地形測繪數(shù)據(jù)和地質(zhì)資料建立三維幾何模型，如圖2。模型采用局部坐標(biāo)系（X、Y、Z），其X軸正方向?yàn)镹E28°，指向河谷下游，Y軸正方向?yàn)镹W62°，指向左岸山里，Z軸正方向鉛直向上。幾何模型建立范圍為：上下游X方向，以猴子巖山梁山脊線為基線分別向河流上下游方向各延伸330m和500m，總長830m；垂直河流方向，從河谷向山內(nèi)側(cè)延伸550m；高程方向，從1480m高程一直到自然坡面，最大高程至2345m。

三維幾何模型中真實(shí)再現(xiàn)了計算區(qū)域復(fù)雜的地形地貌、控制性結(jié)構(gòu)面。對猴子坡山體而言，控制性結(jié)構(gòu)面包括：f7、f28、g37-1～g37-5、近SN向結(jié)構(gòu)面。f7、f28的產(chǎn)狀及位置均確定，如圖2所示。對于層間擠壓錯動帶而言，該區(qū)發(fā)育g37-1～g37-5等多條層間擠壓帶，而計算模型中不可能全部真實(shí)考慮，只能進(jìn)行一定程度的概化。根據(jù)三維剛體極限平衡分析結(jié)果，g37-2作為底滑面的潛在失穩(wěn)塊體穩(wěn)定性最低、塊體體積最大，最可能失穩(wěn)，是該區(qū)的控制性結(jié)構(gòu)面，因此有限元計算中采用g37-2的產(chǎn)狀建立層間擠壓錯動帶，而將層間擠壓帶的厚度設(shè)為10m來綜合反映其它各條錯動帶對邊坡穩(wěn)定性的影響。近SN向結(jié)構(gòu)面可能構(gòu)成猴子坡的下游側(cè)邊界，其產(chǎn)狀對邊坡失穩(wěn)范圍的確定十分重要，計算中采用遍布節(jié)理模型來模擬，節(jié)理模型的產(chǎn)狀采用近SN向結(jié)構(gòu)面的產(chǎn)狀。

圖2 猴子坡數(shù)值計算模型

在幾何模型的基礎(chǔ)上，建立有限元計算模型。由于猴子坡地形復(fù)雜、結(jié)構(gòu)面眾多、產(chǎn)狀變化大，為較好的模擬這些地質(zhì)因素，模型全部采用四面體單元進(jìn)行離散。計算域共剖分了762612個單元，134073個節(jié)點(diǎn)。剖分后的有限元網(wǎng)格圖如圖3所示。對于猴子坡而言，f7、f28、g37～2等結(jié)構(gòu)面是特別關(guān)心的，因此進(jìn)行了重點(diǎn)離散，共離散為62906個單元，如圖4所示。

圖3 模型網(wǎng)格劃分

圖4 結(jié)構(gòu)面網(wǎng)格

4.2 計算參數(shù)和邊界條件

根據(jù)計算域的地質(zhì)情況，按巖體分級分區(qū)選擇材料，巖體采用遍布節(jié)理模型模擬，節(jié)理的產(chǎn)狀設(shè)為近SN向結(jié)構(gòu)面的產(chǎn)狀，f7、f28斷層、層間擠壓錯動帶g37～2等采用莫爾－庫侖材料模擬。材料參數(shù)見表1。

表1 計算中采用的巖體力學(xué)參數(shù)

近SN向結(jié)構(gòu)面的抗剪強(qiáng)度參數(shù)按照20%取B1類指標(biāo)，80%取IV1類巖體指標(biāo)綜合加權(quán)確定。

根據(jù)模型所處的地形地貌條件、河谷對稱性及邊坡載荷方向，模型邊界順河向邊界(X向)、橫河向邊界(Y向)和底部邊界(Z向)分別取法向支座。

4.3 計算工況

計算分別考慮2種工況：

工況1：天然狀態(tài)下，猴子坡山梁的應(yīng)力狀態(tài)及其穩(wěn)定性分析；

工況2：泄洪霧化狀態(tài)下，猴子坡山梁的應(yīng)力狀態(tài)及其穩(wěn)定性分析；

其中，對工況2計算中，考慮最不利的情況，高程1830m以下的巖體均達(dá)到飽和狀態(tài)，巖體容重取飽和容重(γ＝29kN/m3)，在上游側(cè)邊界(－X)和山內(nèi)側(cè)(－Y)1830m高程以下施加隨高度變化(梯度為10kN/m)的三角形荷載，以模擬霧化雨形成的水力邊界，如圖5所示。

圖5 霧化雨形成的水力邊界示意圖

4.4 計算成果分析

(1) 天然狀態(tài)的計算成果

坡體在自重應(yīng)力場作用下的應(yīng)力場分布如圖6～圖11所示。其中，圖6為猴子巖附近區(qū)域的最大主應(yīng)力等值線（拉為正），圖7為坡體表面的拉應(yīng)力區(qū)，為分析f28、f7及g37～2等結(jié)構(gòu)面上的應(yīng)力特征，圖8示出了邊界面上的最大主應(yīng)力云圖。圖9為坡體表面的最小主應(yīng)力等值線。圖10示出了坡面上的塑性屈服區(qū)。

由圖分析可知，猴子坡潛在失穩(wěn)塊體范圍明顯。在自重作用下，由于猴子巖山梁為一突出山梁，其上游和低高程（河谷處）由f7和g37～2所圍限。在f7和g37～2交界部位，各結(jié)構(gòu)面在坡表面均表現(xiàn)為明顯的拉應(yīng)力區(qū)，即使是從最小主應(yīng)力圖來看，這部分坡體表面也呈現(xiàn)出極小的壓應(yīng)力狀態(tài)，表明猴子巖失穩(wěn)塊體可能沿此界限失穩(wěn)下滑。從結(jié)構(gòu)面的最大主應(yīng)力云圖來看，f7和f28淺層均出現(xiàn)拉應(yīng)力區(qū)，其下部和g37～2上游側(cè)則表現(xiàn)為壓應(yīng)力，但應(yīng)力水平較低，而g37～2下游側(cè)壓應(yīng)力水平較高，等值線呈由河谷向山里的展布狀態(tài)，也揭示了其上塊體的失穩(wěn)趨勢。

從塑性區(qū)分布來看，猴子坡低高程潛在失穩(wěn)塊體存在失穩(wěn)趨勢，安全裕度不高，但失穩(wěn)塊體還沒有完全形成。坡面屈服區(qū)僅出現(xiàn)在f7、f28和g37～2的高程局部處，并沒有貫穿整個g37～2結(jié)構(gòu)面。結(jié)構(gòu)面上的屈服區(qū)分布更能體現(xiàn)這一特征，f7和f28的絕大部分均出現(xiàn)了塑性屈服，表明邊坡塊體有失穩(wěn)滑移趨勢，后緣邊界明顯。但g37～2僅坡面附近出現(xiàn)塑性屈服，其分布范圍并不大，則又表明塊體還沒有完全形成，還不可能沿此底滑面出現(xiàn)整體失穩(wěn)滑移。由于兩側(cè)邊界的約束作用以及坡腳巖體的抗滑作用，邊坡仍能保持整體穩(wěn)定。但是由于猴子巖山梁底滑面前部（層間擠壓錯動帶）和上游側(cè)邊界（f7斷層）基本上已經(jīng)屈服，一旦山梁后緣（近SN向結(jié)構(gòu)面）形成拉裂縫，容易發(fā)生順層的滑移拉裂破壞，安全儲備不高。

圖6 最大主應(yīng)力等值線

圖7 拉應(yīng)力區(qū)

圖8 邊界面上的最大主應(yīng)力

圖9 最小主應(yīng)力等值線

(2) 霧化雨作用下的計算成果

圖11～圖13所示為霧化雨作用下的坡體應(yīng)力，圖14示出了坡體的位移場變化，圖15為坡體結(jié)構(gòu)面上的塑性屈服區(qū)。

在霧化雨作用下，坡體應(yīng)力場進(jìn)一步惡化，主要表現(xiàn)在拉應(yīng)力的量值增大，作用范圍集中而貫通，塊體失穩(wěn)邊界逐步形成。猴子坡山梁上游側(cè)（f7通過區(qū)域）拉應(yīng)力分布明顯，整個區(qū)域均呈現(xiàn)出受拉狀態(tài)，沿g37～2向下游，拉應(yīng)力區(qū)集中而呈條帶狀（如圖12），表明g37～2作為上游側(cè)邊界、g37～2作為底滑面的失穩(wěn)范圍十分明確，塊體已經(jīng)沿此出現(xiàn)了拉裂破壞。而從河谷應(yīng)力來看（如圖13），河谷出現(xiàn)了極大的壓應(yīng)力集中，勢必造成巖體屈服，因此，河谷側(cè)的壓應(yīng)力區(qū)也揭示了猴子坡失穩(wěn)塊體的范圍所在。從結(jié)構(gòu)面的最大主應(yīng)力云圖來看，除f7和f28淺層均出現(xiàn)拉應(yīng)力區(qū)外，較之天然狀態(tài)，結(jié)構(gòu)面上的壓應(yīng)力水平大為降低， f7和f28下部和g37～2上游側(cè)壓應(yīng)力水平大為降低，幾乎接近0應(yīng)力狀態(tài)，考慮到這些結(jié)構(gòu)面位于塊體的上游側(cè)，為失穩(wěn)塊體的后緣位置，這些部位出現(xiàn)0應(yīng)力狀態(tài)表明塊體具有明顯的失穩(wěn)下滑趨勢。

在霧化雨作用下，坡體發(fā)生了明顯變形，變形區(qū)域集中在猴子坡前部，揭示了潛在失穩(wěn)區(qū)域。如圖14所示，猴子巖前部發(fā)生了明顯的位移變化，最大位移達(dá)到0.055m，向山里逐漸遞減，位移等值線方向與近SN向結(jié)構(gòu)面的走向基本一致，表明近SN向的結(jié)構(gòu)面是失穩(wěn)塊體的下游邊界面。若以位移大于0.01m的范圍作為失穩(wěn)塊體范圍，則限定的失穩(wěn)塊體范圍如圖20所示?？梢娛Х€(wěn)塊體的邊界為f7、g37～2及近SN向結(jié)構(gòu)面，而f28由于其出露高程較高，沒有形成失穩(wěn)塊體的邊界。

圖15示出了坡體結(jié)構(gòu)面上的塑性區(qū)分布，可見，結(jié)構(gòu)面上的塑性區(qū)大為發(fā)展，g37～2以上的坡體也沿節(jié)理方向發(fā)生了塑性流動，由f7、g37～2和近SN向結(jié)構(gòu)面所圍限的巖體全部發(fā)生了塑性流動，在臨空面具備的情況下，塊體極有可能出現(xiàn)失穩(wěn)下滑。

圖12 拉應(yīng)力區(qū)

圖13 最小主應(yīng)力等值線

圖14 位移大于0.01m的區(qū)域

圖15 結(jié)構(gòu)面上的屈服區(qū)

5 結(jié) 論

猴子坡巖體可能的失穩(wěn)破壞模式為：f7斷層、層間擠壓錯動帶、近SN向結(jié)構(gòu)面形成楔形體組合，有限元計算中，考慮了這些結(jié)構(gòu)面的影響，f7斷層和層間擠壓錯動帶通過建立單獨(dú)的巖體材料反映，近SN向結(jié)構(gòu)面為隨機(jī)分布，采用遍布節(jié)理模型加以反映，重點(diǎn)計算分析了天然邊坡和霧化雨邊坡的穩(wěn)定性。

由于猴子坡的特殊坡體結(jié)構(gòu)，應(yīng)力場特征揭示，在天然邊坡條件下即具有較明顯的失穩(wěn)塊體范圍，其范圍與地質(zhì)判斷一致。塑性區(qū)分布揭示，猴子坡低高程潛在失穩(wěn)塊體存在失穩(wěn)趨勢，安全裕度不高，但失穩(wěn)塊體還沒有完全形成。

受霧化雨作用，坡體應(yīng)力場進(jìn)一步惡化，應(yīng)力場特征揭示，失穩(wěn)塊體邊界已經(jīng)逐步形成，通過位移場的判斷，可進(jìn)一步明確失穩(wěn)塊體的范圍，該范圍內(nèi)的巖體變形在進(jìn)一步發(fā)展中。結(jié)構(gòu)面和坡體節(jié)理方向塑性區(qū)大為擴(kuò)展，在臨空面具備的情況下，塊體具有明顯的失穩(wěn)下滑趨勢。

數(shù)值計算表明，猴子坡在霧化雨作用下可能會出現(xiàn)局部失穩(wěn)，需要對該部分巖體加強(qiáng)支護(hù)，做好防排水措施。

[1] 陶振宇. 試論巖石力學(xué)的最新進(jìn)展[J]. 力學(xué)進(jìn)展,1992 ,22 (2) :161～172.

[2] 常士驃.我國巖土工程回顧與展望[J].水文地質(zhì)工程地質(zhì),1993,20(1): 23～26.

[3] 哈秋聆.三峽工程中的若干力學(xué)問題[J].力學(xué)進(jìn)展,1994,24(4): 433～439.

[4] 楊志法.關(guān)于巖石力學(xué)當(dāng)前發(fā)展戰(zhàn)略的一些看法[J].巖土工程學(xué)報, 1994 ,16 (1) :112～113.

[5] 謝和平,劉夕才,王金安. 關(guān)于21 世紀(jì)巖石力學(xué)發(fā)展戰(zhàn)略的思考[J].巖土工程學(xué)報,1996 ,18 (4) :98～102.

[6] 傅冰駿. 國際巖石力學(xué)發(fā)展動向[J].巖石力學(xué)與工程學(xué)報,1997,16 (2) :195～196.

[7] Zienkiewicz O C. The finite element method.London: Mc-Graw-Hill, 1997.

[8] 王泳嘉,邢紀(jì)波. 離散單元法同拉格朗日元法及其在巖土力學(xué)中的應(yīng)用[J].巖土力學(xué), 1995, 16 (2):10～15.

[9] 王泳嘉.邊界元法在巖石力學(xué)中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報,1986,5 (2) :205～222.

[10] 應(yīng)隆安著.無限元方法[M]. 北京:北京大學(xué)出版社,1992.

[11] Shi G H. Numerical manifold method[M].Proceedings of the First International Conference on Analysis of Discontinuous Deformation, Chungli,Taiwan,1995:187～222.

[12] Shi G H. Discontinuous deformation analysis: A new numerical model for the static and dynamics of block systems. A dissertation submitted in partial satisfaction of the requirements for the degree of Doctor of Philosophy, Berkeley: University of California, 1988.

[13] Margulies M. Combination of the boundary element and finite element methods. Proceeding in Boundary Element Methods 1981(1): 258～288.

[14] Chang C C. Computational simulation of progressive fracture in fibre composites. Journal of Composite Materials, 1987 (21): 834～885.

[15] FLAC3D Version2.0 Online Manual[Z].Minneapolis: Itasca Consulting Group, Inc, 1997.