基于等效零时刻点的水平层状场地黏弹性边界方法研究

贾忠明; 赵伯明; 王子珺; 赵天次; 冯俊

doi:10.11779/CJGE202212020

基于等效零时刻点的水平层状场地黏弹性边界方法研究 English Version

贾忠明^{1, 2,},
赵伯明^{1, 2, ,},
王子珺^{1, 2},
赵天次^{1, 2},
冯俊^{1, 2}

1.
北京交通大学城市地下工程教育部重点实验室，北京 100044
2.
北京交通大学土木建筑工程学院，北京 100044

基金项目:

国家自然科学基金项目 51778046

国家自然科学基金项目 51808035

国家自然科学基金重点项目 U1434210

详细信息

作者简介:
贾忠明(1996—)，男，博士研究生，主要从事结构抗震方面的研究。E-mail：zhmjia@bjtu.edu.cn

通讯作者:
赵伯明, E-mail: bmzhao@bjtu.edu.cn

中图分类号: TU435
计量
- 文章访问数: 118
- HTML全文浏览量: 34
- PDF下载量: 21
出版历程
- 收稿日期: 2021-10-25
- 网络出版日期: 2022-12-13
- 刊出日期: 2022-11-30

Viscoelastic boundary method of horizontal layered site based on equivalent zero-time point

JIA Zhong-ming^{1, 2,},
ZHAO Bo-ming^{1, 2, ,},
WANG Zi-jun^{1, 2},
ZHAO Tian-ci^{1, 2},
FENG Jun^{1, 2}

1.
Key Laboratory of Urban Underground Engineering, Ministry of Education, Beijing Jiaotong University, Beijing 100044, China
2.
School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China

摘要

摘要: 针对黏弹性边界下水平层状场地地震动输入中的时间延迟计算问题，提出了水平层状场地的等效零时刻点迭代时间延迟计算方法，并依托斯涅尔方程及地震波在分层界面上的连续性条件推导了人工边界上等效节点力的计算公式，得出了黏弹性边界下的水平层状场地斜入射地震动输入方法。建立两层地层三维有限元计算模型，进行动力响应分析，并将计算结果与理论值进行对比，验证了方法的有效性和实用性。应用方法进行等效节点力计算的过程中，可根据所求目标地震波的具体信息及其所在地层的具体条件直接求解，无需对所求地震波的全部来源依次进行分析，使得公式更为统一，在简化计算的同时保证计算精度及准确性，并拓展了黏弹性边界理论在水平层状场地斜入射地震动上的应用。
- 地震波斜入射 /
- 水平层状场地 /
- 黏弹性边界 /
- 时间延迟
Abstract: An equivalent zero-time point iterative time delay method for horizontal layered strata is proposed to solve the problem of time delay calculation in ground motion input of layered strata under viscoelastic boundary. It assignes an equivalent zero-time plane to each generated seismic wave. An appropriate point on the plane can be selected as the equivalent zero-time point, and the time delay can be calculated. Based on the Snell equation and the continuity condition of seismic waves at the layered interface, the formula for calculating the equivalent nodal force of the artificial boundary is derived, and the viscoelastic boundary input method for oblique incidence ground motion of the layered strata under the viscoelastic boundary is deduced. The finite element model for a double-layer horizontal layered stratum is established for dynamic response analysis. The calculated results are compared with the theoretical values. The deviation between them is relatively small, which shows the effectiveness and practicability of the method. By using the method it can be solved directly according to the specific information of the seismic wave and the stratum, without analyzing all the sources of the desired seismic wave in turn, which makes the formula unified. It is a way to simplify the calculation without affecting the accuracy. Meanwhile, it can extend the application of the viscoelastic boundary to the oblique incidence ground motion of layered stratum.
- oblique incidence of seismic wave /
- layered stratum /
- viscoelastic boundary /
- time delay

HTML全文

0. 引言

中国地处环太平洋地震带以及欧亚地震带之间，其地震活动呈现出高频性、高强度、浅震源等特点，当地震发生在人口稠密的地区时，会带来惨重的经济和社会损失。对工程进行抗震研究，减轻地震灾害的影响将变得十分重要。

在对工程进行地震响应分析过程中，考虑到现实情况的复杂性，常采用数值方法进行计算^[1]。动力时程分析法可以较好地揭示结构在地震全时段内的响应规律^[2]，随着计算机硬件设备的不断发展，该方法得到了一定的认可。应用有限元方法进行动力时程分析过程中，为降低计算成本，提高计算效率，往往要将地层取出规定大小的有限区域来进行计算，并在此区域外施加人工边界来模拟无限地基的情况^[3]。其中人工边界包括无限元边界、黏性边界和黏弹性边界等形式，而黏弹性人工边界因具有好的低频和高频稳定性，应用方便，被广泛使用^[4]。

