徐曉東

（內(nèi)蒙古包鋼鋼聯(lián)股份有限公司巴潤(rùn)礦業(yè)分公司）

隨著露天礦開(kāi)采規(guī)模的擴(kuò)大及開(kāi)采程度的深入，采礦邊坡滑坡儼然已成為一種主要的地質(zhì)災(zāi)害［1］，對(duì)礦產(chǎn)開(kāi)采、采礦設(shè)備和人員安全造成嚴(yán)重影響。目前的研究對(duì)邊坡穩(wěn)定性和失穩(wěn)原因的探究較為成熟，但對(duì)于邊坡失穩(wěn)后所造成的沖擊影響研究較少［2］。因此，明確邊坡失穩(wěn)后的滑動(dòng)位移和沖擊效應(yīng)是采礦工程領(lǐng)域重要的研究課題。

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

針對(duì)滑體沖擊影響的問(wèn)題，國(guó)內(nèi)外也進(jìn)行了一定的研究。在理論方面，通常采用速滑計(jì)算方法計(jì)算滑體的沖擊能，如能量法［11］和條分法［12］等。根據(jù)相似性原理，利用物理模型實(shí)驗(yàn)?zāi)M滑體的沖擊過(guò)程，評(píng)估沖擊的影響也有涉及［13-15］。然而，巖土體是1 種散體，理論方法在經(jīng)過(guò)過(guò)多的簡(jiǎn)化和假設(shè)后，不能準(zhǔn)確描述滑體滑動(dòng)過(guò)程及沖擊破壞過(guò)程［16］?；谖锢砟Ｐ偷南嗨菩栽O(shè)計(jì)相對(duì)復(fù)雜麻煩，其制作周期長(zhǎng)、費(fèi)用高［17］。

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

1 SPH基本理論

1.1 基本控制方程

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

圖1 中h為光滑長(zhǎng)度；k值為影響系數(shù)，表示影響域大小，一般根據(jù)所研究問(wèn)題來(lái)確定k值，隨著k值的增大，計(jì)算精度會(huì)提高，但同時(shí)也會(huì)降低計(jì)算效率；W為核函數(shù)；i為計(jì)算粒子；j表示相鄰粒子；rij表示粒子i到j(luò)的距離。

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

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

經(jīng)過(guò)離散化處理，式（1）和式（2）在SPH 方法下的計(jì)算公式：

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

1.2 HBP本構(gòu)模型

邊坡從失穩(wěn)到流滑的動(dòng)力過(guò)程十分復(fù)雜，一般認(rèn)為其經(jīng)歷了失穩(wěn)—滑動(dòng)—流滑3 個(gè)過(guò)程。為準(zhǔn)確地描述邊坡在大應(yīng)變狀態(tài)下剪應(yīng)力-剪應(yīng)變非線性關(guān)系，本文將變形體視為具有可變黏度的黏性材料，采用HBP 模型表征變形體的流變特性。HBP 模型由Papanastasiou［18］提出，剪切力是通過(guò)與應(yīng)變張量相關(guān)的剪應(yīng)力張量表示，其本構(gòu)方程如下：

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

其中，有效黏滯系數(shù)計(jì)算公式：

式中，γ為剪切應(yīng)變率；μ為表觀動(dòng)態(tài)黏度系數(shù)；τy表示屈服應(yīng)力，Pa；m，n為計(jì)算常數(shù)，可以模擬剪切變薄或剪切增稠行為，當(dāng)m=0 和n=1 時(shí)，模型變?yōu)榕ｎD模型，而m→∞且n=1時(shí)，模型簡(jiǎn)化為Bingham模型。

因此，最終SPH動(dòng)量方程可表示如下：

1.3 時(shí)間步長(zhǎng)確定

在SPH 求解過(guò)程中，需要對(duì)控制方程和本構(gòu)方程進(jìn)行離散求解，求解方法有蛙跳算法、預(yù)測(cè)矯正法以及龍格-庫(kù)塔法等，在求解過(guò)程中一個(gè)關(guān)鍵參數(shù)是時(shí)間步長(zhǎng)的確定。一般情況下，蛙跳算法所需存儲(chǔ)量低，在每一次計(jì)算中只需要進(jìn)行一次優(yōu)化估值等優(yōu)點(diǎn)，其迭代過(guò)程如下：

式中，ci代表聲速，m/s；vi代表粒子速度，m/s；CFL是穩(wěn)定性常數(shù)，小于1.0。

2 HPB本構(gòu)模型數(shù)值驗(yàn)證

2.1 Poiseuille流動(dòng)

