基于FLAC3D方法對某邊坡的穩定性分析

李勇輝，孫 勇，朱伯文，唐江濤

(1.貴州大學 國土資源部喀斯特環境與地質災害重點實驗室，貴陽 5500033； 2.成都理工大學 地質災害與防治與地質環境保護國家重點實驗室，成都 610059)

西南地區主要是山區，在修建公路鐵路的過程中會遭遇大量的邊坡。而邊坡的穩定性分析是重中之重，首先要分析該邊坡的穩定與否，如若需要進行開挖，則需要對該開挖方案進行穩定性分析，從而判斷該方案的可靠度與可行性。近年來，數值模擬在邊坡穩定性分析中得到大量的運用，其中FLAC3D軟件因其模擬結果與實際情況更加相符，從而得到更廣泛的應用。

1 FLAC3D原理

FLAC3D是由美國itasca公司開發的一種模擬仿真計算軟件，是對二維有限差分程序FLAC2D的擴展，能夠對土體、巖石和其它材料進行三維結構的受力特性模擬和塑性流動分析。通過對三維網格中的多面體單元進行不斷調整來擬合實際的結構單元。單元材料主要有線性和非線性兩種本構模型，在外力作用下，當材料達到屈服流動后，網格可以產生相應的變形和移動(大變形模式)。FLAC3D主要采用顯式拉格朗日算法和混合-離散分區技術，能夠十分準確模擬材料的塑性流動和破壞。因為不用形成剛度矩陣，所以該軟件用較小的內存空間就可以求解相當大范圍的三維問題。FLAC3D主要有三大類模型：開挖模型、彈性模型和塑性模型。FLAC3D網格中的每個區域都可以賦予不同的材料模型，并且可以指定材料參數的統計分布以及變化梯度。而且模型中還包含有節理單元，也稱之為界面單元，可以模擬兩種以及多種材料界面，不同材料性質的間斷特性。節理允許發生滑動或者分離，因此能夠用于模擬巖體中的斷層、節理和摩擦邊界。FLAC3D中的網格生成器gen，可以通過匹配、連接在網格生成器中生成局部網格，從而方便地生成模擬所需要的三維結構網格。還可以自動產生交叉結構網格(比如說相交的巷道)，三維網格由整體坐標系X、Y、Z系統所確定，不同于FLAC程序是由行列方式確定。這就提供了比較靈活的產生和定義三維空間參數。FLAC3D計算步驟為建立網格，劃分網格，選擇本構模型，賦參數，選擇初始邊界條件達到原始平衡，模擬開挖監測坡內響應[1-2]。

2 研究區三維模型建立

2.1 模型的建立

正確的地質模型是計算的基礎，在計算中不可能也無必要將所有細節考慮進去，因此對工程地質條件的深入認識與抽象是建立合理地質概化模型的前提和關鍵。研究區邊坡順擬建公路走向約121 m，垂直公路走向約151 m，最大垂直開挖深度約28 m。邊坡上堆積的殘坡積物較薄，基本為一巖質邊坡，巖性主要為白云巖、泥質白云巖、泥巖等三大類，根據現場調查，研究區邊坡并無大型斷層，結構面面等，因此不考慮斷層及結構面等的影響，為模型長度取121 m，寬度取151 m，底部高程為800 m，參照壩址區工程地質平面圖、剖面圖以及邊坡巖級分帶圖等，建立研究區數值模型。FLAC3D數值模型建立方法大體上有兩種，一種是通過前處理較為方便的其他軟件建立數值模型并劃分網格，將網格節點和單元信息導出，經處理后再導入FLAC3D中[3]；另一種是直接在FLAC3D程序中采用堆砌法生成模型,本文采用第一種方法，利用cad-surfer-ansys-flac聯合建模方法，建立其數值模型。見圖1。

圖1 研究區三維模型Fig.1 Three dimensional model of research area

2.2 模型本構關系的選取

數值計算中，將巖土體看作理想彈塑性材料，屈服準則采用內置的Mohr-Coulomb準則，即：

(1)

式中：NΦ=(1+sinφ)(1-sinφ)；σ1為最大主應力；σ3為最小主應力；φ為內摩擦角；C為內聚力。

若fs<0，表示巖體發生剪性屈服；fs≥0，表示巖體未發生剪性屈服。

當法向應力為張應力時，則超出了莫爾-庫侖準則的力學有效性范圍，此時要求最小主應力不得超過巖體的抗拉強度，否則巖體將出現張性屈服，其判別公式為[4]：

ft=σ3-σt

(2)

若ft≥0，表示巖體張性屈服；若ft<0，表示巖體未發生張性屈服。

對于FLAC-3D數值計算所要求的參數體積模量(K)和剪切模量(G)可以通過下式進行換算[5]：

