细粒铁尾矿的沉积特性与基本物理力学性质试验研究

郭晓霞; 陈之祥; 邵龙潭; 田筱剑

doi:10.11779/CJGE202007005

细粒铁尾矿的沉积特性与基本物理力学性质试验研究 English Version

郭晓霞^{1, 2,},
陈之祥^{1, 2},
邵龙潭^{1, 2, ,},
田筱剑^{1, 2}

1.
大连理工大学工程力学系,辽宁　大连116085
2.
大连理工大学工业装备结构分析国家重点实验室,辽宁　大连　116085

基金项目:

国家重点实验室自主研究课题项目 S18406

国家自然科学基金项目 51479023

国家自然科学基金项目 51309047

国家自然科学基金项目 41877251

详细信息

作者简介:
郭晓霞(1978—),女,博士,高级工程师,主要从事岩土与环境力学等方面的教学和科研工作。E-mail：hanyuer@dlut.edu.cn

通讯作者:
邵龙潭, E-mail：shaolt@hotmail.com

中图分类号: TU43
计量
- 文章访问数: 333
- HTML全文浏览量: 77
- PDF下载量: 193
出版历程
- 收稿日期: 2019-09-03
- 网络出版日期: 2022-12-05
- 刊出日期: 2020-06-30

Experimental study on sedimentary behavior and basic physical mechanical properties of fine iron tailings

GUO Xiao-xia^{1, 2,},
CHEN Zhi-xiang^{1, 2},
SHAO Long-tan^{1, 2, ,},
TIAN Xiao-jian^{1, 2}

1.
Department of Engineering Mechanics, Dalian University of Technology, Dalian 116085, China
2.
State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116085, China

摘要

摘要: 尾矿的沉积与固结过程是形成强度和确定变形量的重要依据。为探求尾矿沉积过程中应力和孔隙水压力演变规律,揭示尾矿料的强度形成机理和变形特征,结合现场溜槽试验,对细粒铁尾矿的沉积坡度、以及沉积过程中的应力、孔隙水压力和有效应力的演变规律进行了监测。同时,采用十字板剪切仪确定了沉积完成后尾矿的剪切强度。在此基础上,对不同沉积断面尾矿料的颗粒级配和渗透特性进行了试验研究,分析尾矿沉积速率的演变机理。试验结果表明：颗粒研磨极其均匀的尾矿料的堆积形式和渗透特性与尾矿料的水力特性存在差别；受尾矿颗粒相对密实度较大、颗粒分布均匀等因素影响,增加尾矿料的可透水边界较施加外荷载,更能提升尾矿料中有效应力的形成速率；尾矿料表层硬化和封闭微孔隙引发了上层滞水入渗困难,应采用必要的引水导流措施降低高势能流体对筑坝稳定性的影响。
- 细粒铁尾矿 /
- 溜槽试验 /
- 沉积特性 /
- 孔隙水压力 /
- 剪切强度
Abstract: The process of deposition and consolidation of fragmented multiphase materials is the basis and prerequisite for studying their strength and deformation. To find out the evolution laws of stress and pore water pressure during the deposition of fine iron tailings and reveal the strength formation mechanism and deformation characteristics of tailings reservoir, the sedimentary slope of tailings, stress, pore water pressure and effective stress during the deposition process are monitored by field large-scale flume tests. Meanwhile, the shear strength of tailings after deposition is determined by the vane shear apparatus. On this basis, grain-size distribution and permeability characteristics of tailings from different sedimentary sections are tested to analyze the evolution mechanism of tailings sedimentation rate. The experimental results show that the accumulation form and permeability characteristics of tailings with extremely uniform particle grinding are different from those of soils. Due to the large proportion of tailings particles and uniform distribution of particles, increasing the permeable boundary of tailings can improve the formation rate of effective stress of tailings better than applying external loads. The surface hardening of tailings and the sealing of micro-pore trigger the formation of effective stress in tailings. It is difficult to infiltrate the stagnant water in the upper layer. Necessary diversion measures should be adopted to reduce the influences of high potential fluid on the stability of dam construction.
- fine iron tailings /
- flume test /
- deposition characteristic /
- pore water pressure /
- shear strength

HTML全文

0. 引言

结构–地基相互作用问题是地基动力学、地震学和抗震工程领域的重要问题,它对于评估施加在体系上的动荷载而引发的响应具有重要意义。由于问题的复杂性,早期有关结构–地基相互作用的计算模型与算法,例如有限元法、边界元法等,大多采用均质半无限地基的假定以满足工程设计施工的需要。实际上,许多大型土木建筑以及海洋工程建筑都坐落在层状地基之上,并且层状地基表现出不同程度的各向异性属性。随着计算机技术水平的提高,横观各向同性层状地基的研究越来越受到学者们的关注。此外,以往大多相关工作主要针对频域内地基动力属性的研究,而较少针对时域内的动力响应进行求解分析。鉴于此,本文针对横观各向同性层状地基上任意形状基础进行时域内动力响应的求解,此项研究对于实际工程应用具有特殊意义,复杂且充满挑战性。

早期的分析方法利用波函数展开来获得二维、三维圆形以及刚性基础在平面SH波激发时的平面外响应^[1-2]。随后,许多学者采用分析方法研究结构–地基动力相互作用,包括地基系统与不同基础相互作用的动力响应,例如半椭圆、半球形基础等^[3-6]。可见,早期分析方法仅限于在均质半空间中具有单一形状的基础的动力响应。然而,近年来,计算机技术的发展促进了数值方法的应用,如有限元法和边界元法。利用有限元法研究地基与结构动力相互作用的关键是应用适当的人工传递边界。Lysmer^[7]首先利用有限元方法并通过设置黏性吸收边界条件计算埋置均匀半空间中的三维基础的动力响应。Tassoulas等^[8]使用半离散方法计算了埋置在分层半空间中的基础的所有刚度系数,并设置了刚度的透射边界,还详细讨论了基础深度的影响。随后,其他学者提出了各种边界条件,包括黏弹性人工边界^[9]、刚度传递边界^[10]、迭代传递边界^[11]、史密斯边界^[12]、近轴边界^[13-14]和Higdon边界^[15]等,这些学者对地基与结构之间的动力相互作用进行了深入研究。边界元法（BEM）不需要人工边界,它使用格林函数作为基本解,能够自动满足无穷远处的辐射条件,Wolf等^[16]利用格林函数求解频域内埋置基础的阻抗函数。Karabalis等^[17-18]利用边界元法求解时域内半无限地基上明置或埋置基础动力刚度系数及其在动荷载作用下的响应。Barros等^[19]使用格林函数计算了埋置层状半空间中的二维半圆形基础的动力响应。付佳等^[20]和Liang等^[21-23]利用边界元法结合分布在斜坡上的格林函数,研究了二维、半圆柱以及刚性基础—剪力墙系统在单层或多层土层中的动力响应问题。艾智勇等^[24-25]提出了层元法给出了三维层状地基频域内格林函数的半解析解,随后韩泽军等^[26-27]利用精细积分算法也给出了二维和三维层状地基任意点频域内格林函数的半解析解。王雨等^[28]基于弹性层状理论,利用传递矩阵法和有限差分法建立了能够考虑土体非均质条件下的地下管道变形计算方法,并进行了参数分析。

为了简化模型和计算,大多数学者把地基假设为各向同性属性,然而,这都是对实际情况的简化。层状地基由于其形成或施工原因,水平向和竖直向材料属性存在不同程度的差异性。因此,应进一步考虑层状地基的横观各向同性属性。Liu等^[29]利用有限元方法结合Green函数研究了埋置横观各向同性层状介质中的单桩的动力响应。Lu等^[30]对横向各向同性体的弹性动力学进行了详细研究,分别讨论了二维问题、轴对称问题和静态问题,并且将有限层元法推广到材料具有黏性的情况。李肇胤等^[31]采用针对横观各向同性层状黏弹性模型模拟半空间之上的层状场地,推导得出环形荷载作用下横观各向同性层状场地位移响应公式,并举例说明一些参数及底面边界条件对这些响应的影响。尽管考虑到地基的横观各向同性属性,但大多只得到了频域解,而较少针对时域内的解进行求解分析。

本文旨在研究时域中三维横观各向同性层状地基–基础相互作用的算法。利用傅里叶变换将广义平面应变问题的动力方程转化到频率–波数域中进行求解,然后采用逆变换得到频域内的位移响应,通过改进相减模型^[34]求解得出频率域内开挖基础动力刚度矩阵,最后通过混合变量技术^[35],构建时域内基础–地基动力方程。混合变量技术与传统的傅里叶变换相比,不仅可以处理线性体系,对非线性体系同样适用。此外,混合变量算法还能处理多自由度体系,在计算得到频域内动力刚度矩阵的基础上可以直接构建多自由度体系的时域方程,从而使动力问题得到求解。该时域方程为一阶常微分方程,可以通过精细积分时程算法^[36]进行求解。通过与已有算法的结果进行对比验证了本文算法的准确性,同时针对横观各向同性属性进行参数分析,重点关注水平向和竖直向杨氏模量的不同对时域内基础动力响应的影响。

1. 问题的描述