選取典型的Poiseuille 流為模擬算例，流體模型為在2個(gè)固定不動(dòng)的平板之間充滿水，并且初始狀態(tài)下水處于靜止?fàn)顟B(tài)。接著，當(dāng)受到水平方向的外力作用時(shí)，水開(kāi)始流動(dòng)。采用SPH 方法模擬瞬態(tài)流動(dòng)過(guò)程時(shí)，建立如圖2 所示的Poiseuille 流模型，計(jì)算區(qū)域?yàn)? m×1 m 的矩形。在此計(jì)算中，主要的計(jì)算參數(shù)：粒子個(gè)數(shù)1 640 個(gè)，粒子的初始間距r0=0.025 m，光滑長(zhǎng)度h=1.2r0，水平外力F=1 N，初始參考密度ρ0=1 000 g/m3，時(shí)間步長(zhǎng)為1.0 s。

該流體模型在初始狀態(tài)下（t=0）位于2 塊固定的無(wú)限大平板之間且處于靜止?fàn)顟B(tài)。然后，在t＞0 時(shí)，一個(gè)恒定的水平作用力作用于整個(gè)流場(chǎng)，導(dǎo)致最初靜止的流體開(kāi)始逐漸流動(dòng)。流動(dòng)速度如圖3所示。

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

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

2.2 Bingham模型對(duì)比分析

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

3 支擋變形數(shù)值模擬

3.1 一維流體與支擋

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

如圖7（a）所示，驗(yàn)證模型測(cè)點(diǎn)A 位置位移隨時(shí)間增加而增加，當(dāng)計(jì)算時(shí)間達(dá)到0.4 s 時(shí)，監(jiān)測(cè)點(diǎn)位移大約為0.16 m。如圖7（b）所示，應(yīng)力也具有同樣的規(guī)律，即隨著時(shí)間的增加，支擋結(jié)構(gòu)所承受的壓應(yīng)力也隨之增大。通過(guò)分析該模型，可知上述方法在變形和支擋結(jié)構(gòu)方面的可行性和正確性。

3.2 柔性支擋

在模型初始高度54 m，初始長(zhǎng)度180 m，支護(hù)結(jié)構(gòu)長(zhǎng)度14 m 的條件下，邊坡變形失穩(wěn)運(yùn)動(dòng)演化過(guò)程及支護(hù)沖擊相互作用如圖8 所示。數(shù)值模擬得到了各瞬態(tài)下粒子的分布情況以及沖擊效應(yīng)分布情況，如圖9所示。

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

為了能夠進(jìn)一步證明該方法的可行性，基于上述模型，對(duì)支護(hù)結(jié)構(gòu)不同錨固深度時(shí)，邊坡失穩(wěn)沖擊效應(yīng)進(jìn)行了研究，得到邊坡在滑動(dòng)過(guò)程中，支護(hù)結(jié)構(gòu)在不同錨固深度的壓力和位移變化數(shù)據(jù)，如圖10 所示。

通過(guò)圖10 可以看出，隨著支護(hù)結(jié)構(gòu)錨固深度的增加，邊坡在滑動(dòng)過(guò)程中支護(hù)結(jié)構(gòu)受到?jīng)_擊作用力呈現(xiàn)增長(zhǎng)的趨勢(shì)，并且最終由于邊坡停止失穩(wěn)，沖擊力均趨于非零穩(wěn)定值；當(dāng)錨固深度3.0 m時(shí)，支護(hù)結(jié)構(gòu)豎向位移變化平穩(wěn)，表明邊坡滑動(dòng)過(guò)程中，支護(hù)結(jié)構(gòu)錨固深度越大，支護(hù)結(jié)構(gòu)在受到?jīng)_擊作用時(shí)結(jié)構(gòu)越穩(wěn)定。

4 結(jié) 論

本文基于光滑流體粒子動(dòng)力方法，結(jié)合HBP 本構(gòu)模型，構(gòu)建了模擬邊坡大變形與支擋結(jié)構(gòu)相互作用過(guò)程的數(shù)值模型，得出以下結(jié)論：

（1）選取HBP 模型作為邊坡變形本構(gòu)模型，有效地克服了土體變形在塑性應(yīng)變過(guò)渡段及大變形情況下可能呈現(xiàn)出非線性流變特性的問(wèn)題，其對(duì)于分析變形較大和兩相物質(zhì)相互作用的問(wèn)題具有一定的可行性。

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

（3）露天礦邊坡在地震荷載作用下的穩(wěn)定性也是工程設(shè)計(jì)關(guān)注的重點(diǎn)問(wèn)題，后續(xù)工作將基于SPH深入開(kāi)展邊坡動(dòng)力失穩(wěn)機(jī)理和滑動(dòng)過(guò)程研究。