周　利，蘇　凱，周亞峰，文喜雨，伍鶴皋

（1.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢　430072；2.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究院，湖北 武漢　430010；3.長(zhǎng)江巖土工程總公司，湖北 武漢　430010）

1　研究背景

隨著社會(huì)經(jīng)濟(jì)與工程技術(shù)的發(fā)展，水工隧洞無(wú)論在數(shù)量上還是在規(guī)模上都達(dá)到了前所未有的高度，帶來(lái)了巨大的社會(huì)與經(jīng)濟(jì)效益，同時(shí)也面臨諸多問(wèn)題。例如：隧洞在內(nèi)水壓力作用下，襯砌混凝土受拉開(kāi)裂，并在施工缺陷以及溫度微裂紋等影響下，形成性態(tài)不一的宏觀裂縫，進(jìn)而導(dǎo)致內(nèi)水外滲［1-2］，此類襯砌結(jié)構(gòu)均稱為透水襯砌。在高內(nèi)壓下襯砌環(huán)向拉應(yīng)力超過(guò)混凝土抗拉強(qiáng)度后，襯砌混凝土損傷破壞，表現(xiàn)為微裂紋萌生、擴(kuò)展并貫通，直至產(chǎn)生宏觀裂縫，洞內(nèi)高壓水沿裂縫通道外滲使得滲透水壓增大，在影響結(jié)構(gòu)應(yīng)力場(chǎng)分布的同時(shí)又加劇了損傷，如此相互影響，是一個(gè)復(fù)雜的多場(chǎng)耦合過(guò)程。卞康等［3］采用增量變彈性損傷剛度計(jì)算迭代的方法分析了襯砌水壓致裂過(guò)程中的耦合效應(yīng)；Simanjuntak等［4］探討了襯砌開(kāi)裂過(guò)程中水力耦合效應(yīng)對(duì)隧洞滲流量的影響；Bian等［5］在分析實(shí)際工程時(shí)，提出了一種滲流-應(yīng)力耦合方法；Fahimifar等［6］在平面軸對(duì)稱假定下，提出了一種解析方法以分析襯砌開(kāi)裂過(guò)程中的耦合效應(yīng)。然而從我國(guó)進(jìn)行的6次較大規(guī)模的隧洞襯砌直接壓水試驗(yàn)結(jié)果來(lái)看（詳見(jiàn)表1），當(dāng)隧洞內(nèi)壓增大到某一數(shù)值時(shí)，裂縫出現(xiàn)，此后無(wú)新的裂縫產(chǎn)生，只是縫寬隨內(nèi)壓增大，裂縫呈少而寬分布［7］，而這一現(xiàn)象在現(xiàn)有的研究中卻鮮有體現(xiàn)?？紤]到水工隧洞襯砌裂縫分布的這一特點(diǎn)，劉秀珍［8］按照彈性地基曲梁結(jié)構(gòu)推導(dǎo)了圓形有壓隧洞襯砌裂縫寬度的計(jì)算方法；沈威［9］采用室內(nèi)水工隧洞結(jié)構(gòu)模型試驗(yàn)研究了襯砌開(kāi)裂過(guò)程，同樣得出了裂縫少而寬的分布特征。

表1　隧洞壓水試驗(yàn)參數(shù)及結(jié)果［7］

襯砌裂縫出現(xiàn)后，洞內(nèi)高壓水勢(shì)必會(huì)沿著裂縫向外流動(dòng)，使得圍巖變形增大，當(dāng)襯砌因鋼筋的約束作用徑向變形與圍巖不一致時(shí)，兩者可能出現(xiàn)脫離現(xiàn)象，單獨(dú)承擔(dān)自身范圍內(nèi)的水壓力［10-11］，體現(xiàn)了襯砌與圍巖典型的有條件聯(lián)合承載機(jī)理［12-13］。此外，工程監(jiān)測(cè)資料顯示，隧洞發(fā)電運(yùn)行時(shí)，襯砌鋼筋應(yīng)力普遍小于設(shè)計(jì)預(yù)期，并沒(méi)有達(dá)到既安全又經(jīng)濟(jì)的效果，說(shuō)明現(xiàn)有的襯砌設(shè)計(jì)方法仍存在一定的問(wèn)題，以致高估了鋼筋應(yīng)力值。例如：在廣蓄電站高壓水工隧洞充水加壓初期，鋼筋應(yīng)力隨著內(nèi)壓增大；襯砌開(kāi)裂后，鋼筋應(yīng)力值變化明顯放緩甚至出現(xiàn)下降的趨勢(shì)，最大值僅為49.5 MPa［14］。類似的情況也發(fā)生在周寧水電站水工隧洞充排水試驗(yàn)中［15］。然而目前的襯砌設(shè)計(jì)理論、解析及數(shù)值方法尚不能很好地解釋這一現(xiàn)象［11］。

