基于Hoek-Brown準(zhǔn)則的巖體彈塑性損傷模型及其應(yīng)力回映算法研究

許夢(mèng)飛，姜諳男，史洪濤，李興盛

(1.大連海事大學(xué)道路與橋梁工程研究所，遼寧，大連 116026；2.中鐵建大橋工程局集團(tuán)第一工程有限公司，遼寧，大連 116026；3.中鐵一局集團(tuán)第二工程有限公司，河北，唐山063000)

在巖體工程開(kāi)挖過(guò)程中，巖體介質(zhì)內(nèi)部隨著裂隙的發(fā)育，貫通產(chǎn)生損傷，導(dǎo)致其強(qiáng)度和剛度發(fā)生劣化，并伴隨有塑性流動(dòng)變形。相當(dāng)部分巖體的力學(xué)行為不但受巖體結(jié)構(gòu)、巖塊強(qiáng)度、應(yīng)力狀態(tài)等因素影響，并且具有強(qiáng)度非線性和應(yīng)變軟化特征。同時(shí)地下水在巖體內(nèi)部將產(chǎn)生孔隙水壓力，會(huì)影響巖體的力學(xué)性質(zhì)和破壞模式。因此，為了更好地反映巖體的力學(xué)行為規(guī)律，建立合理的巖體彈塑性損傷本構(gòu)模型具有重要的意義。

近年來(lái)，國(guó)內(nèi)外學(xué)者根據(jù)彈塑性力學(xué)、連續(xù)損傷力學(xué)理論對(duì)巖石彈塑性損傷模型進(jìn)行了研究，包括引入合理的屈服函數(shù)、損傷變量演化方程和建立穩(wěn)健的數(shù)值求解算法。Luccioni和Armero等[1-2]基于熱力學(xué)框架闡述了塑性與損傷的耦合機(jī)理，并給出了耦合模型的數(shù)值積分算法。Shao、Salari、王軍祥和袁小平等[3-6]建立了基于Druker-Prager (DP)的巖石彈塑性損傷模型，較好地反映了損傷對(duì)粘聚力和剛度的弱化作用。王永亮等[7-8]推導(dǎo)了DP準(zhǔn)則下非均勻巖石損傷本構(gòu)模型，并與流固耦合理論相結(jié)合，開(kāi)發(fā)了損傷巖石的滲流求解程序。劉楊等[9]建立了基于Mohr-Coulomb(MC)準(zhǔn)則的彈塑性損傷本構(gòu)模型，并推導(dǎo)了模型在主應(yīng)力空間中的應(yīng)力回映算法，解決了MC準(zhǔn)則在數(shù)值計(jì)算中的奇異點(diǎn)問(wèn)題。賈善坡等[10]針對(duì)泥巖的力學(xué)試驗(yàn)特性，以等效塑性應(yīng)變?yōu)閾p傷內(nèi)變量，建立了基于MC準(zhǔn)則的塑性損傷模型和蠕變損傷模型，并在ABAQUS平臺(tái)上實(shí)現(xiàn)了模型的數(shù)值求解。楊強(qiáng)等[11]利用細(xì)觀力學(xué)方法建立了巖土材料的彈塑性損傷模型，該模型能夠模擬巖土破壞中的局部化問(wèn)題。杜修力等[12]結(jié)合非線性統(tǒng)一強(qiáng)度模型和考慮圍壓作用的損傷演化方程建立了一種巖石三維彈塑性損傷模型。

已有巖石彈塑性損傷模型大多采用MC或DP等線性強(qiáng)度準(zhǔn)則。同DP和MC準(zhǔn)則相比，廣義Hoek-Brown(HB)屈服準(zhǔn)則[13]更能夠反映巖體的非線性特征及結(jié)構(gòu)面、開(kāi)挖擾動(dòng)等因素對(duì)巖體強(qiáng)度的影響，因此被廣泛應(yīng)用于巖體工程穩(wěn)定性評(píng)價(jià)當(dāng)中。朱合華等[14]闡述了HB強(qiáng)度準(zhǔn)則的研究進(jìn)展，并提出一種廣義三維HB強(qiáng)度準(zhǔn)則。吳順順等[15]研究了HB準(zhǔn)則下的隧道縱向變形曲線。劉立鵬、宗全兵等[16-17]分析了HB參數(shù)對(duì)邊坡穩(wěn)定性的影響。孫闖等[18]提出了簡(jiǎn)化HB應(yīng)變軟化模型，并在此模型基礎(chǔ)上，采用收斂-約束法對(duì)隧道支護(hù)結(jié)構(gòu)的穩(wěn)定性進(jìn)行分析。

