趙煉恒，董瀟陽，胡世紅，左仕，潘秋景

(1.中南大學 土木工程學院，湖南 長沙 410075；2.湖南鐵院土木工程檢測有限公司，湖南 長沙 410075）

巖土體與鋼筋混凝土材料有著本質(zhì)的不同，由于土體在沉積過程及后沉積過程中經(jīng)受各種物理、化學作用，使得土體參數(shù)在空間上存在一定的變異性。工程實際中，從業(yè)者通常將影響區(qū)域內(nèi)的巖土體參數(shù)假定為常數(shù)，并將安全系數(shù)作為評價邊坡穩(wěn)定性的指標，但該種確定性分析方法并不能有效地考慮土體自身固有的變異特征。相比于傳統(tǒng)的確定性分析方法，概率分析法能夠通過計算失效概率或可靠度并以此來量化邊坡的安全程度。近年來，越來越多的學者在進行可靠度計算時將參數(shù)的空間變化特征考慮其中。GRIFFITHS等[1]分別使用隨機變量法和隨機有限元法對黏性土坡的失穩(wěn)概率進行研究；CHO[2]通過KL展開生成二維隨機場模型，將隨機場理論與極限平衡法相結(jié)合來考慮參數(shù)的空間變化特征；HUANG等[3]分析了楊氏模量的空間變異性對隧道襯砌收斂性的作用效應(yīng)；LI等[4]提出多重響應(yīng)面法，并研究了不同自相關(guān)函數(shù)對邊坡可靠度的作用效應(yīng)；GUO等[5]考慮土體的空間變化特征對路堤進行概率分析。盡管眾多研究者在巖土體空間變異性方面進行了大量有價值的研究工作，但還是存在不足：以上研究均是研究單一參數(shù)的空間變化特征，即參數(shù)的自相關(guān)性，并沒有考慮土體參數(shù)間的互相關(guān)性；另一方面研究對象過于單一，大部分研究僅考慮單一土層的層內(nèi)參數(shù)空間變異性，忽視了層間參數(shù)的變異；此外，在使用隨機有限元并結(jié)合蒙特卡洛模擬分析可靠度時，需要耗費大量的時間才能得到足夠精度的結(jié)果?；诖?，本文考慮土體參數(shù)間的互相關(guān)性和空間變異特征，利用Cholesky分解技術(shù)生成二維互相關(guān)對數(shù)正態(tài)隨機場，通過對ABAQUS的二次開發(fā)和引入失效概率變異系數(shù)實現(xiàn)對分層邊坡的高效隨機有限元分析，并結(jié)合蒙特卡洛法進行可靠度和參數(shù)敏感性分析。

1 隨機場模型

1.1 土體參數(shù)空間變異性

土體為天然材料，在經(jīng)過漫長的地質(zhì)、環(huán)境、物理、化學等綜合作用后，其性質(zhì)在空間中存在變異性[6]。為更真實地表征土體在空間中所呈現(xiàn)出的變異特征，需在分析中引入隨機場理論。隨機場應(yīng)用的前提是滿足平穩(wěn)性假設(shè)，其特征通過均值、變異系數(shù)和波動范圍來表征[7]。波動范圍是與空間變化特征的強弱程度相關(guān)的指標，該指標數(shù)值越大表明土體參數(shù)在該范圍內(nèi)的波動較為平緩，取值越小表明參數(shù)波動劇烈。受沉積作用的影響，土體在水平和豎直方向上所表征出的波動范圍存在較大差異。通常豎直方向上波動范圍要小于水平方向上的波動范圍，波動范圍的取值由自相關(guān)函數(shù)確定。本文參數(shù)隨機場假設(shè)滿足平穩(wěn)性假設(shè)，考慮自相關(guān)結(jié)構(gòu)的各向異性，選取被廣泛使用的指數(shù)型自相關(guān)函數(shù)如式(1)所示，用于計算二維平面內(nèi)任意2點間的自相關(guān)系數(shù)并組成自相關(guān)矩陣K。

式中：τx和τy分別為二維平面內(nèi)任意2網(wǎng)格中心點間的水平和豎向相對距離，τx=|xi-xj|，τy=|yiyj|，(xi,yi)，(xj,yj)分別表示第i個和第j個網(wǎng)格中心點的絕對坐標，其中i,j=1,2,…,n，n為劃分的單元數(shù)目，δh和δv分別為水平方向和豎向方向的波動范圍。

為將隨機場與有限元相結(jié)合，需對隨機場進行離散。隨機場可通過中心點[8]、最優(yōu)線性估計[9]、局部平均細分[10]、級數(shù)展開等方法生成，中心點法計算過程簡單便于編程實現(xiàn)，因此本文采用Cholesky中點法對隨機場進行離散。土體參數(shù)取值不可能為負[1]，將任一點的參數(shù)k視作隨機變量并服從對數(shù)正態(tài)分布。確定k的均值uk，變異系數(shù)COVk和波動范圍，將對數(shù)正態(tài)隨機場表示如式(2)所示：

式中：x,y是參數(shù)k的有限元網(wǎng)格中心點的坐標；(x,y)是參數(shù)k的標準高斯隨機場；ulnk和σlnk的值通過對數(shù)正態(tài)變換式(3)～(4)求出：

式中：uk和σk分別是對數(shù)正態(tài)變量k的均值和標準差，ulnk和σlnk分別為正態(tài)變量lnk的均值和標準差；COVk為對數(shù)正態(tài)變量k的變異系數(shù)。

通過Cholesky分解技術(shù)將自相關(guān)矩陣K分解為上三角矩陣ST和下三角矩陣S的乘積如下式：

對于給定的有限元網(wǎng)格和波動范圍，K矩陣中的值是確定的，即ST是確定的，此時參數(shù)k的標準高斯隨機場可表示為：

式中：X={x1,x2,…,xn}為1組服從標準正態(tài)分布的隨機變量所組成的向量。

1.2 互相關(guān)對數(shù)正態(tài)隨機場

在巖土工程中，參數(shù)間相互關(guān)聯(lián)，例如黏聚力和內(nèi)摩擦角呈現(xiàn)負相關(guān)性[5]。實際工程中，通常存在著雙參數(shù)甚至多參數(shù)隨機場的耦合效應(yīng)，本文僅考慮黏聚力c，內(nèi)摩擦角φ的相互作用，且c和φ隨機場均采用指數(shù)型自相關(guān)函數(shù)。蔣水華等[11]對生成互相關(guān)對數(shù)正態(tài)隨機場的方法已做過詳細的研究和描述，這里不做贅述。

2 基于隨機場理論的邊坡可靠度分析

邊坡穩(wěn)定性分析常用的方法有3種：極限平衡、極限分析和有限元。極限平衡和極限分析在求解安全系數(shù)時需假定潛在滑動面，如果在變異性強的土體中，采用這2種方法求解安全系數(shù)可能存在誤差。為將隨機場理論更好地應(yīng)用于邊坡工程，求得更為精確的計算結(jié)果，本文采用基于隨機場理論的隨機有限元法對邊坡的可靠度進行分析。

2.1 蒙特卡羅法

蒙特卡羅法以概率論為基礎(chǔ)，通過對定義的變量進行隨機抽樣并計算出每次抽樣樣本的功能函數(shù)。對計算結(jié)果進行統(tǒng)計和處理，即可得到失效概率及可靠度。盡管該方法需進行大量的計算工作，但數(shù)學公式相對簡單，并且不管模型的復(fù)雜性如何，該方法具有處理幾乎所有可能情況的能力[2]。

首先隨機抽取N個隨機變量樣本，再逐一計算每個樣本下的功能函數(shù)g(X)，并統(tǒng)計g(X)<0的個數(shù)Nf，對于邊坡工程，依據(jù)邊坡的破壞機理建立功能函數(shù)為：g(X)=Fs-1。計算邊坡失效概率如式(7)，失效概率與可靠度間的關(guān)系如式(8)：

可靠度分析的精度會隨樣本數(shù)量的增加而提高，但過多的樣本數(shù)量會增加不必要的工作量，為平衡計算精度與求解時間，需確定出合理的樣本數(shù)量。本文依據(jù)中心極限定理，引入失效概率Pf的變異系數(shù)COVPf來確定蒙特卡洛模擬次數(shù)[12-13]，COVPf定義為Pf的標準差與均值的比值：

蒙特卡洛模擬次數(shù)增加的同時COVPf會以的收斂速率逐漸減小，因此可通過選取合理的收斂標準來確定需要的模擬次數(shù)。本文將收斂標準設(shè)定為0.1，將100組隨機場作為一個分析組，在求解出失效概率后計算相應(yīng)的COVPf，與收斂標準進行對比確定是否需要增加分析組。

2.2 基于ABAQUS平臺邊坡可靠度分析流程

在ABAQUS中實現(xiàn)隨機有限元分析需進行二次開發(fā)，ABAQUS為用戶提供了2種二次開發(fā)路徑：一種是利用前后處理層次的Python語言；另一種是利用求解器層次的FORTRAN子程序。在隨機有限元分析時涉及到前處理過程中的單元參數(shù)賦值；在進行可靠度分析時涉及到后處理過程中的結(jié)果數(shù)據(jù)信息提取及分析；利用Python語言進行前后處理，通過ABAQUS COMMAND將INP文件提交至求解器能減少大量的工作量。此外，Python語言具備簡潔、跨平臺、免費等特點。因此，利用Python語言進行二次開發(fā)更具優(yōu)勢?；诖?，本文利用Python語言基于ABAQUS平臺實現(xiàn)邊坡的可靠度分析。通過Python腳本生成N組互相關(guān)對數(shù)正態(tài)隨機場，并結(jié)合ABAQUS對邊坡進行強度折減確定邊坡安全系數(shù)。統(tǒng)計功能函數(shù)小于0的個數(shù)Nf，并通過式(7)～(8)分別計算失效概率Pf及反算可靠度β，實現(xiàn)流程如圖1所示。該過程由5個步驟所組成：

圖1 邊坡可靠度分析流程Fig.1 Slope reliability analysis process

1)建立邊坡模型及坐標信息提取。首先在ABAQUS中建立邊坡模型并劃分網(wǎng)格，在生成的INP文件中提取各單元對應(yīng)的節(jié)點信息并生成相應(yīng)的節(jié)點信息文件和單元信息文件。

