3

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

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

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

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

1 反演模型

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

(1)

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

(2)

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

(3)

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

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

(4)

2 反演步驟與程序?qū)崿F(xiàn)

由于巖土工程反演模型較為復(fù)雜，無(wú)法準(zhǔn)確寫(xiě)出其解析表達(dá)式，用計(jì)算目標(biāo)函數(shù)的導(dǎo)數(shù)來(lái)計(jì)算出目標(biāo)函數(shù)極小值就無(wú)法實(shí)現(xiàn)，因此本文采用單純形-模擬退火優(yōu)化算法來(lái)解決這個(gè)問(wèn)題。

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

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

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

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

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

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

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

2.2 程序研制

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

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

(2) 用MATLAB語(yǔ)言編寫(xiě)目標(biāo)函數(shù)子程序，調(diào)用ABAQUS的計(jì)算結(jié)果dat文件，獲取測(cè)點(diǎn)的位移計(jì)算值，以此計(jì)算目標(biāo)函數(shù)。

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

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

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

3 材料本構(gòu)模型

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

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

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

(5)

式中,[K′]為局部對(duì)稱剛度矩陣，寫(xiě)入對(duì)應(yīng)的應(yīng)力應(yīng)變對(duì)號(hào)，則可表示為

(6)

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

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

(7)

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

(8)

(9)

(10)

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

(11)

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

4 工程應(yīng)用

4.1 工程概況

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

圖2 圍巖支護(hù)結(jié)構(gòu)布置Fig.2 Rock support structure layout

圖3 隧道開(kāi)挖步驟Fig.3 Excavation steps

4.2 數(shù)值模型

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

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

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

表2 測(cè)點(diǎn)位移實(shí)測(cè)值Tab.2 Measured displacements of selected points mm

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

表3 待反演參數(shù)設(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ù)值模擬過(guò)程中，復(fù)制開(kāi)挖土體并改變其材料參數(shù)，在開(kāi)挖原始土體后添加彈性模量較小的土體，以此達(dá)到彈性模量衰減的目的。

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

4.3 反演結(jié)果

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

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

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

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

測(cè)點(diǎn)位移的反演結(jié)果和與實(shí)測(cè)值比較的結(jié)果如表5所示?？芍?，實(shí)測(cè)值與反演值誤差較小，30測(cè)點(diǎn)的Ux位移值誤差較大，但整體位移值較小，其它測(cè)點(diǎn)的位移相對(duì)誤差均控制在5%以內(nèi)。

表5 反演結(jié)果與實(shí)測(cè)值的比較Tab.5 Comparison of measured results and numerical results for measuring point displacement

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

5 結(jié) 論

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

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