对于层状地基,假定每层都为均质横观各向同性属性,并且层与层之间的交界面水平,如图1所示,任意形状刚性基础坐落在或埋置于层状地基之上/内,层状地基包含 $m + 1$ 层,其中第 $m + 1$ 为刚性基岩或均质半无限空间。每层均有5个弹性控制参数,分别为水平向和竖直向的杨氏模量 $E_{HH i}$ 和 $E_{HV i}$ ,水平面和竖直面泊松比 $ν_{HH i}$ 和 $ν_{HV i}$ 以及竖直面剪切模量 $G_{HV i}$ ,此外定义每层的阻尼比为 $ζ_{i}$ ,密度为 $ρ_{i}$ ,厚度为 $h_{i} = z_{i} - z_{i - 1}$ $(i = 1, 2, 3, \dots, m)$ 。建立如图1所示三维空间直角坐标系,坐标原点O设置在地基表面,x轴和y轴为水平方向,z轴竖直向下。本文针对两种情况进行求解分析：明置基础（h=0）和埋置基础（h>0）,求解其在动荷载作用下的动力位移响应问题。

图 1 任意形状基础与横观各向同性层状地基体系

Figure 1. System of arbitrarily shaped foundation and transversely isotropic multi-layered soil

下载: 全尺寸图片幻灯片

针对明置基础情况,只需要离散基础与地基交界面,求解离散节点频率域内格林影响函数矩阵 ${\tilde{F}}_{g}$ 后便可得到基础动力方程为

$\tilde{P} = \tilde{S} \tilde{U}$ ,

(1)

式中, $\tilde{P}$ 和 $\tilde{U}$ 分别为离散节点的荷载和位移向量； $\tilde{S} = {\tilde{F}}_{g}^{- 1}$ 为基础动力刚度矩阵。

而对于开挖基础,需要将开挖基础进行离散,通过改进相减模型求解开挖基础动力刚度矩阵。对于改进相减模型,其求解理念可参见图2所示。将基础进行有限元离散,离散节点分为两部分：e表示基础外围表面节点,其包括基础与地基交界面以及自由表面；s表示基础内部节点,则此时基础动力刚度矩阵 $\tilde{S}$ 可表示为

图 2 改进相减模型示意图

Figure 2. Schematic diagram of modified subtraction model

下载: 全尺寸图片幻灯片

$\tilde{S} = [\begin{matrix} {\tilde{C}}_{ss}^{f} - {\tilde{C}}_{ss}^{c} & {\tilde{C}}_{s e}^{f} - {\tilde{C}}_{se}^{c} \\ {\tilde{C}}_{es}^{f} - {\tilde{C}}_{es}^{c} & {\tilde{C}}_{ee}^{f} - {\tilde{C}}_{ee}^{c} + {\tilde{S}}_{ee}^{g} \end{matrix}]$ 。

(2)

式中 ${\tilde{C}}^{f} = [\begin{matrix} {\tilde{C}}_{ss}^{f} & {\tilde{C}}_{se}^{f} \end{matrix}; \begin{matrix} {\tilde{C}}_{es}^{f} & {\tilde{C}}_{ee}^{f} \end{matrix}]$ 为基础有限元刚度矩阵； ${\tilde{C}}^{c} = [\begin{matrix} {\tilde{C}}_{ss}^{c} & {\tilde{C}}_{se}^{c} \end{matrix}; \begin{matrix} {\tilde{C}}_{es}^{c} & {\tilde{C}}_{ee}^{c} \end{matrix}]$ 为开挖地基有限元刚度矩阵； ${\tilde{S}}_{ee}^{g}$ 为未开挖横观各向同性层状地基在基础外围节点处的动力刚度矩阵； ${\tilde{C}}^{f}$ 和 ${\tilde{C}}^{c}$ 均可通过有限元求解其静力刚度矩阵和质量矩阵按照下式组合而成：

$\tilde{C} = \tilde{K} - ω^{2} \tilde{M}$ ,

(3)

式中, $\tilde{K}$ 和 $\tilde{M}$ 分别为有限元法计算得到的静力刚度矩阵和质量矩阵, $ω$ 为频率。

对于横观各向同性层状地基动力格林函数 ${\tilde{S}}_{ee}^{g}$ 的求解,本文在第2节内将进行简单的介绍,其详细求解过程可参见文献[27]。

2. 横观各向同性层状地基格林函数求解

2.1 状态方程的建立

由于频率–空间域中的波动方程难以直接求解,故本文利用Fourier变换将空间域中的波动方程转换到波数域进行求解,在得到波数域内位移解之后,通过Fourier逆变换得到空间域的解。设空间域和波数域内的函数分别为 $\tilde{Π} (x, y, z, ω)$ 和 $Π (α, β, z, ω)$ ,则两者之间存在Fourier变换关系为

$Π (α, β, z, ω) = \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} \tilde{Π} (x, y, z, ω) e^{- i (α x + β y)} d x d y$ ,

(4)

$\tilde{Π} (x, y, z, ω) = \frac{1}{4 π^{2}} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} Π (α, β, z, ω) e^{i (α x + β y)} d α d β$ ,

(5)

式中, $α$ 和 $β$ 分别表示沿x轴和y轴的波数, $i = \sqrt{- 1}$ 。

对于横观各向同性地基,其应力–应变关系为

(6)

式中 $D$ 为弹性矩阵； $\tilde{σ}$ 和 $\tilde{τ}$ 分别为下标标注方向的正应力和切应力； $\tilde{ε}$ 和 $\tilde{γ}$ 分别为下标标注方向的正应变和切应变,且

$\begin{array}{l} d_{11} = \frac{E_{HH} (E_{HV} - E_{HH} ν_{HV}^{2})}{(1 + ν_{HH}) (E_{HV} (1 - ν_{HV}) - 2 E_{HH} ν_{HV}^{2})}, \\ d_{12} = \frac{E_{HH} (E_{HV} ν_{HH} + E_{HH} ν_{HV}^{2})}{(1 + ν_{H H}) (E_{HV} (1 - ν_{HV}) - 2 E_{HH} ν_{H V}^{2})}, \\ d_{13} = \frac{E_{HH} E_{HV} ν_{HV}}{E_{HV} (1 - ν_{HV}) - 2 E_{HH} ν_{HV}^{2}}, \\ d_{33} = \frac{E_{HV}^{2} (1 - ν_{HH})}{E_{HV} (1 - ν_{HV}) - 2 E_{HH} ν_{HV}^{2}}, d_{44} = G_{HV} 。 \end{array}}$

(7)

位移–应变关系为

