許夢飛，姜諳男，史洪濤，李興盛

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

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

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

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

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

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

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

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

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

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

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

塑性勢函數(shù)gi與對應(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á)式為：

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

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

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

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

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

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

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

由式(9)可知，損傷最終導(dǎo)致了彈性模量的弱化。如圖2所示，當(dāng)應(yīng)變達(dá)到一定的閾值時，彈性模量值隨著損傷的累積開始減小，材料應(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分別為三個方向的主塑性應(yīng)變。

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

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

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

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

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

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

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

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

已知tn時刻應(yīng)變增量Δεn+1、應(yīng)力σn、硬化內(nèi)變量κn和損傷變量Dn，則tn+1時刻的彈性預(yù)測應(yīng)力為：

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

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

步驟2.塑性修正

塑性修正過程中，保持應(yīng)變增量Δεn+1及損傷變量Dn不變，塑性狀態(tài)下應(yīng)變增量中包含彈性應(yīng)變增量Δεe和塑性應(yīng)變增量Δεp，塑性應(yīng)變增量表達(dá)式為：

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

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

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

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

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

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

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

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

圖4對應(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上取一點σl1=(σ1,σ1,σ3,),l6上取一點σl6=(σ1,σ3,σ3)，若同時滿足：

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

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

由圖4可以看到，由h1、h2組成的平面將尖點位置與l1分割開來，該平面的法向量為：

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

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

圖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ān)鍵是對方程組式(27)進(jìn)行求解，即：

式中，R為殘差值。

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

求解時首先對式(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。求解完成后，對γp進(jìn)行更新，由式(6)獲取更新硬化參數(shù)K。

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

式中：

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

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

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

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

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

步驟3.損傷修正

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

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

2.2 一致切線模量

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

式中：n為屈服函數(shù)值大于零的屈服函數(shù)個數(shù)；Δεp為塑性應(yīng)變增量；Δλi為塑性因子增量；A為n階方陣：

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

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

式中，I為單位矩陣。

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

3 彈塑性損傷模型驗證

3.1 試驗驗證

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

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

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

圖6 試驗曲線與彈塑性損傷模型計算曲線對比Fig.6 Comparison of elastoplastic damage simulations and experimental results

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

3.2 模型收斂特性驗證

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

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

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

4 工程應(yīng)用

4.1 工程概況

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

4.2 有限元計算模型

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

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

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

圖9 有限元計算模型Fig.9 FEM calculation model

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

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

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

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

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

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

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

不同孔隙水壓下，計算開挖損傷區(qū)分布情況如圖11所示。隨著外部孔隙水壓的增大，隧道最終開挖損傷區(qū)逐漸增大，最大損傷值由2.992×10-3增至3.00×10-3。計算結(jié)果表明，孔隙水的存在威脅了隧道開挖穩(wěn)定性，尤其對于海底隧道，漲落潮過程中，地下水位易發(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的開挖過程中監(jiān)測斷面節(jié)點上(參照圖9中的標(biāo)注)的損傷值和位移值變化情況。

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

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

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

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

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

5 結(jié)論

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

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

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

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