基于耐震时程分析法的高效地下结构地震易损性分析

甄立斌; 史跃波; 钟紫蓝; 杜修力; 罗文林

doi:10.11779/CJGE20220188

基于耐震时程分析法的高效地下结构地震易损性分析 English Version

1.
北京市勘察设计研究院有限公司, 北京 100038
2.
北京工业大学城市与工程安全减灾教育部重点实验室, 北京 100124

基金项目:

国家自然科学基金面上项目 51978020

国家自然科学基金联合基金项目 U1839201

详细信息

作者简介:
甄立斌(1994—)，男，硕士，工程师，主要从事地下结构地震工程方面的研究工作。E-mail: bjutzlb666@163.com

通讯作者:
钟紫蓝, E-mail: zilanzhong@bjut.edu.cn

中图分类号: TU311
计量
- 文章访问数: 315
- HTML全文浏览量: 77
- PDF下载量: 73
出版历程
- 收稿日期: 2022-02-21
- 网络出版日期: 2023-04-16
- 刊出日期: 2023-03-31

Efficient seismic fragility of underground structures using endurance time analysis method

1.
BGI Engineering Consultants Ltd., Beijing 100038, China
2.
Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing 100124, China

摘要

摘要: 为研究基于耐震时程分析法对地铁车站进行高效地震易损性分析的可行性，选取Ⅲ类工程场地中典型两层三跨地铁车站作为研究对象进行地震易损性分析。为反映所构造的耐震时程曲线具有所选天然地震记录的离散特性，借鉴IDA分析中分位数回归统计方法，分别以所选实际天然地震记录反应谱的16%，50%，84%分位值作为目标反应谱，各自构造3条人工耐震时程曲线作为有限元模型输入，并基于耐震时程分析结果与IDA结果绘制各级性能水准下的地铁车站地震易损性曲线。结果表明：在Ⅲ类工程场地中两种方法计算得到的地下结构各级性能水准地震易损性曲线吻合良好。因此，考虑输入耐震时程曲线离散特性的耐震时程分析法能够显著降低地下结构进行地震易损性分析过程中的非线性动力时程计算次数，并保证较高的计算精度；同时也验证了基于地震记录反应谱的分位值构造一系列具有离散特性耐震时程曲线进行易损性分析的思路是可行的。
- 耐震时程分析法 /
- 地铁车站 /
- 地震易损性 /
- 离散性 /
- 反应谱分位值
Abstract: The feasibility of efficient seismic fragility analysis of underground structures using the endurance time analysis (ETA) method is investigated. A typical two-story three-span subway station embedded in the site class Ⅲ is used as the prototype. To effectively consider the inherent dispersion of as-recorded earthquake ground motions during the simulation of endurance time acceleration functions, the 16%, 50% and 84% quantiles of response spectra of the as-recorded ground motions at the engineering bedrock are used as the target response spectra to generate three endurance time acceleration functions using the quantile regression method in the incremental dynamic analysis (IDA) for reference as the input motions for the nonlinear soil-structure interaction system, respectively. The seismic fragility evaluation of the underground structures is performed based on the ETA results and the IDA results. It can be seen from the numerical results that the seismic fragility curves obtained by the ETA are in good agreement with those obtained by the IDA. Therefore, the proposed ETA with full consideration of the dispersion of ground motions can significantly reduce the number of nonlinear dynamic analyses, and ensure high accuracy of the computed results. Moreover, the proposed ETA using a series of discrete ETAFs by the quantiles of response spectra of the as-recorded ground motions is a feasible alternative for the seismic fragility analysis of underground structures.
- endurance time analysis method /
- subway station /
- seismic fragility /
- dispersion /
- quantile of response spectrum

HTML全文

0. 引言

在波浪荷载作用下，海床内孔隙水压力增大以及垂直有效应力降低，进而导致海床土体失稳甚至液化^[1]。在人类探索海洋资源的近几十年内，海床的液化对人类在海洋中的作业安全产生了很大的影响^[2]。因此，针对波浪荷载对海床动力响应的研究具有重要的理论意义和工程价值。

自20世纪中叶起，学者们开始分别基于刚性非耦合模型^[3]、Biot流固耦合模型^[4]和多孔介质弹塑性模型^[5]对波浪荷载引起的海床动力响应问题进行了理论研究。其中基于可考虑饱和介质压缩性的Biot三维准静态固结理论^[6]和动力固结理论应用^[7]最为广泛。Yamamoto等^[8]基于Biot准静态固结理论得到了线性波作用下无限厚度海床内孔隙水压力、有效应力和土骨架位移的闭合形式解，为此后的研究奠定了基础。Madsen^[9]在此基础上，引入了Verruijt^[10]提出的考虑海床土体饱和度的储存方程，得到了线性波作用下无限厚度准饱和海床动力响应的解析解，发现海床的动力响应受到海床土体饱和度的影响显著。基于Madsen的研究，Okusa^[11]发现海床土体的参数对无限厚度海床动力响应的影响可以忽略不计。Mei等^[12]采用边界层近似法，即通过将海床分为上下两部分，求解无限厚度海床动力响应的方法求解上半区域海床，并得到下半区域海床的简化解，最终得到波浪荷载作用下有限厚度多孔弹性海床的动力响应的近似解。Hsu等^[13]则得到了波浪荷载作用下有限厚度海床动力响应的解析解，发现当海床厚度趋于无穷大时，有限厚度海床的动力响应与无限厚度海床的解析解结果基本一致。

