横观各向同性层状地基上任意形状刚性基础动力响应求解与分析

杨林青; 韩泽军; 林皋; 周小文; 潘宗泽

doi:10.11779/CJGE202007009

横观各向同性层状地基上任意形状刚性基础动力响应求解与分析 English Version

杨林青^{1, 2,},
韩泽军^1, ,,
林皋³,
周小文¹,
潘宗泽¹

1.
华南理工大学土木与交通学院,广东　广州　510641
2.
广东技术师范大学天河学院,广东　广州　510540
3.
大连理工大学建设工程学部水利工程学院,辽宁　大连　116024

基金项目:

国家自然科学基金青年基金项目 51508203

中央高校基本科研业务费面上项目 2018MS61

详细信息

作者简介:
杨林青(1986—),讲师,主要从事结构–复杂层状地基动力相互作用研究工作。E-mail：ylq0313@126.com

通讯作者:
韩泽军, E-mail：ctzjhan@scut.edu.cn

中图分类号: TU435;O302
计量
- 文章访问数: 276
- HTML全文浏览量: 23
- PDF下载量: 123
出版历程
- 收稿日期: 2019-06-02
- 网络出版日期: 2022-12-05
- 刊出日期: 2020-06-30

Solution and analysis for dynamic response of arbitrarily shaped rigid foundation on transversely isotropic layered soil

1.
School of Civil Engineering and Transportation, South China University of Technology, Guangzhou 510641, China
2.
Tianhe College of Guangdong Polytechnic Normal University, Guangzhou 510540, China
3.
Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, China

摘要

摘要: 基础动力响应的求解对于更好地理解结构–地基相互作用具有重要的意义,而以往的相关研究主要求解各向同性层状地基内的频域解,对于横观各向同性层状地基,尤其是时域解目前较少有相关报道。针对横观各向同性层状地基,利用积分变换方法、混合变量技术以及改进相减模型求解得到了任意形状基础动力响应的时域解。数值算例验证了算法的准确性,同时针对横观各向同性属性中水平向和竖直向材料参数的不同对时域内基础动力响应的影响进行了参数分析。结果表明,地基的非均质特性对于基础动力响应有显著的影响,在实际工程中应予以考虑,以求得到更加精确的结果。
- 横观各向同性层状地基 /
- 任意形状基础 /
- 积分变换方法 /
- 混合变量技术 /
- 时域解
Abstract: The solution of dynamic response of foundation is of great significance for a better understanding of the structure-soil interaction. However, the previous studies have mainly solved the frequency-domain solution for the isotropic layered soil, and few reports have been presented on the transverse isotropic layered soil, especially the solution in time domain. In the study, the dynamic response solution for the arbitrarily shaped foundation resting on (or) embedded in transversely isotropic multi-layered soil in time domain is obtained by using the integral transformation method, mixed variable technique and modified subtraction model. Numerical examples verify the accuracy of the proposed algorithm. Then a parameter study is carried out to analyze the influences of the difference between the material parameters in horizontal and vertical directions on the dynamic response of foundation in time domain. The results show that the heterogeneity of the soil has a significant effect on the dynamic response of the foundation, which should be considered in practical projects to obtain more accurate results.
- transversely isotropic layered soil /
- arbitrarily shaped foundation /
- integral transformation method /
- mixed variable technology /
- time-domain solution

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. System of arbitrarily shaped foundation and transversely isotropic multi-layered soil

下载: 全尺寸图片幻灯片

图 2 改进相减模型示意图

Figure 2. Schematic diagram of modified subtraction model

下载: 全尺寸图片幻灯片

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

Figure 3. Discretization of embedded foundation in horizontal plane

下载: 全尺寸图片幻灯片

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

Figure 4. Convergence analysis of order for rational function fitting

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

图 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

下载: 全尺寸图片幻灯片

表 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

表 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

参考文献(36)