由于HB準(zhǔn)則在棱線和尖點(diǎn)處存在不連續(xù)性，導(dǎo)致其有限元數(shù)值求解過(guò)程十分困難。對(duì)此，Hoek、Sofianos和Priest等[19-21]提出了HB準(zhǔn)則參數(shù)與MC準(zhǔn)則參數(shù)的等價(jià)方法。Pan、Wan和Merifield等[22-24]通過(guò)將HB準(zhǔn)則的角點(diǎn)進(jìn)行圓滑處理或修改屈服函數(shù)的方法，避免了數(shù)值求解中的奇異點(diǎn)問(wèn)題。然而，等價(jià)參數(shù)法存在一定的使用范圍(應(yīng)力值在最大圍壓值與抗拉強(qiáng)度之間)，角點(diǎn)圓滑法本質(zhì)上修改了屈服函數(shù)形式，使其在求解一些經(jīng)典問(wèn)題時(shí)(地基承載力、邊坡安全系數(shù))會(huì)產(chǎn)生偏差。Clausen、陳陪帥等[25-26]建立了主應(yīng)力空間中HB準(zhǔn)則的完全隱式應(yīng)力回映算法，該方法在處理奇異點(diǎn)問(wèn)題時(shí)具有一定的優(yōu)勢(shì)，是當(dāng)前本構(gòu)積分算法的研究熱點(diǎn)。

數(shù)值積分的困難限制了HB模型在有限元數(shù)值模擬中的應(yīng)用，基于HB模型的損傷本構(gòu)研究則更為少見(jiàn)。為了更好地反映巖體介質(zhì)的損傷演化特性，本文引入修正有效應(yīng)力原理，考慮孔隙水壓力作用，建立了基于HB準(zhǔn)則的彈塑性損傷耦合模型，從損傷和塑性兩個(gè)方面反映巖石材料的劣化機(jī)制；推導(dǎo)了主應(yīng)力空間彈塑性損傷本構(gòu)方程的完全隱式返回映射求解算法；基于ABAQUS的用戶子程序接口Umat，建立三維數(shù)值程序；最后，通過(guò)室內(nèi)單軸、三軸壓縮試驗(yàn)和工程案例對(duì)模型進(jìn)行了驗(yàn)證和應(yīng)用。

1 基于HB準(zhǔn)則的彈塑性損傷模型

1.1 考慮損傷的HB彈塑性本構(gòu)模型

基于HB模型的彈塑性損傷屈服函數(shù)和塑性勢(shì)函數(shù)表達(dá)式分別為：

式中：p為靜水壓力；J2、J3分別為第二、第三偏應(yīng)力不變量；θ為羅德角；mb、s、σci和a為HB準(zhǔn)則參數(shù)，mbg、sg、ag為對(duì)應(yīng)的塑性勢(shì)函數(shù)參數(shù)；D為損傷變量。參考已有研究，巖體損傷對(duì)HB參數(shù)的影響主要體現(xiàn)在對(duì)mb和s[27]的弱化作用。因此本文假定損傷只對(duì)參數(shù)mb、s產(chǎn)生影響，受荷過(guò)程中σci和a保持不變。

針對(duì)3個(gè)主應(yīng)力之間的大小關(guān)系，式(1)在主應(yīng)力空間中可以寫成多個(gè)屈服函數(shù)的形式：

塑性勢(shì)函數(shù)gi與對(duì)應(yīng)的屈服函數(shù)形式相同，由參數(shù)mbg、sg、ag分別替換mb、s、a即可，圖1為H-B準(zhǔn)則在主應(yīng)力空間中的圖形。

圖1 HB準(zhǔn)則在主應(yīng)力空間的圖形Fig.1 HB criterion in principal stress space

多屈服面本構(gòu)模型的流動(dòng)法則表達(dá)式為：

式中：dεp為塑性應(yīng)變?cè)隽?；dλi為塑性因子增量；m為式(4)中大于零的屈服函數(shù)的個(gè)數(shù)。

硬化規(guī)律由式(6)進(jìn)行控制：

