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

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

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

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

1 FLAC3D內(nèi)置算法的機(jī)制與問題分析

1.1 內(nèi)置算法原理

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

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

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

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

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

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

圖1 內(nèi)置算法典型流程

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

1.2 內(nèi)置算法存在的問題分析

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

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

2 自適應(yīng)算法的機(jī)制與實(shí)現(xiàn)

2.1 自適應(yīng)算法的求解機(jī)制

根據(jù)上述分析可知FLAC3D內(nèi)置算法求解效率較低是因?yàn)橛嬎氵^程中對系統(tǒng)失穩(wěn)狀態(tài)的判定比較保守，不能充分利用強(qiáng)度折減迭代計算結(jié)果中蘊(yùn)含的信息，常存在已經(jīng)可以判斷系統(tǒng)處于失穩(wěn)狀態(tài)但軟件仍繼續(xù)計算的情況。本文在進(jìn)行這種數(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 某例求解不平衡力比率時變曲線

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

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

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

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

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

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

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

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

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

2.2 自適應(yīng)算法實(shí)現(xiàn)與計算流程

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

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

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

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

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

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

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

圖4 自適應(yīng)算法典型流程圖

3 自適應(yīng)算法的驗(yàn)證與分析

3.1 自適應(yīng)算法的可靠性驗(yàn)證

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

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

表1 材料參數(shù)

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

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

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

表3 計算時步對比

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

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

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

表4 不同網(wǎng)格密度下兩種算法結(jié)果差異對比

圖6 安全系數(shù)對比圖

圖7 求解時步對比圖

圖8 求解時間對比圖

4 結(jié) 論

1）通過分析FLAC3D內(nèi)置安全系數(shù)算法的求解流程和機(jī)制，發(fā)現(xiàn)判定機(jī)制過于保守、無法靈敏及時識別系統(tǒng)失穩(wěn)狀態(tài)，因而計算效率低。此外內(nèi)置算法還存在不能在大變形模式下求解安全系數(shù)、整個求解流程并不對用戶開放等不足。

2）根據(jù)失穩(wěn)狀態(tài)下求解不平衡比率時步曲線的變化規(guī)律，采用比特征響應(yīng)時步更小的分析單位實(shí)現(xiàn)對數(shù)值計算不收斂情況的超前預(yù)測和高效識別，提出了基于改進(jìn)的系統(tǒng)失穩(wěn)狀態(tài)判定機(jī)制的自適應(yīng)算法。

3）針對4 種典型邊坡對比了兩種算法求解的安全系數(shù)和時步，驗(yàn)證了自適應(yīng)算法的正確性和高效性，而且自適應(yīng)算法的求解效率隨著網(wǎng)格密度的增加而明顯增加。