近年来，环境振动问题日益突出^[1-2]，利用桩-土周期结构设计的地基隔振屏障得到了成功应用^[3-4]. 早期学者们沿用传统的周期结构理论，将土体视为理想线弹性材料，对单层^[5]或多层^[6]土体排桩周期结构的禁振带隙和振动衰减规律进行研究. 最近几年，Shi等^[7-10]借助COMSOL有限元软件中的多孔弹性波模块，研究两相饱和土中周期性排桩的带隙和振动传输规律. 简单的散射型声子晶体产生的带隙位于高频段，很难覆盖实际中地震和交通荷载引起的低频振动. Pu等^[11-12]通过填充软材料，研究具有局域共振性质的套管桩隔振结构. 陈晓斌等^[13]利用有限元软件，模拟更复杂的四组元开孔套管桩周期结构的隔振性能. 实际土体由于松弛和蠕变效应，应力作用下土颗粒的重新排列和骨架错动具有明显的时间效应^[14-15]. 蔡袁强等^[16]基于Terzaghi一维固结理论，分析循环荷载作用下土体应力-应变在时间效应下的非线性渐变过程. 曾庆有等^[17]采用Mesri蠕变模型考虑土体的流变性，计算黏弹性土中桩基的长期沉降. 艾志勇等^[18-19]采用有限元-边界元耦合的方法，研究层状黏弹性地基与桩以及梁共同作用的时变行为. 何利军等^[20]提出黏弹性土的蠕变模型，对土体应力-应变的时间效应进行精确分析. 汪磊等^[21]研究半透水边界下分数阶黏弹性饱和土的固结时效行为. 吴奎等^[22]研究流变性岩体中深埋隧道的力学响应问题. 张婉洁等^[23]基于分数阶Bingham模型的线性刚度系统，研究磁流变液阻尼器的时滞问题. 上述研究表明，流变性土的应力-应变对于时间具有依赖性.

为了探究流变性土排桩周期结构的带隙特征以及土体流变性的影响规律，本文采用时间依赖性模量描述土体流变效应，构建桩-土周期结构的连续介质动力学模型. 采用多重散射法推导周期系统的体波弥散方程，通过搜索布里渊域得出周期性排桩地基的频率域弥散曲线，利用模型退化验证算法的正确性. 分析流变性土排桩地基的禁振带隙及核心参数的影响规律.

如图1(a)所示为排桩隔振地基. 假设桩体足够长（桩长远大于元胞尺寸）且与周围土体完全接触，在小应变情形下，元胞平面上各点的位移都与截面平行，所以桩-土结构可以简化为平面应变问题^[7-10]. 针对桩基定义局部极坐标系r-θ（见图1(b)），对Bragg单元定义二维直角坐标系x-y（见图1(c)）. 图中，a为元胞常数，r₀为桩基半径，Γ₁、Γ₂、Γ₃、Γ₄为元胞的4个边界.

采用具有时间依赖性的模量描述土体的流变效应，积分形式的本构方程^[24]为

(1) $ {\sigma ^{\rm{s}}} = \int_{ - \infty }^t {\left[ {{\lambda _{\rm{s}}}\left( {t - \tau } \right){\delta _{ij}}\frac{{{\rm{d}}{\zeta _{\rm{s}}}}}{{{\rm{d}}\tau }}+2{\mu _{\rm{s}}}\left( {t - \tau } \right)\frac{{{\rm{d}}{\varepsilon ^{\rm{s}}}}}{{{\rm{d}}\tau }}} \right]} {\rm{d}}\tau . $

式中：σ^s、ε^s分别为土体的应力和应变，δ_ij为Kronecker delta函数；ζ_s为体应变；λ_s = λ_s(t)、μ_s = μ_s(t)为Lamé常数，由于考虑土体模量的时间依赖性， λ_s、μ_s均为时间t的函数；τ为松弛时间.

根据Lamé常数的定义，可知

(2) $ {\lambda _{\rm{s}}}\left( t \right) = \frac{{{E_{\rm{s}}}\left( t \right){v_{\rm{s}}}}}{{\left( {1+{v_{\rm{s}}}} \right)\left( {1 - 2{v_{\rm{s}}}} \right)}} \text{，} {\mu _{\rm{s}}}\left( t \right) = \frac{{{E_{\rm{s}}}\left( t \right)}}{{2\left( {1+{v_{\rm{s}}}} \right)}} . $

