闞露 張發(fā)明

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

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

中圖分類號：文獻(xiàn)標(biāo)志碼：A

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

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

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

2計(jì)算原理

21有限差分原理

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

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

∫snids=∫Afxi （1）

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

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

=1A∫AfxidA （2）

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

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

=1A∫sni fds （3）

對一個(gè)三角形子單元

=1A∑sniΔs （4）

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

應(yīng)力應(yīng)變：用每一邊速度矢量均值ui代替（4）式中的f，ui取各條邊兩端點(diǎn)的結(jié)點(diǎn)（即差分網(wǎng)格的角點(diǎn)）a和b的速度平均值，則

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

同理可以求出的值。由幾何方程就可以求出其應(yīng)變率

eij=12uixj+ujxi （6）

由材料的本構(gòu)方程和相應(yīng)的邊界條件，就可以求得應(yīng)力增量。對各向同性的材料，有

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

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

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

22強(qiáng)度折減法

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

Cf=CFtrial （8）

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

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

3堆積體基本特征

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

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

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

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

4Flac3D數(shù)值模擬

41幾何模型

模型邊坡的主滑方向?yàn)镾14°E，走向N76°E。模型以主滑方向的反方向?yàn)閄軸正方向，豎直向上為Y軸正方向，走向NE方向?yàn)閆軸正方向。豎直方向坐標(biāo)采用實(shí)際高程坐標(biāo)，建立直角坐標(biāo)系，如圖3 所示。根據(jù)現(xiàn)場提供的地質(zhì)資料，Ⅱ區(qū)堆積體模型選取了13個(gè)剖面，分別是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ū)堆積體三維模型網(wǎng)格圖

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

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

表1巖、土體物理力學(xué)指標(biāo)

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

基巖26526802445110//200

座滑層狀巖體230232022330401//05

堆積體19522002626572024232001

j1泥巖夾層195230025167003//035

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

42計(jì)算結(jié)果分析

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

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

圖5坡體大主應(yīng)力σ1分布圖

圖6坡體小主應(yīng)力σ3分布圖

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

圖7Z=650 m剖面主應(yīng)力分布圖

圖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ū)位于表面堆積體內(nèi)，厚度大約5～10 m。需要強(qiáng)調(diào)的是，計(jì)算結(jié)果顯示的是以Mohr-Coulomb屈服準(zhǔn)則為依據(jù)的塑性區(qū)分布情況。該屈服準(zhǔn)則則認(rèn)為材料進(jìn)入屈服即破壞。

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

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

圖11整體剪應(yīng)變增量展示圖

圖12Z=730 m剪應(yīng)變增量剖面展示圖

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

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

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

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

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

5結(jié)語

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

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

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

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

參考文獻(xiàn)：

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

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

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

[4]陳昌彥， 王思敬，沈小克. 邊坡巖體穩(wěn)定性的人工神經(jīng)網(wǎng)絡(luò)預(yù)測模型[J]. 巖土工程學(xué)報(bào)， 2001， 23（2）：157-161.

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

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

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

[8]黃潤秋，許強(qiáng). 顯式拉格朗日差分分析在巖石邊坡工程中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào)， 1995， 14（4）：346-354.

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

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

（責(zé)任編輯：）endprint