基于光滑粒子流體動力學方法的土坡滑面確定與分析

胡嫚， 劉秋強， 鮑安紅， 楊嘯宇， 高濤

1. 西南大學 工程技術學院， 重慶 400715； 2. 中華人民共和國應急管理部， 北京 100054；3. 西南交通大學 土木工程學院， 成都 610031

在巖土工程領域， 邊坡穩定性問題是最基本的問題之一， 廣泛存在于天然邊坡、 人工邊坡、 擋土結構、 水庫大壩等結構的整體及局部穩定性分析當中， 其中， 確定邊坡潛在滑動面的位置和形狀對邊坡穩定性分析結果具有決定性和控制性作用[1]．

關于邊坡危險滑動面的獲取與確定問題目前已經開展了廣泛的研究， 取得了較多的成果， 邊坡危險滑動面的獲取分析理論大致有極限平衡理論、 滑移線場理論、 極限分析理論、 數值分析方法和差分理論等5種方法[2]． 其中以極限平衡法[3]和數值分析方法[4]應用最為廣泛． 極限平衡法是通過對邊坡內某一給定形狀確定的滑面， 計算該滑面對應的安全系數， 進而搜索邊坡不同位置所有可能的滑面， 找出安全系數最小的滑面作為潛在最危險滑面[5]． 數值分析方法與強度折減概念相結合在目前也開展了許多研究， 包括參數折減模型[6-7]、 失穩判據[8]、 滑面推算[9]、 安全系數[10]、 計算精度[11]、 各種特殊工況下的應用[12-13]等， 已經成為邊坡滑移面推算與搜索的主要數值分析手段． 強度折減理論通過逐步折減邊坡巖土強度參數取值， 最終促使邊坡逼近極限平衡狀態， 極限平衡狀態下強度參數折減倍數即為安全系數， 此時滑移面可通過相應判據規律確定． 由于基于有限元法或有限差分法的強度折減法不需要提前設定滑動面的位置和形狀， 從而顯著地減少了推算過程中的人工干預， 所得滑面結果更加精確． 因此， 部分學者采用此類方法進行邊坡次級滑動面的搜索與確定， 其中包括了趙尚毅等[14]采用最大剪應變增量原則確定邊坡滑動面及求解最小安全系數； 孫冠華、 聶治豹等[15-16]以等效塑性應變為判據計算邊坡滑動面． 此外， 基于應力影響系數法的邊坡滑面確定方法及基于極限應變的判斷原則， 以能量突變為依據的判定方法等方法也得到了發展． 但有限元強度折減法也存在一些難以避免的不足， 如極限平衡狀態的判據問題、 從塑性區中提取滑面的算法問題、 有限元計算中的網格畸變問題等． 由于有限元方法在模擬計算大變形區域過程中會產生數值不穩定性， 從初始狀態到后續大變形的整個變形過程難以連續求解， 因此， 計算流體動力學(CFD)建模和離散建模(如DEM)等方法逐漸發展了起來． 而CFD計算要求材料是特殊類型的流體[17]， 難以解決靜態變形問題， DEM離散建模也不適合處理基于連續介質近似的土體本構模型[18]．

綜上所述， 確定土坡滑面位置的演化， 進而分析土坡的穩定性與動態滑動， 是目前傳統的極限平衡法和有限元方法難以解決的問題， 然而這個過程的求解顯然十分重要， 整個變形過程中包含了許多信息(如破壞區、 滑面演化、 破壞過程、 破壞時間、 推移距離、 沖擊力等)． 因此， 以土坡從小變形到大變形連續求解整個變形過程的視角出發， 探究土坡滑面的位置和演化過程具有重要意義．

無網格法作為求解整個變形過程的一種方法， 是一種有潛力、 有效的大變形求解方法． SPH方法屬于拉格朗日方法， 可以解決大變形問題而不產生網格畸變， 在巖土工程問題中有較好的發展潛力． 此外， 由于SPH方法基于連續介質近似， 能夠較好地表達控制方程和現有巖土材料的本構模型， 進而求解巖土材料的整個變形過程． 目前SPH方法已被用于解決許多相關巖土問題并取得了一些顯著的成果[19-22]， 但目前的研究中尚缺少基于SPH從整個變形過程的視角出發， 確定土坡的滑面位置的相關研究．

