基于無網格方法的邊坡支擋結構變形特征研究

徐曉東

（內蒙古包鋼鋼聯股份有限公司巴潤礦業分公司）

隨著露天礦開采規模的擴大及開采程度的深入，采礦邊坡滑坡儼然已成為一種主要的地質災害［1］，對礦產開采、采礦設備和人員安全造成嚴重影響。目前的研究對邊坡穩定性和失穩原因的探究較為成熟，但對于邊坡失穩后所造成的沖擊影響研究較少［2］。因此，明確邊坡失穩后的滑動位移和沖擊效應是采礦工程領域重要的研究課題。

邊坡失穩時通常會出現大變形特征，當發生流動大變形時，變形不再符合小變形理論，而是進入了非線性的大變形狀態［3-4］。此時，傳統的數值模擬方法，如有限元法（FEM）和有限差分法（FDM），由于其網格屬性，容易在計算大變形時產生網格扭曲，導致對邊坡失穩后的滑移大變形分析存在一定困難［5］。為了克服目前常用的數值模擬方法中的不足，離散元法（DEM）、非連續變形分析法（DDA）、無網格伽遼金法（EFG）以及細胞自動機法（CA）等離散介質力學方法和無網格方法得到了廣泛應用。然而，這些方法在研究邊坡變形方面還存在許多缺點，如難以確定非物理參數和直接描述土體的應力-應變關系等［6-10］。

針對滑體沖擊影響的問題，國內外也進行了一定的研究。在理論方面，通常采用速滑計算方法計算滑體的沖擊能，如能量法［11］和條分法［12］等。根據相似性原理，利用物理模型實驗模擬滑體的沖擊過程，評估沖擊的影響也有涉及［13-15］。然而，巖土體是1 種散體，理論方法在經過過多的簡化和假設后，不能準確描述滑體滑動過程及沖擊破壞過程［16］。基于物理模型的相似性設計相對復雜麻煩，其制作周期長、費用高［17］。

基于此，本文采用1種無網格方法—光滑粒子流體動力學法（SPH）對邊坡變形特征及沖擊效應進行數值模擬研究。結合HBP（Herschel-Bulkley-Papanastasiou）本構模型，模擬了經典的Poiseuille 流問題，同時對比分析了Bingham 模型的流體運動問題。進一步建立露天礦邊坡數值模型，對邊坡失穩演化過程及滑體對支護結構的受力變形規律以及支護效果進行模擬分析，分析柔性支護結構在邊坡失穩過程中的沖擊作用及受力特征，為后期同類工程的施工提供一定的借鑒。

1 SPH基本理論

1.1 基本控制方程

SPH方法是一種具有拉格朗日特征的粒子法，該方法的基本思想是對研究材料進行離散化處理，將其離散為能夠代表研究材料本身的粒子（圖1）。每個粒子攜帶自身的材料屬性，如質量、密度、速度和能量等，并通過跟蹤這些粒子特征反映研究材料的狀態，同時遵循物質的運動守恒控制方程。

圖1 中h為光滑長度；k值為影響系數，表示影響域大小，一般根據所研究問題來確定k值，隨著k值的增大，計算精度會提高，但同時也會降低計算效率；W為核函數；i為計算粒子；j表示相鄰粒子；rij表示粒子i到j的距離。

SPH計算方法中，研究具有材料強度的控制方程有粒子密度、速度和位置變化方程，如下：

式中，ρi為第i個粒子密度，g/m3；xα為空間坐標；vα為速度矢量；σαβ總應力張量，Pa，其拉應力為正，壓應力為負，一般由各向同性壓力P和剪應力τ構成，計算公式為σαβ= -pδαβ+ταβ；α，β為坐標方向。

經過離散化處理，式（1）和式（2）在SPH 方法下的計算公式：

式中，mj為第j個粒子的質量，g；P為正應力，Pa；g為重力，N；Wij為核函數，通常情況下，3次B樣條函數在計算精度和計算效率均比較好。

1.2 HBP本構模型