(3)

(4)

式中：E為彈性模量；γ為泊松比。

實際上，在FLAC3D中使用M-C本構模型時也可直接輸入彈模和泊松比，程序可實現自動換算。

2.3 模型邊界的選取

FLAC3D提供兩種應力邊界條件，一種是位移(速度)邊界條件，一種是應力邊界條件。應力邊界一般用于需考慮原巖應力的情況下。對模型底部采用完全約束，約束其3個方向的自由度；而對模型的前后左右均只作單向約束，即僅約束水平方向速度而不約束其豎向速度。除此之外，由于上部房屋的影響，因此需要考慮施加均布載荷，就需要使用應力邊界荷載條件。這樣，巖土體在生成自重應力階段可發生自由沉降，符合實際情況。

3 邊坡穩定性分析

3.1 自然斜坡穩定性分析

為了解自然斜坡在天然條件下的穩定性情況，考慮上部房屋荷載情況下斜坡在自身重力作用下的穩定性，分別進行斜坡應力分析、斜坡變形分析以及斜坡穩定性分析

3.1.1 斜坡應力分析

斜坡應力分析主要從斜坡的重力場以及最大最小主應力場3個方面進行，斜坡重力場模擬圖見圖2-圖7。

由圖2-圖3可知，坡體應力量值范圍在-2.04～0.13 MPa，坡體應力分布均勻，自上而下逐步增加，應力最大值分布在斜坡地面，坡面由于施加有房屋荷載，所以在一定范圍內，坡表的應力值較大，極少數地方出現拉應力。

圖2 坡體自重應力云圖Fig.2 The slope self-weight stress nephogram

圖3 坡表自重應力云圖Fig.3 The slope surface self-weight stress nephogram

由圖4-圖5可知，邊坡最大主應力的量值范圍為-0.85～0.58 MPa,坡體的最大主應力位于底面，與重力方向的應力相似，坡體中大主應力分布均勻，由上往下逐漸增大，坡體中無應力集中現象。在坡表最大位移云圖中可以知道，位于右側的局部地區出現了拉應力。主要是地形起伏大以及層面間的接觸面處。

圖4坡體最大主應力云圖Fig.4 The slope max principal stress nephogram

圖5 坡表最大主應力云圖Fig.5 The slope surface max principal stress nephogram

由圖6-圖7可知，斜坡的最小主應力的量值范圍為-2.04～0.09 MPa,坡體中最小主應力分布與重力分布類似，由上到下逐漸增大，均勻分布。坡表的最小主應力云圖可知，只有作用均布荷載的地方應力出現一點異常，其他都正常。

圖6 坡體最小主應力云圖Fig.6 The slope mini principal stress nephogram

圖7 坡表最小主應力云圖Fig.7 The slope surface mini principal stress nephogram

3.1.2 斜坡變形分析

坡體變形一定程度上可以反映邊坡的穩定狀態，斜坡的變形分析主要從X方向以及Y方向的位移入手，模擬結果圖見圖8-圖9。

圖8 斜坡X方向位移云圖Fig.8 X direction displacement nephogram of slope

圖9 斜坡Y方向位移云圖Fig.9 Y direction displacement nephogram of slope

圖8-圖9為數值模擬結果圖，從作圖可以知道，坡表X方向的最大位移為2.2 mm，主要集中不同巖性相交位置處，坡表Y方向最大位移為0.6 mm，主要集中在施加均布荷載的前段。自然斜坡的變形量非常小

3.1.3 斜坡穩定性分析

塑性區以及剪應變的增量可以反映邊坡的穩定狀態，一般來說，邊坡塑性區越大，深度越深，則對邊坡的不離程度也相即增大，對模擬結果進行切片，分別沿斜坡橫向以及縱向切片，得到橫向以及縱向的最大剪應變速率。見圖10、圖11。

圖10 斜坡橫向剪應變速率云圖Fig.10 Transversal shear strain rate nephogram of slope

圖11 斜坡豎向剪應變速率云圖Fig.11 Longitudinal shear strain rate nephogram of slope

由圖10、圖11可知，橫向以及豎向的剪應變速率都非常小，幾乎可以忽略不計，說明斜坡處于一個穩定的狀態。

3.2 開挖條件下邊坡穩定性研究

開挖方案為沿著開挖線垂直開挖至指定高程843.0 m處，最大的開挖深度局部可以達到28.0 m，開挖將會形成陡立邊坡，模擬開挖之后邊坡的穩定性以及變形。

首先在自然斜坡達到上述自重平衡的條件下，對原本存在的位移進行清零處理，以便可以清楚掌握在開挖情況下邊坡產生的位移情況。

3.2.1 邊坡應力分析

