引言

土-结构相互作用是地震工程领域的重要问题，以往对该问题的研究集中在地基土为干土的情形，对地基土为饱和土情形的研究相对较少.而实际情形中，地下水位以上可视为干土，水位以下可视为饱和土，对该情形的研究更是鲜有报导.地震波输入下土-结构相互作用分析一般采用子结构法^[1-4]、集中参数法^[5-8]和直接法^[9-12].子结构法和集中参数法将土-结构相互作用问题分解为如下可独立求解的子问题：(1) 自由场问题，即结构和基础不存在时土层在地震波入射时的响应问题；(2) 基础的输入问题，即结构不存在时，无质量刚性基础在地震波入射时的响应问题；(3) 刚性基础的动力刚度或动力柔度矩阵；(4) 结构响应问题.下面仅就饱和土层自由场分析和饱和地基动力刚度研究的相关工作做一简单介绍.

自Biot^[13-15]建立多孔弹性介质波的传播理论以来，Deresiewicz等^[16-18]对于波在饱和土界面上的传播特性进行了系统的研究, 分析了平面波在自由表面、分层多孔介质界面中的反射和透射. Jocker等^[19]通过推广Thomson-Haskell传递矩阵法，以得到一种封闭形式的解析表达式，用于计算水平成层多孔介质的反射率和透射率. Rajapakse等^[20]、Degrande等^[21]、Liang等^[22]分别给出了饱和土层和半空间的精确动力刚度矩阵，采用直接刚度法可求得饱和水平层状场地的自由场响应.李伟华等^[23]采用波函数展开法，求得饱和弹性均匀半空间的自由场，进一步分析了饱和沉积谷场地对SV波的散射问题.赵宇昕等^[24]采用传递矩阵法分析了成层饱和土层自由场响应，并探讨了数值结果中的非因果原因.

1986年，Halpern等^[25]通过对基础所覆盖的单元域上Green函数数值积分给出了三维饱和弹性半空间表面上刚性板的动力柔度系数. Bougacha等^[26]采用基础空间离散振动模态的有限元方法得到了饱和弹性土层上长条形基础和圆形基础的动力刚度. Japon等^[27]应用边界元方法分析了长条形基础的动力刚度，并考虑了渗透力、附加密度、地层深度、基础刚度对基础振动的影响.陈少林等^[28-29]提出了一种时域求解基础阻抗函数的方法，通过给基础输入脉冲位移，应用集中质量显示有限元方法结合局部透射人工边界，得到地基施加于基础的反力时程，然后根据阻抗函数定义，应用傅里叶变换求得基础阻抗函数.

子结构法以频域内的动力刚度为基础，集中参数法通过拟合动力刚度得到参数模型，理论上都只适用于线性系统，不能考虑土体的非线性.直接法将地基、基础和上部结构一并计算，可以方便考虑土体的非均匀性和非线性，但计算量一般较大，尤其是饱和土情形.梁建文等^[30-31]将单相土-结构相互作用的研究拓展到饱和土领域^[32]，采用间接边界元法分析了二维线性饱和土-结构相互作用问题，分析了P-SV波斜入射时，孔隙率、土层深度、土层与基岩间的刚度比等参数的影响.本文拟将显-隐式结合的土-结构相互作用分析方法^[10]，推广到三维饱和土-结构动力相互作用情形，采用集中质量显示有限元结合透射人工边界的方法分析饱和土体^[33-35]，采用隐式方法对上部结构进行分析，实现三维饱和土-基础-结构动力相互作用分析.

根据干土是饱和土的特殊情形这一事实，将两者情形的土-结构相互作用分析统一到饱和土-结构相互作用分析框架，进而可分析考虑地下水位情形的土-结构相互作用分析.

1 基本理论 1.1 土体内节点运动

根据Biot^[1]理论，饱和弹性多孔介质运动方程可以表示为

