ANSYS有限元數值模擬在某開挖高邊坡穩定性分析中的應用

許容 王輝

摘要：ANSYS有限元數值可以很好地模擬巖土體的力學性能，甚至通過選取適當的輸入參數和計算模型，也能模擬出節理裂隙、軟弱夾層、活動性斷裂等地質情形。文章結合某開挖高邊坡工程實例，利用ANSYS有限元軟件，通過應力、變形以及能量計算分析邊坡的穩定性。

關鍵詞：ANSYS有限元;數值模擬;穩定性

0 引言

邊坡穩定分析方法很多，常用的方法主要是通過應力、變形以及能量計算分析邊坡穩定性。其采用矩陣分析模式，更加方便編程分析。應當注意的是巖土工程在做有限元分析時，不應當過度注重其具體量值的大小，應更加注重應力應變的分布情況與相對的變化情況。

1 邊坡工程地質概況

該邊坡位于云南省境內，構造主要受北西向和近南北向構造控制，區域地質構造背景十分復雜。該邊坡區2#傾倒蠕變巖體在邊坡沖溝下游側1 570～1 840 m高程處發育。變形體表層覆蓋有塊石、碎石混粉土的坡積層;傾倒巖體主要為板巖（J2h）;下伏基巖為變質火山角礫巖夾片理化變質凝灰巖（T3xd）。

2 模型建立及數值模擬

2.1 有限元模型建立

選取開挖邊坡橫剖面，該剖面位于正面開挖邊坡中段，剖面線總長535 m，前緣高程1 540 m，后緣高程1 960 m。右岸邊坡開挖至纜機平臺，高程為1 678 m（見下頁圖1）。

本文通過將剖面的二維極限平衡的滲流場通過插值的方式賦予網格中的各個節點，具體步驟如下：

首先將對各個工況下的滲流場進行模擬，得出某個時刻的孔隙水應力圖（見下頁圖2）。提取圖中各個坐標和孔隙水壓力值，通過插值的方法將坐標輸入ANSYS FLUENT模塊中，使用TB、PM命令制定的滲透系數和BIOT系數定義孔隙介質屬性。然后，通過fluent輸出.cdb格式文件，最后再導入ANSYS中。

ANSYS有由頂向下和由底向上兩種建模方式。本論文主要使用由底向上方式建立邊坡的主體模型：首先提取壩址區CAD地質圖中的等高線數據，存為XYZ坐標文件，然后導入到Surfer，選取建模區域內的坐標，通過克里格插值法將這些不規則高程點離散為均勻的網格點，導出這些網格數據為文本，最后在Excel中處理為ANSYS可以識別的APDL格式文件。再根據此初步模型，將上述APDL文本輸入ANSYS，生成地表關鍵點，再通過“由底向上”的方式建立包括下切地層的完整邊坡模型，然后依次分割出各個下切地層、分級開挖體和不同卸荷層等模型分區，最后對幾何模型劃分網格。采用智能剖分（SmartSize）和限定單元尺寸的方法優化網格，多次嘗試后劃分出合適的網格模型。

一般情況下，ANSYS模型越大、越精細，其求解精度就會越高，但是求解的時間也越多。根據本論文的需要和求解時間的限制，本論文建立模型的大小為沿河流方向（X軸）為364 m，垂直于河流方向（Y軸）為363 m，由于邊坡在1 678 m高程，山頂高程為1 906 m（見圖3），考慮到計算的合理性，建模的最低高程應在邊坡高程以下，以防止在模型底部產生的應力集中影響計算結果。模型的頂部應至坡體頂部，能更合理地計算巖土體自重產生的應力應變。綜上所述，建模高程（Z軸）取1 500～1 886 m中間的386 m高程。考慮模型開挖的過程，在模型開挖面以及開挖高程的同時建模。

模型使用SOLID95單元四面體剖分巖土體部分，SOLID95單元具有塑性、蠕變、單元生死、初始應力輸入等適合巖土計算的特性。設置法向接觸剛度為1.0，防止過大的接觸面嵌入并且保證數值計算的穩定性;使用拉格朗日乘子法計算，設置侵入容差為0.05，以定義允許的最大侵入值;設置力學表面作用模式為法向單相接觸，使得模型在計算時，如果接觸處于分離狀態，法向壓力為零;單元網格劃分時，依然綜合考慮了計算精度與計算效率的問題，在坡體表面、開挖區域以及斷層區域使用較小網格剖分，在模型的底部，由于形狀規則，使用較大的網格剖分，模型單元數為10 103，節點數為37 460（見圖4）。

模型分析對于力的容差為0.005，對于位移的容差為0.05，設置使用L2范數檢查SRSS值，位移允許的最小值為1.0。

2.2 數值計算參數

綜合前期勘察報告和作者現場跟蹤調查成果，確定了不同巖層的數值計算參數取值，具體見表1。

3 有限元模擬計算分析

3.1 天然邊坡應力應變分析

邊坡的應力重分布即二次應力是天然狀態邊坡應力場特征的研究關鍵，對于高地應力下的變形破壞具有很高的研究意義。原因在于天然應力場對邊坡開挖后應力重分布的趨勢和結果起決定作用;同時在天然應力場的影響下形成的邊坡內部應力增高帶和應力集中也是影響邊坡穩定性的因素。

3.1.1 應力場分析

