趙勐，肖明?，陳俊濤，金漢城

（1.水資源與水電工程科學國家重點實驗室（武漢大學），湖北 武漢 430072；2.水工巖石力學教育部重點實驗室（武漢大學），湖北 武漢 430072）

為緩解巴基斯坦國家電網(wǎng)嚴重缺電的局面，中國投資并幫助巴基斯坦規(guī)劃建設了一大批水利水電工程，如卡洛特水電站、SK 水電站和正在規(guī)劃設計的阿扎德帕坦水電站，從而形成了為數(shù)眾多的地下水工隧洞群.水工隧洞往往深覆于山體中，大多具有大尺度、大埋深、洞線長等特點，不可避免地要穿越包括強震區(qū)和層狀巖體區(qū)等在內(nèi)的復雜地質(zhì)區(qū)域，面臨突出的抗震安全穩(wěn)定問題[1].在層狀巖體中進行地下工程建設，尤其對于軟硬互層狀巖體而言，當巖層傾角較大，結(jié)構(gòu)面發(fā)育明顯時，層間剪切、擠壓破碎帶較為常見，是水工隧洞抗震的薄弱部分.一旦發(fā)生地震，巖體層間極易發(fā)生剪切滑移破壞，對隧洞結(jié)構(gòu)造成不可逆的損傷破壞.“5.12”汶川大地震震害調(diào)查表明，深埋于地質(zhì)條件較差部位的隧洞結(jié)構(gòu)極易發(fā)生襯砌開裂及錯位等損傷破壞[2].因此，研究層狀巖體中水工隧洞地震響應特性和破壞機理具有重要的現(xiàn)實意義.

層狀巖體水工隧洞地震響應分析主要包括兩方面內(nèi)容：1）地震動的輸入方法；2）層間動力相互作用模擬.在已有的針對地下隧洞群進行的動力時程分析中，多是假設地震動為從模型底部豎直向上入射.但是根據(jù)近年來強震動觀測記錄的統(tǒng)計，發(fā)現(xiàn)基巖場地的地震波入射角平均為60°，從而引起結(jié)構(gòu)的非一致性變形[3].杜修力等[4]研究了地震波斜入射條件下隧洞洞身段的地震響應特征，結(jié)果表明地震波斜入射時隧洞地震響應規(guī)律與豎直入射時明顯不同.李山有等[5]研究了地震波斜入射條件下豎直、傾斜臺階地形引起的波形轉(zhuǎn)換，分析說明了研究斜入射的必要性.Heymsfield[6]分析了二維條件下斜入射SH 波對傾斜基巖自由面位移的放大效應.Stamos等[7]采用一種新的頻域內(nèi)邊界元法研究斜入射體波作用下長大隧道的地震響應.Naggar 等[8]研究了地震波斜入射角度對隧道襯砌彎矩和軸力的影響.從上述文獻可看出，目前對地震動斜入射的研究已取得一些成果，但均未考慮地震動在層狀巖體中傳播時其幅值受巖體阻尼影響而衰減這一情況.

層狀巖體層間動力相互作用屬典型的動接觸問題.在有限元分析中，進行動接觸迭代的數(shù)值模擬方法很多，主要有Lagrange 乘子法[9]、罰函數(shù)法[10]、線性補償法[11]和動接觸力法[12]等.其中，劉晶波等[12]提出的動接觸力法以其計算效率高和穩(wěn)定性好而被廣泛運用，但是該方法忽略了接觸介質(zhì)界面的黏聚力.本文針對地震荷載作用下層狀巖體層間循環(huán)往復相互作用特點，建立一種考慮界面黏結(jié)滑移特性的動接觸力算法，可以反映動力作用下層間非線性接觸特性.

綜上所述，本文建立了一種層狀巖體水工隧洞地震動力響應數(shù)值模擬方法，該方法有效考慮了地震動在層狀巖體中的斜入射特性和層狀巖體層間非線性接觸特性.將本文分析方法應用于巴基斯坦阿扎德帕坦水電站輸水隧洞抗震穩(wěn)定計算，分析軟硬互層狀巖體對隧洞地震響應的影響，以期對復雜層狀巖體中水工隧洞抗震設計進行有益探索.

1 層狀巖體地震動斜入射方法