$ \left. \begin{array}{l} N{\nabla ^2}\mathit{\boldsymbol{u}} + \nabla \left[ {(D + N)\nabla \cdot \mathit{\boldsymbol{u}} + Q\nabla \cdot \mathit{\boldsymbol{U}}} \right] = \\ \;\;\;\;\;\;\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {{\rho _{11}}\mathit{\boldsymbol{u}} + {\rho _{12}}\mathit{\boldsymbol{U}}} \right) + \eta \frac{\partial }{{\partial t}}\left( {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{U}}} \right)\\ \nabla \left( {Q\nabla \cdot \mathit{\boldsymbol{u}} + R\nabla \cdot \mathit{\boldsymbol{U}}} \right) = \\ \;\;\;\;\;\;\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {{\rho _{12}}\mathit{\boldsymbol{u}} + {\rho _{22}}\mathit{\boldsymbol{U}}} \right) - \eta \frac{\partial }{{\partial t}}\left( {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{U}}} \right) \end{array} \right\} $

(1)

其中$\boldsymbol{u}$和$\boldsymbol{U}$为固液相位移矢量，D，N，Q，R为非负弹性常数，$\rho _{11} = \rho _1 + \rho _\alpha $，$\rho _{22} = \rho _2 + \rho _\alpha $，$\rho _{12} = - \rho _\alpha $. $\rho _1 = (1 - \phi )\rho _{\rm s} $，$\rho _2 = \phi \rho _{\rm w} $，$\rho_{\rm s}$为固相质量密度，$\rho_{\rm w}$为液相质量密度，$\rho _{\alpha }$为固、液两相惯性耦合质量密度，$\phi $为孔隙率.衰减系数$\eta = {\overline \mu \phi ^2} /{k_0 }$，$\overline \mu $为流体动力黏度，k₀为渗透率.

从地基土中取一有限土体，采用八节点三维实体元对其进行空间离散.将土体划分成两个区域，即人工边界区和内部计算区(如图 1)，将土体单元节点划分为人工边界点、与基础相连的点和内节点三类.对于内点采用集中质量显式有限单元法.在人工边界上，采用多次透射公式模拟无限域的影响.对于内节点而言，对方程(1)进行有限元空间离散，得到第i个节点的固相运动方程为^[33]

$ \boldsymbol{\ddot u}_i \boldsymbol{M}_{si} + \sum\limits_{e = 1}^N \sum\limits_{j = 1}^J \left( {\boldsymbol{C}_{ssk(i)j}^e \boldsymbol{\dot u}_j^e - \boldsymbol{C}_{swk(i)j}^e \boldsymbol{\dot U}_j^e + }\right.\\ \qquad \left.{ \boldsymbol{K}_{ssk(i)j}^e \boldsymbol{u}_j^e + \boldsymbol{K}_{swk(i)j}^e \boldsymbol{u}_j^e } \right) = \sum\limits_{e = 1}^N \boldsymbol{F}_{si}^e $

(2)

液相运动方程为

$ \begin{array}{l} {{\mathit{\boldsymbol{\ddot U}}}_i}{\mathit{\boldsymbol{M}}_{wi}} + \sum\limits_{e = 1}^N {\sum\limits_{j = 1}^J {\left( { - \mathit{\boldsymbol{C}}_{wsk(i)j}^e\mathit{\boldsymbol{\dot u}}_j^e + \mathit{\boldsymbol{C}}_{wwk(i)j}^e\mathit{\boldsymbol{\dot U}}_j^e + } \right.} } \\ \left. {\mathit{\boldsymbol{K}}_{wsk(i)j}^e\mathit{\boldsymbol{u}}_j^e + \mathit{\boldsymbol{K}}_{wwk(i)j}^e\mathit{\boldsymbol{U}}_j^e} \right) = \sum\limits_{e = 1}^N {\mathit{\boldsymbol{F}}_{wi}^e} \end{array} $

(3)