鑒于此，本文基于混凝土塑性損傷模型和透水襯砌理論，考慮在高內(nèi)壓作用下襯砌損傷開(kāi)裂、滲流與應(yīng)力的耦合作用以及襯砌與圍巖有條件聯(lián)合承載特性，基于大型通用有限元軟件平臺(tái)ABAQUS進(jìn)行二次開(kāi)發(fā)，通過(guò)實(shí)用程序GETVRM獲取材料損傷，利用子程序USDFLD實(shí)現(xiàn)材料滲透系數(shù)隨損傷的動(dòng)態(tài)更新，提出并實(shí)現(xiàn)了透水襯砌滲流-應(yīng)力-損傷耦合算法，模擬分析了水工隧洞初次充水加壓過(guò)程，以此研究襯砌損傷開(kāi)裂過(guò)程中的多場(chǎng)耦合作用，探討鋼筋應(yīng)力的演化特征，為水工隧洞襯砌設(shè)計(jì)問(wèn)題提供一定參考。

2　理論基礎(chǔ)

水工隧洞水壓力由襯砌和圍巖承擔(dān)。本文采用彈塑性本構(gòu)關(guān)系模擬圍巖，以Mohr-Coulomb準(zhǔn)則作為巖體的破壞準(zhǔn)則，由于研究重點(diǎn)在于透水襯砌的耦合分析，限于篇幅，在此不做贅述。對(duì)于襯砌，本文采用間接耦合的方法對(duì)隧洞初次充水期襯砌進(jìn)行滲流-應(yīng)力-損傷耦合分析，考慮滲流場(chǎng)通過(guò)等效節(jié)點(diǎn)荷載影響應(yīng)力場(chǎng)分布，在結(jié)構(gòu)自身所承擔(dān)的荷載以及等效節(jié)點(diǎn)荷載的作用下計(jì)算應(yīng)力場(chǎng)后即可獲取結(jié)構(gòu)應(yīng)力、應(yīng)變以及材料損傷，并建立損傷與滲流之間的耦合關(guān)系，根據(jù)損傷實(shí)時(shí)調(diào)整材料滲透系數(shù)并返回滲流場(chǎng)進(jìn)行計(jì)算，直至計(jì)算收斂、加載結(jié)束，詳見(jiàn)圖1。

圖1　滲流-應(yīng)力-損傷耦合關(guān)系

2.1滲流場(chǎng)與等效節(jié)點(diǎn)荷載的數(shù)值解法假設(shè)水和材料不可壓縮，水在材料孔隙間的滲流滿足微分方程［16-17］：

式中：θ為含水率；分別是x、y、z方向上的滲透系數(shù)；t為時(shí)間變量；水力勢(shì)H=p/γ+z，其中p為孔隙水壓力，z為位置水頭，γ為水容重。

計(jì)算空間域經(jīng)過(guò)離散、插值和積分，應(yīng)用Galerkin加權(quán)余量法及格林公式可得求解滲流場(chǎng)的矩陣方程：

式中：［Ks］為總滲透矩陣；{h}為節(jié)點(diǎn)水頭列向量；{A}為對(duì)滲流邊界積分得到的節(jié)點(diǎn)荷載。

滲流場(chǎng)計(jì)算完成后，水力梯度所產(chǎn)生的水荷載以滲透體積力的方式作用在結(jié)構(gòu)內(nèi)部，按式（3）換算成作用在模型節(jié)點(diǎn)上的等效荷載［3］：

式中：［N］為插值函數(shù)；Ω為節(jié)點(diǎn)積分域。

2.2應(yīng)力場(chǎng)的計(jì)算將計(jì)算空間域離散化，得到滲流與應(yīng)力的耦合方程：

式中：［Km］為結(jié)構(gòu)剛度矩陣；{u}為節(jié)點(diǎn)位移列向量；{F}為結(jié)構(gòu)的節(jié)點(diǎn)荷載。

2.3混凝土塑性損傷模型及鋼筋模型相較于預(yù)裂縫模型與分布裂縫模型，塑性損傷模型能夠較好地模擬混凝土材料軟化過(guò)程、卸載剛度退化及反復(fù)加載過(guò)程［18-19］，并且不必預(yù)先假定混凝土開(kāi)裂位置，程序會(huì)根據(jù)結(jié)構(gòu)的受力狀況判定并計(jì)算結(jié)構(gòu)各個(gè)部位的損傷開(kāi)裂情況，亦不會(huì)出現(xiàn)大面積開(kāi)裂的情形，在工程分析中應(yīng)用廣泛［20］。多軸應(yīng)力狀態(tài)下，混凝土應(yīng)力σ計(jì)算如下：

式中：d為單元損傷；為初始彈性剛度；ε為單元應(yīng)變；εpl為塑性應(yīng)變。

單元損傷d由等效塑性應(yīng)變表示如下：

式中：t、c代表拉、壓狀態(tài)。