${\begin{cases} \begin{matrix} {\tilde{ε}}_{x} = \frac{\partial \tilde{u}}{\partial x}, & {\tilde{ε}}_{y} = \frac{\partial \tilde{v}}{\partial y}, & {\tilde{ε}}_{z} = \frac{\partial \tilde{w}}{\partial z}, \end{matrix} \\ \begin{matrix} {\tilde{γ}}_{y z} = \frac{\partial \tilde{w}}{\partial y} + \frac{\partial \tilde{v}}{\partial z}, & {\tilde{γ}}_{x z} = \frac{\partial \tilde{w}}{\partial x} + \frac{\partial \tilde{u}}{\partial z}, & {\tilde{γ}}_{x y} = \frac{\partial \tilde{v}}{\partial x} + \frac{\partial \tilde{u}}{\partial y} 。 \end{matrix} \end{cases}$

(8)

式中, $\tilde{u}$ , $\tilde{v}$ 和 $\tilde{w}$ 分别为沿x,y和z方向的位移幅值。

频率–空间域内三维波动问题的控制方程为

${\begin{cases} \frac{\partial {\tilde{σ}}_{x}}{\partial x} + \frac{\partial {\tilde{τ}}_{x y}}{\partial y} + \frac{\partial {\tilde{τ}}_{x z}}{\partial z} = ρ \frac{\partial^{2} \tilde{u}}{\partial t^{2}}, \\ \frac{\partial {\tilde{τ}}_{x y}}{\partial x} + \frac{\partial {\tilde{σ}}_{y}}{\partial y} + \frac{\partial {\tilde{τ}}_{y z}}{\partial z} = ρ \frac{\partial^{2} \tilde{v}}{\partial t^{2}}, \\ \frac{\partial {\tilde{τ}}_{x z}}{\partial x} + \frac{\partial {\tilde{τ}}_{y z}}{\partial y} + \frac{\partial {\tilde{σ}}_{z}}{\partial z} = ρ \frac{\partial^{2} \tilde{w}}{\partial t^{2}} 。 \end{cases}$

(9)

利用式（4）的Fourier变换可将频率–空间域内横观各向同性地基的控制方程（式（9））转换到频率–波数域中：

$K_{22} u^{″} + (K_{21} - K_{12}) u^{'} - (K_{11} - ρ ω^{2} I) u = 0$ 。

(10)

式中, $I$ 为单位矩阵, $u^{'}$ 和 $u^{″}$ 分别为位移向量对竖向坐标z的一阶和二阶偏微分导数,系数矩阵分别为

$\begin{array}{l} K_{22} = D_{z z}, K_{21} = K_{12}^{H} = i α D_{x z} + i β D_{y z}, \\ K_{11} = α^{2} D_{x x} + β^{2} D_{y y} + α β (D_{x y} + D_{y z}), \end{array}$

且

$\begin{array}{l} D_{x x} = [\begin{matrix} d_{11} & 0 & 0 \\ 0 & d_{66} & 0 \\ 0 & 0 & d_{55} \end{matrix}], D_{y y} = [\begin{matrix} d_{66} & 0 & 0 \\ 0 & d_{22} & 0 \\ 0 & 0 & d_{44} \end{matrix}], \\ D_{z z} = [\begin{matrix} d_{55} & 0 & 0 \\ 0 & d_{44} & 0 \\ 0 & 0 & d_{33} \end{matrix}], D_{x y} = D_{y x}^{Τ} = [\begin{matrix} 0 & d_{12} & 0 \\ d_{66} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], \\ D_{x z} = D_{z x}^{Τ} = [\begin{matrix} 0 & 0 & d_{13} \\ 0 & 0 & 0 \\ d_{55} & 0 & 0 \end{matrix}], D_{y z} = D_{z y}^{Τ} [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & d_{23} \\ 0 & d_{44} & 0 \end{matrix}] 。 \end{array}$

式（10）为二阶常微分方程,为方便计算,引入对偶应力变量 $p$ ,且令其满足 $p = - (K_{22} u^{'} + K_{21} u)$ ,则式（10）可以转化为一阶常微分方程

$X^{'} = M X$ ,

(11)

式中, $X = {\begin{matrix} u & p \end{matrix}}^{T}$ 为混合向量,系数矩阵 $M$ 写成相应的分块矩阵的形式为

$\begin{matrix} M = [\begin{matrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix}], \\ {\begin{cases} M_{11} = - K_{22}^{- 1} K_{21}, M_{12} = - K_{22}^{- 1}, \\ M_{21} = - K_{11} + K_{12} K_{22}^{- 1} K_{21} + ρ ω^{2} I, \\ M_{22} = K_{12} K_{22}^{- 1} 。 \end{cases} \end{matrix}$ 。

(12)

状态方程中,对于固定的波数 $α$ 和 $β$ ,方程中混合向量 $X$ 仅是竖直坐标 $z$ 的函数。Zhong等^[32]利用精细积分算法对保守体系的方程（11）进行求解,可以得到任意层地基上下表面的位移和应力关系为

${\begin{matrix} p_{a} \\ - p_{b} \end{matrix}} = [\begin{matrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{matrix}] {\begin{matrix} u_{a} \\ u_{b} \end{matrix}}$ ,

(13)

式中, $u_{a}$ 和 $p_{a}$ 分别为上表面的位移和对偶应力向量, $u_{b}$ 和 $p_{b}$ 分别为下表面的位移和对偶应力向量；系数矩阵 $K_{aa}$ , $K_{ab}$ , $K_{ba}$ 和 $K_{bb}$ 为矩阵 $M_{11}$ , $M_{12}$ , $M_{21}$ 和 $M_{22}$ 的函数,其具体求解过程可参见文献[32]。

2.2 边界条件

对于多层地基,相邻层的交界面处需满足位移和应力的连续条件,即

$u_{k}^{+} = u_{k}^{-}, p_{k}^{+} = p_{k}^{-}, (k = 1, 2, 3, \dots, l - 1, l)$ 。

(14)

对于层状地基底部（第 $m + 1$ 层）为刚性基岩的情况,需要满足m层底部位移为 $0$ 的边界条件：

$u_{m} = 0$ 。

(15)

而对于层状地基下卧半无限空间的情况,需要考虑无限远处的辐射边界条件。首先求解均质半无限地基特征矩阵 $M$ 的特征值问题：

$M Φ = Φ Λ$ ,

(16)

式中, $Λ$ 为特征值矩阵, $Φ$ 为特征向量矩阵,二者写成分块矩阵的形式为

$Λ = [\begin{matrix} λ_{i} \\ - λ_{i} \end{matrix}], Φ = [\begin{matrix} Φ_{11} & Φ_{12} \\ Φ_{21} & Φ_{22} \end{matrix}]$ ,

(17)

式中, $λ_{i}$ 包含3个实部为正值的特征值。

引入

$b = Φ^{- 1} X$ ,

(18)

此时方程（11）可以写成

$\dot{b} = Λ b, b (z) = [\begin{matrix} exp(λ_{i} z) \\ exp(- λ_{i} z) \end{matrix}] {\begin{matrix} c_{1} \\ c_{2} \end{matrix}}$ 。

(19)

式中, $c_{1}$ 和 $c_{2}$ 为积分常数,向量 $b (z)$ 为一个 $(6 \times 1)$ 列向量。由式（19）可以看出, $b (z)$ 的分量表示向上或者向下传播的波。由于矩阵 $Λ$ 和 $Φ$ 已经按照式（17）的特定方式分块, $b (z)$ 的前3个分量表示向上传播的波,而后3个分量表示向下传播的波。对于均质半无限地基,因为没有向上传播的波,位移向量 $u$ 和应力向量 $p$ 必须保证为有限值,因此需要有 $c_{1} = 0$ 。

将方程（19）代入方程（18）中可得

$X = [\begin{matrix} Φ_{11} & Φ_{12} \\ Φ_{21} & Φ_{22} \end{matrix}] {\begin{cases} c_{1} exp(λ_{i} z) \\ c_{2} exp(- λ_{i} z) \end{cases}}$ 。

(20)

应用方程（20）可以得到半无限地基表面的边界条件为

$p_{m} = K_{\infty} u_{m}, K_{\infty} = Φ_{22} Φ_{12}^{- 1}$ 。

(21)

2.3 频率波数域内动力刚度矩阵

假设基础外围离散节点共有k层,将基础深度范围内的地基划分为 $k - 1$ 层,利用2.2节求解任意相邻两层节点之间的系数矩阵 $K_{aa}^{i}$ , $K_{ab}^{i}$ , $K_{ba}^{i}$ 和 $K_{bb}^{i} (i = 1, 2, 3, \dots, k - 1)$ ,结合2.2节的边界条件,便可组装得到k层节点在频率–波数域内的动力刚度矩阵：

${\begin{matrix} P_{1} \\ P_{2} \\ ⋮ \\ P_{k - 1} \\ P_{k} \end{matrix}} = [\begin{matrix} K_{aa}^{1} & K_{ab}^{1} & 0 & 0 & 0 \\ K_{ba}^{1} & K_{bb}^{1} + K_{aa}^{2} & K_{ab}^{2} & 0 & 0 \\ 0 & K_{ba}^{2} & ⋱ & K_{ab}^{k - 2} & 0 \\ 0 & 0 & K_{ab}^{k - 2} & K_{bb}^{k - 2} + K_{aa}^{k - 1} & K_{ab}^{k - 1} \\ 0 & 0 & 0 & K_{ba}^{k - 1} & K_{bb}^{k - 1} + K_{\infty} \end{matrix}] {\begin{matrix} U_{1} \\ U_{2} \\ ⋮ \\ U_{k - 1} \\ U_{k} \end{matrix}} 。$

(22)

式中 $P_{i} (i = 1, 2, 3, \dots, k)$ 表示第i层节点处荷载向量,且 $P_{i} = - p_{b}^{i - 1} + p_{a}^{i}$ ,特殊情况 $P_{1} = p_{a}^{1}$ ； $U_{i} (i = 1, 2, 3, \dots, k)$ 表示第i层节点处位移向量,且 $U_{i} = u_{b}^{i - 1} = u_{a}^{i}$ ； $K_{\infty}$ 表示第k层节点底部层状地基在该层面的动力刚度矩阵。特殊情况,对于基础底部即为刚性基岩的情况,可直接令 $U_{k} = 0$ 即可得到频率–波数域内的动力刚度矩阵。

2.4 位移格林函数求解

由于集中荷载作用下位移响应会出现奇异情况,本文格林函数的求解均转化为均布荷载作用。对于明置基础,可以直接将基础–地基交界面离散为圆盘单元进行求解计算,每个圆盘中心即为求解节点。对于开挖基础,在引用改进相减模型时,为与有限元模型进行耦合,将基础在水平面方向按照正方形单元 $(Δ a \times Δ a)$ 进行离散,如图3所示,将离散节点的控制面积转化为半径为 $Δ$ 的圆盘,每个圆盘的圆心即为计算节点。所有节点的控制面积共包括4种情况,即A/4,2A/4,3A/4和A（A为正方形单元面积）,并将该控制面积等面积转化为半径为 $Δ_{1}$ , $Δ_{2}$ , $Δ_{3}$ 和 $Δ_{4}$ 的圆盘进行计算。最大圆盘半径的选取由计算截止频率 $ω_{\max}$ 决定,即需满足在波长 $2π c_{s} / ω_{\max}$ 范围内至少需要6个圆盘单元,每个圆盘单元的直径应不大于 $π c_{s} / (3 ω_{\max})$ ,则最大圆盘半径应为 $π c_{s} / (6 ω_{\max})$ 。

图 3 开挖基础水平面单元离散

Figure 3. Discretization of embedded foundation in horizontal plane

下载: 全尺寸图片幻灯片

在得到频率–波数域内k层节点动力刚度矩阵（式（22））之后,利用式（5）的Fourier逆变换求解空间域内任意点的位移响应,即

${\tilde{U}}_{i j} (x, y, z, ω) = \frac{1}{{4π}^{2}} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} U_{r s} (α, β, z, ω) e^{i (α x + β y)} d α d β,$

(23)

式中, ${\tilde{U}}_{i j} (i, j = 1, 2, 3, \dots, n)$ （n为需要计算格林函数的节点数）为频率空间域内在j点施加荷载在i点产生的位移；第i点处于开挖基础的第r层节点内,第j点处于开挖基础的第s层节点内, $U_{r s}$ 表示频率波数域内s层点施加荷载在r层点产生的位移,其可通过式（22）求逆获得。显而易见, ${\tilde{U}}_{i j}$ 和 $U_{r s}$ 均为 $3 \times 3$ 的矩阵,可分别记作

${\tilde{U}}_{i j} = {[\begin{matrix} {\tilde{U}}_{x x} & {\tilde{U}}_{x y} & {\tilde{U}}_{x z} \\ {\tilde{U}}_{y x} & {\tilde{U}}_{y y} & {\tilde{U}}_{y z} \\ {\tilde{U}}_{z x} & {\tilde{U}}_{z y} & {\tilde{U}}_{z z} \end{matrix}]}_{i j}, U_{r s} = {[\begin{matrix} U_{r r} & U_{r θ} & U_{r z} \\ U_{θ r} & U_{θ θ} & U_{θ z} \\ U_{z r} & U_{z θ} & U_{z z} \end{matrix}]}_{r s} 。$

(24)

式中, ${\tilde{U}}_{ξ η} (ξ, η = x, y, z)$ 和 $U_{ξ η} (ξ, η = r, θ, z)$ 表示在 $η$ 方向施加荷载产生的 $ξ$ 方向的位移。

值得注意的是,式（23）是一个双重无限积分,计算量比较大,本文采用文献[33]的方法,利用横观各向同性介质的轴对称特性,将双重积分简化为一重积分进行求解,即令 $α = 0$ 且 $β = κ$ ,则均布荷载作用下的位移格林函数可表示为

${\tilde{U}}_{i j} (x, y, z, ω) = \frac{{\tilde{P}}_{j}}{π Δ_{k}} R (θ) \int_{0}^{+ \infty} {\tilde{G}}_{i j} R^{Τ} (θ) J_{1} (κ Δ_{k}) d κ$ ,

(25)

式中 $Δ_{j}$ 表示施加均布荷载圆盘的半径； $J_{n} (#)$ 表示第一类n阶贝塞尔函数； $R (θ)$ 为

$R (θ) = [\begin{matrix} \sin θ & \cos θ & 0 \\ - \cos θ & \sin θ & 0 \\ 0 & 0 & 1 \end{matrix}]$ ,

(26)

$θ$ 为i点相对于j点在柱坐标系下的方位角,如图3所示。 ${\tilde{G}}_{i j}$ 中的元素分别为

$\begin{array}{l} {\tilde{G}}_{11} = \frac{1}{2} ([J_{0} (κ r) - J_{2} (κ r)] U_{r r} + [J_{0} (κ r) + J_{2} (κ r)] U_{θ θ}), \\ {\tilde{G}}_{22} = \frac{1}{2} ([J_{0} (κ r) - J_{2} (κ r)] U_{θ θ} + [J_{0} (κ r) + J_{2} (κ r)] U_{r r}), \\ {\tilde{G}}_{23} = {iJ}_{1} (κ r) U_{r z}, G_{32} = {iJ}_{1} (κ r) U_{z r}, \\ {\tilde{G}}_{33} = J_{0} (κ r) U_{z z}, {\tilde{G}}_{12} = {\tilde{G}}_{13} = {\tilde{G}}_{21} = {\tilde{G}}_{31} = 0 。 \end{array}}$

(27)

式（25）虽然是个无限积分,但是由于被积函数具有良好的收敛性,可采用有限积分代替无限积分进行求解,得到精度较高的计算结果。

至此,式（1）中的 ${\tilde{F}}_{g}^{- 1}$ 和（2）中的 ${\tilde{S}}_{aa}^{g}$ 可以直接通过格林函数组装而成,进而可以得到柔性基础的动力刚度矩阵。

3. 刚性基础动力刚度矩阵

对于刚性基础,假设基础和地基完全接触,则节点位移与刚体位移必须满足条件：

$\tilde{U} = Ξ {\tilde{U}}_{0}$ 。

(28)

式中 ${\tilde{U}}_{0}$ 为刚性基础的6个刚体自由度位移幅值,包括3个平动分量 $(\begin{matrix} {\tilde{U}}_{x} & {\tilde{U}}_{y} & {\tilde{U}}_{z} \end{matrix})$ 和3个转动分量 $(\begin{matrix} {\tilde{U}}_{θ x} & {\tilde{U}}_{θ y} & {\tilde{U}}_{θ z} \end{matrix})$ ；矩阵 $Ξ$ 为相应的形函数,并且有

${\begin{cases} {\tilde{U}}_{0} = {\begin{matrix} \begin{matrix} {\tilde{U}}_{x} & {\tilde{U}}_{y} & {\tilde{U}}_{z} \end{matrix} & {\tilde{U}}_{θ x} & {\tilde{U}}_{θ y} & {\tilde{U}}_{θ z} \end{matrix}}^{T} \\ Ξ = {[\begin{matrix} Ξ_{1} & Ξ_{2} & \dots & Ξ_{n} \end{matrix}]}^{T} \\ Ξ_{i} = {[\begin{matrix} 1 & 0 & 0 & 0 & z_{i} - z_{0} & - (y_{i} - y_{0}) \\ 0 & 1 & 0 & - (z_{i} - z_{0}) & 0 & x_{i} - x_{0} \\ 0 & 0 & 1 & y_{i} - y_{0} & - (x_{i} - x_{0}) & 0 \end{matrix}]}^{T} (i = 1, 2, 3, \dots, n) \end{cases} 。$

(29)

式中, $(x_{0}, y_{0}, z_{0})$ 为基础参考点的坐标,对于明置基础选取为基础形心位置,开挖基础选取为基础地面的形心位置； $(x_{i}, y_{i}, z_{i}) (i = 1, 2, 3, \dots, n)$ 为节点坐标。

根据力的平衡条件,作用在刚性基础上的集中荷载和离散节点处施加的均布荷载幅值满足

${\tilde{P}}_{0} = Ξ^{Τ} \tilde{P}$ ,

(30)

式中, ${\tilde{P}}_{0} = {\begin{matrix} \begin{matrix} {\tilde{F}}_{x} & {\tilde{F}}_{y} & {\tilde{F}}_{z} \end{matrix} & {\tilde{M}}_{x} & {\tilde{M}}_{y} & {\tilde{M}}_{z} \end{matrix}}^{Τ}$ 。

联立式（1）,（2）,（28）和（30）可得

${\tilde{P}}_{0} = {\tilde{S}}_{f} {\tilde{U}}_{0} 且 {\tilde{S}}_{f} = Ξ^{Τ} \tilde{S} Ξ$ 。

(31)

矩阵 ${\tilde{S}}_{f}$ 即为刚性基础的动力刚度矩阵,其可以写成如下形式：

${\tilde{S}}_{f} = [\begin{matrix} S_{11} & 0 & 0 & 0 & S_{15} & 0 \\ 0 & S_{22} & 0 & S_{24} & 0 & 0 \\ 0 & 0 & S_{33} & 0 & 0 & 0 \\ 0 & S_{42} & 0 & S_{44} & 0 & 0 \\ S_{51} & 0 & 0 & 0 & S_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{66} \end{matrix}]$ ,

(32)

式中 $S_{i i}$ $(i = 1, 2, 3)$ 为刚性基础平动刚度系数； $S_{i i}$ $(i = 4, 5, 6)$ 为刚性基础转动刚度系数； $S_{15} = S_{51}$ 和 $S_{24} = S_{42}$ 为平动与摇摆的耦合项。

4. 动刚度的有理拟合

在得到频域–空间域内基础动力刚度的基础上,采用混合变量方法求解得到基础在外荷载作用下的时域解。假设通过1～3节求解得到了w个一系列频率下基础离散的动力刚度矩阵 $\tilde{S} (ω)$ 或 ${\tilde{S}}_{f} (ω)$ ,利用Ruge等^[35]提出的有理函数进行拟合,即

$\tilde{S} (ω) 或 {\tilde{S}}_{f} (ω) = \frac{Y (i ω)}{T (i ω)} = \frac{Y_{0} + i ω Y_{1} + \dots + {(i ω)}^{δ + 1} Y_{δ + 1}}{T_{0} + i ω T_{1} + \dots + {(i ω)}^{δ} T_{δ}}$ 。

(33)

式中 $δ$ 为有理函数的阶次；系数矩阵 $Y_{i} (i = 1, 2, \dots,$ $δ + 1)$ 和 $T_{i} (i = 1, 2, \dots, δ)$ 可以通过离散动力刚度矩阵求得,具体求解过程可参见文献[35]。公式（33）中 $δ$ 阶的有理函数分子中 $i ω$ 的最高次为 $δ + 1$ ,分母中为 $i ω$ 的最高次为 $δ$ ,因此总是可以将该有理函数分解为 $i ω$ 的线性函数项和另一个 $δ - 1$ 阶的有理函数项,如此每一步分解得到的有理函数项都可以产生一个新的内部变量,从而引入一个既包含荷载,又包含位移的混合变量,详细的分解过程可参见附录A。

最终频域和时域内混合变量方程为

$\begin{matrix} A Z + i ω B Z = P & （频域） & , \end{matrix}$

(34)

$\begin{matrix} Α Z (t) + B \dot{Z} (t) = P (t) & （时域） & 。 \end{matrix}$

(35)

式中 A和B为通过频率–空间域内基础动力刚度矩阵计算得到的系数矩阵； $Z$ 为混合变量,其中既含有位移向量,又包括虚拟荷载向量； $P$ 为外荷载向量。方程（35）是一阶常微分方程,本文采用精细积分时程法对其进行求解,不仅可以保证较高的精度,而且有较高的求解效率。

5. 数值算例

为了验证本文算法的准确性,本文选择了两个算例加以说明：第1个算例选取横观各向同性均质半无限地基情况进行计算；第2个算例选取多层横观各向同性半无限地基上埋置方形基础的情况进行验证。在验证本文算法准确的基础上,针对横观各向同性层状地基的各向异性属性（水平向和竖直向杨氏模量的不同）进行参数分析,研究其对基础动力响应时域解的影响。

5.1 刚性圆形基础置于横观各向同性均质半无限地基之上

为验证本文算法的准确性,现选取刚性圆形基础置于均质半无限地基情况进行动力响应计算,圆形基础的半径为r₀。地基材料属性为横观各向同性属性,其弹性参数为水平向和竖直向的杨氏模量 $E_{HH} =$ $3.864 G_{HV}$ 和 $E_{HV} = 2 .863 G_{HV}$ ,水平面和竖直面泊松比 $ν_{HH} = 0.185$ 和 $ν_{HV} = 0.301$ 以及竖直面剪切模量 $G_{HV}$ ,密度为 $ρ$ ,阻尼比为0.001。基础与地基交界面离散成200个圆盘单元,每个圆盘单元的圆心即为一个计算节点,圆盘单元的半径为0.0705r₀。利用本文提出算法求解水平、竖直、摇摆和扭转四个模态的动力刚度系数,并利用有理函数进行拟合,为说明有理函数拟合随阶数M的收敛性,选取竖直向动力刚度系数拟合结果进行说明。针对M=1～6分别进行拟合,并将计算得到的动力刚度矩阵系数按照下式进行无量纲化：

$\begin{matrix} S_{i i} = K_{static} [k (a_{0}) + a_{0} c (a_{0})] (i = 1, 3, 4, 6), \\ a_{0} = ω \sqrt{G_{HV} / ρ} / r_{0} \end{matrix}$ 。

(36)

式中, $K_{static}$ 为基础对应的静力刚度矩阵系数,如图4所示。

图 4 有理函数拟合阶数收敛性分析

Figure 4. Convergence analysis of order for rational function fitting

下载: 全尺寸图片幻灯片

由图4可见,有理函数拟合随阶数M有良好的收敛性,在M=4之后拟合结果基本相同。其余模态动力刚度系数拟合结果如图5所示,由图5可见,动力刚度系数的实部和虚部均完美拟合。采用拟合好的动力刚度矩阵有理函数表达式建立时域内的运动方程,分别在刚性基础中心施加水平、竖直、摇摆和扭转4个方向的动荷载,外部激励荷载定义为

图 5 频域–空间域动力刚度矩阵拟合曲线

Figure 5. Fitted curves of dynamic stiffness matrices in frequency-spatial domain

下载: 全尺寸图片幻灯片

$F (t) = F_{0} π \sin (\frac{π}{d} \sqrt{\frac{G_{HV}}{ρ}} t) \exp [- \frac{1}{4} {(\frac{t - 10}{2})}^{2}]$ 。

(37)

计算时间总长为20r₀/c_s,时间步长为0.01r₀/c_s,计算得到的动力响应如图6所示。

图 6 均质横观各向同性半无限地基上圆形基础时域内动力响应

Figure 6. Dynamic responses versus time for a circular foundation on a transversely isotropic half-space

下载: 全尺寸图片幻灯片

目前为止,尚未找到任何文献时域计算结果予以对比,本文采用有限单元法进行求解,选取地基尺寸水平向半径为20r₀,竖直向厚度也为20r₀,从而保证没有任何运动波经有限元模型边界反射回到基础。由图6可见,两种结果基本完全拟合,从而验证了本文针对横观各向同性层状地基上刚性基础动力响应时域解的准确性。值得注意的是,本文算法只需要在基础范围内进行离散,自由度数较少,因此有较高的求解效率,以本算例而言,通过本文算法进行计算只需要30 min,而有限元算法需要12 h以上。此外,本文算法的计算耗时主要集中在频域内动力刚度矩阵的求解上,时域内的计算耗时很少,对于相同模型不同荷载情况的计算,总耗时会更少。

5.2 刚性方形基础埋置于多层横观各向同性地基

为进一步说明本文对于多层横观各向同性地基的适用性,选取宽度为2a的刚性方形基础埋置于多层横观各向同性地基内的情况予以验证。本算例中多层地基的材料参数见表1所示,地基的阻尼比均为0.001。埋置基础的深度h/a=0.5,基础在水平面方向离散为12×12共计144个方形单元,竖直向离散为5层节点。基础中心参考点坐标选为（0,0,h）,在基础参考点施加公式（37）所示的激励荷载,利用混合变量方法求解得到水平、竖直、摇摆和扭转4个模态的动力响应,时间步长取为0.01a/c_s $(c_{s} = \sqrt{G_{0} / ρ_{0}})$ ,计算时长为20a/c_s。同样将计算结果与有限单元法计算结果进行对比（见图7）。由图7可见,两种结果基本完全拟合,从而进一步验证了本文算法对于多层地基情况的适用性。

表 1 横观各向同性层状地基材料属性

Table 1. Material properties of transversely isotropic layered soil

层号	$\frac{E_{HH}}{G_{0}}$	$\frac{E_{HV}}{G_{0}}$	$ν_{HH}$	$ν_{HV}$	$\frac{ρ}{ρ_{0}}$	$\frac{G_{HV}}{G_{0}}$	H (a)
1	3.864	2.863	0.185	0.301	1.0	1.0	2.0
2	7.728	5.726	0.250	0.330	1.4	1.1	1.0
3	9.660	7.1575	0.250	0.250	2.0	1.2	半无限

下载: 导出CSV

| 显示表格

图 7 多层地基上埋置方形基础时域内动力响应

Figure 7. Dynamic responses versus time for a square foundation embedded in a multi-layered soil

下载: 全尺寸图片幻灯片

5.3 各向异性属性对基础动力响应的影响

以往研究主要针对各向同性层状地基进行求解分析,而对于横观各向同性层状地基的研究主要局限于频率域内,对于时域内动力响应的求解较少。而层状地基由于其形成原因以及施工工艺等因素,水平向和竖直向的弹性参数往往不同,而在实际工程和研究中往往忽略这一点。为研究地基各向异性特性对基础动力响应的影响,选取刚性圆盘基础置于单层半无限地基的算例进行说明。基础半径为1.0 m,地基材料参数见表2所示,在保持水平向杨氏模量不变的情况下,按照比值χ=E_HH/E_HV更改竖直向杨氏模量E_HV的值,分别对χ=0.5,1,2三种情况进行计算。分别在刚性基础上施加水平、竖直、摇摆和扭转4个方向的激励荷载,外部激励荷载定义见图8所示,时间步长 $Δ t$ 为5.79×10^-5 s,计算得到的4种模态的动力响应结果见图9所示。由图9可见,比值χ对于基础脉冲响应有显著的影响,尤其对于动力响应的幅值,随着χ的增大,4个位移响应的幅值均明显的增大,反之则减小。对于水平向和扭转弯矩荷载作用下,位移响应幅值发生的时刻均与χ无关,均在同一个时间点产生最大幅值；不同的是,对于竖直向和摇摆向,不仅幅值随χ有较大的变化,对于幅值发生时间点也有变换,随着χ的增大,该时间点均有一定的延迟。此外,随着χ的增大,即随着竖直向杨氏模量的减小,竖直向和摇摆向位移在基础回归到起始位置后逐渐发生负向位移。

表 2 横观各向同性层状半无限地基材料属性

Table 2. Material properties of transversely soil isotropic layered half-space

层号	E_HH/GPa	$ν_{HH}$	$ν_{HV}$	ρ/(kg·m^-3)	H/m
1	0.50	0.33	0.37	2100	1.0
2	1.00	0.32	0.35	2300	半无限

下载: 导出CSV

| 显示表格

图 8 三角形荷载时程

Figure 8. Triangular load versus time

下载: 全尺寸图片幻灯片

图 9 各向异性特性对基础动力响应的影响

Figure 9. Influences of horizontal and vertical Young’s moduli on impulse response for a massless foundation

下载: 全尺寸图片幻灯片

6. 结论

时域中横观各向同性层状地基–结构相互作用求解是近年来国内外学者十分关注的课题,与各向同性地基相比更符合实际工程的需要。本文提出了一种横观各向同性多层地基上任意形状基础动力响应时域解的求解方法。首先波动方程的解在频率–波数域中给出,随后采用Fourier逆变换和改进相减模型得到频率–空间域内的动力刚度矩阵,最后通过使用混合变量公式实现从频域到时域的转换,不仅将标量有理拟合扩展到多自由度问题,并且频域中的线性公式很容易转换为时域一阶常微分方程。

（1）数值算例验证了本文算法的准确性,该算法不仅可以求解横观各向同性层状半无限地基上明置基础动力响应时域解,而且可以求解开挖基础的情况。

（2）参数分析表明了,地基水平向和竖直向材料属性的差异对于基础动力响应有着显著的影响,尤其是对于在相同荷载作用下的位移幅值,对于四种模态具有相同的影响规律,均随着竖直向杨氏模量的减小而增大。

（3）此外,地基水平向和竖直向材料属性的差异对于位移响应幅值对应的时间点,四种不同的模态影响有所差别,其中水平向和扭转向没有影响,而对于竖直向和摇摆向,随着n的增大均有所延迟。

附录：

以 $δ = 3$ 为例说明详细分解过程以及动力方程的构建（假设 $λ = i ω$ ）：

由式（33）可得初始方程为

$\begin{array}{l} \frac{Y^{(0)}}{T^{(0)}} {\tilde{U}}_{0} = {\tilde{P}}_{0}, \\ Y^{(0)} = Y_{0} + λ Y_{1} + λ^{2} Y_{2} + λ^{3} Y_{3} + λ^{4} Y_{4}, \\ T^{(0)} = T_{0} + λ T_{1} + λ^{2} T_{2} + λ^{3} T_{3} 。 \end{array}}$

(38)

第一次分解结果为

$\begin{array}{l} \frac{Y^{(0)}}{T^{(0)}} \overset{!}{=} Q^{(0)} + \frac{R^{(0)}}{T^{(0)}}, \\ R^{(0)} = R_{0}^{(0)} + λ R_{1}^{(0)} + λ^{2} R_{2}^{(0)}, \\ Q^{(0)} = Q_{0}^{(0)} + λ Q_{1}^{(0)} 。 \end{array}}$

(39)

其中,系数矩阵 $Q_{k}^{(0)}$ 和 $R_{k}^{(0)}$ 可以通过针对 $λ^{k}$ 的待定系数法计算得到：

$Y^{(0)} = Q^{(0)} T^{(0)} + R^{(0)}$ ,

(40)

$\begin{array}{l} λ^{4} : Q_{1}^{(0)} T_{3}^{(0)} = Y_{4}^{(0)} \to Q_{1}^{(0)}, \\ λ^{3} : Q_{0}^{(0)} T_{3}^{(0)} + Q_{1}^{(0)} T_{2}^{(0)} = Y_{3}^{(0)} \to Q_{0}^{(0)}, \\ λ^{2} : Q_{0}^{(0)} T_{2}^{(0)} + Q_{1}^{(0)} T_{1}^{(0)} + R_{2}^{(0)} = Y_{2}^{(0)} \to R_{2}^{(0)}, \\ λ^{1} : Q_{0}^{(0)} T_{1}^{(0)} + Q_{1}^{(0)} T_{0}^{(0)} + R_{1}^{(0)} = Y_{1}^{(0)} \to R_{1}^{(0)}, \\ λ^{0} : Q_{0}^{(0)} T_{0}^{(0)} + R_{0}^{(0)} = Y_{0}^{(0)} \to R_{0}^{(0)} 。 \end{array}}$

(41)

将分解得到的有理函数方程（2）代入方程（1）中可得

$\begin{matrix} [Q^{(0)} + \frac{R^{(0)}}{T^{(0)}}] {\tilde{U}}_{0} = {\tilde{P}}_{0} \to Q^{(0)} {\tilde{U}}_{0} + P_{1} = {\tilde{P}}_{0} 且 \\ P_{1} = \frac{R^{(0)}}{T^{(0)}} {\tilde{U}}_{0} \end{matrix}$ 。

(42)

式中, ${\tilde{P}}_{1}$ 是一个新的变量,为了进行下一步的分解,构成一个新的有理函数方程,并使其满足以下条件：

$\begin{array}{l} \frac{Y^{(1)}}{T^{(1)}} P_{1} = {\tilde{U}}_{0}, \\ Y^{(1)} \equiv T^{(0)} = T_{0}^{(0)} + λ T_{1}^{(0)} + λ^{2} T_{2}^{(0)} + λ^{3} T_{3}^{(0)}, Y_{k}^{(1)} = T_{k}^{(0)}, \\ T^{(1)} \equiv R^{(0)} = R_{0}^{(0)} + λ R_{1}^{(0)} + λ^{2} R_{2}^{(0)}, T_{k}^{(1)} = R_{k}^{(0)} 。 \end{array}}$

(43)

第二次分解结果为

$\begin{array}{l} \frac{Y^{(1)}}{T^{(1)}} = Q^{(1)} + \frac{R^{(1)}}{T^{(1)}}, \\ \begin{array}{l} Q^{(1)} = Q_{0}^{(1)} + λ Q_{1}^{(1)}, \\ R^{(1)} = R_{0}^{(1)} + λ R_{1}^{(1)} 。 \end{array} \end{array}}$

(44)

系数矩阵 $Q_{k}^{(1)}$ 和 $R_{k}^{(1)}$ 可以通过针对 $λ^{k}$ 的待定系数法求解得到

$Q^{(1)} T^{(1)} + R^{(1)} = Y^{(1)}$ ,

(45)

$\begin{array}{l} λ^{3} : Q_{1}^{(1)} T_{2}^{(1)} = Y_{3}^{(1)} \to Q_{1}^{(1)}, \\ λ^{2} : Q_{0}^{(1)} T_{2}^{(1)} + Q_{1}^{(1)} T_{1}^{(1)} = Y_{2}^{(1)} \to Q_{0}^{(1)}, \\ λ^{1} : Q_{0}^{(1)} T_{2}^{(1)} + Q_{1}^{(1)} T_{0}^{(1)} + R_{1}^{(1)} = Y_{1}^{(1)} \to R_{1}^{(1)}, \\ λ^{0} : Q_{0}^{(1)} T_{0}^{(1)} + R_{0}^{(1)} = Y_{0}^{(1)} \to R_{0}^{(1)} 。 \end{array}}$

(46)

将式（7）代入式（6）中,并定义新的变量 $U_{1}$ 可得

$[Q^{(1)} + \frac{R^{(1)}}{T^{(1)}}] P_{1} = {\tilde{U}}_{0} \to Q^{(1)} P_{1} + U_{1} = {\tilde{U}}_{0}$ 。

此时有

$U_{1} = \frac{R^{(1)}}{T^{(1)}} P_{1}$ 。

(47)

为进行下一步的分解,再次建立一个新的有理函数表达式,使其满足如下条件：

$\begin{array}{l} \frac{Y^{(2)}}{T^{(2)}} U_{1} = P_{1}, \\ Y^{(2)} \equiv T^{(1)} = R_{0}^{(0)} + λ R_{1}^{(0)} + λ^{2} R_{2}^{(0)}, Y_{k}^{(2)} = R_{k}^{(0)}, \\ T^{(2)} \equiv R^{(1)} = R_{0}^{(1)} + λ R_{1}^{(1)}, T_{k}^{(2)} = R_{k}^{(1)} 。 \end{array}}$

(48)

第三次分解结果为

$\begin{array}{l} \frac{Y^{(2)}}{T^{(2)}} = Q^{(2)} + \frac{R^{(2)}}{T^{(2)}}, \\ \begin{array}{l} Q^{(2)} = Q_{0}^{(2)} + λ Q_{1}^{(2)}, \\ R^{(2)} = R_{0}^{(2)} 。 \end{array} \end{array}}$

(49)

类似前两步分解,系数矩阵 $Q_{k}^{(2)}$ 和 $R_{0}^{(2)}$ 可以通过针对 $λ^{k}$ 的待定系数法求解得到

$Y^{(2)} = Q^{(2)} T^{(2)} + R^{(2)}$ ,

(50)

$\begin{array}{l} λ^{2} : Q_{1}^{(2)} T_{1}^{(2)} = Y_{2}^{(2)} \to Q_{1}^{(2)}, \\ λ^{1} : Q_{0}^{(2)} T_{1}^{(2)} + Q_{1}^{(2)} T_{0}^{(2)} = Y_{1}^{(2)} \to Q_{0}^{(2)}, \\ λ^{0} : Q_{0}^{(2)} T_{0}^{(2)} + R_{0}^{(2)} = Y_{0}^{(2)} \to R_{0}^{(2)} 。 \end{array}}$