式中：K為硬化參數(shù)向量，K=｛mb,s｝ ;Kin為mb、s的初始強(qiáng)度；Kfin為mb、s的殘余強(qiáng)度；當(dāng)Kin＞Kfin時(shí)，材料發(fā)生軟化，Kin＜Kfin時(shí)，材料發(fā)生硬化，Kin=Kfin時(shí)，材料為理想彈塑性[28]；κ為塑性內(nèi)變量，本文取κ=γp，γp為塑性剪切應(yīng)變，dγp=dε1-dε3；κfin為巖石發(fā)生破壞時(shí)的塑性內(nèi)變量。

由Kuhn-Tucker加卸載準(zhǔn)則可得：

考慮損傷的應(yīng)力-應(yīng)變本構(gòu)關(guān)系為：

式中：εe為彈性應(yīng)變張量；為損傷彈性矩陣，為四階對(duì)稱張量，分別為損傷剪切模量和體積模量。

G(D)、K(D)可以用初始剪切模量G0和初始體積模量K0表示：

由式(9)可知，損傷最終導(dǎo)致了彈性模量的弱化。如圖2所示，當(dāng)應(yīng)變達(dá)到一定的閾值時(shí)，彈性模量值隨著損傷的累積開(kāi)始減小，材料應(yīng)力-應(yīng)變關(guān)系顯現(xiàn)出非線性特征。

圖2 損傷引起剛度退化Fig.2 Stiffness degeneration caused by damage

1.2 損傷變量D的演化方程

本文將損傷變量D表示為等效塑性應(yīng)變的冪指數(shù)函數(shù)形式：

式中：α取值范圍為[0, +∞]，決定了損傷后巖石材料軟化曲線的初始斜率；β取值范圍為[0,1]，決定了巖石最大損傷值。不同α、β值下的損傷變量變化曲線如圖3所示。從圖3可以看出：α越大，損傷演化速率越慢；β越大，巖石的最終損傷值越小。

為等效塑性應(yīng)變，其表達(dá)式為：

式中，εp1、εp2、εp3分別為三個(gè)方向的主塑性應(yīng)變。

圖3 損傷參數(shù)對(duì)損傷演化規(guī)律的影響Fig.3 The influence of damage parameters on damage evolution

1.3 修正有效應(yīng)力公式

巖體是由巖石骨架和相互連通的孔隙以及其中儲(chǔ)存的流體組成的多孔介質(zhì)，在流體運(yùn)動(dòng)作用下，巖石力學(xué)性質(zhì)會(huì)發(fā)生改變。根據(jù)Biot理論，巖體有效應(yīng)力表達(dá)式為：

式中：為有效應(yīng)力；σ為總應(yīng)力；δ為ijijKronecker符號(hào)；α0為Biot系數(shù)，通常取α=1；pw和pa分別為孔隙水壓力和孔隙氣壓力；χ與飽和度和表面張力有關(guān)，一般取χ=sw，sw為孔隙水飽和度。

2 彈塑性損傷模型數(shù)值積分算法

式(1)不變量的形式在塑性計(jì)算中經(jīng)常使用，但其形式較為復(fù)雜[9]，且HB模型存在棱線和尖點(diǎn)處的奇異點(diǎn)，在這些奇異點(diǎn)處屈服函數(shù)外法線方向不唯一，導(dǎo)數(shù)不連續(xù)，使得數(shù)值實(shí)施存在一定的困難。因此本文從主應(yīng)力空間角度出發(fā)，對(duì)HB彈塑性損傷模型完全隱式的返回映射求解算法進(jìn)行推導(dǎo)：在塑性狀態(tài)求解過(guò)程中，首先對(duì)應(yīng)力空間進(jìn)行劃分，判斷應(yīng)力回映點(diǎn)的位置(單一屈服面、雙屈服面相交的棱線或者多屈服面相交的尖點(diǎn)處)，根據(jù)回映位置的不同，建立更新應(yīng)力及多個(gè)或單一塑性因子的Newton-Raphson求解式，很好地解決了空間角點(diǎn)問(wèn)題。通過(guò)ABAQUS軟件的用戶子程序接口Umat，實(shí)現(xiàn)考慮孔隙水壓作用的彈塑性損傷模型三維數(shù)值求解過(guò)程。HB彈塑性損傷模型返回映射求解算法分為彈性預(yù)測(cè)、塑性修正和損傷修正三個(gè)部分，其中，彈性預(yù)測(cè)與塑性修正均在有效應(yīng)力空間進(jìn)行，最后通過(guò)損傷修正得到最終的名義應(yīng)力。具體求解步驟如下。

