馬云峰，徐林生

(1.重慶交通大學(xué) 土木工程學(xué)院，重慶400074；2.重慶水利電力職業(yè)技術(shù)學(xué)院 市政工程系，重慶402160)

截止到2016年年底我國公路隧道共計15 181處、總長度14 037.9 km，其中特長隧道共計815處、總長度3 622.7 km，長隧道共計3 520處、總長度6 045.5 km[1]。伴隨著隧道建設(shè)事業(yè)的磅礴發(fā)展，運營隧道的病害問題日益突出，尤其是襯砌裂損給隧道的安全運營帶來巨大隱患。隧道襯砌裂損的實質(zhì)主要是由于作用于隧道結(jié)構(gòu)的荷載超出其承載能力產(chǎn)生的，但最終都會反映在隧道的二襯結(jié)構(gòu)上。近幾年國內(nèi)外諸多學(xué)者對隧道襯砌病害的探索研究取得了一系列的研究成效[2]。但是對于襯砌裂縫的大多數(shù)研究還是依托工程實例，尚未能對裂縫的開裂機理、裂紋的動態(tài)擴展、裂縫的影響因素和分布規(guī)律形成完整系統(tǒng)的認(rèn)識[3]。盡管目前數(shù)值分析對襯砌裂縫的研究取得了較大的進(jìn)展，但大部分還是針對襯砌裂縫部位及裂損部位的應(yīng)力和變形研究展開的，而且在模擬精度上參數(shù)的選取和模型精確性還需進(jìn)一步研究[4]。有限單元法(FEM)由傳統(tǒng)連續(xù)介質(zhì)力學(xué)發(fā)展而來，對于求解斷裂損傷問題需要預(yù)判開裂位置及擴展方向，而且對于網(wǎng)格的依賴性較強，計算過程較為復(fù)雜，計算結(jié)果往往與實際情況差異性較大。近場動力學(xué)理論(PD)在求解非連續(xù)性問題時基于非局部作用的積分思想，避免了由于微分方程求解而產(chǎn)生的奇異性問題，但其求解的計算效率遠(yuǎn)低于有限單元法[5]。近場動力學(xué)和傳統(tǒng)有限元混合建模的方法在分析襯砌損傷和開裂方面能夠有效提高計算精度和計算效率，加快收斂速度。將計算模型中的FEM域和PD域根據(jù)應(yīng)變能密度等效原則進(jìn)行耦合。在驗證模型的可行性基礎(chǔ)上對襯砌裂縫的擴展和分布規(guī)律進(jìn)行分析探討。

1 近場動力學(xué)與有限元混合模型

1.1 ABAQUS計算模型離散化

將模型的PD域離散為具有密度和體積信息的物質(zhì)點，得到控制方程在時間和空間上的離散化形式[6]：

(1)

式中：ρ為物質(zhì)點的物理密度；xi和xj為近場范圍內(nèi)任兩物質(zhì)點的位置；ui和uj為兩物質(zhì)點的位移；n為時間步長；Vj為j處物質(zhì)點體積；bi為物質(zhì)點外荷載的密度；f為本構(gòu)力函數(shù)。

1.2 本構(gòu)關(guān)系和破壞準(zhǔn)則

本構(gòu)模型采用修正的PMB模型，物質(zhì)點對視為cohesive單元(鍵)連接。當(dāng)‖ξ‖δ時，物質(zhì)點對相互作用的本構(gòu)力函數(shù)f表示形式為[7]：

(2)

當(dāng)‖ξ‖>δ時，f=0。式中：ξ為物質(zhì)點對相對位置；δ為近場范圍；η為物質(zhì)點對的相對位移；E為彈性模量；c為材料的微模量；μ為物質(zhì)點產(chǎn)生的位移；s為物質(zhì)點對的相對伸長率。

PD理論中對于損傷斷裂的描述引入了標(biāo)量函數(shù)σ和φ(x,t)來表示[8]：

(3)