在进行工程抗震分析时，根据波动理论，考虑到地壳密度随地层深度的增加而增加，普遍认为地震波在地表附近的传播方向接近竖向^[5]，故而常常探讨地震波在垂直入射下情况。但现实情况下，根据地层条件的不同，地震波往往存在一定的入射角度。Jin等^[6]和Takahiro等^[7]通过对各地震动记录进行反演分析，得出了不同的地震波入射角度。因此，若仅考虑垂直入射的情况对结构进行地震响应分析，或将带来不能忽视的误差^[8]。因此，考虑地震波斜入射也具有一定的现实意义。而在此方面，许多学者^[9-11]同样进行了很多相关的研究。

而对于水平层状场地地震波斜入射问题，刘晶波等^[12-13]提出了一维化时域有限元计算方法；赵密等^[14-15]和尹侯权等^[16]利用弹性介质的应力–位移本构关系建立了人工边界处应力与速度的阻抗边界条件，替代黏性边界条件，改进P-SV波斜入射时成层半空间自由场的时域算法；王笃国等^[17-18]应用频域传递矩阵结合等效线性化方法，求解了地震波斜入射下水平层状场地的地震响应。而直接运用黏弹性边界理论，因时间延迟计算较为繁琐，缺乏深入研究。

对于地震波斜入射中时间延迟的计算，目前普遍采用立体几何的方式对地震波路径逐一进行求解，或采用坐标变换^[19]的方式对计算过程进行一定成度的简化。然而，当面对水平层状场地等特殊情况时，存在地震波多次反射和透射的情况，此方法不易迭代，计算成本较高。

对此，本文尝试建立一种新的时间延迟计算方法，以当前地震波属性表征其生成路径全过程上的时间延迟；在此基础上，推导黏弹性边界下水平层状场地地震波斜入射方法，进而实现合理有效的地震动数值计算。最后建立有限元计算模型，利用此方法对单脉冲作用下的数值解与理论解进行对比，验证方法的准确性和适用性。

1. 黏弹性边界地震动输入

黏弹性人工边界的思想就是将弹簧和黏滞阻尼器并联，形成一种弹簧-阻尼器物理元件，并将其施加在边界上。要使用这种边界，首先要将地震波转化成人工边界上各处所对应的等效荷载，其中包括产生内行场反应所需要抵抗的人工边界物理元件上的结点力，以及因产生内行场反应而需抵抗近场介质所需的结点力^[20-21]。当黏弹性边界可以将计算区域内向外传播的散射波完全吸收时，人工边界的节点上便会形成因地震作用而产生的自由场运动。

这样，地震动的输入问题便转化为在人工边界的节点上所作用的自由场运动问题，通等效节点力的方式施加在所求模型的边界上。

等效节点力的计算公式为

${F_{{\text{b}}}} = ({{{\boldsymbol{ K}}}_{{\text{b}}}}{{\boldsymbol{u }}_{{\text{b}}}} + {{{\boldsymbol{C }}}_{{\text{b}}}}{\dot {\boldsymbol{ u}}_{{\text{b}}}} + {\sigma _{{\text{b}}}}n){A_{{\text{l}}}} 。$

(1)

在弹簧刚度矩阵 ${{\boldsymbol{K }}_{\text{b}}}$ 和阻尼系数矩阵 ${{\boldsymbol{ C}}_{\text{b}}}$ 的取值方面，许多学者都进行了研究，同时给出了一定的建议，但不同方法对结果影响不大。对此，本文选用文献[22]的方法进行计算。

在水平层状场地的计算过程中，式(1)的自由场位移向量 ${{\boldsymbol{u }}_{{\text{b}}}}$ 和自由场速度向量 ${\dot {\boldsymbol{ u}}_{{\text{b}}}}$ 需要根据斯涅尔方程得到地震波在成层界面上的转换情况，并分析其与入射波之间的时空关系，进而利用波的叠加原理进行求解，这其中便涉及到时间延迟的计算问题。在此方面，当前的计算方法较为复杂，求解较为困难。下面针对水平层状场地时间延迟的计算进行分析。

2. 时间延迟计算方法

针对水平层状场地，本文假定各层地层均匀弹性，在水平方向上无限延伸。根据地震波波动理论可知，平面P波斜入射时，在地层界面处生成反射P波、反射SV波、透射P波和透射SV波；平面SV波斜入射时，如入射角度不超过临界角，将生成反射P波、反射SV波、透射P波和透射SV波。因此，为将时间延迟计算公式进行统一，在这里假设入射波为m，反射波n，透射波l，其中m，n，l为平面P波或平面SV波，入射波、反射波及透射波的方向角分别为 $\alpha$ ， $\beta$ ， $\gamma$ ， ${T_0}$ 时刻入射波波前与入射波传播方向垂直，仅区分反射和透射两种情况进行讨论。

对反射波，如所示， $\alpha$ 和 $\beta$ 满足：

$\frac{{{{\text{sin}}}\alpha }}{{{v_1}}} = \frac{{{{\text{sin}}}\beta }}{{{v_2}}}{\text{，}}$