2.1 主應(yīng)力空間應(yīng)力回映算法

步驟1.彈性預(yù)測(cè)

已知tn時(shí)刻應(yīng)變?cè)隽喀う舗+1、應(yīng)力σn、硬化內(nèi)變量κn和損傷變量Dn，則tn+1時(shí)刻的彈性預(yù)測(cè)應(yīng)力為：

若f1＜0，材料仍處于彈性階段，對(duì)變量進(jìn)行更新：

若f1＞0，則進(jìn)入塑性修正階段。

步驟2.塑性修正

塑性修正過(guò)程中，保持應(yīng)變?cè)隽喀う舗+1及損傷變量Dn不變，塑性狀態(tài)下應(yīng)變?cè)隽恐邪瑥椥詰?yīng)變?cè)隽喀う舉和塑性應(yīng)變?cè)隽喀う舙，塑性應(yīng)變?cè)隽勘磉_(dá)式為：

式中：m為式(4)中函數(shù)值大于零的方程個(gè)數(shù)；n+1指的是tn+1時(shí)刻。基于返回應(yīng)力的位置，式(15)具有不同的形式，當(dāng)返回應(yīng)力位于屈服面、棱線或尖點(diǎn)處時(shí)，分別需要求解1個(gè)、2個(gè)或3個(gè)塑性因子。

此時(shí)的更新應(yīng)力位于屈服面以外，需要對(duì)其進(jìn)行修正，修正量Δσp的表達(dá)式為：

更新后的應(yīng)力和內(nèi)變量表達(dá)式為：

1) 應(yīng)力返回位置的判斷

在f1與f2的交線l1上有σ1=σ2(為方便描述，對(duì)下標(biāo)n進(jìn)行省略)，因此直線方程為：

同理，f1和f6的交線l6上有σ2=σ3，直線方程為：

當(dāng)更新應(yīng)力只位于屈服面f1上時(shí)，由式(16)可得此時(shí)塑性修正應(yīng)力增量的方向h1：

式中：ν為泊松比；k的表達(dá)式如下：

圖4對(duì)應(yīng)力空間進(jìn)行了區(qū)域劃分。由圖4可以看出，h1與f1的兩條邊界線組成的邊界面確定了應(yīng)返回至f1面的應(yīng)力區(qū)域；h1、h2和f1、f2的交線確定了返回至l1線的區(qū)域，h1、h6和f1、f2的交線確定了返回值l6線的區(qū)域。

在l1上取一點(diǎn)σl1=(σ1,σ1,σ3,),l6上取一點(diǎn)σl6=(σ1,σ3,σ3)，若同時(shí)滿足：

則更新應(yīng)力位于f1面上。式中：

式中，rl1、rl6為棱線l1和l6的方向向量。

由圖4可以看到，由h1、h2組成的平面將尖點(diǎn)位置與l1分割開(kāi)來(lái)，該平面的法向量為：

由h1、h6組成的平面將尖點(diǎn)位置與l6區(qū)域分割開(kāi)來(lái)，平面法向量為：

尖點(diǎn)處應(yīng)力值為σapex=(σa,σa,σa)，σa=sσci/mb。若有：同時(shí)成立，則更新應(yīng)力位于尖點(diǎn)處。

圖4 HB準(zhǔn)則在主應(yīng)力空間的區(qū)域劃分Fig.4 Regional division for HB criterion in principal stress space

2) 應(yīng)力及內(nèi)變量求解

應(yīng)力求解過(guò)程的關(guān)鍵是對(duì)方程組式(27)進(jìn)行求解，即：

式中，R為殘差值。

當(dāng)應(yīng)力回映到f1上時(shí)，式(27)可以寫為:

求解時(shí)首先對(duì)式(28a)進(jìn)行求解，建立關(guān)于σ1,n+1和Δλ1的Newton-Raphson式：

式中，k為迭代步數(shù)。

為確保Δλ1非負(fù)，應(yīng)按照下式進(jìn)行更新：