水工隧洞相比地表、地基和城市地下建筑物的特殊性在于其完全埋置于無限域的山體介質(zhì)之中.在對深埋水工隧洞進行動力時程分析時，由于其所在的工程區(qū)相對整個地震區(qū)域是微小的，遠域地震的波動對工程區(qū)域影響較小，因此可將工程區(qū)外無限域地震區(qū)假設為彈性無限介質(zhì)體，從中截取隧洞所在的有限區(qū)域進行有限化模擬.在時域分析中首先需對計算模型各邊界設置人工邊界，以模擬外行波的透射、內(nèi)行波的入射、邊界處無限域波動場及彈性位移場，且需考慮地表自由面反射波對工程區(qū)的影響.

基于三維黏彈性人工邊界[13]，可將無限域地震波動場問題轉(zhuǎn)化為求解作用于人工邊界節(jié)點的等效節(jié)點力問題，以實現(xiàn)模型內(nèi)外波動場的交互.本文基于波場分解原理，假定計算模型外無限域為均勻彈性介質(zhì)體，地震波為傾斜入射的彈性平面波，將波動場分解為內(nèi)行場和外行場，外行場主要由自模型內(nèi)部向無限域透射的外行波動場構(gòu)成，可由模型內(nèi)顯式有限元逐步積分計算求得，故地下隧洞地震動斜入射實現(xiàn)的關(guān)鍵是求解無限域的斜向內(nèi)行場以及相應的地表反射內(nèi)行場.

對于空間任意入射角度的地震波，可分解為質(zhì)點振動方向與波的傳播方向相一致的P 波（壓縮波）和質(zhì)點振動方向與波的傳播方向相垂直的S 波（剪切波）.下面詳細闡述P 波三維斜入射下地震荷載的計算方法.

1.1 P 波三維斜入射

平面P 波在半空間自由表面經(jīng)反射后會發(fā)生波形轉(zhuǎn)換，產(chǎn)生反射P 波和反射SV 波（如圖1 所示），此時側(cè)向邊界區(qū)內(nèi)的內(nèi)行場位移（t）（i＝1，2，3）和應力（t）（i＝1，2，3）應由入射P 波、反射P 波和反射SV 波各自內(nèi)行場的位移和應力疊加而成；在底部邊界，內(nèi)行場為P 波入射場.即：

假設入射P 波在零時刻的位移時程為u0（t），入射波零時刻波陣面與水平面夾角為入射角α，節(jié)點l（x0，y0，z0）為模型人工邊界上某一節(jié)點，L 為模型底部至自由面高度.需注意本文研究基于入射波波前平行于隧洞軸線，故根據(jù)入射P 波波前和人工邊界節(jié)點l 的幾何空間位置關(guān)系，可求得模型各人工邊界處內(nèi)行場位移時程：

圖1 P 波和SV 波三維斜入射示意圖Fig.1 3D diagram of oblique incidence of plane P and SV waves

式中：β1為反射SV 波在半空間自由面的反射角，β1=arcsin（cssin α/cp），cp和cs分別為P 波和SV 波的波速；A1和A2分別為反射P 波和反射SV 波的幅值放大系數(shù)，其值參考文獻[14]取得；η（·）為考慮巖體阻尼情況下，地震動的幅值沿傳播距離的衰減系數(shù)，胡進軍等[15]研究表明地震動在復雜巖層中傳播時可認為其幅值受巖體阻尼的影響而呈線性衰減，張志國等[16]進一步將其表述為：

式中：η0是一與圍巖阻尼和彈性模量相關(guān)的經(jīng)驗常數(shù)，而與地震動的輸入方式無關(guān)，張志國等[16]對各因素對η0值的影響做了大量數(shù)值試驗，可參考取值；Δli為地震波i 傳播至邊界節(jié)點l 所經(jīng)過的路程；Δt為入射P 波從零時刻波陣面?zhèn)鞑サ竭吔绻?jié)點l 所需的時間；Δl1～Δl3、Δt1～Δt3可表示為：

側(cè)向人工邊界節(jié)點A 處的應力為：

底部人工邊界節(jié)點A 處的應力為：

式中：λ、G 為介質(zhì)的拉梅常數(shù)和剪切模量.由于各人工邊界上y 方向的位移時程為零，故相應的y 方向應力也為零.

由式（2）～（6），可求解出P 波三維斜入射條件下人工邊界處內(nèi)行場的位移時程（t）、速度時程（t）及其相應的應力時程（t），以此實現(xiàn)三維斜入射條件下P 波荷載輸入：