天然狀態下邊坡大主應力場主要來源于巖體自重力，隨著巖體深度的變化，分布較均勻，主要為壓應力，局部位置為拉應力。邊坡的坡面為臨空面，初始的最大主應力為沿坡面切線方向受重力影響指向坡面以下，表現為拉應力集中在邊坡頂部附近下游部位，此處為J2h千枚狀泥質板巖，且后緣部位為2#傾倒變形巖體，最大主應力與最小主應力都表現為較小應力，與現實相符合。2#傾倒變形巖體附近最大主應力大小為0.82～1.33 MPa，坡底部表現出由重力引起的較大的壓應力，大小為1.25～3.23 MPa;坡面最小主應力為拉應力，大小為0～0.27 MPa，在邊坡底部有較大的壓應力。

3.1.2 位移場分析

天然狀態下，邊坡模型在施加自重應力后會產生初始位移（見圖5～7）。后面的工況中，ANSYS計算出的位移包含了初始位移，因此在繪制位移云圖時從現位移中減去該初始位移，即可得到因荷載改變產生的增量位移，即真實位移。

3.2 開挖邊坡應力應變分析

3.2.1 應力場分析

天然邊坡坡面主應力較小，在開挖后J2h板巖傾倒變形體開挖面部位在水平方向上位移較大，開挖后的臨空面應力釋放，最大主應力由壓應力變為拉應力，大小為0～0.29 MPa，最小主應力為0～0.21 MPa，邊坡下游圓弧開挖面應力集中現象不明顯，上游折線開挖面在坡腳與開挖折線處產生較大的應力集中，最大主應力大小為3.30 MPa，此開挖腳處產生較大的塑性變形（見圖8和圖9）。坡體內部最大主應力在F12斷層附近出現應力松弛效應，在F12斷層附近，斷層上盤的主應力大于斷層下盤的主應力。邊坡水平向位移最大位置在開挖面處J2h板巖傾倒變形體F12斷層上盤附近，最大主應力在斷層上盤的應力值大于下盤，連續的應力等值線在此處錯開。在J2h千枚狀泥質板巖與T3xd火山角礫巖、碎屑巖夾凝灰巖巖層分界面處，T3xd巖層最大主應力應力值大于J2h巖層應力。

3.2.2 位移場分析

水平向的位移仍然以上部J2h板巖區域傾倒變形體最大。在坡體內部斷層與巖層界面之間的巖體主應力值相較兩側巖體小，斷層具有應力松弛效應，傾倒變形體與下伏基巖界面上部變形體應力值小于下部基巖。2#傾倒變形體的后部J2h千枚狀泥質板巖在表層有較大的指向坡外的位移，基本與現場監測情況一致。開挖期間，邊坡向臨空方向卸荷回彈變形以及沿各緩傾斷層滑移變形繼續增大（見圖10）。

4 結語

本文對云南省某開挖邊坡進行了三維模擬，計算邊坡天然應力，對邊坡進行開挖模擬，計算開挖對邊坡巖體的應力應變影響，并對比開挖后的應力應變、位移，得出開挖邊坡的穩定性分析結論：

（1）針對邊坡開挖的卸荷變形問題，介紹了使用有限單元法模擬邊坡開挖的方法。

（2）自然邊坡在人工開挖后，邊坡開挖面在卸荷作用下，開挖面附近坡體原有應力場將發生變化，應力場的變化導致了在開挖面附近巖土體產生變形。變形在開挖面上最大，向坡體內部逐漸減小。在開挖過程中，開挖面附近巖體卸荷，由于坡面沒有約束，應力釋放明顯。隨著向坡體內部的深入，隨著巖土體間的內摩擦力與粘聚力的作用越來越強，應力釋放就越來越小，但是這種變化是有限制的，開挖邊坡存在一個卸荷影響范圍，當超出此范圍以后開挖卸荷影響就會減小，此處邊坡巖土體的應力就基本仍為原邊坡的應力狀態。

（3）邊坡的人工開挖對開挖范圍附近的應力分布會有較大影響，其程度會因開挖面的深度、開挖面的形狀與開挖方式而不同。開挖后開挖面會發生卸荷，隨著初始應力由坡頂向坡底逐漸增加，開挖的卸荷量也逐漸增加。在開挖的臺階處發生了較大的應力集中，會導致在開挖的過程中就有巖體的卸荷變形現象。順坡向開挖坡腳出現向外的位移，同時豎直方向在開挖面上出現指向坡外的位移，上游側折線開挖面在坡腳處常常會發生應力集中，具有較大的張拉應力，這些都是由于巖體傾倒蠕變所導致的，如果不及時加固會發生巖體卸荷變形。所以開挖時應避免較大折角的產生。

（4）現場監測分析表明該開挖邊坡穩定，理論分析與實際結論一致，可為今后類似工程模型的建立提供依據，具有一定的理論意義及工程運用價值。

參考文獻：

[1]王敏強，許 原. 有限元分析中開挖釋放荷載計算的討論[J]. 武漢大學學報（工學版），2001，34（1）：56-59.

[2]章 青.有限元分析中開挖釋放荷載的正確計算[J].河海大學學報，1999，27（3）：112-115.

[3]李 圍. ANSYS土木工程應用實例[M]. 北京：中國水利水電出版社，2006.

[4]王 勇，殷宗澤.有限元計算深開挖中挖方等效荷載的分析[J].河海大學學報，1998，26（5）：71-74.

收稿日期：2020-06-02