基于土體理想彈塑性本構模型， 根據最大剪切應變率原則選取最大剪切應變率SPH粒子， 多組最大剪切應變率粒子貫通可確定滑面位置與形態． 文中采用了Fortran編程進行SPH理論計算與滑面確定， 所確定的滑面與極限平衡法確定滑面進行了比較， 并分析了土坡穩定性與土坡形變間的關系， 此方法為土坡滑面的確定與分析提供了一種新的方式．

1 理論與方法

1.1 SPH的基本理論

在SPH方法中， 計算對象用一系列粒子的集合表示． 一方面， 粒子的運動可以單獨求解； 另一方面， SPH方法通過采用了獨特的插值理論， 將一系列粒子結合起來以表示連續介質的變形． 上述插值理論包括兩個關鍵步驟： 核近似和粒子近似． 使用核近似法， 對于參考粒子α， 基于支持域內的相鄰粒子β， 引用光滑函數W插值計算， 如圖1所示，h是光滑域半徑， 影響域半徑的h由粒子間的初始距離r乘以參數kd得到。

圖1 SPH方法光滑域與光滑函數

核近似場函數表示為[23]

(1)

光滑函數W用來表示臨近的β粒子對α粒子的影響． 本文中三次樣條函數取為光滑函數[24]．

(2)

其中，R=|x-x′|/h， 在二維和三維空間，αd分別取15/(7πh2)和2πh3．

另外， 粒子近似將粒子的信息用光滑域內的其他粒子表示， 其特征如下：

(3)

其中Wij=W(xi-xj，h)，m和ρ分別表示粒子的質量和密度．

1.2 固體力學的SPH方法

本文中使用的控制方程是基于固體力學的方法， 連續性方程和運動方程可以表示為

(4)

(5)

其中ui是速度矢量，ρ是密度，σij是應力張量和Fi是外力向量． 將SPH插值理論應用于狀態方程， 上式可以寫為

(6)

(7)

其中Cij為人工黏性項與人工應力項， 表示為

(8)

本文研究中將巖土體假設為理想彈塑性材料， 采用Drucker-Prager屈服準確定巖土體的塑性流動狀態， 其抗剪強度使用黏聚力和內摩擦角來表示． 當巖土體的應力狀態處于塑性屈服狀態時， 它隨屈服函數的變化而變化． 然而， 數值誤差會導致應力狀態超過屈服函數． 因此， 模擬計算過程中需要及時糾正計算粒子的應力狀態， 本文采用Chen等[25]提出的應力映射調整算法． 圖2為拉伸開裂處理示意圖． 如果計算粒子的應力狀態在計算中出現誤差， 處在了拉伸狀態， 即對應式(9)條件， 將靜水壓力應力分量進行調整， 修正到臨界拉伸狀態．

圖2 拉伸狀態處理

(9)

圖3表明了在計算模擬過程中出現計算誤差、 粒子應力狀態超出了屈服表面、 采取應力調整的過程． 按比例系數R減小偏切應力分量， 靜水應力分量I1保持不變， 如式(10)所示．

圖3 超出屈服面狀態處理

(10)

(11)

除了上述應力調整之外， 本文還采用了Jaumann應力速率來考慮剛體自旋運動的影響．

1.3 基于SPH的強度折減方法

在邊坡穩定性分析中(極限平衡法)， 安全系數通常定義為抗滑力與下滑力的比值， 因此， 以彈塑性巖土體為例， 極限平衡臨界狀態意味著滑面上巖土體均處于塑性屈服狀態． 與基于有限元法或有限差分法的強度折減法相似， 基于SPH的邊坡穩定性分析的強度折減法通過對整個土坡的巖土體抗剪強度參數c和φ進行一系列折減得到ct和φt， 如式(12)和式(13)所示．

(12)

(13)