式中：v_s为泊松比（本文中假定v_s不随时间变化）；E(t)为弹性模量，松弛函数形式^[24]为

(3) $ {E_{\rm{s}}}\left( t \right) = {E_\infty }+\left( {{E_0} - {E_\infty }} \right){{\rm{exp}}} \;\left( { - t/\tau } \right). $

式中：E₀、E_∞分别为初始和最终稳定状态的弹性模量.

对式（3）引入时间t的Fourier变换，可以得到频率域的表达式：

(4) $ {E_{\rm{s}}}\left( \omega \right) = E_{\rm{s}}^{\rm{R}}\left( \omega \right)+{\rm{i}}E_{\rm{s}}^{\rm{I}}\left( \omega \right) . $

式中：ω为角频率；i为虚数单位； $E_{\rm{s}}^{\rm{R}}\left( \omega \right)$、 $E_{\rm{s}}^{\rm{I}}\left( \omega \right)$表示实部和虚部，

$ E_{\rm{s}}^{\rm{R}}\left( \omega \right) = \frac{{{E_\infty }+{E_0}{\omega ^2}{\tau ^2}}}{{1+{\omega ^2}{\tau ^2}}} \text{，} E_{\rm{s}}^{\rm{I}}\left( \omega \right) = \frac{{\left( {{E_\infty } - {E_0}} \right)\omega \tau }}{{1+{\omega ^2}{\tau ^2}}} . $

对式（1）进行时间t的Fourier变换，可以得到频率域的本构方程：

(5) $ {\sigma ^{\rm{s}}} = {\lambda _{\rm{s}}}\left( \omega \right){\delta _{ij}}{\zeta _{\rm{s}}}+2{\mu _{\rm{s}}}\left( \omega \right){\varepsilon ^{\rm{s}}} . $

式中：λ_s(ω)、μ_s(ω)为频率域中对应的Lamé常数，表达式为

(6) $ {\lambda _{\rm{s}}}\left( \omega \right) = \frac{{{\lambda _0}\left( {\eta +{\omega ^2}{\tau ^2}} \right)}}{{1+{\omega ^2}{\tau ^2}}}+{\rm{i}}\frac{{{\lambda _0}\omega \tau \left( {\eta - 1} \right)}}{{1+{\omega ^2}{\tau ^2}}} , _{ } $

(7) $ {\mu _{\rm{s}}}\left( \omega \right) = \frac{{{\mu _0}\left( {\eta +{\omega ^2}{\tau ^2}} \right)}}{{1+{\omega ^2}{\tau ^2}}}+{\rm{i}}\frac{{{\mu _0}\omega \tau \left( {\eta - 1} \right)}}{{1+{\omega ^2}{\tau ^2}}} . $

其中η为初始和最终状态的模量比，η = E_∞/E₀；λ₀、μ₀为初始Lamé常数，

$ {\lambda _0} = \frac{{{v_{\rm{s}}}{E_0}}}{{\left( {1+{v_{\rm{s}}}} \right)\left( {1 - 2{v_{\rm{s}}}} \right)}} \text{，} {\mu _0} = \frac{{{E_0}}}{{2\left( {1+{v_{\rm{s}}}} \right)}} . $

根据黏性理论，可得滞回阻尼比为

(8) $ {\xi _{\rm{d}}} = \frac{{{E^{\rm{I}}}}}{{{E^{\rm{R}}}}} = \frac{{\left( {\eta - 1} \right)\omega \tau }}{{\eta +{\omega ^2}{\tau ^2}}} . $

从式（6）~（8）可以看出，土体的变形模量λ_s、μ_s及ξ_d均为频率的函数. 图2给出不同的η和τ下，ξ_d随频率的变化曲线. 可知，阻尼比随频率发生非单调性的变化. 从图1（a）可知，当η > 1（即 E_∞ > E₀）时，阻尼比为正值，且先增大而后减小，最终趋于收敛. 当η < 1（即 E_∞ < E₀）时，阻尼比为负值，且先负向增大而后收敛. 图1(a)说明，阻尼比随着频率的增大而出现非单调性变化. 从图1(b)可知，τ影响阻尼比的峰值位置.

