深厚复杂海峡场地二维非线性地震反应特性

王彦臻; 范宏飞; 赵凯; 陈国兴; 江志杰

doi:10.11779/CJGE20221307

深厚复杂海峡场地二维非线性地震反应特性 English Version

王彦臻^{1, 2,},
范宏飞¹,
赵凯^{1, 2},
陈国兴^{1, 2, ,},
江志杰^{1, 3}

1.
南京工业大学岩土工程研究所, 江苏南京 210009
2.
江苏省土木工程防震技术研究中心, 江苏南京 210009
3.
中国地震局兰州地震研究所, 甘肃兰州 730000

基金项目:

国家自然科学基金项目 51978334

详细信息

作者简介:
王彦臻(1990—)，男，博士研究生，主要从事岩土地震工程方面的研究工作。E-mail: soledadwang@163.com

通讯作者:
陈国兴, E-mail: gxc6307@163.com

中图分类号: TU435
计量
- 文章访问数: 326
- HTML全文浏览量: 35
- PDF下载量: 107
出版历程
- 收稿日期: 2022-10-22
- 网络出版日期: 2023-06-05
- 刊出日期: 2024-01-31

2D nonlinear seismic response characteristics of a strait site with deep inhomogeneous soil deposits

WANG Yanzhen^{1, 2,},
FAN Hongfei¹,
ZHAO Kai^{1, 2},
CHEN Guoxing^{1, 2, ,},
JIANG Zhijie^{1, 3}

1.
Institute of Geotechnical Engineering, Nanjing Tech University, Nanjing 210009, China
2.
Civil Engineering and Earthquake Disaster Prevention Center of Jiangsu Province, Nanjing 210009, China
3.
Lanzhou Institute of Seismology, CEA, Lanzhou 730000, China

摘要

摘要: 海峡地层常为深厚沉积盆地构造，海床微地形类型多样，海峡场地地震反应受微地形的影响显著。以跨琼州海峡通道剖面场地为研究对象，综合考虑海峡盆地构造、海床微地形和海床沉积土异质性，建立海床剖面二维精细化有限元模型；依据场址区域地震环境与历史地震活动性，选择Kik-net台网井下台阵基岩地震记录作为海床基岩输入地震动，考虑土的非线性滞回特性，模拟海峡场地的非线性地震反应。微地形效应和海床沉积土异质性对跨海峡场地地震反应的耦合影响显著，表现为盆地类共振现象、海床凹陷边缘和海床表面凸起处的反应加剧；与海峡场地基本频率接近的基岩地震动频率成分易于传播与放大；与海床水平向地震反应相比，海床地表微地形对其竖向地震反应的影响更为显著；水平向和竖向地震波传播的叠加效应会在海床表面产生面波，导致海床地表水平向地震反应放大。
- 海峡场地 /
- 琼州海峡 /
- 盆地构造 /
- 海床微地形 /
- 沉积土异质性 /
- 土的非线性
Abstract: The strait seabed is usually a deep sediment basin with various types of micro-topographies, and its seismic responses are strongly affected by the micro-topographies. A nonlinear seismic response analysis of 2D refined finite element simulated model for a strait section crossing Qiongzhou Strait is established. In the response analyses of the cross-strait site, the basin structure, surface micro-topography and spatial inhomogeneous soil with variable S- and P- wave velocities are considered. Considering the regional tectonic setting and historical seismicity around the strait site, three bedrock recordings of moderate-far field earthquake scenarios are justified for determining the bedrock shakings. The combined influences of seabed cross-section topography and spatial inhomogeneity of soils on the seismic site responses are complicated and significant, including the basin resonance-like effects and the aggravating effect of the local concave seabed surface edges and seabed convex areas on soil amplification. The simulated results imply that it is easier to propagate upward earthquake waves through the soil deposits for the bedrock motion components close to the fundamental frequency of strait seabed. Compared with the horizontal seabed site responses, the influences of the surface micro-topography on the vertical seismic responses are more significant. The generated surface waves at seabed surface due to the wave constructive interference of the horizontal and vertical components upward earthquake waves aggravate the horizontal responses of seabed site surface.
- strait site /
- Qiongzhou strait /
- basin structure /
- surface micro-topography /
- inhomogeneous soil deposit /
- soil nonlinearity

HTML全文

0. 引言

DDA(discontinuous deformation analysis)即非连续变形分析，是一种基于非连续介质力学的数值分析方法。该方法以块体作为基本对象进行分析，块体边界可表示实际岩体裂隙，将各块体的位移视作基本未知量，通过块体间的接触关系将独立块体组合成一个块体系统，在处理块体大变形、大位移问题上具有一定的优势。因此DDA被广泛地应用在岩体这一块体离散介质的研究当中。众多学者将DDA应用于边坡稳定、岩体锚固、矿山开挖、地震等众多工程问题当中^[1-6]。

