樊玉超

(國家能源集團新疆能源有限責(zé)任公司，新疆維吾爾自治區(qū)烏魯木齊市，830000)

0 引言

確定邊坡臨界滑面的位置是露天煤礦邊坡穩(wěn)定性分析的關(guān)鍵問題，目前工程中常用的方法主要分為極限平衡法和有限元法兩類。極限平衡法一般采用滑面形狀為弧形或者多段型的假設(shè)，對整個露天煤礦邊坡區(qū)域進行遍歷搜索，計算各個假設(shè)滑面上塊體的靜力平衡方程，遴選出安全系數(shù)最小的滑面即為該邊坡的臨界滑面[1-3]。有限元法采用彈塑性本構(gòu)模型計算露天煤礦邊坡的應(yīng)力場和位移場，通過強度折減、過載等方式使邊坡達(dá)到臨界狀態(tài)，對有限元計算結(jié)果進行后處理，估算出滑面的位置。

目前，用于判斷露天煤礦邊坡失穩(wěn)破壞的判據(jù)主要有3種：數(shù)值解非收斂判據(jù)[4]，以靜力平衡計算不收斂作為邊坡整體失穩(wěn)的標(biāo)志；塑性區(qū)貫通判據(jù)[5]，以塑性區(qū)或某一幅值的等效塑性應(yīng)變由坡腳至坡頂貫通作為邊坡整體失穩(wěn)的標(biāo)志；位移突變判據(jù)[6]，以土體內(nèi)某些點的應(yīng)變和位移發(fā)生突變且無限發(fā)展作為邊坡整體失穩(wěn)的標(biāo)志。

確定臨界滑動面空間位置的方法包括：基于畸變網(wǎng)格或貫通的塑性應(yīng)變等色圖[4，7]；沿深度方向搜索最大塑性應(yīng)變、最大剪應(yīng)變增量、最大位移變化率位置[8-10]；求解滑動面控制方程[11-12]；對邊坡位移及節(jié)點安全系數(shù)聚類[13]。

與極限平衡法相比，有限元法考慮巖土體變形協(xié)調(diào)條件與本構(gòu)關(guān)系，可以模擬露天煤礦邊坡漸進破壞過程，能夠更好地適用于各種復(fù)雜工況條件。然而，有限元法中邊坡失穩(wěn)狀態(tài)的判識與臨界滑面的確定通常互相割裂，彼此之間缺乏內(nèi)在關(guān)聯(lián)。為此，提出一種局部追蹤算法[14]，通過追蹤局部化帶形成、擴展、匯合直至貫通的過程，在判斷邊坡的極限平衡狀態(tài)的同時，確定臨界滑動面的位置。局部追蹤算法在處理單條局部化帶擴展時十分簡單，但在處理多條局部化帶擴展時非常繁雜，缺少魯棒性，并且構(gòu)造的局部化帶路徑也不光滑。

為解決上述問題，引入OLIVER J等[15-16]提出的全局追蹤算法(Global Tracking Algorithm)，結(jié)合局部化帶貫通判據(jù)，發(fā)展一種基于全局追蹤算法的露天煤礦邊坡臨界滑面確定方法，并通過算例分析驗證了臨界滑面確定方法的有效性。

1 材料失效理論與邊坡失穩(wěn)判據(jù)

1.1 材料失效理論

局部化帶理論[17-18]認(rèn)為，土體材料失效表現(xiàn)為由原先光滑、連續(xù)化的變形模式轉(zhuǎn)變?yōu)橐詣×易冃为M窄條帶形式出現(xiàn)的、高度局部化的變形模式。如圖1所示，當(dāng)土體材料出現(xiàn)劇烈變形的局部化帶時，位移率在跨越局部化帶時保持連續(xù)，而位移梯度率出現(xiàn)間斷，即：

圖1 土體材料應(yīng)變局部化現(xiàn)象

(1)

“+”“-”——帶狀區(qū)域的正表面和負(fù)表面；

n——局部化帶的單位法向量；

m——局部化帶的單位切向量；

1.2 邊坡失穩(wěn)判據(jù)

在局部化帶理論框架下，邊坡失穩(wěn)破壞過程可以看作土體單元不斷發(fā)生局部化失效的過程。坡體內(nèi)部局部化帶演化如圖2所示。

圖2 坡體內(nèi)部局部化帶演化