通过对式（8）求极值，可得阻尼比的峰值和峰值对应的频率：

(9) $ {\left( {{\xi _{\rm{d}}}} \right)_{\max }} = \frac{{\eta - 1}}{{\sqrt 2 \left( {\eta +{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)}} \text{，} {\omega _{{\rm{cr}}}} = \frac{1}{{\sqrt 2 \tau }} . $

从图2及式（9）可知，η决定阻尼比的幅值，τ决定阻尼比随频率的变化速率. 在实际工程中，土体受载后被压密，变形模量将增大（E_∞ > E₀），因此不会出现η < 1的情况，不再讨论该情况.

根据弹性波动理论可知，土体和桩基的频率域Navier波动方程为

(10) $ \left( {{\lambda _{\rm{s}}}+{\mu _{\rm{s}}}} \right)\nabla \left( {\nabla \cdot {{\boldsymbol{u}}^{\rm{s}}}} \right)+{\mu _{\rm{s}}}{\nabla ^2}{{\boldsymbol{u}}^{\rm{s}}}+{\rho _{\rm{s}}}{\omega ^2}{{\boldsymbol{u}}^{\rm{s}}} = {\boldsymbol{0}} , $

(11) $ \left( {{\lambda _{\rm{p}}}+{\mu _{\rm{p}}}} \right)\nabla \left( {\nabla \cdot {{\boldsymbol{u}}^{\rm{p}}}} \right)+{\mu _{\rm{p}}}{\nabla ^2}{{\boldsymbol{u}}^{\rm{p}}}+{\rho _{\rm{p}}}{\omega ^2}{{\boldsymbol{u}}^{\rm{p}}} = {\boldsymbol{0}} . $

式中：u^s、u^p、ρ_s、ρ_p分别为土体和桩基的位移矢量和密度；与流变性土体不同，桩基为均质弹性材料，其Lamé常数λ_p、μ_p为常数.

引入势函数 ${\phi _{\rm{s}}}$、ψ_s，对u^s在极坐标系下进行Helmholtz分解，得到

(12) $ u_r^{\rm{s}} = \frac{{\partial {\phi _{\rm{s}}}}}{{\partial r}}+\frac{1}{r}\frac{{\partial {\psi _{\rm{s}}}}}{{\partial \theta }} \text{，} u_\theta ^{\rm{s}} = \frac{1}{r}\frac{{\partial {\phi _{\rm{s}}}}}{{\partial \theta }} - \frac{{\partial {\psi _{\rm{s}}}}}{{\partial r}} . $

将式（12）代入式（10），可得

(13) $ \left. \begin{gathered} {\nabla ^2}{\phi _{\rm{s}}}+k_{{\rm{l}}1}^2{\phi _{\rm{s}}} = 0\begin{array}{*{20}{c}} , \end{array} \\ {\nabla ^2}{\psi _{\rm{s}}}+k_{{\rm{t}}1}^2{\psi _{\rm{s}}} = 0\begin{array}{*{20}{c}} . \end{array} \\ \end{gathered} \right\} $

式中：k_l1、k_t1分别为土体中压缩波和剪切波的波数，k_l1 = ω/c_l1、k_t1= ω/c_t1，其中c_l1、c_t1为2种波的波速.

将式（12）代入本构方程（5），可得

(14) $ \left. \begin{gathered} \sigma _r^{\rm{s}} = {\lambda _{\rm{s}}}\left( \omega \right){\nabla ^2}{\phi _{\rm{s}}}+2{\mu _{\rm{s}}}\left( \omega \right)\frac{\partial }{{\partial r}}\left( {\frac{{\partial {\phi _{\rm{s}}}}}{{\partial r}}+\frac{1}{r}\frac{{\partial {\psi _{\rm{s}}}}}{{\partial \theta }}} \right)\begin{array}{*{20}{c}} , \end{array} \\ \sigma _{r\theta }^{\rm{s}} = {\mu _{\rm{s}}}\left( \omega \right)\left( {\frac{2}{r}\frac{{{\partial ^2}{\phi _{\rm{s}}}}}{{\partial r\partial \theta }} - \frac{2}{{{r^2}}}\frac{{\partial {\phi _{\rm{s}}}}}{{\partial \theta }} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{\psi _{\rm{s}}}}}{{\partial {\theta ^2}}} - \frac{{{\partial ^2}{\psi _{\rm{s}}}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {\psi _{\rm{s}}}}}{{\partial r}}} \right). \end{gathered} \right\} $

