基于Python語(yǔ)言和Abaqus平臺(tái)的邊坡可靠度計(jì)算自動(dòng)化算法開(kāi)發(fā)

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

? ? ?

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

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

中圖分類(lèi)號(hào)：TU431 ? 文獻(xiàn)標(biāo)志碼：A ? 文章編號(hào)：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 ?可靠度分析法采用可靠度指標(biāo)（失效概率）代替安全系數(shù)進(jìn)行邊坡穩(wěn)定性分析[1-4]，是一種非確定性方法，更加符合邊坡巖土體的非均勻性及失穩(wěn)破壞不確定性。但由于可靠度分析法的計(jì)算相對(duì)復(fù)雜，目前在巖土工程分析和設(shè)計(jì)中的應(yīng)用尚處于研究和探索階段。

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

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

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

1 土體參數(shù)的非平穩(wěn)隨機(jī)場(chǎng)

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

2 自動(dòng)化計(jì)算流程

2.1 前期處理

非平穩(wěn)隨機(jī)場(chǎng)的形成通過(guò)以下4步來(lái)實(shí)現(xiàn)，如圖1所示。

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

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

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

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

2.2 求解過(guò)程

求解過(guò)程包括7步，如圖2所示。

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

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

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

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

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

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

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

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

3 算例驗(yàn)證

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

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

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

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

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

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

4 結(jié)論?

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

1）該程序能夠使用隨機(jī)場(chǎng)理論自動(dòng)計(jì)算邊坡的失效概率。

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

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

參考文獻(xiàn)：

[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] 宋永東.考慮巖土體強(qiáng)度參數(shù)空間變異性的邊坡可靠度研究[D].北京：中國(guó)科學(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] 胡金政， 張潔， 陳宏智， 等. 基于隨機(jī)場(chǎng)理論的強(qiáng)降雨條件下花崗巖殘積土邊坡的穩(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ù)空間分布對(duì)邊坡穩(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] 蔣水華. 水電工程邊坡可靠度非侵入式隨機(jī)分析方法[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] 王新. 基于二維隨機(jī)場(chǎng)理論的非均質(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] 袁葳， 常曉林， 段寅， 等. 考慮抗剪強(qiáng)度參數(shù)空間變異性的邊坡穩(wěn)定性分析[J]. 中南大學(xué)學(xué)報(bào)（自然科學(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é)報(bào)（自然科學(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）

（編輯 鄧云）