式中：Kli和Cli分別表示節(jié)點l（x0，y0，z0）在方向i 上的人工邊界彈簧和阻尼系數(shù)，可參考文獻[14]取值；Al為人工邊界上節(jié)點l 的控制面積.

SV 波和SH 波三維斜入射的計算公式可根據(jù)上述P 波斜入射條件下地震荷載的計算公式進行類似推導，不再贅述.由此本文建立了一種與三維黏彈性人工邊界相適應的層狀巖體中地震動三維斜入射的輸入方法.

1.2 算例驗證

為驗證本文層狀巖體中地震動斜入射輸入方法的合理性，建立了一個有限元模型來分析斜入射條件下半無限域三維彈性介質(zhì)體的動力響應問題.三維有限元模型的尺寸為800 m×800 m×800 m，共計46 656 個六面體單元，如圖2 所示.模型彈性模量為10 GPa，泊松比為0.3，密度為2 000 kg/m3.計算模型底部和側(cè)向邊界施加三維黏彈性人工邊界，頂部為自由面，取頂部自由面中點A（400，400，800）為監(jiān)測點.入射P 波的位移時程曲線如圖3 所示.

圖2 三維有限元模型Fig.2 3D finite element model

圖3 入射P 波位移時程Fig.3 Displacement time-history of incident P waves

圖4 所示為當平面P 波入射角分別為15°和30°時監(jiān)測點A 的豎直向（z 向）位移時程的理論解與數(shù)值解.需要說明的是本文入射角α 是指入射波零時刻波陣面與水平面夾角，故P 波入射角分別為15°和30°時波陣面的法向向量分別為（0.259，0，0.966）和（0.5，0，0.866）.從圖中可看出不同入射角下監(jiān)測點A 的豎直向位移時程的計算結(jié)果與理論值比較符合.

圖4 監(jiān)測點A 的豎直向位移時程Fig.4 Displacement time-history of the observation point A

2 層狀巖體層間動接觸系統(tǒng)分析模型

在構(gòu)造應力和地震荷載聯(lián)合作用下層狀巖體層間易發(fā)生相對錯動，從而對隧洞結(jié)構(gòu)造成嚴重的損傷破壞.根據(jù)巴基斯坦阿扎德帕坦水電站現(xiàn)場觀測資料，工程區(qū)層間剪切、擠壓破碎帶較為常見.本文針對層狀巖體層間循環(huán)往復作用特點，建立一種考慮界面黏結(jié)與滑移特性的層狀巖體接觸系統(tǒng)動力響應分析模型，用來模擬層間非線性滑移破壞.

2.1 動接觸力法的基本方程

對層狀巖體層間接觸系統(tǒng)進行有限元離散后，可得到包含動接觸力的節(jié)點運動方程：

式中：M、C 和K 分別表示層間接觸節(jié)點的質(zhì)量、阻尼和剛度矩陣；、和u 分別表示接觸節(jié)點的加速度、速度和位移向量；F 為接觸節(jié)點的外荷載向量；R為接觸節(jié)點的動接觸力向量，且R=N+T，N 和T分別為R 的法向和切向分量.

采用中心差分法求解接觸節(jié)點的運動方程，相應的時域積分格式為：

聯(lián)立式（8）、式（9）和式（10），可得：

由式（11）～（13）可看出，由t 時刻接觸節(jié)點的運動狀態(tài)和動接觸力，可求解出t+Δt 時刻接觸節(jié)點的運動狀態(tài).t 時刻接觸節(jié)點的運動狀態(tài)是已知的，而動接觸力Rt是未知量，需根據(jù)t～t+Δt 時刻的接觸狀態(tài)基于相應的接觸條件計算.

2.2 考慮界面黏結(jié)滑移特性的動接觸力算法

假定在地震作用前層狀巖體層間接觸良好，考慮界面黏聚力，則接觸點對處于黏結(jié)接觸狀態(tài)，如圖5 所示.在強震過程中，層間可能會發(fā)生相對滑動，導致接觸節(jié)點與其對應的單元某一面發(fā)生接觸，接觸面進入滑動接觸狀態(tài)，此時不再考慮黏聚力.

圖5 接觸模型和點對上動接觸力示意圖Fig.5 Contact model and dynamic contact force on node pairs