原始DDA程序主要考虑块体应力应变场，难以直接应用到岩体渗流问题模拟当中。考虑复杂渗流场-应力应变场耦合最终建立能够真实模拟岩体赋存条件的DDA模型将大幅提高该算法的模拟精度和实用性，成为国内外学者关注的热点^[7]。例如，张国新等^[8]基于能量原理将水压荷载加入到原始DDA程序中，分析了水位变化条件下岩质边坡的倾倒变形机理。虞松等^[9]提出了一种矩阵搜索方法，识别并删除DDA生成的裂隙岩体中的不连通裂隙，以生成更精确的裂隙水通网络，并对地下油库水封条件开挖的工程案例进行模拟，发现考虑流固耦合所获得的裂隙水压相对较小而涌水量则有所增加。王知深^[10]在DDARF基础上推导达西-非达西流裂隙岩体渗流场和应力场控制方程，并建立DDARF平台裂隙岩体渗流应力耦合模型，将其应用在碳酸盐岩水力压裂过程分析中。江巍等^[11]采用DDA对藕塘滑坡进行了失稳模式分析，采用块体在饱和状态下的强度参数来进行降雨工况的模拟，但在分析过程中没能考虑流固耦合的影响。Hu等^[12]利用离散裂隙网络模型模拟页岩的自然裂隙，结合DDA提出复杂水力裂隙网络耦合流体力学模型，并将其应用在页岩气压裂开采设计中。Morgan等^[13]利用有限体积裂隙网络模型模拟裂缝中的可压缩流体，研究了在内部有开口加压的格里菲斯裂缝在无渗漏无补给条件下压强和裂缝开度的变化，通过与半解析解进行对比验证了算法的正确性，并对高黏性压裂液在不渗透介质中的传播进行了数值分析。

上述研究成果对于DDA流固耦合的处理多是基于饱和稳态渗流理论。在实际工程中降雨入渗、库水位变化等工程案例均属于非饱和瞬态渗流问题，而目前运用DDA处理这一类问题通常是将降雨工况下材料强度参数进行折减，用饱和强度参数进行模拟，该方法并不能反映降雨入渗过程。针对这一问题，本文对DDA进行二次开发，推导DDA瞬态整体渗流矩阵，并将其加入到总体平衡方程中，实现DDA非饱和瞬态渗流流固耦合计算，通过土柱入渗、砂槽入渗试验验证该计算方法的正确性，为实际工程提供新的数值模拟手段和借鉴参考。

1. DDA稳态渗流计算

DDA以岩块和裂隙组成的岩体作为研究对象，由于岩块渗透能力较弱，岩体渗流通常只考虑裂隙渗流。采用立方定律来描述岩体单裂隙渗流：

。

$q = \frac{{g{e^3}}}{{12\mu }}J 。$

(1)

式中： $e$ 为岩体裂隙开度；q为岩体裂隙流量；μ为流体的动力黏滞系数；J为光滑平板两端的水力梯度。

在实际工程当中，岩体裂隙具有一定的起伏和粗糙度。裂隙的起伏和粗糙度会使流体在裂隙中的流线变长，因此在相同情况下单位时间内通过粗糙裂隙的渗流量也会小于光滑裂隙。为了使计算结果更接近于真实值，引入粗糙度修正系数C，可得

。

$q = \frac{{g{e^3}}}{{12C\mu }}J 。$

(2)

根据相关的水力学试验，粗糙度修正系数C可由岩体裂隙内表面凹凸粗糙程度Δ与裂隙开度e计算得到^[14]

。

$C = 1 + {\left( {\frac{\Delta }{e}} \right)^{1.5}} 。$

(3)

为将其写入DDA程序当中，将其改写成类似于DDA平衡方程的系数矩阵和未知数矩阵乘积的形式，即

，

$q = T\Delta H \text{，}$

(4)

式中，T为岩体裂隙的导流系数。

对于岩体裂隙网络中的某一节点，其总流量为各条单裂隙流量对该节点流量的影响之和。对一个由n个节点组成的岩体模型来说，其整体渗流方程可写为

$\left[ {\begin{array}{*{20}{c}} {{T_{11}}}&{{T_{12}}}& \ldots &{{T_{1n}}} \\ {{T_{21}}}&{{T_{22}}}& \ldots &{{T_{2n}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{T_{n1}}}&{{T_{n2}}}& \ldots &{{T_{nn}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{H_1}} \\ {{H_2}} \\ \vdots \\ {{H_n}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{Q_1}} \\ {{Q_2}} \\ \vdots \\ {{Q_n}} \end{array}} \right\} 。$

(5)

式中：裂隙导流系数矩阵T的对角线元素T_ii为与节点i相连接的所有裂隙的导流系数之和；T_ij为由节点i和节点j组成裂隙的导流系数的相反数；水头矩阵H中的H_i为节点i对应的水头；流量矩阵Q中的Q_i为节点i对应的总流量，由连接该点的所有裂隙流量叠加得到，且以流体流入该点为负流出该点为正。

在渗流的模拟过程当中，为了求解式（5）还需要在迭代前给出边界条件，即流量边界和水头边界。当某一节点的水头为已知量时，在迭代时将该点导流系数矩阵T的对角线元素T_ii设置为1，其余元素T_ij设置为0，以保证水头不变。当某一节点的流量为已知量时，将流量矩阵Q中对应的元素Q_i设为已知量即可^[15]。

2. DDA瞬态渗流计算

2.1 瞬态整体渗流矩阵推导

为了运用DDA进行瞬态渗流场模拟，将用体积含水量表示的二维非饱和Richards方程引入到DDA程序当中：

$\frac{\partial }{{\partial x}}\left( {{D_x}(\theta )\frac{{\partial \theta }}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{D_y}(\theta )\frac{{\partial \theta }}{{\partial y}}} \right) = \frac{{\partial \theta }}{{\partial t}} 。$

(6)

式中：θ为体积含水量，D_x(θ)，D_y(θ)分别为x方向和y方向上的水力扩散系数。

推导符合DDA形式的瞬态整体渗流矩阵，需要将二维Richards方程进行离散化处理。DDA中的块体边界代表了岩体中的裂隙，具有实际的物理意义，因此直接用DDA中块体划分形成的网格作为基础，将 ${D_x}(\theta )\frac{{\partial \theta }}{{\partial x}}$ 视为整体采用中心差分将式（6）中的第一项偏导数离散化，可得

$\frac{\partial }{\partial x}\left({D}_{x}(\theta )\frac{\partial \theta }{\partial x}\right)=\frac{D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i+\frac{1}{2}hx}^{j+1}}{\partial x}-D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i-\frac{1}{2}hx}^{j+1}}{\partial x}}{\Delta x}。$

