基于近場(chǎng)動(dòng)力學(xué)的隧道襯砌結(jié)構(gòu)漸進(jìn)裂損數(shù)值模擬分析

馬云峰，徐林生

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

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

1 近場(chǎng)動(dòng)力學(xué)與有限元混合模型

1.1 ABAQUS計(jì)算模型離散化

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

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

(5)

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

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

(6)

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

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

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

(7)

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

(8)

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

2 襯砌裂縫近場(chǎng)動(dòng)力學(xué)可行性研究

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

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

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

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

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

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

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

本文主要以天水市八盤(pán)山隧道的維護(hù)完善工程為主要實(shí)例，進(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ù)模塊實(shí)現(xiàn)有限元域和近場(chǎng)動(dòng)力學(xué)域的混合建模，將隧道襯砌開(kāi)裂部位劃分為FEM域和PD域，見(jiàn)圖4。根據(jù)式(2)cohesive單元的剛度為：

(9)

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

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

3.1 局部偏壓

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

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

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

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

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

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

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

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

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

3.2 背后空洞

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

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

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

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

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

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

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

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

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

3.3 襯砌厚度不足

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

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

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

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

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

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

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

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

4 結(jié) 論

1) 通過(guò)與原型試驗(yàn)對(duì)比分析，近場(chǎng)動(dòng)力學(xué)結(jié)合有限元的方法可以適用于模擬襯砌裂縫的開(kāi)裂與擴(kuò)展過(guò)程量化分析研究。

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

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

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