假定接觸節(jié)點對在地震作用下沒有相對滑動，則在t+Δt 時刻接觸點對仍可認為處于黏結(jié)接觸狀態(tài).此時相應的接觸邊界上的節(jié)點對i 和i′須滿足接觸邊界條件，包括位移接觸條件和接觸力邊界條件：

式中：ni為接觸節(jié)點對的單位法向矢量，由節(jié)點i′指向i；ti為單位切向矢量.

聯(lián)立式（11）、式（13）和式（14），可得：

式中：Mi和Mi′分別為節(jié)點i 和i′的集中質(zhì)量.

顯然，式（15）中法向和切向動接觸力是在層狀巖體層間處于黏結(jié)狀態(tài)下求得的，而層間接觸面是水工隧洞中薄弱部分，其在地震循環(huán)作用下的損傷是不可忽略的.實際上，在動力作用下，層狀巖體層間接觸面存在黏結(jié)接觸、滑動接觸和分離等多種接觸狀態(tài)，因此在每一時步計算完畢后，需要對接觸節(jié)點對的接觸狀態(tài)進行判別，并對動接觸力進行修正[17].層間接觸面的破壞形式主要包括沿切向的剪切滑移和沿法向的張拉開裂，具體方法如下：

1）若Δ1i<0 且‖‖＞cAi，表明接觸面受拉，并且發(fā)生了張拉開裂破壞，接觸點對i 和i′處于分離狀態(tài).此時，動接觸力修正如下：

2）若Δ1i<0 且‖‖≤cAi，或者Δ1i>0 且‖‖≤μs‖‖，表明接觸節(jié)點對處于黏結(jié)接觸狀態(tài)，此時無需對和進行修正.

3）若Δ1i>0 且‖‖≥μs‖‖，表明接觸面受壓，且接觸點對處于滑動接觸狀態(tài)，即接觸面發(fā)生剪切滑移破壞.此時，須校核切向動接觸力，如下所示：

式中：μs和μd分別為接觸面靜、動摩擦因數(shù)；c 為接觸面的黏聚力.若接觸點對一直處于黏結(jié)接觸狀態(tài)，則c>0.若接觸點對從黏結(jié)接觸轉(zhuǎn)變到滑動接觸或分離狀態(tài)，則c=0；Ai為接觸節(jié)點i 的控制面積.由式（15）～（17）即可求得t 時刻動接觸力Rt，將Rt代入式（13），計算t+Δt 時刻附加位移場Δut+Δt，更新接觸節(jié)點的總位移ut+Δt.

3 工程實例

3.1 工程概況和計算模型

阿扎德帕坦水電站位于巴基斯坦的Jhelum 河上，為該河段水電開發(fā)中的一級，以發(fā)電為主.引水發(fā)電系統(tǒng)位于河流左岸，導流隧洞布置在右岸.導流隧洞工程區(qū)地質(zhì)條件復雜，隧洞所穿越的基巖巖性為砂巖與非砂巖類呈互層狀分布，地層主要為單斜構(gòu)造，巖層傾角∠71°～∠81°，層間剪切、擠壓破碎帶較為常見.

阿扎德帕坦工程區(qū)屬于地震活動區(qū)，主要受印度板塊持續(xù)向歐亞板塊俯沖運動影響.巴基斯坦北部和阿扎德地區(qū)為地震強烈活動地區(qū)，受多個地震板塊構(gòu)造影響，地震活動多由本地區(qū)活動斷層運動引起，阿扎德帕坦水電站工程區(qū)就在這一區(qū)域.時間最近的大地震為2005 年10 月8 日的7.6 級地震.根據(jù)中國地震局地質(zhì)研究所研究成果，阿扎德帕坦水電站工程區(qū)50 年超越概率10%（DBE）的峰值加速度為0.315g，對應的地震基本烈度為Ⅷ度.2#導流隧洞洞長667 m，洞身段埋深在140 m 左右.采用圓形斷面結(jié)構(gòu)，開挖洞徑10.0 m，圍巖以Ⅲ、Ⅳ類為主，襯砌采用C25 鋼筋混凝土結(jié)構(gòu)，厚度為60 cm.

選取如圖6 所示洞身段含泥、砂巖互層的隧洞區(qū)域建立水工隧洞三維有限元模型，考慮到若模型建至地表，則單元數(shù)量過多，動力計算耗時將呈指數(shù)級增長.為了提高計算效率，隧洞頂部取50 m.層狀巖層走向與洞軸線垂直，傾角為70°.模型共剖分了48 608 個八節(jié)點六面體單元和52 041 個節(jié)點，其中混凝土襯砌單元3 888 個.模型范圍及坐標系：x 方向從-60.0 m 到60.0 m，與洞軸線垂直；y 方向從-85.0 m 到85.0 m，與洞軸線重合，順水流為正；z 方向從400.0 m 到520.0 m，與大地坐標系平行，豎直向上為正.