2)確定統(tǒng)計參數(shù)生成隨機場。確定強度參數(shù)的均值、變異系數(shù)及互相關(guān)系數(shù)等統(tǒng)計參數(shù)，選取指數(shù)型自相關(guān)函數(shù)和波動范圍，利用以上信息根據(jù)第2節(jié)步驟編制生成互相關(guān)對數(shù)正態(tài)隨機場的Python程序，并利用第1步生成的節(jié)點、單元信息生成N組c，φ互相關(guān)對數(shù)正態(tài)隨機場。

3)模型建立、參數(shù)輸入及獲取初始應(yīng)力場。利用Python程序?qū)崿F(xiàn)模型建立、參數(shù)賦值直至計算全過程并生成初始INP文件，再通過Python腳本更改INP文件中的強度參數(shù)批量生成N個賦有隨機場參數(shù)的INP文件，再通過ABAQUS COMMAND提交批量計算程序得到N個初始應(yīng)力場。

4)批量地應(yīng)力平衡。同步驟3，在Python腳本中添加以下語句：場變量及相應(yīng)參數(shù)值、導入第(3)步獲取的初始應(yīng)力場、修改模型關(guān)鍵字，實現(xiàn)地應(yīng)力平衡并得到最終的INP文件，再通過Python腳本更改此INP文件中的強度參數(shù)快速生成N個賦有隨機場參數(shù)信息的INP文件，通過ABAQUS COMMAND提交批量計算程序得到N結(jié)果文件。

