3

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

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

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

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

1 反演模型

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

(1)

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

(2)

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

(3)

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

采用精確罰函數(shù)法將上述約束優(yōu)化問題變成無約束優(yōu)化問題，即：

(4)

2 反演步驟與程序實現(xiàn)

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

2.1 單純形-模擬退火優(yōu)化算法

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

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

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

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

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

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

2.2 程序研制

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

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

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

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

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

優(yōu)化反分析程序如圖1所示。

3 材料本構模型

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

圖1 巖土力學參數(shù)優(yōu)化反演程序框圖Fig.1 Flow chart of inversion analysis for mechanical parameters of geotechnical materials

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

(5)

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

(6)

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

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

(7)

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

(8)

(9)

(10)

當5個彈性參數(shù)確定時，即可確定任一節(jié)理面傾角的彈性參數(shù)，即：

(11)

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

4 工程應用

4.1 工程概況

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

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

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

4.2 數(shù)值模型

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

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

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

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

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

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

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

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

4.3 反演結果

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

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

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

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

測點位移的反演結果和與實測值比較的結果如表5所示?？芍瑢崪y值與反演值誤差較小，30測點的Ux位移值誤差較大，但整體位移值較小，其它測點的位移相對誤差均控制在5%以內(nèi)。

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

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

5 結 論

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

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