基于單純形-模擬退火算法的巖土力學參數反演

3

(1.長江大學 城市建設學院，湖北 荊州 434023； 2.中國一冶集團有限公司，湖北 武漢 430081； 3.東北石油大學 非常規油氣研究院，黑龍江 大慶 163318)

近年來，隨著計算機技術的飛速發展和各種數值計算方法的廣泛應用，運用數值方法研究巖土工程實際問題已成為常態。為了使數值計算結果更加準確、客觀、實用，計算參數的選取尤為重要，但因為巖土介質的復雜性，在實際中如何選取計算參數值就成為了一個難題。目前確定巖土參數的途徑主要有試驗法和反演分析法。由于室內試驗或現場原位試驗確定的相應巖土力學參數都與實際結果有著較大的偏差，因此得到的數值計算結果存在一定的誤差，這在一定程度上也限制了數值方法在實際工程中的應用。而基于實測信息的反演分析方法則為巖土力學參數的獲取提供了新的途徑，且表現出了獨特的優勢[1-2]。由于巖體結構在開挖或變形過程中位移是最容易測定的，利用現場位移量測信息反演巖土力學參數的位移反分析法在巖土工程中的應用較為普遍。

位移反分析法是將參數反演問題轉化為目標函數尋優問題，但是目標函數是一個很復雜的非線性多峰函數，傳統優化方法存在著優化結果依賴于初值選擇、容易陷入局部極值等不足，采用全局優化算法成為較理想的選擇。一些學者將遺傳算法、粒子群算法等智能優化算法應用到巖土參數的反演計算中，取得了較為理想的結果，但是存在著計算量大、收斂速度慢以及計算時間長等不足，而且這些方法也并沒有真正將有限元程序嵌入到相應優化算法中[3-7]。一些學者將自己開發的有限元程序嵌入到優化算法中，編寫了優化反分析程序，但自編的有限元程序其功能不如大型商業有限元軟件優越，所以也在一定程度上限制了反演分析的功能，對復雜本構模型以及多場耦合參數的反演往往無法實現[2，8-9]。還有一些學者將節理巖體視為橫觀各向同性介質對巖土參數進行反演，忽略了節理面的影響，反演結果與實際存在著較大誤差[10-11]。

本文將單純形算法和模擬退火算法相結合組成混合算法，并基于該算法編寫優化反分析程序，將ABAQUS作為一個單獨模塊嵌入到優化算法程序中，實現MATLAB和ABAQUS的聯合反演，加快優化算法的收斂速度，減少反演計算的迭代次數。基于ABAQUS自帶的節理材料模型，編寫了節理巖體各向異性彈塑性本構模型子程序USDFLD，并應用到實際工程數值模型中，提高了反演參數精度。

1 反演模型

位移反分析是利用現場實測的變形作為觀測信息反演巖土力學參數。待反演參數可用向量x表示，即x=[x1,x2,…,xm]T，其中m表示待反演參數總個數。以計算值與測點實測值之間的誤差建立目標函數如下：

(1)

在參數優化反演過程中，尋找一組材料參數{x*}，構造目標函數的均值為新的目標函數，使新目標函數得到最小值時的反演參數即視為現場巖體材料參數，此時目標函數值為

(2)

巖土工程反演問題通常較為復雜，為保證解的穩定性和唯一性需要加一些約束條件作為約束優化問題來處理，則有限元優化反演模型可表示為

(3)

式中，Mi(x)為第i個等式約束函數；Nj(x)為第j個不等式約束。

采用精確罰函數法將上述約束優化問題變成無約束優化問題，即：

(4)

2 反演步驟與程序實現

由于巖土工程反演模型較為復雜，無法準確寫出其解析表達式，用計算目標函數的導數來計算出目標函數極小值就無法實現，因此本文采用單純形-模擬退火優化算法來解決這個問題。

2.1 單純形-模擬退火優化算法

單純形優化算法是先求出一個基本可行解，然后逐步改善可行解，使目標函數值逐步增大(或減小)，直到目標函數達到極值時，就可得到最優解。模擬退火優化算法是基于Monte-Carlo迭代求解策略的一種隨機尋優算法，是源于組合優化與物理中固體物質退火過程之間的相似性。模擬退火優化算法是在某一較高初溫下，伴隨溫度參數的不斷下降，結合概率突跳特性在解空間中隨機尋找目標函數的全局最優解，即在局部最優解中能概率性地跳出并最終趨于全局最優解。