5)計算失效概率及可靠度。在ABAQUS中通過運行腳本的方式調(diào)用已編寫好的Python腳本文件，批量打開結(jié)果數(shù)據(jù)文件并獲取到N個邊坡的安全系數(shù)值，使用EXCEL對安全系數(shù)值進行統(tǒng)計，計算出邊坡的失效概率及可靠度。

3 算例分析

本文的算例為分層土坡，坡高10 m，坡比1:1.5。幾何尺寸及網(wǎng)格劃分如圖2所示。使用摩爾-庫倫破壞準則和非關(guān)聯(lián)流動法則的理想彈塑性本構(gòu)模型。模型兩側(cè)約束水平位移，模型底部約束水平位移和豎向位移，荷載僅考慮土體自重。有限元模型劃分有2 071個4節(jié)點平面應(yīng)變單元和2 179個節(jié)點，上部土層951個單元，下部土層1 120個單元，單元尺寸為0.5 m×0.5 m。HUANG等[14]建議在隨機有限元分析中所采用的單元尺寸與空間相關(guān)距離之比應(yīng)小于0.5，指數(shù)型自相關(guān)函數(shù)空間相關(guān)距離是波動范圍的一半，因此本文單元尺寸/水平相關(guān)距離之比=0.5/1.5=0.33<0.5，滿足要求。分別對以下2種情況進行分析：確定性和隨機性。基于ABAQUS平臺，采用有限元分析不再收斂作為失穩(wěn)判據(jù)。