其中，$\boldsymbol{\ddot u}_i $和$\boldsymbol{\ddot U}_i $为节点i(全局编号)的固相加速度矢量和液相加速度矢量，$\boldsymbol{ u}_j^e $和$ \boldsymbol{\ddot u}_j^e $为单元e中节点j (单元局部编号)的固相位移和速度矢量，$\boldsymbol{U}_j^e $和$ \boldsymbol{\ddot U}_j^e $为单元e中节点j的液相位移和速度矢量. $\boldsymbol{M}_{si}$和$\boldsymbol{ M}_{wi}$分别表示集中在i节点上的固相和液相质量阵；$\boldsymbol{C}_{ssk(i)j}^e $, $\boldsymbol{C}_{wwk(i)j}^e $, $\boldsymbol{ C}_{wsk(i)j}^e $, $\boldsymbol{C}_{swk(i)j}^e $是考虑固液两相互相运动所产生的单元阻尼阵的子矩阵，下标$k(i)$表示节点$i $(全局编号)在单元中对应的局部节点编号为k，即全局节点编号和单元内的局部节点编号之间的对应关系. $\boldsymbol{ K}_{ssk(i)j}^e $, $\boldsymbol{K}_{wwk(i)j}^e $, $\boldsymbol{K}_{swk(i)j}^e $, $\boldsymbol{K}_{wsk(i)j}^e $表示单元刚度阵的子矩阵，$\boldsymbol{F}_{si}^e $和$\boldsymbol{F}_{wi}^e $为单元e中分配给节点$i $的固相载荷矢量和液相载荷矢量. N为包含节点i的单元个数，J为单元e的节点个数.

此外，我们可以注意到，上述阻尼矩阵仅考虑固液两相互相运动所产生的阻尼，对固相骨架考虑瑞雷阻尼

$ \boldsymbol{C}_{Rss}^e = \alpha \boldsymbol{M}_{sk(i)}^e + \beta \boldsymbol{K}_{ss}^e $

(4)

那么式(2)可修正为

$ \boldsymbol{\ddot u}_i \boldsymbol{M}_{si} + \sum\limits_e \sum\limits_{j = 1}^J \left( {\boldsymbol{C}_{ssk(i)j}^e \boldsymbol{\ddot u}_j^e - \boldsymbol{C}_{swk(i)j}^e \boldsymbol{\ddot U}_j^e + \boldsymbol{ C}_{Rssk(i)j}^e \boldsymbol{\ddot u}_j^e } \right. + \\ \qquad \left.{ \boldsymbol{K}_{ssk(i)j}^e \boldsymbol{u}_j^e + \boldsymbol{K}_{swk(i)j}^e \boldsymbol{U}_j^e } \right) = \sum\limits_e \boldsymbol{F}_{si}^e $

(5)

对内节点方程(5)和方程(3)采用中心差分与单边结合的积分格式进行时域积分.假定在p+1及其以前时刻节点i的位移矢量已知，时间离散步长用$\Delta t$表示，那么p时刻的加速度和速度可以分别表示为

$ \boldsymbol{\ddot W}^p = \dfrac{\boldsymbol{W}^{p + 1} - 2\boldsymbol{W}^p + \boldsymbol{W}^{p - 1}}{\Delta t^2} , \quad \boldsymbol{\ddot W}^p = \dfrac{\boldsymbol{W}^p - \boldsymbol{W}^{p - 1}}{\Delta t} $

(6)

其中$\boldsymbol{W}=\boldsymbol{u}$, $\boldsymbol{U}$.上式代入式(5)和式(3)可得节点i的位移递推公式

$ \begin{array}{l} \mathit{\boldsymbol{u}}_i^{p + 1} = 2\mathit{\boldsymbol{u}}_i^p - \mathit{\boldsymbol{u}}_i^{p - 1} + \Delta {t^2}\mathit{\boldsymbol{M}}_{si}^{ - 1}\left\{ {\sum\limits_e {\mathit{\boldsymbol{F}}_{si}^e} - } \right.\qquad \\ \sum\limits_e {\sum\limits_{j = 1}^J {\left[ {\left( {\mathit{\boldsymbol{C}}_{ssk(i)j}^e + \mathit{\boldsymbol{C}}_{Rssk(i)j}^e} \right)\frac{{\mathit{\boldsymbol{u}}_j^{ep} - \mathit{\boldsymbol{u}}_j^{e\left( {p - 1} \right)}}}{{\Delta t}} - } \right.} } \qquad \\ \left. {\left. {\mathit{\boldsymbol{C}}_{swk(i)j}^e\frac{{\mathit{\boldsymbol{U}}_j^{ep} - \mathit{\boldsymbol{U}}_j^{e\left( {p - 1} \right)}}}{{\Delta t}} + \mathit{\boldsymbol{K}}_{ssk(i)j}^e\mathit{\boldsymbol{u}}_j^e + \mathit{\boldsymbol{K}}_{swk(i)j}^e\mathit{\boldsymbol{U}}_j^e} \right]} \right\} \end{array} $

