任斌斌，蘇立君，張崇磊，謝奇峻

(1.中國科學(xué)院 山地災(zāi)害與地表過程重點實驗室；中國科學(xué)院、水利部 成都山地災(zāi)害與環(huán)境研究所，成都 610041；2．中國科學(xué)院 青藏高原地球科學(xué)卓越創(chuàng)新中心，北京 100101；3．中國科學(xué)院大學(xué) 人工智能學(xué)院，北京 100049)

可靠度分析法采用可靠度指標(biāo)(失效概率)代替安全系數(shù)進(jìn)行邊坡穩(wěn)定性分析[1-4]，是一種非確定性方法，更加符合邊坡巖土體的非均勻性及失穩(wěn)破壞不確定性。但由于可靠度分析法的計算相對復(fù)雜，目前在巖土工程分析和設(shè)計中的應(yīng)用尚處于研究和探索階段。

目前的可靠度方法研究和應(yīng)用中，蒙特卡洛模擬[5-6]應(yīng)用較為廣泛。已有學(xué)者基于隨機場理論，使用蒙特卡洛模擬進(jìn)行邊坡穩(wěn)定性研究。宋永東[7]運用Matlab離散隨機場，利用Excel作為銜接手段，將離散后的土體強度參數(shù)導(dǎo)入有限差分軟件Flac3D，計算邊坡的穩(wěn)定性；胡金政等[8]用Flac3D建模，利用Fish語言將離散場與網(wǎng)格單元一一對應(yīng)，反復(fù)N次，進(jìn)行邊坡穩(wěn)定性計算，最后使用Matlab讀取計算結(jié)果；曹少剛[9]使用Matlab 編寫程序，得到能夠表現(xiàn)土體參數(shù)空間變異性的一系列隨機變量，然后使用Flac3D計算邊坡的安全系數(shù)；蔣水華[10]利用有限元軟件Abaqus和GeoStudio編寫接口程序，計算邊坡的可靠度指標(biāo)；Griffths等[11]使用Fortran語言編寫耦合隨機場理論與邊坡可靠度分析軟件；王新[12]使用Matlab獲取離散隨機場，然后與Abaqus模型相結(jié)合進(jìn)行邊坡的可靠度計算；袁葳等[13]以隨機場理論為基礎(chǔ)，使用Abaqus提供的用戶子程序接口編寫隨機有限元程序，使用Python腳本進(jìn)行后期處理。

上述研究在“隨機有限元程序”應(yīng)用方面取得了一定的進(jìn)展，但仍然存在不足。首先，使用Flac3D與Matlab計算時，對不同的土體參數(shù)進(jìn)行敏感性或影響程度分析時，需要在兩者之間進(jìn)行數(shù)以萬次反復(fù)轉(zhuǎn)換，計算量較大且耗時較長；其次，調(diào)用Abaqus內(nèi)核進(jìn)行批處理時，并沒有涉及地應(yīng)力迭代過程，這將導(dǎo)致計算結(jié)果存在一定的誤差；最后，使用GeoStudio、Abaqus與Matlab相結(jié)合時，有限元軟件與編寫程序所使用的語言不一致，會降低原程序的計算效率。

本文利用Abaqus腳本建模使用的Python語言編寫程序，將有限元建模、隨機場賦值和強度折減計算有機結(jié)合起來，進(jìn)行批量自動化運算，實現(xiàn)高效精確的邊坡可靠度分析。

1 土體參數(shù)的非平穩(wěn)隨機場

Phoon等[14]指出，土體的強度是各向異性的，可以從兩個方向?qū)ν馏w的強度進(jìn)行分析。對正常固結(jié)土來說，在垂直方向上，從地表開始，土體的強度隨深度加深逐漸增加；在水平方向上，土體的強度是與深度無關(guān)的隨機波動量。由于非平穩(wěn)隨機場能夠合理地模擬土體強度參數(shù)隨埋深逐漸增加的特征，因此，許多學(xué)者進(jìn)行了相關(guān)研究。Li等[15]通過式(1)研究了土體黏聚力隨有效應(yīng)力及固結(jié)比的變化關(guān)系。

su/σv′=(0.23±0.04)OCR0.8

(1)

式中：su、σv′和OCR分別表示土體的黏聚力、有效應(yīng)力及固結(jié)比。Jiang等[16]考慮土體抗剪強度參數(shù)的隨機波動量服從對數(shù)正態(tài)分布，將其隨深度的變化關(guān)系表示為式(2)。

(2)

式中:su(x,z)為某位置處的土體黏聚力；t為比例因子；w(x,z)服從正態(tài)分布。Griffiths等[11,17]、Der Kiureghian等[18]等利用式(3)建立非平穩(wěn)隨機場，研究邊坡的失效概率。