邊坡從失穩到流滑的動力過程十分復雜，一般認為其經歷了失穩—滑動—流滑3 個過程。為準確地描述邊坡在大應變狀態下剪應力-剪應變非線性關系，本文將變形體視為具有可變黏度的黏性材料，采用HBP 模型表征變形體的流變特性。HBP 模型由Papanastasiou［18］提出，剪切力是通過與應變張量相關的剪應力張量表示，其本構方程如下：

式中，τij為剪應力，Pa；μeff為有效黏滯系數，描述流體的流變特性；εαβ為應變張量，Pa。

其中，有效黏滯系數計算公式：

式中，γ為剪切應變率；μ為表觀動態黏度系數；τy表示屈服應力，Pa；m，n為計算常數，可以模擬剪切變薄或剪切增稠行為，當m=0 和n=1 時，模型變為牛頓模型，而m→∞且n=1時，模型簡化為Bingham模型。

因此，最終SPH動量方程可表示如下：

1.3 時間步長確定

在SPH 求解過程中，需要對控制方程和本構方程進行離散求解，求解方法有蛙跳算法、預測矯正法以及龍格-庫塔法等，在求解過程中一個關鍵參數是時間步長的確定。一般情況下，蛙跳算法所需存儲量低，在每一次計算中只需要進行一次優化估值等優點，其迭代過程如下：

式中，ci代表聲速，m/s；vi代表粒子速度，m/s；CFL是穩定性常數，小于1.0。

2 HPB本構模型數值驗證

2.1 Poiseuille流動

選取典型的Poiseuille 流為模擬算例，流體模型為在2個固定不動的平板之間充滿水，并且初始狀態下水處于靜止狀態。接著，當受到水平方向的外力作用時，水開始流動。采用SPH 方法模擬瞬態流動過程時，建立如圖2 所示的Poiseuille 流模型，計算區域為2 m×1 m 的矩形。在此計算中，主要的計算參數：粒子個數1 640 個，粒子的初始間距r0=0.025 m，光滑長度h=1.2r0，水平外力F=1 N，初始參考密度ρ0=1 000 g/m3，時間步長為1.0 s。

該流體模型在初始狀態下（t=0）位于2 塊固定的無限大平板之間且處于靜止狀態。然后，在t＞0 時，一個恒定的水平作用力作用于整個流場，導致最初靜止的流體開始逐漸流動。流動速度如圖3所示。

從圖3中可以看出，流體的速度沿著水平軸中心線具有對稱性。即靠近平板附近的流體由于黏附在靜止的平板上，因此它們的速度較小且近似為零；而流道中心的速度最大，這與日常物理現象相一致。在經歷0.8 s 后，流動達到最大狀態，其最大流速約為0.006 m/s，如圖3（c）所示。為了進一步驗證方法的正確性，在流場中選取了20個空間測量點（圖2）采用有限容積法（FVM 方法）與所提方法進行對比，得到HBP 本構模型條件下Poiseuille 流的速度對比圖，如圖4所示。

由圖4 可看出，在時間步長分別取1.0，2.0 和5.0 s 時，平行板兩側粒子速度幾乎相同，靠近兩側粒子速度逐漸減小，各時刻流速分布圖趨勢均是拋物線型，并且左右對稱。SPH 方法與FVM 方法的所得的結果具有較好的吻合性，驗證了SPH 方法在模擬黏彈性流體瞬態流動問題的準確性和有效性。

2.2 Bingham模型對比分析

基于HBP 和Bingham 2 種本構模型，建立三維水槽流動模型，模型兩側為Bingham 本構模型，中間為HBP 本構模型。粒子初始間距0.008 5 m，光滑長度為粒子初始間距的1.2 倍，初始時間步長為1.0×10-4s，對比不同本構關系下的SPH數值模擬結果，得到了各瞬態下粒子的分布情況以及速度場分布情況，如圖5所示。如圖5（a）所示，水流運動過程中最大傳播速度出現于HBP 本構模型的流體前端部位，并逐漸向后端遞減，且相對于Bingham 流體模型更為明顯。隨著時間的推移，流體前端速度繼續加大，并且運動距離明顯大于上下兩側流體，如圖5（b）所示。如圖5（c）所示，t=60.0 s 時，流體速度達到峰值后逐漸減小至零，達到靜止平衡狀態，從橫向速度場分布及運動距離上來看，HBP 本構模型可以更好地模擬變形流變特征，最接近水體自由運動。