(7)

$ \begin{array}{l} \mathit{\boldsymbol{U}}_i^{p + 1} = 2\mathit{\boldsymbol{U}}_i^p - \mathit{\boldsymbol{U}}_i^{p - 1} + \Delta {t^2}\mathit{\boldsymbol{M}}_{wi}^{ - 1}\left[ {\sum\limits_e {\mathit{\boldsymbol{F}}_{wi}^e} } \right. - \\ \sum\limits_e {\sum\limits_{j = 1}^J {\left( { - \mathit{\boldsymbol{C}}_{wsk(i)j}^e\frac{{\mathit{\boldsymbol{u}}_j^{ep} - \mathit{\boldsymbol{u}}_j^{e\left( {p - 1} \right)}}}{{\Delta t}} + } \right.} } \\ \left. {\left. {\mathit{\boldsymbol{C}}_{wwk(i)j}^e\frac{{\mathit{\boldsymbol{U}}_j^{ep} - \mathit{\boldsymbol{U}}_j^{e\left( {p - 1} \right)}}}{{\Delta t}} + \mathit{\boldsymbol{K}}_{wsk(i)j}^e\mathit{\boldsymbol{u}}_j^e + \mathit{\boldsymbol{K}}_{wwk(i)j}^e\mathit{\boldsymbol{U}}_j^e} \right)} \right\} \end{array} $

(8)

1.2 土体边界节点运动

对于人工边界上的节点，其运动状态可通过多次透射公式^[36]进行确定

$ u_0^{p + 1} = \sum\limits_{n = 1}^N {\left( { - 1} \right)^{n + 1}C_n^N \boldsymbol{T}_n \boldsymbol{u}_n^{p + 1 - n} } $

(9)

其中，$u_0^{p + 1} $为边界节点0在p+1时刻的外传波位移(散射位移)，N为透射阶数，二项式$C_n^N $可表达为

$ C_n^N = \dfrac{N!}{\left( {N - n} \right)!n!} $

(10)

$\boldsymbol{u}_n^{p + 1 - n} $为由人工边界节点及内部节点的散射位移构成的列向量

$ \boldsymbol{u}_n^{p + 1 - n} = \left( {u_0^{p + 1 - n}, u_1^{p + 1 - n}, \cdots, u_{2n}^{p + 1 - n} } \right)^{\rm T} $

(11)

行向量$\boldsymbol{T}_n $为

$ \boldsymbol{T}_n = \left( {t_{n, 1}, t_{n, 2}, \cdots, t_{n, 2n + 1} } \right) $

(12)

其中

$ t_{n, m} = \sum {t_{1, k_1 } } t_{1, k_2 } \cdots t_{1, k_j } , \ \ m = 1, 2, \cdots, 2n + 1 $

(13)

上式中记号"$\Sigma $"表示对所有满足以下条件的项求和

$ k_1 + k_2 + \cdots + k_j = m + n-1 , \ \ k_1, k_2, \cdots, k_j = 1, 2, 3 $

(14)

以一阶Modulation Transfer Function (MTF)为例，那么可知

$ u_0^{p + 1} = \boldsymbol{T}_1 \boldsymbol{u}_1^p = t_{1, 1} u_0^p + t_{1, 2} u_1^p + t_{1, 3} u_2^p $

(15)

$ \left. {\begin{array}{*{20}{l}} \begin{array}{l} {t_{1,1}} = \left( {2 - S} \right)\left( {1 - S} \right)/2\\ {t_{1,2}} = S\left( {2 - S} \right)\\ {t_{1,3}} = S\left( {S - 1} \right)/2 \end{array} \end{array}} \right\} $