上式中c和φ為巖土體實際抗剪強度參數；ct為強度折減后的黏聚力值；φt為強度折減后內摩擦角值；FOSt為試算的安全系數(折減倍數)． 一般初始時刻FOSt設置得足夠小， 以便系統穩定． 隨后，FOSt值逐漸增加， 直到出現貫穿的滑面并持續呈現， 邊坡失穩， 使邊坡出現臨界失穩狀態的FOSt的最終值被定義為邊坡的安全系數．

基于有限元法或有限差分法的強度折減法在確定邊坡的極限平衡狀態時有多種判別標準， 如計算收斂、 塑性區貫通、 最大位移等[26]． 與此不同的是基于SPH的強度折減法在模擬計算過程中， 可以連續求解從小應變到大變形的整個變形過程， 而且可以求解破壞區、 破壞時間、 移動距離等．

1.4 滑面分析方法

應用強度折減法確定滑動面的形態計算時， 確定極限平衡狀態與判定滑動面位置是其中的關鍵． 從連續求解由小應變到大變形的整個變形過程的視角出發， 不同于有限元方法， SPH框架下土坡在整個變形過程中確定極限平衡狀態臨界點的意義不大， 因為土坡變形的計算過程必定會趨于收斂， 故采用SPH方法重要的一步是判斷滑動面位置．

當土體承受的剪切力達到最大抗剪強度后， 隨著土體的連續變形， 土體抗剪強度便降到殘余強度， 本文的研究假設不考慮滑面上巖土體材料塑性屈服后剪切強度的減弱， 也就是說， 切面殘余剪切強度仍等于峰值剪切強度． 在判斷土坡滑面位置的分析中， 時間間隔Δt內， 每個滑面搜索區間內最大剪切應變率的粒子發生的剪切位移是最大的， 因此， 粒子最大剪切應變率可以作為土坡中滑動面判定的標準． 具體而言， 若任一土坡其穩定性不足， 即會發生變形， 進而失穩滑動． 在土體滑動發生的過程中， 首先在全局搜索區域(土坡全長或滑面可能出現范圍)中劃分單次搜索區間， 單次搜索區間的最優長度是1.1～1.2倍初始粒子間距， 即S=(1.1～1.2)d， 根據單次搜索區間的最大剪切應變率原則選取最大剪切應變率SPH粒子， 在土坡滑動進行的某一時刻開始， 多組最大剪切應變率粒子會形成貫通， 即確定了土坡的滑面位置與形態， 如圖4所示．

圖4 土坡滑面SPH粒子搜索

上述SPH方法的土坡滑面確定流程圖如圖5所示， 依次確定滑面全局搜索區域、 單次搜索區域、 開始計算和時間步增加， 直至某一時刻下多組最大剪切應變率粒子形成貫通， 即確定出土坡滑面． 值得注意的是： (1) 當土坡強度參數足夠保持穩定時， 顯然土坡處于穩定狀態， 此時通過SPH方法進行土坡模擬計算， 土坡不會發生變形與運動， 因而不會形成滑面， 若引入強度折減法降低土坡強度參數， 某一狀態下土坡會開始變形與運動， 此時才能形成滑面； (2) 某一狀態下土坡會開始變形與運動， 當多組最大剪切應變率粒子形成貫通后， 隨著土坡運動的發展， 不同時刻下最大剪切應變率粒子貫通的滑面并不一致， 他們是位置相互靠近且十分接近， 這表明某一狀態下的土坡運動過程中滑面位置并非絲毫不變， 而是有一定變化， 并且滑面位置均十分靠近而形成“滑帶”， 這與土坡失穩運動后的實際情況是一致的．

圖5 SPH方法的土坡滑面確定流程

2 算例

2.1 土坡滑面的確定

基于彈塑性本構的土體SPH模型精度在前期的研究中已經得到驗證[27]， 在這個算例中， 為了滿足經典瑞典條分法的分析條件， 計算中將剩余強度設置為峰值強度． 圖6是數值算例的尺寸示意圖， 邊坡由均質土體構成， 為了盡量消除豎直邊界條件對模擬結果的影響， 模型加大了底面尺寸建立了如圖形狀簡單的模型． 模型左側是90°垂直邊界， 右側是坡角45°的斜坡， 表1列出了本次模擬使用的參數． 模擬使用了Drucker-Prager屈服準則， 模型底部的邊界采取了無滑移邊界處理， 模擬考慮了Jaumann應力速率， 沒有考慮孔隙水壓力的影響， 在靜止土壓力的基礎上對邊坡進行重力加載．

