疏干条件下隧洞涌水量及孔压分布解析研究

吴建; 周志芳; 庄超

doi:10.11779/CJGE202010015

疏干条件下隧洞涌水量及孔压分布解析研究 English Version

河海大学地球科学与工程学院,江苏　南京　211100

基金项目:

国家重大研究计划重点支持项目 91747204

详细信息

作者简介:
吴建(1994—),男,博士研究生,主要从事地下水数值模拟相关研究。E-mail：wujian1994@hhu.edu.cn

通讯作者:
周志芳, E-mail：zhouzf@hhu.edu.cn

中图分类号: U45
计量
- 文章访问数: 289
- HTML全文浏览量: 46
- PDF下载量: 152
出版历程
- 收稿日期: 2019-01-06
- 网络出版日期: 2022-12-07
- 刊出日期: 2020-09-30

Analytical study on water inflow and pore water pressure of tunnels under drainage

School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China

摘要

摘要: 隧洞涌水和孔隙水压力计算一直是地下工程领域最为关注的焦点问题,然而对于疏干条件下隧洞涌水量及衬孔隙水压力的变化规律仍未得到很好的解决。针对浅埋、深埋及衬砌隧洞分别进行解析分析,然后利用数值计算方法反算出关键参数（影响距离、等效地下水位及衬砌等效地下水位）的计算公式,进一步提出了隧洞涌水量的解析公式,探讨了孔隙水压力的分布规律。结果表明,疏干条件下隧洞涌水量取决于隧洞半径、初始地下水位高度及围岩与衬砌相对渗透性,其中围岩与衬砌相对渗透性起到决定性因素；衬砌圈显著改变了隧洞涌水量大小与外水压力分布规律,随着相对渗透性的增大,孔压解析解与数值解的拟合程度逐渐提高；通过与传统解析公式及数值的相互对比,提出的解析公式真实刻画了浅埋、深埋和衬砌隧洞渗流场特征,计算结果更加合理和准确。
- 隧洞 /
- 涌水量 /
- 孔隙水压力 /
- 解析公式 /
- 数值模拟
Abstract: The calculation of water inflow and pore water pressure of tunnels is always the focus of underground engineering. However, the effect of water inflow and pore water pressure of tunnels under drainage has not been well solved. The analytical analysis of shallow, deep and lined tunnels is carried out respectively. Then, the expressions for the key parameters (influencing distance, equivalent underground water level and equivalent underground water level of linings) are obtained by using the numerical method, and the analytical formulae for water inflow into tunnels are put forward. Finally, the distribution rules of pore water pressure are discussed. The results show that the water inflow into tunnels depends on the tunnel radius, the initial groundwater level and the relative permeability of the surrounding rock and linings, among which the relative permeability of the surrounding rock and linings plays a key role. The lining ring significantly changes the amount of water inflow and distribution of pore water pressure of tunnels. With the increase of the relative permeability, the fitting degree between the analytical solution and the numerical solution of pore water pressure increases gradually. Compared with traditional analytical solution and numerical solution, the proposed analytical formula more truly describes the seepage field characteristics of shallow, deep and lined tunnels, and the calculated results are more reasonable and accurate.
- tunnel /
- water inflow /
- pore water pressure /
- analytical formula /
- numerical simulation

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. Schematic representation of shallow tunnel

下载: 全尺寸图片幻灯片

图 2 深埋隧洞渗流场示意图

Figure 2. Schematic representation of deep tunnel

下载: 全尺寸图片幻灯片

图 3 衬砌隧洞渗流场示意图

Figure 3. Schematic representation of lined tunnel

下载: 全尺寸图片幻灯片

图 4 有限元模型网格剖分示意图

Figure 4. Finite element mesh of models

下载: 全尺寸图片幻灯片

图 5 隧洞涌水量数值计算结果

Figure 5. Water inflows of tunnels in different cases

下载: 全尺寸图片幻灯片

图 6 不同初始地下水位高度下疏干隧洞水平向与垂向的孔隙水压力分布

Figure 6. Distribution of pore water pressure along horizontal and vertical directions of drained tunnel under various initial groundwater levels

下载: 全尺寸图片幻灯片

图 7 不同围岩与衬砌相对渗透性下衬砌隧洞水平向与垂向的.孔隙水压力分布

Figure 7. Distribution of pore water pressure along horizontal and vertical directions of lined tunnel with different k_m/k_l

下载: 全尺寸图片幻灯片

图 8 d/D关于r/h的散点分布及拟合曲线

Figure 8. Scatter distribution and fitting curves of d/D with respect to r/h

下载: 全尺寸图片幻灯片

