深埋互层岩体地下洞室地震响应数值模拟研究

赵勐; 肖明; 陈俊涛; 杨步云

doi:10.11779/CJGE202112002

深埋互层岩体地下洞室地震响应数值模拟研究 English Version

赵勐^{1, 2,},
肖明^{1, 2, ,},
陈俊涛^{1, 2},
杨步云^{1, 2}

1.
武汉大学水资源与水电工程科学国家重点实验室,湖北　武汉　430072
2.
武汉大学水工岩石力学教育部重点实验室,湖北　武汉　430072

基金项目:

国家自然科学基金项目 52079097

国家重点基础研究发展计划（“973”计划）项目 2015CB057904

国家自然科学基金项目 51579191

详细信息

作者简介:
赵勐(1994— ),男,博士研究生,主要从事地下结构稳定数值分析研究。E-mail：zhaomeng@whu.edu.cn

通讯作者:
肖明, E-mail：mxiao@whu.edu.cn

中图分类号: TU45
计量
- 文章访问数: 249
- HTML全文浏览量: 49
- PDF下载量: 351
出版历程
- 收稿日期: 2021-03-02
- 网络出版日期: 2022-11-30
- 刊出日期: 2021-11-30

Numerical simulation of seismic response of a deeply-buried underground cavern in interbedded rock mass

ZHAO Meng^{1, 2,},
XIAO Ming^{1, 2, ,},
CHEN Jun-tao^{1, 2},
YANG Bu-yun^{1, 2}

1.
State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
2.
Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education, Wuhan University, Wuhan 430072, China

摘要

摘要: 针对深埋大型地下洞室地震波场特性,考虑近场斜入射地震动的方向性、多面性和非一致性,通过将场地地震反应转化为人工边界上等效荷载实现了深埋地下洞室地震波斜入射。针对地震作用下互层岩体层间的动力相互作用特点,建立考虑层面震动劣化效应和黏结滑移特性的动接触力算法。由此构建地震动斜入射下深埋互层岩体地下洞室地震响应分析方法,将该方法应用于阿扎德帕坦水电站地下厂房地震损伤演化分析中,研究结果表明：斜入射地震动加剧了衬砌结构的位移和应力响应,主要体现在波动幅值上,厂房上部边墙和顶拱损伤破坏程度最大；考虑动接触后,层面附近洞室的地震响应增大,岩层间产生明显的地震劣化现象和剪切滑移破坏,层间错动更加明显,最大错动位移在5.9 cm处趋于稳定；并从横向和纵轴向两个角度归纳总结了互层岩体地下洞室结构的震损特征和破坏模式。
- 地震动斜入射 /
- 互层岩体 /
- 地下洞室 /
- 震动劣化 /
- 地震响应 /
- 数值分析
Abstract: In order to study the seismic wave field characteristics of large deeply-buried underground caverns, the input methods for obliquely incident earthquakes in deeply-buried underground caverns are proposed by converting the site seismic response into the equivalent load acting on the artificial boundaries. It is suitable to reflect the incident direction, the multi-incident surfaces and the inconsistency of near-field obliquely incident earthquakes. According to the dynamic interaction characteristics between interlayers in interbedded rock mass under seismic action, a dynamic contact analysis method considering the seismic deterioration effect and the bond-slip characteristics of interface is also established. Consequently, a dynamic response analysis method for a deeply-buried underground cavern in interbedded rock mass under obliquely incident seismic waves is constructed and applied to the seismic damage evolution process analysis of the underground powerhouse of Azad Pattan Hydropower Station. The simulated results reveal that the obliquely incident earthquake aggravates the seismic reaction of lining structure, which mainly lies in the amplitudes of the displacement and stress fluctuations. Some areas at the upper side wall and the arch of the linings in the main powerhouse are damaged to the most severe degree. After considering the dynamic contact, the seismic response of caverns near interface further increases, and the evident seismic deterioration effect and shear slip failure occur at interface. The dislocation between interlayers is more obvious, and the maximum dislocation displacement tends to be stable at 5.9 cm. Besides, the seismic damage characteristics and failure modes of the lining structures of the underground powerhouse in interbedded rock mass are discussed from the transverse and longitudinal angles.
- obliquely incident earthquake /
- interbedded rock mass /
- underground cavern /
- vibration deterioration /
- seismic response /
- numerical simulation

HTML全文

0. 引言

随着“一带一路”合作深入发展,为加快巴基斯坦水能资源的开发和缓解其电网严重缺电的局面,中国投资帮助巴基斯坦规划兴建了众多大型水电站,进而形成了大量的水电站地下洞室群。这些水电站工程区位于喜马拉雅山脉西段南麓,毗邻中国西南高山峡谷地区。区内强震活动频发,且多种地层岩性分布广泛,软硬互层岩体或者泥化软弱夹层发育明显,工程地质条件极其复杂^[1]。此外,西南高地震烈度区地下洞室一般埋深较大,故研究深埋互层岩体地下洞室地震响应特性和破坏机理对工程建设安全具有重要意义。

深埋互层岩体地下洞室地震响应时程分析主要包括两方面的内容：①深埋近场地震动输入方法；②互层岩体岩层间动力接触分析。在已有的针对地下厂房洞室群的动力时程分析中,多是假定地震动从模型底部竖直向上入射,且模型顶部不作处理,直接建至地表^[2]。这对浅埋地下洞室群的动力时程计算是可行的,但是对于深埋近场地震,即震源距工程场地较近时,一方面有限元单元数量成倍增加,计算模型规模过于庞大,动力计算效率低下；另一方面地震动可能从任意方向传播至地下洞室群工程区域^[3]。根据近年来强震动观测记录的统计,发现基岩场地的地震动入射角平均在60°,从而引起地下结构的非一致性变形^[4]。杜修力等^[5]研究了地震动斜入射下岩体隧道的时域地震反应规律,研究成果表明斜入射下隧道结构的动力响应远大于竖直入射；赵宝友等^[6]研究了不同P波入射角度下大型地下岩体洞室群地震反应,研究发现斜入射地震动会对洞室群的稳定性产生更大的不利影响；赵密等^[7]和黄景琦等^[8]分别研究地震动斜入射对跨断层隧道和地铁车站结构地震响应的影响；张志国等^[9]考虑地震波入射的非一致性,提出一种适用于地下洞室群围岩动力时程分析的地震波空间斜入射法。Kouretzis等^[10]基于三维薄壳应变分析法研究斜入射S波作用下长隧道的地震响应。总体看来,目前的研究重点是斜入射地震动下浅埋地下结构的地震响应,但在中国西南高山峡谷地区多为深埋地下洞室群,关于深埋近场斜入射地震动对地下洞室结构地震响应的研究相对较少。