將兩種算法結合形成一種混合算法，即單純形-模擬退火優化算法，該算法的主要計算步驟如下。

(1) 給出待反演參數的初值x0，初解為G(x0,u0)，設定算法控制求解精度。

(4) 利用Metropolis準則判別函數收斂性。

(5) 若優化后目標函數值大于上步的迭代函數值，則將上步的反演參數替代該反演參數，返回步驟(2)。

(6) 設置線性溫度降低函數，初溫設為某一固定值。

2.2 程序研制

基于本文提出的反演思路和方法，在ABAQUS的基礎上，采用MATLAB語言為平臺并結合單純形-模擬退火算法編制了優化反演子程序SMSAInverse.m。優化反分析程序的具體實施步驟如下。

(1) 用MATLAB語言編寫反演主程序，給出待反演參數的初值和測點的實測位移值，并用MATLAB語言調用ABAQUS的命令流文件進行計算。

(2) 用MATLAB語言編寫目標函數子程序，調用ABAQUS的計算結果dat文件，獲取測點的位移計算值，以此計算目標函數。

(3) 編寫單純形-模擬退火算法，在算法中定義相應的邊界條件并初始化其參數。

(4) 先用單純形算法優化出待反演參數，再用strrep命令修正命令流文件中的待反演參數，當目標函數值滿足精度要求時，以反演后的參數作為模擬退火算法的初值進行優化迭代，并同樣通過修改命令流文件更新反演參數，直到滿足迭代次數或計算精度要求為止，停止迭代。

優化反分析程序如圖1所示。

3 材料本構模型

目前理論上多將節理巖體視為橫觀各向同性介質[12]，其應力-應變關系包含5個彈性參數，但由于巖體中存在節理面，節理巖體和巖石的力學特征存在較大的差別，因此節理巖體表現出各向異性特征。由橫觀各向同性理論[13]可知，局部坐標系中橫觀各向同性增量的彈性應力-應變關系表達式如下：

圖1 巖土力學參數優化反演程序框圖Fig.1 Flow chart of inversion analysis for mechanical parameters of geotechnical materials

[Δσ′]=[K′][Δε′e]

(5)

式中,[K′]為局部對稱剛度矩陣，寫入對應的應力應變對號，則可表示為

(6)

巖石應力應變的關系不僅與這5個彈性參數有關，還與節理面傾角有關[14]，由各向異性理論可知，軸向應力σy與應變(εx,εy,εz)的關系如下：

εx=a12σy,εy=a22σy,εz=a23σy

(7)

式中,a12,a22,a23是由5個彈性參數和節理面傾角θ決定的函數關系式，表達式如下：

(8)

(9)

(10)

當5個彈性參數確定時，即可確定任一節理面傾角的彈性參數，即：

(11)

基于上述理論分析，本文以ABAQUS自帶的節理材料模型為基礎，結合節理巖體彈性參數的各向異性，建立了節理巖體各向異性彈塑性本構關系，并用FORTRAN語言開發了材料本構模型子程序USDFLD，將其應用到本文數值模擬的模型當中。

4 工程應用

4.1 工程概況

武漢花山大道寶蓋山隧道工程采用小凈距分離式隧道，隧道樁號為K0+960～K1+542，隧道全長582 m，穿越寶蓋山與缽孟峰山谷，隧道線型為直線，具體概況詳見武漢花山大道工程報告[15]。本文選取寶蓋山隧道左線樁號K1+360為研究對象，隧道圍巖屬于Ⅲ級圍巖，采用上下臺階開挖法，該隧道支護方式如圖2所示，采用復合式襯砌設計，并釘入注漿錨桿，混凝土強度等級為C25，系統錨桿采用Φ22/Φ25的中空注漿錨桿，二次模筑襯砌厚度為450 mm。隧道開挖的主要施工步驟如圖3所示。

圖2 圍巖支護結構布置Fig.2 Rock support structure layout

圖3 隧道開挖步驟Fig.3 Excavation steps

4.2 數值模型

由寶蓋山隧道左線K1+360處地質資料可知，隧洞方向長度約為243 m，隧道高度方向為10.157 m，跨度方向長度為13.75 m。以此建立隧洞斷面二維模型，模型的豎直方向高度為99.5 m，水平方向長度為160 m。襯砌厚度為450 mm，單元類型為CPE4，且錨桿采用半徑為25mm的B21桿單元，襯砌和錨桿均采用節點公用單元。模型兩側施加水平方向位移約束，底邊施加豎直方向固定約束，單元類型為CPE4(四節點平面應變單元)，有限元數值模型及測點分布如圖4所示。

