闞露 張發(fā)明

摘要：文章以位于某水電站壩址區(qū)右岸古滑坡Ⅱ區(qū)堆積體為研究對象，建立了一個較為復雜的三維模型，并從主應力、塑性區(qū)、剪應力和位移四個方面對古滑坡堆積體在蓄水條件下的變形和穩(wěn)定性情況進行了詳細的模擬和分析。計算得到的四個方面的結果都相互吻合相互印證。文章以FLAC3D中的單元界面能自動依附于指定范圍內模型表面生成的特性為思路來模擬地下水位。研究結果表明，堆積體雖然未出現(xiàn)整體性滑動破壞，但要注意坡腳局部地區(qū)坍塌和土層內部局部滑弧滑動，具體直觀地在模型上指出堆積體上的不穩(wěn)定區(qū)域，為后期針對性的加固治理方案提供參考，為工程優(yōu)化設計提供借鑒意義。

關鍵詞：堆積體；三維數(shù)值模擬； 穩(wěn)定性

中圖分類號：文獻標志碼：A

文章編號：1672-1098（2013）03-0000-00

崩塌、滑坡發(fā)生后形成的堆積體變形及穩(wěn)定性問題一直是工程地質領域研究的重點問題之一，滑坡形成的堆積體結構松散，強度低，在受到擾動后容易失穩(wěn)，對水電站的水庫正常運營會造成很大的影響，尤其是大型堆積體如果發(fā)生破壞，很可能會堵塞河道，后果不堪設想。迄今為止，研究者們已提出了如極限平衡分析法、塊體理論法[1]、數(shù)值方法[2-3]等計算條件相對簡單的傳統(tǒng)方法，也提出了人工神經網絡理論[4]、模糊數(shù)學理論[5]、灰色系統(tǒng)理論[6]、可靠性理論[7]等基于未確知性理論的計算方法。

本文采用FLAC3D數(shù)值模擬的方法對某水電站壩址區(qū)右岸古滑坡Ⅱ區(qū)堆積體進行穩(wěn)定性分析， 由于FLAC3D采用的是顯式差分法求解微分方程，不形成剛度矩陣，每一步計算僅需要很小的內存。在求解過程中通過疊加每一時步的小變形獲得大變形求解，從而僅占用計算機很少的內存就可以模擬大變形，并且在模擬材料的塑性破壞和塑性流動方面具有優(yōu)勢。但是由于FLAC前處理功能較弱，建立復雜的三維模型非常困難，故本文選用ANSYS進行堆積體前期建模，并劃分計算網格，再導入FLAC3D程序中計算，充分利用了兩種軟件的優(yōu)點，優(yōu)化分析過程。

2計算原理

21有限差分原理

有限差分法就是在利用數(shù)值計算方法求解偏微分方程時，用有限差分近似公式代替每一處導數(shù)，從而將求解偏微分方程的問題轉化成求解代數(shù)方程的問題。

FLAC有限差分公式：由高斯散度定理有

∫snids=∫Afxi （1）

式中：∫s為封閉曲面上沿邊界的積分； ni為曲面s的單位法向量；f為標量、向量或張量； xi為位置向量； ds為增量弧長；∫A為對表面積A積分。

定義f在面A上的梯度平均值為：

=1A∫AfxidA （2）

式中：<>表示求平均值。

將（2）代入（1）后得到

=1A∫sni fds （3）

對一個三角形子單元

=1A∑sniΔs （4）

式中：s為三角形某一條邊的邊長，等式右邊的求和在三角形三條邊上加總。取為相應邊上的平均值[8]。

應力應變：用每一邊速度矢量均值ui代替（4）式中的f，ui取各條邊兩端點的結點（即差分網格的角點）a和b的速度平均值，則

=12A∑N[（uai+ubi）niΔs]≈uixj （5）

同理可以求出的值。由幾何方程就可以求出其應變率

eij=12uixj+ujxi （6）

由材料的本構方程和相應的邊界條件，就可以求得應力增量。對各向同性的材料，有

σij=λδijθ+3μeij （7）

式中：λ，μ為拉梅常數(shù)；θ為體積應變，當時i=j，δij為1，否則，δij為零。

這樣通過上述各式的迭代求解，就可以得出每一迭代時步對應的單元的應力應變值[9]。

22強度折減法