地震作用下岩层间动力循环相互作用具有复杂的接触非线性特性,并且伴随层面的震动劣化过程,因此一直是互层岩体地下洞室地震响应分析的难点问题。刘晶波等^[11]通过引入动接触界面的约束条件,基于显式中心差分法求解动接触力,该方法步骤清晰,计算效率高且收敛性较好,被广泛应用于复杂接触系统的动接触计算。但是该法没有考虑接触介质之间的黏聚力和地震循环荷载对接触面震动劣化效应的影响。Lee等^[12]、刘博等^[13]、倪卫达等^[14]通过大量循环剪切试验发现了岩石节理的峰值抗剪强度在循环剪切过程中的劣化效应。

综合上述分析,本文首先基于波场分解原理和三维黏弹性人工边界条件^[15],建立1种深埋地下洞室近场地震动斜入射的输入方法；针对互层岩体层间的动力相互作用特点,建立1种考虑接触面震动劣化效应和多种接触状态的动接触力算法。将本文分析方法应用于巴基斯坦阿扎德帕坦水电站地下厂房抗震稳定计算中,分析地震动斜入射下互层岩体对地下厂房地震响应的影响,以期对地下洞室的抗震减震设计提供有益的参考。

1. 深埋地下洞室地震动斜入射方法

1.1 计算原理

针对动力时程分析中深埋地下洞室群有限元模型不便于建立到山顶地表的现实情况,建立1种考虑半无限空间自由波动场传播特性的斜入射地震动输入方法。采用该方法,对于深埋地下洞室,模型可以不必建立至地表,在模型顶部进行人工截断,设置黏弹性人工边界^[15]。同时采用解析法求解模型顶部至地表自由面空间范围内的地震波动场,并通过黏弹性人工边界将求解的反射波场入射到有限元模型内部。模型的四周边界采用波场分解技术,将总波动场分解为外行场和内行场,外行场为自模型内部通过人工边界向无限域透射的波场,可通过显式有限元逐步积分计算求解。故深埋地下洞室地震动斜入射实现的关键是求解无限域的斜向内行场以及相应的地表反射内行场,并通过人工边界将无限域地震波动场问题转化为求解作用于边界节点的等效节点力问题,以实现模型内外波动场的波场交换。

1.2 深埋地下洞室近场地震动斜入射输入方法

深埋近场地震动的斜入射与传统的竖直入射区别在于3方面：①斜入射地震波会在半空间自由表面发生复杂的波形转换现象,如图1所示,入射P波会经自由面反射分别形成反射P波和反射SV波,故模型人工边界处内行场位移和应力则是由入射P波、反射P波和SV波各自内行场的位移和应力叠加而成,SV波斜入射时也具有相同的现象^[16]；②深埋地下洞室在地震动斜入射时,其模型人工边界具有多个迎波面,每个迎波面上不同节点到震源的距离各不相同,因此在同一时刻人工边界上每个节点接受到的入射地震波会存在相位差；③大型地下洞室的顶拱和高边墙由于具有较大的临空面,斜入射地震波会在此区域发生复杂的反射和衍射现象,产生放大效应。假定计算模型人工边界以外的无限域为均匀弹性介质,地震波为弹性平面波,地表自由面为水平面,则计算模型各人工边界上节点的内行场位移时程可推导如下：

图 1 深埋地下洞室地震动斜入射示意图

Figure 1. 3D diagram of oblique incident earthquake for a deeply-buried underground cavern

下载: 全尺寸图片幻灯片

对于斜入射P波：

$u_{l i}^{R} (t) = {\begin{cases} u_{P_{2}}^{R} (t) + u_{{SV}_{2}}^{R} (t) (顶部边界) \\ u_{P_{1}}^{R} (t) + u_{P_{2}}^{R} (t) + u_{{SV}_{2}}^{R} (t) (侧向边界) \\ u_{P_{1}}^{R} (t) (底部边界) \end{cases}$ 。

(1)

对于斜入射SV波：

$u_{l i}^{R} (t) = {\begin{cases} u_{{SV}_{2}}^{R} (t) + u_{P_{2}}^{R} (t) (顶部边界) \\ u_{{SV}_{1}}^{R} (t) + u_{{SV}_{2}}^{R} (t) + u_{P_{2}}^{R} (t) (侧向边界) \\ u_{{SV}_{1}}^{R} (t) (底部边界) \end{cases}$ 。

(2)

式中, $u_{P_{1}}^{R} (t)$ , $u_{{SV}_{1}}^{R} (t)$ , $u_{P_{2}}^{R} (t)$ , $u_{{SV}_{2}}^{R} (t)$ 为入射P波、入射SV波、反射P波和反射SV波的位移时程。

假设入射P波和SV波在零时刻的位移时程分别为 $u_{P 0} (t)$ , $u_{SV0} (t)$ ,入射波波前与水平面的夹角定义为入射角 $α$ ,且入射波波前平行于洞室的纵轴线。有限元模型高度为L,模型顶部至地表深度为H, $l (x_{l}, y_{l}, z_{l})$ 为模型人工边界上某一节点。故根据入射波波前和人工边界节点l的空间位置关系,采用行波公式可求得模型各人工边界处内行场位移时程：

对于斜入射P波：

$\begin{array}{l} u_{P 1}^{R} (t) = u_{P 0} (t - Δ t_{1}) {[\begin{matrix} \sin α & 0 & \cos α \end{matrix}]}^{T} ， \\ u_{P 2}^{R} (t) = A_{1} u_{P 0} (t - Δ t_{2}) {[\begin{matrix} - \sin α & 0 & \cos α \end{matrix}]}^{T} ， \\ u_{SV 2}^{R} (t) = A_{2} u_{P 0} (t - Δ t_{3}) {[\begin{matrix} \cos β_{1} & 0 & \sin β_{1} \end{matrix}]}^{T} 。 \end{array}}$

(3)

对于斜入射SV波：

$\begin{array}{l} u_{SV 1}^{R} (t) = u_{SV 0} (t - Δ t_{4}) {[\begin{matrix} \cos α & 0 & - \sin α \end{matrix}]}^{T} ， \\ u_{SV 2}^{R} (t) = A_{3} u_{SV 0} (t - Δ t_{5}) {[\begin{matrix} - \cos α & 0 & - \sin α \end{matrix}]}^{T} ， \\ u_{P 2}^{R} (t) = A_{4} u_{SV 0} (t - Δ t_{6}) {[\begin{matrix} \sin β_{2} & 0 & - \cos β_{2} \end{matrix}]}^{T} 。 \end{array}}$

(4)

式中,β₁,β₂为斜入射P波和SV波在地表自由面的反射角,β₁=arcsin(c_ssinα/c_p),β₂=arcsin(c_psinα/c_s), $c P$ , $c SV$ 为P波和SV波波速； $A_{1}$ , $A_{2}$ 为反射P波和SV波幅值与入射P波幅值的比值, $A_{3}$ , $A_{4}$ 为反射SV波和P波幅值与入射SV波幅值的比值,其值可参考文献[16]取得； $Δ t$ 为入射波从零时刻波阵面传播到边界节点l所需的时间：