由于邊坡內(nèi)材料的非均勻性、初始弱面的存在以及非均勻分布的荷載與邊界約束，邊坡內(nèi)部的應(yīng)力場始終是不均勻的，一些應(yīng)力集中單元率先發(fā)生實現(xiàn)，形成局部化帶。隨著外部荷載的增加、材料強度的減小以及土體單元實現(xiàn)后釋放載荷在邊坡內(nèi)部的轉(zhuǎn)移和調(diào)整，新的局部化帶不斷形成，已有的局部化帶繼續(xù)擴展、匯合，直至在坡體內(nèi)部完全貫通。在此極限狀態(tài)，處于局部化的土體單元沿其局部化帶滑移方向除承擔(dān)與其強度相適應(yīng)的載荷之外，已無能力承擔(dān)附加荷載，任何微小的擾動都可能導(dǎo)致滑體沿連續(xù)貫通的局部化帶發(fā)生大變形，最終引發(fā)邊坡模型的整體失穩(wěn)。

因此，可將局部化帶由坡底至坡頂?shù)耐耆B續(xù)貫通時刻定義為邊坡力學(xué)模型的極限平衡狀態(tài)，將局部化帶的貫通路徑定義為邊坡模型的臨界滑動面。只要準(zhǔn)確地追蹤到局部化帶自萌生至貫通的完整路徑，即可同時確定邊坡模型的極限狀態(tài)和滑動面位置。

2 局部化帶追蹤技術(shù)

2.1 局部化信息

為準(zhǔn)確追蹤局部化帶路徑、確定邊坡臨界滑面，需要發(fā)展有效的追蹤技術(shù)，即根據(jù)每個時間步邊坡有限元計算結(jié)果，通過局部化帶追蹤算法，以后處理的方式再現(xiàn)邊坡內(nèi)部局部化帶的演化狀態(tài)與擴展路徑。在局部化帶理論框架下，從邊坡有限元計算結(jié)果中提取追蹤算法所需的兩類信息。

(2)

(3)

考慮到滿足式(2)的解不唯一，當(dāng)有限元單元發(fā)生局部化失穩(wěn)時，至少存在2組潛在的局部化擴展方向。本篇采用“沿實際局部化帶切線方向，絕對位移率分量最大”原則，從多個候選局部化向量中選取局部化帶的實際擴展方向，其數(shù)學(xué)表達(dá)式為：

(4)

2.2 追蹤算法

在數(shù)值實現(xiàn)時，邊坡有限元將整個邊坡離散為一組單元的組合體，將連續(xù)的無限自由度問題變成離散的有限自由度問題。由于單元內(nèi)部的應(yīng)力場被平均化處理，由單元積分點(高斯點)的應(yīng)力近似代表，因此對于滿足式(2)的局部化單元，局部化帶出現(xiàn)的位置是不確定的，可以出現(xiàn)在局部化單元的任意區(qū)域。

為此，發(fā)展了局部追蹤算法(Local Tracking Algorithms)和全局追蹤算法(Global Tracking Algorithms)來確定局部化帶的具體擴展路徑。這兩類算法均遵循如下假定條件：局部化帶以離散線段形式在每個局部化單元中擴展；局部化帶在相鄰單元擴展時，要求保持連續(xù)性。

2.2.1 局部追蹤算法

根據(jù)前文給出的局部化信息，局部追蹤算法從根單元出發(fā)，逐個單元地追蹤局部化帶擴展路徑。局部追蹤算法示意如圖3所示。

圖3 局部追蹤算法示意

局部追蹤算法基本流程如下所示。

(1)根單元挑選。從局部化單元集合中挑選出最早進入局部化的單元，將其定義為根單元，如圖3(a)所示。

(2)局部化帶啟動。從根單元的中心點出發(fā)，局部化帶沿局部化擴展方向擴展至單元邊界，如圖3(b)所示。

(3)局部化帶擴展。核查局部化帶端點兩側(cè)單元是否為局部化單元，如果是，從局部化帶端點出發(fā)，繼續(xù)沿對應(yīng)單元的局部化擴展方向擴展局部化帶，直至到達(dá)邊界或局部化帶端點兩側(cè)單元為非局部化單元，如圖3(c)所示。

局部追蹤算法在處理單條局部化帶擴展時十分簡單，但在處理多條局部化帶擴展時非常繁雜[14-16]，缺少魯棒性，并且構(gòu)造的局部化帶路徑也不光滑。

2.2.2 全局追蹤算法

全局追蹤算法將局部化帶路徑追蹤問題轉(zhuǎn)換為求解穩(wěn)態(tài)熱傳導(dǎo)類問題，由等溫線獲得所有潛在局部化帶，從中確定局部化帶擴展路徑，如圖4所示。