引入势函数 ${\phi _{\rm{p}}}$、ψ_p，对桩基波动方程进行Helmholtz分解，可得

(15) $ \left. \begin{gathered} {\nabla ^2}{\phi _{\rm{p}}}+k_{{\rm{l}}2}^2{\phi _{\rm{p}}} = 0\begin{array}{*{20}{c}} , \end{array} \\ {\nabla ^2}{\psi _{\rm{p}}}+k_{{\rm{t}}2}^2{\psi _{\rm{p}}} = 0\begin{array}{*{20}{c}} . \end{array} \\ \end{gathered} \right\} $

(16) $ \left. \begin{gathered} \sigma _r^{\rm{p}} = {\lambda _{\rm{p}}}{\nabla ^2}{\phi _{\rm{p}}}+2{\mu _{\rm{p}}}\frac{\partial }{{\partial r}}\left( {\frac{{\partial {\phi _{\rm{p}}}}}{{\partial r}}+\frac{1}{r}\frac{{\partial {\psi _{\rm{p}}}}}{{\partial \theta }}} \right)\begin{array}{*{20}{c}} , \end{array} \\ \sigma _{r\theta }^{\rm{p}} = {\mu _{\rm{p}}}\left( {\frac{2}{r}\frac{{{\partial ^2}{\phi _{\rm{p}}}}}{{\partial r\partial \theta }} - \frac{2}{{{r^2}}}\frac{{\partial {\phi _{\rm{p}}}}}{{\partial \theta }} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{\psi _{\rm{p}}}}}{{\partial {\theta ^2}}} - \frac{{{\partial ^2}{\psi _{\rm{p}}}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {\psi _{\rm{p}}}}}{{\partial r}}} \right). \end{gathered} \right\} $

式中：k_l2、k_t2分别为桩基中压缩波和剪切波的波数.

根据弹性波散射理论^[25]，将土体和桩体的势函数 ${\phi _{\rm{s}}}$、ψ_s、 ${\phi _{\rm{p}}}$、ψ_p用柱波函数展开，可得

(17) $ \left. \begin{gathered} {\phi _{\rm{s}}} = \sum\limits_{m = - \infty }^\infty {\left[ {{A_{m1}}{J_m}\left( {{k_{{\rm{l}}1}}r} \right)+{B_{m1}}{Y_m}\left( {{k_{{\rm{l}}1}}r} \right)} \right]{{\rm{exp}}}\; \left( {{\rm{i}}m\theta } \right) , \begin{array}{*{20}{c}} \end{array}} \\ {\psi _{\rm{s}}} = \sum\limits_{m = - \infty }^\infty {\left[ {{A_{m2}}{J_m}\left( {{k_{{\rm{t}}1}}r} \right)+{B_{m2}}{Y_m}\left( {{k_{{\rm{t}}1}}r} \right)} \right]{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right) , \begin{array}{*{20}{c}} \end{array}} \\ {\phi _{\rm{p}}} = \sum\limits_{m = - \infty }^\infty {{C_{m1}}{J_m}\left( {{k_{{\rm{l}}2}}r} \right){{\rm{exp}}} \left( {{\rm{i}}m\theta } \right) , \begin{array}{*{20}{c}} \end{array}} \\ {\psi _{\rm{p}}} = \sum\limits_{m = - \infty }^\infty {{C_{m2}}{J_m}\left( {{k_{{\rm{t}}2}}r} \right){{\rm{exp}}} \left( {{\rm{i}}m\theta } \right).\begin{array}{*{20}{c}} \end{array}} \\ \end{gathered} \right\} $

式中：A_m1、A_m2、B_m1、B_m2、C_m1、C_m2为待定积分常数，通过桩-土界面条件确定；J_m、Y_m分别为m阶第1类和第2类Bessel函数，其中m为整数.

将式（17）代入式（12）、（14）、（16），可得土体和桩基位移场及应力场的级数表达式：

