0.   引言

近年来，随着经济全球化的进一步发展，越来越多的海上桩基础结构在建设。此外，中国处于欧亚地震带和环太平洋地震带之间，这些海上结构不同于陆地上的桩支撑结构^[1-2]，桩基础支撑的海上结构往往需要考虑地震作用下水-桩-土3者之间耦合作用对结构的影响^[3-5]，如跨海大桥、海上风机、石油平台和高桩码头等海上结构。因此，评估考虑水-桩-土相互作用的海上桩基础结构在地震作用下的动力响应对海上桩基础结构设计有一定的参考价值。

目前，关于桩-土相互作用的耦合振动解析研究方法大概可以分为3类：①第一种是Winkler模型方法^[6-8]，Winkler模型简单易用，促进了桩基理论水平的研究，但也受限于一些模型参数的取值；②第二种是平面应变方法^[9-10]，这种方法忽略了周围土体的三维波效应；③第三种是考虑波散射效应的三维轴对称连续性模型方法^[11-12]。在上述研究中，土体底部和桩端土都是基于地基刚性假定，这种简化方法一般适用于端承桩，对于浮动桩来说存在一定程度上的误差。

不同于土体底部或桩端土的刚性边界条件，后续又建立了支撑土体或桩端土的黏弹性边界条件（由线性弹簧和阻尼器并联而成）。例如，Wang等^[7]给出了嵌在层状土体中的黏弹性桩体在轴向荷载作用下动力响应的解析解；Deng等^[13]和Lü等^[14]通过引入虚拟桩-土模型，分别研究了在单层土体和多层土体中桩-土相互作用对桩体竖向振动响应的影响；吕述晖等^[15]推导了在土体弹性半空间情况下管桩在竖向动力荷载作用下动力响应的解析解，同时也给出了桩上土体剪切复合刚度，对管桩的工程应用具有一定的参考价值。通过上述研究发现，土体弹性半空间情况下的弹簧和阻尼值对桩体动力响应有一定的影响。同时也发现针对土体弹性半空间条件下桩-土的解析研究都是基于在桩顶施加竖向荷载的情况，对地震作用下考虑土体弹性半空间的桩-土相互作用系统对桩体动力响应的影响研究目前还没有。

海上桩基础结构一般都是处于深水环境中，因此在动力分析时往往需要考虑桩-水相互作用对桩体动力响应的影响。桩-水相互作用不仅对结构的动力特性有一定的影响，而且水体还会对结构产生额外的动水压力作用。基于刚性地基，各学者对桩-水动力相互作用进行了一系列的研究^[16-18]，如Wei等^[19]建立了一种等截面轴对称柱体与水体相互作用的简化抗震设计分析方法；Jiang等^[20]通过精细的数据拟合，推导出了圆柱体动水压力的简化公式，其主要参数为柱体半径和高度。

以上研究分析发现桩-水相互作用大多基于刚性基础，考虑地基特性影响的研究较少。而关于桩-土相互作用的研究，对于近海结构来说，往往需要考虑水体的影响，从而建立一个系统模型。目前，关于水-桩-土相互作用的数值研究较为丰富^[21-23]。从上述水-桩-土相互作用的数值研究中可以发现，几乎所有的假设都是基于刚性基础，很少考虑地基为弹性半空间的情况，特别是针对水平动荷载的研究更是少见。虽然数值模拟可以更方便、有效地解决水-桩-土相互作用问题，但对于某些参数规律的研究和讨论，解析解有时较为方便。此外，根据笔者目前的认识，考虑土体弹性半空间的水-桩-土相互作用系统在地震作用下的解析研究仍是一个亟待解决的问题。为此，本文建立了考虑土体弹性半空间的水-桩-土相互作用在地震作用下解析解，重点分析了土体弹性半空间和水-土参数对桩体位移响应的影响。

1.   基本假定及方程的建立