(51)

定义一个新的变量 $P_{2}$ ,同时将式（12）代入式（11）可得

$[Q^{(2)} + \frac{R^{(2)}}{T^{(2)}}] U_{1} = P_{1} \to Q^{(2)} U_{1} + P_{2} = P_{1}$ 。

此时有

$P_{2} = \frac{R^{(2)}}{T^{(2)}} U_{1}$ 。

(52)

为进行下一步的分解,定义新的变量对 $Y^{(3)}$ 和 $T^{(3)}$ 建立新的有理函数方程,并使其满足如下条件：

$\begin{array}{l} \frac{Y^{(3)}}{T^{(3)}} P_{2} = U_{1}, \\ Y^{(3)} = T^{(2)} = R_{0}^{(1)} + λ R_{1}^{(1)}, Y_{k}^{(3)} = R_{k}^{(1)}, \\ T^{(3)} = R^{(2)} = R_{0}^{(2)}, T_{0}^{(3)} = R_{0}^{(2)} 。 \end{array}}$

(53)

其中, $Y^{(3)} / T^{(3)}$ 只包含两项,即

$\frac{Y^{(3)}}{T^{(3)}} = Q_{0}^{(3)} + λ Q_{1}^{(3)}$

,此时有

$Q_{0}^{(3)} = R_{0}^{(1)} / R_{0}^{(2)}, Q_{1}^{(3)} = R_{1}^{(1)} / R_{0}^{(2)}$ ,