$\begin{array}{l} Δ t_{1} = | x_{l} \sin α + z_{l} \cos α | / c_{p} ， \\ Δ t_{2} = | x_{l} \sin α + [2 (H + L) - z_{l}] \cos α | / c_{p} ， \\ Δ t_{3} = | H + L - z_{l} | / (c_{s} \cos β_{1}) + \\ | x_{l} \sin α + (H + L) \cos α - (H + L - z_{l}) \sin α \tan β_{1} | / c_{p} ， \\ Δ t_{4} = | x_{l} \sin α + z_{l} \cos α | / c_{s} ， \\ Δ t_{5} = | x_{l} \sin α + [2 (H + L) - z_{l}] \cos α | / c_{s} ， \\ Δ t_{6} = | H + L - z_{l} | / (c_{p} \cos β_{2}) + \\ | x_{l} \sin α + (H + L) \cos α - (H + L - z_{l}) \sin α \tan β_{2} | / c_{s} 。 \end{array}}$

(5)

根据人工边界处内行场位移时程 $u_{l i}^{R} (t)$ ,可通过对时间求导或差分获得相应的速度场 ${\dot{u}}_{l i}^{R} (t)$ 。根据广义胡克定律,可求得人工边界处的应力场 $σ_{l i}^{R} (t)$ 。将上述求得的结果带入下式,由此实现深埋地下洞室近场地震动斜入射的等效荷载输入：

$f_{l i} = (K_{l i} u_{l i}^{R} + C_{l i} {\dot{u}}_{l i}^{R} + σ_{l i}^{R}) A_{l}$ 。

(6)

式中 $f_{l i}$ 为斜入射地震动作用下人工边界处节点l在i方向的等效节点力； $A_{l}$ 为人工边界面上节点 $l$ 的控制面积； $K_{l i}$ , $C_{l i}$ 为三维黏弹性人工边界处节点 $l (x_{l}, y_{l}, z_{l})$ 在方向i上的弹簧和阻尼元件参数,可参考文献[15]取值。

2. 地震荷载下层间接触面的劣化效应

地震荷载对岩体结构面抗剪强度的劣化作用是一个复杂的动态变化过程。本文采用震动劣化系数 $D (t)$ 来定量描述层间接触面在地震过程中的劣化衰减程度,则任意t时刻下层间接触面的抗剪强度可表示为

$τ_{p} (t) = τ_{p0} D (t) = σ (t) \tan φ_{0} D (t) + c_{0} D (t)$ 。

(7)

式中 $τ_{p0}$ , $τ_{p} (t)$ 为初始时刻和t时刻接触面峰值抗剪强度； $σ (t)$ 为t时刻接触面的法向应力； $φ_{0}$ , $c_{0}$ 为接触面初始内摩擦角和初始黏聚力。

总结前人研究成果^[12-14],地震荷载既是循环荷载也是动荷载,因此深埋互层岩体层间接触面在地震作用下的力学特性同时受到循环荷载的循环剪切作用和动荷载的剪切变形速率这两个因素影响^[14]。基于此,本文分别采用震动磨损影响系数 $η (t)$ 和相对速度影响系数 $γ (t)$ 来定量描述二者对接触面的影响。前者主要通过循环剪切次数和循环剪切幅值造成起伏不平的接触面磨损和钝化,降低其抗剪强度。后者主要是通过互层岩体岩层间在地震作用下的相对运动而产生的相对速度导致接触面的动摩擦系数降低。由于这两个因素对接触面强度劣化的作用机理是不同的,可假定在整个地震历时过程中,二者对接触面抗剪强度的劣化作用是相互独立的,则震动劣化系数 $D (t)$ 可表示为

$D (t) = η (t) γ (t)$ 。

(8)

2.1 震动磨损影响系数数学表达式

在地震过程中,层间接触面的循环剪切作用会造成界面的疲劳磨损和钝化,并随着地震持时不断累积,导致接触面抗剪强度参数逐渐降低。假定在地震过程中由循环剪切作用引起的接触面劣化程度按负指数关系衰减演化^{[14, 17]},则接触面由于地震动荷载的循环剪切作用造成的震动磨损影响系数 $η (t)$ 可表示为

$η (t) = R_{0} + (1 - R_{0}) e^{- λ (t / t_{s})}$ ,

(9)

式中,t为当前地震荷载作用时刻, $t_{s}$ 为地震总持时, $λ$ 为待定参数, $R_{0}$ 为震后接触面弱化程度收敛值,均由试验拟合。

2.2 相对速度影响系数数学表达式

针对相对剪切运动速度对层间接触面抗剪强度劣化的影响,李海波等^[18]的研究表明,岩石节理的峰值剪切强度与相对剪切速度关系呈现明显的负指数衰减形式：

$γ (t) = m {(| Δ {\dot{δ}}_{t} | + n)}^{- p}$ ,

(10)

式中, $Δ {\dot{δ}}_{t}$ 为t时刻下接触点对间的相对剪切速度,m,n,p均为待定参数,由试验拟合。

2.3 震动劣化系数表达式

将式（9）,（10）代入式（8）,得到t时刻下层间接触面震动劣化系数 $D (t)$ ：

$D (t) = [R_{0} + (1 - R_{0}) e^{- λ (t / t_{s})}] [m {(| Δ {\dot{δ}}_{t} | + n)}^{- p}]$ ,

(11)

式中,待定参数 $R_{0}$ , $λ$ ,m,n,p的选取均应根据实际工程的岩石材料剪切试验确定,本文从工程偏安全的角度出发,参考文献[14,18]的试验成果,选取 $R_{0}$ =0.75, $λ$ =5,m=0.883,n=0.02,p=0.032。

式（11）中接触点对间的相对剪切速度可由显式中心差分法求解。在地震动力响应计算的每一时步内,由此可显式求解接触面震动劣化系数 $D (t)$ ,实时刷新接触面的抗剪强度参数,从而动态获得整个地震历时过程中层间接触面的劣化衰减过程,并运用到下节接触面动力算法中。

3. 互层岩体层间接触系统动力响应分析模型

互层岩体层间地震动响应下的受力机理是一个复杂的动接触问题。地震荷载作用下,层间接触面可能发生脱开、滑移,从而对洞室结构造成严重的损伤破坏。针对互层岩体层间存在的复杂动接触行为,建立一种考虑接触面震动劣化效应和多种接触状态的层间接触系统动力响应分析模型。

对岩层间的接触系统采用有限元离散后,可得到地震作用下包含动接触力的节点运动微分方程：

$M \ddot{u} + C \dot{u} + K u = F + R$ 。

(12)

式中, $M$ , $C$ , $K$ 为层间接触节点系统的质量、阻尼和刚度矩阵； $\ddot{u}$ , $\dot{u}$ , $u$ 为接触节点的加速度、速度和位移向量； $F$ 为自重、地震等外荷载向量； $R$ 为动接触力向量, $R = N + T$ , $N$ , $T$ 为 $R$ 的法向和切向动接触力分量。

