一種基于強度折減法的自適應安全系數(shù)算法研究①

2020-03-25 01:46:34伍禮杰鄧紅衛(wèi)張亞南

礦冶工程 2020年1期

伍禮杰，鄧紅衛(wèi)，張亞南

（中南大學資源與安全工程學院，湖南長沙410083）

穩(wěn)定性分析是邊坡工程建設、評估、治理中的重要內容，目前主要通過理論計算、試驗研究和數(shù)值模擬等方法進行研究，其中利用強度折減法計算安全系數(shù)是目前最常用的數(shù)值模擬方法之一。保證安全系數(shù)計算準確的關鍵是精準判定邊坡的狀態(tài)，目前在邊坡穩(wěn)定性分析中常使用潛在破壞面貫通［1－2］、特征點位移突變［3］、能量突變［4－6］和數(shù)值計算不收斂［7］作為邊坡失穩(wěn)的判據(jù)，但前3 種判據(jù)對邊坡臨界狀態(tài)的判定依賴在折減系數(shù)等增量變化下對判據(jù)變化的觀察，且貫通塑性區(qū)發(fā)展部位的判斷依賴于人的經驗，在事先不能確定安全系數(shù)范圍的情況下要想獲得高精度的解，計算過程繁雜、效率低，即使采取增量細化搜索策略［8］進行優(yōu)化，其計算量也非常龐大。數(shù)值計算不收斂判據(jù)的生效并不依賴折減系數(shù)等增量變化下某些觀測量表現(xiàn)出的突變性和漸變貫通性，并能通過交叉二分法提高安全系數(shù)的搜索效率，使用方便且計算效率較高，因此在大型數(shù)值計算中具有良好的實用性。基于數(shù)值計算不收斂判據(jù)，許多研究通過各種方法來提高算法效率［9－12］，這在一定程度上對強度折減安全系數(shù)算法的搜索區(qū)間、搜索機制、初始應力狀態(tài)回歸等方面進行了優(yōu)化，但對數(shù)值計算收斂判據(jù)的標準依舊采用人為劃定的收斂限值和迭代計算上限，收斂限值過大會影響安全系數(shù)的準確性，過小則會提高計算收斂難度導致計算量過大甚至無法收斂；迭代計算上限過大或過小會導致邊坡狀態(tài)的誤判，且過大的迭代計算上限還會導致計算量過大。尋找位于安全系數(shù)求解準確性和效率兩者平衡點間的迭代計算上限值需要豐富的經驗，所得結果會因人而異，故在保證安全系數(shù)計算準確性前提下提高求解效率的關鍵是采用合理數(shù)值計算不收斂判據(jù)。

本文從FLAC3D內置安全系數(shù)算法（以下簡稱內置算法）出發(fā)，針對內置算法的不收斂判據(jù)存在的一些不足，通過對安全系數(shù)搜索和計算流程中的系統(tǒng)失穩(wěn)狀態(tài)判定的優(yōu)化研究，建立了一種原理簡單、方便易用、準確高效的自適應安全系數(shù)算法（以下簡稱自適應算法）。

1 FLAC3D內置算法的機制與問題分析

1.1 內置算法原理

FLAC3D內置算法采用強度折減法計算邊坡的安全系數(shù)，即不斷折減材料的剪切強度以獲取邊坡臨界平衡狀態(tài)，此時的材料折減系數(shù)即邊坡穩(wěn)定性安全系數(shù)FS，其原理的表達式為：

式中c、φ 分別為材料的粘聚力和內摩擦角；Ftrial為折減系數(shù)；ctrial、φtrial分別為材料強度折減后的粘聚力和內摩擦角。

大多數(shù)工程邊坡特別是失穩(wěn)后亟需治理和處治后，要求邊坡的安全系數(shù)介于0.5～2.0 之間，因此本文討論的安全系數(shù)范圍為［0.5，2.0］。FLAC3D5.0 軟件的內置算法機制可以視為運用Ftrial搜索區(qū)間為［0.5，2.0］且首次Ftrial為1.0 的特殊二分法來逼近邊坡臨界狀態(tài)的折減系數(shù)，其典型的求解流程如圖1 所示，具體步驟如下：

1）初始計算條件設置，包括材料本構關系、性質和模型邊界條件等參數(shù)的設置，重力加速度的施加等。

2）初始平衡計算，生成初始地應力，并對數(shù)值計算產生的變形位移、變形速度和塑性屈服區(qū)進行歸零處置，獲得模型初始應力狀態(tài)文件。