式中，s0為物質(zhì)點對的臨界伸長率。當(dāng)μ=0時，表示物質(zhì)點對之間的相互作用消失，微觀的鍵斷裂映射在宏觀上產(chǎn)生了裂紋，當(dāng)物質(zhì)點的破壞累積成面時，就形成了宏觀的斷裂，該點近場范圍內(nèi)裂損度用φ(x,t)表示[9]：

(4)

式中，H為空間域內(nèi)物質(zhì)點x的進(jìn)場范圍。

近場動力學(xué)應(yīng)變能密度的表達(dá)式為[10]：

(5)

式中：α、k′ 為待求參數(shù)；ω為影響函數(shù)；θ為體積膨脹率；ed為伸長狀態(tài)的偏張量；W為應(yīng)變能密度。

由傳統(tǒng)應(yīng)變能密度和近場動力學(xué)應(yīng)變能密度相等可以推出：

(6)

式中：k、G分別為體積模量和剪切模量；m為加權(quán)體積標(biāo)量。

為了模擬隧道襯砌開裂狀態(tài)建立基于能量的損傷耦合方程[11]：

Ψ(e,φ)=W[(1-φ)e]

(7)

由式(1),(4)～(7)可得：

(8)

式中：M為沿著變形鍵方向的單位矢量。

2 襯砌裂縫近場動力學(xué)可行性研究

日本隧道研究中心曾于2000年開展過1:1襯砌模型室內(nèi)試驗，在襯砌外表面周身布設(shè)千斤頂來模擬圍巖抗力和外加荷載，主要針對襯砌拱頂集中荷載和拱頂30°范圍的襯砌背后空洞等工況作用下襯砌的開裂情況進(jìn)行研究。具體參數(shù)詳見文獻(xiàn)[12]。為了驗證近場動力學(xué)模擬襯砌裂縫開裂過程，選取30°范圍的襯砌背后空洞工況進(jìn)行對比分析。為了貼合原型試驗，本文的數(shù)值模擬不考慮隧道仰拱的作用，同時忽略隧道埋深對其襯砌結(jié)構(gòu)開裂的影響，選取0°～180°范圍內(nèi)襯砌進(jìn)行研究(見圖1)。

圖1 隧道模型示意圖Fig.1 Sketch of tunnel model

分析過程中主要考慮了三種破壞模式，一種是由于混凝土達(dá)到極限壓應(yīng)變而產(chǎn)生的壓潰； 一種是由于混凝土達(dá)到極限拉應(yīng)變而產(chǎn)生的拉裂；一種是由于裂縫的延展而導(dǎo)致的隧道襯砌承載力失效。基于應(yīng)變能密度等效的原則，裂縫的擴展準(zhǔn)則采用應(yīng)變能密度因子準(zhǔn)則。襯砌的損傷采取最大應(yīng)力準(zhǔn)則，損傷演化選取基于能量的、線性軟化的、混合模式的指數(shù)損傷演化準(zhǔn)則。

圖2～3分別是近場動力學(xué)方法和室內(nèi)原型試驗得到的荷載-拱頂位移對比曲線和裂縫分布狀態(tài)對比圖。對比結(jié)果顯示，二者在荷載-拱頂位移曲線中趨勢相同，曲線擬合較好，數(shù)值模擬前期加載產(chǎn)生的拱頂位移大于原型試驗；數(shù)值模擬后期加載產(chǎn)生的拱頂位移小于原型試驗。對原型試驗過程分析，千斤頂模擬外加荷載和圍巖抗力的不均勻性是造成結(jié)果差異的主要原因。裂縫分布狀態(tài)圖中顯示二者在開裂位置和裂縫分布上的結(jié)果基本吻合。因此，采用近場動力學(xué)方法來模擬襯砌裂縫漸進(jìn)裂損過程是可行的。

圖2 數(shù)值模擬與原型試驗荷載-拱頂位移對比曲線Fig.2 Comparing load-vault displacement curves between numerical simulation and prototype test

圖3 數(shù)值分析與原型試驗裂縫分布狀態(tài)對比圖Fig.3 Comparison of fracture distribution between numerical analysis and prototype test

3 隧道襯砌漸進(jìn)裂損近場動力學(xué)分析