在ABAQUS中，混凝土損傷后，應(yīng)力、損傷與塑性應(yīng)變的關(guān)系可以按照列表形式輸入，本文襯砌材料采用C25混凝土，抗拉、抗壓強(qiáng)度取標(biāo)準(zhǔn)值1.78 MPa、16.7 MPa，結(jié)合《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》（GB50010-2010）附錄C確定相關(guān)參數(shù)，具體應(yīng)力、損傷與塑性應(yīng)變之間的關(guān)系如圖2所示。

采用埋藏式模型模擬鋼筋，在迭代計(jì)算中鋼筋由線性桿單元模擬，選取理想彈塑性本構(gòu)關(guān)系，應(yīng)力σs與應(yīng)變?chǔ)舠關(guān)系如下：

圖2　塑性損傷演化曲線

式中：Es為鋼筋彈性模量；fy為屈服應(yīng)力。

2.4損傷開(kāi)裂與滲流的耦合關(guān)系采用塑性損傷模型時(shí)，損傷將導(dǎo)致混凝土結(jié)構(gòu)剛度矩陣的退化，以此模擬混凝土損傷開(kāi)裂對(duì)結(jié)構(gòu)受力狀態(tài)的影響，同時(shí)混凝土損傷開(kāi)裂對(duì)結(jié)構(gòu)的滲透特性亦會(huì)帶來(lái)顯著的影響。材料單元由損傷相與未損傷相兩部分組成，單元滲透系數(shù)計(jì)算如下［21］：

式中：km為未損傷相滲透系數(shù)；kd為損傷相滲透系數(shù)；損傷相的塑性體積應(yīng)變?yōu)閱卧苄泽w積應(yīng)變；若單元同時(shí)出現(xiàn)拉損傷與壓損傷，分別計(jì)算對(duì)應(yīng)狀態(tài)下的單元滲透系數(shù)，k取較大值。

考慮到脆性材料出現(xiàn)宏觀裂縫后，其滲透特性較開(kāi)裂前會(huì)出現(xiàn)突增現(xiàn)象，因此引入突跳因子ξ以計(jì)算損傷相滲透系數(shù)［22］：

式中：對(duì)于壓剪損傷ξ=100，對(duì)于拉損傷

3　技術(shù)路線

在耦合計(jì)算時(shí)，首先根據(jù)式（1）—（3）計(jì)算滲流場(chǎng)與等效節(jié)點(diǎn)荷載，然后根據(jù)式（4）計(jì)算應(yīng)力場(chǎng)，并依據(jù)塑性損傷演化曲線獲取材料損傷，根據(jù)式（11）和式（12）調(diào)整材料滲透系數(shù)后重新計(jì)算滲流場(chǎng)，如此循環(huán)直至計(jì)算收斂、加載結(jié)束。

為提高計(jì)算效率并保證計(jì)算精度［3］，本文內(nèi)壓分3個(gè)階段進(jìn)行施加以模擬隧洞初次充水運(yùn)行工況：（1）通過(guò)試算確定混凝土即將出現(xiàn)損傷時(shí)所對(duì)應(yīng)的內(nèi)壓，該內(nèi)壓采取一次性加載的方式進(jìn)行施加；（2）襯砌損傷后，內(nèi)壓水頭增量以較小值逐級(jí)緩慢施加；（3）貫穿縫出現(xiàn)后，內(nèi)水外滲通道形成，適當(dāng)增大內(nèi)壓水頭增量。每施加一次內(nèi)壓水頭增量，須保證在該級(jí)內(nèi)壓作用下耦合迭代完全，滲流場(chǎng)達(dá)到穩(wěn)定，以保證耦合精度。

本文基于ABAQUS平臺(tái)，通過(guò)Fortran語(yǔ)言調(diào)用子程序USDFLD進(jìn)行二次開(kāi)發(fā)，使用USDFLD時(shí)，必須調(diào)用實(shí)用程序GETVRM以獲取材料點(diǎn)的數(shù)據(jù)，以此完成耦合迭代過(guò)程，具體步驟如下（詳見(jiàn)圖3）：

①施加上覆壓力以模擬初始應(yīng)力場(chǎng)的影響，計(jì)算隧洞開(kāi)挖后應(yīng)力場(chǎng)，考慮隧洞開(kāi)挖后，應(yīng)力釋放完全后進(jìn)行襯砌支護(hù)，即襯砌不承擔(dān)開(kāi)挖釋放荷載，計(jì)算襯砌支護(hù)后應(yīng)力場(chǎng)；

②襯砌內(nèi)壁施加增量水頭計(jì)算滲流場(chǎng)，將各單元滲透體積力轉(zhuǎn)換為等效節(jié)點(diǎn)荷載；

③將等效節(jié)點(diǎn)荷載施加于模型節(jié)點(diǎn)計(jì)算得到應(yīng)力、應(yīng)變以及材料損傷，調(diào)用實(shí)用程序GETVRM獲取單元損傷值，對(duì)損傷單元根據(jù)式（11）和式（12）調(diào)用子程序USDFLD調(diào)整其滲透系數(shù)；