圖6 層狀巖體模型Fig.6 Interbed layered rock mass model

三維初始地應力場根據(jù)設計院提供的實測地應力反演分析得到，側(cè)壓力系數(shù)取kx=1.1，ky=0.85，kz=1.0.層狀巖體砂巖、泥巖、接觸面和襯砌的材料力學參數(shù)取值見表1.圍巖臨界阻尼比取5%，則參考文獻[16]砂巖層的地震動幅值衰減系數(shù)取0.04%，非砂巖層的地震動幅值衰減系數(shù)取0.07%.動力計算之前，采用三維彈塑性損傷有限元法進行隧洞的靜力開挖與支護計算，其相應的計算結(jié)果作為動力計算的初始條件.

表1 材料力學參數(shù)Tab.1 Physico-mechanical parameters of materials

3.2 計算條件

計算程序采用課題組自主開發(fā)的大型地下洞群抗震穩(wěn)定動力時程分析平臺[18]，并將本文地震動輸入方法和動接觸力算法嵌入其中.接觸面模擬的部位分別位于非砂巖層（即軟巖）與砂巖層（即硬巖）相交界面處.動力加載前，首先基于節(jié)點分離技術(shù)[19]，通過增加砂巖層與非砂巖層接觸面兩側(cè)的共用節(jié)點，并設置一一對應的接觸節(jié)點對，以完成軟巖與硬巖單元的節(jié)點分離.圍巖和襯砌采用基于M-C 屈服準則的動力彈塑性損傷本構(gòu)模型[18]，三維損傷演化方程如下所示：

式中：R 為材料的損傷常數(shù)；Di、分別為第i 個主應變方向上的損傷系數(shù)和塑性偏應變，為便于后處理顯示材料的損傷狀態(tài)，采用標量形式的損傷系數(shù)D=.

模型的底部、四周和頂部均采用黏彈性人工邊界，以吸收斜入射條件下的地震波及其在地表自由面的反射波.地震波采用美國強震記錄的El-Centro波，并根據(jù)阿扎德帕坦水電站工程區(qū)抗震設防烈度，將峰值加速度調(diào)整為3.15 m/s2，截取其中變化劇烈、幅值較大的20 s 時段作為入射波，經(jīng)濾波和基線校正處理后加速度時程曲線如圖7 所示.計算時考慮斜入射（∠30°入射角）和豎直入射的地震動對圍巖和襯砌的作用.動力計算同時考慮P 波和SV 波對水工隧洞的作用.其中，SV 波采用如圖7 所示的入射波，P 波加速度時程取為SV 波的2/3[20].

圖7 入射波加速度時程曲線Fig.7 Time-history curve of input wave acceleration

選取軟巖中間斷面為監(jiān)測斷面，布置如圖8 所示的監(jiān)測方案，監(jiān)測點A、B、C 分別位于監(jiān)測斷面上襯砌的頂拱、左腰部和底部幾個關(guān)鍵部位，用以監(jiān)測地震加載過程中襯砌的位移和應力等指標特性，監(jiān)測點D、E 分別位于層狀巖體硬巖與軟巖接觸面兩側(cè)，用以監(jiān)測層狀巖體層間相對運動特征.動力計算分3 種工況：①地震動豎直入射，不考慮動接觸；②地震動斜入射，不考慮動接觸；③地震動斜入射，考慮動接觸.需要注意的是，工況①地震動豎直入射是指地震動自模型底部垂直入射，水平面內(nèi)沿垂直水流向（x 向）振動，豎直向振動取水平向振動的2/3.地震動斜入射的入射方向矢量為（0.5，0，0.866）.

圖8 監(jiān)測點布置Fig.8 Layout of the monitoring points

3.3 計算結(jié)果及分析

3.3.1 圍巖破壞區(qū)分布

在地震循環(huán)荷載作用下，水工隧洞洞周圍巖循環(huán)加卸載，圍巖應力一直處于波動狀態(tài)，且塑性變形不斷累積，使得圍巖總應變逐漸增加，導致圍巖總破壞區(qū)體積相應增大.3 種工況下震后洞周圍巖破壞區(qū)分布如圖9 所示.