3）系統(tǒng)特征響應時步Nr的計算，以摩爾-庫侖本構模型為例，此計算過程的實質是給模型粘聚力和抗拉強度一個極大值，使系統(tǒng)內部應力狀態(tài)產生巨大變化，默認給系統(tǒng)應力狀態(tài)一擾動因子為2 的擾動，計算系統(tǒng)將歷經多少時步回歸平衡，這個時步值即為Nr，反映系統(tǒng)的求解計算難度。

圖1 內置算法典型流程

4）采用交叉法猜測強度折減系數(shù)Ftrial的上下限值［v1，v2］，并用二分法不斷試算逼近模型臨界狀態(tài)的FS。每經過Nr時步的強度折減迭代計算都會進行系統(tǒng)狀態(tài)判定，當不平衡力比率小于收斂限值Rd（軟件默認為1.0×10－5）時，則認為數(shù)值計算結果收斂，系統(tǒng)狀態(tài)是穩(wěn)定的，更新穩(wěn)定狀態(tài)文件FOSStable.f3sav。當某次迭代計算Nr時步后系統(tǒng)的平均力比率的變化量小于10%或者其變化量一直大于10%但迭代時步超過6Nr，都認為數(shù)值計算結果不收斂，系統(tǒng)處于失穩(wěn)狀態(tài)，更新失穩(wěn)狀態(tài)文件FOSUnstable.f3sav。以上條件都不符合時，系統(tǒng)處于未確認狀態(tài)，進入下一輪Nr時步的迭代計算，直到滿足條件后才會退出強度折減迭代計算。

1.2 內置算法存在的問題分析

FLAC3D5.0 內置算法的本質是傳統(tǒng)強度折減法，雖然有著眾多優(yōu)點，但在計算邊坡安全系數(shù)時依舊存在有待改進之處，比如對粘聚力和內摩擦角的同比例折減策略并不完全符合實際邊坡失穩(wěn)破壞過程中粘聚力和內摩擦角的變化［13－14］，邊坡失穩(wěn)狀態(tài)難以精準判定等。圖2 為內置算法計算過程中不平衡力比率和折減參數(shù)的時變曲線，可以看出在時步（6.5～13）×104計算階段的三次強度折減中，計算結果不再收斂，但算法并不能及時精準識別這一系統(tǒng)狀態(tài)。結合內置算法的求解流程分析可知，以特征響應時步Nr作為分析單位過于保守，對系統(tǒng)失穩(wěn)狀態(tài)的識別度較低，在采用不平衡力比率判定系統(tǒng)狀態(tài)時，無法及時有效識別系統(tǒng)失穩(wěn)狀態(tài)使得求解時間變長，效率變低，尤其對于尺寸大、網(wǎng)格密度高和工況復雜的大型數(shù)值模型計算，安全系數(shù)求解時間會隨著求解難度的提高明顯變長，而且限于軟件權限，內置算法無法對邊坡折減過程進行連續(xù)監(jiān)測以及在大變形模式下計算安全系數(shù)。

圖2 某例內置算法求解不平衡力比率及折減系數(shù)時變曲線

2 自適應算法的機制與實現(xiàn)

2.1 自適應算法的求解機制

根據(jù)上述分析可知FLAC3D內置算法求解效率較低是因為計算過程中對系統(tǒng)失穩(wěn)狀態(tài)的判定比較保守，不能充分利用強度折減迭代計算結果中蘊含的信息，常存在已經可以判斷系統(tǒng)處于失穩(wěn)狀態(tài)但軟件仍繼續(xù)計算的情況。本文在進行這種數(shù)值計算時發(fā)現(xiàn)所有的失穩(wěn)狀態(tài)下的不平衡力比率時變曲線都符合前期迅速下降、后期平穩(wěn)波動的規(guī)律，具體如圖3 所示。可以看出Ⅰ區(qū)域曲線下降速度快、幅度明顯，計算收斂迅速，Ⅱ區(qū)域曲線在某一穩(wěn)定值上下波動，但變化幅值極小，一般而言整個計算周期不平衡力比率最小值會出現(xiàn)在這一區(qū)域，Ⅲ區(qū)域曲線以極緩極小幅度上升，表明數(shù)值計算不再收斂，系統(tǒng)失穩(wěn)。

圖3 某例求解不平衡力比率時變曲線

為解決內置算法計算效率低的問題，本文利用不平衡力比率變化量來表征強度折減計算過程中失穩(wěn)狀態(tài)下的不平衡比率時變曲線的變化規(guī)律，通過比特征響應時步更小的分析單位來觀測折減計算的收斂情況，提出一種改進的安全系數(shù)求解算法——自適應算法，其主要機制如下：