3 支擋變形數值模擬

3.1 一維流體與支擋

考慮自由流體與支護結構相互作用模型，其上部粒子數為12 160 個，表征自由流體，受自重作用，高度4 m，寬度2 m，如圖6（a）所示。下部粒子數352模擬固體支擋結構，對其進行封閉。如圖6（b）所示，底部支擋結構收到壓力主要集中在兩側，支擋結構未產生變形。隨著時間增加，如圖6（c）所示，底部支擋結構應力集中轉移至支擋結構中心位置，且在壓應力的作用下，支擋結構發生彎曲，但是未達到損傷、斷裂現象。為了更清楚地分析流體與支護結構之間的相互作用，選取驗證模型底部中點為監測點，其位移、壓力與時間關系如圖7所示。

如圖7（a）所示，驗證模型測點A 位置位移隨時間增加而增加，當計算時間達到0.4 s 時，監測點位移大約為0.16 m。如圖7（b）所示，應力也具有同樣的規律，即隨著時間的增加，支擋結構所承受的壓應力也隨之增大。通過分析該模型，可知上述方法在變形和支擋結構方面的可行性和正確性。

3.2 柔性支擋

在模型初始高度54 m，初始長度180 m，支護結構長度14 m 的條件下，邊坡變形失穩運動演化過程及支護沖擊相互作用如圖8 所示。數值模擬得到了各瞬態下粒子的分布情況以及沖擊效應分布情況，如圖9所示。

由圖8（a）可以看出，當邊坡失穩后，滑體與支護結構接觸時，支護結構受到的最大沖擊力主要作用在土體與接觸部位，即支護結構的底部，但是由于滑動速度較小，沖擊不足以使支護結構發生變形。隨著邊坡滑動時間的增加，邊坡加速滑動向前，如圖8（b）所示，支護結構阻擋了絕大數失穩土體，且受到的沖擊力逐漸上移，支護結構發生彎曲傾倒變形。當邊坡失穩速度達到最大時（圖9（a）），此時滑體體積也是最大，部分失穩土體越過支護結構，此時支護結構變形達到最大值，如圖8（c）所示。隨著邊坡失穩速度逐漸減小直至邊坡停止，如圖9（b）所示，此時達到新的平衡狀態，支護結構所受到的沖擊力突然驟降，柔性支護結構恢復自身變形，但是其底部出現較小的殘余應力作用（圖8（d）），主要是由于失穩土體的土壓力造成的。

為了能夠進一步證明該方法的可行性，基于上述模型，對支護結構不同錨固深度時，邊坡失穩沖擊效應進行了研究，得到邊坡在滑動過程中，支護結構在不同錨固深度的壓力和位移變化數據，如圖10 所示。

通過圖10 可以看出，隨著支護結構錨固深度的增加，邊坡在滑動過程中支護結構受到沖擊作用力呈現增長的趨勢，并且最終由于邊坡停止失穩，沖擊力均趨于非零穩定值；當錨固深度3.0 m時，支護結構豎向位移變化平穩，表明邊坡滑動過程中，支護結構錨固深度越大，支護結構在受到沖擊作用時結構越穩定。

4 結 論

本文基于光滑流體粒子動力方法，結合HBP 本構模型，構建了模擬邊坡大變形與支擋結構相互作用過程的數值模型，得出以下結論：

（1）選取HBP 模型作為邊坡變形本構模型，有效地克服了土體變形在塑性應變過渡段及大變形情況下可能呈現出非線性流變特性的問題，其對于分析變形較大和兩相物質相互作用的問題具有一定的可行性。

（2）由于SPH 方法不依賴于網格及其純拉格朗日特性，自然而然在處理具有變形或者不連續問題上有比較好的優勢，為研究該類問題提供了一種新的方法和研究思路。

（3）露天礦邊坡在地震荷載作用下的穩定性也是工程設計關注的重點問題，后續工作將基于SPH深入開展邊坡動力失穩機理和滑動過程研究。