考虑土体弹性半空间的水-桩-土相互作用系统模型如图 1所示。模型底部受到SV波${u_{\text{g}}}{{\text{e}}^{{\text{i}}\omega t}}$作用（其中$\omega$，$t$分别为圆频率和时间），(r, z)为桩体柱坐标系，桩体半径、桩体出水面高度、水深、土体厚度和桩体总高度分别为r₀，h₁，h₂，h₃，H。

图  1  水-桩-土相互作用系统模型

Figure  1.  Model for water-pile-soil interaction system

下载: 全尺寸图片 幻灯片

本文在建立方程时需要采用如下假定条件：①水体的水平方向是无限域，且是无旋、无黏性，初始时刻不发生扰动；②桩体采用线黏弹性Euler梁模拟，桩体截面为圆形。桩体与水体、土体和土体弹性地基接触良好，水-桩-土体系为小变形振动，接触面不发生滑移；③土体为黏弹性、均质、各向同性介质，土体材料阻尼采用与频率无关的滞回阻尼；④桩体和土体主要变形在水平振动方向，忽略与水平振动相关的竖向位移^[9-10]。

参考竖向动荷载作用下桩-土相互作用的弹簧和阻尼系数^[13-15]，并进一步结合由Gazetas^[24]提出的水平弹簧和阻尼系数，桩体和土体底部弹性半空间可以简化为弹簧-阻尼系统。此外，桩体弹簧、土体弹簧以及桩体阻尼和土体阻尼系数的表达式可表示为

式中：${G_{\text{b}}}$，${\nu _{{\text{sb}}}}$，${c_{{\text{sb}}}}$，${\rho _{\text{b}}}$，${A_{\text{p}}}$分别为土体底部弹性半空间的剪切模量、泊松比、剪切波速、密度以及桩体截面面积，${c_{{\text{sb}}}} = \sqrt {{G_{\text{b}}}/{\rho _{\text{b}}}}$。

1.1   水体控制方程

根据势流理论，水体动水压力在柱坐标下可以表示为$p(r,\theta ,z){{\text{e}}^{{\text{i}}\omega t}}$，以水体动水压力表示的控制方程为

式中：$k = \omega /{c_{\text{w}}}$；${c_{\text{w}}}$为水中波传播的速度。

当水体波速${c_{\text{w}}} \to \infty$时，$k \to 0$，此时水体表示不可压缩。

1.2   土体控制方程

土体介质在SV波作用下，可以分为自由场响应和散射场响应，其中自由场位移响应将在下一章节介绍，而散射场位移响应可以表示为$u_r^{}(r,\theta ,z){{\text{e}}^{{\text{i}}\omega t}}$和$u_\theta ^{}{\text{(}}r,\theta ,z{\text{)}}{{\text{e}}^{{\text{i}}\omega t}}$。土体径向和环向位移表示的耦合动力方程可以表示为

式中：$\lambda _{\text{s}}^ * = 2G_{\text{s}}^ * {\nu _{\text{s}}}/(1 - 2{\nu _{\text{s}}})$为土体复Lame常数；${\nu _{\text{s}}}$为土体泊松比；$G_{\text{s}}^ * = {G_{\text{s}}}(1 + 2{\text{i}}{\beta _{\text{s}}})$为土体复剪切模量；${G_{\text{s}}}$为土体剪切模量；${\beta _{\text{s}}}$为土体滞后阻尼；${\rho _{\text{s}}}$为土体质量密度。

1.3   桩体控制方程

基于一维欧拉-伯努利梁理论，桩体在空气、水体和土体中的水平位移${w_{{p_i}}}{{\text{e}}^{{\text{i}}\omega t}}(z){\text{(}}i = 1,2,3{\text{)}}$表示的动力控制方程为

式中：${E_{\text{p}}}$，${I_{\text{p}}}$，m_p分别为桩体弹性模量、截面惯性矩和单位长度的质量，${m_{\text{p}}} = {\rho _{\text{p}}}{\text{π }}r_0^2$；$E_{\text{p}}^ * = {E_{\text{p}}}(1 + 2i{\beta _{\text{p}}})$，${\beta _{\text{p}}}$为桩体滞后阻尼；${f_i}(i = 1,2,3)$分别为空气$({f_1} = 0)$、水体和土体对桩体的阻抗力。需要注意的是关于桩-土地底部弹性半空间和上部土体-土体底部弹性半空间相互作用简化为弹簧-阻尼模型时，地震动是在弹簧的末端输入，通过水平应力平衡来体现水平SV波作用下桩体和上部土体与土体弹性半空间的相互作用。