cz=c0(μcu0+ρz)/μcu0

(3)

式中:cz為某深度處的土體黏聚力；c0為由非平穩(wěn)隨機場得到的土體黏聚力值；μcu0為地表處黏聚力均值；ρ為比例因子；z為特定的土體深度。

相對于式(1)和式(2)，式(3)的認(rèn)可度較高，因此，使用式(3)將平穩(wěn)隨機場轉(zhuǎn)化為非平穩(wěn)隨機場，研究土體參數(shù)的空間變異性對邊坡可靠度影響。

2 自動化計算流程

2.1 前期處理

非平穩(wěn)隨機場的形成通過以下4步來實現(xiàn)，如圖1所示。

1)Abaqus平臺模型網(wǎng)格劃分。首先，給定邊坡，在Abaqus平臺下劃分網(wǎng)格，得到各個單元所對應(yīng)的初始節(jié)點及節(jié)點坐標(biāo)，并將其導(dǎo)出。

2)Python讀取數(shù)據(jù)并對單元排序。將步驟1)導(dǎo)出的單元重新排序，保持節(jié)點序號不變，目的是使離散后的隨機場變量能夠批量賦值給對應(yīng)的邊坡單元。

3)平穩(wěn)隨機場。主要包括隨機場的離散和有限元的結(jié)合，利用中心點法離散隨機場，得到一系列隨機變量，然后按照邊坡的實際空間位置，將隨機變量映射到步驟2)得到的有限元單元中。

4)非平穩(wěn)隨機場。非平穩(wěn)隨機場與土體參數(shù)實際分布比較接近，將步驟3)得到的平穩(wěn)隨機場轉(zhuǎn)化為非平穩(wěn)隨機場。

圖1 自動化程序的前處理部分Fig.1 The preprocessing part of the automation

2.2 求解過程

求解過程包括7步，如圖2所示。

1)Python形成初始Inp文件。將邊坡的幾何參數(shù)、材料信息及隨機場數(shù)據(jù)寫入Inp文件，該文件稱為初始Inp文件。

2)Inp文件進(jìn)行初次運算。該步驟的主要目的是平衡地應(yīng)力。在初始Inp文件中，有施加土體重力的分析步驟。邊坡在初始狀態(tài)下，由于自重作用，存在與重力相平衡的應(yīng)力狀態(tài)，因此，在進(jìn)行數(shù)值模擬時，需要在邊坡開始運算之前建立相應(yīng)的應(yīng)力場。

3)提取初始應(yīng)力生成Rpt文件。經(jīng)過步驟2)的初始運算，得到一系列包含各個單元應(yīng)力的Job文件，然后使用Python編寫的腳本文件，提取各個單元的內(nèi)力，并生成包含各個應(yīng)力提取代碼的Rpt文件。

4)地應(yīng)力平衡的Csv文件。在Abaqus平臺上運行步驟3)得到Rpt文件，得到與每種情況相對應(yīng)的Csv文件，以便平衡地應(yīng)力。

5)地應(yīng)力平衡。模型的地應(yīng)力平衡結(jié)果滿足要求后，程序自動調(diào)用提前編寫的命令讀取Csv文件。

6)得到最終的Inp文件。在初始Inp文件中加入強度折減法的分析步，得到最終的Inp文件。

7)Abaqus強度折減法運算。調(diào)用Abaqus求解器得到最終包含邊坡變形、應(yīng)力和場變量等信息的Job文件。

8)計算邊坡的失效概率。根據(jù)Pf=Nfs<1/N(Pf表示邊坡的失效概率;Nfs<1表示安全系數(shù)小于1的數(shù)量;N表示總的計算次數(shù))輸出邊坡的失效概率。

圖2 自動化程序的求解過程Fig.2 The solving part of the automation

3 算例驗證

為驗證編寫的自動化算法程序的精度，采用經(jīng)典邊坡算例。邊坡尺寸如圖3所示，坡高比為1∶2。黏聚力均值為15 kPa，標(biāo)準(zhǔn)差為4 kPa，水平相關(guān)距離為38 m，豎向相關(guān)距離為3.8 m。在隨機場理論中，相關(guān)距離是指土體中任意兩點性質(zhì)不相關(guān)的最小距離，是土體的天然特性。土體天然密度ρ為2 000 kg/m3，變形模量為10 MPa，泊松比v=0.3[17]。為簡化計算，只考慮黏聚力生成的隨機場，內(nèi)摩擦角為0°。

圖3 典型邊坡幾何尺寸