④根據(jù)新的滲透系數(shù)計(jì)算滲流場(chǎng)，檢查滲流場(chǎng)是否穩(wěn)定。當(dāng)前后兩次滲流場(chǎng)計(jì)算結(jié)果滿足式（13）時(shí)，滲流場(chǎng)達(dá)到穩(wěn)定狀態(tài)；如不滿足，則繼續(xù)求解等效節(jié)點(diǎn)荷載施加于模型節(jié)點(diǎn)進(jìn)行耦合迭代，直至滲流場(chǎng)穩(wěn)定；

式中：Pn+1、Pn為迭代n+1和n次后節(jié)點(diǎn)孔壓；δ為迭代控制閥值，為提高計(jì)算效率并保證計(jì)算精度，將δ取為 1%［3］。

⑤施加下一級(jí)增量水頭，在同一級(jí)內(nèi)壓水頭下重復(fù)②—④步，直至加載完成、計(jì)算結(jié)束。

圖3　耦合迭代流程圖

4　耦合算例分析

4.1模型參數(shù)及邊界條件如圖4為一豎井段圓形隧洞軸對(duì)稱模型，X、Y軸平行于隧洞斷面，位于水平面內(nèi)，豎直方向?yàn)閆軸，向上為正，襯砌內(nèi)半徑r1為4.2 m，外半徑r2為5 m，圍巖沿徑向取30倍開(kāi)挖洞徑［23］，為300 m，模型厚10 m，襯砌與圍巖均由八節(jié)點(diǎn)實(shí)體單元模擬，共計(jì)9865個(gè)節(jié)點(diǎn)和7888個(gè)單元，襯砌與圍巖之間共節(jié)點(diǎn)。在襯砌內(nèi)外兩側(cè)鋪設(shè)直徑為28 mm的鋼筋，鋼筋到隧洞中心距離R1、R2分別為4.3 m、4.9 m。模型材料參數(shù)如表2所示。

圖4　有限元模型

隧洞埋深100 m，沿Z軸負(fù)向施加100 m的巖體上覆壓力于模型上表面以模擬初始應(yīng)力場(chǎng)的影響。在初始時(shí)刻，襯砌及圍巖處于飽和狀態(tài)，結(jié)構(gòu)內(nèi)部為飽和滲流，襯砌及圍巖孔隙率設(shè)為1.5%、10%。采用增量水頭法計(jì)算滲流場(chǎng)及應(yīng)力場(chǎng)［24］，忽略重力對(duì)增量水頭影響范圍即水力半徑的影響，針對(duì)增量水頭的影響范圍將圍巖外圈設(shè)為零水頭邊界［25］。模型頂面以及底面設(shè)定為不透水邊界，襯砌內(nèi)壁設(shè)為定水頭邊界，計(jì)算過(guò)程中增大襯砌內(nèi)壁的定水頭邊界得到增量水頭以模擬充水加壓過(guò)程。計(jì)算滲流場(chǎng)時(shí)，模型施加全位移約束；計(jì)算應(yīng)力場(chǎng)時(shí)，模型底部和外圈節(jié)點(diǎn)施加法向位移約束。不考慮模型自身以及洞內(nèi)水體的重力，只考慮襯砌及圍巖中水體的重力。內(nèi)壓水頭從0.5 MPa逐級(jí)加載至1.1 MPa，每級(jí)增量水頭控制在0.025～0.1 MPa，保證每一級(jí)內(nèi)壓水頭下耦合迭代完全。

表2　材料參數(shù)

4.2襯砌損傷開(kāi)裂特征通過(guò)試算，內(nèi)壓小于0.5 MPa時(shí)，襯砌不會(huì)出現(xiàn)損傷。內(nèi)壓增至0.65 MPa時(shí)，損傷分布均勻，從內(nèi)向外損傷值由大變小，最大值為0.108，詳見(jiàn)圖5（a）。在0.67 MPa下迭代4次達(dá)到穩(wěn)定，襯砌環(huán)向出現(xiàn)3條貫穿縫，等間距分布于0°（隧洞右腰）、120°和240°位置（逆時(shí)針），裂縫分布少而寬，詳見(jiàn)圖5（b）。此后直至加載完畢，不再出現(xiàn)新的裂縫，裂縫處損傷值隨內(nèi)壓持續(xù)增大至0.999，且裂縫間混凝土中部出現(xiàn)拉損傷，原因在于裂縫出現(xiàn)后，襯砌結(jié)構(gòu)類似于彈性地基（曲）梁，使得裂縫間襯砌外側(cè)出現(xiàn)反彎點(diǎn)，拉應(yīng)力增大，出現(xiàn)明顯的拉損傷區(qū)［26］，詳見(jiàn)圖5（c）。通過(guò)耦合迭代計(jì)算得到的襯砌損傷開(kāi)裂特征與隧洞壓水試驗(yàn)結(jié)果（表1）吻合較好［7］，說(shuō)明本文建立的耦合模型能夠較好地反映襯砌的開(kāi)裂過(guò)程。