開挖邊坡主要從最大主應力、最小主應力進行初步分析，分析是否存在應力集中部位，模擬結果見圖12-圖16。

由圖12-圖13可以看出，開挖情況下，坡體的最大主應力總體分布較均勻，自上而下逐步增加，最大主應力的量值范圍為0.055～-0.85 MPa，但在開挖的坡腳出現了應力集中的情況，應力的大小約為-0.2 MPa。此外，在坡表作用于均布荷載的地方，應力也出現“異常”，其他地方都應力分布形式都為正常情況。

圖12 坡體最大應力分布圖Fig.12 Max stress distribution in slope

圖13 坡表最大應力分布圖Fig.13 Max stress distribution in slope surface

圖14 坡體最小應力分布圖Fig.14 Mini stress distribution in slope

圖15 坡表最小應力分布圖Fig.15 Mini stress distribution in slope surface

坡體最小主應力圖的分布形式與坡體最大主應力分布形式類似，應力分布總體均勻，自上而下應力逐漸增大，應力分布的量值為0.008～-2.054 MPa，最小主應力的最大值分布在坡地，最小值分布于坡表面，在開挖的坡腳段也出現了應力集中的現象，作用于均布載荷的上部應力“異常”，邊坡其他部位都處于正常狀況。

從圖16可以看出，邊坡的最大剪應力分布在開挖段坡腳部分，剪應力的量值范圍為0.005～0.682 MPa，圖16紅色區域就是最大剪應力分布區域，值的大小約為0.67 MPa。綜合以上分析可知，在開挖坡腳部分出現應力集中的情況，可能對邊坡的穩定性造成不小的危害，需要進行進一步的分析論證。

圖16 邊坡最大剪應力分布圖Fig.16 Max shear stress distribution of slope

3.2.2 邊坡變形分析

對邊坡的變形分析主要從開挖之后坡體在X以及Y方向的位移及其位移值，布設的監測點進行分析。見圖18-圖19。

從邊坡X方向的位移云圖可知，X方向最大位移量值方位為-13.6～12.5 mm，最大位移主要是層面接觸方向以及均布荷載作用的邊緣，如圖左上角紅色區域為均布荷載作用邊緣處。在Y方向的位移云圖中可知，位移的量值范圍為3.2～-93.2 mm，最大位移集在相對巖體力學性質軟弱處，在邊坡的右側，由于其開挖量比較小，邊坡Y方向的位移值比較小，呈現mm級別的，但是在邊坡左側由于上部房屋的作用以及開挖高度大的影響，邊坡Y方向的位移呈cm級別，最大位移為9.3 cm，邊坡處于不穩定狀態。

3.2.3 邊坡穩定性分析

塑性區以及剪應變的增量可以反映邊坡的穩定狀態，一般來說，邊坡塑性區越大，深度越深，則對邊坡的不離程度也相繼增大，對模擬結果進行切片，分別沿斜坡橫向以及縱向切片，得到橫向以及縱向的最大剪應變速率。見圖19-圖22。

圖17 邊坡X方向位移云圖Fig.17 X direction displacement nephogram of slope

圖18 邊坡Y方向位移云圖Fig.18 Y direction displacement nephogram of slope

由圖19-圖22可以看到，X方向的最大剪切應變處于開挖坡腳線附近，其量值范圍為1.00×10-6～8.48×10-6，Y方向的最大剪切應變位于坡表，量值范圍為1.00×10-7～9.50×10-7。從Y方向塑性區的分布云圖可以知道，開挖線后約20 m范圍內巖體處于受剪切的狀態，受剪區與開挖線的夾角大約為45°。在X方向塑性區分布云圖可知，在高程為850.0 m以上的巖體基本處于受剪切狀態，表面巖體受拉張作用。

圖19 X方向最大剪應變云圖Fig.19 X direction max shear strain nephogram

圖20 Y方向最大剪應變云圖Fig.20 Y direction max shear strain nephogram

圖21 X方向塑性區分布云圖Fig.21 X direction plastic zone distribution nephogram

圖22 Y方向塑性區分布云圖Fig.22 Y direction plastic zone distribution nephogram

4 結 論

由FLAC3D數值模擬軟件做出的以上兩種分別在自然工況和開挖工況的穩定性分析，得到以下結論：

1) 自然條件下斜坡處于一個穩定的狀態。

2) 在該方案開挖條件下，邊坡的位移量值非常大，邊坡開挖，房屋的穩定性會受到極大的影響，容易造成房屋的開裂，建筑物發生拉裂縫，開挖線后20 m以及高程850.0 m以上的巖體都會因開挖受到非常大的影響，在此種方案下，邊坡處于不穩定狀態，極易產生失穩變形破壞。