(18) $ \begin{gathered} \left[ {\begin{array}{*{20}{c}} {u_r^{\rm{s}}} \\ {u_\theta ^{\rm{s}}} \end{array}} \right] = \sum\limits_{m = - \infty }^{+\infty } {{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right)\left\{ {{{\boldsymbol{F}}_m}\left[ {\begin{array}{*{20}{c}} {{A_{m1}}} \\ {{A_{m2}}} \end{array}} \right]+{{\boldsymbol{G}}_m}\left[ {\begin{array}{*{20}{c}} {{B_{m1}}} \\ {{B_{m2}}} \end{array}} \right]} \right\}} \begin{array}{*{20}{c}} , \end{array} \\ \left[ {\begin{array}{*{20}{c}} {\sigma _r^{\rm{s}}} \\ {\sigma _{r\theta }^{\rm{s}}} \end{array}} \right] = \sum\limits_{m = - \infty }^{+\infty } {{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right)\left\{ {{{\boldsymbol{f}}_m}\left[ {\begin{array}{*{20}{c}} {{A_{m1}}} \\ {{A_{m2}}} \end{array}} \right]+{{\boldsymbol{g}}_m}\left[ {\begin{array}{*{20}{c}} {{B_{m1}}} \\ {{B_{m2}}} \end{array}} \right]} \right\}} \begin{array}{*{20}{c}} , \end{array} \\ \left[ {\begin{array}{*{20}{c}} {u_r^{\rm{p}}} \\ {u_\theta ^{\rm{p}}} \end{array}} \right] = \sum\limits_{m = - \infty }^{+\infty } {{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right){{\boldsymbol{D}}_m}\left[ {\begin{array}{*{20}{c}} {{C_{m1}}} \\ {{C_{m2}}} \end{array}} \right]} \begin{array}{*{20}{c}} , \end{array} \\ \left[ {\begin{array}{*{20}{c}} {\sigma _r^{\rm{p}}} \\ {\sigma _{r\theta }^{\rm{p}}} \end{array}} \right] = \sum\limits_{m = - \infty }^{+\infty } {{{\boldsymbol{d}}_m}{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right)\left[ {\begin{array}{*{20}{c}} {{C_{m1}}} \\ {{C_{m2}}} \end{array}} \right]} \begin{array}{*{20}{c}} . \end{array} \\ \end{gathered} $

式中：系数矩阵F_m、G_m、f_m、g_m、D_m、d_m均为已知的二阶矩阵，具体表达式见附录.

根据桩-土界面处（r = r₀）的位移和应力连续，得到界面条件如下：

(19) $ u_r^{\rm{s}} = u_r^{\rm{p}} \text{，} u_\theta ^{\rm{s}} = u_\theta ^{\rm{p}} \text{，} \sigma _r^{\rm{s}} = \sigma _r^{\rm{p}} \text{，} \sigma _{r\theta }^{\rm{s}} = \sigma _{r\theta }^{\rm{p}} . $

将式（18）代入式（19），可得待定系数A_m1、A_m2与B_m1、B_m2的关系：

(20) $ \left[ {\begin{array}{*{20}{c}} {{A_{m1}}} \\ {{A_{m2}}} \end{array}} \right] = {{\boldsymbol{P}}_m}\left[ {\begin{array}{*{20}{c}} {{B_{m1}}} \\ {{B_{m2}}} \end{array}} \right] . $

将式（20）代入式（18），可以将土体位移写为

(21) $ \left[ {\begin{array}{*{20}{c}} {u_r^{\rm{s}}} \\ {u_\theta ^{\rm{s}}} \end{array}} \right] = \sum\limits_{m = - \infty }^{+\infty } {{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right){{\boldsymbol{Q}}_m}\left[ {\begin{array}{*{20}{c}} {{B_{m1}}} \\ {{B_{m2}}} \end{array}} \right]} . $

式中：P_m、Q_m为已知的二阶矩阵，表达式见附录.

为了推导桩-土周期结构的弥散方程，须将极坐标系中的土体位移转换到元胞系统的直角坐标系（x-y）中，转换关系如下：