邊坡采用平面應(yīng)變單元CPE4，共劃分910個單元，971個單元節(jié)點，土體失效模式采用Mohr-Coulomb 屈服準(zhǔn)則。邊界條件為約束邊界的側(cè)向位移及底部的水平及豎向位移。Der Kiureghian等[18]和Huang等[19]指出單元尺寸與相關(guān)距離之比應(yīng)小于0.25。單元水平長度為2 m，高度為0.5 m，其中，單元水平長度/水平相關(guān)距離=2/38=0.05<0.25，單元高度/豎向相關(guān)距離=0.5/3.8=0.13<0.25，單元尺寸符合要求。

地應(yīng)力平衡是巖土工程數(shù)值模擬過程中的重要內(nèi)容，根據(jù)一般巖土工程對地應(yīng)力平衡的要求，土體變形小于10-4m即可滿足工程實際要求[20]，非均勻隨機場下自動化程序計算的邊坡地應(yīng)力平衡結(jié)果如圖4(a)所示，土體變形最大值所在的量級為10-5m，滿足要求。在非均勻隨機場下邊坡的失效變形模式如圖4(b)所示，為圓弧形面，符合規(guī)律。

圖4 邊坡變形圖

圖5為黏聚力在非平穩(wěn)隨機場下的變化規(guī)律及均值線性趨勢圖。由于圖5可知，本程序計算的結(jié)果與Jiang等[17]的非平穩(wěn)隨機場下土體的黏聚力值均在其均值線性趨勢線的右邊。這是因為將平穩(wěn)隨機場轉(zhuǎn)化為非平穩(wěn)隨機場，并沒有將土體強度參數(shù)隨深度增加的趨勢分量與波動分量分開，波動分量并不明顯。而假定的土體強度參數(shù)分布為對數(shù)正態(tài)分布，由對數(shù)正態(tài)函數(shù)的頻率分布圖可知，均值右側(cè)的隨機變量遠(yuǎn)多于左側(cè)數(shù)據(jù)，因此，由隨機場得到的隨機變量大多浮動在線性趨勢的右側(cè)。

圖5 黏聚力隨深度的變化關(guān)系Fig.5 The relationship between cohesion and

根據(jù)自動化程序得到10 000組安全系數(shù)的散點圖，如圖6所示，安全系數(shù)大多分布在1～3之間。經(jīng)統(tǒng)計得出，安全系數(shù)的均值為1.967，標(biāo)準(zhǔn)差為0.479，最小值為0.716，最大值為3.86。根據(jù)安全系數(shù)的分布直方圖以及擬合的正態(tài)分布，可知安全系數(shù)服從正態(tài)分布。

圖6 安全系數(shù)分布圖

圖7為不同模擬次數(shù)下對應(yīng)的邊坡失效概率，可知，當(dāng)模擬次數(shù)在1 000～10 000之間時，由自動化程序得到的該邊坡失效概率曲線趨于穩(wěn)定，此時，所對應(yīng)的邊坡的失效概率為7.2‰。Jiang等[17]采用同樣土體參數(shù)計算得到邊坡失效概率為5.28‰，誤差來源主要為失效概率計算方法的偏差，Jiang等[17]采用子集模擬法計算邊坡的失效概率，而本文使用蒙特卡洛模擬法。子集模擬是一種求解失效概率的近似方法，所得到的結(jié)果是近似結(jié)果，而蒙特卡洛模擬是檢驗其他方法的依據(jù)，并且兩者結(jié)果僅相差1.92‰， 可以認(rèn)為本文的計算結(jié)果準(zhǔn)確。

圖7 邊坡的失效概率趨勢圖Fig.7 Trend Chart of Slope Failure

4 結(jié)論

基于Abaqus開放接口，使用Python語言進(jìn)行二次開發(fā)，編寫隨機有限元腳本文件，用以計算邊坡的可靠度。當(dāng)邊坡的幾何形狀確定后，只需運行幾個特定的腳本文件，便可利用該程序求解基于隨機場理論的邊坡可靠度，使用方便。主要結(jié)論如下：

1)該程序能夠使用隨機場理論自動計算邊坡的失效概率。

2)當(dāng)土體黏聚力服從對數(shù)正態(tài)分布時，邊坡的安全系數(shù)分布比較集中，服從正態(tài)分布。使用隨機場理論計算邊坡的穩(wěn)定性時，邊坡失效時的滑動面為圓弧面，符合規(guī)律。

3)計算邊坡的失效概率時，蒙特卡洛模擬次數(shù)較大時，計算得到的失效概率逐漸趨于穩(wěn)定，在計算未知邊坡的失效概率時，不必過多設(shè)置模擬次數(shù)，以免耗時過長，可以近似認(rèn)為邊坡失效概率曲線穩(wěn)定時，對應(yīng)的失效概率為邊坡的實際失效概率。