(54)

$[Q_{0}^{(3)} + λ Q_{1}^{(3)}] {\tilde{P}}_{2} = {\tilde{U}}_{1}$ 。

(55)

此时通过方程（5）,（10）,（15）和（18）便可得到频域和时域内的混合变量方程（34）和（35）。其中的系数矩阵表达式为

$\begin{matrix} A = (\begin{matrix} Q_{0}^{(0)} & I & 0 & 0 \\ - I & Q_{0}^{(1)} & I & 0 \\ 0 & - I & Q_{0}^{(2)} & I \\ 0 & 0 & - I & Q_{0}^{(3)} \end{matrix}), Z = (\begin{matrix} {\tilde{U}}_{0} \\ P_{1} \\ U_{1} \\ P_{2} \end{matrix}), P = (\begin{matrix} {\tilde{P}}_{0} \\ 0 \\ 0 \\ 0 \end{matrix}) \\ B = diag {\begin{matrix} Q_{1}^{(0)} & Q_{1}^{(1)} & Q_{1}^{(2)} & Q_{1}^{(3)} \end{matrix}} \end{matrix}$ 。

(56)

图 1 试验布置

Figure 1. Test arrangement

下载: 全尺寸图片幻灯片

图 2 现场流槽试验结果

Figure 2. Results of flume tests on tailings