(22) $ \begin{split} &\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos\;\theta }&{ - \sin\;\theta } \\ {\sin\;\theta }&{\cos\;\theta } \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {u_r^{\rm{s}}} \\ {u_\theta ^{\rm{s}}} \end{array}} \right] = \\ & \sum\limits_{m = - \infty }^{+\infty } {{{\rm{exp}}} \left( {{\rm{i}}m\theta } \right){{\boldsymbol{{\bar {\boldsymbol{Q}}}}}_m}\left[ {\begin{array}{*{20}{c}} {{B_{m1}}} \\ {{B_{m2}}} \end{array}} \right]} . \end{split} $

式中：

$ {{\boldsymbol{\bar Q}}_m} = \left[ {\begin{array}{*{20}{c}} {\cos\;\theta }&{ - \sin\;\theta } \\ {\sin\;\theta }&{\cos\;\theta } \end{array}} \right]{{\boldsymbol{Q}}_m} . $

将图1(c)中元胞4个边界上的法向位移表示为

(23) $ \frac{\partial }{{\partial {{n}}}}{\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma _i} } = {\boldsymbol{\varOmega}} {\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma_ i}} . $

式中：n= [n_x, n_y, −n_x, −n_y]为4个边界对应的单位法向量；系数矩阵Ω为8N×8N阶矩阵，矩阵元素由 ${{\boldsymbol{\bar Q}}_m}$的矩阵元素组成，其中N为每个边界上选取的计算点数. Ω与 ${{\boldsymbol{\bar Q}}_m}$的矩阵元素对应关系如下：

(24) $ {\boldsymbol{\varOmega}} = {{\boldsymbol{\varOmega}} _2}{\boldsymbol{\varOmega}} _1^{ - 1} . $

式中：

$ {\left( {{{{\varOmega}} _1}} \right)_{jk}} = \left[ {{{{\boldsymbol{\bar Q}}}_{k - 1 - 2N}}\left( {{r_j},{\theta _j}} \right)} \right] \text{，} $

${\left( {{{{\varOmega}} _2}} \right)_{jk}} = \dfrac{\partial }{{\partial {{n}}}}\left[ {{{{\boldsymbol{\bar Q}}}_{k - 1 - 2N}}\left( {{r_j},{\theta _j}} \right)} \right]$；j, k = 1, 2, ···, 4N.

根据Bloch–Floquet定理可知，周期结构中的所有场量Φ均满足如下周期条件：

(25) $\left. \begin{array}{l} \varPhi \left( {x+a} \right) = \varPhi \left( x \right){{\rm{exp}}\;({{\rm{i}}{k_x}a})} ，\\ \varPhi (y + a) = \varPhi (y){\rm{exp}}\;({\rm{i}}{k_y}a).\end{array} \right\}$

式中：X为坐标矢量，X = [x, y]；k为波矢，k = [k_x, k_y].

对元胞边界上的位移运用周期条件（25），可得

(26) $\left. \begin{split} &{\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma_3}} = {{\rm{exp}}\;({{\rm{i}}{k_y}a})}{\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma_1}},\\ & {\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma_4}} = {{\rm{exp}}\;({{\rm{i}}{k_x}a})}{\left[ {\begin{array}{*{20}{c}} {u_x^{\rm{s}}} \\ {u_y^{\rm{s}}} \end{array}} \right]_{\varGamma_2}} . \end{split} \right\}$

将式（26）代入式（23），可得桩-土系统的特征方程：

(27) $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varOmega}} _{31}}}&{{{\boldsymbol{\varOmega}} _{32}}} \\ {{{\boldsymbol{\varOmega}} _{41}}}&{{{\boldsymbol{\varOmega}} _{42}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\left. {{{\boldsymbol{u}}^{\rm{s}}}} \right|}_{\varGamma_1}}} \\ {{{\left. {{{\boldsymbol{u}}^{\rm{s}}}} \right|}_{\varGamma_2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol M}_{11}}}&{{{\boldsymbol M}_{12}}} \\ {{{\boldsymbol M}_{21}}}&{{{\boldsymbol M}_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\left. {{{\boldsymbol{u}}^{\rm{s}}}} \right|}_{\varGamma_1}}} \\ {{{\left. {{{\boldsymbol{u}}^{\rm{s}}}} \right|}_{\varGamma_2}}} \end{array}} \right] . $

式中：Ω₃₁、Ω₃₂、Ω₄₁、Ω₄₂为Ω中的二维分块矩阵；