1.4   系统边界条件

桩体顶部边界条件：

水体表面边界条件：

水-土交界面边界条件：

弹性半空间表面边界条件：

式中：$u_x^{\text{s}}$为土体散射场在直角坐标系下的x方向位移，与笛卡尔坐标系下位移关系为$u_x^{\text{s}} = u_r^{\text{s}}\cos \theta - u_\theta ^{\text{s}}\sin \theta$。

桩段连接面边界条件：

桩-土交界面边界条件：

桩-水交界面边界条件：

水体动水压力无穷远处边界条件：

土体位移无穷远处边界条件：

2.   方程的求解

2.1   水-土-弹性半空间自由场相互作用

SV波垂直入射时，水-土自由场模型如图 2所示，其中直角坐标系原点在水-土交界面处。模型中弹性半空间固体介质的波系${u_{\text{b}}}(z,t)$由入射SV波${u_{{\text{bi}}}}(z,t)$和反射SV波${u_{{\text{br}}}}(z,t)$构成，弹性半空间上层土体介质由透射SV波${u_{{\text{st}}}}(z,t)$和反射SV波${u_{{\text{sr}}}}(z,t)$构成。与P波垂直入射时不同，SV波垂直入射时海域场地自由场模型中固、液介质不发生能量交换，水体介质中没有波的输入。

图  2  水-土自由场模型

Figure  2.  Model for water-soil free field

下载: 全尺寸图片 幻灯片

根据Helmholtz定理，引入位移势函数$\phi (z,t)$，则波系中各成分的位移可以表示为

对于频率为$\omega$的简谐入射波，稳态解的势函数可以表示为

式中：${c_{{\text{ss}}}} = \sqrt {{G_{\text{s}}}/{\rho _{\text{s}}}}$为弹性半空间上层土体介质的SV波波速${c_{ss}}{\text{ = }}\sqrt {{G_s}/{\rho _s}}$；${A_{{\text{bi}}}}$，${A_{{\text{br}}}}$分别为弹性半空间固体介质入射波和反射波的势函数幅值；${A_{{\text{st}}}}$，${A_{{\text{sr}}}}$分别为弹性半空间上层土体介质和反射波的势函数幅值。

在土体顶部($z = 0$)和弹性半空间与上层土体介质耦合界面处($z = - {h_3}$)，满足土体顶部切向应力为零、法向位移连续和法向应力平衡的边界条件，

将式（27）~（30）代入平衡条件式（31）~（33），进一步可以得到：

式中：

将式（34）~（36）代入式（26），并作Fourier变换，得到波场频域解：

对于任一给定的输入地震动${u_{\text{g}}}(t)$，可以是简谐波或者实际地震动。考虑到实际地震动的传播时间较长，在傅里叶变换到频域下的频谱时较为繁琐，因此本文算例考虑使用简谐波时程，同时也利于一些参数规律化研究。本文研究采用的简谐波时程是狄拉克脉冲，位移幅值为1 cm，脉冲激励时间为0.2 s，脉冲时程如图 3所示。对时域下位移时程作Fourier变换，然后由${\left. {U_{{\text{si}}}^{\text{f}}(z,\omega )} \right|_{z = - {h_3}}} = {U_{\text{g}}}(\omega )$得到：

图  3  狄拉克脉冲

Figure  3.  Dirac impulse

下载: 全尺寸图片 幻灯片

将式（40）代入式（39），得到弹性半空间固体域介质和弹性半空间上层土体介质的频域解，再作Fourier逆变换后，即可得到弹性半空间波场时域解。