(16)

其中${S = c_{\rm a} \Delta t} / {\Delta x}$.

式(9)分别用于饱和多孔介质的固相位移和液相位移，可得到边界点在p+1时刻的散射波场位移，加上边界点的自由场位移，即可得到边界点在p+1时刻的总位移.式(11)中的散射位移向量可通过式(7)和式(8)得到的总位移减去自由场位移得到.对于饱和成层土体的自由场计算可参见文献[24].对于上覆干土，下卧饱和土的情形，可通过满足干土和饱和土分界面上的4个应力位移连续条件，得到相应的传递矩阵.

具体计算流程可参见图 2.

1.3 上部结构的运动

对于一般情况下的结构或系统，其动力学方程均可写成

$ \boldsymbol{M}\boldsymbol{\ddot u} + \boldsymbol{c}\boldsymbol{\ddot u} + \boldsymbol{k}\boldsymbol{u} = \boldsymbol{ F} $

(17)

在$p\sim p+1$时步的间隔$\Delta t$内，Newmark积分格式对位移、速度、加速度采用如下的假设关系

$ \left. \begin{array}{l} {\mathit{\boldsymbol{u}}^{p + 1}} = {\mathit{\boldsymbol{u}}^p} + {{\mathit{\boldsymbol{\dot u}}}^p}\Delta t + \left[ {\left( {\frac{1}{2} - \beta } \right){{\mathit{\boldsymbol{\ddot u}}}^p} + \beta {{\mathit{\boldsymbol{\ddot u}}}^{p + 1}}} \right]\Delta {t^2}\\ {{\mathit{\boldsymbol{\dot u}}}^{p + 1}} = {{\mathit{\boldsymbol{\dot u}}}^p} + \left[ {\left( {1 - \gamma } \right){{\mathit{\boldsymbol{\ddot u}}}^p} + \gamma {{\mathit{\boldsymbol{\ddot u}}}^{p + 1}}} \right]\Delta t\\ {{\mathit{\boldsymbol{\ddot u}}}^{p + 1}} = \frac{1}{{\beta \Delta {t^2}}}\left( {{\mathit{\boldsymbol{u}}^{p + 1}} - {\mathit{\boldsymbol{u}}^p}} \right) - \frac{1}{{\beta \Delta t}}{{\mathit{\boldsymbol{\dot u}}}^p} - \left( {\frac{1}{{2\beta }} - 1} \right){{\mathit{\boldsymbol{\ddot u}}}^p} \end{array} \right\} $

(18)

式中，$\beta $和$\gamma $为积分参数.因此将式(18)代入式(17)可以得到系统p+1时刻的运动平衡方程

$ \begin{array}{l} \left( {\mathit{\boldsymbol{K}} + \frac{1}{{\beta \Delta {t^2}}}\mathit{\boldsymbol{M}} + \frac{\gamma }{{\beta \Delta t}}\mathit{\boldsymbol{C}}} \right){\mathit{\boldsymbol{u}}^{p + 1}} = {\mathit{\boldsymbol{F}}^{p + 1}} + \mathit{\boldsymbol{M}}\left[ {\frac{1}{{\beta \Delta {t^2}}}{\mathit{\boldsymbol{u}}^p} + } \right.\\ \left. {\frac{1}{{\beta \Delta t}}{{\mathit{\boldsymbol{\dot u}}}^p} + \left( {\frac{1}{{2\beta }} - 1} \right){{\mathit{\boldsymbol{\ddot u}}}^p}} \right] + \mathit{\boldsymbol{C}}\left[ {\frac{\gamma }{{\beta \Delta t}}{\mathit{\boldsymbol{u}}^p} + } \right.\\ \left. {\left( {\frac{\gamma }{\beta } - 1} \right){{\mathit{\boldsymbol{\dot u}}}^p} + \left( {\frac{\gamma }{{2\beta }} - 1} \right)\Delta t{{\mathit{\boldsymbol{\ddot u}}}^p}} \right] \end{array} $

(19)