$ \begin{split} &{{\boldsymbol M}_{11}} = {{\boldsymbol{\varOmega}} _{13}}{\rm{exp}}\;(2{\rm{i}}{k_y}a)+\left( {{{\boldsymbol{\varOmega}} _{11}} - {{\boldsymbol{\varOmega}} _{33}}} \right){\rm{exp}}\;({\rm{i}}{k_y}a) \text{，} \\ &{{\boldsymbol M}_{12}} = {{\boldsymbol{\varOmega}} _{12}}{\rm{exp}}\;({\rm{i}}{k_y}a) - {{\boldsymbol{\varOmega}} _{34}}{\rm{exp}}\;({\rm{i}}{k_x}a) + {{\boldsymbol{\varOmega}} _{14}}{{\rm{exp}}\;({{\rm{i}}\left( {{k_x} + {k_y}} \right)a)}} \text{，} \\ &{{\boldsymbol M}_{21}} = {{\boldsymbol{\varOmega}} _{21}}{\rm{exp}}\;({\rm{i}}{k_x}a) + {{\boldsymbol{\varOmega}} _{23}}{{\rm{exp}}\;({{\rm{i}}\left( {{k_x} + {k_y}} \right)a)}} - {{\boldsymbol{\varOmega}} _{43}}{\rm{exp}}\;({\rm{i}}{k_y}a) \text{，} \\ &{{\boldsymbol M}_{22}} = {{\boldsymbol{\varOmega}} _{24}}{\rm{exp}}\;(2{\rm{i}}{k_x}a)+\left( {{{\boldsymbol{\varOmega}} _{22}} - {{\boldsymbol{\varOmega}} _{44}}} \right){\rm{exp}}\;({\rm{i}}{k_x}a) . \end{split} $

将式（27）进一步整理为标准方程形式：

(28) $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol T}_1}{{\rm{exp}}\;{({\rm{i}}{k_{x,y}}a)}}+{{\boldsymbol T}_2}}&{{{\boldsymbol T}_3}} \\ { - {\boldsymbol{I}}}&{{\boldsymbol{I}}{{\rm{exp}}\;{({\rm{i}}{k_{x,y}}a)}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\boldsymbol{V}} \\ {\boldsymbol{U}} \end{array}} \right] = {\bf{{0}}} . $

式（28）为8N×8N阶矩阵方程，其中I为4N×4N阶单位矩阵； ${\boldsymbol{V}} = {\boldsymbol{U}}{{\rm{exp}}\;{({\rm{i}}{k_{x,y}}a)}}$，其中U为元胞边界上的位移矩阵， ${\boldsymbol{U}} = {[ {{{ {u_x^{\rm{s}}} |}_{\varGamma_1}},{{ {u_y^{\rm{s}}} |}_{\varGamma_1}},{{ {u_x^{\rm{s}}} |}_{\varGamma_2}},{{ {u_y^{\rm{s}}} |}_{\varGamma_2}}} ]^{\rm{T}}}$；矩阵T₁、T₂、T₃的表达式见附录.

利用非凡解条件，可得

(29) $ \det \left[ {\begin{array}{*{20}{c}} {{{ {\boldsymbol{T}}}_1}{{\rm{exp}}\;{({\rm{i}}{k_{x,y}}a)}}+{{ {\boldsymbol{T}}}_2}}&{{{ {\boldsymbol{T}}}_3}} \\ { - {\boldsymbol{I}}}&{{\boldsymbol{I}}{{\rm{exp}}\;{({\rm{i}}{k_{x,y}}a)}}} \end{array}} \right] = 0 . $

式（29）为桩-土周期结构的弥散方程，det表示取行列式. 对方程给定角频率ω，在第一不可约布里渊域Γ-X-M（见图3）中搜索满足方程的波矢k_x、k_y，可得周期系统中压缩波和剪切波的弥散曲线（频率-波数关系）. 当搜索Γ-X边界时（0 ≤ k_x ≤ π/a），方程中的指数项取 ${{\rm{exp}}\;{({\rm{i}}{k_x}a)}}$；当搜索X-M（0 ≤ k_y ≤ π/a）边界时，指数项取 ${{\rm{exp}}\;{({\rm{i}}{k_y}a)}}$.