圖2 邊坡模型尺寸及網(wǎng)格劃分Fig.2 Slope model size and mesh generation

確定性分析時，不考慮空間變異特征，黏聚力、內(nèi)摩擦角、楊氏模量、重度、泊松比等參數(shù)采用均值如表1所示，利用ABAQUS軟件求得安全系數(shù)為1.210，利用FLAC3D進行有限差分求得安全系數(shù)為1.216，兩者基本一致。在隨機分析時，僅考慮c和φ的空間變異特征，其余參數(shù)視為定值，楊氏模量、泊松比、剪脹角等參數(shù)對邊坡穩(wěn)定性狀態(tài)分析結(jié)果的影響可忽略不計[15]。本文分析時均假設(shè)在不同的土層中土體的參數(shù)相互獨立，上下2層有著相同的自相關(guān)函數(shù)、相關(guān)系數(shù)、變異系數(shù)和波動范圍，c和φ均服從對數(shù)正態(tài)分布，均值、變異系數(shù)、相關(guān)系數(shù)及波動范圍等參數(shù)參考文獻[4,16]選取，初步設(shè)立隨機場參數(shù)取值如下：

ρ(c,φ)=0，COVc=0.3，COVφ=0.2，δh=30 m，δv=3 m，用于檢驗所設(shè)置的收斂標準是否合理，土體的參數(shù)值如表1所示。

表1 有限元模型土體參數(shù)Table 1 Soil parameters of finite element model

3.1 收斂標準的檢驗

可靠度的計算結(jié)果是否精確取決于所設(shè)置的收斂標準。為檢驗設(shè)定的收斂標準是否合理，采用圖2所示的邊坡模型，隨機場參數(shù)如前所述，邊坡的物理力學參數(shù)取值如表1所示。按第3小節(jié)所述步驟計算并統(tǒng)計出不同蒙特卡洛模擬次數(shù)下失效概率對應(yīng)的變異系數(shù)COVPf和相應(yīng)的可靠度β，統(tǒng)計結(jié)果如圖3所示。從圖3中可看出，當COVPf小于0.1時，β逐漸收斂于1.29附近，因此將0.1作為收斂標準能得到足夠精度的可靠度結(jié)果。