圖4 全局追蹤算法示意

全局追蹤算法分為3個部分：

(5)

(2)求解向量場T(x)的包絡(luò)線(潛在局部化帶)。邊坡內(nèi)部的局部化帶可以看作無厚度的曲線集合{Si}，曲線Si上任意坐標(biāo)x處的切線方向為局部化擴展方向T(x)。因此，從數(shù)學(xué)角度看，追蹤局部化帶路徑等同于求解向量場T(x)的包絡(luò)線。

在邊坡區(qū)域Ω內(nèi)定義標(biāo)量函數(shù)θ(x)，使其梯度與T(x)始終正交，即：

T(x)·?θ(x)=0x∈Ω

(6)

考慮到函數(shù)θ(x)的等值線(θ(x)為常數(shù))垂直于梯度?θ(x)，可知θ(x)的等值線即為向量場T(x)的包絡(luò)線，物理上代表所有潛在局部化帶，如圖4(b)所示。

為了求解θ(x)，在式(6)等號兩邊同乘T(x)，將式(6)擴展為一組熱傳導(dǎo)類偏微分方程組，即：

(7)

式中：υ——垂直于?qΩ的單位向量；

θd——在Dirichlet邊界條件上的值。

為獲得θ(x)的非零解，應(yīng)至少對位于不同等值線上的2個節(jié)點設(shè)置θd。并且，考慮到K(x)存在奇異性，需要將式(7)中的K(x)=T(x)?T(x)修正為：

K(x)=T(x)?T(x)+∈I

(8)

式中：I——單位張量；

∈——數(shù)值很小的各項同性熱導(dǎo)率。

式(7)表明，標(biāo)量函數(shù)θ(x)可以類比為溫度場，-(K(x)?θ(x))為熱通量，K(x)為各向異性熱傳導(dǎo)張量。因此，追蹤局部化帶路徑問題轉(zhuǎn)變?yōu)榍蠼鉄o內(nèi)部熱源和零熱通量輸入的穩(wěn)態(tài)熱傳導(dǎo)類問題，通過求解溫度場的等溫線，獲得所有潛在局部化帶。

(3)確定局部化帶擴展路徑。在獲得所有潛在局部化帶后，需要從中確定局部化帶擴展路徑。從未追蹤的局部化單元集合中挑選出最早進入局部化的單元(根單元)，將通過根單元中心點的等溫線標(biāo)記為潛在局部化帶，如圖4(c)所示。追蹤與潛在局部化帶相交的所有局部化單元，將其標(biāo)記為已追蹤局部化單元，將潛在局部化帶與這些局部化單元的相交線段標(biāo)記為局部化帶擴展路徑。

如果局部化單元集合中仍然存在未被追蹤的局部化單元，則繼續(xù)從未追蹤的局部化單元中挑選下一個根單元，開始下一條局部化帶路徑追蹤，直至所有局部化單元均已追蹤完畢。

3 邊坡臨界滑面判識流程

與局部追蹤算法相比，全局追蹤算法理論推導(dǎo)嚴(yán)密，更適合于多局部化帶路徑追蹤問題。因此，可以基于全局追蹤算法來判識邊坡臨界滑面，具體實施流程如下。

(1)作為后處理程序，邊坡臨界滑面判識程序在每個或每幾個時間步結(jié)束后調(diào)用，通過追蹤當(dāng)前荷載下局部化帶擴展路徑來判識邊坡滑面狀態(tài)。

(6)核查局部化單元是否全部標(biāo)記為已追蹤狀態(tài)。如果不是，轉(zhuǎn)步驟(5)，開始下一條局部化帶路徑追蹤。如果是，則完成全部局部化帶路徑追蹤。

(7)根據(jù)各條局部化帶路徑，判斷局部化帶貫通情況。如果某條局部化帶由坡底至坡頂完全連續(xù)貫通，則整個邊坡處于極限平衡狀態(tài)，該條局部化帶擴展路徑為臨界滑面，邊坡穩(wěn)定性計算結(jié)束；如果沒有，則轉(zhuǎn)步驟(1)，開始下一個時間步長加載及邊坡穩(wěn)定性計算。

4 有限元數(shù)值實現(xiàn)

全局追蹤算法的數(shù)值實現(xiàn)流程：使用邊坡有限元模型網(wǎng)格，將給定的邊坡區(qū)域Ω離散為nelem個有限單元和nnode個有限節(jié)點。采用有限元方法將穩(wěn)態(tài)熱傳導(dǎo)控制方程(7)離散化[19]，導(dǎo)出其對應(yīng)的離散形式，即：