下载: 全尺寸图片幻灯片

图 3 尾矿流槽中的非均匀流态

Figure 3. Non-uniform flow patterns in tailings flume

下载: 全尺寸图片幻灯片

图 4 不同断面处尾矿料的颗粒级配

Figure 4. Grain-size distribution curves of different sections

下载: 全尺寸图片幻灯片

图 5 现场流槽十字板剪切试验结果

Figure 5. Results of vane shear tests of tailings flume

下载: 全尺寸图片幻灯片

图 6 基于十字板抗剪强度的内摩擦角

Figure 6. Friction angles based on vane shear strength tests

下载: 全尺寸图片幻灯片

表 1 试验过程中各阶段

Table 1 Schedule of flume tests on tailings

序号	时间	事项
1	2011年10月13日	开始排放尾矿
2	2011年10月13日	初始阶段放矿
3	2011年10月14日—15日	第二阶段放矿
4	2011年10月16日—17日	第三阶段放矿
5	2011年10月23日—24日	第四阶段放矿
6	2011年11月08日	降雨
7	2011年11月20日	监测到尾矿冻结
8	2011年12月10日	停止数据监测

下载: 导出CSV

表 2 不同断面处尾矿料的颗粒组成

Table 2 Grain composition of different sections

断面/m	土粒相对密实度	颗粒组成/%					不均匀系数C_u	曲率系数C_c	液限含水率/%	塑限含水率/%	塑性指数I_p	土质类型
断面/m	土粒相对密实度	2~0.5 mm	0.5~0.25 mm	0.25~0.075 mm	0.075~0.005 mm	<0.005 mm	不均匀系数C_u	曲率系数C_c	液限含水率/%	塑限含水率/%	塑性指数I_p	土质类型
5	2.877	0.24	6.69	28.13	64.70	0.25	4.47	0.57	11.68	21.55	9.87	低液限尾粉土
10	2.851	0.36	9.36	35.86	53.92	0.50	4.56	1.07	13.26	23.17	9.91	低液限尾粉土
20	2.846	0.06	4.61	39.56	54.68	1.10	4.42	1.00	13.30	22.77	9.47	低液限尾粉土
40	2.838	0.02	2.54	25.34	71.90	0.21	3.53	1.07	14.49	25.80	11.31	低液限尾粉质黏土
50	2.861	0.02	0.45	23.28	75.86	0.40	3.94	1.08	15.76	27.80	12.04	低液限尾粉质黏土
60	2.839	0.01	0.24	14.60	84.40	0.75	2.87	0.89	14.57	25.73	11.16	低液限尾粉质黏土
70	2.806	0.02	0.25	15.74	82.50	1.50	3.57	1.29	12.32	25.99	13.67	低液限尾粉质黏土
80	2.868	0	0.37	12.27	86.66	0.70	3.50	1.59	14.54	25.85	11.31	低液限尾粉质黏土
90	2.864	0	0.22	13.89	84.19	1.70	3.15	0.95	15.36	26.23	10.87	低液限尾粉质黏土
100	2.877	0	0.16	11.31	86.04	2.50	4.11	0.97	15.57	27.42	11.85	低液限尾粉质黏土
120	2.884	0	0.42	5.72	91.66	2.20	3.08	0.90	17.09	27.23	10.14	低液限尾粉质黏土
140	2.878	0	0.10	5.84	91.57	2.50	3.61	0.79	16.61	29.05	12.44	低液限尾粉质黏土