Zienkiewicz等^[14]基于Biot准静态固结理论和Biot动力固结理论，整理总结得到了适应于饱和土体的准静态理论、部分动态理论和完整动态理论。Ulker等^[15]通过对准静态解、部分动态解和完整动态解的对比分析研究，发现在一定的海床和波浪条件下三者得到的结果基本一致。考虑到准静态理论的推导简单，便于用于实际工程分析，于是学者们在准静态理论的基础上进行了一系列深入的研究工作。海洋中波浪特性往往相当复杂^[16]，而且波浪非线性特性对海床动力响应的研究尚不够深入。Madsen^[9]虽然得到了无限厚度海床动力响应的解答，但是该研究没有考虑海床底部下卧基岩的情况以及波浪非线性的影响。王忠涛等^[17]则通过数值模拟的方法给出了波浪非线性对海床动力响应的影响。郭秀军等^[18]研究基于有限差分法研究了海床的分层对动力响应的影响。Zhou等^[19]和Zhang等^[20]考虑了海床分层的影响得到了二阶Stokes波作用下多层海床的半解析解。周晓智等^[21]研究了驻波作用下海床应力路径特性。但是前述学者的研究均没有考虑波浪参数对海床动力响应的影响，且没有得到二阶Stokes波作用下有限厚度海床动力响应的准静态解析解。

因此本文基于准静态理论求解均质各向同性海床在二阶Stokes波作用下动力响应的解析解，通过与既有解进行对比验证，并结合数值计算分析了波浪和海床土体的参数对海床动力响应的影响规律。

1. 控制方程和边界条件

如图 1所示，本文考虑了下卧刚性不透水基岩的多孔弹性海床在二阶Stokes波作用下动力响应问题。其中二阶Stokes波的传播方向为x轴的正方向。与Li等^[1]的研究一致，取海床和海洋的接触面为笛卡尔坐标系的零点，z轴垂直向上时为正值，负值在海床和海洋的接触面以下。此外，图中H和L分别为二阶Stokes波的波高和波长；d为水深；h为海床的厚度。

图 1 波浪和多孔弹性海床之间的动力相互作用

Figure 1. Dynamic interaction between ocean waves and poroelastic seabed

下载: 全尺寸图片幻灯片

1.1 控制方程

基于Biot三维准静态固结理论^[6]与Verruijt储能方程^[10]，多孔弹性海床的控制方程可以表示为

${\sigma _{ij,j}} = 0 \text{，}$

(1)

$- {p_{,i}} = {{{\rho _{\text{f}}}g{{\dot w}_i}} \mathord{\left/ {\vphantom {{{\rho _{\text{f}}}g{{\dot w}_i}} K}} \right. } K} 。$

(2)

式中：下标i和j表示该场量沿x或z方向的分量；下标‘, i’和‘, j’表示该场量对x或z方向的一阶偏导；字母上的‘•’表示该场量对时间t的一阶导数； ${\sigma _{ij}}$ 为总应力张量；p为孔隙水压力； ${\rho _{\text{f}}}$ 为孔隙流体密度；g为重力加速度；K为渗透系数（水力传导系数）。

流体的质量守恒方程为

${\dot u_{i,i}} + {\dot w_{i,i}} = - n\beta \dot p = \partial {\rho _{\text{f}}}/\partial t 。$

(3)

式中：u为土骨架的位移；w为孔隙流体相对土骨架的位移；n为海床土体的孔隙率；β为孔隙流体压缩性系数。当海床土体中含有少量气体时，β可以通过下式计算：

$\beta = 1/{K_{\text{w}}} + (1 - {S_{\text{r}}})/{P_{{\text{w0}}}} 。$

(4)

式中： ${K_{\text{w}}}$ 为孔隙流体的实际体积弹性模量，一般取 ${K_{\text{w}}}$ =2×10⁹ Pa； ${S_{\text{r}}}$ 为海床土体的饱和度； ${P_{{\text{w0}}}}$ 为由水深d确定的绝对压力， ${P_{{\text{w0}}}}$ = ${\rho _{\text{f}}}$ gd。

将式（2）代入式（3），消去w项，可以将控制方程改写为

$K\left( {\frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{{\partial ^2}p}}{{\partial {z^2}}}} \right) - {\gamma _{\text{w}}}n\beta \frac{{\partial p}}{{\partial t}} = {\gamma _{\text{w}}}\frac{\partial }{{\partial t}}\left( {\frac{{\partial {u_x}}}{{\partial x}} + \frac{{\partial {u_z}}}{{\partial z}}} \right) 。$

(5)

式中： ${\gamma _{\text{w}}}$ 为孔隙流体的重度。

根据式（1），采用 ${\sigma^{\prime}_{xx}}$ ， ${\sigma{zz}^{\prime}_}$ 和τ_xz分别表示沿x方向和z方向的有效应力以及剪切应力，那么多孔弹性海床的整体平衡方程可以表示为

$\frac{{\partial {{\sigma^{\prime}}_{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{xz}}}}{{\partial z}} = - \frac{{\partial p}}{{\partial x}} \text{，}$

(6)

$\frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {{\sigma^{\prime}}_{zz}}}}{{\partial z}} = - \frac{{\partial p}}{{\partial z}} 。$

(7)

多孔弹性海床的本构关系可以表示为

$\left.\begin{array}{l}{{\sigma }^{\prime }}_{xx}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\partial {u}_{x}}{\partial x}+\frac{\mu }{1-\mu }\frac{\partial {u}_{z}}{\partial z}\right)\text{ }\text{，}\\ {{\sigma }^{\prime }}_{zz}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\mu }{(1-\mu )}\frac{\partial {u}_{x}}{\partial x}+\frac{\partial {u}_{z}}{\partial z}\right)\text{ }\text{，}\\ {\tau }_{xz}=-G\left(\frac{\partial {u}_{x}}{\partial z}+\frac{\partial {u}_{z}}{\partial x}\right)\text{ }。\end{array}\right\}$

(8)

式中：μ为泊松比；G为剪切模量。

1.2 边界条件