对任意一接触节点对i和 $i^{'}$ ,如图2所示,根据t时刻接触节点的位移 $u^{t}$ 、速度 ${\dot{u}}^{t}$ 和动接触力 $R^{t}$ ,可采用中心差分法求解 $t + Δ t$ 时刻的位移 $u^{t + Δ t}$ 和速度 ${\dot{u}}^{t + Δ t}$ ：

$u^{t + Δ t} = {\bar{u}}^{t + Δ t} + Δ t^{2} M^{- 1} R^{t} / 2$ ,

(13)

$\begin{matrix} {\bar{u}}^{t + Δ t} = Δ t^{2} M^{- 1} F^{t} / 2 + (I - Δ t^{2} M^{- 1} K / 2) u^{t} + \\ (Δ t I - Δ t^{2} M^{- 1} C / 2) {\dot{u}}^{t} \end{matrix}$ ,

(14)

${\dot{u}}^{t + Δ t} = 2 (u^{t + Δ t} - u^{t}) / Δ t - {\dot{u}}^{t}$ 。

(15)

图 2 接触模型和点对上动接触力示意图

Figure 2. Contact model and dynamic contact force on node pairs

下载: 全尺寸图片幻灯片

在上述求解迭代过程中,只有动接触力 $R^{t}$ 是未知量,需根据t— $t + Δ t$ 时刻的接触状态基于相应的接触条件求解。假定 $t + Δ t$ 时刻接触节点对处于黏结接触状态,则可引入接触边界约束条件,包括位移接触条件和接触力边界条件：

$\begin{array}{l} n_{i}^{T} (u_{i'}^{t + Δ t} - u_{i}^{t + Δ t}) = 0 ， \\ t_{i}^{T} (u_{i'}^{t + Δ t} - u_{i}^{t + Δ t}) = t_{i}^{T} (u_{i'}^{t} - u_{i}^{t}) ， \\ N_{i}^{t} = - N_{i'}^{t}, T_{i}^{t} = - T_{i'}^{t} ， \end{array}}$

(16)

式中, $n_{i}$ 为接触面在接触节点处的单位法向矢量,由节点 $i'$ 指向 $i$ , $t_{i}$ 为单位切向矢量。

将式（13）,（14）代入式（16）可得

$\begin{array}{l} N_{i}^{t} = \frac{2 M_{i} M_{i'}}{(M_{i} + M_{i'}) Δ t^{2}} Δ_{1 i} n_{i} ， \\ T_{i}^{t} = \frac{2 M_{i} M_{i'}}{(M_{i} + M_{i'}) Δ t^{2}} Δ_{2 i} t_{i} ， \\ Δ_{1 i} = n_{i}^{T} ({\bar{u}}_{i'}^{t + Δ t} - {\bar{u}}_{i}^{t + Δ t}) ， \\ Δ_{2 i} = t_{i}^{T} [({\bar{u}}_{i'}^{t + Δ t} - {\bar{u}}_{i}^{t + Δ t}) - (u_{i'}^{t} - u_{i}^{t})] 。 \end{array}}$

(17)

式中 $N_{i}^{t}$ , $T_{i}^{t}$ 为动接触力 $R_{i}^{t}$ 的法向和切向分量； $M_{i}$ , $M_{i'}$ 为节点 $i$ 和 $i^{'}$ 的集中质量； $Δ_{1 i}$ , $Δ_{2 i}$ 为接触点对的法向和切向接触间隙。

式（17）中动接触力 $N_{i}^{t}$ , $T_{i}^{t}$ 的计算公式是在接触节点对处于黏结接触的假定下求得。实际上,在强震过程中,当动接触力超过接触面黏结强度或抗剪强度后,层间界面还存在滑移或脱开等多种接触状态。因此每一时步计算完毕后,需根据接触面极限承载强度,对接触状态重新判别并修正动接触力。接触面法向极限抗拉力和切向剪切力由下式计算：

$\begin{array}{l} P_{N} = c D (t) A_{i} ， \\ P_{T} = μ_{s} ‖ N ‖ D (t) + c D (t) A_{i} 。 \end{array}}$

(18)

式中 $A_{i}$ 为接触节点i的控制面积； $c$ , $μ_{s}$ 为接触面的黏聚力和静摩擦系数； $D (t)$ 为本文第二节所求得的t时刻接触面震动劣化系数；若接触点对处于黏结接触状态,则需考虑界面黏聚力即 $c > 0$ 。若接触点对在地震加载过程中发生相对滑动,界面进入滑动接触状态,则不再考虑黏聚力即 $c = 0$ 。

（1）如果 $Δ_{1 i} < 0$ 且 $‖ N_{i}^{t} ‖ > P_{N}$ ,接触节点对处于分离状态,表明接触面受拉且发生拉裂破坏,则动接触力按照下式修正：

$N_{i}^{t} = 0, T_{i}^{t} = 0$ 。

(19)

（2）如果 $Δ_{1 i} < 0$ 且 $‖ N_{i}^{t} ‖ \leq P_{N}$ ,或者 $Δ_{1 i} > 0$ 且 $‖ N_{i}^{t} ‖ \leq P_{T}$ ,表明接触节点对处于黏结接触状态,此时接触面的损伤特性主要体现在层面的震动劣化上,无需对动接触力 $N_{i}^{t}$ 和 $T_{i}^{t}$ 进行修正。

（3）如果 $Δ_{1 i} > 0$ 且 $‖ N_{i}^{t} ‖ > P_{T}$ ,接触节点对处于压剪状态,表明接触面进入滑动接触且发生剪切滑移破坏。此时,须按照下式修正切向动接触力：

$T_{i}^{t} = μ_{d} ‖ N_{i}^{t} ‖ D (t) T_{i}^{t} / ‖ T_{i}^{t} ‖$ ,

(20)

式中, $μ_{d}$ 为接触面动摩擦系数。由式（17）～（20）可求解t时刻动接触力 $R^{t}$ ,将其回代入式（13）～（15）即可得到接触系统在下一时步的运动状态。

4. 工程实例

4.1 工程概况与计算模型

阿扎德帕坦水电站位于巴基斯坦巴控克什米尔地区的Jhelum河上,为该河梯级水电开发中的一级。地下厂房位于坝下游左岸山体内。工程区地质形态特征主要受地层岩性控制,岩性为60%的砂岩岩层与40%的非砂岩类岩层呈互层分布,地层主要为单斜构造,多呈陡倾、垂直展布,倾向NE,倾角∠65°～75°,层间剪切、挤压破碎带较为常见,对地下洞室高边墙的围岩稳定性影响较大。

阿扎德帕坦水电站工程区为地震强烈活动地区,主要受印度板块向欧亚板块持续俯冲运动影响,该区域表现为强烈的现今地震活动性。根据中国地震局地质研究所研究成果,工程区50 a超越概率10%（DBE）的峰值加速度为0.32g,对应的地震基本烈度为Ⅷ度。该区地震地质背景总体较为复杂。地下厂房采用单洞室方案,埋深一般超过300 m,主厂房长×宽×高为175.6 m×24.9 m×60.15 m,具有大跨度和高边墙的特点。主厂房内部衬砌厚度设置为顶拱1.25 m,边墙1.0 m。主厂房下部还包括机墩等大体积混凝土结构。

选取如图3所示含泥、砂岩互层的地下厂房区域建立数值模型。有限元模型采用8节点6面体单元离散,模型中包括砂岩层、非砂岩层、衬砌和厂房内部混凝土结构。计算模型共剖分了个138972实体单元和147975个节点,其中混凝土衬砌单元18194个。三维计算坐标x轴正向为N148ºE,与厂房纵轴线方向垂直,指向顺水流方向为正；y轴向为N58ºE,与厂房纵轴线方向重合,从^#1机组指向^#4机组为正；z轴与大地坐标重合,竖直向上为正。

图 3 三维有限元模型

Figure 3. 3D finite element model for underground powerhouse

下载: 全尺寸图片幻灯片

三维初始地应力场根据设计院提供的地应力资料反演得到,侧压力系数取k_x=1.4,k_y=1.2,k_z=1.15。计算采用的材料物理力学参数见表1。在进行动力分析前,首先对地下洞室进行静力开挖与支护模拟,将喷锚支护开挖完毕后得到的围岩应力场作为动力计算前的初始应力场。

表 1 模型材料力学参数

Table 1. Mechanical parameters of model materials

材料	变形模量/GPa	泊松比	黏聚力/MPa	内摩擦角/(°)	抗拉强度/MPa
砂岩	8.0	0.240	1.10	42.0	1.95
非砂岩	2.0	0.300	0.35	26.0	0.25
衬砌	28.0	0.167	1.80	46.0	1.27
接触面	—	—	1.00	30.0	—

下载: 导出CSV

| 显示表格

4.2 计算条件

计算程序采用课题组自主研发的大型地下洞室群地震灾变模拟系统^[19],并将本文的深埋互层岩体地下洞室动力响应分析方法嵌入其中。地震加载前,首先在模型砂岩层单元与非砂岩层单元相交界面处设置接触面,通过增加界面两侧的共用接触节点对,以完成软岩与硬岩单元的节点分离。围岩采用基于M-C准则的动力弹塑性损伤本构模型^[19],混凝土衬砌采用塑性损伤本构模型模拟^[20],标量损伤系数d为

$\begin{array}{l} d = 1 - (1 - d_{c}) (1 - s d_{t}) ， \\ d_{t} = 1 - \frac{(1 - A_{t}) ε_{t}^{f}}{\bar{ε}} - \frac{A_{t}}{\exp [B_{t} (\bar{ε} - ε_{t}^{f})]} ， \\ d_{c} = 1 - \frac{(1 - A_{c}) ε_{c}^{f}}{\bar{ε}} - \frac{A_{c}}{\exp [B_{c} (\bar{ε} - ε_{c}^{f})]} 。 \end{array}}$

(21)

式中 $d_{t}$ , $d_{c}$ 为单元受拉和受压损伤系数；s为单元由受拉转变为受压状态的刚度恢复系数, $0 \leq$ $s \leq 1$ ； $A_{t}$ , $B_{t}$ 和 $A_{c}$ , $B_{c}$ 分别为单元相应的拉、剪损伤常数,可根据试验结果拟合得到； $ε_{t}^{f}$ , $ε_{c}^{f}$ 为拉、剪损伤阈值； $\bar{ε} = \sqrt{{〈 ε_{1} 〉}^{2} + {〈 ε_{2} 〉}^{2} + {〈 ε_{3} 〉}^{2}}$ , $ε_{1}$ , $ε_{2}$ , $ε_{3}$ 为单元的3个主应变, $〈 ε_{i} 〉$ 的函数式定义为 $〈 ε_{i} 〉 = (ε_{i} + | ε_{i} |) / 2 ， i = 1, 2, 3$ 。在实际计算中,应保证拉伸准则具有更高的优先级,若混凝土单元受力状态满足最大拉应变准则,则产生拉损,且不需再判断该单元是否满足Mohr-Coulomb屈服准则；否则的话,采用Mohr-Coulomb屈服准则判别其压应变状态。

采用1940 El-Centro波作为动力计算输入地震波,并根据阿扎德帕坦水电站工程区抗震设防烈度,将峰值加速度调整为0.32g,截取其中地震动较为强烈的前20 s时段作为入射波,经滤波、基线校正后得到如图4所示的加速度时程曲线。鉴于SV波（沿x方向振动）对地下洞室的影响要大于SH波（沿y方向振动）,动力计算同时考虑P波和SV波对洞室的作用。其中,SV波采用如图4所示的入射波,P波加速度取为SV波的2/3。

图 4 入射波加速度时程曲线

Figure 4. Time-history curves of input wave acceleration

下载: 全尺寸图片幻灯片

为研究地下厂房结构不同部位的地震动响应时程特性,选取Y=0软岩层中间断面为监测断面,取该断面主厂房衬砌结构上4个监测点,用以监测衬砌结构的位移和应力等响应特性,另取监测点分别位于软岩与硬岩层接触面两侧围岩内,用以监测层间相对运动特征,监测点布置如图5所示。抗震计算分3种工况：①地震动竖直入射,不考虑层面劣化和动接触；②地震动斜入射,不考虑层面劣化和动接触；③地震动斜入射,且考虑层面劣化和动接触。地震动斜入射的入射方向矢量为 $(1 / 2, 0, \sqrt{3} / 2)$ ,入射基准点为模型坐标原点。

图 5 监测点布置

Figure 5. Layout of monitoring points

下载: 全尺寸图片幻灯片

5. 地震响应时程分析

5.1 互层岩体层间相对运动及层面地震劣化分析

在强震加载过程中,软、硬岩层间循环往复相互作用,近场斜入射地震动加剧了层间复杂的动接触行为,产生相对位移,进而导致接触面产生震动劣化效应。3种工况下围岩层间相对位移时程曲线如图6所示。工况①下,围岩层间相对位移在0线上下波动变化,层间最大相对位移为1.28 cm,震后基本为0。工况②下,受斜入射地震动输入非一致性影响,围岩层间相对位移相比工况①波动剧烈,波动范围主要为-1.20～1.50 cm,最大相对位移为2.15 cm,但震后逐渐减小为0。工况③下,考虑动接触作用后,在2～7 s时间段内,层面震动劣化效应加剧了层间相对运动,且层间发生了明显的错动,最大相对位移为-6.01 cm；后期该相对位移值基本在-5.0～-6.0 cm上下波动,震后为-5.93 cm,表明考虑层面劣化和动接触后,强震作用下软、硬岩层震动不同步,层间发生了明显的剪切滑移破坏,震后接触面处于滑动状态。

图 6 层间接触面相对位移时程曲线

Figure 6. Time-history curves of relative displacement of interface

下载: 全尺寸图片幻灯片

工况③下主厂房边墙和顶拱部位软岩与硬岩层面震动劣化系数时程曲线如图7所示。由图可知,边墙与顶拱的软、硬岩层间接触面震动劣化系数 $D (t)$ 总体上随地震持时呈现明显的负指数衰减规律。在0～7 s, $D (t)$ 迅速衰减,震动劣化效应明显；7 s之后, $D (t)$ 不再继续衰减,只随接触面相对剪切运动速度的影响而在收敛值上下波动。主厂房高边墙相比顶拱部位的层面震动劣化效应明显,尤其在入射波加速度峰值过后时间段内,边墙处 $D (t)$ 受相对速度的影响波动幅度相比顶拱较大,前者主要在0.71～0.78上下波动,后者主要在0.75附近小幅度波动。震动劣化系数 $D (t)$ 的时程波动曲线直观地反映了整个强震历时过程中软、硬岩层间接触面的动态劣化过程。

图 7 震动劣化系数时程曲线

Figure 7. Time-history curves of vibration deterioration coefficient

下载: 全尺寸图片幻灯片

5.2 地下厂房衬砌结构位移时程分析

鉴于地震过程中地下洞室高边墙和顶拱相比底部破坏严重^{[17, 21]},故本节主要分析3种工况下边墙（B点）和顶拱（A点）位移时程波动规律,如图8所示。本文采用3个监测点（A,B,C）与其洞室底部测点（D）的峰值位移差来表征相对位移,如图9所示。由图8,9可以看出3种工况下：①地下厂房衬砌结构各监测点合位移时程曲线的波形和波动规律较为相似,均出现多个明显的波峰,不同部位几乎同时出现峰值和谷值,表明衬砌结构不同部位处于同步震动状态；②在输入地震动较为强烈的前5 s内,各监测点位移时程曲线大幅波动,且洞室边墙位移明显大于顶拱和底部,易引起洞室结构相对应区域的非一致性变形,这也与文献[21]的研究结果相符。

图 8 监测点合位移时程曲线

Figure 8. Time-history curves of displacement of monitoring points

下载: 全尺寸图片幻灯片

图 9 监测点峰值位移差

Figure 9. Peak values of relative displacements

下载: 全尺寸图片幻灯片

工况①下,主厂房衬砌结构各监测点最大位移为6.83 cm,左、右边墙和顶拱与底部的峰值位移差分别为3.7,3.5,2.1 cm。工况②下,考虑地震动的三维斜入射特性后,其对洞室结构位移时程响应的影响相比工况①主要表现在波动幅值上,衬砌各监测点最大位移为8.14 cm,左、右边墙和顶拱的峰值位移差分别为4.3,3.9,2.3 cm,表明洞室衬砌结构的位移响应受地震动入射角影响较大,且洞室结构产生了较大的非一致性变形。这主要是因为斜入射考虑了地震动入射在模型人工边界上的方向性和非一致性,同时考虑了地表反射波对深埋洞室工程区波动场的影响。斜入射地震动会在地表自由面发生复杂的波形转换现象,致使地震动在洞室周围发生衍射效应,故模型人工边界处地震波场是由不同入射波和反射波叠加而成,使得人工边界上各节点具有不同的振动方向和振动波形,产生放大效应。这是地震动竖直入射时所没有的波场特性。

工况③下,考虑层面劣化和动接触后,各监测点最大位移为10.06 cm,左、右边墙和顶拱的峰值位移差分别为5.9,4.8,3.0 cm。表明考虑软、硬岩层间接触面震动劣化条件下,界面抗剪强度参数在强震历时过程中逐步劣化,无法提供足够的抗滑力,且斜入射地震动加剧了震动劣化效应,层间的剪切滑动导致了厂房边墙产生了较大的相对变形,易引起洞室围岩的垮塌和大变形。

5.3 地下厂房衬砌结构应力时程分析

由于拉裂破坏是地震作用下地下厂房衬砌结构的主要破坏模式^[19],因此本文主要分析3种工况下厂房衬砌边墙（B点）和顶拱（A点）最大主应力时程波动规律,如图10所示。3种工况下厂房顶拱和边墙最大主应力时程曲线变化规律基本一致,均表现为在0～6 s,最大主应力波动剧烈,增幅明显；6 s之后随着输入地震加速度的衰减,波动趋于平缓。这种现象与输入地震波的波动规律较为相似,表明厂房衬砌结构的最大主应力波动规律主要受到输入地震波的影响。

图 10 监测点最大主应力时程曲线

Figure 10. Time-history curves of maximum principal stress of monitoring points

下载: 全尺寸图片幻灯片

3种工况下衬砌顶拱的最大主应力分别为1.18,1.30,1.49 MPa,边墙的最大主应力分别为1.31,1.53,1.87 MPa。工况②和工况③下边墙和顶拱的最大拉应力明显超过了混凝土的抗拉强度,表明厂房边墙和顶拱发生拉裂损伤破坏。当考虑地震动斜入射时,边墙和顶拱最大拉应力幅值相比竖直入射工况增加16.8%和10.2%。当考虑层面劣化和动接触后,工况③下边墙和顶拱最大拉应力幅值相比工况②增加22.2%和14.6%。表明受地表自由面反射波场叠加效应的影响,斜入射时厂房衬砌结构应力水平的地震响应远大于竖直入射,互层岩体层面震动劣化和层间相对滑移对厂房衬砌结构受力有重要影响,加剧了顶拱和高边墙拉应力的波动,使其发生严重的疲劳破坏。

5.4 地下厂房衬砌结构震损分析

采用衬砌结构损伤系数可以直观地对厂房震后损伤破坏程度进行量化,选取软岩穿过部位的机组段进行震损分析,图11为不同工况下软岩穿过部位机组段衬砌结构损伤系数分布。工况①下,衬砌结构震损区域主要分布在顶拱、高边墙、发电机层楼板和水轮机层楼板等结构,大部分区域损伤系数量值较小,在0.3左右。考虑斜入射后,损伤区范围沿着洞室纵轴向和竖直向进一步延伸,主厂房顶拱损伤区与上、下游高边墙损伤区有连通的趋势。考虑层面劣化和动接触后,地下洞室衬砌结构损伤程度达到最大,主厂房结构大部分区域被损伤区覆盖,局部区域如顶拱、高边墙中心部位和风罩损伤严重,其最大损伤系数达到0.9。

图 11 不同工况下衬砌结构损伤系数分布

Figure 11. Distribution of damage coefficient of underground structures under different cases

下载: 全尺寸图片幻灯片

3种工况下地下厂房衬砌结构损伤类型分布如图12所示,横向上看,厂房衬砌上部结构相比下部结构更易发生损伤破坏,其中拉裂损伤区主要分布在厂房顶拱顶部、上部高边墙和楼板层,剪切损伤区多分布在厂房顶拱与整个边墙,拉裂损伤区的损伤系数比剪切损伤区的更大。主要由于上部结构厚度小,只受围岩单向约束,具有较大的临空面,斜入射地震波在临空面产生衍射和放大效应,导致这些部位损伤破坏更为严重。纵轴向上看,拉裂损伤主要分布在互层岩体软岩层穿过的中心部位,剪切损伤主要分布在层间错动面上,并呈现从错动面向两侧扩展分布的趋势,形成一个层间剪切、挤压破碎带。

图 12 不同工况下衬砌结构损伤类型分布

Figure 12. Distribution of damage types of underground structures under different cases

下载: 全尺寸图片幻灯片

根据数值模拟的结果,可以归纳出地震作用下互层岩体地下厂房2种主要的破坏模式,由结构变形和层间错动引起的结构破坏。如图13（a）所示,地下厂房主要可以分为由顶拱、上部边墙和楼板层组成的上部结构和由机墩大体积混凝土组成的下部结构。斜入射地震动在传播过程中会对结构产生放大的水平力,引起上、下部结构的相对变形,如图13（b）,震后厂房结构x向位移由顶拱至机墩底部逐渐减小,上部结构变形大于下部,上、下部结构水平横向相对位移在3.1～4.2 cm,洞室经历了较大的结构变形,最终导致顶拱和上部边墙损伤破坏严重。

图 13 地下厂房结构破坏模式与数值结果对比

Figure 13. Comparison between failure modes and numerical results

下载: 全尺寸图片幻灯片

层间错动下厂房结构综合响应示意图如图13（c）所示。地震作用下,互层岩体软岩穿过部位衬砌变形与层间错动面两侧衬砌变形差异显著（图13（d））,易导致厂房衬砌结构在该处发生剪切破坏。

6. 结论

（1）基于深埋近场地震动斜入射传播特性和考虑层面震动劣化和黏结滑移特性的动接触力分析方法,建立一种深埋互层岩体地下洞室地震动力响应分析模型,并将该模型应用于巴基斯坦阿扎德帕坦水电站地下厂房地震动响应分析中。

（2）地震作用下厂房衬砌结构不同部位监测点的位移时程曲线和应力时程曲线较为相似。衬砌结构的位移和应力响应受地震动入射角影响较大,主要体现在波动幅值上,斜入射地震动输入的非一致性加剧了厂房的地震响应。考虑层面震动劣化和动接触作用后,层面附近衬砌位移和应力进一步加大,厂房上部边墙和顶拱拉应力超过衬砌的抗拉强度。

（3）地震作用下软、硬岩层震动不同步,层间易发生相互错动,最大错动位移达5.93 cm,层间发生明显的剪切滑移破坏和震动劣化效应,厂房上部边墙相比顶拱劣化明显。层面震动劣化系数随地震持时呈负指数衰减,其时程变化曲线可直观反映层面的动态劣化过程。

（4）横向上看,厂房衬砌结构损伤区主要分布在上下游边墙、顶拱及其与边墙交界处和发电机层楼板,其中拉裂损伤主要分布在顶拱顶部和上部高边墙,该区域损伤系数较大,破坏严重；纵轴向上看,拉裂损伤主要分布在互层岩体软岩层穿过的中心部位,剪切损伤主要分布在层间错动面上,并呈现从层面向两侧扩展分布的趋势,形成一个层间剪切、挤压破碎带。归纳总结地震作用下互层岩体地下厂房2种主要的破坏模式,地震荷载对结构产生的放大水平力导致结构变形和层间错动导致接触面处衬砌结构发生剪切变形。

图 1 深埋地下洞室地震动斜入射示意图

Figure 1. 3D diagram of oblique incident earthquake for a deeply-buried underground cavern

下载: 全尺寸图片幻灯片

图 2 接触模型和点对上动接触力示意图

Figure 2. Contact model and dynamic contact force on node pairs

下载: 全尺寸图片幻灯片

图 3 三维有限元模型

Figure 3. 3D finite element model for underground powerhouse

下载: 全尺寸图片幻灯片

图 4 入射波加速度时程曲线

Figure 4. Time-history curves of input wave acceleration

下载: 全尺寸图片幻灯片

图 5 监测点布置

Figure 5. Layout of monitoring points

下载: 全尺寸图片幻灯片

图 6 层间接触面相对位移时程曲线

Figure 6. Time-history curves of relative displacement of interface

下载: 全尺寸图片幻灯片

图 7 震动劣化系数时程曲线

Figure 7. Time-history curves of vibration deterioration coefficient

下载: 全尺寸图片幻灯片

图 8 监测点合位移时程曲线

Figure 8. Time-history curves of displacement of monitoring points

下载: 全尺寸图片幻灯片

图 9 监测点峰值位移差

Figure 9. Peak values of relative displacements

下载: 全尺寸图片幻灯片

图 10 监测点最大主应力时程曲线

Figure 10. Time-history curves of maximum principal stress of monitoring points

下载: 全尺寸图片幻灯片

图 11 不同工况下衬砌结构损伤系数分布

Figure 11. Distribution of damage coefficient of underground structures under different cases

下载: 全尺寸图片幻灯片

图 12 不同工况下衬砌结构损伤类型分布

Figure 12. Distribution of damage types of underground structures under different cases

下载: 全尺寸图片幻灯片

图 13 地下厂房结构破坏模式与数值结果对比

Figure 13. Comparison between failure modes and numerical results

下载: 全尺寸图片幻灯片

表 1 模型材料力学参数

Table 1 Mechanical parameters of model materials

材料变形模量/GPa 泊松比黏聚力/MPa 内摩擦角/(°) 抗拉强度/MPa

砂岩 8.0 0.240 1.10 42.0 1.95
非砂岩 2.0 0.300 0.35 26.0 0.25
衬砌 28.0 0.167 1.80 46.0 1.27
接触面 — — 1.00 30.0 —

下载: 导出CSV

参考文献(21)

[1]	SHEN Y S, GAO B, YANG X M, et al. Seismic damage mechanism and dynamic deformation characteristic analysis of mountain tunnel after Wenchuan earthquake[J]. Engineering Geology, 2014, 180: 85-98. doi: 10.1016/j.enggeo.2014.07.017
[2]	李小军, 卢滔. 水电站地下厂房洞室群地震反应显式有限元分析[J]. 水力发电学报, 2009, 28(5): 41-46. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB200905008.htm LI Xiao-jun, LU Tao. Explicit finite element analysis of earthquake response for underground caverns of hydropower stations[J]. Journal of Hydroelectric Engineering, 2009, 28(5): 41-46. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB200905008.htm
[3]	杜修力, 陈维, 李亮, 等. 斜入射条件下地下结构时域地震反应分析初探[J]. 震灾防御技术, 2007, 2(3): 290-296. doi: 10.3969/j.issn.1673-5722.2007.03.009 DU Xiu-li, CHEN Wei, LI Liang, et al. Preliminary study of time-domain seismic response for underground structures to obliquely incident seismic waves[J]. Technology for Earthquake Disaster Prevention, 2007, 2(3): 290-296. (in Chinese) doi: 10.3969/j.issn.1673-5722.2007.03.009
[4]	尤红兵, 赵凤新, 荣棉水. 地震波斜入射时水平层状场地的非线性地震反应[J]. 岩土工程学报, 2009, 31(2): 234-240. doi: 10.3321/j.issn:1000-4548.2009.02.014 YOU Hong-bing, ZHAO Feng-xin, RONG Mian-shui. Nonlinear seismic response of horizontal layered site due to inclined wave[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(2): 234-240. (in Chinese) doi: 10.3321/j.issn:1000-4548.2009.02.014
[5]	杜修力, 黄景琦, 赵密, 等. SV 波斜入射对岩体隧道洞身段地震响应影响研究[J]. 岩土工程学报, 2014, 36(8): 1400-1406. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201408005.htm DU Xiu-li, HUANG Jing-qi, ZHAO Mi, et al. Effect of oblique incidence of SV waves on seismic response of portal sections of rock tunnels[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(8): 1400-1406. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201408005.htm
[6]	赵宝友, 马震岳, 丁秀丽. 不同地震动输入方向下的大型地下岩体洞室群地震反应分析[J]. 岩石力学与工程学报, 2010, 29(增刊1): 3395-3402. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2010S1117.htm ZHAO Bao-you, MA Zhen-yue, DING Xiu-li. Seismic response of a large underground rock cavern groups considering different incident angles of earthquake waves[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(S1): 3395-3402. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2010S1117.htm
[7]	赵密, 欧阳文龙, 黄景琦, 等. P波作用下跨断层隧道轴线地震响应分析[J]. 岩土力学, 2019, 40(9): 3645-3655. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201909039.htm ZHAO Mi, OUYANG Wen-long, HUANG Jing-qi, et al. Analysis of axis dynamic response of rock tunnels through fault fracture zone under P waves of earthquake[J]. Rock and Soil Mechanics, 2019, 40(9): 3645-3655. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201909039.htm
[8]	黄景琦, 杜修力, 田志敏, 等. 斜入射SV波对地铁车站地震响应的影响[J]. 工程力学, 2014, 31(9): 81-88, 103. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201409013.htm HUANG Jing-qi, DU Xiu-li, TIAN Zhi-min, et al. Effect of the oblique incidence of seismic SV waves on the seismic response of subway station structure[J]. Engineering Mechanics, 2014, 31(9): 81-88, 103. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201409013.htm
[9]	张志国. 地下洞室群地震响应数值分析方法研究[D]. 武汉: 武汉大学, 2012: 76-87. ZHANG Zhi-guo. Study on Numerical Simulation Methods for Seismic Response of Underground Cavern Complexes[D]. Wuhan: Wuhan University, 2012: 76-87. (in Chinese)
[10]	KOURETZIS G P, BOUCKOVALAS G D, GANTES C J. 3-D shell analysis of cylindrical underground structures under seismic shear (S) wave action[J]. Soil Dynamics and Earthquake Engineering, 2006, 26(10): 909-921. doi: 10.1016/j.soildyn.2006.02.002
[11]	LIU J B, SHARAN S K. Analysis of dynamic contact of cracks in viscoelastic media[J]. Computer Methods in Applied Mechanics and Engineering, 1995, 121(1/2/3/4): 187-200.
[12]	LEE H S, PARK Y J, CHO T F, et al. Influence of asperity degradation on the mechanical behavior of rough rock joints under cyclic shear loading[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(7): 967-980.
[13]	刘博, 李海波, 朱小明. 循环剪切荷载作用下岩石节理强度劣化规律试验模拟研究[J]. 岩石力学与工程学报, 2011, 30(10): 2033-2039. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201110011.htm LIU Bo, LI Hai-bo, ZHU Xiao-ming. Experiment simulation study of strength degradation of rock joints under cyclic shear loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 2033-2039. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201110011.htm
[14]	倪卫达, 唐辉明, 刘晓, 等. 考虑结构面震动劣化的岩质边坡动力稳定分析[J]. 岩石力学与工程学报, 2013, 32(3): 492-500. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201303008.htm NI Wei-da, TANG Hui-ming, LIU Xiao, et al. Dynamic stability analysis of rock slope considering vibration deterioration of structural planes under seismic loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(3): 492-500. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201303008.htm
[15]	刘晶波, 王振宇, 杜修力, 等. 波动问题中的三维时域粘弹性人工边界[J]. 工程力学, 2005, 22(6): 46-51. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX200506007.htm LIU Jing-bo, WANG Zhen-yu, DU Xiu-li, et al. Three- dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics, 2005, 22(6): 46-51. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX200506007.htm
[16]	杜修力. 工程波动理论与方法[M]. 北京: 科学出版社, 2009: 15-21. DU Xiu-li. Theories and Methods of Wave Motion for Engineering[M]. Beijing: Science Press, 2009: 15-21. (in Chinese)
[17]	ZHOU Hao, XIAO Ming, YANG Yang, et al. Seismic response analysis method for lining structure in underground cavern of hydropower station[J]. KSCE Journal of Civil Engineering, 2019, 23(3): 1236-1247.
[18]	李海波, 冯海鹏, 刘博. 不同剪切速率下岩石节理的强度特性研究[J]. 岩石力学与工程学报, 2006, 25(12): 2435-2440. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200612011.htm LI Hai-bo, FENG Hai-peng, LIU Bo. Study on strength behaviors of rock joints under different shearing deformation velocities[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(12): 2435-2440. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200612011.htm
[19]	张志国, 肖明, 陈俊涛. 大型地下洞室地震灾变过程三维动力有限元模拟[J]. 岩石力学与工程学报, 2011, 30(3): 509-523. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201103011.htm ZHANG Zhi-guo, XIAO Ming, CHEN Jun-tao. Simulation of earthquake disaster process of large-scale underground caverns using three-dimensional dynamic finite element method[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(3): 509-523. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201103011.htm
[20]	MAZARS J. A description of micro- and macro-scale damage of concrete structures[J]. Engineering Fracture Mechanics, 1986, 25(5/6): 729-737.
[21]	隋斌, 朱维申, 李晓静. 地震荷载作用下大型地下洞室群的动态响应模拟[J]. 岩土工程学报, 2008, 30(12): 1877-1882. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200812024.htm SUI Bin, ZHU Wei-shen, LI Xiao-jing. Simulation on dynamic response of large underground opening complex under seismic loads[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(12): 1877-1882. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200812024.htm