圖9 震后圍巖破壞區(qū)分布Fig.9 Distribution of surrounding rock failure zone after earthquake

當?shù)卣饎迂Q直入射時，洞周圍巖破壞區(qū)分布較少，從橫向上看，洞周塑性破壞區(qū)在腰拱處有逐漸向深部擴展的趨勢，且塑性區(qū)深度為2.52 m，開裂區(qū)分布較少，深度為0.97 m；從縱向上看，開裂區(qū)主要分布在軟巖穿過的區(qū)域，其余處較少，僅在腰拱處出現(xiàn).考慮地震動斜入射時，軟巖及層間破碎帶處的塑性破壞區(qū)顯著變大，而開裂區(qū)增大不明顯，塑性區(qū)深度達4.96 m，開裂區(qū)深度為1.65 m，表明考慮地震動在層狀巖體中的傳播特性后，斜入射地震動極大削弱了洞周圍巖的穩(wěn)定性，有可能導致隧洞結(jié)構(gòu)的局部失穩(wěn)和損傷破壞.考慮地震動的斜入射和層間動接觸相互作用時，層間破壞區(qū)進一步擴展，主要表現(xiàn)為塑性區(qū)和開裂區(qū)均明顯增大.與僅考慮地震動斜入射工況相比，塑性區(qū)深度增加了2.15 m，開裂區(qū)深度增加了0.86 m，且軟巖穿越的洞周幾乎被開裂區(qū)包圍.從圖9 中還可看出震后開裂區(qū)主要分布在隧洞的腰拱及其上部，這與Wang 等[21]的研究結(jié)果相符.

3.3.2 襯砌結(jié)構(gòu)位移時程分析

3 種工況下襯砌結(jié)構(gòu)不同監(jiān)測點處合位移時程如圖10 所示.由圖10 可看出3 種工況下，1）頂拱、腰部和底部位移時程曲線的波形和波動規(guī)律基本一致，均出現(xiàn)了多個明顯的波峰.監(jiān)測點位移同時出現(xiàn)波峰和波谷，表明水工隧洞襯砌結(jié)構(gòu)各部位處于同步震動狀態(tài)；2）在0～5 s 時間段內(nèi)，襯砌結(jié)構(gòu)各部位位移時程曲線呈現(xiàn)大幅度波動，腰部最大位移要比頂拱和底部大.本文采用腰部與頂拱的合位移差值來表征相對位移，相對位移可以更好地表征隧洞襯砌結(jié)構(gòu)的變形特征，圖11 所示為3 種工況下腰部及頂拱位移動力響應差值時程曲線.

工況①下，襯砌結(jié)構(gòu)各監(jiān)測點處最大位移為7.4 cm，腰部與頂拱最大相對位移為0.67 cm，發(fā)生在4.95 s，但震后腰、拱相對位移為0.21 cm，隧洞襯砌結(jié)構(gòu)相對變形量值較小.

圖10 監(jiān)測點合位移時程曲線Fig.10 Displacement time history curves of lining

圖11 腰部與頂拱相對位移時程曲線Fig.11 Relative displacement time history curve of the haunch and top arch

工況②下，當考慮地震動的空間斜入射特性時，各監(jiān)測點位移及腰、拱相對位移時程曲線相比工況①的差別主要表現(xiàn)在波動幅值上，各監(jiān)測點處最大位移為8.7 cm，腰、拱最大相對位移為0.85 cm，發(fā)生在5.5 s，震后相對位移為0.52 cm，表明襯砌變形受地震動入射角影響較大，地震動的三維斜入射特性和入射邊界的非一致特性對隧洞襯砌結(jié)構(gòu)位移響應影響較大.這主要是因為考慮地震動斜入射時，地震波場與豎直入射時具有明顯差異.地震動斜入射時，入射P 波和SV 波會在自由面發(fā)生波形轉(zhuǎn)換，各自分別形成反射P 波和反射SV 波，故模型人工邊界處波場是由不同入射波和反射波疊加而成的，使得人工邊界上各節(jié)點具有不同的振動波形，產(chǎn)生了放大效應.而地震動從模型底部豎直入射時，經(jīng)地表自由面反射后地震波仍具有相同的反射方向和振動幅值，模型人工邊界上節(jié)點也具有相同的振動方向和振動波形.因此，斜入射時襯砌結(jié)構(gòu)各部位位移響應相比豎直入射時要大.