下载: 导出CSV

表 3 不同断面处尾矿料的渗透系数

Table 3 Permeability coefficient of different sections

断面位置/m	土粒相对密实度	干密度/(g·cm^-3)	最大干密度/(g·cm^-3)	最优含水率/%	渗透系数/(cm·s^-1)
5	2.88	1.513	1.878	14.5	8.85×10^-5
15	2.84	1.517	1.868	15.5	2.63×10^-4
25	2.84	1.553	1.915	13.7	6.02×10^-5
40	2.84	1.498	1.795	15.9	3.28×10^-4
80	2.87	1.475	1.828	15.8	1.27×10^-4
120	2.89	1.516	1.819	16.7	7.75×10^-5

下载: 导出CSV

表 4 各测试点土样的基本物性参数

Table 4 Basic physical properties of soil samples at different test points

取土点/m	颗粒相对密实度	天然含水率/%	天然密度/(g·cm^-3)
10	2.83	27.4	2.032
20	2.84	27.8	2.027
40	2.84	27.0	1.985
60	2.84	27.7	1.965
70	2.85	23.1	1.939
80	2.87	32.9	1.960
90	2.86	26.6	1.984
100	2.87	27.8	1.922

下载: 导出CSV

参考文献(25)

[1]	张力霆, 齐清兰, 李强, 等. 尾矿库坝体溃决演进规律的模型试验研究[J]. 水利学报, 2016, 47(2): 229-235. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201602013.htm ZHANG Li-ting, QI Qing-lan, LI Qiang, et al. Experimental model study on dam break and evolution law of tailings pond[J]. Journal of Hydraulic Engineering, 2016, 47(2): 229-235. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201602013.htm
[2]	张力霆. 尾矿库溃坝研究综述[J]. 水利学报, 2013, 44(5): 594-600. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201305015.htm ZHANG Li-ting. Summary on the dam-break of tailing pond[J]. Journal of Hydraulic Engineering, 2013, 44(5): 594-600. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201305015.htm
[3]	陈生水. 尾矿库安全评价存在的问题与对策[J]. 岩土工程学报, 2016, 38(10): 1869-1873. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201610020.htm CHEN Sheng-shui. Problems and countermeasures of safety evaluation of tailing pond[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(10): 1869-1873. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201610020.htm
[4]	孙从露, 邵龙潭, 郭晓霞. 尾矿砂不同频率动三轴试验研究[J]. 矿业研究与开发, 2015, 35(11): 69-73. https://www.cnki.com.cn/Article/CJFDTOTAL-KYYK201511017.htm SUN Cong-lu, SHAO Long-tan, GUO Xiao-xia. Study on dynamic triaxial tests of tailings sand under different vibration frequencies[J]. Mining Research and Development, 2015, 35(11): 69-73. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-KYYK201511017.htm
[5]	于斯滢, 邵龙潭, 刘士乙. 基于有限元极限平衡法的尾矿坝坝体稳定分析[J]. 岩土力学, 2013, 34(4): 1185-1190. YU Si-ying, SHAO Long-tan, LIU Shi-yi. Stability analysis of tailings dam based on finite element limit equilibrium method[J]. Rock and Soil Mechanics, 2013, 34(4): 1185-1190. (in Chinese)
[6]	刘建民, 邱月, 郭婷婷, 等. 饱和重塑黄土液化孔压增长模型的试验研究[J]. 科学技术与工程, 2019, 19(19): 226-232. LIU Jian-min, QIU Yue, GUO Ting-ting, et al. Experimental study on pore pressure growth model of saturated remolded loess sample during liquefaction[J]. Science Technology and Engineering, 2019, 19(19): 226-232. (in Chinese)
[7]	GARING C, CHALENDAR J A D, VOLTOLINI M, et al. Pore-scale capillary pressure analysis using multi-scale X-ray micromotography[J]. Advances in Water Resources, 2017, 104: 223-241.
[8]	刘叔灼, 李慧子, 单毅, 等. 基于能量法的尾矿土动孔压模型研究[J]. 岩土工程学报, 2016, 38(11): 2051-2058. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201611019.htm LIU Shu-zhuo, LI Hui-zi, SHAN Yi, et al. Energy method for analyzing dynamic pore water pressure model for tailing soil[J]. Chinese Journal of Geotechnical Engineering, 2017, 104: 223-241. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201611019.htm
[9]	巫尚蔚, 杨春和, 张超, 等. 尾矿浆沉积室内模拟试验[J]. 工程科学学报, 2017, 39(10): 1485-1492. https://www.cnki.com.cn/Article/CJFDTOTAL-BJKD201710004.htm WU Shang-wei, YANG Chun-he, ZHANG Chao, et al. Indoor scale-down test of tailings[J]. Chinese Journal of Engineering, 2017, 39(10): 1485-1492. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BJKD201710004.htm
[10]	刘洋, 赵学同, 吴顺川. 快速冲填尾矿库静力液化分析与数值模拟[J]. 岩石力学与工程学报, 2014, 33(6): 1158-1168. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201406009.htm LIU Yang, ZHAO Xue-tong, WU Shun-chuan. Analysis of static liquefaction and numerical simulation for tailings pond under high deposting rates[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(6): 1158-1168. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201406009.htm
[11]	LIANG W, XIONG D S. Research on the solution to coefficient of earth pressure at rest using effective stress spade[J]. Journal of Railway Engineering Society, 2017, 34(7): 5-9, 35.
[12]	罗战友, 夏建中, 龚晓南, 等. 考虑孔压消散的静压单桩挤土位移场研究[J]. 岩石力学与工程学报, 2014, 33(增刊1): 2765-2772. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2014S1025.htm LUO Zhan-you, XIA Jian-zhong, GONG Xiao-nan. Study of compacting soil displacements around jacked single pile based on excess pore pressure dissipation[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S1): 2765-2772. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2014S1025.htm
[13]	郑彬彬. 高浓度尾矿上游式堆坝基础性问题研究及坝体稳定性分析[D]. 重庆: 重庆大学, 2017. ZHENG Bin-bin. Research on the Basic Issues for Upstream Method of High Thickened Tailings and Stability Analysis of Tailings Dam[D]. Chongqing: Chongqing University, 2017. (in Chinese)
[14]	张鹏伟, 吴辉, 胡黎明, 等. 铁矿尾矿料力学特性及坝体变形稳定性研究[J]. 工程地质学报, 2015, 23(6): 1189-1195. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201506025.htm ZHANG Peng-wei, WU Hui, HU Li-ming, et al. Mechanical characteristics of iron mine tailing materials and analysis on deformation and stability of tailing dam[J]. Journal of Engineering Geology, 2015, 23(6): 1189-1195. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201506025.htm
[15]	成词峰, 徐颖, 郑庭. 混合式堆坝尾矿料的基本性能及沉积规律分析[J]. 青岛理工大学学报, 2018, 39(2): 36-41. https://www.cnki.com.cn/Article/CJFDTOTAL-QDJG201802008.htm CHENG Ci-feng, XU Ying, ZHENG Ting. Properties of tailings and the law of deposition about tailings dam constructed by upstream and centerline methods[J]. Journal of Qingdao University of Technology, 2018, 39(2): 36-41. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-QDJG201802008.htm
[16]	尹光志, 张千贵, 魏作安, 等. 孔隙水运移特性及对尾矿细观结构作用机制试验研究[J]. 岩石力学与工程学报, 2012, 31(1): 71-79. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201201011.htm YIN Guang-zhi, ZHANG Qian-gui, WEI Zuo-an, et al. Experimental study of migration characteristics of pore water and its effect on meso-structure of tailings[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(1): 71-79. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201201011.htm
[17]	魏作安, 杨永浩, 赵怀军, 等. 小打鹅尾矿库尾矿堆积坝稳定性研究[J]. 东北大学学报(自然科学版), 2016, 37(4): 589-593. https://www.cnki.com.cn/Article/CJFDTOTAL-DBDX201604028.htm WEI Zuo-an, YANG Yong-hao, ZHAO Huai-jun, et al. Stability of tailings dam of Xiaodae Tailings pond[J]. Journal of Northeastern University(Natural Science), 2016, 37(4): 589-593. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DBDX201604028.htm
[18]	冯彦芳, 李顺群, 陈之祥, 等. 基于土体各向异性的雨水入渗渗井试验研究与验证[J]. 长江科学院院报, 2019, 36(3): 110-115. https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201903022.htm FENG Yan-fang, LI Shun-qun, CHEN Zhi-xiang, et al. Experimental study and verification of rainwater infiltration well based on soil anisotropy[J]. Journal of Yangtze River Scientific Research Institute, 2019, 36(3): 110-115. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201903022.htm
[19]	李顺群, 贾红晶, 王杏杏, 等. 轴平移技术在基质吸力测控中的局限性和误差分析[J]. 岩土力学, 2016, 37(11): 3089-3095, 3252. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201611007.htm LI Shun-qun, JIA Hong-jing, WANG Xing-xing, et al. Limitation and error analysis of axis translation technique for measuring and controlling matric suction[J]. Rock and Soil Mechanics, 2016, 37(11): 3089-3095, 3252. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201611007.htm
[20]	王建州, 刘书幸, 周国庆, 等. 深季节冻土地区基坑工程水平冻胀力试验研究[J]. 中国矿业大学学报, 2018, 47(4): 815-821. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGKD201804015.htm WANG Jian-zhou, LIU Shu-xing, ZHOU Guo-qing, et al. Model experiment on frost-heave force of foundation pit at deepseasonal frozen regions[J]. Journal of China University of Mining & Technology, 2018, 47(4): 815-821. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGKD201804015.htm
[21]	蔡正银, 李小梅, 韩林, 等. 考虑级配和颗粒破碎影响的堆石料临界状态研究[J]. 岩土工程学报, 2016, 38(8): 1357-1364. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201608001.htm CAI Zheng-yin, LI Xiao-mei, HAN Lin, et al. Critical state of rockfill materials considering particle gradation and breakage[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(8): 1357-1364. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201608001.htm
[22]	ONO K, YU Y, SAWADA Y, et al. Lateral force–displacement prediction for buried pipe under different effective stress condition[J]. International Journal of Geotechnical Engineering, 2017(2): 1-9.
[23]	CHUN L I, ZHAO F L, XIU Z G, et al. Testing investigation on effective stress increment of unsaturated sandy soils[J]. Journal of Northeastern University, 2017, 38(8): 1158-1162.
[24]	张常光, 赵均海, 朱倩. 非饱和土抗剪强度公式分类及总结[J]. 建筑科学与工程学报, 2012, 29(2): 74-82. https://www.cnki.com.cn/Article/CJFDTOTAL-XBJG201202016.htm ZHANG Chang-guang, ZHAO Jun-hai, ZHU Qian. Classification and summary of shear strength formulae for unsaturated soils[J]. Journal of Architecture and Civil Engineering, 2012, 29(2): 74-82. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-XBJG201202016.htm
[25]	陈星欣, 白冰, 于涛, 等. 粒径和渗流速度对多孔介质中悬浮颗粒迁移和沉积特性的耦合影响[J]. 岩石力学与工程学报, 2013, 32(增刊1): 2840-2845. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2013S1033.htm CHEN Xing-xin, BAI Bing, YU Tao, et al. Coupled effects of particle size and flow rate on characteristics of particle transportation and deposition in porous media[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(S1): 2840-2845. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2013S1033.htm

施引文献(9)

期刊类型引用(4)

1.	朱宝龙，李凯，林其，于时恩. 层状地基对脉冲风洞天平基础振动特征影响分析. 地震工程与工程振动. 2024(06): 125-137 . 百度学术
2.	张聪，冯忠居，王富春，马晓谦，陈慧芸. 强震作用下嵌岩群桩时程响应振动台试验. 应用基础与工程科学学报. 2023(03): 703-714 . 百度学术
3.	王志宇，唐贞云，杜修力. 时域稳定的基础频响离散有理近似参数识别方法. 岩土工程学报. 2021(09): 1708-1714 . 本站查看
4.	冯忠居，张聪，何静斌，董芸秀，袁枫斌. 强震作用下嵌岩单桩时程响应振动台试验. 岩土力学. 2021(12): 3227-3237 . 百度学术