[1]	LUCO J E. Dynamic interaction of a shear wall with the soil[J]. Journal of the Engineering Mechanics Division, 1969, 95(2): 333-346. doi: 10.1061/JMCEA3.0001101
[2]	TRIFUNAC M D. Interaction of a shear wall with the soil for incident plane SH waves[J]. Bulletin of the Seismological Society of America, 1972, 62(1): 63-83. doi: 10.1785/BSSA0620010063
[3]	WONG H L, TRIFUNAC M D. Interaction of a shear wall with the soil for incident plane SH waves: elliptical rigid foundation[J]. Bulletin of the Seismological Society of America, 1974, 64(6): 1825-1842. doi: 10.1785/BSSA0640061825
[4]	LUCO J E. Torsional response of structures for SH waves: the case of hemispherical foundations[J]. Bulletin of the Seismological Society of America, 1976, 66(1): 109-123. doi: 10.1785/BSSA0660010109
[5]	APSEL R J, LUCO J E. Torsional response of rigid embedded foundation[J]. Journal of the Engineering Mechanics Division, 1976, 102(6): 957-970. doi: 10.1061/JMCEA3.0002193
[6]	LEE V W. Investigation of Three-Dimensional Soil-Structure Interaction[R]. Los Angeles: Department of Civil Engineering, University of Southern California, 1979.
[7]	LYSMER J. Finite dynamic model for infinite media[J]. Journal of the Engineering Mechanics Division, 1969, 95(4): 859-878. doi: 10.1061/JMCEA3.0001144
[8]	TASSOULAS J L, KAUSEL E. On the effect of the rigid sidewall on the dynamic stiffness of embedded circular footings[J]. Earthquake Engineering & Structural Dynamics, 1983, 11(3): 403-414.
[9]	周兴涛, 盛谦, 冷先伦, 等. 颗粒离散单元法地震动力时程计算黏性人工边界及其应用[J]. 岩石力学与工程学报, 2017, 36(4): 928-939. doi: 10.13722/j.cnki.jrme.2016.0426 ZHOU Xing-tao, SHENG Qian, LENG Xian-lun, et al. Viscous artificial boundary for seismic dynamic time-history analysis with granular discrete element method and its application[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(4): 928-939. (in Chinese) doi: 10.13722/j.cnki.jrme.2016.0426
[10]	WAAS G. Earthquake Vibration Effects and Abatement for Military Facilities[R]. Miss: Vicksburg, 1972.
[11]	LIAO Z F, HUANG K, YANG B, et al. A transmitting boundary for transient wave analyses[J]. Science in China (Series A—Mathematics, Physics, Astronomy & Technological Science), 1969, 27(10): 1063-1076.
[12]	SMITH W D. A nonreflecting plane boundary for wave propagation problems[J]. Journal of Computational Physics, 1974, 15(4): 492-503. doi: 10.1016/0021-9991(74)90075-8
[13]	CLAYTON R, ENGQUIST B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bulletin of the Seismological Society of America, 1977, 67(6): 1529-1540. doi: 10.1785/BSSA0670061529
[14]	ENGQUIST B, MAJDA A. Absorbing boundary conditions for the numerical simulation of waves[J]. Mathematics of Computation, 1977, 31(139): 629-651. doi: 10.1090/S0025-5718-1977-0436612-4
[15]	HIGDON R L. Absorbing boundary conditions for difference approximations to the multidimensional wave equation[J]. Mathematics of Computation, 1986, 47(176): 437-459.
[16]	WOLF J P, DARBRE G R. Dynamic-stiffness matrix of soil by the boundary-element method: embedded foundation[J]. Earthquake Engineering & Structural Dynamics, 1984, 12(3): 385-400.
[17]	KARABALIS D, BESKOS D. Dynamic response of 3-D rigid surface foundations by time domain boundary element method[J]. Earthquake Engineering & Structural Dynamics, 1984, 12(1): 73-93.
[18]	KARABALIS D, BESKOS D. Dynamic response of 3-D embedded foundations by the boundary element method[J]. Computer Methods in Applied Mechanics and Engineering, 1986, 56(1): 91-119. doi: 10.1016/0045-7825(86)90138-6
[19]	BARROS F, LUCO J E. Dynamic response of a two-dimensional semi-circular foundation embedded in a layered viscoelastic half-space[J]. Soil Dynamics and Earthquake Engineering, 1995, 14(1): 45-57. doi: 10.1016/0267-7261(94)00021-8
[20]	付佳, 梁建文. 场地动力特性对土–基础相互作用的影响[J]. 地震工程与工程振动, 2016, 36(1): 141-147. https://www.cnki.com.cn/Article/CJFDTOTAL-DGGC201601018.htm FU Jia, LIANG Jian-wen. Influence of site dynamic characteristics on soil-foundation interaction[J]. Earthquake Engineering Engineering Dynamics, 2016, 36(1): 141-147. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DGGC201601018.htm
[21]	LIANG Jian-wen, FU Jia, TODOROVSKA M I, et al. Effects of the site dynamic characteristics on soil-structure interaction: I incident SH-waves[J]. Soil Dynamics and Earthquake Engineering, 2013, 44: 27-37. doi: 10.1016/j.soildyn.2012.08.013
[22]	LIANG Jian-wen, FU Jia, TODOROVSKA M I, et al. In-plane soil-structure interaction in layered, fluid-saturated, poroelastic half-space: I structural response[J]. Soil Dynamics and Earthquake Engineering, 2016, 81: 84-111. doi: 10.1016/j.soildyn.2015.10.018
[23]	LIANG Jian-wen, FU Jia, TODOROVSKA M I, et al. In-plane soil-structure interaction in layered, fluid-saturated, poroelastic half-space: II pore pressure and volumetric strain[J]. Soil Dynamics and Earthquake Engineering, 2017, 92: 585-595. doi: 10.1016/j.soildyn.2016.10.012
[24]	艾智勇, 曾凯, 曾文泽. 层状地基三维问题的解析层元解[J]. 岩土工程学报, 2012, 34(6): 1154-1158. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201206031.htm AI Zhi-yong, ZENG Kai, ZENG Wen-ze. Analytical layer-element solution for three-dimensional problem of multilayered foundation[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(6): 1154-1158. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201206031.htm
[25]	艾智勇, 胡亚东. 3D横观各向同性地基非耦合解析层元[J]. 岩土工程学报, 2013, 35(2): 717-720. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2013S2131.htm AI Zhi-yong, HU Ya-dong. Uncoupled analytical layer-element for 3D transversely isotropic multilayered foundation[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(2): 717-720. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2013S2131.htm
[26]	韩泽军, 林皋, 周小文, 等. 横观各向同性层状地基上埋置刚性条带基础动力刚度矩阵求解[J]. 岩土工程学报, 2016, 38(6): 1117-1124. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201606019.htm HAN Ze-jun, LIN Gao, ZHOU Xiao-wen, et al. Solution of dynamic stiffness matrix for rigid strip foundations embedded in layered transversely isotropic soil[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(6): 1117-1124. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201606019.htm
[27]	韩泽军, 林皋, 周小文. 三维横观各向同性层状地基任意点格林函数求解[J]. 岩土工程学报, 2016, 38(12): 2218-2225. doi: 10.11779/CJGE201612010 HAN Ze-jun, LIN Gao, ZHOU Xiao-wen. Solution to Green's functions for arbitrary points in 3D cross-anisotropic multi-layered soil[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(12): 2218-2225. (in Chinese) doi: 10.11779/CJGE201612010
[28]	王雨, 陈文化, 王凯旋. 地层非均质成层性对地下管道的变形影响[J]. 岩土工程学报, 2018, 40(5): 900-909. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805018.htm WANG Yu, CHEN Wen-hua, WANG Kai-xuan. Effect of soil stratification on deformation of buried pipelines[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(5): 900-909. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805018.htm
[29]	LIU W, NOVAK M. Dynamic response of single piles embedded in transversely isotropic layered media[J]. Earthquake Engineering & Structural Dynamics, 1994, 23(11): 1239-1257.
[30]	LU W D, YANG Z W. Finite layer analysis for semi-infinite soils (II): special cases and numerical examples[J]. Applied Mathematics and Mechanics, 1992, 13(11): 987-993.
[31]	李肇胤, 周锋, 陈镕, 等. 横观各向同性层状场地受垂直圆盘状简谐荷载时的响应分析[J]. 岩石力学与工程学报, 2000(5): 608-612. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200005015.htm LI Zhao-yin, ZHOU Feng, CHEN Rong, et al. Respon se analysis on transversely isotropic elastic strata to harmonic disk loads[J]. Chinese Journal of Rock Mechanics and Engineering, 2000(5): 608-612. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200005015.htm
[32]	GAO Q, ZHONG W X, HOWSON W P. A precise method for solving wave propagation problems in layered anisotropic media[J]. Wave Motion, 2004, 40: 191-207.
[33]	ANDERSEN L, CLAUSEN J. Impedance of surface footings on layered ground[J]. Computers and Structures, 2008, 86: 72-87.
[34]	TABATABAIE, M. The Flexible Volume Method for Dynamic Soil-Structure Interaction Analysis[D]. Berkeley: University of California.
[35]	RUGE P, TRINKS C, WITTE S. Time-domain analysis of unbounded media using mixed-variable formulations[J]. Earthquake Engineering & Structural Dynamics, 2001, 30(6): 899-925.
[36]	钟万勰. 结构动力方程的精细时程积分法[J]. 大连理工大学学报, 1994, 34(2): 131-136. https://www.cnki.com.cn/Article/CJFDTOTAL-DLLG199402003.htm ZHONG Wan-xie. On precise time-integration method for structural dynamics[J]. Journal of Dlian University of Technology, 1994, 34(2): 131-136. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DLLG199402003.htm

施引文献(9)

期刊类型引用(4)

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