工況③下，當?shù)卣饎有比肷淝铱紤]層間動力相互作用后，各監(jiān)測點處最大位移為9.8 cm，腰、拱最大相對位移達到1.41 cm，發(fā)生在5.15 s，在前期地震動波動較為劇烈的0～6 s 時間段內(nèi)，其腰、拱相對位移時程曲線波動幅度相比前2 個工況較大，波動范圍為-0.55～1.41 cm，震后相對位移為1.04 cm.從理論上說，襯砌結(jié)構(gòu)的相對位移在一次地震結(jié)束后應當回到0[3].然而，工況③中襯砌腰部和頂拱的相對位移值在震后為1.04 cm，表明隧洞結(jié)構(gòu)在地震荷載作用下發(fā)生永久變形.可見在構(gòu)造應力和地震荷載聯(lián)合作用下層狀巖體層間易發(fā)生剪切滑移破壞，且斜入射地震動加劇了襯砌腰部的相對變形，致使襯砌腰部抗剪段安全問題突出.

3.3.3 層狀巖體層間相對運動分析

在層狀巖體層間接觸系統(tǒng)地震動響應過程中，硬巖與軟巖層間循環(huán)往復作用，發(fā)生了復雜的動接觸行為，如：黏結(jié)接觸、滑動接觸和分離等多種接觸狀態(tài)，進而產(chǎn)生了層間錯動位移，對隧洞結(jié)構(gòu)造成嚴重的損傷破壞.3 種工況下層狀巖體層間相對位移時程曲線如圖12 所示，進一步說明了地震動斜入射對層狀巖體水工隧洞接觸響應的影響.

圖12 層間接觸面相對位移時程曲線Fig.12 Relative displacement time history curve of interface

當?shù)卣饎迂Q直入射時，層狀巖體層間相對位移在0 線附近波動變化，層間最大相對位移為1.77 cm，在0～6 s 時間段內(nèi)，層間相對位移波動幅度較大，主要在-1.60～1.77 cm 范圍內(nèi)上下波動，后期波動幅度逐漸減小，震后基本為0.當考慮地震動斜入射后，在0～6 s 時間段內(nèi)，層間相對位移在0 線上下劇烈波動，波動范圍為-1.96～2.16 cm，層間最大相對位移為2.16 cm，發(fā)生在5.75 s，受斜入射下地震動輸入非一致性的影響，后期波動幅度雖較0～6 s 時間段有所減小，但相比工況①仍較大，震后逐漸減小到0 線附近.當考慮層間動接觸力后，在前期0～1.5 s 時間段內(nèi)，層間相對位移在0 線附近波動極小，表明此時間段內(nèi)層狀巖體層間接觸面接觸良好，處于黏結(jié)接觸狀態(tài)，或從黏結(jié)接觸狀態(tài)向滑動接觸或分離接觸狀態(tài)轉(zhuǎn)變的過渡階段；在1.5～6 s 時間段內(nèi)，層狀巖體層間相對位移波動較為明顯，且出現(xiàn)了明顯的錯動位移，最大錯動位移為-3.91 cm，發(fā)生在5.2 s；后期層間相對位移主要在-3.00～-4.00 cm 范圍內(nèi)上下波動，震后為-3.94 cm，表明考慮層間動接觸力后，層狀巖體層間發(fā)生了明顯的剪切滑移破壞，且受地震動斜入射的影響，層間相對位移在-3.50 cm 上下波動明顯.

3.3.4 襯砌結(jié)構(gòu)應力時程分析

由于混凝土的抗壓強度遠大于其抗拉強度，地震作用下水工隧洞襯砌結(jié)構(gòu)的損傷破壞主要是拉裂破壞.因此本文主要分析地震作用下襯砌結(jié)構(gòu)的最大主應力變化規(guī)律.3 種工況下襯砌結(jié)構(gòu)不同監(jiān)測點處最大主應力時程曲線如圖13 所示.

圖13 襯砌監(jiān)測點最大主應力時程曲線Fig.13 Maximum principle stress time history curves of lining

當?shù)卣饎迂Q直入射時，在前5 s 內(nèi)，頂拱、腰部和底部的最大主應力變化劇烈，波動范圍主要為-0.50～1.47 MPa，在峰值處腰部、頂拱和底部的拉應力量值分別達到1.47 MPa、1.20 MPa、0.96 MPa，腰部的拉應力量值超過了混凝土的抗拉強度，且腰部的應力水平明顯大于頂拱和底部的應力水平.