以NS取100 時步為分析單位，提取n·NS時步后的不平衡力比率Rn，并與（n－1）·NS時步后的不平衡力比率Rn－1做差值，得到當前第n 次分析計算的不平衡力比率變化量Vn：

Vn反映出第n 次分析計算過程中不平衡力比率的收斂速度，分析Vn在不平衡力比率時變曲線不同階段的特性，易知：

1）第n 次分析計算處于Ⅰ區(qū)域時，Vn＜0，且其值隨著時步增加而增加，其絕對值隨著時步增加而減小；

2）第n 次分析計算處于Ⅱ區(qū)域時，Vn值有正有負，其值在0 上下以極小幅值波動；

3）第n 次分析計算處于Ⅲ區(qū)域時，Vn＞0，Vn值穩(wěn)定且極小。

從強度折減迭代計算開始到不平衡力比率最小值出現(xiàn)結束，Vn的絕對值逐漸減小且逼近0，而失穩(wěn)狀態(tài)的實質就是系統(tǒng)的不平衡力比率無法在實際可行的有限時步內達到收斂限值以下。因此為識別系統(tǒng)失穩(wěn)狀態(tài)，理論上可以通過Vn來預測系統(tǒng)當前達到收斂限值需要的最少計算時步Ni，人為設定Ni的限值（取1.0×105），當Ni超過限值時，認為系統(tǒng)不具備在實際可行的有限時步內達到穩(wěn)定狀態(tài)的可能性，即系統(tǒng)處于失穩(wěn)狀態(tài)。

定義一個新的變量目標收斂特征比Tn，用以表征當前循環(huán)不平衡力比率變化量占比目標差值Rn－Rd的多少，其計算如式（4）所示，相應的Tn的限值D1可由式（5）確定，當Tn＞D1時，說明系統(tǒng)處于失穩(wěn)狀態(tài)。

根據(jù)Vn的變化規(guī)律可知，當前觀測到繼續(xù)計算Ni時步后無法收斂，意味著實際這個預測的判斷區(qū)間長度是不止1.0×105計算時步的。與內置算法的最大判斷區(qū)間長度3.0×105計算時步（實際絕大多數(shù)情況下，這個判斷區(qū)間是達不到的）相比，兩者處于同一數(shù)量級，因此出于計算效率和實用性的考慮以1.0×105為判斷區(qū)間是合理的。由于本算法的分析單位NS相較于內置算法很小，其更新頻率高、能更快識別出系統(tǒng)失穩(wěn)狀態(tài)，故這種判斷方法效率更高。所有采用數(shù)值計算收斂判據(jù)的強度折減安全系數(shù)算法結構都是折減參數(shù)循環(huán)包含著折減計算循環(huán)，自適應算法能提高每一個系統(tǒng)狀態(tài)失穩(wěn)的計算循環(huán)，從而明顯提高安全系數(shù)求解效率。

2.2 自適應算法實現(xiàn)與計算流程

自適應算法采用NS作為折減計算循環(huán)單位，由于其計算時步較短，使得其失穩(wěn)狀態(tài)識別功能的適應性較弱，受限于網(wǎng)格數(shù)量、模型條件和工況條件等，可能會出現(xiàn)計算奇異點對這個識別功能帶來干擾，計算結果會失準。為了避免出現(xiàn)上述情況，本文設計了多重判斷機制來避免識別功能陷入計算奇異點的陷阱。執(zhí)行的失穩(wěn)狀態(tài)判定條件有：

1）判定的前置條件：①某次計算循環(huán)的預測迭代次數(shù)超過Ni且其不平衡力比率變化量是負值，即D＜Tn＜0。②連續(xù)5 次計算循環(huán)的不平衡力比率變化量是正值，即min [ Tn－4，Tn－3，Tn－2，Tn－1，Tn]＞0。③以10NS時步為間隔，某次最大不平衡力比率Un不再減小，即Un－Un－9＞0。

當條件①、②或③累計出現(xiàn)次數(shù)超過1 次后，視為前置條件已滿足。

2）單個分析單位的最終判定條件：當前置條件達到后，連續(xù)3 次小循環(huán)的預測迭代次數(shù)超過Ni，即min [ Tn－2，Tn－1，Tn]＞0，認為當前系統(tǒng)為失穩(wěn)狀態(tài)。