由式(19)可以得到$\boldsymbol{u}^{p + 1}$，进而通过式(18)可得到p+1时刻速度和加速度.

1.4 基础的运动

刚性基础作为土体和上部结构的连接部分，起到了承上启下的作用，传递结构与土体之间的作用力.其运动可由6个分量描述，即3个平动分量和3个转动分量土和结构作用在基础上的合力使基础产生刚体运动

$ \boldsymbol{M}_{\rm F} \boldsymbol{\ddot u}_{\rm F} = \boldsymbol{F} = \boldsymbol{F}_{\rm S} + \boldsymbol{F}_{\rm D} $

(20)

式中, $\boldsymbol{M}_{\rm F}$为基础的集中质量阵，对角线元素依次为三方向质量$M_{Fx}$, $M_{Fy}$和$M_{Fz}$，三方向绕质心的转动惯量$I_{Fx}$, $I_{Fy}$和$I_{Fz}$；$\boldsymbol{F}_{\rm S}$和$\boldsymbol{F}_{\rm D}$分别为结构和土体作用于基础上的广义力矢量.

首先考虑p时刻土体对刚性基础作用力$\boldsymbol{F}_{\rm D}^p $.饱和土体与基础接触任一节点k处，没有力矩的传递($M_{kx}^p = M_{ky}^p = M_{kz}^p = 0)$，只有3个方向集中力的传递，这三方向集中力分别由每一时刻固液两相本构力合成

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{c}} {F_{kx}^p}\\ {F_{ky}^p}\\ {F_{kz}^p} \end{array}} \right\} = \frac{1}{{\Delta t}}\sum\limits_L^N {\left[ {\left( {{C_{ssL}} + {C_{RssL}}} \right)\left( {u_k^p - u_k^{p - 1}} \right) - } \right.} \\ \left. {\;\;\;\;\;\;{C_{swL}}\left( {U_k^p - U_k^{p - 1}} \right)} \right] + \sum\limits_L^N {\left( {{K_{ssL}}u_k^p + {K_{swL}}U_k^p} \right)} + \\ \frac{1}{{\Delta t}}\sum\limits_L^N {\left[ { - {C_{wsL}}\left( {u_k^p - u_k^{p - 1}} \right) + {C_{wwL}}\left( {U_k^p - U_k^{p - 1}} \right)} \right]} + \\ \sum\limits_L^N {\left( {{K_{wsL}}u_k^p + {K_{wwL}}U_k^p} \right)} \end{array} $

(21)

因此饱和土体对刚性基础的合力可以表示为

$ \boldsymbol{F}_{\rm D}^p = \sum\limits_{k = 1}^m \boldsymbol{A }_k^{\rm T} \boldsymbol{F }_k^p $

(22)

$ \boldsymbol{F }_k^p = \left( {F_{kx}^p, F_{ky}^p, F_{kz}^p, 0, 0, 0} \right)^{\rm T} $

(23)

$ \boldsymbol{F }_D^p = \left( {F_{{\rm D}x}^p, F_{{\rm D}y}^p, F_{{\rm D}z}^p, M_{{\rm D}x}^p, M_{{\rm D}y}^p, M_{{\rm D}z}^p } \right)^{\rm T} $

(24)

$ {\mathit{\boldsymbol{A}}_k} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&{\Delta {Z_k}}&{ - \Delta {Y_k}}\\ 0&1&0&{ - \Delta {Z_k}}&0&{\Delta {X_k}}\\ 0&0&1&{\Delta {Y_k}}&{ - \Delta {X_k}}&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1 \end{array}} \right] $

(25)

其中，m表示土体与基础接触点总数；$\boldsymbol{F}_k^p $表示节点k对基础作用力；$F_{{\rm D}x}^p $, $F_{{\rm D}y}^p $和$F_{{\rm D}z}^p $是土体对基础三方向的集中作用合力；$M_{{\rm D}x}^p $, $M_{{\rm D}y}^p $和$M_{{\rm D}z}^p $是土体对基础三方向的转动合力矩；$\Delta X_{k}$, $\Delta Y_{k}$和$\Delta Z_{k}$是节点k对基础质心的相对坐标；矩阵$\boldsymbol{A}$是一几何转换矩阵，代表土体接触点或结构接触点和基础质心的相对坐标关系.