本文通過強度折減法來計算堆積體的穩(wěn)定系數(shù)。強度折減法的原理是通過對巖土體的強度指標C和Φ值不斷地折減，反復計算，直到塑性區(qū)貫通整個坡體，即邊坡達到了臨界的破壞狀態(tài)，此時的折減系數(shù)就作為坡體的穩(wěn)定系數(shù)Fs。

Cf=CFtrial （8）

Φf=tan-1（tanΦFtrial） （9）

式中：Cf為折減后的粘結力；Φf為折減后的摩擦角；Ftrial為折減系數(shù)。

3堆積體基本特征

古滑坡堆積體位于貴州某水電站選壩河段，分為四個區(qū)域，其中Ⅱ區(qū)巖土體結構最為復雜，控制底界以上物質主要以崩坡積塊碎石土組成，故以Ⅱ區(qū)作為本文研究對象。堆積體所在河段河谷為不對稱的“V”型河谷，總體為順向坡，呈明顯圈椅狀地貌，順河長600～700 m，橫河寬約730 m，前緣臨河，后緣高程為740～790 m。順河方向上呈上升臺階狀，上游側地勢低于下游側40～70 m?？傮w上地形坡度約25°。后緣緊靠灰?guī)r陡壁，區(qū)域內發(fā)育兩條淺沖溝。覆蓋層物質組成為崩坡積碎塊石土為主。下伏基巖主要為J1夾層以下的二疊系下統(tǒng)棲霞組（P1q）中厚層夾薄層含炭質灰?guī)r。

圖1Ⅱ區(qū)地質分區(qū)平面示意圖及13條剖面線位置

J1夾層為該區(qū)早期層狀巖體失穩(wěn)破壞的控制底界，夾層以上層狀座滑體已經完全剝蝕搬運走，目前該區(qū)上部覆蓋物質來源應為后緣陡壁后期崩塌堆積形成。平面示意圖和地質橫剖面示意圖分別見圖1和圖2。

圖2Ⅱ區(qū)典型地質橫剖面示意圖

4Flac3D數(shù)值模擬

41幾何模型

模型邊坡的主滑方向為S14°E，走向N76°E。模型以主滑方向的反方向為X軸正方向，豎直向上為Y軸正方向，走向NE方向為Z軸正方向。豎直方向坐標采用實際高程坐標，建立直角坐標系，如圖3 所示。根據(jù)現(xiàn)場提供的地質資料，Ⅱ區(qū)堆積體模型選取了13個剖面，分別是Z=0，Z=57，Z=128，Z=215，Z=272，Z=325，Z=380，Z=438，Z=501，Z=561，Z=627，Z=697，Z=786（剖面線位置見圖1），左邊界略包括了Ⅰ區(qū)部分巖體。模型邊坡底面高程190 m，坡頂高程850 m，沿走向方向長度為786 m，沿傾向方向長度為1400 m。模型在X軸上的范圍是-200～1200 m，在Y軸上的范圍是190～850 m，在Z軸上的范圍是0～786 m。

圖3Ⅱ區(qū)堆積體三維模型網格圖

圖4z=750 m剖面靜水壓力云圖

模型分為5個部分，如圖3，對應圖上顏色，1，5部分為基巖，2為j1夾層，3為Ⅰ區(qū)座滑層狀巖體，4為Ⅱ區(qū)表面堆積體，其中1，3，4，5部分選用solid45單元類型，夾層選用shell99單元類型。Ⅱ區(qū)堆積體模型共有15569個單元，4969個節(jié)點。計算模型邊界約束形式為：側邊界只對水平方向進行約束，底邊界在水平和豎直方向都進行約束，模型的上部邊界取為自由面。FLAC程序中，節(jié)點速度是主要變量，所以選取模型的邊界條件是通過約束模型邊界的節(jié)點速度實現(xiàn)的，即模型底部邊界的水平、豎直方向的速度約束和四周邊界水平方向的速度約束。在程序中表達為邊界上xvel，yvel設置0. 根據(jù)實際地勘報告，計算采用的各層物理力學參數(shù)如表1所示。

表1巖、土體物理力學指標

地層巖性密度ρ（g/cm3）天然飽和泊松比抗剪斷（巖/巖）（°）c′/MPa飽和（°）飽和c′彈性模量/GPa

基巖26526802445110//200

座滑層狀巖體230232022330401//05

堆積體19522002626572024232001

j1泥巖夾層195230025167003//035