圖6 計算模型尺寸圖

表1 數值算例模型參數

根據前述模型尺寸與參數建模計算， 在初始強度參數(c=12.31 kPa，φ=31.30°)下， 計算結果中模型隨時間推移并沒有發生明顯的形變， 說明邊坡處于穩定狀態． 按照強度折減的理念， 逐漸增加折減系數依次進行計算， 當折減倍數為0.85時， 土坡開始在重力作用下發生失穩破壞、 滑動． 圖7-圖11中橫坐標與縱坐標均表示模型尺寸．

圖7 瞬時剪切應變率隨時間的變化

圖7展示了算例在SPH框架下從小變形到大變形連續求解的過程， 顯示了土坡SPH粒子在不同時刻的瞬時剪切應變率分布， 圖7(a)-7(b)可見土坡的剪切應變率逐漸增大， 分布由分散狀態逐漸聚集為條狀， 但土體塑性屈服區尚未貫穿； 圖7(c)-7(e)土坡變形速率增大， 剪切應變率繼續增大到最大值， 形成貫穿的塑性屈服區， 即塑性剪切帶； 圖7(f)-7(g)中土坡變形速率逐漸減小， 剪切應變率也逐漸減小， 部分剪切區瞬時剪切應變率降低為零； 圖7(h)中時間到4.4 s以后土坡運動停止， 瞬時剪切應變率均為零．

圖8顯示了與圖7相對應時刻滑面判定的最大剪切應變率SPH粒子位置的分布， 圖8(a)-8(b)可見滑面尚未貫穿， 圖8(c)-8(g)顯示土坡運動過程中可確定其滑面， 圖8(g)可見土坡運動停止后， 最大剪切應變率SPH粒子位置分布呈無序狀態． 圖9提取了計算過程中最大剪切應變率出現的時間(1.7 s)， 此時最大剪切應變率值為0.298 0， 圖10與圖11分別顯示了算例中滑坡變形運動終止時的位移分布與累計塑性應變分布．

圖8 滑面位置的出現及分布

圖9 最大剪切應變率出現時間及分布(t=1.7 s， 最大剪切應變率值0.298 0)

圖10 變形終止位移圖

圖11 累計塑性應變

2.2 與極限平衡法滑面位置的對比

目前極限平衡法是邊坡穩定性分析中最成熟、 最方便、 應用最廣的方法， 常用的極限平衡法有瑞典條分法、 畢肖普條分法等， 由于幾種方法所算得的最危險滑面位置與安全系數較為接近， 本文用瑞典條分法求解滑面， 并與SPH方法進行比較分析．

圖12中， 黑色實線圓弧表示用瑞典條分法獲得的最小安全系數的圓弧， 紅色實線代表本文SPH最大剪切應變率粒子滑面， 可見本文方法與瑞典條分發得到的滑面位置均通過土坡坡角， 由于滑面形狀預先的假設， 瑞典條分法滑面是完美的圓弧形狀， 而本文方法求解的滑面為近似平滑圓弧曲線， 其圓心O1距邊坡外側更遠， 半徑更大． 本文方法得到的滑面位置更接近于使用作圖法得到的滑面位置， 即破壞面與大主應力作用方向(坡面)夾角θ=45°-φ/2．

圖12 兩種方法滑面對比

2.3 土坡形變與穩定性分析