上述3 種最終判定條件只要滿足其中一種就判定系統(tǒng)處于失穩(wěn)狀態(tài)，系統(tǒng)穩(wěn)定狀態(tài)的判定則與內置算法一樣，以不平衡力比率小于收斂限值為條件，一旦失穩(wěn)或穩(wěn)定判定條件生效，算法就會退出當前的強度折減計算。其他情況則屬于系統(tǒng)狀態(tài)未確定，算法繼續(xù)強度折減迭代計算，直至系統(tǒng)狀態(tài)可以判定。

這種判定機制的建立是基于不平衡力比率時變曲線的變化規(guī)律，前置條件的判斷是為了保證計算階段已處于平穩(wěn)波動區(qū)，條件2）則是通過多次識別系統(tǒng)的失穩(wěn)狀態(tài)，進一步精確判斷系統(tǒng)狀態(tài)，降低誤識別的機率；條件3）和4）則是為了充分利用已經計算得到的不平衡力比率數(shù)據(jù)，在保證失穩(wěn)狀態(tài)識別功能完備前提下彈性調整判定限值，避免進入計算奇異點，陷入死循環(huán)。

自適應算法采用軟件提供的FISH 語言編寫實現(xiàn)，系統(tǒng)狀態(tài)判定采用前述的判定機制，具體流程如圖4 所示。

圖4 自適應算法典型流程圖

3 自適應算法的驗證與分析

3.1 自適應算法的可靠性驗證

為對比內置算法與自適應算法在求解結果上的差異和求解效率上的優(yōu)劣，參考文獻［15］中的模型數(shù)據(jù)和參數(shù)設置了4 個具有代表性的邊坡模型來進行對比分析，其邊坡具體尺寸如圖5 所示，模型寬度統(tǒng)一取5.0 m，其材料參數(shù)如表1 所示，安全系數(shù)都處于［0.5，2.0］區(qū)間內。

圖5 4 種邊坡網(wǎng)格構形

表1 材料參數(shù)

建模時采用Midas GTS NX 進行網(wǎng)格劃分，網(wǎng)格類型為四面體＋六面體的混合網(wǎng)格，網(wǎng)格尺寸為1.0 m，模型底部施加3 方向位移約束的邊界條件，四周施加垂直面的方向位移約束的邊界條件，上部施加自由邊界。為保證計算效果的可靠性，兩種算法對同一邊坡進行安全系數(shù)計算時采用的網(wǎng)格模型、重力加速度、材料參數(shù)和軟硬件環(huán)境保持相同，求解精度均取0.001。

表2 和表3 分別為自適應算法與內置算法求解的安全系數(shù)和所需的計算時步，可以看出兩種算法求解結果差異非常小，4 個邊坡自適應算法結果與內置算法結果的差值比分別為0.52‰，1.86‰，0.00‰和0.00‰，但在計算時步上自適應算法與內置算法相比，4個邊坡的計算時步分別減少了42.24%，65.75%，39.90%和18.98%。這說明自適應算法不僅求解結果是準確的，而且在求解效率上是大于內置算法的，因此自適應算法一定程度上改進了內置算法存在的系統(tǒng)失穩(wěn)狀態(tài)識別效率低的缺點。

表2 安全系數(shù)對比

表3 計算時步對比

3.2 網(wǎng)格密度對自適應算法求解結果的影響分析

為進一步研究在大型數(shù)值計算中網(wǎng)格密度對兩種算法的計算結果差異的影響，以邊坡（c）為基礎，劃分0.3～1.2 m 共10 種不同網(wǎng)格尺寸的模型，獲得了不同網(wǎng)格密度下兩種算法所得的安全系數(shù)、求解時步和求解時間，具體結果如表4 和圖6～8 所示。

由表4 和圖6 可知，2 種算法求解的安全系數(shù)曲線的變化規(guī)律基本相同，數(shù)據(jù)差異很小，始終保持在1.1%以下，說明網(wǎng)格密度對自適應算法計算結果的準確性基本上沒有影響。由表4 和圖7 可知，隨著網(wǎng)格密度增大，內置算法的求解時步總體上不斷增大，而自適應算法所需的求解時步能夠始終保持在一定的數(shù)值范圍內，對比內置算法平均能夠減少55.62%的求解時步，自適應算法在計算時步上的優(yōu)化效果顯著。由表4 和圖8 看出，隨著網(wǎng)格密度增加，2 種算法的求解時間都在增加，但自適應算法的增幅小于內置算法，相比內置算法，自適應算法平均能夠減少50.10%的求解時間。因此自適應算法在計算效率上比內置算法具有更加明顯的優(yōu)勢。