求解出σ1,n+1和Δλ1后，代入f1=0中求解σ3,n+1，最后由式(28b)求出σ2,n+1。求解完成后，對(duì)γp進(jìn)行更新，由式(6)獲取更新硬化參數(shù)K。

當(dāng)應(yīng)力回映至l1上，滿足σ1=σ2，式(27)可以寫為：

式中：

此時(shí)，式(31)有3個(gè)未知數(shù)σ1,n+1、Δλ1和Δλ2。

若應(yīng)力回映至l6線上，滿足σ2=σ3，式(27)可以寫為：

式(33)有3個(gè)未知數(shù)σ1,n+1、Δλ1和Δλ6。

若上述結(jié)果均不滿足，當(dāng)應(yīng)力回映到尖點(diǎn)處時(shí)，更新應(yīng)力狀態(tài)滿足σ1,n+1=σ2,n+1=σ3,n+1=σa，此時(shí)需要求解方程組：

按照式(28)的形式建立Newton-Raphson求解式，分別對(duì)式(31)、式(33)和式(35)進(jìn)行求解，即可求得更新應(yīng)力值σn+1。

步驟3.損傷修正

由式(10)對(duì)損傷變量進(jìn)行更新：

因此，tn+1時(shí)刻的所對(duì)應(yīng)的名義應(yīng)力張量為：

2.2 一致切線模量

為了保證有限元方程組整體迭代求解過(guò)程中具有二階收斂速度，需要給出一致切線模量的表達(dá)式：

式中：n為屈服函數(shù)值大于零的屈服函數(shù)個(gè)數(shù)；Δεp為塑性應(yīng)變?cè)隽浚沪う薸為塑性因子增量；A為n階方陣：

式中：δij為Kronecker符號(hào)；αi的表達(dá)式為：

()c為修正彈性矩陣，為修正矩陣，表達(dá)式為：

式中，I為單位矩陣。

圖5 彈塑性損傷模型應(yīng)力回映算法求解過(guò)程Fig.5 Stress mapping algorithm solution process for elastoplastic damage

3 彈塑性損傷模型驗(yàn)證

3.1 試驗(yàn)驗(yàn)證

為驗(yàn)證本文模型的合理性，對(duì)室內(nèi)常規(guī)巖石力學(xué)試驗(yàn)進(jìn)行有限元模擬和結(jié)果對(duì)比。巖石試樣取自吉林省輝白隧道，其圍巖主要組成為混合片麻巖，包含長(zhǎng)石、石英和各種暗色礦物(云母、角閃石、輝石等)，其中長(zhǎng)石和石英含量大于50%。將同一掌子面的巖塊取回至實(shí)驗(yàn)室，然后采用姜堰市星光機(jī)電廠生產(chǎn)的ZS-200型自動(dòng)取芯機(jī)取芯，再用SHM-200雙端面磨石機(jī)打磨加工制成直徑50 mm、高100 mm的圓柱狀標(biāo)準(zhǔn)試樣。試驗(yàn)裝置采用大連海事大學(xué)和長(zhǎng)春朝陽(yáng)試驗(yàn)機(jī)廠聯(lián)合研制的多功能RLW-2000巖石三軸儀，試驗(yàn)過(guò)程采用位移加載模式進(jìn)行，加載速率為0.002 mm/s。根據(jù)文獻(xiàn)[29]的研究，由常規(guī)三軸和單軸試驗(yàn)獲取數(shù)值計(jì)算參數(shù)如下：彈性模量E=21.34 GPa，泊松比ν=0.2，重度γ=24 kN/m3，單軸抗壓強(qiáng)度σci=63 MPa，采用關(guān)聯(lián)流動(dòng)法則對(duì)塑性流動(dòng)特性進(jìn)行描述，強(qiáng)度參數(shù)mbini=mg=3，sini=sg=0.15，a=ag=0.5，殘余強(qiáng)度mbfin=0.05mbini，sfin=0.05sini，破壞時(shí)的塑性剪切應(yīng)變?chǔ)胒in=1×10-3。根據(jù)峰后段擬合損傷參數(shù)α=6×10-4，β=0.45。計(jì)算不同圍壓下巖石應(yīng)力-應(yīng)變曲線如圖6所示。