（1）海床表面的边界条件（z=0）

如图 1所示，海床表面的垂直有效应力和剪切应力趋于零，且海洋表面产生波浪运动引起的周期荷载：

${\sigma{zz}^{\prime}_} = {\tau _{xz}} = 0 \text{，} {p_{\text{f}}} = {p_{\text{b}}}(x,t) 。$

(9)

海床表面的波浪压力可以进一步表示为

${p_{\text{b}}}(x,t) = \sum\limits_{m = 1}^2 {{P_m}\cos } (mkx - m\omega t) \text{，}$

(10)

其中

${P_1} = \frac{{{\rho _{\text{f}}}gH}}{{2\cosh (kd)}} \text{，}$

(11)

${P_2} = \frac{3}{4}{\gamma _{\text{w}}}H\left( {\frac{{{\text{π }}H}}{L}} \right)\frac{1}{{\sinh (2kd)}}\left[ {\frac{1}{{{{\sinh }^2}(kd)}} - \frac{1}{3}} \right] 。$

(12)

式中： ${p_{\text{b}}}$ 为海床表面的周期荷载；m=1, 2，P₁和P₂分别为海床表面的一阶和二阶压力；k=2π/L为波数，由波长L确定；ω=2π/T为波的角频率，由周期T确定。

（2）海床底面的边界条件（z=-h）

如图 1所示，假定有限厚度海床的底面为刚性、不透水基岩，底部边界条件可以表示为

${u_x} = {u_z} = 0 \text{，} {{\partial p} \mathord{\left/ {\vphantom {{\partial p} {\partial z}}} \right. } {\partial z}} = 0 。$

(13)

2. 海床动力响应求解

由于作用在海床表面的波浪荷载呈现周期性变化，因此未知场量同样随着波数k与角频率ω的变化产生周期性变化，可以基于Madsen^[14]求解过程中使用的复变量对该问题进行求解。同时，考虑到二阶Stokes波由不同阶次的波浪组成，因此可以将所有场量改写成由不同阶次求和的复变量形式：

$\left\{ {\begin{array}{*{20}{c}} {{u_x}(x,z,t)} \\ {{u_z}(x,z,t)} \\ {p(x,z,t)} \end{array}} \right\} = \sum\limits_{m = 1}^2 {{Re} \left\{ {\begin{array}{*{20}{c}} {U_x^{\left( m \right)}(z)} \\ {U_z^{\left( m \right)}(z)} \\ {{P^{\left( m \right)}}(z)} \end{array}} \right\}} {{\text{e}}^{{\text{i}}m\left( {kx - \omega t} \right)}} 。$

(14)

式中：i=(-1)^1/2为虚数；Re代表对结果取实部； $U_x^{\left( m \right)}$ 、 $U_z^{\left( m \right)}$ 和 ${P^{\left( m \right)}}$ 为波浪荷载引起的海床动力响应的幅值。值得注意的是，由于将场量拓展到复数域进行求解，因此仅取最后得到结果的实部进行分析。

波浪作用在海床表面的荷载在复数域可以表示为

${p_{\text{b}}}(x,t) = \sum\limits_{m = 1}^2 {{P_m}{{\text{e}}^{{\text{i}}m\left( {kx - \omega t} \right)}}} 。$

(15)

通过将场量表示为复变量形式，式（5），（6）和（8）可以分别改写为

$\frac{{{\partial ^2}{P^{\left( m \right)}}}}{{\partial {z^2}}} - \left( {{m^2}{k^2} - \frac{{{\gamma _{\text{w}}}n\beta }}{K}{\text{i}}m\omega } \right){P^{\left( m \right)}}$

$\begin{array}{*{20}{c}} {}&{}&{}&{} \end{array} = - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}}}{K}\left( {{\text{i}}mkU_x^{\left( m \right)} + \frac{{\partial U_z^{\left( m \right)}}}{{\partial z}}} \right) \text{，}$

(16)

${P^{\left( m \right)}} = - \left( {{{\sigma^{\prime}}_{xx}} + \frac{1}{{{\text{i}}mk}}\frac{{\partial {\tau _{xz}}}}{{\partial z}}} \right) \text{，}$

(17)

$\left.\begin{array}{l}{{\sigma }^{\prime }}_{xx}=-\frac{2(1-\mu )G}{1-2\mu }\left(\text{i}mk{U}_{x}^{\left(m\right)}+\frac{\mu }{1-\mu }\frac{\partial {U}_{z}^{\left(m\right)}}{\partial z}\right)\text{ }\text{，}\\ {{\sigma }^{\prime }}_{zz}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\mu }{1-\mu }\text{i}mk{U}_{x}^{\left(m\right)}+\frac{\partial {U}_{z}^{\left(m\right)}}{\partial z}\right)\text{ }\text{，}\\ {\tau }_{xz}=-G\left(\frac{\partial {U}_{x}^{\left(m\right)}}{\partial z}+\text{i}mk{U}_{z}^{\left(m\right)}\right)\text{ }。\end{array}\right\}$

(18)

将式（17）代入式（7）和式（16）可以分别得到

$\frac{{\partial {{\sigma^{\prime}}_{zz}}}}{{\partial z}} + {\text{i}}mk{\tau _{xz}} = \frac{{\partial {{\sigma^{\prime}}_{xx}}}}{{\partial z}} + \frac{1}{{{\text{i}}mk}}\frac{{{\partial ^2}{\tau _{xz}}}}{{\partial {z^2}}} \text{，}$

(19)