尋找：

以致：

(10)

式中：Ni(x)——標(biāo)準(zhǔn)形函數(shù)；

K——相應(yīng)的剛度矩陣。

作為后處理程序，式(10)需要在邊坡穩(wěn)定性計算的每個或每幾個時間步完成后調(diào)用，并且根據(jù)前文描述的追蹤算法判識邊坡臨界滑面來確定每個局部化單元中局部化帶的具體擴展位置，見步驟(4)～(7)。

需要注意的是，如果未指定邊界溫度(?θΩ=???qΩ=?Ω)，則K矩陣的秩為nnode-1。因此，必須在至少一個節(jié)點處指定溫度，以確保有限元方程(10)具有唯一解。此外，為了避免溫度場θ出現(xiàn)常數(shù)解(這將無法區(qū)分不同的等溫線)，應(yīng)在另一個節(jié)點處指定溫度。如何指定溫度值不影響局部化帶位置的判定，只要不在同一條等溫線上的2個點上施加2種不同的溫度值即可。

5 算例驗證

為驗證全局追蹤算法的有效性，分別采用強度折減和位移加載方式，開展兩類邊坡臨界滑面追蹤數(shù)值分析。土體材料均為理想彈塑性材料，屈服準(zhǔn)則為平面應(yīng)變條件下摩爾-庫侖結(jié)合DP準(zhǔn)則(DP4準(zhǔn)則)。

5.1 算例1：強度折減

采用文獻[4，14]中的邊坡算例：邊坡坡高H=20 m，坡角β=26.57°，楊氏模量E=100 MPa，泊松比v=0.3，容重γ=20 kN/m3，黏聚力c=10 kPa，內(nèi)摩擦角φ=20°。土體材料為理想彈塑性材料，強度準(zhǔn)則采用平面應(yīng)變條件下摩爾-庫侖結(jié)合DP準(zhǔn)則(DP4準(zhǔn)則)。算例1邊坡尺寸及網(wǎng)格劃分如圖5所示，邊坡左、右兩側(cè)邊界為法向約束，底邊為雙向固定約束，共劃分四邊形單元5 686個、節(jié)點5 878個。采用有限元強度折減方法，按照c/ω，tanφ/ω方式不斷折減土體材料強度參數(shù)，模擬邊坡漸進破壞過程及其位移場演化信息，直至邊坡發(fā)生整體失穩(wěn)。其中，ω為折減系數(shù)，其初始值ω0=1.0，步長△ω=0.002。

圖5 算例1邊坡尺寸及網(wǎng)格劃分

折減系數(shù)ω=1.384時等效塑性應(yīng)變云圖如圖6所示。由圖6可以大致確定滑動面的分布范圍，但無法精確定位。

圖6 折減系數(shù)ω=1.384時等效塑性應(yīng)變云圖

采用文獻[12]中局部追蹤技術(shù)，折減系數(shù)ω=1.384時局部化帶擴展路徑如圖7所示[12]。由圖7可以看出，局部追蹤技術(shù)給出的局部化區(qū)域要比塑性云圖的范圍小得多，主要集中在邊坡的局部破壞區(qū)域。然而在坡角及坡肩附近，多條局部化帶相互交織，并不光滑，這給臨界局部化帶位置的判識帶來了困難。此外，這些局部化滑帶寬度與有限元網(wǎng)格尺寸相關(guān)，通常占據(jù)2～3個單元。

圖7 折減系數(shù)ω=1.384時局部化帶擴展路徑

折減系數(shù)ω=1.384時，采用全局追蹤技術(shù)獲得的溫度場等值線和局部化帶擴展路徑如圖8所示。與局部追蹤技術(shù)相比，全局追蹤技術(shù)同時利用局部化單元的擴展方向和非局部化單元可能的擴展方向求解邊坡向量場的包絡(luò)線，從而從所有潛在局部化帶中確定局部化帶擴展路徑。這種方法追蹤滑面的魯棒性更高，不僅主級滑面(紅色線條)為無厚度的光滑曲線，次級滑面(藍(lán)色線條)清晰，而且各級滑面之間以層狀方式平行分布，不再出現(xiàn)滑面彼此交叉的情況。

圖8 折減系數(shù)ω=1.384時采用全局追蹤技術(shù)追蹤效果

通過以上對比可以看出，借助塑性應(yīng)變云圖的傳統(tǒng)臨界滑動面確定方法具有很大的局限性，僅能大致確定滑動面的分布范圍；局部追蹤技術(shù)雖然能夠給出更具體的局部化區(qū)域，但在坡角及坡肩附近的局部化帶相互交織，影響了臨界局部化帶位置的判識；而全局追蹤技術(shù)則能夠克服這些問題，提供更準(zhǔn)確、更全面的局部化帶擴展路徑信息。