由圖6可知，不同圍壓下，數(shù)值計(jì)算所得巖體峰值偏應(yīng)力分別為68.96 MPa、72.41 MPa、78.65 MPa和90.60 MPa，與試驗(yàn)所得數(shù)值68.90 MPa、70.89 MPa、77.87 MPa和90.01 MPa十分接近，峰后軟化段的應(yīng)力大小與殘余強(qiáng)度值的選取有關(guān)。

在不考慮圍壓的情況下，本文模型對(duì)不同mb、s值下計(jì)算的巖石應(yīng)力-應(yīng)變曲線的變化規(guī)律研究見(jiàn)圖7。由圖7可知，隨著s的增大，巖石應(yīng)力應(yīng)變曲線發(fā)生上移，峰值強(qiáng)度隨之增大；改變mb的取值對(duì)巖石應(yīng)力-應(yīng)變曲線影響較小。這與文獻(xiàn)[30]的研究結(jié)論具有較好的一致性。

圖6 試驗(yàn)曲線與彈塑性損傷模型計(jì)算曲線對(duì)比Fig.6 Comparison of elastoplastic damage simulations and experimental results

圖7 不同s、mb值下單軸抗壓應(yīng)力-應(yīng)變計(jì)算曲線Fig.7 Uniaxial compression curve under different s, mb values

3.2 模型收斂特性驗(yàn)證

為分析所建模型在求解過(guò)程中的收斂效率問(wèn)題，通過(guò)Message file文件對(duì)圖6中圍壓為零的試件壓縮計(jì)算過(guò)程進(jìn)行記錄。計(jì)算過(guò)程分為11個(gè)增量步，平均每個(gè)增量步中含有2.90個(gè)迭代步，由于篇幅所限，只記錄8、9和10增量步迭代過(guò)程中的最大殘余力變化規(guī)律，如表1所示(收斂標(biāo)準(zhǔn)為默認(rèn)值5.0×10-3)。

表1 最大殘余力變化規(guī)律Table 1 Largest residual force change law

由表1可知，一致切線模量的引入使得計(jì)算過(guò)程具有二階或接近二階的收斂速度，保證了模型在實(shí)際應(yīng)用過(guò)程中的數(shù)值穩(wěn)定性。

4 工程應(yīng)用

4.1 工程概況

大連地鐵五號(hào)線04標(biāo)段工程起止里程為YK10+061.992～YK12+932.454，其中海域段長(zhǎng)度為2310 m。本文所選研究段位于K12+400～K12+753，由地勘報(bào)告可知，該段主要穿過(guò)中風(fēng)化白云質(zhì)灰?guī)r、少量中風(fēng)化板巖，海水深度約9 m～14 m，海底距隧道頂部約為29.5 m。區(qū)間采用土壓平衡式盾構(gòu)機(jī)，實(shí)行單洞雙線雙層襯砌的開(kāi)挖方案，管片內(nèi)徑10.8 m、外徑11.8 m、環(huán)寬2.0 m、管片厚度50 cm。

4.2 有限元計(jì)算模型

盾構(gòu)施工是一個(gè)復(fù)雜的三維問(wèn)題，包括地應(yīng)力初始平衡、土體開(kāi)挖、開(kāi)挖面土體應(yīng)力釋放、盾構(gòu)機(jī)行進(jìn)、管片襯砌安裝、盾尾注漿壓力、注漿層硬化等施工過(guò)程。每一掘進(jìn)步的施工流程如圖8所示。

圖8 盾構(gòu)施工流程Fig.8 Shield construction process

根據(jù)實(shí)際工程概況建立有限元模型如圖9所示。模型大小為60 m×30 m×60 m，共劃分為17640個(gè)單元和17568個(gè)節(jié)點(diǎn)。管片厚度0.5 m，注漿層厚度0.1 m，開(kāi)挖半徑5.9 m。數(shù)值模擬段主要穿過(guò)中風(fēng)化白云巖，依據(jù)地勘情況并查閱文獻(xiàn)[13]中所提供的表格確定巖體力學(xué)參數(shù)和主要支護(hù)參數(shù)如表2所示，mb與s的殘余強(qiáng)度均取為峰值強(qiáng)度的1/10，破壞時(shí)的塑性剪切應(yīng)變?chǔ)胒in=1×10-3。模型兩側(cè)施加沿z軸方向成線性變化的孔隙水，設(shè)置最大水頭壓力分別為1.5 MPa、3.5 MPa(實(shí)際工況)和4.5 MPa。盾構(gòu)機(jī)經(jīng)過(guò)27個(gè)開(kāi)挖步，由y=0 m推進(jìn)至y=45 m。