本文主要以天水市八盤山隧道的維護完善工程為主要實例，進(jìn)行襯砌結(jié)構(gòu)漸進(jìn)裂損的數(shù)值模擬。模型相關(guān)參數(shù)：隧道半徑為5 250 mm，襯砌厚度為55 cm厚C20混凝土，斷裂能為G1C=G2C=G3C=60 N/m，Young模量為21 GPa，Poisson比為0.3，圍巖壓力為22 kPa，采用彈簧單元模擬圍巖抗力。利用ABAQUS中的鑲嵌單元技術(shù)模塊實現(xiàn)有限元域和近場動力學(xué)域的混合建模，將隧道襯砌開裂部位劃分為FEM域和PD域，見圖4。根據(jù)式(2)cohesive單元的剛度為：

(9)

圖4 FEM節(jié)點和PD物質(zhì)點連接示意圖Fig.4 Diagram for connection between FEM node and PD material point

選取對我們所關(guān)心的裂縫擴展區(qū)域采用四結(jié)點二維cohesive單元，其他區(qū)域采用四結(jié)點雙線性平面應(yīng)變四邊形單元，并對開裂區(qū)域網(wǎng)格加密。

3.1 局部偏壓

偏壓現(xiàn)象主要是由于荷載的不對稱性產(chǎn)生的，本文不考慮地形的不對稱產(chǎn)生的偏壓，僅分析荷載局部偏壓對隧道襯砌開裂的影響，計算模型見圖5。以隧道左側(cè)80°為荷載偏壓基線，選取10°(80°～70°)、20°(80°～60°)、30°(80°～50°)、40°(80°～40°)四種局部單偏壓工況進(jìn)行分析研究。

圖5 偏壓計算模型圖示Fig.5 Graphics of bias pressure computing model

圖6反映了裂縫的擴展?fàn)顟B(tài)，給出襯砌四種工況下裂紋水平層次集的等值面。隨著偏壓荷載的不斷增大，在襯砌內(nèi)表面偏壓荷載的中心附近出現(xiàn)一條剪拉狀態(tài)的主裂縫，并伴隨荷載的增加主裂縫繼續(xù)擴展由剪拉狀態(tài)轉(zhuǎn)變?yōu)榧魤籂顟B(tài)，其兩側(cè)出現(xiàn)新裂縫的發(fā)育并呈現(xiàn)剪壓狀。從初始圍巖壓力作用到襯砌承載力失效的過程中，外表面的混凝土首先被壓潰，內(nèi)表面混凝土由張裂轉(zhuǎn)變?yōu)閴簼?，最終隨著裂縫貫通導(dǎo)致襯砌結(jié)構(gòu)失效。

圖6 偏壓作用裂縫分布狀態(tài)圖Fig.6 Fracture distribution state diagram under bias pressure

伴隨著襯砌裂縫的擴展，在偏壓范圍內(nèi)出現(xiàn)不同程度的位移，根據(jù)內(nèi)表面主裂縫的擴展?fàn)顟B(tài)，以裂縫的徑向深度來分析位移和裂縫擴展的關(guān)系見圖7。

圖7 偏壓作用主裂縫擴展法向深度與位移關(guān)系圖Fig.7 Relation diagram of normal depth and displacement of main crack propagation under bias pressure

伴隨著裂縫深度的不斷增加產(chǎn)生的環(huán)向位移和徑向位移也逐漸增大。隨著偏壓范圍的增加，襯砌的位移隨裂縫深度擴展變化明顯，直到裂縫貫穿襯砌，位移的變化沒有平緩的過程，可見襯砌的破壞主要是由于裂縫擴展而導(dǎo)致的襯砌承載力失效。

圖8給出了襯砌內(nèi)表面主裂縫尖端單元應(yīng)變能和時間歷程關(guān)系圖。隨著偏壓范圍增加，裂縫尖端的應(yīng)變能也隨之減小，說明襯砌承載力失效所需的破壞荷載隨偏壓范圍增加而減小。在加荷的中前期裂縫主要經(jīng)歷成核狀態(tài)并呈現(xiàn)穩(wěn)態(tài)擴展，隨荷載的逐步增加裂縫由剪拉狀態(tài)向剪壓狀態(tài)轉(zhuǎn)變，應(yīng)變能呈現(xiàn)減少趨勢，而后出現(xiàn)失穩(wěn)擴展?fàn)顟B(tài)，最終達(dá)到斷裂。