在无其他特殊说明情况下采用表 1所示的水-桩-土参数，为了验证上述推导的SV波作用下水-土-弹性半空间自由场解析解的合理性，将本自由场解析解退化为只有水-土情况下，然后分别与张国良的时域一维有限元法^[25]、宝鑫的解析解法^[26]以及刚度矩阵法^[27]进行对比，结果如图 4所示。从图 4可以看出，几种方法吻合较好，也反映了本方法的合理性。

表  1  桩体、土体和水体参数

Table  1.  Parameters of pile, soil and water

参数	幅值		参数	幅值
桩长H	60 m		杨氏模量Es	50 MPa
桩出水面长度h1	0 m		密度${\rho _{\text{s}}}$	2000${\text{kg}}/{{\text{m}}^3}$
直径2r0	4 m		滞后阻尼${\beta _{\text{s}}}$	0.02
杨氏模量Ep	30 GPa		泊松比${\nu _{\text{s}}}$	0.35
剪切模量Gp	—		弹性半空间剪切模量Gb	100 MPa
密度${\rho _{\text{p}}}$	2500${\text{kg}}/{{\text{m}}^3}$		弹性半空间密度${\rho _{\text{b}}}$	2200${\text{kg}}/{{\text{m}}^3}$
滞后阻尼${\beta _{\text{p}}}$	0.02		水深h2	30 m
泊松比${\nu _{\text{p}}}$	0.25		密度${\rho _{\text{w}}}$	1000${\text{kg}}/{{\text{m}}^3}$
土体厚度h3	30 m		波速${c_{\text{w}}}$	1438 m/s

下载: 导出CSV 

| 显示表格

图  4  几种方法的位移时程对比

Figure  4.  Comparison of time histories of displacement by several methods

下载: 全尺寸图片 幻灯片

图 5，6分别表示土体顶部位移随着土体底部弹性半空间剪切模量和土体弹性模量的变化趋势。

图  5  不同土体弹性半空间剪切模量情况下土体顶部位移时程

Figure  5.  Time histories of displacement of soil top under different elastic half-space shear moduli

下载: 全尺寸图片 幻灯片

图  6  土体顶部最大位移随着土体底部弹性半空间剪切模量变化

Figure  6.  Variation of maximum displacement of soil top with elastic half-space shear modulus at soil bottom

下载: 全尺寸图片 幻灯片

从图 5，6可以看出，软弱土层的土体位移偏小以及相比于土体底部弹性地基情况，土体底部刚性地基情况下的位移偏大，也就是说土体底部地基不能简单的视为刚性地基。

2.2   桩-土相互作用

根据亥姆霍兹分解、分离变量法和一系列边界条件（12），（13），（21），（22），（24），（25）可以得到土体径向位移和环向位移为^[28]

式中：$q_i^2{\text{ = }}\frac{{G_{\text{s}}^ * g_i^2 - {\rho _{\text{s}}}{\omega ^2}}}{{\lambda _{\text{s}}^ * + 2G_{\text{s}}^ * }}$；$s_i^2 = \frac{{G_{\text{s}}^{\text{*}}g_i^2 - {\rho _{\text{s}}}{\omega ^2}}}{{G_{\text{s}}^ * }}$；${A_i}$，${B_i}$为未知系数；${z_3} = z - {h_1} - {h_2}$。

${g_i}$可以通过边界条件（14）求解，同时还需要满足在任意角度下等式都成立。因此，整理得到以下含有虚数的超越方程：

式中：${\gamma _{\text{s}}} = \frac{{{k_{{\text{sh}}}} + {\text{i}}\omega {c_{{\text{sh}}}}}}{{G_{\text{s}}^ * }}$。

根据桩-土位移连续边界条件，且桩完全附着在土体上，因此可以假设桩体位移为

根据式（39），（41），（42）可以得到：

式中：

式中：$U_{1i}^{}$为与z相关的未知系数。

再根据力和位移的关系，可以得到土体阻抗力${f_1}$：

式中：${b_{{\text{h}}i}}$为第i阶水平阻抗因子。

当$i = 3$时，将式（48）代入式（7），并利用通解和特解求得桩体在土中的位移表达式为^[28]

式中：${Q_1}$，${Q_2}$，${Q_3}$，${Q_4}$为由未知常数。