圖3 COVPf和β隨蒙特卡洛模擬次數(shù)的變化曲線圖Fig.3 Curve of COVPf and β with simulation times

3.2 可靠度分析

在3.1節(jié)的基礎(chǔ)上，將計算所得1 000組安全系數(shù)進行統(tǒng)計分析并繪制出安全系數(shù)散點圖和頻率分布直方圖如圖4(a)和4(b)所示。由圖4可知安全系數(shù)的分布于0.80～1.45之間，安全系數(shù)均值為1.13，最大值為1.43，最小值為0.84。1.43對應(yīng)的臨界滑面及參數(shù)隨機場如圖5所示，0.84對應(yīng)的臨界滑面及參數(shù)隨機場如圖6所示，圖中灰色的梯度變化代表著參數(shù)值大小的變化。由圖4(b)中的擬合曲線可知，安全系數(shù)服從對數(shù)正態(tài)分布。由1 000組安全系數(shù)求得邊坡的可靠度指標為1.30。

圖5 最大安全系數(shù)1.43對應(yīng)的臨界滑面及隨機場分布Fig.5 Distribution of critical slip surface and random field corresponding to the maximum safety factor of 1.43

圖6 最小安全系數(shù)0.84對應(yīng)的臨界滑面及隨機場分布Fig.6 Distribution of critical slip surface and random field corresponding to the maximum safety factor of 0.84

由圖4可知，考慮c，φ的空間變異性時，區(qū)別于確定性分析，安全系數(shù)會受到顯著的影響，同一隨機場參數(shù)下求得的安全系數(shù)既可能大于1，也可能小于1，這表明以安全系數(shù)作為邊坡的穩(wěn)定性狀態(tài)評判指標的傳統(tǒng)方法存在爭議，也表明可靠度分析的必要性。對比隨機分析計算出的安全系數(shù)均值與確定性分析計算出的安全系數(shù)，可看出確定性分析計算出的安全系數(shù)要大于隨機分析計算出的安全系數(shù)均值，這表明傳統(tǒng)的確定性分析不能科學合理的考慮土體的變異性，用確定性分析獲取的安全系數(shù)評價邊坡的穩(wěn)定性可能造成不保守估計。

圖4 安全系數(shù)分布Fig.4 Safety factor distribution

3.3 參數(shù)敏感性分析

本小節(jié)進行參數(shù)敏感性分析，邊坡的物理力學參數(shù)取值為定值如表1所示，改變2個參數(shù)間的相關(guān)系數(shù)ρ(c,φ)，黏聚力變異系數(shù)COVc，內(nèi)摩擦角變異系數(shù)COVφ，水平波動范圍δh和豎向波動范圍δv，研究隨機場參數(shù)的變化對邊坡的可靠度和安全系數(shù)均值的作用規(guī)律。

研究表明c，φ間呈負相關(guān)性，基于已有文獻[2,4]，本文ρ(c,φ)的取值范圍為[-0.7,0]，COVc取值范圍為[0.1,0.7]，COVφ取值范圍為[0.05,0.2]，δh取值范圍為[10,50]，δv取值范圍為[1,5]，計算不同參數(shù)下的可靠度及安全系數(shù)均值，結(jié)果見圖7～11。

從圖7中可看出，當ρ(c,φ)從-0.7增加至0時，相應(yīng)的β從2.457減小至1.298，uF從1.160減小至1.134，β降幅為47.17%，uF降幅為2.27%。這種趨勢表明相關(guān)系數(shù)的改變會對邊坡的可靠度和安全系數(shù)均值產(chǎn)生顯著的作用效應(yīng)，呈現(xiàn)出負相關(guān)性，即β和uF會隨ρ(c,φ)的增加而減小，且相比于uF，β對ρ(c,φ)的變化更加敏感，負相關(guān)性越強，邊坡的失效概率會減小，相應(yīng)的可靠度會增加，這與文獻[17-18]得出的規(guī)律相一致。