圖8 偏壓作用主裂縫應(yīng)變能時程關(guān)系曲線Fig.8 Strain energy time-history curve of main crack under bias pressure

3.2 背后空洞

本文采用彈簧單元的失效來模擬襯砌背后空洞(見圖9)，初始狀態(tài)圍巖壓力均勻分布于襯砌周身，每一個分析步開始移除空洞范圍襯砌圍壓作用，重新分布圍巖壓力，將空洞圍壓在每一個分析步中按時間增量均勻加載在其他襯砌上。逐漸增加圍巖壓力直至裂縫擴展導(dǎo)致襯砌承載力失效。在此過程不考慮拱腰范圍空洞大小，主要分析拱頂10°、20°、30°、45°對稱范圍空洞大小(4種工況)對襯砌開裂的影響。

圖9 背后空洞計算模型圖示Fig.9 Diagram of calculation model for back cavity

圖10給出了襯砌四種工況下裂紋水平層次集的等值面。隨著圍壓的增加，拱頂10°、20°工況下，拱肩和拱腳出現(xiàn)了較為明顯的裂縫擴展，拱肩外表面出現(xiàn)主裂縫有向拱頂延展趨勢并最終貫穿襯砌，拱頂20°工況下，拱頂開始出現(xiàn)明顯裂縫擴展；到拱頂30°、45°工況時襯砌周身均出現(xiàn)了明顯的裂縫擴展，裂縫的增擴現(xiàn)象較為明顯，其中在拱腰外表面出現(xiàn)了主裂縫貫穿襯砌。在此過程中空洞范圍的邊界附近最先出現(xiàn)裂縫的擴展，裂縫多呈剪拉狀態(tài)，拱腳位置的裂縫主要呈剪壓狀態(tài)。拱頂10°、20°工況下拱腳外表面混凝土首先壓潰，隨著空洞范圍的不斷增大，拱頂內(nèi)表面混凝土也出現(xiàn)壓潰，最終隨著裂縫的貫通，整個襯砌承載力失效。

圖10 背后空洞裂縫分布狀態(tài)圖Fig.10 State graph of fracture distribution in back cavity

在加載破壞的過程中，空洞 20°工況時拱頂外表面出現(xiàn)較為明顯的裂縫擴展，隨著空洞范圍增大，裂縫由剪拉狀態(tài)向剪壓狀態(tài)轉(zhuǎn)變，并產(chǎn)生不同程度位移，按照時間增量分析拱頂內(nèi)表面位移發(fā)展關(guān)系，見圖11。

圖11 背后空洞作用下主裂縫擴展法向深度與位移關(guān)系Fig.11 Normal depth and displacement of main crack propagation under back cavity

初期加載過程位移曲線出現(xiàn)了明顯上揚，主要是材料的接觸變形和自壓縮產(chǎn)生的形變，隨著荷載的持續(xù)增加，在加載過程的中后期拱頂內(nèi)表面的混凝土出現(xiàn)較為明顯的壓潰，位移也出現(xiàn)了小幅度的回落。整個時間歷程中，環(huán)向位移明顯大于徑向位移，隨著空洞范圍增大，產(chǎn)生的位移量也逐步增加。

以貫穿襯砌的主裂縫來分析在整個時間歷程中應(yīng)變能的變化規(guī)律見圖12。不同拱頂空洞范圍的應(yīng)變能隨時程的發(fā)展呈現(xiàn)逐級增大趨勢，可見裂縫的擴展過程沒有出現(xiàn)受力狀態(tài)的轉(zhuǎn)變。空洞范圍較大的工況下襯砌裂縫產(chǎn)生應(yīng)變能較大，裂縫擴展速率也較大。

圖12 背后空洞作用下主裂縫應(yīng)變能時程關(guān)系曲線Fig.12 Time-history curve of strain energy of main crack under back cavity