结合式（44），（50），进一步利用三角函数正交性，可以得到${U_{1i}}$的表示式。然后再带回到式（51），整理得到桩体在土中的位移表达式为

式中：

2.3   桩-水相互作用

根据分离变量法和一系列边界条件（10），（11），（23），（24），可以得到水体水平方向对桩体的阻抗力${f_2}$^[28]：

式中：${h_1} \leqslant {z_2} \leqslant {h_1} + {h_2}$；${a_{hj}}$为第j阶水平阻抗因子；${r_1} = - {\text{i}}r{\lambda _j}\sqrt {{\text{1}} - {{(\omega /{\lambda _j}c)}^2}}$，${\lambda _j} = \frac{{(2j - 1){\text{π }}}}{{2{h_2}}}$，$j{\text{ = }}1,2,3 \cdots$。

当$i = 2$时，将式（54）代入式（7），解法与桩-土相互作用求解类似，整理得到桩体在水中的位移表达式为

式中：${Q_5}$，${Q_6}$，${Q_7}$，${Q_8}$为未知常数。且

2.4   空气中桩

当$i = 1$时，桩体在空气中的位移表达式为

式中：${Q_9}$，${Q_{10}}$，${Q_{11}}$，${Q_{12}}$为由未知常数。

2.5   桩体总位移求解

将式（52），（58），（60）代入边界条件（8），（9），（15）~（20），经过转化可得

将式（62）代入式（63），并进一步代入式（64）再结合式（61）可得

未知系数${Q_{{w_1}}}$，${Q_{{w_2}}}$，${Q_{{w_3}}}$分别可以通过式（62），（63），（65）得到。基于这些系数，即可求得频域下桩体总位移，再通过Fourier逆变换，可得到时域下桩体总位移。

3.   模型验证

由式（48），（54）可以看出土体阻抗和水体阻抗的表达式是一种无穷级数形式，因此有必要分析无穷级数的收敛性。从图 7中可以看出，无水和有水两种情况下将无穷级数截止到第5项之后，计算结果趋于平稳，说明无穷级数趋于收敛。因此，考虑到计算精确和计算效率，i和j都取为5。

图  7  收敛性分析

Figure  7.  Convergence tests

下载: 全尺寸图片 幻灯片

同时，为了验证本文解析解的合理性，将本解析解退化为两种情况：①退化为无水有土情况（弹性地基），然后与子结构法进行对比^[29]，结果如图 8所示；②退化为有水有土情况（刚性地基），然后与已有的刚性解析解进行对比^[28]，结果如图 9所示，其中弹簧取大值表示在本算例中土体底部弹性半空间的剪切模量取值为${G_{\text{b}}} = 100\;{\text{GPa}}$。从图 8，9中可以看出，3种方法吻合度较好，也即说明本解析解推导的合理性。

图  8  本解析法与子结构方法的对比(h₂=0 m)

Figure  8.  Comparison between analytical and substructure methods (h₂=0 m)

下载: 全尺寸图片 幻灯片

图  9  本解析法与刚性解析解法的对比

Figure  9.  Comparison between analytical and rigid analytical methods

下载: 全尺寸图片 幻灯片

4.   数值算例

4.1   土体底部弹性半空间刚度的影响

为了进一步了解土体底部和桩端土弹簧和阻尼的取值范围对桩体位移响应的影响，在其他条件不变以及无水情况下，分析了桩体最大位移随着土体底部弹性半空间剪切模量以及土体弹性模量的变化趋势，结果如图 10所示。从图 10可以看出，桩体最大位移随着土体底部弹性半空间剪切模量的增大而增大，也就是说相比于弹性地基情况，刚性地基情况下对桩体位移响应有一定的放大作用，且土体弹性模量越大，放大作用越明显。需要注意的是图 10研究的是桩顶最大位移随着土体底部弹性半空间剪切模量和土体弹性模量的变化趋势，而图 6研究的是水-土-土体弹性半空间的土体自由场顶部最大位移随着土体底部弹性半空间剪切模量和土体弹性模量的变化趋势。