$\frac{{{\partial ^2}{{\sigma^{\prime}}_{xx}}}}{{\partial {z^2}}} + \frac{1}{{{\text{i}}mk}}\frac{{{\partial ^3}{\tau _{xz}}}}{{\partial {z^3}}} - \left( {{m^2}{k^2} - \frac{{{\gamma _{\text{w}}}n\beta }}{K}{\text{i}}m\omega } \right)\left( {{{\sigma^{\prime}}_{xx}} + \frac{1}{{{\text{i}}mk}}\frac{{\partial {\tau _{xz}}}}{{\partial z}}} \right) \\ \;\;\;\;\;\;\;\;\;\; = \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}}}{K}\left( {{\text{i}}mkU_x^{(m)} + \frac{{\partial U_z^{(m)}}}{{\partial z}}} \right) 。$

(20)

将多孔弹性介质本构关系式（18）代入式（19），可以得到

$\frac{{{\partial ^2}U_z^{(m)}}}{{\partial {z^2}}} - {m^2}{k^2}U_z^{(m)} = \frac{1}{{{\text{i}}mk}}\left( {\frac{{{\partial ^3}U_x^{(m)}}}{{\partial {z^3}}} - {m^2}{k^2}\frac{{\partial U_x^{(m)}}}{{\partial z}}} \right) 。$

(21)

将式（7）代入式（20），可以得到

$\begin{array}{l} \left\{ {\frac{{{\partial ^3}U_z^{\left( m \right)}}}{{\partial {z^3}}} - {m^2}{k^2}\frac{{\partial U_z^{\left( m \right)}}}{{\partial z}} + } \right. \\ \left. {\frac{{{\text{i}}m\omega {\gamma _{\text{w}}}\left[ {n\beta + (1 - 2\mu )/G} \right]}}{K}\frac{{\partial U_z^{\left( m \right)}}}{{\partial z}}} \right\} \\ = - \frac{1}{{{\text{i}}mk}}\left\{ {(1 - 2\mu )\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} - \left[ {2(1 - \mu ){m^2}{k^2} + (1 - 2\mu ) \cdot } \right.} \right. \\ {\text{ }}\left. {\left( {{m^2}{k^2} - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}n\beta }}{K}} \right)} \right]\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} + \\ {m^2}{k^2}\left[ {2(1 - \mu ){m^2}{k^2} - {\text{i}}m\omega {\gamma _{\text{w}}} \cdot } \right. \\ \end{array}$

$\left. {\left. {\frac{{2n\beta (1 - \mu ) + (1 - 2\mu )/G}}{K}} \right]U_x^{\left( m \right)}} \right\} 。$

(22)

将式（21）代入式（22），可以得到

$\begin{array}{l} \frac{{\partial U_z^{\left( m \right)}}}{{\partial z}} = \frac{{2(1 - \mu )K}}{{{m^2}\omega k{\gamma _{\text{w}}}\left[ {n\beta + (1 - 2\mu )/G} \right]}} \times \\ \;\;\;\;\;\;\;\left\{ {\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} - \left[ {2{m^2}{k^2} - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}n\beta }}{K}\frac{{(1 - 2\mu )}}{{2(1 - \mu )}}} \right]\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} + } \right. \end{array} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{m^2}{k^2}({m^2}{k^2} - {\kappa ^2})U_x^{\left( m \right)}\left. \begin{gathered} \\ \\ \end{gathered} \right\} 。$

(23)

式中：

${\kappa ^2} = \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}\left( {n\beta + \frac{{(1 - 2\mu )}}{{2G(1 - \mu )}}} \right)}}{K} 。$

(24)

将式（23）代入式（22），可以得到