图 9 $h^{'} / h$ 关于 $r / h$ 的散点分布及拟合曲线

Figure 9. Scatter distribution and fitting curves of $h^{'} / h$ with respect to $r / h$

下载: 全尺寸图片幻灯片

图 10 $h^{'} / h$ 关于 $r / h$ 的散点分布及拟合曲线

Figure 10. Scatter distribution and fitting curves of $h^{'} / h$ with respect to $r / h$

下载: 全尺寸图片幻灯片

图 11 不同计算方法下隧洞中心右侧水平向与顶拱垂向的水压力分布曲线

Figure 11. Distribution of pore water pressure along horizontal and vertical directions of tunnel by different methods

下载: 全尺寸图片幻灯片

图 12 折减系数与相对渗透性之间的关系曲线

Figure 12. Relationship between discount coefficient and k_m/k_l

下载: 全尺寸图片幻灯片

图 13 不同计算方法隧洞涌水量结果对比

Figure 13. Comparison of water inflows into tunnel by different methods

下载: 全尺寸图片幻灯片

参考文献(17)

[1]	刘志春, 王梦恕. 隧道工程因素对地下水环境影响研究[J]. 岩土力学, 2015, 36(增刊2): 281-288. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2038.htm LIU Zhi-chun, WANG Meng-shu. Research on impact of tunnel engineering factors on groundwater environment[J]. Rock and Soil Mechanics, 2015, 36(S2): 281-288. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2038.htm
[2]	潘东东, 李术才, 许振浩, 等. 岩溶隧道承压隐伏溶洞突水模型试验与数值分析[J]. 岩土工程学报, 2018, 40(5): 828-836. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805009.htm PAN Dong-dong, LI Shu-cai, XU Zhen-hao, et al. Model tests and numerical analysis for water inrush caused by karst caves filled with confined water in tunnels Chinese[J]. Journal of Geotechnical Engineering, 2018, 40(5): 828-836. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805009.htm
[3]	HARR M E, Groundwater and Seepage[M]. New York: Dover Publications, 1962: 249-264.
[4]	吴建, 周志芳, 李鸣威, 等. 隧洞涌水量预测计算方法研究进展[J]. 工程地质学报, 2019, 27(4): 890-902. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201904023.htm WU Jian, ZHOU Zhi-fang, LI Ming-wei, et al. Advance on the methods for predicting water inflow into tunnels[J]. Journal of Engineering Geology, 2019, 27(4): 890-902. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201904023.htm
[5]	GOODMAN R E, MOYE D G, VAN SCHALKWYK A, et al. Ground Water Inflows During Tunnel Driving[M]. California: College of Engineering, University of California, 1964.
[6]	LEI S. An analytical solution for steady flow into a Ttonnel[J]. Groundwater, 1999, 37(1): 23-26. doi: 10.1111/j.1745-6584.1999.tb00953.x
[7]	RAYMER J H. Predicting groundwater inflow into hard-rock tunnels: estimating the high-end of the permeability distribution[C]//2001 Rapid Excavation and Tunneling Conference, 2001, Littleton: 1027-1038.
[8]	EL TANI M. Circular tunnel in a semi-infinite aquifer[J]. Tunnelling and Underground Space Technology, 2003, 18(1): 49-55. doi: 10.1016/S0886-7798(02)00102-5
[9]	MOON J, FERNANDEZ G. Effect of excavation-induced groundwater level drawdown on tunnel inflow in a jointed rock mass[J]. Engineering Geology, 2010, 110(3/4): 33-42.
[10]	苏凯, 陈锐, 周亚峰, 等. 圆形隧洞渗流量计算的解析方法及其应用[J]. 岩石力学与工程学报, 2017, 36(增刊1): 3332-3341. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2017S1024.htm SU Kai, CHEN Rui, ZHOU Ya-feng, et al. Analytical method for calculating the groundwater inflow of circular tunnel and its application[J]. Rock and Soil Mechanics, 2017, 36(S1): 3332-3341. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2017S1024.htm
[11]	EL TANI M, KAMALI A, GHOLAMI M A. Analytic assessment of the water table drawdown, seepage, and back pressure at rudbar PSPP[J]. Rock Mechanics and Rock Engineering, 2019, 52(7): 2227-2243.
[12]	HOOGHOUDT S B. Bijdragen Tot De Kennis Van Eenige Natuurkundige Grootheden Van Den Grond: Algemeene Beschouwing Van Het Probleem Van De Detailontwatering En De Infiltratie Door Middel Van Parallel Loopende Drains, Greppels, Slooten En Kanalen[R]. Hague: Algemeene Landsdrukkerij, 1940: 515-707.
[13]	FERNANDEZ G, MOON J. Excavation-induced hydraulic conductivity reduction around a tunnel-part 1: Guideline for estimate of ground water inflow rate[J]. Tunnelling and Underground Space Technology, 2010, 25(5): 560-566.
[14]	FERNANDEZ G, MOON J. Excavation-induced hydraulic conductivity reduction around a tunnel-part 2: Verification of proposed method using numerical modeling[J]. Tunnelling and Underground Space Technology, 2010, 25(5): 567-574.
[15]	周亚峰, 苏凯, 伍鹤皋. 水工隧洞钢筋混凝土衬砌外水压力取值方法研究[J]. 岩土力学, 2014, 35(增刊2): 198-203. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2014S2027.htm ZHOU Ya-feng, SU Kai, WU He-gao. Study of external water pressure estimation method for reinforced concrete lining of hydraulic tunnels[J]. Rock and Soil Mechanics, 2014, 35(S2): 198-203. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2014S2027.htm
[16]	DAI Y, ZHOU Z. Steady seepage simulation of underground oil storage caverns based on Signorini type variational inequality formulation[J]. Geosciences Journal, 2015, 19(2): 341-355.
[17]	郑宏, 刘德富, 李焯, 等. 一个新的有自由面渗流问题的变分不等式提法[J]. 应用数学和力学, 2005, 26(3): 363-371. https://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200503016.htm ZHENG Hong, LIU De-fu, LEE C F, et al. New variational inequality formulation for seepage problems with free surfaces[J]. Applied Mathematics and Mechanics, 2005, 26(3): 363-371. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200503016.htm

施引文献(2)

期刊类型引用(2)

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