图  10  桩体最大位移随着土体底部弹性半空间剪切模量变化

Figure  10.  Variation of maximum displacement of pile with elastic half-space shear modulus at soil bottom

下载: 全尺寸图片 幻灯片

从模型简化的过程中，可以看出造成上述桩体位移响应差异的主要原因有两个：①土体底部弹簧和阻尼取值对土体模态${g_i}$和阻抗${b_{{\text{h}}i}}\cos ({g_i}{z_3})$的影响，进而影响桩体位移响应，分析结果如图 11，12所示；②桩端土弹簧和阻尼取值对桩体位移响应的影响，分析结果如图 13所示。

图  11  土体模态随着土体底部弹性半空间剪切模量的变化

Figure  11.  Variation of soil mode with elastic half-space shear modulus at soil bottom

下载: 全尺寸图片 幻灯片

图  12  土体阻抗随土体底部弹性半空间剪切模量的变化

Figure  12.  Variation of soil impedance with elastic half-space shear modulus at soil bottom

下载: 全尺寸图片 幻灯片

图  13  T₈和T₉矩阵随着土体底部弹性半空间剪切模量的变化

Figure  13.  Variation of T₈ and T₉ matrices with elastic half space shear modulus at soil bottom

下载: 全尺寸图片 幻灯片

从图 11，12可以看出，随着土体底部弹性半空间剪切模量的变化，土体模态和阻抗的变化趋势较小，说明地震作用下土体底部地基可以视为刚性地基。从图 13可以看出，随着土体底部弹性半空间剪切模量的变化（桩端土弹簧和阻尼值随之变化），对${T_{\mathbf{8}}}(2 \times 4)$和${T_9}(2 \times 1)$矩阵造成了一定程度上的影响，而这也是造成桩体顶部最大位移差异的主要原因。

4.2   水深的影响

从工程应用的角度来看，探讨水体深度变化对桩体位移响应造成的影响有一定意义。因此，定义了一个无量纲参数E_r（误差），

式中：${u}_{有水}$为不同水深情况下的桩体位移峰值；${u}_{无水}$为无水情况下的桩体位移峰值。

图 14，15表示在不同桩体半径以及土体弹性模量情况下，研究桩体响应的误差E_r与水深(桩长随着水深的变化而变化)之间的关系。从图 14，15中可以看出，不同桩体半径和土体弹性模量情况下，水深的变化都对桩体位移响应有造成一定的影响。从工程角度来看，对于软弱土层且处于低水位情况下，水深的变化对桩体位移响应造成的影响（8%左右）可以忽略不计。

图  14  不同桩体半径情况下${E_{\text{r}}}$随着水深的变化

Figure  14.  Variation of ${E_{\text{r}}}$ with water depth under different pile radii

下载: 全尺寸图片 幻灯片

图  15  不同土体弹性模量情况下${E_{\text{r}}}$随着水深的变化

Figure  15.  Variation of ${E_{\text{r}}}$ with water depth under different soil elastic moduli

下载: 全尺寸图片 幻灯片

5.   结论

本文建立了一种考虑土体弹性半空间的水-桩-土相互作用系统在SV波作用下的解析解。为了深入了解土体底部弹性地基和刚性地基对桩体位移响应造成的差异，重点讨论了桩体位移响应随着土体底部弹性半空间剪切模量的变化趋势，同时也研究了桩体、土体参数以及水深对桩体位移响应造成的影响。得出3点结论。

（1）在本文所研究的参数范围内，相比于土体底部弹性地基情况，土体底部是刚性地基情况下对土体自由场位移响应有一定的影响，而对土体散射场位移响应造成的影响可以忽略。

（2）同样，在本文所研究的参数范围内，场地是中硬土以上情况下土体底部地基可以视为刚性地基，而对桩体位移造成的影响主要取决于桩端土地基类型的选择。

（3）随着桩体和土体参数变化，水体对桩体位移响应有一定影响。然而从工程角度来看，水深小于15 m时对软弱土层中桩体位移响应造成的影响可以忽略。