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

? ? ?

摘 要：可靠度方法因其更符合邊坡巖土體的非均勻性及失穩(wěn)破壞的不確定性特征，受到科研及工程設(shè)計人員的重視。然而，目前沒有成熟且能夠用于可靠度分析的隨機有限元軟件，而跨平臺的隨機場生成和穩(wěn)定性分析增加了可靠度分析的難度，從而限制了其推廣應(yīng)用?；赑ython語言和Abaqus平臺，開發(fā)了一套能自動計算邊坡可靠度的隨機有限元算法。在給定邊坡幾何參數(shù)及土體抗剪強度參數(shù)的均值、相關(guān)距離和變異系數(shù)前提下，利用該算法可自動實現(xiàn)非平穩(wěn)隨機場的離散及邊坡失效概率的計算;該程序有效地解決了多種軟件交互使用的穩(wěn)定對接和子程序編寫等難題;與經(jīng)典的邊坡算例進行對比，計算結(jié)果驗證了該方法的可靠性。

關(guān)鍵詞：非平穩(wěn)隨機場;自動化程序;邊坡可靠度;失效概率

中圖分類號：TU431 ? 文獻標志碼：A ? 文章編號：2096-6717（2019）05-0067-06

Abstract：Based on Python and Abaqus platform， a stochastic finite element algorithm program was developed. The program can automatically discretize the non-stationary random field and calculate the failure probability when the geometric parameters of the slope and the mean， correlation distance and variation coefficient of the soil shear strength parameters are provided. The algorithm is used to calculate the benckmark slope examples， and the results verify the reliability of the proposed method. Moreover， the program can also effectively solve the complex problems of multiple software interaction and subroutine compilation. The program is then redeveloped on Abaqus platform， which is a widely applicable finite element software. It was beneficial to the popularization of the program in practical engineering applications.

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

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

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

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

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

? ? ?相對于式（1）和式（2），式（3）的認可度較高，因此，使用式（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é)點坐標，并將其導(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)隨機場。

2.2 求解過程

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

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

2）Inp文件進行初次運算。該步驟的主要目的是平衡地應(yīng)力。在初始Inp文件中，有施加土體重力的分析步驟。邊坡在初始狀態(tài)下，由于自重作用，存在與重力相平衡的應(yīng)力狀態(tài)，因此，在進行數(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ù)）輸出邊坡的失效概率。

3 算例驗證

為驗證編寫的自動化算法程序的精度，采用經(jīng)典邊坡算例。邊坡尺寸如圖3所示，坡高比為1∶2。黏聚力均值為15 kPa，標準差為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°。

邊坡采用平面應(yīng)變單元CPE4，共劃分910個單元，971個單元節(jié)點，土體失效模式采用Mohr-Coulomb 屈服準則。邊界條件為約束邊界的側(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-4 m即可滿足工程實際要求[20]，非均勻隨機場下自動化程序計算的邊坡地應(yīng)力平衡結(jié)果如圖4（a）所示，土體變形最大值所在的量級為10-5 m，滿足要求。在非均勻隨機場下邊坡的失效變形模式如圖4（b）所示，為圓弧形面，符合規(guī)律。

圖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è)的隨機變量遠多于左側(cè)數(shù)據(jù)，因此，由隨機場得到的隨機變量大多浮動在線性趨勢的右側(cè)。

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

圖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‰， 可以認為本文的計算結(jié)果準確。

4 結(jié)論?

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

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

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

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

參考文獻：

[1] HUANG J， GRIFFITHS D V. Determining an appropriate finite element size for modelling the strength of undrained random soils [J]. Computers and Geotechnics， 2015， 69： 506-513.

[2] LOW B K. Reliability analysis of rock slopes involving correlated nonnormals [J]. International Journal of Rock Mechanics and Mining Sciences， 2007， 44（6）：922-935.

[3] GRIFFITHS D V， FENTON G A. Probabilistic slope stability analysis by finite elements[J]. Journal of Geotechnical and Geoenvironmental Engineering， 2004， 130（5）： 507-518.

[4] LI D Q， CHEN Y F， LU W B， et al. Stochastic response surface method for reliability analysis of rock slopes involving correlated non-normal variables[J]. Computers and Geotechnics， 2011， 38（1）： 58-68.

[5] LOW B K， LACASSE S， NADIM F. Slope reliability analysis accounting for spatial variation[J]. Georisk： Assessment and Management of Risk for Engineered Systems and Geohazards， 2007， 1（4）： 177-189.