3.3 襯砌厚度不足

本文以八盤山隧道襯砌厚度55 cm為基本模型，襯砌厚度30 cm為襯砌厚度不足模型，隧道的外徑不變，分析從初始圍壓加載至30 cm厚度襯砌發(fā)生裂縫的延展而導(dǎo)致的襯砌承載力失效時整個過程襯砌裂縫開裂狀態(tài)和分布規(guī)律，計算模型見圖13。

圖13 襯砌厚度不足計算模型圖示Fig.13 Diagram of calculation model for insufficient lining thickness

圖14為30 cm和55 cm兩種厚度襯砌模型破壞時襯砌裂縫的分布狀態(tài)圖。在持續(xù)加載的過程中，兩種結(jié)構(gòu)開裂部位較為相似，30 cm厚度襯砌模型發(fā)生破壞時，裂縫多以張開型裂縫為主，拱頂部位產(chǎn)生貫穿裂縫，55 cm厚度襯砌模型在拱腰部位產(chǎn)生最大長度裂縫為8 cm，裂縫寬度0.25 cm。

圖15為兩種模型中最大裂縫的擴展深度隨時間歷程的關(guān)系曲線。當(dāng)最大裂縫擴展至相同長度，30 cm厚度襯砌需要荷載較小。從兩條曲線斜率的發(fā)展趨勢可以看出，30 cm厚度襯砌裂縫的擴展速率遠(yuǎn)大于55 cm厚度襯砌。

圖14 不同的襯砌厚度裂縫分布狀態(tài)圖Fig.14 State diagram of crack distribution in different lining thicknesses

圖15 最大裂縫擴展深度的時程曲線Fig.15 Time-history curve of maximum fracture propagation depth

為了更好反應(yīng)襯砌厚度不足情況下襯砌的漸進(jìn)裂損狀態(tài)，輸出整個模型應(yīng)變能隨時間歷程的變化，見圖16。55 cm厚度襯砌在整個加荷過程中應(yīng)變能變化較為平緩，裂紋擴展基本處于穩(wěn)態(tài)擴展階段。30 cm厚度襯砌在加荷中期后期，曲線出現(xiàn)斷層式跳躍并以更大的斜率不斷增長，說明在此階段裂縫已經(jīng)進(jìn)入失穩(wěn)擴展?fàn)顟B(tài)，并且裂縫在襯砌周身大面積擴展，直至裂縫的擴展而導(dǎo)致襯砌承載能力失效。

圖16 整個襯砌應(yīng)變能釋放的時程曲線Fig.16 Time-history curve of strain energy release for whole lining

4 結(jié) 論

1) 通過與原型試驗對比分析，近場動力學(xué)結(jié)合有限元的方法可以適用于模擬襯砌裂縫的開裂與擴展過程量化分析研究。

2) 局部偏壓作用下，隧道襯砌結(jié)構(gòu)在偏壓范圍內(nèi)產(chǎn)生的裂縫較為集中，在偏壓荷載中心附近出現(xiàn)一條主裂縫，隨著主裂縫的擴展，裂縫的受力狀態(tài)由剪拉狀向剪壓狀轉(zhuǎn)變，在此階段應(yīng)變能的釋放出現(xiàn)明顯的降低。整個過程中產(chǎn)生的環(huán)向位移遠(yuǎn)大于徑向位移。

3) 拱頂背后空洞作用下，在10°、20°工況下，襯砌結(jié)構(gòu)周身產(chǎn)生的裂縫較少，隨著空洞范圍的擴大，裂縫的增擴現(xiàn)象逐漸明顯。在模擬空洞過程中，其應(yīng)力產(chǎn)生重分布；由于材料的接觸變形和自壓縮，襯砌在加荷初期產(chǎn)生較為明顯的環(huán)向位移和徑向位移。

4) 對比分析襯砌厚度不同情況，在加荷過程中產(chǎn)生裂縫的位置大致相同，但是襯砌厚度不足情況下裂縫的擴展速率和裂縫的增擴速率較大，并且整個襯砌的產(chǎn)生的應(yīng)變能遠(yuǎn)大于正常襯砌模型。