基于Python語言和Abaqus平臺的邊坡可靠度計算自動化算法開發

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

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

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

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

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

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

1 土體參數的非平穩隨機場

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

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

(1)

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

(2)

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

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

(3)

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

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

2 自動化計算流程

2.1 前期處理

非平穩隨機場的形成通過以下4步來實現，如圖1所示。

1)Abaqus平臺模型網格劃分。首先，給定邊坡，在Abaqus平臺下劃分網格，得到各個單元所對應的初始節點及節點坐標，并將其導出。

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

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

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

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

2.2 求解過程

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

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

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

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

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

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

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

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

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

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

3 算例驗證

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

圖3 典型邊坡幾何尺寸

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

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

圖4 邊坡變形圖

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

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

根據自動化程序得到10 000組安全系數的散點圖，如圖6所示，安全系數大多分布在1～3之間。經統計得出，安全系數的均值為1.967，標準差為0.479，最小值為0.716，最大值為3.86。根據安全系數的分布直方圖以及擬合的正態分布，可知安全系數服從正態分布。

圖6 安全系數分布圖

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

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

4 結論

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

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

2)當土體黏聚力服從對數正態分布時，邊坡的安全系數分布比較集中，服從正態分布。使用隨機場理論計算邊坡的穩定性時，邊坡失效時的滑動面為圓弧面，符合規律。

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