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

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

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

可靠度分析法采用可靠度指標(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)

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

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

(1)

式中：su、σv′和OCR分別表示土體的黏聚力、有效應(yīng)力及固結(jié)比。Jiang等[16]考慮土體抗剪強(qiáng)度參數(shù)的隨機(jī)波動(dòng)量服從對(duì)數(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)隨機(jī)場(chǎng)，研究邊坡的失效概率。

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

(3)

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

相對(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)。

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

2.2 求解過(guò)程

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

1)Python形成初始Inp文件。將邊坡的幾何參數(shù)、材料信息及隨機(jī)場(chǎng)數(shù)據(jù)寫(xiě)入Inp文件，該文件稱為初始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ù))輸出邊坡的失效概率。

圖2 自動(dòng)化程序的求解過(guò)程Fig.2 The solving part of the automation

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°。

圖3 典型邊坡幾何尺寸

邊坡采用平面應(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-4m即可滿足工程實(shí)際要求[20]，非均勻隨機(jī)場(chǎng)下自動(dòng)化程序計(jì)算的邊坡地應(yīng)力平衡結(jié)果如圖4(a)所示，土體變形最大值所在的量級(jí)為10-5m，滿足要求。在非均勻隨機(jī)場(chǎng)下邊坡的失效變形模式如圖4(b)所示，為圓弧形面，符合規(guī)律。

圖4 邊坡變形圖

圖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è)。

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

根據(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)分布。

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

圖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)確。

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

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í)際失效概率。