(7)

式中：θ的上标j+1为时间差分，表示在第j+1时间步对应的体积含水量，下标 $i + \frac{1}{2}hx$ 为空间差分，表示对第i个节点在x方向上取1/2个空间步长对应的体积含水量；Δx为x方向上两个节点的间距。

对式（7）中的偏微分式继续离散化处理，可得

$\left.\begin{array}{l}D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i+\frac{1}{2}hx}^{j+1}}{\partial x}=D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{{\theta }_{i+hx}^{j+1}-{\theta }_{i}^{j+1}}{\Delta x}\text{ }\text{，}\\ D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i-\frac{1}{2}hx}^{j+1}}{\partial x}=D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{{\theta }_{i}^{j+1}-{\theta }_{i-hx}^{j+1}}{\Delta x}\text{ }。\end{array}\right\}$

(8)

对于y方向，也有类似于式（8）的计算结果，只需要将空间差分的下标和节点间距分别记为y方向上的对应值即可。

对式（6）的等式右侧同样进行离散化处理，有

$\frac{{\partial \theta }}{{\partial t}} = \frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} 。$

(9)

由于差分在渗流模拟过程中具有一定的滞后性，因此对于水力扩散系数D(θ)可以近似的用第j时步的计算结果去表达第j+1时步^[16]。为了简化算式表达进行如下处理：

$\left.\begin{array}{l}\alpha {x}_{i}^{j}=\frac{D\left({\theta }_{i+\frac{1}{2}hx}^{j}\right)}{\Delta {x}^{2}}\text{ }\text{，}\\ \beta {x}_{i}^{j}=\frac{D\left({\theta }_{i-\frac{1}{2}hx}^{j}\right)}{\Delta {x}^{2}}\text{ }\text{，}\\ \alpha {y}_{i}^{j}=\frac{D\left({\theta }_{i+\frac{1}{2}hy}^{j}\right)}{\Delta {y}^{2}}\text{ }\text{，}\\ \beta {y}_{i}^{j}=\frac{D\left({\theta }_{i-\frac{1}{2}hy}^{j}\right)}{\Delta {y}^{2}}\text{ }。\end{array}\right\}$

(10)

对式（7）~（10）进行整理可以得对某一节点体积含水量离散形式的表达式：

$\begin{aligned} \frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} =& \alpha x_i^j\theta _{i + hx}^{j + 1} + \beta x_i^j\theta _{i - hx}^{j + 1} + \alpha y_i^j\theta _{i + hy}^{j + 1} +\\ &\beta y_i^j\theta _{i - hy}^{j + 1} - (\alpha x_i^j + \beta x_i^j + \alpha y_i^j + \beta y_i^j)\theta _i^{j + 1} 。 \end{aligned}$

(11)

式中：Δt为瞬态渗流的时间步长，上标为j表示本时步的体积含水量，均为已知量；j+1表示下一时步的体积含水量，通过本时步的状态求解；下标i+hx表示在x方向上与第i个节点间隔hx的下一个节点对应的体积含水量。

由此可知式（11）中每节点在x方向和y方向上均有两个相邻节点连接，对应了采用较为规则的四边形网格划分的情况。在DDA中块体可以是任意形状，因此网格的划分也会随之变得不规则，对网格中某一点i来说，其体积含水量仍然受到与之相连的其他节点的影响，将式（11）整理成矩阵形式，可得