當?shù)卣饎有比肷鋾r，受自由面反射波場疊加效應的影響，整體上看隧洞結(jié)構(gòu)應力水平的地震反應大于豎直入射時的地震反應.斜入射時襯砌腰部最大主應力相比頂拱和底部一直波動較為劇烈，在峰值處襯砌腰部、頂拱和底部的拉應力量值分別為1.61 MPa、1.27 MPa、1.23 MPa，其中腰部和頂拱最大拉應力達到了混凝土的抗拉強度，表明襯砌的腰部為襯砌結(jié)構(gòu)受力的不利部位.

當?shù)卣饎有比肷淝铱紤]層間動接觸力后，斜入射的地震動加劇了層狀巖體硬巖和軟巖間循環(huán)往復相互作用，腰部最大主應力波動范圍增大，主要在0.5～1.81 MPa，在峰值處襯砌腰部、頂拱和底部的拉應力量值分別為1.90 MPa、1.69 MPa、1.45 MPa，均超過了混凝土的抗拉強度.表明在考慮層間動接觸力后，腰部的應力響應更為明顯，層狀巖體層間相對滑移對襯砌結(jié)構(gòu)受力具有重要影響，加劇了隧洞腰部拉應力的變化，使得襯砌結(jié)構(gòu)發(fā)生拉裂損傷破壞.

3.3.5 襯砌結(jié)構(gòu)損傷分析

圖14 所示為考慮地震動斜入射和層間動接觸相互作用后，襯砌結(jié)構(gòu)震后損傷系數(shù)分布圖.由圖可知，在地震動輸入完成后，襯砌損傷區(qū)主要分布在層狀巖體層間接觸部位和軟巖穿過部位，且向兩側(cè)延伸約5 m 范圍內(nèi)，最大損傷系數(shù)接近于1，且主要位于襯砌結(jié)構(gòu)腰部.表明在進行層狀巖體水工隧洞抗震設計時，需要采用抗斷技術(shù)限制襯砌結(jié)構(gòu)發(fā)生開裂破壞.

圖14 襯砌結(jié)構(gòu)損傷系數(shù)分布Fig.14 Damage coefficient distribution of lining structure

4 結(jié)論

基于地震動三維空間斜入射輸入方法和考慮層狀巖體層間黏結(jié)滑移特性的動接觸力算法，建立了一種層狀巖體水工隧洞地震動力響應分析模型.結(jié)合阿扎德帕坦水電站水工隧洞工程實例，對地震作用下襯砌結(jié)構(gòu)動力響應及損傷破壞進行分析，得到如下結(jié)論：

1）層狀巖體中軟巖穿過部位的破碎帶加劇了隧洞結(jié)構(gòu)的地震響應，表現(xiàn)為圍巖塑性破壞區(qū)及開裂區(qū)的顯著增大.斜入射地震動削弱了水工隧洞襯砌結(jié)構(gòu)的穩(wěn)定性，在考慮層間動力相互作用后接觸面附近破壞區(qū)進一步發(fā)展且軟巖穿越的洞周幾乎被開裂區(qū)包圍.

2）地震作用下隧洞襯砌結(jié)構(gòu)不同部位監(jiān)測點的位移時程曲線與輸入地震動位移時程曲線相類似，表現(xiàn)為一種同步震動趨勢.對比地震動斜入射和豎直入射，隧洞襯砌結(jié)構(gòu)的變形響應受地震動入射角影響較大.當同時考慮地震動斜入射和層間動接觸力時，層狀巖體層間易表現(xiàn)為剪切滑移破壞.

3）與地震動豎直入射相比，斜入射時隧洞襯砌結(jié)構(gòu)的應力響應更大.當考慮層間動接觸作用后，接觸面附近處襯砌應力進一步加大，襯砌腰部應力響應相比頂拱較大，因而會首先發(fā)生開裂破壞，故腰部是水工隧洞襯砌結(jié)構(gòu)抗震設計的薄弱部位.

4）地震作用下層狀巖體層間震動不同步，極易發(fā)生相互錯動，進而使襯砌損傷破壞.襯砌損傷區(qū)主要分布于層間接觸部位和軟巖穿過部位，且向兩側(cè)延伸約5 m.