不同于判定滑面位置問題計算收斂的確定性， 在SPH框架下土坡的穩定性分析需要首先確定整個變形過程中的極限平衡臨界狀態． 然而， SPH框架下土坡的穩定性分析的極限平衡臨界狀態， 不是上述最大剪切應變率粒子貫通形成滑面時， 而是在土坡塑性剪切區域形成過程中， 故難以明確極限平衡臨界狀態． 本文從土坡安全系數與坡體位移量的關系出發， 對土坡形變與穩定性關系作一定的探討． 選取坡面上具有代表性兩個位置的位移量： 坡角水平位移(X)和坡頂豎向位移(Y)如圖13． 通過前述基于SPH的強度折減法， 不斷折減算例中巖土體材料的強度參數黏聚力c與內摩擦角φ， 同時記錄位移量X與Y， 如表2． 不同強度參數對應安全系數下橫向位移與豎向位移水平如圖14所示．

圖13 位移量選取點位置示意圖

表2 強度折減與位移量

由圖14可知， 當安全系數接近1.1時，ΔX/δ與ΔY/δ都小于0.01， 沒有觀察到較明顯變形． 在安全系數逐漸減小的過程中， 坡頂的豎向位移和坡角的水平位移呈現逐漸增大的趨勢， 與工程邊坡中滑體后緣形成裂縫和前緣向前滑動的特征相對應起來． 反之， 當安全系數小于1.0時， 剪切帶發育明顯， 邊坡變形明顯． 當安全系數接近1.0時， 水平位移ΔX與坡面高度δ之比接近0.02， 而豎向位移ΔY與坡面長度δ之比接近0.05， 較前者大， 造成二者差異的原因可能是在邊坡初始狀態， 在粒子自重的影響下， SPH粒子間相互作用力與相互距離在計算初始階段出現一定的調整引起的． 若認為極限平衡法的安全系數是可信的， 則當安全系數為1.0時， 土坡出現明顯變形，ΔX/δ為0.018 8(<0.02)， 表明土坡在極限平衡狀態下時， 出現不明顯變形； 當安全系數小于1.0時， 土坡開始出現明顯變形(安全系數0.95，ΔX/δ為0.055 20)．

圖14 不同安全系數下橫向位移與豎向位移

綜上所述， 本文依據SPH對土坡變形過程連續求解的方式， 可以確定土坡滑面的位置與形態， 并且評估土坡穩定性與形變的關系． 傳統的極限平衡法預先設定的滑面形狀， 考慮了滿足力和力矩的平衡條件， 忽略了巖土體應力—應變的關系， 特別是變形的發展． 同時， 由于忽略了失穩過程中巖土體應力—應變的關系， 因而極限平衡法缺失了邊坡的變形信息， 但從實際工程的邊坡監測中來看， 邊坡發生較大變形、 出現較大位移時， 邊坡的平衡條件發生了變化， 穩定性系數也隨之改變； 另一方面， 極限平衡法的土坡滑動面是提前假設的， 假設的滑動面是設定的而非精確計算的． 因此， 與傳統的滑面確定方法相比， 基于SPH的滑面確定方法更能準確地預測在安全與危險邊界處巖土材料的變形．

3 結論

土坡滑面的確定對其穩定性分析起著決定性和控制性的作用． SPH方法作為一種能夠連續求解巖土體整體變形過程的方法， 在處理土體大變形方面具有較大的優勢． 本文基于SPH算法框架下土體的理想彈塑性本構提出了土坡滑面判定方法， 得出如下結論：

1) 從連續求解整個變形過程的視角出發， 在全局搜索區域中劃分單次搜索區間， 依據強度折減法與最大剪切應變率原則， 提出了滑面判定方法， 得到了理想彈塑性本構土體下土坡的滑面， 并與極限平衡條件下條分法所確定滑面進行了對比分析， 本文方法求解的滑面為近似平滑圓弧曲線， 其圓心距邊坡外側更遠， 半徑更大；

2) 依據強度折減法， 可通過SPH土坡全過程模擬方法初步分析土坡穩定性與土坡形變間的關系， 與傳統的滑面確定方法相比， 基于SPH的滑面確定方法更能準確地預測在安全與危險邊界處巖土材料的變形． 本文為土坡滑面確定與穩定性分析提供了一種新的、 有效的方式， 較極限平衡法與有限元法， 本文方法能夠從土坡滑面確定與穩定性分析中提取出土坡全過程變形信息， 是此類問題分析較好的補充．