$\frac{{{\partial ^6}U_x^{\left( m \right)}}}{{\partial {z^6}}} - {m^2}{k^2}\left( {3 - \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} +$

${m^4}{k^4}\left( {3 - 2\frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} - {m^6}{k^6}\left( {1 - \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)U_x^{\left( m \right)} = 0 。$

(25)

式（25）为六阶偏微分方程，该方程的解可以假设为

$U_x^{\left( m \right)} = A{{\text{e}}^{\gamma mkz}} 。$

(26)

将所假设的解式（26）代入式（25）中，可以得到

$({\gamma ^2} - 1)({\gamma ^2} - 1)\left( {{\gamma ^2} - 1 + \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right) = 0 。$

(27)

由式（27），可以得到 $U_x^{\left( m \right)}$ 的解析解为

${U}_{x}^{\left(m\right)}=({A}_{1}+{A}_{2}z){\text{e}}^{mkz}+({A}_{3}+{A}_{4}z){\text{e}}^{-mkz}+{A}_{5}{\text{e}}^{{k}_{*}mkz}+{A}_{6}{\text{e}}^{-{k}_{*}mkz}。$

(28)

式中： ${k_*}$ =[1-κ²/(mk)²]。

考虑到后续计算中矩阵求逆的正指数溢出问题，通过将式（28）重构为

$U_x^{\left( m \right)} = ({A_1} + {A_2}z){{\text{e}}^{mkz}} + {A_5}{{\text{e}}^{{k_*}mkz}} + ({A_3} + {A_4}z){{\text{e}}^{ - mk\left( {z + h} \right)}} +$

${A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}} 。$

(29)

将式（29）代入式（23），并对z求积分，可以得到 $U_z^{\left( m \right)}$ 的解析解为

$U_z^{\left( m \right)} = - {\text{i}}\left\{ {\left( {{A_1} + {A_2}z - \frac{{1 + 2\delta }}{{mk}}{A_2}} \right){{\text{e}}^{mkz}}{\text{ + }}{k_*}{A_5}{{\text{e}}^{{k_*}mkz}} - } \right.$

$\left. {\left( {{A_3} + {A_4}z + \frac{{1 + 2\delta }}{{mk}}{A_4}} \right){{\text{e}}^{ - mk\left( {z + h} \right)}} - {k_*}{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} 。$

(30)

式中：

$\delta = \frac{{(1 - 2\mu )n\beta }}{{n\beta + {{(1 - 2\mu )} \mathord{\left/ {\vphantom {{(1 - 2\mu )} G}} \right. } G}}} 。$

(31)

同理可得孔隙水压力和有效应力的解析解为

${P^{\left( m \right)}} = 2G{\text{i}}\left\{ {\frac{{(\delta + 2\mu - 1)}}{{1 - 2\mu }}{A_2}{{\text{e}}^{mkz}} + \left[ {\frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }}} \right]mk{A_5}{{\text{e}}^{{k_*}mkz}} + } \right. \\ \left. {\frac{{1 - 2\mu - \delta }}{{1 - 2\mu }}{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(32)

$\sigma_{x x}^{\prime(m)} = - 2G{\text{i}}\left\{ {mk{A_1}{{\text{e}}^{mkz}} + \frac{{(1 - 2\mu )mkz + 2\mu \delta }}{{1 - 2\mu }}{A_2}{{\text{e}}^{mkz}} + } \right.$

${\text{ }}\frac{{1 - \mu - \mu k_*^2}}{{1 - 2\mu }}mk{A_5}{{\text{e}}^{{k_*}mkz}} + mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} +$

${\text{ }}\frac{{(1 - 2\mu )mkz - 2\mu \delta }}{{1 - 2\mu }}{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} +$

$\left. {\frac{{1 - \mu - \mu k_*^2}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(33)

$\begin{gathered} \sigma_{z z}^{\prime(m)} = 2G{\text{i}}\left\{ {mk{A_1}{{\text{e}}^{mkz}} - \left[ { - mkz + \frac{{2(1 - \mu )}}{{1 - 2\mu }}\delta } \right]{A_2}{{\text{e}}^{mkz}} - } \right. \\ {\text{ }}\frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }}mk{A_5}{{\text{e}}^{{k_*}mkz}} + mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \\ {\text{ }}\left[ {mkz + \frac{{2\delta (1 - \mu )}}{{1 - 2\mu }}} \right]{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} - \\ \end{gathered}$

$\left. {\frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(34)

$\begin{gathered} \tau _{xz}^{\left( m \right)} = 2G\left\{ { - mk{A_1}{{\text{e}}^{mkz}} - (mkz - \delta ){A_2}{{\text{e}}^{mkz}} - {k_*}mk{A_5}{{\text{e}}^{{k_*}mkz}}} \right. + \\ {\text{ }}mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} + (mkz + \delta ){A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \end{gathered} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \left. {{k_*}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} 。$

(35)

式中：未知参数A₁~A₆可以代入边界条件（式（9），（13））确定，如下式所示：

$\left[\begin{array}{c}{A}_{1}\\ {A}_{2}\\ {A}_{3}\\ {A}_{4}\\ {A}_{5}\\ {A}_{6}\end{array}\right]={\left[\begin{array}{cccccc}{B}_{1}& -{B}_{2}& {B}_{1}{B}_{5}& {B}_{2}{B}_{5}& -{B}_{3}{B}_{1}& -{B}_{3}{B}_{1}{B}_{10}\\ -{B}_{1}& \delta & {B}_{1}{B}_{5}& \delta {B}_{5}& -{B}_{4}& {B}_{4}{B}_{10}\\ {B}_{5}& -{B}_{5}h& 1& -h& {B}_{10}& 1\\ {B}_{5}& -{B}_{5}h-{B}_{7}{B}_{5}& -1& h-{B}_{7}& {k}_{*}{B}_{10}& -{k}_{*}\\ 0& {B}_{8}{B}_{5}& 0& {B}_{8}& {B}_{9}{B}_{4}{B}_{10}& -{B}_{9}{B}_{4}\\ 0& 2G\text{i}{B}_{8}& 0& -2G\text{i}{B}_{8}{B}_{5}& 2G\text{i}{B}_{9}{B}_{1}& 2G\text{i}{B}_{9}{B}_{1}{B}_{10}\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ {p}_{b}\end{array}\right]。$

(36)

式中：

$\begin{array}{l} {B_1} = mk;{B_2} = \frac{{2\delta (1 - \mu )}}{{1 - 2\mu }};{B_3} = \frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }};{B_4} = {k_*}mk; \\ {B_5} = {{\text{e}}^{ - mkh}};{B_6} = {{\text{e}}^{mkh}};{B_7} = \frac{{1 + 2\delta }}{{mk}};{B_8} = \frac{{\delta + 2\mu - 1}}{{1 - 2\mu }}; \end{array} \\ \;\;\;\; {B_9} = \frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }};{B_{10}} = {{\text{e}}^{ - {k_*}mkh}};{B_{11}} = {{\text{e}}^{{k_*}mkh}} 。$

(37)

3. 二阶Stokes波作用下海床动力响应分析

3.1 与既有解的对比

为了对比验证本文得到的二阶Stokes波作用下海床动力响应的解析解，本文分别给出了线性波和二阶Stokes波作用下的对比结果，对比结果中实线代表本文解析解，点代表既有解。首先通过选取与Hsu等^[13]一致的参数：T=15 s，d=70 m，L=324 m，h=25 m，G=10⁷ Pa，n=0.3，μ=1/3以及K=10^-2 m/s，将本文解析解与Hsu等^[13]得到的线性波作用下的解进行对比。由图 2可知本文解析解得到的孔隙水压力和有效应力曲线与Hsu等的曲线完全一致，从而可以证明本文解析解的正确与可靠。

图 2 线性波作用下本文解析解与Hsu等的解之间的对比

Figure 2. Comparison of seabed response between present solution and solution by Hsu et al. under linear waves

下载: 全尺寸图片幻灯片

进一步地，将本文解析解与Zhang等^[20]二阶Stokes波作用下的解进行对比如图 3所示。对比中选用的参数与Zhang等一致：T=12 s，H=8 m，d=0.125L，h=24 m，S_r=0.975，n=0.35，μ=0.4，K=10^-4 m/s以及G=5×10⁶ Pa。对比结果表明本文解析解与Zhang等^[20]的解有良好的匹配性，进一步表明了本文解的正确性和可靠性。

图 3 二阶斯托克斯波作用下本文解析解与Zhang等的解之间的对比

Figure 3. Comparison of seabed response between present solution and solution by Zhang et al. under second-order Stokes waves

下载: 全尺寸图片幻灯片

3.2 波浪和海床参数对动力响应的影响

由二阶Stokes波作用在海床表面周期压力的表达式（10）~（12）可知，波浪的周期T、波高H以及水深d对二阶Stokes波的影响较大；而在对海床动力响应的求解过程中可知，海床土体的剪切模量G和渗透系数K同样会对海床的动力响应产生影响。为了探究波浪参数和海床土体参数对海床动力响应的影响，给出了如表 1所示的波浪和海床土体参数并进行算例计算和分析。在下列算例计算和分析中，没有提及的参数的取值均来源于表 1。

表 1 多孔弹性海床和波浪的基本参数

Table 1. Basic poroelastic properties and wave conditions

波浪和海床土参数	数值
波的周期T	10 s
波高H	5 m
水深d	0.125L^*
海床厚度h	20 m
孔隙流体的实际体积弹性模量K_w	2×10⁹ Pa
孔隙流体密度ρ_f	1000 kg/m³
孔隙率n	0.35
泊松比μ	0.33
渗透系数K	10^-4 m/s
剪切模量G	1×10⁷ Pa
海床土饱和度S_r	0.975
^*表中L代表波长。

下载: 导出CSV

| 显示表格

（1）波浪参数的影响

首先对波浪参数对海床动力响应的影响进行研究。图 4是二阶Stokes波作用下波浪周期T对海床动力响应的影响。其中波浪周期T的取值分别为8，9，10，11，12 s。从图 4中可以发现不同的波浪周期T对海床动力响应随深度方向z/h的变化趋势是一致的，还可以从图中发现：在海床相同深度z/h时，最大孔隙水压力|p|随着波浪周期T的增加而增加；最大水平有效应力|σ′_xx|在较浅区域随着波浪周期T的增加而减小，在较深区域则呈现相反的趋势；最大垂直有效应力|σ′_zz|和最大剪切应力|τ_xz|在较浅区域随着波浪周期T的增加而减小，在较深区域则存在转折点。同时还可以看出在相对水深d/L不变的情况下，波浪周期T越小，海床的动力响应受到二阶Stokes波的二阶项的影响越明显。

图 4 二阶Stokes波作用下波浪周期T对海床动力响应的影响

Figure 4. Influences of wave period on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 5是二阶Stokes波作用下水深d对海床动力响应的影响。Le Méhauté^[22]认为当水深d小于0.1L时，Stokes波不能严格满足自由水面的边界条件，故水深d的取值分别为：0.1L，0.12L，0.125L，0.13L，0.14L。从图 5可以得知水深d对海床的动力响应有明显的影响，同时水深d还对二阶Stokes波的二阶项引起的海床动力响应的大小有相当显著的影响。与波浪周期T的影响类似，海床动力响应随深度方向z/h的变化趋势是一致的。最大孔隙水压力|p|、最大水平有效应力|σ′_xx|、最大垂直有效应力|σ′_zz|和最大剪切应力|τ_xz|随着水深d的增加与波浪周期T增加时均表现出相似的规律。当波浪周期T不变时，随着相对水深d/L较小时（例如d/L=0.1），海床的动力响应受二阶Stokes波二阶项的影响有显著的影响。

图 5 二阶Stokes波作用下水深d对海床动力响应的影响

Figure 5. Influences of water depth on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 6是二阶Stokes波作用下波高H对海床动力响应的影响。其中波高H的取值分别为3，4，5，6，7 m。由图 6可以得知波高H对海床动力响应的影响较小。在海床的相同深度z/h处，最大孔隙水压力|p|、最大垂直有效应|σ′_zz|和最大剪切应力|τ_xz|均随着波高H的增大而增大。同时，波高H对二阶Stokes波二阶项引起的海床动力响应基本没有影响。

图 6 二阶Stokes波作用下波高H对海床动力响应的影响

Figure 6. Influences of seabed thickness on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

（2）海床土体参数对海床动力响应的影响

图 7是二阶Stokes波作用下渗透系数K对海床动力响应的影响。其中渗透系数K的取值分别为10^-2，10^-3，10^-4，10^-5，10^-6 m/s。由图 7可得：渗透系数K对海床的动力响应有显著的影响。渗透系数K对海床动力响应的影响体现在沿深度方向z/h的变化率上，随着渗透系数K的减小在深度较浅处的最大孔隙水压力|p|减小速率越快，最大垂直有效应力|σ′_zz|增大速率越快，而在较深区域这些场量的变化速率与渗透系数K的相关性减小。当渗透系数K较大时，最大水平有效应力|σ′_xx|和最大剪切应力|τ_xz|受到的影响较大，渗透系数K较小时对最大剪切应力|τ_xz|基本没有影响，且对最大水平有效应力|σ′_xx|的影响较小。

图 7 二阶Stokes波作用下渗透系数K对海床动力响应的影响

Figure 7. Influences of permeability on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 8是二阶Stokes波作用下剪切模量G对海床动力响应的影响。其中剪切模量G的取值分别为1×10⁶，5×10⁶，1×10⁷，5×10⁷，1×10⁸ Pa。从图 8中可以得出：剪切模量G对海床的动力响应有显著的影响。在海床的相同深度z/h处，最大孔隙水压力|p|和最大垂直有效应力|σ′_zz|随着剪切模量G的增大而分别减小和增大。而最大剪切应力|τ_xz|在深度z/h大约为0.5处存在一个转折点，在转折点以上的相同深度z/h处随着剪切模量G的增大而增大，在转折点以下则呈现相反的变化趋势，此现象可能是由于竖向和水平向有效应力状态改变引起的，即在z/h=0.5以下部分，水平有效应力出现了快速的增大而竖直有效应力出现了减小，在竖向偏应力状态改变的情况下，剪切应力状态呈现了相反的变化趋势。

图 8 二阶Stokes波作用下剪切模量G对海床动力响应的影响

Figure 8. Influences of shear modulus on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

4. 结论

本文研究引入了二阶Stokes波理论模拟波浪荷载，考虑了波浪非线性的影响，最终建立了二阶Stokes波作用下有限厚度海床动力响应的解析解。并将本文得到的解析解与线性波既有解进行了对比，验证了本文解的正确性和可靠性。最后分析了波浪和海床土体参数对海床动力响应的影响。通过本文的研究主要得出以下3点结论。

（1）在相对水深d/L确定，波浪周期T较小时，以及波浪周期T确定，相对水深较小时，海床的动力响应受到二阶Stokes波的二阶项影响均增加。

（2）波浪周期T、水深d和波高H均会对海床的动力响应产生影响，但波高H对二阶Stokes波的二阶项引起的海床动力响应基本没有影响。

（3）海床土体的渗透系数K和剪切模量G均会对海床的动力响应有显著的影响。渗透系数K对海床动力响应的影响表现为沿深度方向的变化速率，而剪切模量G对海床动力响应体现为幅值大小的影响。

图 1 地铁车站二维整体有限元模型示意图

Figure 1. 2D finite element model for soil-underground structure

下载: 全尺寸图片幻灯片

图 2 基岩场地地震动加速度反应谱(ξ=5%)

Figure 2. Response spectra of acceleration of bedrock ground motions (ξ=5%)

下载: 全尺寸图片幻灯片

图 3 地铁车站IDA分析结果

Figure 3. IDA results of subway station

下载: 全尺寸图片幻灯片

图 4 地铁车站IDR-PGA线性回归

Figure 4. Linear regression PGA for IDR of subway station based on IDA

下载: 全尺寸图片幻灯片

图 5 Ⅲ类场地中地铁车站易损性曲线

Figure 5. Curves of seismic fragility of subway station in class Ⅲ site

下载: 全尺寸图片幻灯片

图 6 所选地震动反应谱的分位值

Figure 6. Quantile values of response spectra of selected ground motions

下载: 全尺寸图片幻灯片

图 7 基于反应谱分位值构造的ETAFs

Figure 7. ETAFs based on quantile values of response spectra of ground motions

下载: 全尺寸图片幻灯片

图 8 ETA分位值结果与IDA结果对比

Figure 8. Comparison of quantile value results between ETA and IDA

下载: 全尺寸图片幻灯片

图 9 ETA与IDA各分位值结果相关性

Figure 9. Correlation of quantile value results between ETA and IDA

下载: 全尺寸图片幻灯片

图 10 地铁车站ETA分析结果

Figure 10. ETA results of subway station

下载: 全尺寸图片幻灯片

图 11 地铁车站IDR-PGA线性回归

Figure 11. Linear regression PGA for IDR of subway station based on ETA

下载: 全尺寸图片幻灯片

图 12 ETA易损性分析结果验证

Figure 12. Curves of seismic fragility obtained by ETA and IDA

下载: 全尺寸图片幻灯片

表 1 基于层间位移角的性能指标限值

Table 1 Performance index system based on interstory drift ratio

指标性能水准

1 2 3 4

层间位移角 1/1223 1/343 1/161 1/105

下载: 导出CSV

参考文献(22)

[1]	韩建平, 吕西林, 李慧. 基于性能的地震工程研究的新进展及对结构非线性分析的要求[J]. 地震工程与工程振动, 2007, 27(4): 15-23. doi: 10.3969/j.issn.1000-1301.2007.04.003 HAN Jianping, LÜ Xilin, LI Hui. State-of-the-art of performance-based earthquake engineering and need for structural nonlinear analysis[J]. Journal of Earthquake Engineering and Engineering Vibration, 2007, 27(4): 15-23. (in Chinese) doi: 10.3969/j.issn.1000-1301.2007.04.003
[2]	崔臻, 盛谦, 冷先伦, 等. 基于增量动力分析的大型地下洞室群性能化地震动力稳定性评估[J]. 岩石力学与工程学报, 2012, 31(4): 703-712. doi: 10.3969/j.issn.1000-6915.2012.04.009 CUI Zhen, SHENG Qian, LENG Xianlun, et al. Performance-based seismic stability assessment of large underground cavern group with incremental dynamic analysis[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(4): 703-712. (in Chinese) doi: 10.3969/j.issn.1000-6915.2012.04.009
[3]	LENG X L, CUI Z, LIU J J, et al. Performance-based seismic stability assessment of large underground carven complex with incremental dynamic analysis[J]. Disaster Advances, 2012, 5(4): 686-692.
[4]	LIU T, CHEN Z Y, YUAN Y, et al. Fragility analysis of a subway station structure by incremental dynamic analysis[J]. Advances in Structural Engineering, 2017, 20(7): 1111-1124. doi: 10.1177/1369433216671319
[5]	ZHONG Z L, SHEN Y Y, ZHAO M, et al. Seismic fragility assessment of the Daikai subway station in layered soil[J]. Soil Dynamics and Earthquake Engineering, 2020, 132: 106044. doi: 10.1016/j.soildyn.2020.106044
[6]	钟紫蓝, 甄立斌, 申轶尧, 等. 基于耐震时程分析法的地下结构抗震性能评价[J]. 岩土工程学报, 2020, 42(8): 1482-1490. doi: 10.11779/CJGE202008013 ZHONG Zilan, ZHEN Libin, SHEN Yiyao, et al. Seismic performance evaluation of underground structures using endurance time analysis[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(8): 1482-1490. (in Chinese) doi: 10.11779/CJGE202008013
[7]	钟紫蓝, 甄立斌, 陈群, 等. 浅埋地下结构耐震时程分析法最优持时研究[J]. 岩土工程学报, 2021, 43(4): 698-705. doi: 10.11779/CJGE202104012 ZHONG Zilan, ZHEN Libin, CHEN Qun, et al. Optimal duration of endurance time acceleration functions for shallow buried underground structures[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(4): 698-705. (in Chinese) doi: 10.11779/CJGE202104012
[8]	RAHIMI E, ESTEKANCHI H E. Collapse assessment of steel moment frames using endurance time method[J]. Earthquake Engineering and Engineering Vibration, 2015, 14(2): 347-360. doi: 10.1007/s11803-015-0027-0
[9]	TAVAZO H A, RANJBARAN A. Fragility analysis of 3D reinforced concrete frames based on endurance time method with derived standard deviation[J]. Journal of Earthquake and Tsunami, 2017, 11(4): 1750011. doi: 10.1142/S1793431117500117
[10]	HARIRI-ARDEBILI M A, SATTAR S, ESTEKANCHI H E. Performance-based seismic assessment of steel frames using endurance time analysis[J]. Engineering Structures, 2014, 69: 216-234. doi: 10.1016/j.engstruct.2014.03.019
[11]	何海峰. 基于耐震时程法的冲刷环境下桥梁结构地震易损性分析[D]. 成都: 西南交通大学, 2020. HE Haifeng. Seismic Fragility Analysis of Bridges Considering Scour Effect Based on Endurance Time Method[D]. Chengdu: Southwest Jiaotong University, 2020. (in Chinese)
[12]	徐舒桐. 基于耐震时程法的混凝土坝易损性分析[D]. 大连: 大连理工大学, 2019. XU Shutong. Seismic Fragility Analysis of Concrete Dam Based on Endurance Time Analysis Method[D]. Dalian: Dalian University of Technology, 2019. (in Chinese
[13]	朱少坤. 基于耐震时程分析方法的水电站厂房易损性分析[D]. 西安: 西安理工大学, 2020. ZHU Shaokun. Study on Seismic Behavior of Hydropower Plant Based on Seismic Time-History Analysis[D]. Xi'an: Xi'an University of Technology, 2020. (in Chinese)
[14]	VAMVATSIKOS D, ALLIN CORNELL C. Incremental dynamic analysis[J]. Earthquake Engineering & Structural Dynamics, 2002, 31(3): 491-514.
[15]	VAMVATSIKOS D, ALLIN CORNELL C. Applied incremental dynamic analysis[J]. Earthquake Spectra, 2004, 20(2): 523-553. doi: 10.1193/1.1737737
[16]	陈建宝, 丁军军. 分位数回归技术综述[J]. 统计与信息论坛, 2008, 23(3): 89-96. https://www.cnki.com.cn/Article/CJFDTOTAL-TJLT200803019.htm CHEN Jianbao, DING Junjun. A review of technologies on quantile regression[J]. Statistics & Information Forum, 2008, 23(3): 89-96. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TJLT200803019.htm
[17]	陈国兴, 庄海洋, 杜修力, 等. 土-地铁车站结构动力相互作用大型振动台模型试验研究[J]. 地震工程与工程振动, 2007, 27(2): 171-176. https://www.cnki.com.cn/Article/CJFDTOTAL-DGGC200702026.htm CHEN Guoxing, ZHUANG Haiyang, DU Xiuli, et al. Analysis of large-scale shaking table test of dynamic soil-subway station interaction[J]. Journal of Earthquake Engineering and Engineering Vibration, 2007, 27(2): 171-176. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DGGC200702026.htm
[18]	钟紫蓝, 申轶尧, 甄立斌, 等. 地震动强度参数与地铁车站结构动力响应指标分析[J]. 岩土工程学报, 2020, 42(3): 486-494. doi: 10.11779/CJGE202003010 ZHONG Zilan, SHEN Yiyao, ZHEN Libin, et al. Ground motion intensity measures and dynamic response indexes of metro station structures[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(3): 486-494. (in Chinese) doi: 10.11779/CJGE202003010
[19]	钟紫蓝, 申轶尧, 郝亚茹, 等. 基于IDA方法的两层三跨地铁地下结构地震易损性分析[J]. 岩土工程学报, 2020, 42(5): 916-924. doi: 10.11779/CJGE202005014 ZHONG Zilan, SHEN Yiyao, HAO Yaru, et al. Seismic fragility analysis of two-story and three-span metro station structures based on IDA method[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(5): 916-924. (in Chinese) doi: 10.11779/CJGE202005014
[20]	杜修力, 蒋家卫, 许紫刚, 等. 浅埋矩形框架地铁车站结构抗震性能指标标定研究[J]. 土木工程学报, 2019, 52(10): 111-119, 128. https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201910012.htm DU Xiuli, JIANG Jiawei, XU Zigang, et al. Study on quantification of seismic performance index for rectangular frame subway station structure[J]. China Civil Engineering Journal, 2019, 52(10): 111-119, 128. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201910012.htm
[21]	ARGYROUDIS S A, PITILAKIS K D. Seismic fragility curves of shallow tunnels in alluvial deposits[J]. Soil Dynamics and Earthquake Engineering, 2012, 35: 1-12.
[22]	NIBS. HAZUS-MH: Technical Manuals[M]. Washington D C: Federal Emergency Management Agency and National Institute of Building Science, 2004.