圖4 有限元網格模型及測點布置Fig.4 Measuring points and finite element mesh model

如圖4所示，參與反演的測點為A、B、C、D、E五個測點，其中A、B、C、D為水平方向測點，E為豎直方向測點，測點對應的有限元模型中單元節點分別為991、699、28、30、48，有部分測點與單元節點無法對應，采用插入法取測點位移值，則該5個測點實測位移值如表2所示。

表2 測點位移實測值Tab.2 Measured displacements of selected points mm

筆者用正交設計試驗方法分析了巖土力學參數對于圍巖變形的敏感性，最后選取了節理巖體凝聚力、內摩擦角和節理面凝聚力作為待反演參數，具體過程在此不再贅述。將USDFLD材料本構子程序應用到隧道數值模型中，且模型中參與反演的材料參數的初值及取值范圍如表3所示。

表3 待反演參數設置Tab.3 Initial setting of undetermined parameters for inversion model

隧洞圍巖參數中彈性模量為0.244 GPa，泊松比取0.18，密度為2 500 kg/m3，節理巖體凝聚力為3.2 MPa，節理巖體內摩擦角為24.001°，節理面凝聚力為0.3 MPa，節理面內摩擦角為10°；襯砌的彈性模量為9 GPa，泊松比為0.25，密度為2 500 kg/m3；錨桿彈性模量為180 GPa，泊松比為0.22，密度為7 800 kg/m3。在數值模擬過程中，復制開挖土體并改變其材料參數，在開挖原始土體后添加彈性模量較小的土體，以此達到彈性模量衰減的目的。

數值模擬開挖的過程如下：創建二維數值模型，設定材料參數和荷載邊界條件，進行地應力平衡；復制①、②兩部分土體，且每部分復制兩份，首先開挖原始土體①部分并添加模量較小的土體以達到模量衰減的目的，然后施加上部襯砌；同樣的方法開挖②部分土體，然后施加下部襯砌，開挖結束。

4.3 反演結果

對上述二維數值模型進行有限元計算，將實際測點對應的模型測點的位移值作為理論計算值，依據該計算值可優化反演出相關待反演參數，并通過比較實測點的計算位移值與實測位移值對本文方法的合理性和精度做出判斷。

利用單純形-模擬退火算法進行優化迭代，首先以待反演參數的反演區間上下限的平均值作為初始值，由單純形算法進行優化迭代，得出反演參數值。然后以單純形反演值作為初始值，進行模擬退火算法優化計算，最終得到待反演參數的反演結果，如表4所示。最終優化程序迭代了65次，調用有限元軟件ABAQUS計算177次。反演求解迭代次數曲線如圖5所示。

表4 參數反演結果Tab.4 Back analysis results for mechanical parameters

圖5 目標函數值變化趨勢Fig.5 Iteration process curve of back analysis

測點位移的反演結果和與實測值比較的結果如表5所示。可知，實測值與反演值誤差較小，30測點的Ux位移值誤差較大，但整體位移值較小，其它測點的位移相對誤差均控制在5%以內。

表5 反演結果與實測值的比較Tab.5 Comparison of measured results and numerical results for measuring point displacement

由上可知，由單純形局部優化算法與模擬退火全局優化算法相結合的混合算法，解決了單純形算法易陷入局部極小值和模擬退火算法搜索效率低的問題，提高了搜索效率和反演參數精度，為確定地下工程圍巖力學參數提供了有效的途徑。將自行開發的節理本構模型應用到花山大道寶蓋山隧道工程當中，能夠很好地反映隧洞圍巖變形及受力特性，為指導實際工程起到了一定的預見作用。

5 結 論

(1) 將單純形優化算法和模擬退火算法相結合，構造一種混合算法并結合合適的精確罰函數，建立了基于計算值和測點實測值之間的目標函數，并將有限元軟件ABAQUS作為單獨模塊嵌入到單純形-模擬退火算法中，實現ABAQUS與MATLAB聯合反演，可保證該算法不依賴于初始值就能搜索到全局最優解，且可加快搜索速度。

(2) 建立花山大道寶蓋山隧道開挖的二維數值模型，將開發的節理巖體各向異性本構子程序應用到數值模型中，并嵌入到編寫的優化反演程序中，對節理巖體凝聚力、節理巖體內摩擦角和節理面凝聚力三個力學參數進行反演，得到了合理的參數值，說明本文提出的方法可用于實際工程，為確定圍巖力學參數提供了一套有效的方法。