[6] CHO S E. Probabilistic assessment of slope stability that considers the spatial variability of soil properties[J]. Journal of Geotechnical and Geoenvironmental Engineering， 2010， 136（7）： 975-984.

[7] 宋永東.考慮巖土體強度參數(shù)空間變異性的邊坡可靠度研究[D].北京：中國科學(xué)院大學(xué)，2017.

SONG Y D. Study on slope reliability considering spatial variability of rock and soil strength parameters[D]. Beijing： University of Chinese Academy of Sciences，2017.（in Chinese）

[8] 胡金政， 張潔， 陳宏智， 等. 基于隨機場理論的強降雨條件下花崗巖殘積土邊坡的穩(wěn)定性分析[J]. 南方能源建設(shè)， 2017， 4（3）： 107-114.

HU J Z， ZHANG J， CHEN H Z， et al. Stability of completely decomposed granite slopes under intense rainfall infiltration based on the random field theory [J]. Southern Energy Construction， 2017， 4（3）： 107-114.（in Chinese）

[9] 曹少剛.土體參數(shù)空間分布對邊坡穩(wěn)定性影響研究[D].西安：西安理工大學(xué)，2017.

CAO S G. Study on the influence of spatial distribution of soil on slope stability [D]. Xian： Xian University of Technolog，2017. （in Chinese）

[10] 蔣水華. 水電工程邊坡可靠度非侵入式隨機分析方法[D]. 武漢： 武漢大學(xué)， 2014.

JIANG S H. A non-intrusive stochastic method for slope reliability in hydroelectricity engineering[D]. Wuhan： Wuhan University， 2014.（in Chinese）

[11] GRIFFTHS D V， HUANG J， FENTON G. Probabilistic slope stability analysis using RFEM with non-stationary random fields[C]//Geotechnical Safety & Risk V， 2015：704-709.

[12] 王新. 基于二維隨機場理論的非均質(zhì)邊坡穩(wěn)定分析[D]. 昆明： 昆明理工大學(xué)， 2018.

WANG X. Stability analysis of heterogeneous slope based on two-dimensional random field theory[D]. Kunming： Kunming University of Science and Technology， 2018.（in Chinese）

[13] 袁葳， 常曉林， 段寅， 等. 考慮抗剪強度參數(shù)空間變異性的邊坡穩(wěn)定性分析[J]. 中南大學(xué)學(xué)報（自然科學(xué)版）， 2016， 47（11）： 3899-3908.

YUAN W， CHANG X L， DUAN Y， et al. Stability analysis of slope considering spatial variation of shear strength parameters[J]. Journal of Central South University（Science and Technology）， 2016， 47（11）： 3899-3908.（in Chinese）

[14] PHOON K K， KULHAWY F H. Characterization of geotechnical variability[J]. Canadian Geotechnical Journal， 1999， 36（4）： 612-624.

[15] LI D Q， QI X H， PHOON K K， et al. Effect of spatially variable shear strength parameters with linearly increasing mean trend on reliability of infinite slopes [J]. Structural Safety， 2014， 49： 45-55.

[16] JIANG S H， HUANG J S. Modeling of non-stationary random field of undrained shear strength of soil for slope reliability analysis[J]. Soils and Foundations， 2018， 58（1）： 185-198.

[17] GRIFFTHS D V， ZHU D S， HUANG J S， et al. 6th Asian-Pacific[C]//Symposium on Structural Reliability and its Applications， May 28-30， 2016， Shanghai， China.

[18] DER KIUREGHIAN A， KE J B. The stochastic finite element method in structural reliability[J]. Probabilistic Engineering Mechanics， 1988， 3（2）： 83-91.

[19] HUANG J， GRIFFLTHS D V. Determining an appropriate finite element size for modelling the strength of undrained random soils [J]. Computers and Geotechnics， 2015， 69：506-513.

[20] 代汝林， 李忠芳， 王姣. 基于ABAQUS的初始地應(yīng)力平衡方法研究[J]. 重慶工商大學(xué)學(xué)報（自然科學(xué)版）， 2012， 29（9）： 76-81.

DAI R L， LI Z F， WANG J. Research on initial geo-stress balance method based on ABAQUS [J]. Journal of Chongqing Technology and Business University（Natural Science Edition）， 2012， 29（9）： 76-81.（in Chinese）

（編輯 鄧云）