地下水面的生成：由于在FLAC3D中直接生成符合勘察資料的地下水面比較困難，所以本文中建立模型時充分利用FLAC3D中的單元界面能自動依附于指定范圍內模型表面生成的特性，生成水面。正常河水位為383 m，水庫正常蓄水位為400 m。根據(jù)勘察報告，本文采用地下水位從坡腳處高程400 m開始，以5°向上傾斜，來近似模擬實際地下水位形態(tài)，如圖4。沿著水位線，對三維模型采用切割命令，得到地下水面及其下部的幾何實體，并剖分網格，導入FLAC3D，并以之為輔助單元，生成與地下水面空間形態(tài)一致的單元界面。采用FISH語言[10]遍歷界面單元節(jié)點，生成水面，同時生成靜水壓力。

42計算結果分析

421主應力場規(guī)律分析

計算表明，坡體的最大主應力及最小主應力基本為層狀分布，并在坡腳高程約400～460 m處出現(xiàn)應力集中，坡體大主應力σ1最大值為18139 MPa，小主應力值σ3最大值為10585 MPa。具體應力分布見圖5及圖6。

圖5坡體大主應力σ1分布圖

圖6坡體小主應力σ3分布圖

從剖面Z=650 m的應力分布圖（如圖7）可以看出，剖面附近的最大主應力（壓應力），基本順著坡面方向，并一直延伸到坡腳。而往邊坡內部，最大主應力方向與水平軸的夾角逐步變大，直至鉛直；由于巖層分界面的存在，使得其附近區(qū)域的最大主應力方向要比其他區(qū)域最大主應力方向變化大而且迅速，但并未影響主應力分布的總體走勢。這些都表明邊坡深部土體主要受鉛垂方向的壓應力作用，體現(xiàn)為受壓屈服。

圖7Z=650 m剖面主應力分布圖

圖8坡體塑性區(qū)分布圖

422塑性區(qū)分布規(guī)律分析

從圖8來看，坡體的屈服區(qū)域較集中，范圍較小。坡體后緣（高程約780～850 m）分布有小塊拉張塑性區(qū)，坡腳處（高程約380～450 m）分布零散的剪切塑性區(qū)。堆積體發(fā)生張拉剪切破壞的可能性很小，即使發(fā)生，也僅是局部區(qū)域，不會對堆積體整體穩(wěn)定性造成重大影響。從剖面塑性區(qū)分布圖（圖9，圖10）來看，僅僅在坡體的后緣和坡腳處出現(xiàn)了零星的塑性區(qū)，這表明堆積體處于正常的工作狀態(tài)。

圖9Z=210 m剖面塑性區(qū)分布圖

圖10Z=750 m剖面塑性區(qū)分布圖

坡腳處的剪切塑性區(qū)位于表面堆積體內，厚度大約5～10 m。需要強調的是，計算結果顯示的是以Mohr-Coulomb屈服準則為依據(jù)的塑性區(qū)分布情況。該屈服準則則認為材料進入屈服即破壞。

423剪應力（剪應變增量）規(guī)律分析

判斷堆積體的（潛在）滑動面（帶），可根據(jù)其剪應力（應變增量）來判斷；剪應力較為集中或剪應變增量較大（絕對值）的部位，則為其（潛在）滑動面（帶）。變形破壞也都沿此處發(fā)生；剪應力分布較為分散（均勻）或者剪應變增量較小或基本上沒發(fā)生變化的部位，一般不會有潛在滑動面產生，因此，這些部位也不會發(fā)生較大的變形和破壞。應力集中帶主要出現(xiàn)在3個范圍內，坡體上下游側坡腳處和j1夾層附近，它們是堆積體最有可能發(fā)生破壞的部位，見圖11。J1夾層附近出現(xiàn)剪應變增量增高帶，說明坡體有可能的破壞模式是沿著以j1夾層附近為滑面進行滑動，但由于計算得到的堆積體整體穩(wěn)定系數(shù)Fs為110，所以發(fā)生整體性滑動的可能性不大；坡體上下游側坡腳處的兩個剪應變增量增高帶，說明坡體可能沿土層內部局部滑弧滑動。Z=730剪應變增量剖面見圖12。

圖11整體剪應變增量展示圖

圖12Z=730 m剪應變增量剖面展示圖

424位移場規(guī)律分析

坡體變形主要限于堆積體中、下部（600 m高程以下）區(qū)域。豎直方向位移總體表現(xiàn)為下沉，最大下沉區(qū)與水平方向最大位移區(qū)一致，局部沉降最大量值為51368 cm，其他下沉部位量值為0～1 cm（如圖13）。