对本文的求解方法和计算程序进行验证. 可知，当τ = 0时，土体模型退化为线弹性模型. Huang等^[26]利用平面波展开法，计算线弹性土排桩地基的禁振带隙. 图4给出本文退化模型Γ-X边界上的剪切波弥散曲线，所需的计算参数取值如表1所示，并与Huang的结果进行对比. 图中，f^*为归一化无量纲频率，f^*= ωa/(2πc_t2). 对比发现，本文解与Huang的结果能够吻合. 图4中的阴影区域为周期结构的禁振带隙，当0 ≤ ω ≤ 500 rad/s时出现了2个带隙，f_u、f_l分别为带隙的上界频率和下界频率，带隙宽度W^* = f_u−f_l.

为了对比流变性对带隙的影响，利用Rayleigh黏弹性模型^[27]计算阻尼比为常数时的带隙，与流变性模型进行结果对比. 如图5所示为Rayleigh模型的常阻尼比ξ = 0.177与流变模型(ξ_d)_max = 0.177的结果对比，(ξ_d)_max = 0.177对应于η = 1.5. 对比发现，Rayleigh模型的计算结果与线弹性无阻尼模型的结果无明显差异. 流变模型的计算结果显示，在考虑土体流变性后，带隙明显上移，带宽小幅减小.

通过改变Rayleigh模型的常阻尼比ξ，考察ξ对计算结果的影响（计算结果见表2）. 从表2可知，只有当阻尼比非常大（接近0.7）时，才会对带隙产生影响，但实际土体的阻尼比不可能达到0.7. 由此可见，常阻尼模型不能反映土体流变性对排桩地基带隙的影响. 相反，流变性模型的分析结果表明，土体流变性会引起带隙位置上移，带宽减小. 即便在阻尼比很小的情况下，这种现象依然很明显.

如图3所示为η增大时带隙的变化曲线（即阻尼比幅值增大）. 可知，前半段曲线（0 < η < 6）带隙小幅上移、带宽减小，在 η > 6后（对应于( ξ_d)_max > 0.544），带隙上移加快，由于下界频率增速显著慢于上界频率，带宽显著增大. 实际土体的模量比不会大于6，阻尼比达不到0.5，因此 图6中后半段曲线在实际中不会出现.

η决定阻尼比的幅值，τ决定阻尼比随频率的变化速率. 如图7所示为带隙随τ的变化曲线，由此分析阻尼比的变化速率对带隙的影响. 可知，随着松弛时间的增大，带隙大幅上移，带宽显著减小，当松弛时间达到0.01 s时趋于收敛，带隙不再随松弛时间变化. 对比图6、7可以看出，松弛时间对带隙的影响明显大于模量比. 由此证明，流变性影响带隙的因素主要是阻尼的时变效应，而非阻尼值.

如图8所示为带隙随桩基填充率F_s的变化曲线， ${F_{\rm{s}}} = \text{π} r{}_0^2/\mathop a\nolimits^2$. 可知，尽管流变性导致带隙位置上移且带宽减小，但带隙随F_s的变化规律与线弹性土地基一致，即随着填充率的增大，带隙宽度先增大后减小. 当松弛时间较大时，带隙宽度随F_s的变化率减小（见图8(b)中τ=10⁻² s对应的曲线），说明流变性能够降低填充率对带隙的影响.

（1）流变性土的阻尼比随频率发生非单调性变化，始终态模量比决定阻尼比的幅值，而松弛时间决定阻尼比随频率的变化速率. 阻尼比的幅值随着模量比的增大而增大，松弛时间越大，阻尼比随频率的变化速率越大.

（2）土体流变性引起排桩地基的带隙位置上移，且带宽减小. 模量比和松弛时间越大，则带隙位置上移，带宽减小越明显. 松弛时间对带隙的影响程度显著大于模量比，这是土体流变性影响带隙的主要因素.

（3）流变性减弱了桩基填充率对带隙的影响，土体流变性越强，带隙随填充率的变化幅度越小.

（4）受土体流变性的影响，实际工程中排桩地基的带隙频率高于理论值，且带宽减小，减弱了排桩的隔振效果，消除桩周土的流变性将有利于排桩发挥隔振作用.