圖7 β和uF隨ρ(c,φ)變化曲線Fig.7 Curve of β and uF with ρ(c,φ)

圖9 β和uF隨COVφ變化曲線Fig.9 Curve of β and uF with COVφ

圖8～9分別為COVc和COVφ的變化對β和uF的作用曲線。c，φ負相關(guān)時，2幅圖中的曲線有相同的走勢，β和uF隨變異系數(shù)的增加而減小。圖8中COVc從0.1增加至0.7，相應(yīng)的β從2.256減小至1.270，uF從1.152減小至1.126，β降幅為43.70%，uF降幅為2.22%。圖8中COVφ從0.05增加至0.2，相應(yīng)的β從2.878降低至2.097，uF從1.184減小至1.150，β降幅為27.14%，uF降幅為2.96%。對比數(shù)據(jù)，在同一的增幅下，β比uF降幅更大，說明可靠度對變異系數(shù)更為敏感。計算圖8～9中起末2點斜率的絕對值，分別為1.643和5.208，可看出內(nèi)摩擦角變異系數(shù)對邊坡的可靠度影響更為顯著。

圖8 β和uF隨COVc變化曲線Fig.8 Curve of β and uF with COVc

圖11 β和uF隨δv的變化曲線Fig.11 Curve of β and uF with δv

圖10～11分別為δh和δv的變化對β和uF的作用曲線。C，φ負相關(guān)時，2幅圖的曲線呈現(xiàn)相同走勢，隨波動范圍的增大，參數(shù)的變異性減弱，β和uF呈現(xiàn)相反走勢，β隨波動范圍增大而減小，uF隨波動范圍增大而增大。δh從10增加至50，相應(yīng)的β從2.409減小至1.927，uF從1.148增加至1.152，β降幅為20.01%，uF增幅為0.34%；δv從1增加至5，相應(yīng)的β從2.366減小至1.995，uF從1.143增加至1.158，β降幅為15.6%，uF增幅為1.23%。對比分析各項數(shù)據(jù)，δh和δv均會對可靠度產(chǎn)生顯著的影響，對安全系數(shù)均值的影響較弱，即可靠度對水平和豎向波動范圍敏感。分別計算圖10～11中起末2點間斜率的絕對值，分別為0.012和0.092，可看出豎向波動范圍對邊坡的可靠度影響更為顯著。

圖10 β和uF隨δh的變化曲線Fig.10 Curve of β and uF with δh

4 結(jié)論

1)采用傳統(tǒng)的確定性分析方法評價邊坡的穩(wěn)定狀態(tài)時，僅依賴于安全系數(shù)，未考慮到巖土體的空間變異特征，缺乏合理性，可能造成不保守估計。

2)c，φ間的相關(guān)性會影響邊坡的穩(wěn)定性狀態(tài)，可靠度和安全系數(shù)均值隨相關(guān)系數(shù)的增加而減小，可靠度對相關(guān)系數(shù)更為敏感，隨著相關(guān)系數(shù)絕對值的增大，失效概率隨之減小，相應(yīng)的可靠度增加。

3)變異系數(shù)會影響邊坡的穩(wěn)定性狀態(tài)，c，φ負相關(guān)時，可靠度和安全系數(shù)均值隨變異系數(shù)的增加而減小，可靠度對變異系數(shù)更為敏感，并且內(nèi)摩擦角變異系數(shù)對可靠度的影響更為顯著。

4)波動范圍會影響邊坡的穩(wěn)定性狀態(tài)，c，φ負相關(guān)時，隨著波動范圍的逐漸增大，參數(shù)的變異性減弱，可靠度隨之減小，安全系數(shù)均值增大，可靠度對波動范圍更為敏感，并且豎向波動范圍對可靠度的影響更為顯著。