圖5　襯砌損傷開(kāi)裂情況

4.3滲透孔壓演化特征為避開(kāi)模型頂面以及底面邊界條件的影響，選取模型中部位置對(duì)滲透孔壓演化特征進(jìn)行分析。內(nèi)壓為0.5 MPa時(shí)，襯砌未進(jìn)入塑性狀態(tài)，無(wú)損傷，襯砌承擔(dān)了大部分的水壓力，圍巖內(nèi)孔壓相對(duì)較小，不超過(guò)0.115 MPa，詳見(jiàn)圖6（a）。在內(nèi)壓達(dá)到0.67 MPa，襯砌開(kāi)裂前，襯砌內(nèi)側(cè)出現(xiàn)損傷，損傷區(qū)域孔壓較大，圍巖內(nèi)孔壓有所增大，最大值為0.257 MPa，詳見(jiàn)圖6（b）；3條貫穿縫出現(xiàn)后，內(nèi)水外滲，滲流場(chǎng)出現(xiàn)明顯的重分布現(xiàn)象，圍巖內(nèi)孔壓升高，在襯砌環(huán)向0°、120°及240°裂縫位置滲流場(chǎng)變化明顯，孔壓梯度明顯降低，裂縫處圍巖承受的水壓接近內(nèi)壓水頭，孔壓自裂縫位置向圍巖深部逐漸消散，滲流場(chǎng)分布與襯砌損傷演化情況相對(duì)應(yīng)，基本滿足等間距分布，詳見(jiàn)圖 6（c）。

圖6　隧洞滲流場(chǎng)分布（單位：Pa）

隨著內(nèi)壓的增大，襯砌損傷峰值與滲透系數(shù)峰值同步變化，表示材料損傷與滲透系數(shù)實(shí)現(xiàn)了同步更新。裂縫處混凝土損傷與滲透系數(shù)在襯砌開(kāi)裂時(shí)出現(xiàn)陡增現(xiàn)象，此后隨內(nèi)壓持續(xù)增大，加載完成時(shí)，滲透系數(shù)相較于初始1×10-9m/s提高了約3個(gè)量級(jí)，裂縫成為內(nèi)水外滲的主要通道，詳見(jiàn)圖7。

4.4鋼筋應(yīng)力演化特征充水加壓過(guò)程中，裂縫處鋼筋應(yīng)力隨內(nèi)壓增大，而未開(kāi)裂位置鋼筋全程處于低應(yīng)力狀態(tài)，應(yīng)力值小于14 MPa。當(dāng)內(nèi)壓小于0.67 MPa時(shí)，鋼筋應(yīng)力緩慢增長(zhǎng)；貫穿縫形成后，裂縫位置混凝土應(yīng)力釋放，大部分荷載轉(zhuǎn)由鋼筋承擔(dān)，裂縫處鋼筋應(yīng)力突增至50.120 MPa；當(dāng)內(nèi)壓從0.67 MPa增至0.9 MPa時(shí)，裂縫處鋼筋應(yīng)力緩慢增至62.770 MPa，當(dāng)內(nèi)壓達(dá)到1.0 MPa時(shí)，裂縫處鋼筋應(yīng)力進(jìn)一步增大，在內(nèi)壓1.1 MPa迭代穩(wěn)定后裂縫處鋼筋應(yīng)力達(dá)到228.800 MPa，詳見(jiàn)圖8。原因在于貫穿縫出現(xiàn)后，裂縫處的鋼筋所發(fā)揮的作用類似于錨栓，單獨(dú)承擔(dān)拉力，而未開(kāi)裂位置鋼筋與襯砌混凝土共同承擔(dān)襯砌環(huán)向拉力，混凝土承擔(dān)大部分拉力，鋼筋所承擔(dān)的拉力較小。

圖7　襯砌損傷峰值及滲透系數(shù)峰值隨內(nèi)壓的演化

圖8　鋼筋應(yīng)力演化特征

5　襯砌與圍巖有條件聯(lián)合承載機(jī)理研究

襯砌開(kāi)裂前，透水性遠(yuǎn)弱于圍巖，內(nèi)部水力梯度較大，承擔(dān)大部分的水荷載，而圍巖內(nèi)部水力梯度較小，承擔(dān)的水荷載較小，因此襯砌沿徑向朝外變形的趨勢(shì)更大，襯砌與圍巖處于壓緊狀態(tài)；襯砌開(kāi)裂后，透水性顯著增強(qiáng)，內(nèi)部水力梯度大幅降低，所承擔(dān)的水荷載明顯減小，此時(shí)圍巖的透水性相對(duì)較弱，內(nèi)部水力梯度相對(duì)較大，承擔(dān)大部分的水荷載，因此圍巖沿徑向朝外變形的趨勢(shì)更為顯著，圍巖拉著襯砌向外變形，同時(shí)襯砌受到內(nèi)置鋼筋的約束，在襯砌與圍巖交界面上產(chǎn)生拉應(yīng)力，當(dāng)拉應(yīng)力超過(guò)交界面粘結(jié)強(qiáng)度時(shí)，襯砌與圍巖將出現(xiàn)脫離現(xiàn)象，表現(xiàn)出有條件聯(lián)合承載特性［10，27-29］。