(2)

图 1 反射波时间延迟计算示意图

Figure 1. Time delay calculation of reflected waves

下载: 全尺寸图片幻灯片

式中， ${v_1}$ ， ${v_2}$ 分别为入射波m和反射波n在对应地层内的传播速度。

假设 ${X_{{\text{f}}}}({x_0}, {y_0})$ 为所求模型边界上的一点，AB为零时刻的波前面，则反射波n在该点处的时间延迟为

$\Delta t = {t_1} + {t_2} {\text{，}}$

(3)

${t_1} = \frac{{OD}}{{{v_1}}}{\text{，}}$

(4)

${t_2} = \frac{{O{X_{{\text{f}}}}}}{{{v_2}}}{\text{，}}$

(5)

式中， ${t_1}$ ， ${t_2}$ 分别为计算所需入射波m和反射波n的时间延迟。

在此，过A点做ON垂线，与ON延长线交于点C。由几何关系得

$OD = AO\sin \alpha {\text{，}}$

(6)

$OC = AO\sin \beta {\text{，}}$

(7)

${t_1} = \frac{{OD}}{{{v_1}}} = \frac{{OC}}{{{v_2}}}{\text{，}}$

(8)

$\Delta t = {t_1} + {t_2} = \frac{{OC + O{X_{{\text{f}}}}}}{{{v_2}}} = \frac{{C{X_{{\text{f}}}}}}{{{v_2}}}。$

(9)

设点 $U({x_1}, {y_1}) \in AC$ ，则式（9）可改写为

$\Delta t = \frac{1}{{{v_2}}}\left[ {({x_1} - {x_0})\cos \beta + ({y_1} - {y_0})\sin \beta } \right]。$

(10)

对透射波，如所示， $\alpha$ 和 $\gamma$ 满足：

$\frac{{{{\text{sin}}}\alpha }}{{{v_1}}} = \frac{{{{\text{sin}}}\gamma }}{{{v_3}}}{\text{，}}$

(11)

图 2 透射波时间延迟计算示意图

Figure 2. Time delay calculation of transmission waves

下载: 全尺寸图片幻灯片

式中， ${v_3}$ 为透射波l在对应地层内的传播速度。

假设 ${X_t}({x_2}, {y_2})$ 为所求模型边界上的一点，AB为零时刻的波前面，则透射波l在该点处的时间延迟为

$\Delta t = {t_1} + {t_3}{\text{，}}$

(12)

${t_1} = \frac{{OD}}{{{v_1}}}{\text{，}}$

(13)

${t_3} = \frac{{O{X_{{\text{t}}}}}}{{{v_3}}}{\text{，}}$

(14)

式中，t₃为计算所需透射波的时间延迟。

在此，过A点做ON垂线，与ON延长线交于点 $C'$ 。由几何关系得

$OD = AO\sin \alpha {\text{，}}$

(15)

$OC' = AO\sin \gamma {\text{，}}$

(16)