5.2 算例2：位移加載

采用文獻[14，20]中的邊坡算例：坡高H=10 m，坡角β= 45°，楊氏模量E= 10 MPa，泊松比v= 0.4，容重γ= 20 kN/m3，黏聚力c= 20 kPa，內(nèi)摩擦角φ= 30°。邊坡幾何尺寸及網(wǎng)格劃分如圖9所示，邊坡右側(cè)邊界為法向約束，底邊為雙向固定約束。除重力載荷G外，邊坡在坡頂承受豎直向下位移u作用，F(xiàn)為豎向位移u產(chǎn)生的豎向反力。重力載荷通過初始應(yīng)力場體現(xiàn)，所產(chǎn)生的位移不計入豎向位移u中。

圖9 算例2邊坡幾何尺寸及網(wǎng)格劃分

邊坡在失穩(wěn)臨界狀態(tài)下(豎向位移u=0.14 m)的等效塑性應(yīng)變云圖、采用局部追蹤技術(shù)和全局追蹤技術(shù)獲得的局部化帶擴展路徑如圖10～12所示。由圖10～12可以看出：通過塑性應(yīng)變云圖呈現(xiàn)出2條條帶狀分布的貫通滑動面，無法準(zhǔn)確判斷臨界滑動面的位置；局部追蹤技術(shù)和全局追蹤技術(shù)都能夠判斷臨界滑動面的位置；然而，局部追蹤技術(shù)給出的臨界滑面不夠準(zhǔn)確，其路徑不光滑，且在坡體下半?yún)^(qū)域內(nèi)的局部滑帶也呈現(xiàn)出條帶狀分布，占據(jù)了2～3個單元的厚度；與局部追蹤技術(shù)相比，全局追蹤技術(shù)追蹤的主級滑面(紅色線條)以無厚度的光滑曲線的方式呈現(xiàn)，沒有出現(xiàn)滑面互相交叉以及滑面條帶化的情況，因此，全局追蹤技術(shù)可以為后續(xù)的邊坡工程治理提供更有效地支撐。

圖11 u=0.14 m時局部化帶擴展路徑

圖12 u=0.14 m時采用全局追蹤技術(shù)追蹤效果

通過對比分析，可以得出全局追蹤技術(shù)在追蹤邊坡滑動面方面具有更高的魯棒性和準(zhǔn)確性，能夠以光滑曲線的方式呈現(xiàn)局部化帶路徑，克服了傳統(tǒng)方法和局部追蹤技術(shù)的局限性，從而為后續(xù)的邊坡工程治理提供了重要的分析支撐手段。

6 結(jié)論

(1)在局部化帶理論框架下，將全局追蹤算法與局部化帶貫通判據(jù)相結(jié)合，提出了一種基于全局追蹤技術(shù)的露天煤礦邊坡臨界滑面確定方法。該方法將追蹤局部化帶路徑的問題轉(zhuǎn)化為求解無內(nèi)部熱源和零熱通量輸入的穩(wěn)態(tài)熱傳導(dǎo)問題，通過追蹤經(jīng)過局部化單元中心的等溫線來辨識局部化帶路徑的擴展及貫通狀況。

(2)算例結(jié)果比對分析表明：借助塑性應(yīng)變云圖的傳統(tǒng)臨界滑動面確定方法存在很大的局限性，只能大致確定滑動面的分布范圍，而無法準(zhǔn)確定位；局部追蹤技術(shù)雖然能夠給出更具體的局部化區(qū)域，但局部化帶經(jīng)常相互交織，影響了臨界局部化帶位置的精準(zhǔn)判識，此外，局部追蹤技術(shù)給出的臨界滑面不夠準(zhǔn)確，其路徑不光滑，經(jīng)常呈現(xiàn)出條帶狀分布特點；基于全局追蹤技術(shù)的邊坡臨界滑面確定方法具有更好的魯棒性，它能以光滑曲線的方式呈現(xiàn)滑面，不會出現(xiàn)滑面互相交叉以及滑面條帶化的情況。這種方法更適合多條局部化帶同時擴展情況下的臨界滑面的確定。

(3)提出的全局追蹤技術(shù)結(jié)合局部化帶貫通判據(jù)的方法能夠有效地確定邊坡的臨界滑面，克服了傳統(tǒng)方法的局限性。