在原有模型的基礎(chǔ)上，對(duì)計(jì)算模型進(jìn)行局部修改，考慮襯砌與圍巖之間不共節(jié)點(diǎn)，在襯砌與圍巖之間建立面-面接觸單元以模擬襯砌與圍巖有條件聯(lián)合承載特性，接觸單元切向摩擦系數(shù)取為0.5，局部襯砌-接觸單元-圍巖模型如圖9所示。此時(shí)襯砌與圍巖接觸面只存在法向力與切向摩擦力，當(dāng)法向力為壓力時(shí)，兩者處于壓緊狀態(tài)，一旦法向受拉，兩者脫離，不再傳遞法向應(yīng)力，也不再滿足位移連續(xù)條件，即默認(rèn)襯砌與圍巖交界面粘結(jié)強(qiáng)度為零，接觸壓力與間隙的關(guān)系如圖10所示。在應(yīng)力場(chǎng)計(jì)算中，接觸面的作用主要是傳遞法向應(yīng)力；在滲流場(chǎng)計(jì)算中，接觸可以保證孔隙流體在不同材料之間的滲透作用。

5.1襯砌與圍巖的接觸狀態(tài)內(nèi)壓增加到0.625 MPa前，襯砌未開(kāi)裂，襯砌滲透特性遠(yuǎn)弱于圍巖，襯砌承擔(dān)了大部分水荷載，在0.625 MPa內(nèi)壓下迭代前，10 m厚度范圍內(nèi)襯砌內(nèi)壁內(nèi)水外滲總流量?jī)H為1.109 L/min，此時(shí)襯砌與圍巖不存在脫離現(xiàn)象，詳見(jiàn)圖11（a）；在0.625 MPa內(nèi)壓下迭代完成后，襯砌環(huán)向在290°位置（襯砌右腰為0°位置，逆時(shí)針旋轉(zhuǎn)）形成一條徑向貫穿裂縫，同時(shí)襯砌除貫穿縫外的其他部位的損傷值在迭代過(guò)程中有所增大，滲透性增強(qiáng)，兩部分因素的作用下使得襯砌內(nèi)壁內(nèi)水外滲總流量陡增至6.850 L/min。此后隨著內(nèi)壓增大，不再出現(xiàn)新的裂縫，襯砌損傷開(kāi)裂特征與隧洞壓水試驗(yàn)一致（表1）。由于裂縫處水力梯度的降低以及襯砌向外變形受到內(nèi)置鋼筋的約束，襯砌與圍巖首先在290°開(kāi)裂位置出現(xiàn)脫離現(xiàn)象，產(chǎn)生的空隙由滲透水充滿，脫離范圍在迭代過(guò)程中逐步擴(kuò)大，迭代11次后達(dá)到穩(wěn)定狀態(tài)，脫離范圍達(dá)75%，詳見(jiàn)圖11（b）。隨著內(nèi)壓增至0.635 MPa，迭代4次后達(dá)到穩(wěn)定狀態(tài)，襯砌與圍巖已經(jīng)完全脫離，詳見(jiàn)圖11（c）。

圖9　襯砌-接觸單元-圍巖模型

圖10　接觸壓力與間隙的關(guān)系

圖11　襯砌與圍巖接觸狀態(tài)（深色為脫離部分）

5.2鋼筋應(yīng)力演化特征引入接觸單元后，襯砌鋼筋應(yīng)力表現(xiàn)出明顯的差異。襯砌開(kāi)裂前鋼筋應(yīng)力較小，內(nèi)壓達(dá)到0.625 MPa時(shí)，襯砌開(kāi)裂，鋼筋應(yīng)力呈整體上升趨勢(shì)，且存在較大波動(dòng)，迭代過(guò)程中裂縫處鋼筋應(yīng)力最大值達(dá)到327.500 MPa，此后隨迭代過(guò)程鋼筋應(yīng)力逐漸趨于穩(wěn)定，并在水壓力緩慢增加的過(guò)程中，鋼筋應(yīng)力基本穩(wěn)定并呈緩慢降低的趨勢(shì)。在1.1 MPa內(nèi)壓作用下，鋼筋應(yīng)力最大值為56.560 MPa，整體量值較小，詳見(jiàn)圖12。出現(xiàn)此現(xiàn)象的原因在于襯砌開(kāi)裂后，襯砌與圍巖逐步脫離，圍巖成為水荷載的主要承擔(dān)對(duì)象，而襯砌承擔(dān)的水荷載逐漸減小，導(dǎo)致鋼筋應(yīng)力量值較小。同時(shí)與完全聯(lián)合承載情況不同，完全聯(lián)合承載時(shí)襯砌在開(kāi)裂之后仍然會(huì)隨著圍巖向外擴(kuò)張變形，導(dǎo)致鋼筋應(yīng)力持續(xù)增大。相比之下襯砌與圍巖相互脫離，即有條件聯(lián)合承載極有可能是隧洞發(fā)電運(yùn)行過(guò)程中鋼筋應(yīng)力普遍小于設(shè)計(jì)預(yù)期的原因，并且鋼筋應(yīng)力的這一演化特征和量值很好地解釋了相關(guān)工程的監(jiān)測(cè)結(jié)果［14-15，30］。