圖9 有限元計(jì)算模型Fig.9 FEM calculation model

表2 模型計(jì)算參數(shù)Table 2 Calculating parameters of model

4.3 數(shù)值結(jié)果分析

由2.2節(jié)可知，損傷參數(shù)α和β控制了損傷變化速率和最終損傷值的大小。不同損傷參數(shù)下，計(jì)算出的巖體材料剛度退化程度也不同，為研究損傷參數(shù)對(duì)數(shù)值模擬結(jié)果的影響，計(jì)算不同α、β下，由y=0 m至y=45 m模型開(kāi)挖完成后，地表沉降值如圖10所示。

由圖10(a)可知，地表沉降值隨著β的減小而增大，當(dāng)β=0.5時(shí)，計(jì)算出現(xiàn)不收斂。由圖10(b)可知，隨著α的增大，地表沉降值不斷減小，當(dāng)α＞0.05時(shí)，沉降值幾乎不再發(fā)生變化。計(jì)算結(jié)果表明，損傷參數(shù)的選取對(duì)計(jì)算結(jié)果影響較大，模型進(jìn)行工程應(yīng)用時(shí)宜首先根據(jù)監(jiān)測(cè)數(shù)據(jù)對(duì)損傷參數(shù)進(jìn)行反演，獲得合理的取值，從而保證計(jì)算結(jié)果的準(zhǔn)確性。

圖10 不同損傷參數(shù)下地表沉降值Fig.10 Ground surface settlement under different damage parameter

本文選用差異進(jìn)化算法(DE)對(duì)損傷參數(shù)進(jìn)行反演，首先建立基于位移值的目標(biāo)函數(shù)：

然后通過(guò)MATLAB語(yǔ)言反復(fù)調(diào)用DE算法和ABAQUS求解器對(duì)目標(biāo)函數(shù)式(42)進(jìn)行尋求，最終獲得較為合理的損傷參數(shù)，具體方法見(jiàn)文獻(xiàn)[31]。為保證結(jié)果的唯一性，本文選取了3個(gè)斷面處共9個(gè)監(jiān)測(cè)點(diǎn)的數(shù)據(jù)進(jìn)行反演，最終得到損傷參數(shù)為α=0.002，β=0.51。

不同孔隙水壓下，計(jì)算開(kāi)挖損傷區(qū)分布情況如圖11所示。隨著外部孔隙水壓的增大，隧道最終開(kāi)挖損傷區(qū)逐漸增大，最大損傷值由2.992×10-3增至3.00×10-3。計(jì)算結(jié)果表明，孔隙水的存在威脅了隧道開(kāi)挖穩(wěn)定性，尤其對(duì)于海底隧道，漲落潮過(guò)程中，地下水位易發(fā)生變化，造成圍巖損傷加劇，引發(fā)工程事故。

圖11 不同孔隙水壓下?lián)p傷區(qū)分布圖Fig.11 Damage profile under different pore pressure

圖12為最大孔隙水壓值為1.5 MPa、3.5 MPa、4.5 MPa的開(kāi)挖過(guò)程中監(jiān)測(cè)斷面節(jié)點(diǎn)上(參照?qǐng)D9中的標(biāo)注)的損傷值和位移值變化情況。

由圖12可知，盾構(gòu)施工過(guò)程中，監(jiān)測(cè)斷面處的損傷值逐步增長(zhǎng)，變化速率隨著開(kāi)挖面的臨近不斷增大，當(dāng)開(kāi)挖面遠(yuǎn)離監(jiān)測(cè)斷面時(shí)，損傷累積速率逐步減緩，最終變?yōu)榱?，損傷值也不再發(fā)生變化。同一監(jiān)測(cè)點(diǎn)處的最終損傷值隨著孔隙水壓的增大而加大，不同開(kāi)挖步下的損傷增長(zhǎng)速率不同。隨著孔隙水壓的增大，拱頂點(diǎn)處的最終損傷值分別為0.2988、0.2998和0.2999；拱底點(diǎn)處的最終損傷值分別為0.2980、0.2997和0.2995；拱腰點(diǎn)處的最終損傷值分別為0.2698、0.2861和0.2928，因此可知，孔隙水壓的存在加劇了圍巖的損傷程度。而在相同水壓下，最終損傷值總有拱頂點(diǎn)＞拱底點(diǎn)＞拱腰點(diǎn)，建議在施工監(jiān)測(cè)時(shí)嚴(yán)密關(guān)注拱頂處是否有襯砌開(kāi)裂的情況。