$\left[ {\begin{array}{*{20}{c}} {{A_{11}}}&{{A_{12}}}& \ldots &{{A_{1n}}} \\ {{A_{21}}}&{{A_{22}}}& \ldots &{{A_{2n}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{A_{n1}}}&{{A_{n2}}}& \ldots &{{A_{nn}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\theta _1^{j + 1}} \\ {\theta _2^{j + 1}} \\ \vdots \\ {\theta _n^{j + 1}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{B_1}} \\ {{B_2}} \\ \vdots \\ {{B_n}} \end{array}} \right\} 。$

(12)

与式（5）相似，式（12）的系数矩阵A是一个只在相邻节点连接的位置上有对应数值的稀疏矩阵。其中，非对角元素A_ij表示由i，j两点构成裂隙对应的系数，即式（10）中的4个代数式之一；对角元素A_ii表示在第i行中所有非对角元素A_ij的和再加上时间步长 $\Delta t$ 的倒数；而等式右侧的常数矩阵B为每一个节点本时步的体积含水量θ与时间步长 $\Delta t$ 的比值。由此可知，对于系数矩阵A由于有时间项的存在，其每一行的对角线元素都大于该行非对角线元素的总和，是具有严格对角占优特点的矩阵。为了保证其收敛性，采用Gauss-Seidel迭代法求解。

对式（12）进行求解需要已知初始条件下每个节点的体积含水量，以及计算模型的边界条件。对DDA模型先进行稳态渗流计算，以稳态渗流的结果作为瞬态渗流的初始条件。对体积含水量已知的边界，将对角线元素A_ii设为1，非对角线元素A_ij设为0即可。

得到节点体积含水量后，通过下式转化成有效饱和度：

${S_{\text{e}}} = \frac{{\theta - {\theta _{\text{r}}}}}{{{\theta _{\text{s}}} - {\theta _{\text{r}}}}} 。$

(13)

式中：S_e为有效饱和度；θ_s为饱和体积含水量；θ_r为残余体积含水量。

根据Van Genuchten^[17]模型和Mualem模型^[18]来表示有效饱和度和相对渗透系数、裂隙压力间的关系：

${P_{\text{c}}} = \frac{{{{({S_{\text{e}}}^{ - \frac{1}{m}} - 1)}^{\frac{1}{n}}}}}{\alpha } \text{，}$

(14)

${T_{\text{r}}} = {S_{\text{e}}}^{\frac{1}{2}}{\left[ {1 - {{\left( {1 - {S_{\text{e}}}^{\frac{1}{m}}} \right)}^m}} \right]^2} 。$

(15)

式中：m，n，α均为土水特征曲线经验拟合参数，m可由 $m = 1 - \frac{1}{n}$ 计算得到；P_c为压力水头；T_r为相对渗透系数，加入到导流系数矩阵T中。

2.2 流固耦合计算

将渗流场计算得到的水头矩阵块体外荷载加入到DDA总体平衡方程中。假设有一块体i，该块体由点1和点2组成的边界受到水压力的作用。将水压力视作在两点之间呈梯形分布的荷载，如图 1所示。

图 1 块体水压力

Figure 1. Water pressures on block

下载: 全尺寸图片幻灯片

点1上的水压力h分别向x，y方向分解，并规定与坐标轴同向为正，与坐标轴反向为负。

$\left.\begin{array}{l}{f}_{x}{}_{1}=-{h}_{1}\mathrm{sin}\alpha \text{ }\text{，}\\ {f}_{y}{}_{1}={h}_{1}\mathrm{cos}\alpha \text{ }。\end{array}\right\}$

(16)

式中， $\alpha$ 角为在点1处块体边界线与水平线的夹角。

对于点2上的水压力沿x，y方向分解可得到相似的结果。类似于DDA的体积力荷载子矩阵，点1和点2的水压力对块体i产生的势能为

${\varPi _i} = - \int\limits_0^{{l_{12}}} {({f_{{\text{xs}}}}{u_{\text{s}}} + {f_{{\text{ys}}}}{v_{\text{s}}})} {\text{d}}s 。$

(17)

式中：l₁₂为点1到点2的距离；s为在该裂隙长度方向上选取的某一微元段；u_s和v_s为该微元段在x，y方向上的平移分量。

对式（17）进行势能最小化处理：

$- \frac{{\partial {\Pi _i}(0)}}{{\partial {d_{ri}}}} = {\boldsymbol T_i}^{\text{T}}({x_1},{y_1})\left[ E \right] + {\boldsymbol T_i}^{\text{T}}({x_2},{y_2})\left[ G \right]$

(18)

式中： ${\boldsymbol T_i}^{\text{T}}$ (x₁, y₁)和T_i^T(x₂, y₂)分别为点1和点2的位移转换矩阵，[E]和[G]由下式表示：

$\begin{array}{l}\left[E\right]=\left(\begin{array}{cc}-\frac{{l}_{12}}{3}\mathrm{sin}\alpha & -\frac{{l}_{12}}{6}\mathrm{sin}\alpha \\ \frac{{l}_{12}}{3}\mathrm{cos}\alpha & \frac{{l}_{12}}{6}\mathrm{cos}\alpha \end{array}\right)\left\{\begin{array}{c}{p}_{1}\\ {p}_{2}\end{array}\right\}\text{ }\text{，}\\ \left[G\right]=\left(\begin{array}{cc}-\frac{{l}_{12}}{6}\mathrm{sin}\alpha & -\frac{{l}_{12}}{3}\mathrm{sin}\alpha \\ \frac{{l}_{12}}{6}\mathrm{cos}\alpha & \frac{{l}_{12}}{3}\mathrm{cos}\alpha \end{array}\right)\left\{\begin{array}{c}{p}_{1}\\ {p}_{2}\end{array}\right\}\text{ }。\end{array}$

(19)

式（18）的计算结果为一个6×1的子矩阵，将其加入到DDA总体平衡方程的荷载项F_i中，即式(20)的右侧：

$\left[\begin{array}{cccc}{K}_{11}& {K}_{12}& \dots & {K}_{1n}\\ {K}_{21}& {K}_{22}& \dots & {K}_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ {K}_{n1}& {K}_{n2}& \dots & {K}_{nn}\end{array}\right]\left\{\begin{array}{c}{D}_{1}\\ {D}_{2}\\ ⋮\\ {D}_{n}\end{array}\right\}=\left\{\begin{array}{c}{F}_{1}\\ {F}_{2}\\ ⋮\\ {F}_{n}\end{array}\right\}\text{ }。$

(20)

式中，K_ij为6×6的子矩阵，用于描述块体之间的相互作用，D_i和F_i为6×1的子矩阵用以描述块体上的位移和荷载。

通过式（18）可以将渗流场计算的结果作为外荷载施加到应力应变场计算中，实现渗流场对应力应变场的影响。外荷载的施加则会改变下一时步岩体裂隙的开度，反向施加影响于渗流场，进而实现渗流场与应变场的双向耦合。

3. 算例验证

3.1 土柱入渗

该算例选取长1 m土柱计算水体从其顶部入渗的压力水头分布^[19]。在模型顶部初始施加-7 m的水头，底部施加-8 m的水头。当降雨开始后模型底部水头不发生变化，顶部施加1 m的水头作为边界条件，如图 2所示。

图 2 计算模型及边界条件

Figure 2. Computational model and boundary conditions

下载: 全尺寸图片幻灯片

该模型的具体计算参数如表 1所示，将46800 s的压力水头数值解与解析解进行对比。

表 1 模型参数

Table 1. Model parameters

初始裂隙开度/mm	θ_r	θ_s	α/m^-1	n	Δt/s
0.111	0.363	0.168	1	1.53	180

下载: 导出CSV

| 显示表格

如图 3所示，数值解与解析解的压力水头均呈现随高程下降而逐渐降低的趋势，且模型上下两端压力水头随高程变化较慢而模型中部压力水头随高程变化较快，这是由于模型上下两端设置了固定水头作为边界条件，而模型中部受边界条件影响较小。与解析解相比，数值解模型顶部压力水头上升较快，原因可能为当块体体积含水量接近饱和体积含水量时，压力水头会近似处理成饱和状态压力水头，从而导致压力水头上升较快。数值解与解析解的压力水头具有相似的变化规律，验证了该数值方法的正确性。

图 3 46800 s压力水头分布

Figure 3. Distribution of pressure head at 46800 s

下载: 全尺寸图片幻灯片

3.2 砂槽入渗试验

该算例采用砂槽入渗试验对程序的正确性进行验证，该试验常用于检验饱和-非饱和非稳定渗流算法的性能^[20]。试验区域为一个高2 m宽6 m的矩形区域。边界条件为两侧边界自由排水、底部边界不透水、上部边界降雨。地下水位在0.65 m处，在试验开始后，对该试验区域中心1 m的范围内均匀施加0.148 m/h的降雨强度，降雨时长为8 h。初始裂隙宽度为0.0003 m，饱和渗透系数为9.72×10^-5 m/s，其数值模型如图 4所示。

图 4 算例模型

Figure 4. Example model

下载: 全尺寸图片幻灯片

由于该模型是完全对称的，因此截取模型的右半段进行分析。完成建模后按照表 2参数对模型进行赋值并计算。

表 2 模型参数^[19]

Table 2. Model parameters

初始裂隙开度/mm	θ_r	θ_s	α/m^-1	n	Δt/s
1.29	0.01	0.3	3.3	4.1	3600

下载: 导出CSV

| 显示表格

降雨开始2 h后的压力水头等值线如图 5所示。

图 5 2 h压力水头等值线

Figure 5. Contours of pressure head at 2 h

下载: 全尺寸图片幻灯片

由于在模型的左上角施加降雨边界，在该边界条件下的地层压力水头较大，随着右侧边界条件变为自由边界，压力水头迅速下降。符合在降雨区域地层饱和度上升负孔隙水压消散，在未降雨区域地层状态与初始状态接近这一规律。整体模型压力水头等值线呈左高右低的趋势，且随着深度降低等值线越发趋于平缓。符合浅表地层受降雨影响较大，体积含水量上升，压力水头上升，而深处地层受到影响较小的规律。与NMM的计算结果主要在降雨边界附近存在差异，其主要原因是NMM可以通过数学片对块体进一步细分，而DDA只在块体的边界有计算点。因此DDA模拟结果的等值线分布相NMM要更为稀疏，但总体分布规律与NMM相似。该模型主要研究降雨引起地下水位的变化，对于地下水位以下部分未进行插值。

降雨开始8 h后的压力水头等值线如图 6所示。随着降雨的进行，模型的压力水头整体呈上升趋势，尤其是在降雨边界处，由2 h的-0.49 m上升到了-0.36 m。地下水位线也随降雨进行而有所升高，且呈现左高右低的趋势，符合离降雨边界较近处饱和度上升较快的规律。等值线出现偏折是由于受到边界条件影响，在降雨边界和自由边界的分界处产生计算数值上的突变。DDA计算结果压力水头的分布规律与NMM计算结果相似。选取试验区域中央截面，分析地下水位随时间上升的规律如图 7所示。

图 6 8 h压力水头等值线

Figure 6. Contours of pressure head at 8 h

下载: 全尺寸图片幻灯片

图 7 地下水位线变化

Figure 7. Variation of groundwater level line

下载: 全尺寸图片幻灯片

图 7将Vuaclin的室内试验结果和DDA数值模拟的结果进行对比。DDA数值模拟水位线上升的高度稍小于室内试验得到的结果，与试验结果相比最大误差出现在降雨4 h后，为3.74%。因此可以认为数值模拟结果与室内试验结果较为接近，进而认为该程序可以应用在非饱和瞬态渗流相关的模拟当中。

4. 工程案例

将改进后的DDA程序应用到重庆市巫山县三溪乡大梯子岩危岩体中。根据工程地质剖面图建立如图 8所示的计算模型。模型的力学参数如表 3所示。

图 8 计算模型

Figure 8. Model

下载: 全尺寸图片幻灯片

表 3 模型力学参数

Table 3. Mechanical parameters for model

材料	密度/(kg·m^-3)	弹性模量/GPa	泊松比	抗拉强度/MPa	黏聚力/MPa	摩擦角/(°)
危岩体	2400	40	0.3	10	10	35
基岩	2400	50	0.25	20	50	50

下载: 导出CSV

| 显示表格

在自然工况下对危岩体进行数值计算，结果如图 9所示。如图 9所示，在自然工况下，位移最大值点出现在危岩体的底部，其值为4.35 mm。可能由于差异风化导致危岩体底部出现凹腔，在上部自重等荷载的作用下导致危岩体在底部产生较大位移。其此时危岩体的整体稳定性系数为1.2127，危岩体整体处于相对稳定状态。

图 9 自然工况位移云图

Figure 9. Nephogram of displacement under natural condition

下载: 全尺寸图片幻灯片

对于危岩体的非饱和渗流采用Van Genuchten模型描述，具体取值如表 4所示。

表 4 渗流参数

Table 4. Seepage parameters

初始裂隙开度/mm	Δt/s	n	α/m^-1	${\theta _{\text{s}}}$	${\theta _{\text{r}}}$
0.3	1800	2	0.1	0.39	0.05

下载: 导出CSV

| 显示表格

根据当地降雨统计数据模拟20 h强度为7.07 mm/h的较长时间持续降雨，得到结果如图 10所示。

图 10 降雨工况体积含水量云图

Figure 10. Nephogram of volume water content under rainfall condition

下载: 全尺寸图片幻灯片

经过20 h降雨作用后，危岩体顶部体积含水量接近饱和体积含水量，产生饱和度大于0.95的暂态饱和区^[21]，而危岩体下部体积含水量则几乎不受到降雨的影响。由于降雨边界只在该模型顶部设置，因此底部节点的体积含水量受到降雨影响小，降雨从危岩体顶部入渗且无法渗透整个危岩体。危岩体右侧降雨入渗深度相对于左侧较深：一方面由于危岩体顶部有一定坡度，右侧高度低于左侧更有利于雨水入渗；另一方面，降雨入渗计算网格与切割岩体裂隙网格相同，岩体右侧裂隙长度较左侧更短，因此在相同的降雨条件下，危岩体右侧体积含水量更容易受到影响。而岩体顶部最右侧临空块体由于没有联通的裂隙对其进行切割，属于孤立块体，因此不在降雨入渗计算网格当中。危岩体的总位移云图如图 11所示。

图 11 降雨工况总位移云图

Figure 11. Nephogram of total displacement under rainfall condition

下载: 全尺寸图片幻灯片

危岩体最大位移为2.6 cm，出现在模型顶部的临空块体，发生较大位移的区域位于危岩体底部，其位移为1.8 cm。与自然工况相比，危岩体的位移有明显增长。危岩体顶部受到降雨影响，块体的强度有一定弱化，且顶部受到降雨形成的水头压力使其产生较大的位移。底部块体由于顶部块体的重力和变形增大，也产生了较大的位移。此时危岩体的整体稳定性系数为0.6304，危岩体整体处于不稳定状态。危岩体顶部和底部块体出现裂缝扩展、错动等大变形现象，如图 12所示。

图 12 块体大变形

Figure 12. Large deformations of block

下载: 全尺寸图片幻灯片

该模拟结果符合工程实际规律，且能用于描述块体大变形过程，验证了本文提出的方法在工程应用中的可行性。

5. 结论

针对DDA难以进行瞬态渗流场模拟这一问题对其进行二次开发，将渗流场计算引入到原始DDA当中，形成一种流固耦合的DDA方法并通过相关经典算例和工程案例进行验证。提出了一种可以反映降雨入渗过程的DDA方法，为实际工程提供新的数值模拟手段，得到以下2点结论。

（1）对土柱入渗和砂槽入渗的模拟结果表明，本文建立的非饱和瞬态流固耦合计算模型与试验结果和NMM计算模型相比具有较好的一致性。能够对降雨造成的地下水位抬升、负压水头消散等变化过程进行模拟，拓展了DDA的适用范围。

（2）对大梯子岩危岩体的模拟结果表明，降雨导致危岩体顶部出现暂态饱和区，造成岩体强度参数弱化以及裂隙水压力上升，从而产生较大的位移。在顶部块体的重力和位移作用下，危岩体底部块体产生相对较大的位移，危岩体整体处于欠稳定状态。在危岩体治理中应做好排水措施，防止降雨造成危岩体失稳崩塌。

图 1 区域地震构造与历史地震、工程场址位置及跨海通道剖面示意图

Figure 1. Regional tectonic setting and historical seismicity around sea-crossing project, site location and profile of seabed cross-section

下载: 全尺寸图片幻灯片

图 2 跨海峡通道海床场地S和P波波速结构、钻孔及观测点位置示意

Figure 2. Profiles of S- and P- wave velocities of seabed site of strait-crossing project, boreholes and typical locations at seabed surface

下载: 全尺寸图片幻灯片

图 3 海床土的动剪切模量比与阻尼比曲线

Figure 3. Curves of shear modulus reduction and damping ratio of undisturbed seabed soils

下载: 全尺寸图片幻灯片

图 4 水下深槽左侧细部网格

Figure 4. Detailed grids at the left side of the deep underwater groove

下载: 全尺寸图片幻灯片

图 5 基岩地震动的加速度时程、傅氏谱及频域累积能量比曲线(PGA = 0.1g)

Figure 5. Bedrock seismograms, Fourier spectra and cumulative energy distributions for frequency domain (level of 0.1g)

下载: 全尺寸图片幻灯片

图 6 跨海峡海床面微地形特征及水平向、双向(水平向+竖向)地震动作用下海床面各观测点的加速度反应时程(PBA_H = 0.2g)

Figure 6. Characteristics of micro-topography of seabed cross-section surface (a) and acceleration responses at typical locations of seabed surface for only horizontal shaking (b) and bidirectional (horizontal + vertical) shakings (c) (PBA_H = 0.2g)

下载: 全尺寸图片幻灯片

图 7 海床表面加速度反应传递函数(ATF)

Figure 7. ATFs at seabed surface

下载: 全尺寸图片幻灯片

图 8 跨海峡海床剖面的地表峰值加速度放大系数AF_PGA

Figure 8. Variation of AF_PGA at seabed surface

下载: 全尺寸图片幻灯片

图 9 跨海峡海床剖面的峰值加速度放大系数AF_PA

Figure 9. Variation of AF_PA at seabed section crossing the strait

下载: 全尺寸图片幻灯片

图 10 -130 m高程的峰值加速度

Figure 10. Peak accelerations at altitude of -130 m

下载: 全尺寸图片幻灯片

图 11 海床面的地震动持时延长系数D_5-95-PF

Figure 11. Variation of D_5-95-PF at seabed surface

下载: 全尺寸图片幻灯片

图 12 海床表面峰值加速度比例系数

Figure 12. Proportional coefficient of peak acceleration of seabed surface

下载: 全尺寸图片幻灯片

图 13 基岩输入双向地震动时海床地表加速度的动力系数β

Figure 13. Values of acceleration dynamic coefficient β at seabed surface under bidirectional bedrock motions

下载: 全尺寸图片幻灯片

表 1 基岩地震记录的原始信息

Table 1 Original information of bedrock earthquake recordings

台站	地震或震源位置		地震时间	地震动分量	震级M	震中距/km	基岩峰值加速度/(m·s^-2)	有效持时D_5-95 /s	第一卓越频率f_e1 /Hz	第二卓越频率f_e2 /Hz
台站	经度E/(°)	纬度N/(°)	地震时间	地震动分量	震级M	震中距/km	基岩峰值加速度/(m·s^-2)	有效持时D_5-95 /s	第一卓越频率f_e1 /Hz	第二卓越频率f_e2 /Hz
HYGH13	133.85	35.38	2016	NS	6.6	77	0.089	15.78	0.19	0.84
HYGH13	133.85	35.38	2016	UD	6.6	77	0.045	27.09	0.88	1.20
FKSH18	141.70	37.73	2021	EW	7.3	106	0.109	9.60	1.11	1.79
FKSH18	141.70	37.73	2021	UD	7.3	106	0.121	11.93	1.85	1.93
FKOH01	130.81	32.74	2016	NS	6.5	128	0.046	24.29	0.20	1.66
FKOH01	130.81	32.74	2016	UD	6.5	128	0.036	29.01	0.21	1.81

下载: 导出CSV

表 2 1D模型计算的钻剖面孔ZK-6、ZK-7、ZK-11、ZK-13处地表的AF_PGA和地表PGA_2D/1D比较

Table 2 AF_PGA obtained by 1D model and surface PGA aggregation factor PGA_2D/1D at boreholes ZK-6, ZK-7, ZK-11 and ZK-13

钻孔编号		地震记录HYGH13				地震记录FKSH18				地震记录FKOH01
		PBA_H = 0.075g		PBA_H = 0.2g		PBA_H = 0.075g		PBA_H = 0.2g		PBA_H = 0.075g		PBA_H = 0.2g
		AF_PGA(1D)	地表PGA_2D/1D	AF_PGA(1D)	地表PGA_2D/1D	AF_PGA(1D)	地表PGA_2D/1D	AF_PGA(1D)	地表PGA_2D/1D	AF_PGA(1D)	地表PGA_2D/1D	AF_PGA(1D)	地表PGA_2D/1D
海床地表水平向PGA	ZK6	2.5	1.28	1.4	1.29	2.3	1.30	1.3	1.71	2.2	1.27	1.4	1.36
	ZK7	2.6	1.04	1.6	1.35	2.5	1.00	1.6	1.32	2.1	1.10	1.5	1.13
	ZK11	2.5	0.92	1.9	1.02	2.7	1.19	1.4	1.68	2.4	1.21	1.7	1.47
	ZK13	2.4	0.92	1.6	1.12	2.2	1.27	1.3	1.40	2.3	1.17	1.5	1.20
海床地表竖向PGA	ZK6	3.2	1.41	2.6	1.25	3.1	2.16	2.3	2.18	3.0	1.57	1.9	1.55
	ZK7	3.3	1.61	2.5	2.08	3.4	1.41	2.6	1.60	2.7	2.00	1.8	1.75
	ZK11	2.9	1.00	1.8	1.06	2.7	1.48	1.9	1.51	2.6	1.77	2.0	1.34
	ZK13	3.0	0.97	1.7	1.00	2.9	1.14	1.9	1.32	2.8	1.32	1.9	1.49

下载: 导出CSV

参考文献(22)

[1]	HU J J, TAN J Y, ZHAO J X. New GMPEs for the Sagami Bay region in Japan for moderate magnitude events with emphasis on differences on site amplifications at the seafloor and land seismic stations of K-NET[J]. Bulletin of the Seismological Society of America, 2020, 110(5): 2577-2597. doi: 10.1785/0120190305
[2]	BOORE D M, SMITH C E. Analysis of earthquake recordings obtained from the Seafloor Earthquake Measurement System (SEMS) instruments deployed off the coast of southern California[J]. Bulletin of the Seismological Society of America, 1999, 89(1): 260-274. doi: 10.1785/BSSA0890010260
[3]	CHEN B K, WANG D S, LI H N, et al. Characteristics of Earthquake Ground Motion on the Seafloor[J]. Journal of Earthquake Engineering, 2015, 19(5/6): 874-904.
[4]	NAKAMURA T, NAKANO M, HAYASHIMOTO N, et al. Anomalously large seismic amplifications in the seafloor area off the Kii peninsula[J]. Marine Geophysical Research, 2014, 35(3): 255-270. doi: 10.1007/s11001-014-9211-2
[5]	DHAKAL Y P, AOI S, KUNUGI T, et al. Assessment of nonlinear site response at ocean bottom seismograph sites based on S-wave horizontal-to-vertical spectral ratios: a study at the Sagami Bay area K-NET sites in Japan[J]. Earth, Planets and Space, 2017, 69(1): 1-7. doi: 10.1186/s40623-016-0587-x
[6]	RUAN B, ZHAO K, WANG S Y, et al. Numerical modeling of seismic site effects in a shallow estuarine bay (Suai Bay, Shantou, China)[J]. Engineering Geology, 2019, 260: 105233. doi: 10.1016/j.enggeo.2019.105233
[7]	CHEN G X, RUAN B, ZHAO K, et al. Nonlinear response characteristics of undersea shield tunnel subjected to strong earthquake motions[J]. Journal of Earthquake Engineering, 2020, 24(3): 351-380. doi: 10.1080/13632469.2018.1453416
[8]	陈国兴, 岳文泽, 阮滨, 等. 金塘海峡海床地震反应特征的二维非线性分析[J]. 岩土工程学报, 2021, 43(11): 1967-1975. doi: 10.11779/CJGE202111002 CHEN Guoxing, YUE Wenze, RUAN Bin, et al. Two-dimensional nonlinear seismic response analysis for seabed site effect assessment in Jintang Strait[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(11): 1967-1975. (in Chinese) doi: 10.11779/CJGE202111002
[9]	陈国兴, 夏高旭, 王彦臻, 等. 琼州海峡海床地震反应特性的一维非线性分析[J]. 工程力学, 2022, 39(5): 75-85. CHEN Guoxing, XIA Gaoxu, WANG Yanzhen, et al. One-dimensional nonlinear seismic response analysis for seabed site effect assessment in the Qiongzhou strait[J]. Engineering Mechanics, 2022, 39(5): 75-85. (in Chinese)
[10]	CHEN G X, ZHU J, QIANG M Y, et al. Three-dimensional site characterization with borehole data-A case study of Suzhou area[J]. Engineering Geology, 2018, 234: 65-82. doi: 10.1016/j.enggeo.2017.12.019
[11]	朱姣, 陈国兴, 许汉刚, 等. 苏州第四纪深厚地层剪切波速空间变化特征及其应用[J]. 岩土工程学报, 2018, 40(4): 726-735. doi: 10.11779/CJGE201804017 ZHU Jiao, CHEN Guoxing, XU Hangang, et al. Spatial variation characteristics of shear wave velocity structure and its application to quaternary deep sediment layers in Suzhou region[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(4): 726-735. (in Chinese) doi: 10.11779/CJGE201804017
[12]	ALFELD P. A trivariate clough—tocher scheme for tetrahedral data[J]. Computer Aided Geometric Design, 1984, 1(2): 169-181. doi: 10.1016/0167-8396(84)90029-3
[13]	CHEN G X, WANG Y Z, ZHAO D F, et al. A new effective stress method for nonlinear site response analyses[J]. Earthquake Engineering & Structural Dynamics, 2021, 50(6): 1595-1611.
[14]	赵丁凤, 阮滨, 陈国兴, 等. 基于Davidenkov骨架曲线模型的修正不规则加卸载准则与等效剪应变算法及其验证[J]. 岩土工程学报, 2017, 39(5): 888-895. doi: 10.11779/CJGE201705013 ZHAO Dingfeng, RUAN Bin, CHEN Guoxing, et al. Validation of modified irregular loading-unloading rules based on Davidenkov skeleton curve and its equivalent shear strain algorithm implemented in ABAQUS[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(5): 888-895. (in Chinese) doi: 10.11779/CJGE201705013
[15]	孙田, 陈国兴, 周恩全, 等. 琼州海峡100 m以浅海洋土动剪切模量比和阻尼比试验研究[J]. 岩土工程学报, 2013, 35(增刊2): 375-382. http://www.cgejournal.com/cn/article/id/15412 SUN Tian, CHEN Guoxing, ZHOU Enquan, et al. Experimental study on dynamic shear modulus ratio and damping ratio of marine soils in Qiongzhou Strait with depth less than 100 m[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(S2): 375-382. (in Chinese) http://www.cgejournal.com/cn/article/id/15412
[16]	孙田, 陈国兴, 周恩全, 等. 深层海床粉质黏土动剪切模量和阻尼比试验研究[J]. 土木工程学报, 2012, 45(增刊1): 9-14. SUN Tian, CHEN Guoxing, ZHOU Enquan, et al. Experimental research on the dynamic shear modulus and the damping ratio of deep-seabed marine silty clay[J]. China Civil Engineering Journal, 2012, 45(S1): 9-14. (in Chinese)
[17]	章小龙, 李小军, 陈国兴, 等. 黏弹性人工边界等效荷载计算的改进方法[J]. 力学学报, 2016, 48(5): 1126-1135. ZHANG Xiaolong, LI Xiaojun, CHEN Guoxing, et al. An improved method of the calculation of equivalent nodal forces in viscous-elastic artificial boundary[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5): 1126-1135. (in Chinese)
[18]	NAKAMURA Y. What is the nakamura method?[J]. Seismological Research Letters, 2019: 1437-1443.
[19]	EVANGELIA G, IOANNIS A, GEORGE G, et al. Soil, basin and soil-building-soil interaction effects on motions of Mexico City during seven earthquakes[J]. Géotechnique, 2021: 1-9.
[20]	HOUGH S E, ALTIDOR J R, ANGLADE D, et al. Localized damage caused by topographic amplification during the 2010 M 7.0 Haiti earthquake[J]. Nature Geoscience, 2010, 3(11): 778-782. doi: 10.1038/ngeo988
[21]	CHEN G X, JIN D D, ZHU J A, et al. Nonlinear analysis on seismic site response of Fuzhou Basin, China[J]. Bulletin of the Seismological Society of America, 2015, 105(2A): 928-949. doi: 10.1785/0120140085
[22]	中国地震动参数区划图: GB 18306—2015[S]. 北京: 中国标准出版社, 2016. Seismic Ground Motion Parameters Zonation Map of China: GB 18306—2015[S]. Beijing: Standards Press of China, 2016. (in Chinese)

施引文献(5)

期刊类型引用(2)

1.	邱军领，贾玎，赖金星，唐琨杰，强磊. 含裂隙土质隧道降雨入渗双通道渗流模型. 岩土工程学报. 2025(03): 548-558 . 本站查看
2.	张小东，路喆津，刘敏，杨峰，赵炼恒. 强降雨作用下引入分形理论的浅层边坡入渗及稳定性分析. 哈尔滨工业大学学报. 2024(11): 140-150 . 百度学术