圖12　有條件聯(lián)合承載時(shí)鋼筋應(yīng)力演化特征

6　結(jié)語(yǔ)

本文基于混凝土塑性損傷模型與透水襯砌理論，通過(guò)等效節(jié)點(diǎn)荷載實(shí)現(xiàn)了滲流場(chǎng)對(duì)應(yīng)力場(chǎng)的作用，建立了損傷開(kāi)裂與滲流之間的耦合關(guān)系，以此提出了透水襯砌滲流-應(yīng)力-損傷耦合算法?；贏BAQUS平臺(tái)進(jìn)行二次開(kāi)發(fā)，通過(guò)Fortran語(yǔ)言調(diào)用實(shí)用程序GETVRM獲取材料點(diǎn)的損傷值，利用子程序USDFLD實(shí)現(xiàn)了材料滲透系數(shù)隨損傷的動(dòng)態(tài)更新，完成了隧洞初次充水期間襯砌損傷開(kāi)裂的耦合迭代過(guò)程。研究了襯砌損傷開(kāi)裂、滲透孔壓以及鋼筋應(yīng)力的演化特征，并采用面-面接觸單元探討了襯砌與圍巖有條件聯(lián)合承載機(jī)理，可以發(fā)現(xiàn)：（1）本文提出的耦合算法所得到的襯砌損傷開(kāi)裂特征與襯砌壓水試驗(yàn)結(jié)果吻合，表現(xiàn)為襯砌在一定內(nèi)壓下開(kāi)裂，此后不再出現(xiàn)新的裂縫，只在裂縫間襯砌外側(cè)出現(xiàn)拉損傷，較好地反映了水工隧洞初次充水加壓工況中襯砌損傷開(kāi)裂的過(guò)程。（2）襯砌開(kāi)裂前后，滲流場(chǎng)表現(xiàn)出明顯的差異。開(kāi)裂后，裂縫處孔壓梯度明顯降低，滲透系數(shù)呈現(xiàn)數(shù)量級(jí)的增大，裂縫成為內(nèi)水外滲的主要通道，這與工程普遍認(rèn)知是相符的，說(shuō)明本文所提出的耦合算法是可行的，具有一定的推廣性。（3）完全聯(lián)合承載時(shí)，襯砌與圍巖之間徑向變形協(xié)調(diào)，裂縫處鋼筋應(yīng)力持續(xù)增大，未開(kāi)裂區(qū)鋼筋全程維持在低應(yīng)力水平；考慮有條件聯(lián)合承載特性后，襯砌開(kāi)裂之后與圍巖逐步脫離，圍巖成為承擔(dān)水荷載的主體，襯砌鋼筋應(yīng)力出現(xiàn)應(yīng)力回縮現(xiàn)象，在加載完畢時(shí)應(yīng)力量值維持在較小水平，揭示了有條件聯(lián)合承載是隧洞發(fā)電運(yùn)行過(guò)程中鋼筋應(yīng)力普遍小于設(shè)計(jì)預(yù)期的原因，解釋了相關(guān)工程的監(jiān)測(cè)結(jié)果，說(shuō)明本文所提出的耦合算法可以為隧洞襯砌設(shè)計(jì)問(wèn)題提供一定參考。

參考文獻(xiàn)：

［1］楊林德，丁文其.滲水高壓引水隧洞襯砌的設(shè)計(jì)研究［J］.巖石力學(xué)與工程學(xué)報(bào)，1997，16（2）：17-22.

［2］杜小凱，任青文，陳偉 .有壓引水隧洞內(nèi)水外滲作用研究［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2008，40（5）：63-68.

［3］卞康，肖明.水工隧洞襯砌水壓致裂過(guò)程的滲流-損傷-應(yīng)力耦合分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29（S2）：3769-3776.

［4］SIMANJUNTAK T D Y F，MARENCE M，MYNETT A E，et al.Mechanical-hydraulic interaction in the lining cracking process of pressure tunnels［J］.Int.J.Hydropow.Dams，2013，20（5）：98-105.

［5］BIAN K，XIAO M，CHEN J.Study on coupled seepage and stress fields in the concrete lining of the underground pipe with high water pressure［J］.Tunnelling and Underground Space Technology，2009，24：287-295.