其次考虑p时刻结构对刚性基础作用力$\boldsymbol{F}_{\rm S}^p $.类似于上述的方法，对于任一结构与基础接触节点i处，上部结构对基础作用力可表示为

$ \boldsymbol{F }_i^p = \left( {F_{ix}^p, F_{iy}^p, F_{iz}^p, M_{ix}^p, M_{iy}^p , M_{iz}^p } \right)^{\rm T} $

(26)

可通过ANSYS分析得到各点反力值$\boldsymbol{F }_i^p $.因此上部结构对刚性基础的合力可以表示为

$ \boldsymbol{F }_{\rm S}^p = \sum\limits_{i = 1}^n \boldsymbol{A }_i^{\rm T} \boldsymbol{F }_i^p $

(27)

$ \boldsymbol{F }_{\rm S}^p = \left( {F_{{\rm S}x}^p, F_{{\rm S}y}^p, F_{{\rm S}z}^p, M_{{\rm S}x}^p, M_{{\rm S}y}^p, M_{{\rm S}z}^p } \right)^{\rm T} $

(28)

其中，n表示结构与基础接触点总数；$F_{{\rm S}x}^p $, $F_{{\rm S}y}^p $和$F_{{\rm S}z}^p $是结构对基础三方向的集中作用合力；$M_{{\rm S}x}^p $, $M_{{\rm S}y}^p $和$M_{{\rm S}z}^p $是结构对基础三方向的转动合力矩.

得到p时刻作用力$\boldsymbol{F}_{\rm D}$与$\boldsymbol{F}_{\rm S}$后，对式(20)采用中心差分格式积分，可得

$ \boldsymbol{u}_{\rm F}^{p + 1} = 2\boldsymbol{u}_{\rm F}^p - \boldsymbol{u}_{\rm F}^{p - 1} + \Delta t^2\boldsymbol{M }_{\rm F}^{ - 1} \left( {\boldsymbol{F }_{\rm D}^p + \boldsymbol{ F }_{\rm S}^p } \right) $

(29)

通过刚性基础位移可得到与其接触的土体点或结构点的位移

$ \boldsymbol{u}^{p + 1} = \boldsymbol{A \pmb u}_{\rm F}^{p + 1} $

(30)

1.5 饱和土-基础-结构相互作用分析流程

上面对饱和土-基础-结构相互作用分析的各个部分计算方法做了详细阐述，下面简述其基本分析步骤.假定p及以前时刻土体、基础与结构位移已知，计算p+1时刻各点的位移，基本步骤可以可描述如下：

(1) 采用集中质量显式有限元的方法，由式(7)和式(8)得到饱和土体内节点固、液相位移；

(2) 利用多次透射式(9)和自由场位移得到饱和土体人工边界节点的固、液相位移；

(3) 根据式(29)得到基础整体位移，进而由式(30)得到与刚性基础相连土体节点位移以及结构节点位移；

(4) 得到与基础相联的结构节点位移后，将该位移作为结构约束，根据式(19)可得到结构各点位移, 进一步可得到该时刻结构对基础的作用力；

(5) 重复执行步骤(1)~(4)，即可得到饱和土-基础-结构体系在各时刻的动力反应.

根据上述原理，通过编制相应的FORTRAN程序实现了地震波入射时饱和土-基础-结构相互作用分析.其中结构可通过ANSYS有限元软件进行分析，通过FORTRAN调用ANSYS实现饱和土-基础-结构系统的整体分析.

2 饱和土退化为干土的实现方法及验证

当流体体积模量$K_{f}$以及孔隙率$\phi $均置为零时，饱和土退化为干土^[29].下面采用如下模型对上述退化进行数值验证.