${t_1} = \frac{{OD}}{{{v_1}}} = \frac{{OC'}}{{{v_3}}}{\text{，}}$

(17)

$\Delta t = {t_1} + {t_3} = \frac{{OC' + O{X_{{\text{t}}}}}}{{{v_3}}} = \frac{{C{X_{{\text{t}}}}}}{{{v_3}}}。$

(18)

设点 $V({x_3}, {y_3}) \in AC'$ ，式（18）可改写为

$\Delta t = \frac{1}{{{v_3}}}(({x_3} - {x_2})\cos \gamma + ({y_3} - {y_2})\sin \gamma )。$

(19)

至此，反射波n和透射波l的时间延迟计算公式中已经不含有其对应入射波m的任何参数，所有参数均为当前地震波的参数，且形式较为统一，便于进行迭代计算。

由于点 $U({x_1}, {y_1}) \in AC$ ，点 $V{{\text{(}}}{x_3}, {y_3}{{\text{)}}} \in AC'$ ，不妨将其均取为波前与地层界面相交的A点，作为所求地震波的等效零时刻点，由此便可以得到时间延迟的计算公式为

$\Delta {t_i} = \frac{1}{{{v_i}}}{{\text{(}}}\Delta {x_i}\cos {\theta _i} + \Delta {y_i}\sin {\theta _i}{{\text{)}}}{\text{，}}$

(20)

式中， $\Delta {x_i}$ ， $\Delta {y_i}$ 为所求点与A点坐标差， ${\theta _i}$ 为所求地震波传播方向与竖向的夹角， ${v_i}$ 为所求地震波在对应地层中的传播速度。

将其引申至三维模型，如图 3所示，波前与地层界面相交于A₁A₂，故A点可为A₁A₂上的任意一点，为便于计算可将其设为特殊点。此时时间延迟计算公式：

$\Delta {t_i} = \frac{1}{{{v_i}}}\left[ {\begin{array}{*{20}{c}} {\Delta {x_i}}&{\Delta {y_i}}&{\Delta {z_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos {\alpha _i}} \\ {\cos {\beta _i}} \\ {\cos {\gamma _i}} \end{array}} \right]{\text{，}}$

(21)

图 3 三维模型时间延迟计算示意图

Figure 3. Time delay calculation of 3D model

下载: 全尺寸图片幻灯片

式中， ${\alpha _i}$ ， ${\beta _i}$ ， ${\gamma _i}$ 分别为所求地震波传播方向与3个坐标轴的夹角。

3. 水平层状场地地震波斜入射

对于水平层状场地而言，地震波斜入射边界条件的计算重点在于等效节点力的求解上，而地震波在成层界面的反射透射导致其计算过程较为困难。因此，通过将上述时间延迟计算思路应用于水平层状场地的情况，采用迭代方式实现水平层状场地地震动的计算问题。

利用文献[]中地震波在成层界面上的对应关系，在计算过程中，对任意地震波，可通过将其产生过程中每一次转换所得到的振幅比值进行连乘来记录其与原入射波之间的幅度比值，即 ${E_i} = \prod {A{}_i{B_i}{C_i}{D_i}}$ 。其中A，B，C，D分别表示同类反射波、转换反射波、同类透射波、转换透射波与入射波幅值的比值。同时，由于地震波遇到成层界面时将会发生反射和透射现象，故在实际计算过程中存在终止条件的选择问题。在此取计算精度 $\delta$ ，忽略所得地震波与原入射波的比值小于 $\delta$ 的分量（即 ${E_i} < \delta$ ）及其后的地震波。考虑到地震是扩散衰减的，应用此方法进行计算无需担心收敛问题。 $\delta$ 的大小将在一定程度上影响计算精度和计算时效，故应在对二者进行权衡，但与动力有限元计算相比，此处计算量占比较小，无需过多斟酌。对地震波求解完成后，即可进一步计算黏弹性边界上各个节点的节点力，并施加在模型上。计算流程如图 4所示。

图 4 计算流程图

Figure 4. Flow chart of calculation

下载: 全尺寸图片幻灯片

以第k层地层边界面上的一点 $N({x_n}, {y_n}, {z_n})$ 为例，求解黏弹性边界水平层状场地地震波斜入射下的等效节点力。

假设Ⅰ₁为该层地层中与Z轴正向夹角大于90°的平面P波，Ⅰ₂为该层地层中与Z轴正向夹角小于90°的平面P波，Ⅰ₃为该层地层中与Z轴正向夹角大于90°的平面SV波，Ⅰ₄为该层地层中与Z轴正向夹角小于90°的平面SV波。计算过程中，应注意地震波传播过程中正方向的连续性，以保证公式的准确性。

假设 ${u_0}$ 为基岩入射波的位移时程，则N点处自由场位移：

${u_{Nx}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} {E_i}{u_0}(t - \Delta {t_i}){{\text{sin}}}{\alpha _i} + \mathop \sum \limits_{i \in {{{\text{I}}}_3}, {{{\text{I}}}_4}} {E_i}{u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i} {\text{，}}$

(22)

${u_{Ny}} = 0{\text{，}}$

(23)

${u_{Nz}} = - \mathop \sum \limits_{i \in {{{\text{I}}}_1}} {E_i}{u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i} + \mathop \sum \limits_{i \in {{{\text{I}}}_2}} {E_i}{u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i} +$

$\sum\limits_{i\in {{\text{I}}}_{3}}{E}_{i}{u}_{0}(t-\Delta {t}_{i}){\text{sin}}{\alpha }_{i}-\sum\limits_{i\in {{\text{I}}}_{4}}{E}_{i}{u}_{0}(t-\Delta {t}_{i}){\text{sin}}{\alpha }_{i}{\text{ }}。$

(24)

将式（22）~（24）对时间求导可以得到N点处自由场速度：

${\dot u_{Nx}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} {E_i}{\dot u_0}(t - \Delta {t_i}){{\text{sin}}}{\alpha _i} + \mathop \sum \limits_{i \in {{{\text{I}}}_3}, {{{\text{I}}}_4}} {E_i}{\dot u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i}{\text{，}}$

(25)

${\dot u_{Ny}} = 0{\text{，}}$

(26)

${\dot u_{Nz}} = - \mathop \sum \limits_{i \in {{{\text{I}}}_1}} {E_i}{\dot u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i} + \mathop \sum \limits_{i \in {{{\text{I}}}_2}} {E_i}{\dot u_0}(t - \Delta {t_i}){{\text{cos}}}{\alpha _i} +$

$\sum \limits_{i\in {{\text{I}}}_{3}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{sin}}{\alpha }_{i}-\sum \limits_{i\in {{\text{I}}}_{4}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{sin}}{\alpha }_{i}。$

(27)

在计算自由场应力时，可以根据波的叠加效应分别对模型的每个侧面及底面进行分析。

当对模型任意一点 $(x, y, z)$ ，都有x_n≥x时，N点处自由场应力：

${\sigma _{Nx}} = - \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} \frac{{{\lambda _k} + 2{G_k}{{\text{si}}}{{{\text{n}}}^2}{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) -$

$\sum \limits_{i\in {{\text{I}}}_{1}, {{\text{I}}}_{2}}\frac{{G}_{k}{\text{sin}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}{\text{，}}$

(28)

${\sigma _{Ny}} = 0{\text{，}}$

(29)

${\sigma _{Nz}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) - \mathop \sum \limits_{i \in {{{\text{I}}}_2}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}} \cdot {E_i}{\dot u_0}(t - \Delta {t_i}) + \mathop \sum \limits_{i \in {{{\text{I}}}_3}} \frac{{{G_k}{{\text{cos}}}2{\alpha _i}}}{{{c_{sk}}}}{E_i}{\ddot u_0}(t - \Delta {t_i}) -$

$\sum \limits_{i\in {{\text{I}}}_{4}}\frac{{G}_{k}{\text{cos}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}。$

(30)

当对模型任意一点 $(x, y, z)$ ，都有x_n≤x时，N点处自由场应力：

${\sigma _{Nx}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} \frac{{{\lambda _k} + 2{G_k}{{\text{si}}}{{{\text{n}}}^2}{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) +$

$\sum \limits_{i\in {{\text{I}}}_{1}, {{\text{I}}}_{2}}\frac{{G}_{k}{\text{sin}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}{\text{，}}$

(31)

${\sigma _{Ny}} = 0{\text{，}}$

(32)

${\sigma _{Nz}} = - \mathop \sum \limits_{i \in {{{\text{I}}}_1}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) + \mathop \sum \limits_{i \in {{{\text{I}}}_2}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}} \cdot$

${E_i}{\dot u_0}(t - \Delta {t_i}) - \mathop \sum \limits_{i \in {{{\text{I}}}_3}} \frac{{{G_k}{{\text{cos}}}2{\alpha _i}}}{{{c_{sk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) +$

$\sum \limits_{i\in {{\text{I}}}_{4}}\frac{{G}_{k}{\text{cos}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}。$

(33)

当对模型任意一点 $(x, y, z)$ ，都有y_n≥y时，N点处自由场应力：

${\sigma _{Nx}} = 0{\text{，}}$

(34)

${\sigma _{Ny}} = - \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} \frac{{{\lambda _k}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}){\text{，}}$

(35)

${\sigma _{Nz}} = 0。$

(36)

当对模型任意一点 $(x, y, z)$ ，都有y_n≤y时，N点处自由场应力：

${\sigma _{Nx}} = 0{\text{，}}$

(37)

${\sigma _{Ny}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} \frac{{{\lambda _k}}}{{{c_{pk}}}}{E_i}{\dot u_0}{{\text{(}}}t - \Delta {t_i}{{\text{)}}}{\text{，}}$

(38)

${\sigma _{Nz}} = 0。$

(39)

当对模型任意一点 $(x, y, z)$ ，都有z_n≤z时，N点处自由场应力：

${\sigma _{Nx}} = - \mathop {\sum }\limits_{i \in {{{\text{I}}}_1}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) + \mathop \sum \limits_{i \in {{{\text{I}}}_2}} \frac{{{G_k}{{\text{sin}}}2{\alpha _i}}}{{{c_{pk}}}} \cdot$

${E_i}{\dot u_0}(t - \Delta {t_i}) - \mathop \sum \limits_{i \in {{{\text{I}}}_3}} \frac{{{G_k}{{\text{cos}}}2{\alpha _i}}}{{{c_{sk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) + {{\text{ }}}$

$\sum \limits_{i\in {{\text{I}}}_{4}}\frac{{G}_{k}{\text{cos}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}{\text{，}}$

(40)

${\sigma _{Ny}} = 0{\text{，}}$

(41)

${\sigma _{Nz}} = \mathop \sum \limits_{i \in {{{\text{I}}}_1}, {{{\text{I}}}_2}} \frac{{{\lambda _k} + 2{G_k}{{\text{co}}}{{{\text{s}}}^2}{\alpha _i}}}{{{c_{pk}}}}{E_i}{\dot u_0}(t - \Delta {t_i}) -$

$\sum \limits_{i\in {{\text{I}}}_{3}, {{\text{I}}}_{4}}\frac{{G}_{k}{\text{sin}}2{\alpha }_{i}}{{c}_{sk}}{E}_{i}{\dot{u}}_{0}(t-\Delta {t}_{i}){\text{ }}。$

(42)

式（22）~（42）中的 $\Delta {t_i}$ 可根据式（21）进行计算。将式（22）~（42）代入等效节点力的计算公式，即可求得黏弹性边界下分层地层地震波斜入射过程中施加在模型边界节点上的等效节点力。在此基础上，依托现有的黏弹性边界地震动输入理论，利用常用的有限元计算软件和计算方法，即可对水平层状场地进行地震动有限元分析。

4. 方法验证

为验证本文方法的有效性和准确性，使用相关有限元计算软件建立如所示的计算模型。考虑两层水平层状场地结构进行模拟，模型上层 ${H_1} = 100\;{{\text{ m}}}$ ，下层 ${H_2} = 150\;{{\text{ m}}}$ 。长度L和宽度D均取250 m。采用T=0.25 s的狄拉克脉冲函数作为入射波，其位移的时程关系曲线如所示，先考虑模型两层密度 ${\rho _1} = {\rho _2} =$ 2.0×10³ ${{\text{kg/}}}{{{\text{m}}}^{{\text{3}}}}$ ，弹性模量 ${E_1} = {E_2} = 1{{\text{ GPa}}}$ ，泊松比 ${\mu _1} =$ ${\mu _2} = 0.2$ 的特殊情况，取入射角 $\alpha$ 为10°，20°平面P波和 $\alpha$ 为5°，15°平面SV波作为入射波，输入的时间间隔取0.005 s，对比本文方法与文献[19]方法的异同。

图 5 计算模型示意图

Figure 5. Computational model

下载: 全尺寸图片幻灯片

图 6 输入位移时程曲线

Figure 6. Time histories of input displacement

下载: 全尺寸图片幻灯片

分别采用两种方法得到有限元计算输入文件，图 7给出了文件的MD5计算结果，两种方法得到的MD5完全相同，说明文件一致，因此本文方法在有限次终止时与文献[19]方法的结果保持一致。

图 7 MD5计算结果

Figure 7. Results of MD5

下载: 全尺寸图片幻灯片

考虑一般情况，取上层密度 ${\rho _1} =$ 2.0×10³ ${{\text{kg/}}}{{{\text{m}}}^{{\text{3}}}}$ ，弹性模量 ${E_1} = 1\;{{\text{ GPa}}}$ ，泊松比 ${\mu _1} = 0.2$ ；下层密度 ${\rho _2} =$ 2.0×10³ ${{\text{kg/}}}{{{\text{m}}}^{{\text{3}}}}$ ，泊松比 ${\mu _2} = 0.2$ ，弹性模量 ${E_2} = 1.5$ GPa，分别取入射角 $\alpha$ 为10°，20°平面P波和 $\alpha$ 为5°，15°平面SV波作为入射波。

取 $\delta = 0.01$ ，图 8给出了模型在不同角度及不同类型入射波下地表中点处的水平（X方向）和竖向（Z方向）的位移时程，并与理论解进行对比。理论解采用频域传递矩阵法^[24]结合快速傅里叶变换得到^[14]。从图中可以看出，模型解与理论解吻合较好，相对误差较小，说明此方法在求解水平层状场地地震波斜入射方面具有良好的精度和准确性。

图 8 地表中点处的位移时程

Figure 8. Time histories of displacement at midpoint of surface

下载: 全尺寸图片幻灯片

以入射角 $\alpha =$ 20°的平面P波为例给出了模型在不同时刻下的位移场云图，从图中可以看出，本文的方法所引起地层内部变形规律与边界一致，变化较为均匀合理，且云图中的位移角度与地震波的斜入射的理论值基本上是吻合的，验证了本方法在计算水平层状场地地震波斜入射过程中的可靠性。

图 9 不同时刻位移场云图

Figure 9. Cloud diagram of displacement field at different time

下载: 全尺寸图片幻灯片

5. 结论

（1）提出了新的时间延迟计算方法。使用解析方式推导时间延迟的计算公式，使公式更加简化。在进行时间延迟的求解过程中，可通过所求目标地震波的具体信息及其所在地层的具体条件直接求解，摒弃以往计算方式中对所求地震波的全部来源依次进行计算的方式，使时间延迟公式统一，增强时间延迟计算方法的适用性。

（2）推导了基于等效零时刻点迭代的水平层状场地地震波斜入射方法。对于黏弹性边界下水平层状场地地震波波动输入问题，推导了基于本文时间延迟计算方法的人工边界上等效节点力的计算公式，进而得到了黏弹性边界下水平层状场地地震动输入的计算公式，拓展了黏弹性边界理论在水平层状场地斜入射地震动上的应用。

（3）进行方法验证。使用该方法建立水平层状场地三维有限元模型，求解两层地层结构地震波斜入射下的响应。将结果与理论值进行对比，验证了该方法的准确性和实用性。

图 1 反射波时间延迟计算示意图

Figure 1. Time delay calculation of reflected waves

下载: 全尺寸图片幻灯片

图 2 透射波时间延迟计算示意图

Figure 2. Time delay calculation of transmission waves

下载: 全尺寸图片幻灯片

图 3 三维模型时间延迟计算示意图

Figure 3. Time delay calculation of 3D model

下载: 全尺寸图片幻灯片

图 4 计算流程图

Figure 4. Flow chart of calculation

下载: 全尺寸图片幻灯片

图 5 计算模型示意图

Figure 5. Computational model

下载: 全尺寸图片幻灯片

图 6 输入位移时程曲线

Figure 6. Time histories of input displacement

下载: 全尺寸图片幻灯片

图 7 MD5计算结果

Figure 7. Results of MD5

下载: 全尺寸图片幻灯片

图 8 地表中点处的位移时程

Figure 8. Time histories of displacement at midpoint of surface

下载: 全尺寸图片幻灯片

图 9 不同时刻位移场云图

Figure 9. Cloud diagram of displacement field at different time

下载: 全尺寸图片幻灯片

参考文献(24)

[1]	胡聿贤. 地震工程学[M]. 北京: 地震出版社, 2006. HU Yu-xian. Earthquake Engineering[M]. Beijing: Seismological Press, 2006. (in Chinese)
[2]	耿萍. 铁路隧道抗震计算方法研究[D]. 成都: 西南交通大学, 2012. GENG Ping. Research on the Seismic Calculation Method of Railway Tunnel[D]. Chengdu: Southwest Jiaotong University, 2012. (in Chinese)
[3]	刘晶波, 杜义欣, 闫秋实. 黏弹性人工边界及地震动输入在通用有限元软件中的实现[J]. 防灾减灾工程学报, 2007, 27(增刊1): 37–42. https://cpfd.cnki.com.cn/Article/CPFDTOTAL-OGTY200703001009.htm LIU Jing-bo, DU Yi-xin, YAN Qiu-shi. Implementation of viscoelastic artificial boundary and ground motion input in general finite element software[J]. Journal of Disaster Prevention and Mitigation Engineering, 2007, 27(S1): 37–42. (in Chinese) https://cpfd.cnki.com.cn/Article/CPFDTOTAL-OGTY200703001009.htm
[4]	刘晶波, 谷音, 杜义欣. 一致黏弹性人工边界及黏弹性边界单元[J]. 岩土工程学报, 2006, 28(9): 1070–1075. doi: 10.3321/j.issn:1000-4548.2006.09.004 LIU Jing-bo, GU Yin, DU Yi-xin. Consistent viscous-spring artificial boundaries and viscous-spring boundary elements[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(9): 1070–1075. (in Chinese) doi: 10.3321/j.issn:1000-4548.2006.09.004
[5]	何建涛, 马怀发, 张伯艳, 等. 黏弹性人工边界地震动输入方法及实现[J]. 水利学报, 2010, 41(8): 960–969. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201008014.htm HE Jian-tao, MA Huai-fa, ZHANG Bo-yan, et al. Method and realization of seismic motion input of viscous-spring boundary[J]. Journal of Hydraulic Engineering, 2010, 41(8): 960–969. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201008014.htm
[6]	XING J, MECHANICS L Z O E, SSB, et al. Statistical research on S-wave incident angle[J]. Earthquake Research in China, 1994, 8(1): 121–131.
[7]	TAKAHIRO S, TAKESHI U, RYOICHI T, et al. Estimation of earthquake motion incident angle at rock site[C]//Proc of 12th World Conference Earthquake Engineering, NZ National Society for Earthquake Engineering. Auckland. 2000.
[8]	杜修力, 黄景琦, 赵密, 等. SV波斜入射对岩体隧道洞身段地震响应影响研究[J]. 岩土工程学报, 2014, 36(8): 1400–1406. doi: 10.11779/CJGE201408004 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) doi: 10.11779/CJGE201408004
[9]	王文晖. 地下结构实用抗震分析方法及性能指标研究[D]. 北京: 清华大学, 2013. WANG Wen-hui. Research on Practical Seismic Analysis Methods and Performance Index of Underground Structures[D]. Beijing: Tsinghua University, 2013. (in Chinese)
[10]	刘晶波, 王东洋, 谭辉, 等. 整体式反应位移法的理论推导及一致性证明[J]. 土木工程学报, 2019, 52(8): 18–23. https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201908003.htm LIU Jing-bo, WANG Dong-yang, TAN Hui, et al. Theorectical derivation and consistency proof of the integral response deformation method[J]. China Civil Engineering Journal, 2019, 52(8): 18–23. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201908003.htm
[11]	赵密, 欧阳文龙, 黄景琦, 等. 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
[12]	刘晶波, 王艳. 成层介质中平面内自由波场的一维化时域算法[J]. 工程力学, 2007, 24(7): 16–22. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX200707006.htm LIU Jing-bo, WANG Yan. A 1d time-domain method for in-plane wave motion of free field in layered media[J]. Engineering Mechanics, 2007, 24(7): 16–22. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX200707006.htm
[13]	刘晶波, 王艳. 成层半空间出平面自由波场的一维化时域算法[J]. 力学学报, 2006, 38(2): 219–225. https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB200602010.htm LIU Jing-bo, WANG Yan. A 1-d time-domain method for 2-d wave motion in elastic layered half-space by antiplane wave oblique incidence[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(2): 219–225. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB200602010.htm
[14]	赵密, 杜修力, 刘晶波, 等. P-SV波斜入射时成层半空间自由场的时域算法[J]. 地震工程学报, 2013, 35(1): 84–90. https://www.cnki.com.cn/Article/CJFDTOTAL-ZBDZ201301012.htm ZHAO Mi, DU Xiu-li, LIU Jing-bo, et al. Time-domain method for free field in layered half space under P-SV waves of oblique incidence[J]. China Earthquake Engineering Journal, 2013, 35(1): 84–90. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZBDZ201301012.htm
[15]	ZHAO M, YIN H Q, DU X L, et al. 1D finite element artificial boundary method for layered half space site response from obliquely incident earthquake[J]. Earthquakes and Structures, 2015, 9(1): 173–194.
[16]	尹侯权. 地震波斜入射时成层半空间场地反应分析方法及其应用[D]. 北京: 北京工业大学, 2015. YIN Hou-quan. A Method for Layered Half-Space Site Response Under Obliquely Incident Earthquake and its Application[D]. Beijing: Beijing University of Technology, 2015. (in Chinese)
[17]	王笃国, 赵成刚. 地震波斜入射下考虑场地非线性、地形效应和土结动力相互作用的大跨连续刚构桥地震响应分析[J]. 工程力学, 2017, 34(4): 32–41. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201704005.htm WANG Du-guo, ZHAO Cheng-gang. Seismic analysis of long-span continuous rigid frame bridge considering site nonlinearity, topography effect and soil-structure dynamic interaction under oblique incidence[J]. Engineering Mechanics, 2017, 34(4): 32–41. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201704005.htm
[18]	王笃国, 赵成刚. 地震波斜入射时二维成层介质自由场求解的等效线性化方法[J]. 岩土工程学报, 2016, 38(3): 554–561. doi: 10.11779/CJGE201603020 WANG Du-guo, ZHAO Cheng-gang. Two- dimensional equivalent linear seismic analysis of free field in layered half-space due to oblique incidence[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(3): 554–561. (in Chinese) doi: 10.11779/CJGE201603020
[19]	黄景琦. 岩体隧道非线性地震响应分析[D]. 北京: 北京工业大学, 2015. HUANG Jing-qi. Study on Nonlinear Seismic Response of Rock Tunnels[D]. Beijing: Beijing University of Technology, 2015. (in Chinese)
[20]	徐海滨, 杜修力, 赵密, 等. 地震波斜入射对高拱坝地震反应的影响[J]. 水力发电学报, 2011, 30(6): 159–165. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201106031.htm XU Hai-bin, DU Xiu-li, ZHAO Mi, et al. Effect of oblique incidence of seismic waves on seismic responses of high arch dam[J]. Journal of Hydroelectric Engineering, 2011, 30(6): 159–165. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201106031.htm
[21]	刘云贺, 张伯艳, 陈厚群. 拱坝地震输入模型中黏弹性边界与黏性边界的比较[J]. 水利学报, 2006, 37(6): 758–763. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200606019.htm LIU Yun-he, ZHANG Bo-yan, CHEN Hou-qun. Comparison of spring-viscous boundary with viscous boundary for arch dam seismic input model[J]. Journal of Hydraulic Engineering, 2006, 37(6): 758–763. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200606019.htm
[22]	李明超, 张佳文, 张梦溪, 等. 地震波斜入射下混凝土重力坝的塑性损伤响应分析[J]. 水利学报, 2019, 50(11): 1326–1338, 1349. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201911005.htm LI Ming-chao, ZHANG Jia-wen, ZHANG Meng-xi, et al. Plastic damage response analysis of concrete gravity dam due to obliquely incident seismic waves[J]. Journal of Hydraulic Engineering, 2019, 50(11): 1326–1338, 1349. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201911005.htm
[23]	WEI X C, CHEN T S. Joint PP and PS AVO inversion based on Zoeppritz equations[J]. Earthquake Science, 2011, 24(4): 329–334.
[24]	傅淑芳, 刘宝诚. 地震学教程[M]. 北京: 地震出版社, 1991. FU Shu-fang, LIU Bao-cheng. Seismology Tutoria[M]. Beijing: Scismic Press, 1991. (in Chinese)

施引文献(2)

期刊类型引用(2)

1.	杜闯，高启辉，宋帅. 地震波斜入射下港珠澳大桥沉管隧道地震响应分析. 河北工业大学学报. 2024(03): 76-86 . 百度学术
2.	邓捷程，吕龙，孙海忠，吕玺琳. 基于等效线性土体模型的地基与结构整体地震响应模拟. 结构工程师. 2023(03): 85-91 . 百度学术