［6］FAHIMIFAR A，MOHAMMAD R Z.A theoretical solution for analysis of tunnels below groundwater considering the hydraulic-mechanical coupling［J］.Tunnelling and Underground Space Technology，2009，24：634-646.

［7］張有天.水工隧洞建設(shè)的經(jīng)驗(yàn)和教訓(xùn)（上）［J］.貴州水力發(fā)電，2001，15（4）：76-84.

［8］劉秀珍.關(guān)于圓形有壓水工隧洞襯砌裂縫計(jì)算的探討［J］.水力發(fā)電，1983，12：39-46.

［9］沈威.圓形水工壓力隧洞襯砌變形特性與限裂設(shè)計(jì)研究［D］.大連：大連理工大學(xué)，2012.

［10］蘇凱，伍鶴皋，周創(chuàng)兵 .內(nèi)水壓力下水工隧洞襯砌與圍巖承載特性研究［J］.巖土力學(xué)，2010，31（8）：2407-2412.

［11］ZHOU Y，SU K，WU H.Hydro-mechanical interaction analysis of high pressure hydraulic tunnel［J］.Tunnelling and Underground Space Technology，2015，47：28-34.

［12］GRAZIANI A，BOLDINI D.Influence of hydro-mechanical coupling on tunnel response in clays［J］.Journal of Geotechnical&Geoenvironmental Engineering，2012，138（3）：415-418.

［13］YOO C.Interaction between tunneling and groundwater-numerical investigation using three dimensional stress-pore pressure coupled analysis［J］.Journal of Geotechnical&Geoenvironmental Engineering，2005，131（2）：240-250.

［14］葉冀升.廣蓄電站水工高壓隧洞設(shè)計(jì)施工的若干問(wèn)題［J］.水力發(fā)電學(xué)報(bào)，1998（2）：39-50.

［15］李煊明.周寧水電站透水襯砌高壓鋼筋混凝土岔管研究［J］.水電站設(shè)計(jì)，2007，23（2）：15-17.

［16］吳夢(mèng)喜.飽和—非飽和土中滲流Richards方程有限元算法［J］.水利學(xué)報(bào)，2009，40（10）：1274-1279.

［17］紀(jì)佑軍，劉建軍，程林松 .考慮流-固耦合的隧道開(kāi)挖數(shù)值模擬［J］.巖土力學(xué)，2011，32（4）：1229-1233.

［18］LUBLINER J，OLIVER J，OLLER S，et al.A plastic-damage model for concrete［J］.International Journal of Solids&Structures，1989，25（3）：299-326.

［19］LEE J，F(xiàn)ENVES G L.Plastic-damage model for cyclic loading of concrete structures［J］.Journal of Engineering Mechanics，1998，124（8）：892-900.

［20］蘇凱，張偉，伍鶴皋，等.考慮摩擦接觸特性的鋼襯鋼筋混凝土管道承載機(jī)理研究［J］.水利學(xué)報(bào)，2016，47（8）：1070-1078.

［21］賈善坡，陳衛(wèi)忠，于洪丹，等.泥巖隧道施工過(guò)程中滲流場(chǎng)與應(yīng)力場(chǎng)全耦合損傷模型研究［J］.巖土力學(xué)，2009，30（1）：19-26.

［22］周亞峰.水工隧洞滲流-應(yīng)力-損傷/開(kāi)裂耦合理論及方法研究［D］.武漢：武漢大學(xué)，2016.

［23］周亞峰，蘇凱，伍鶴皋.水工隧洞鋼筋混凝土襯砌外水壓力取值方法研究［J］.巖土力學(xué)，2014，35（S2）：198-203.

［24］張有天，張武功.隧洞水荷載的靜力計(jì)算［J］.水利學(xué)報(bào)，1980（3）：52-62.

［25］潘家錚，唐家驥，張有天，等.隧洞水荷載靜力計(jì)算（學(xué)術(shù)討論）［J］.水利學(xué)報(bào)，1981（5）：75-81.

［26］李寧，劉波 .城門(mén)洞型壓力隧洞襯砌限裂配筋設(shè)計(jì)方法研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2013（12）：2409-2414.

［27］張有天，張武功，王鐳.再論隧洞水荷載的靜力計(jì)算［J］.水利學(xué)報(bào)，1985（3）：21-31.

［28］BOBET A，NAM S W.Stresses around pressure tunnels with semi-permeable liners［J］.Rock Mechanics and Rock Engineering，2007，40（3）：287-315.

［29］FERNANDEZ G.Behavior of pressure tunnels and guidelines for liner design［J］.Journal of Geotechnical Engineering，1994，120（10）：1768-1791.

［30］劉素琴，張秀麗，劉啟釗，等.高水頭鋼筋混凝土岔管結(jié)構(gòu)計(jì)算與觀測(cè)對(duì)比分析［J］.水利水電技術(shù)，2006，37（8）：24-28.