上部结构通过ANSYS进行分析，采用悬臂梁模型，截面尺寸0.1 m×0.1 m，梁高2 m，密度2 500 kg/m³，弹性模量3.00 GPa，泊松比0.2，阻尼比0.05，剪切波速707 m/s，一阶振型周期0.226 s.刚性基础尺寸6 m×6 m×4 m，密度2 500 kg/m³.考虑表 1中前两类地基土，输入相同的单位脉冲SV波：时间步距dt=0.000 2 s，脉冲宽度t₀=0.15 s，步数n=10¹³ (见图 3).土体划分单元尺寸均为1 m×1 m×1 m，满足波动模拟精度要求.结构隐式分析的时间步距dT=0.002 s，即FORTRAN每10步对ANSYS执行一次调用.

图 4和图 5分别给出了基础和结构顶部的位移.实线为干土情形的结果，虚线为饱和土退化情形的结果，两者完全吻合.因此，地基土为干土情形的土-基础-结构相互作用分析可纳入到饱和土-基础-结构相互作用分析框架，进而可以分析考虑地下水位的干土覆盖饱和土情形的土-结构相互作用分析.另外，结构在基础有效输入为零后(因为结构柔且小，不足以带动基础做相应运动)做有阻尼的自由振动，其周期正好为结构一阶振型周期0.226 s.

3 不同地基土类型对土-结构相互作用的影响

实际工程中，地下水位位于地表以下一定深度，若直接考虑为干土或饱和土并不合适.以下就地震波入射时土-基础-结构动力相互作用问题，比较地基土分别为干土情形、饱和土情形、干土覆盖饱和土情形对基础和结构响应的影响.仍采用悬臂梁结构及刚性基础，首先考虑地基土为表 1中的3号、4号、5号的情形，分别称为CASE1, CASE2, CASE3.其中，3号为饱和土，4号为具有与3号饱和土相同剪切波速和密度的干土，在工程中，常常对饱和土进行如此简化；5号为干土，具有与3号饱和土相同的骨架剪切波速和固相密度.最后考虑由上覆4 m厚的5号干土，下卧14 m厚的3号饱和土所组成的地基土，这种场地描述了当地下水位为4 m时的情况，称其为CASE4.

图 6(a)和图 6(b)分别为以上4种地基土情形时基础的位移时程和频谱. 图 7(a)和图 7(b)为结构顶点的位移时程和频谱，为了比较，图中同时给出了不考虑土结相互作用情形时的结果.由图可知，就本文算例而言，若通过两种干土模型(CASE2, CASE3)代替饱和土模型(CASE1)进行分析，基础与结构的位移均将偏小.由于CASE3中的干土剪切波速更大，其到时也比CASE1与CASE2的情况要早. CASE4与前3种情况相比，基础与结构的位移都有较大幅度减小，并不是介于干土和饱和土情形的结果之间，表明一定深度的地下水位对基础和结构响应的影响较大.不考虑土结相互作用时，脉冲位移从结构底部输入，到时较其他情形早(见图 7(a))，最大位移与饱和土-结构相互作用情形CASE1相当，但要大于其他情形(CASE2~CASE4).自由振动阶段的位移幅值小于CASE1，大于其余情形.另外，由图 7可知，由于结构柔且小，在脉冲输入结束后，结构独自做有阻尼的自由振动(此时基础反应很小)，此时土-结构相互作用影响较小，与土体情形基本无关，其周期均为结构一阶振型周期0.226 s.

4 结论

土体中的孔隙水将影响土体的地震反应，进而影响支撑其上的基础和结构的响应，为此，本文发展了一种饱和土-基础-结构动力相互作用分析的高效方法，主要结论如下：

(1) 采用集中质量显式有限元方法结合多次透射人工边界模拟半无限饱和土体，利用隐式积分方法分析上部结构，实现了地震波输入时饱和土-基础-结构动力相互作用的高效分析.

(2) 通过将流体体积模量$K_{\rm f}$，孔隙率$\phi $置为零，将干土统一到饱和土计算框架，通过算例验证了这一方案，并实现了地基为干土覆盖饱和土情形的土-结构相互作用分析.

(3) 对比了饱和土地基、干土地基、干土覆盖饱和土地基(考虑地下水位)情形时，土-结构相互作用对基础和结构响应的影响，算例表明，地下水位对基础和结构响应的影响较大.