圖12 監(jiān)測(cè)點(diǎn)損傷值和位移值變化曲線Fig.12 Monitoring damage and displacement values change curve

由位移計(jì)算結(jié)果可知，盾構(gòu)掘進(jìn)過(guò)程中，位移變化速率在臨近開(kāi)挖面時(shí)會(huì)突然增大，隨著開(kāi)挖面的遠(yuǎn)離速率不斷減小，位移值變化趨于穩(wěn)定。不同位置處的位移變化形式不同，拱頂監(jiān)測(cè)點(diǎn)主要產(chǎn)生沉降，拱底監(jiān)測(cè)點(diǎn)主要產(chǎn)生隆起，而拱腰監(jiān)測(cè)點(diǎn)主要產(chǎn)生收斂現(xiàn)象。隨著孔隙水壓的增大，拱頂點(diǎn)處的最終位移值分別為0.0296 m、0.0420 m和0.0500 m；拱底點(diǎn)處的最終位移值分別為0.0133 m、0.0345 m和0.0451 m，拱腰點(diǎn)處的最終位移值分別為0.0234 m、0.0348 m和0.0436 m。各點(diǎn)處位移值均會(huì)受到孔隙水壓的影響，因此在施工過(guò)程中，應(yīng)緊密關(guān)注監(jiān)測(cè)點(diǎn)位移值變化情況，避免海水漲落引發(fā)的水壓變化對(duì)工程造成危害。

由于監(jiān)測(cè)點(diǎn)布設(shè)技術(shù)的限制，只能對(duì)監(jiān)測(cè)面開(kāi)挖后的監(jiān)測(cè)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，但由圖12可看出，位移值變化規(guī)律與孔隙水壓為3.5 MPa(實(shí)際工況)時(shí)的計(jì)算值較為吻合，表明所建模型的計(jì)算精度能夠滿足工程需要。

本文所建彈塑性損傷模型能夠計(jì)算得到盾構(gòu)掘進(jìn)過(guò)程中巖體的損傷值和位移變化規(guī)律，可以為盾構(gòu)施工參數(shù)(掘進(jìn)速度、注漿壓力和土倉(cāng)壓力等)合理調(diào)整提供依據(jù)。

5 結(jié)論

本文建立了基于H-B準(zhǔn)則的巖體彈塑性損傷模型，并給出了模型的數(shù)值積分算法。對(duì)該模型進(jìn)行了試驗(yàn)驗(yàn)證和工程應(yīng)用，得出以下結(jié)論。

(1) 本文引入修正有效應(yīng)力原理來(lái)考慮孔隙水壓力的作用，建立了基于HB準(zhǔn)則的巖體彈塑性損傷本構(gòu)模型。既具有強(qiáng)度非線性的優(yōu)點(diǎn)，又能考慮損傷引起的剛度退化和塑性導(dǎo)致的流動(dòng)兩種破壞機(jī)制的耦合作用。

(2) 為解決HB彈塑性準(zhǔn)則應(yīng)力空間奇異點(diǎn)導(dǎo)致難以數(shù)值求解的問(wèn)題，本文從主應(yīng)力空間出發(fā)，推導(dǎo)了HB彈塑性損傷模型的包括彈性預(yù)測(cè)、塑性修正和損傷修正三個(gè)步驟的隱式返回映射算法。并通過(guò)ABAQUS軟件的用戶子程序接口Umat成功開(kāi)發(fā)了程序。

(3) 通過(guò)試驗(yàn)驗(yàn)證了本模型的合理性?；诮⒌腍B彈塑性損傷模型和程序，對(duì)大連地鐵工程的海底盾構(gòu)隧道進(jìn)行了三維數(shù)值計(jì)算，能夠合理地反映盾構(gòu)施工工序、不同海水水壓導(dǎo)致的圍巖位移和損傷演化規(guī)律。計(jì)算結(jié)果可以為盾構(gòu)施工參數(shù)的合理調(diào)整提供依據(jù)。