圖13坡體豎直方向位移云圖

出現(xiàn)局部最大沉降的區(qū)域有兩處，第一處位于Ⅱ區(qū)下游側，高程（Y坐標）范圍大約位于420～425 m之間，走向方向（Z坐標）范圍大約處于130～150 m之間，表面積小于100 m2；第二處位于Ⅱ區(qū)上游側，高程（Y坐標）范圍大約位于410～415 m之間，走向方向范圍（Z坐標）范圍大約處于760～770 m之間，表面積小于50 m2。水平方向局部最大位移量值為13.483 cm，出現(xiàn)區(qū)域與豎直方向局部最大沉降區(qū)一致，坡體其他部位的水平方向位移量值大都處于0～2 cm之間（如圖14）。坡體最大變形區(qū)位置與塑性區(qū)位置相吻合。endprint

圖14 坡體水平方向位移云圖

5結語

1） 根據(jù)強度折減法，得到的堆積體整體穩(wěn)定系數(shù)Fs為110，所以發(fā)生整體性滑動的可能性不大。文章對滑坡主應力，塑性區(qū)，剪應力和位移四個方面的計算結果進行了詳細分析，并且四個方面的計算結果都相互吻合相互印證。

2） 可以直觀地在模型上看出在蓄水條件下，上下游坡腳處的碎石土堆積層是Ⅱ區(qū)的主要不穩(wěn)區(qū)域，具體在兩個位置：①下游側坡腳處，高程（Y坐標）范圍大約位于420～425 m之間，走向方向（Z坐標）范圍大約處于130～150 m之間，②上游側坡腳處，高程（Y坐標）范圍大約位于410～415 m之間，走向方向范圍（Z坐標）范圍大約處于760～770 m之間，對后期堆積體的加固治理方案提供了直觀詳細的位置。

3） 本文在模擬蓄水條件下的地下水位時利用了FLAC3D中的單元界面能自動依附于指定范圍內模型表面生成的特性，通過建立輔助單元生成水面并同時生成靜水壓力，思路清晰明了，適用于比較復雜的坡體模型。

4） 利用FLAC3D軟件對堆積體進行穩(wěn)定性分析具有方便，經濟，直觀的優(yōu)勢，但是由于建立的模型不可能完全模擬出現(xiàn)實堆積體的所有特性，選取的計算參數(shù)也不可能完全與實際相符，所以需要結合其他技術手段綜合評價，才能得出更為全面正確的結果。

參考文獻：

[1]李素梅. 塊體理論在巖石路塹邊坡穩(wěn)定分析中的應用[J]. 云南交通科技， 2003， 19（4）：14-16.

[2]張均峰， 丁樺.邊坡穩(wěn)定性分析的三維極限平衡法及應用[J]. 巖石力學與工程學報， 2005，24（3）：365-370.

[3]楊強， 朱玲， 薛利軍. 基于三維多重網格法的極限平衡法在錦屏高邊坡穩(wěn)定性分析中的應用[J]. 巖石力學與工程學報， 2005， 24（增2）：5 313-5 318.

[4]陳昌彥， 王思敬，沈小克. 邊坡巖體穩(wěn)定性的人工神經網絡預測模型[J]. 巖土工程學報， 2001， 23（2）：157-161.

[5]賈厚華， 賀懷建. 邊坡穩(wěn)定模糊隨機可靠度分析[J]. 巖土力學， 2003， 24（4）：657-660.

[6]陳新民， 羅國煜. 基于經驗的邊坡穩(wěn)定性灰色系統(tǒng)分析與評價[J]. 巖土工程學報， 1999， 21（5）：638-641.

[7]祝玉學. 邊坡可靠性分析[M]. 北京：冶金工業(yè)出版社， 1993.

[8]黃潤秋，許強. 顯式拉格朗日差分分析在巖石邊坡工程中的應用[J]. 巖石力學與工程學報， 1995， 14（4）：346-354.

[9]鄒棟， 鄭宏. 快速拉格朗日法及其在邊坡穩(wěn)定性分析中的應用[J]. 礦業(yè)研究與開發(fā)， 2005， 25， （5）：80-83.

[10]陳育民， 徐鼎平. FLAC/FLAC3D基礎與工程實例[M]. 北京： 中國水利水電出版社， 2008.

（責任編輯：）endprint