基于人工弹簧模型的周期排桩计算方法与特性

Calculation Method and Characteristics of Periodic Pile Arrangement Based on Artificial Spring Model

摘要

周期性排桩作为一种新型减振体系，可实现对振动的抑制甚至隔离，在隔振减振领域具有重要意义。土体中周期性排桩存在衰减域，对处于衰减域内的振动具有明显的隔离作用。本文基于周期理论和能量原理引入人工弹簧模型模拟周期单元边界和桩土密接边界，将边界约束转化为弹簧的弹性势能，以克服传统能量法难以构造满足边界条件位移形函数的问题。扫描波数计算衰减域时，只有周期边界的弹性势能对应的刚度矩阵包含波数项，质量矩阵和其余刚度矩阵不含波数项，不需要重复计算，大幅度减少了计算量。对形函数项数和人工弹簧刚度进行收敛性分析，并与有限元法、平面波展开法对比，证明本文方法不仅具有较高的准确性，还具有更高的计算效率。对比六边形和正方形不同排桩形式发现，当填充率相同时，六边形排桩比正方形排桩的衰减域宽度大；对比管桩和实心桩不同桩型发现，当填充率相同时，管桩的衰减域宽度比实心桩更大。此外，分析填充土体对衰减域的影响发现，衰减域宽度随填充土体密度增大而增大，随填充土体弹性模量增大先增大后减小且在弹性模量约16 MPa时衰减域宽度最大。

Abstract

As a new type of vibration reduction system, periodic row piles achieve vibration suppression and even isolation, which is of great significance in the field of vibration isolation and damping. Periodic row piles in soil exhibit an attenuation domain, which has a significant isolation effect on vibrations within this domain. This study introduces an artificial spring model based on periodic theory and the energy principle to simulate the boundary of periodic elements and the closely connected interface between piles and soil. The boundary constraints are transformed into the elastic potential energy of the spring, overcoming the difficulty of constructing displacement shape functions that satisfy boundary conditions using traditional energy methods. When calculating the attenuation domain by scanning the wave number, only the stiffness matrix corresponding to the elastic potential energy of the periodic boundary contains the wave number term. In contrast, the mass matrix and other stiffness matrices do not include this term, eliminating the need for repeated calculations and significantly reducing the computational workload. The convergence analysis of the number of terms in the shape function and the stiffness of the artificial spring, compared to the finite element method, demonstrates that the proposed method not only provides high accuracy but also improves computational efficiency. When comparing different forms of hexagonal row piles and square row piles with the same filling rate, the attenuation domain width of hexagonal row piles is larger than that of square row piles. When comparing different pile types, such as pipe piles and solid piles, with the same filling rate, the attenuation domain width of pipe piles is also larger than that of solid piles. In addition, analyzing the influence of filled soil on the attenuation domain indicates that the width of the attenuation domain increases with the density of the filled soil, and first increases and then decreases with the increase of the elastic modulus of the filled soil. The maximum width of the attenuation domain is obtained when E=16 MPa.

Graphical abstract

关键词

岩土工程 / 减振 / 周期排桩 / 衰减域 / 能量法 / 人工弹簧

Key words

geotechnical engineering / vibration reduction / periodic pile arrangement / attenuation domain / energy method / artificial spring

本刊网刊

随着经济社会的快速发展，城市建筑愈发密集，市内交通不断完善并呈立体化、多样化发展，环境振动问题日益突出。为减少环境有害振动，研究人员对减振、隔振进行了大量研究。工程实践和研究表明，周期性排桩结构可以有效隔离环境振动^[1‒3]。

Richart等^[4]首次提出在土体中设置单排或多排薄壁衬砌圆孔作为隔振屏障，用来减弱环境振动的传播。Kattis等^[5]将单排桩模型等效为连续屏障模型，并对单排桩减振进行了参数分析。Tsai等^[6]应用3维边界元法研究了圆形钢管桩、混凝土空心桩、混凝土实心桩对基础竖向振动的隔离效果。高广运等^[7]通过理论和试验验证了圆截面排桩具有隔振作用。Huang 等^[8‒10]基于周期理论研究了单相土多排桩特性，通过理论计算和数值模拟，解释了多排桩在特定频段内具有隔振作用的机理。孟庆娟^[11]应用COMSOL PDE有限元法研究了饱和土中周期排桩的减振特性。

在对周期性排桩的研究中，衰减域的计算十分重要，相关计算方法主要有有限元法^[12‒15]、能量法^[16‒17]、传递矩阵法^[18‒19]、平面波展开法^[20‒21]、多重散射理论^[22‒23]等。有限元法在计算周期排桩结构衰减域方面适用性强，但其精度取决于网格精细程度，获得高精度结果需要消耗较高的计算资源。当计算频率较高、结构尺寸较大时，会导致分析精度不够、计算效率较低。Guo等^[24]采用人工弹簧模拟周期边界条件，将边界约束转变为人工弹簧的弹性势能，用以分析周期结构的带隙问题。能量法将微分方程边值问题转变为泛函极值问题，有利于复杂边界条件问题的求解，因而被广泛应用于周期结构衰减域的计算。但传统能量法（如瑞利里兹法）在计算排桩结构衰减域时需要构造满足周期性边界条件的位移场形函数，在数学上难度较大；且重构以后的位移场形函数包含波数，这会导致涉及形函数的质量矩阵和刚度矩阵中也含有波数，在计算带隙时，结构的质量和刚度矩阵需随着波数的变化进行反复计算，增加了计算成本。

本文基于能量法，引入人工弹簧模拟周期单元边界条件及桩土密接边界条件，解耦形函数与边界条件，解决了传统能量法构造满足边界条件形函数困难的问题。利用人工弹簧模拟周期边界和桩土密接边界，将边界约束转化为人工弹簧的弹性势能，在扫描波数时，只有周期边界的弹性势能对应的刚度矩阵包含波数项，质量矩阵和其余刚度矩阵不含波数项，不需重复计算，极大地减少了计算量。

1 人工弹簧法计算周期排桩带隙

以图1（a）中的周期性排桩为研究对象，按正方形周期性排布，2a为元胞常数，R为桩半径，取出单个桩土胞元进行研究，如图1（b）所示。正方形周期性排桩第一布里渊区（first Brillouin zone）如图1（c）所示^[25]，正方形区域为第一布里渊区，灰色三角形区域为简约布里渊区，简约布里渊区由几个高对称点控制：\Gamma =(\mathrm{0,0})，X=\frac{\mathrm{2}\mathrm{\pi }}{a}\left(\frac{\mathrm{1}}{\mathrm{2}},\mathrm{0}\right)，M=\frac{\mathrm{2}\mathrm{\pi }}{a}\left(\frac{\mathrm{1}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}}\right)。图1（c）中，k_x\mathrm{、}{k}_y为波矢方向。正方形排桩填充率f_{\mathrm{s}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{e}}=(\mathrm{\pi }{R}^{\mathrm{2}})/{(\mathrm{2}a)}^{\mathrm{2}}。

平面应力问题中，周期单元x、y方向位移u\left(x,y,t\right)和v\left(x,y,t\right)可表示为基函数f_i\left(x,y\right)和未知权重系数d_i(t)的组合。以x方向为例，u\left(x,y,t\right)可表示为：

u(x,y,t)={\mathit{d}}^{\mathrm{T}}\mathit{f}=\mathit{f}{\mathit{d}}^{\mathrm{T}}（1）

其中，

\left\{\begin{array}{l}\mathit{f}=\mathit{\alpha }(x)\otimes \mathit{\beta }(y),\\ \mathit{\alpha }=[{\alpha }_{\mathrm{1}}(x),{\alpha }_{\mathrm{2}}(x),\cdots ,{\alpha }_i(x),\cdots ,{\alpha }_m(x)],\\ \mathit{\beta }=[{\beta }_{\mathrm{1}}(y),{\beta }_{\mathrm{2}}(y),\cdots ,{\beta }_i(y),\cdots ,{\beta }_n(y)]\end{array}\right.（2）

式（1）、（2）中， f 为形函数， d 为d_i(t)矩阵形式，符号\otimes表示克罗内克积， α 、 β 分别为x、y方向的容许形函数，m和n分别为x方向和y方向的截断项数。

本文选用切比雪夫级数作为周期排桩的位移形函数^[26]，具体形式为：

\left\{\begin{array}{c}\begin{array}{l}{\alpha }_{\mathrm{1}}(x)=\mathrm{1},{\alpha }_{\mathrm{2}}(x)=x,\\ {\alpha }_m(x)=\mathrm{2}x\cdot {\alpha }_{m-\mathrm{1}}(x)-{\alpha }_{m-\mathrm{2}}(x);\end{array}\\ \begin{array}{l}{\beta }_{\mathrm{1}}(y)=\mathrm{1},{\beta }_{\mathrm{2}}(y)=y,\\ {\beta }_n(y)=\mathrm{2}y\cdot {\beta }_{n-\mathrm{1}}(y)-{\beta }_{n-\mathrm{2}}(y)\end{array}\end{array}\right.（3）

根据Bloch定理，周期单元在边界处的周期性边界条件为^[27]：

\left\{\begin{array}{l}u(-a,y)=u(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a},\\ v(-a,y)=v(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a}\end{array}\right.（4）

\left\{\begin{array}{l}u(x,-a)=u(x,-a){\mathrm{e}}^{-i{k}_y\mathrm{2}a},\\ v(x,-a)=v(x,-a){\mathrm{e}}^{-i{k}_y\mathrm{2}a}\end{array}\right.（5）

桩土密接边界条件为：

\left\{\begin{array}{l}{u}_{\mathrm{p}}(x_r,y_r)-u_{\mathrm{s}}(x_r,y_r)=\mathrm{0},\\ v_{\mathrm{p}}(x_r,y_r)-v_{\mathrm{s}}(x_r,y_r)=\mathrm{0}\end{array}\right.（6）

式中，(x_r,y_r)为桩土接触面处坐标，u_{\mathrm{p}}\left(x_r,y_r\right)-u_{\mathrm{s}}\left(x_r,y_r\right)为桩和土在(x_r,y_r)处法向位移差值，v_{\mathrm{p}}\left(x_r,y_r\right)-v_{\mathrm{s}}\left(x_r,y_r\right)为桩和土在(x_r,y_r)处切向位移差值。

形函数\mathit{f}需满足Bloch周期性条件和桩土密接条件，构造满足条件的形函数困难较大，限制了传统能量法在周期排桩问题中的使用。故本文采用人工弹簧处理周期性边界和桩土密接边界，以解耦边界条件与位移形函数^[28‒29]。在单元边界上定义x向和y向线弹簧，物理含义是在x和y方向通过人工弹簧将相邻周期单元连接，如图2（a）所示。桩土边界通过压缩弹簧和剪切弹簧连接，如图2（b）所示。

周期边界势能U_{\mathrm{p}}^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}表示为：

U_{\mathrm{p}}^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}=\frac{\mathrm{1}}{\mathrm{2}}\left\{{\int }_{-a}^a[k_{\mathrm{t}}{\left(u(-a,y)-u(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a}\right)}^{\mathrm{2}}\right.+k_{\mathrm{s}}{\left(v(-a,y)-v(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a}\right)}^{\mathrm{2}}]\mathrm{d}y+{\int }_{-a}^a[k_{\mathrm{s}}{\left(u(x,-a)-u(x,a){\mathrm{e}}^{-i{k}_y\mathrm{2}a}\right)}^{\mathrm{2}}+\left.k_{\mathrm{t}}{\left(v(x,-a)-v(x,a){\mathrm{e}}^{-i{k}_y\mathrm{2}a}\right)}^{\mathrm{2}}]\mathrm{d}x\right\}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}(k_x,k_y)\mathit{d}（7）

式中，u(-a,y)-u(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a}、v(-a,y)-v(a,y){\mathrm{e}}^{-i{k}_x\mathrm{2}a}分别为相邻周期首尾端x、y向位移差值，k_t、k_s分别为拉压弹簧和剪切弹簧的刚度，{\mathit{K}}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}(k_x,k_y)为单胞边界刚度矩阵。

由式（4）、（5）可知，x向、y向位移差值趋于0，即

\left\{\begin{array}{c}\mathrm{l}\mathrm{i}\mathrm{m}\left[u\left(-a,y\right)-u\left(a,y\right){\mathrm{e}}^{-i{k}_{x}a}\right]=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{t}}}{k_{\mathrm{t}}}\right)=\mathrm{0},\\ \mathrm{l}\mathrm{i}\mathrm{m}\left[v\left(-a,y\right)-v\left(a,y\right){\mathrm{e}}^{-i{k}_{x}a}\right]=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{s}}}{k_{\mathrm{s}}}\right)=\mathrm{0}\end{array}\right.（8）

\left\{\begin{array}{c}\mathrm{l}\mathrm{i}\mathrm{m}\left[u\left(x,-a\right)-u\left(x,a\right){\mathrm{e}}^{-i{k}_{x}a}\right]=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{s}}}{k_{\mathrm{s}}}\right)=\mathrm{0},\\ \mathrm{l}\mathrm{i}\mathrm{m}\left[v\left(x,-a\right)-v\left(x,a\right){\mathrm{e}}^{-i{k}_{x}a}\right]=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{t}}}{k_{\mathrm{t}}}\right)=\mathrm{0}\end{array}\right.（9）

因为式（8）、（9）中结构压力F_{\mathrm{t}}和剪力F_{\mathrm{s}}都为有限值，所以当k_t、k_s足够大时，式（8）、（9）才能成立。式（7）中，U_{\mathrm{p}}^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}已包含周期性边界条件，故形函数\mathit{f}的构造不需要满足周期性边界。

桩土接触面之间弹性势能U_{\mathrm{p}}^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}可以表示为：

\begin{array}{l}{U}_{\mathrm{p}}^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}\left\{\oint k_{\mathrm{t}}{\left[u_{\mathrm{p}}(x_r,y_r)-u_{\mathrm{s}}(x_r,y_r)\right]}^{\mathrm{2}}\mathrm{d}s\right.+\\ \left.\oint k_{\mathrm{s}}{\left[v_{\mathrm{p}}(x_r,y_r)-v_{\mathrm{s}}(x_r,y_r)\right]}^{\mathrm{2}}\mathrm{d}s\right\}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}\mathit{d}\end{array}（10）

式中，\mathrm{d}s为桩土接触面微分长度，{\mathit{K}}_{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}为桩土接触面刚度矩阵。

根据桩土密接边界条件，有

\left\{\begin{array}{l}{u}_{\mathrm{p}}-u_{\mathrm{s}}=\mathrm{0},\\ v_{\mathrm{p}}-v_{\mathrm{s}}=\mathrm{0}\end{array}\right.（11）

可以得到：

\left\{\begin{array}{l}\mathrm{l}\mathrm{i}\mathrm{m}(u_{\mathrm{p}}-u_{\mathrm{s}})=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{t}}}{k_{\mathrm{t}}}\right)=\mathrm{0},\\ \mathrm{l}\mathrm{i}\mathrm{m}(v_{\mathrm{p}}-v_{\mathrm{s}})=\mathrm{l}\mathrm{i}\mathrm{m}\left(-\frac{F_{\mathrm{s}}}{k_{\mathrm{s}}}\right)=\mathrm{0}\end{array}\right.（12）

式（10）中，U_{\mathrm{p}}^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}已包含桩土密接边界条件，故形函数\mathit{f}的构造不需要满足桩土密接边界。

周期单元中桩体的动能U_{\mathrm{p}}和应变能E_{\mathrm{p}}表示为：

\begin{array}{l}{U}_{\mathrm{p}}=\frac{\mathrm{1}}{\mathrm{2}}{\iint }_{\Omega }^{}\frac{{\widehat{E}}_{\mathrm{p}}}{\left(\mathrm{1}+\nu \text{'}\right)\left(\mathrm{1}-\mathrm{2}\nu \text{'}\right)}[\left(\mathrm{1}-\nu \text{'}\right)\left({\left(\frac{\partial u}{\partial x}\right)}^{\mathrm{2}}+{\left(\frac{\partial v}{\partial y}\right)}^{\mathrm{2}}\right)+\\ \mathrm{2}\nu \text{'}\frac{\partial u}{\partial x}·\frac{\partial v}{\partial y}+\frac{\mathrm{1}-\mathrm{2}\nu \text{'}}{\mathrm{2}}{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^{\mathrm{2}}]\mathrm{d}s=\\ \frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\left\{{\iint }_{\Omega }^{}\frac{{\widehat{E}}_{\mathrm{p}}}{\left(\mathrm{1}+\nu \text{'}\right)\left(\mathrm{1}-\mathrm{2}\nu \text{'}\right)}\left[\left(\mathrm{1}-\nu \text{'}\right)\cdot {\left(\frac{\partial {\mathit{f}}_x}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial x}\right)}^{\mathrm{T}}\right.\right.\\ \left(\frac{\partial {\mathit{f}}_x}{\partial x}·\frac{\partial {\mathit{f}}_{\mathrm{x}}^{\mathrm{T}}}{\partial x}+\frac{\partial {\mathit{f}}_y}{\partial y}·\frac{\partial {\mathit{f}}_{\mathit{y}}^{\mathbf{T}}}{\partial y}\right)+\nu \text{'}\left(\frac{\partial {\mathit{f}}_x}{\partial x}·\frac{\partial {\mathit{f}}_{\mathit{y}}^{\mathbf{T}}}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial y}·\frac{\partial {\mathit{f}}_{\mathrm{x}}^{\mathrm{T}}}{\partial x}\right)+\frac{\mathrm{1}-\mathrm{2}\nu \text{'}}{\mathrm{2}}\left(\frac{\partial {\mathit{f}}_x}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial x}\right)·\\ \left.\left.{\left(\frac{\partial {\mathit{f}}_x}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial x}\right)}^{\mathrm{T}}\right]\mathrm{d}s\right\}\mathit{d}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{p}}\mathit{d}\end{array}（13）

E_{\mathrm{p}}=\frac{\mathrm{1}}{\mathrm{2}}{\iint }_{\Omega }^{}{\rho }_{\mathrm{p}}\left({\dot{u}}^{\mathrm{2}}+{\dot{v}}^{\mathrm{2}}\right)\mathrm{d}s=\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}\left[{\iint }_{\Omega }^{}{\rho }_{\mathrm{p}}\left({\mathit{f}}_{\mathit{x}}^{\mathrm{2}}+{\mathit{f}}_{\mathit{y}}^{\mathrm{2}}\right)\mathrm{d}s\right]\mathit{d}=\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{p}}\mathit{d}（14）

式（13）、（14）中，{\mathit{K}}_{\mathrm{p}}为桩体变形刚度矩阵，{\mathit{M}}_{\mathrm{p}}为桩体质量矩阵，{\widehat{E}}_{\mathrm{p}}为桩体弹性模量，{\rho }_{\mathrm{p}}为桩体密度，\nu \text{'}为土体泊松比。{\widehat{E}}_{\mathrm{p}}和{\rho }_{\mathrm{p}}可分别表示为：

（15）

（16）

周期单元中，土体的动能U_{\mathrm{s}}和应变能E_{\mathrm{s}}表示为：

\begin{array}{l}{U}_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\left\{{\iint }_{\Omega }^{}\frac{{\widehat{E}}_{\mathrm{s}}}{\left(\mathrm{1}+\nu \text{'}\right)\left(\mathrm{1}-\mathrm{2}\nu \text{'}\right)}\left[\left(\mathrm{1}-\nu \text{'}\right)\left(\frac{\partial {\mathit{f}}_x}{\partial x}·\frac{\partial {\mathit{f}}_{\mathrm{x}}^{\mathrm{T}}}{\partial x}+\right.\right.\right.\\ \left.\frac{\partial {\mathit{f}}_y}{\partial y}·\frac{\partial {\mathit{f}}_{\mathit{y}}^{\mathbf{T}}}{\partial y}\right)+\nu \left(\frac{\partial {\mathit{f}}_x}{\partial x}·\frac{\partial {\mathit{f}}_{\mathit{y}}^{\mathbf{T}}}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial y}·\frac{\partial {\mathit{f}}_{\mathrm{x}}^{\mathrm{T}}}{\partial x}\right)+\left.\left.\frac{\mathrm{1}-\mathrm{2}\nu \text{'}}{\mathrm{2}}\left(\frac{\partial {\mathit{f}}_x}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial x}\right){\left(\frac{\partial {\mathit{f}}_x}{\partial y}+\frac{\partial {\mathit{f}}_y}{\partial x}\right)}^{\mathrm{T}}\right]\mathrm{d}s\right\}\mathit{d}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{s}}\mathit{d}\end{array}（17）

E_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}\left[{\iint }_{\Omega }^{}{\rho }_{\mathrm{s}}\left({\mathit{f}}_{\mathit{x}}^{\mathrm{2}}+{\mathit{f}}_{\mathit{y}}^{\mathrm{2}}\right)\mathrm{d}s\right]\mathit{d}=\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{s}}\mathit{d}（18）

式（17）、（18）中，{\mathit{K}}_{\mathrm{s}}为土体变形刚度矩阵，{\mathit{M}}_{\mathrm{s}}为土体质量矩阵，{\widehat{E}}_{\mathrm{s}}为土体弹性模量，{\rho }_{\mathrm{s}}为土体密度。{\widehat{E}}_{\mathrm{s}}和{\rho }_{\mathrm{s}}可分别表示为：

（19）

（20）

周期排桩单元的整体拉格朗日量为：

\begin{array}{l}L=E_{\mathrm{p}}+E_{\mathrm{s}}-U_{\mathrm{p}}-U_{\mathrm{s}}-U_{\mathrm{p}}^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}-U_{\mathrm{p}}^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}=\\ \frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{s}}\mathit{d}+\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{d}}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{p}}\mathit{d}-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{p}}\mathit{d}-\\ \frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{s}}\mathit{d}-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}\left(k_x,k_y\right)\mathit{d}-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}\mathit{d}\end{array}（21）

定义与时间相关向量\mathit{d}\left(t\right)=\widehat{\mathit{D}}{\mathrm{e}}^{\mathrm{i}\omega t}（\omega为频率，i表示复数），根据变分原理得到：

\frac{\partial }{\partial t}\left(\frac{\partial L}{\partial \dot{\mathit{d}}}\right)-\frac{\partial L}{\partial \mathit{d}}=\mathrm{0}（22）

得到周期排桩的运动方程：

\left[\mathit{K}-{\omega }^{\mathrm{2}}\mathit{M}\right]\widehat{\mathit{D}}=\mathrm{0}（23）

式中，\mathit{K}={\mathit{K}}_{\mathrm{p}}+{\mathit{K}}_{\mathrm{s}}+{\mathit{K}}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}(k_x,k_y)+{\mathit{K}}_{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{e}}，\mathit{M}={\mathit{M}}_{\mathrm{s}}+{\mathit{M}}_{\mathrm{p}}。

通过扫描第一布里渊区，有k_x\times k_y=[-\mathrm{\pi }/a,\mathrm{\pi }/a]\times [-\mathrm{\pi }/a,\mathrm{\pi }/a]，可得周期排桩结构的频散曲线。

2 算例验证

2.1 收敛性分析

位移形函数截断项m、n对计算结果有较大影响，人工弹簧刚度k_t、k_s也会影响计算结果的准确性，故针对截断项和弹簧刚度进行收敛性分析。元胞常数2a=2 m，填充率f_{\mathrm{s}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{e}}=\mathrm{0.35}，桩半径R=0.668 m，桩土材料参数见表1。计算正方形周期排桩的频散曲线，在频散曲线中选择波矢k_x=k_y=\mathrm{\pi }/a处的频率做收敛性分析。

首先，分析截断项m、n不同取值对频率的影响，结果见表2。由表2可知：随着截断项数增大，频率趋于稳定；当m=n=13时，频率值基本不再变化，表明随着形函数截断项数增加，结果满足收敛性。故后续计算取截断项m=n=13。

对人工弹簧刚度取值进行收敛性分析，结果如图3所示。

由图3可以看出，当弹簧刚度大于10¹⁰ N/m时，频率值趋于稳定，证明了人工弹簧法在计算周期排桩问题上的收敛性。故后续计算中取弹簧刚度k_{\mathrm{t}}=k_{\mathrm{s}}={\mathrm{10}}^{\mathrm{10}}\mathrm{N}/\mathrm{m}。

2.2 准确性分析

通过计算正方形周期排桩的频散曲线并与COM-SOL有限元软件和IPWE计算结果^[25]对比，验证本文方法计算周期结构带隙的准确性。相关参数取值同第2.1节，周期排桩频散关系曲线如图4所示。由图4可知，在33.3~37.3 Hz范围产生了一个完全的衰减域（阴影部分），本文方法与COMSOL、IPWE计算结果^[25]吻合良好，验证了人工弹簧计算周期排桩频散曲线的准确性。

周期性排桩除了正方形排桩外，还有六边形排桩。图5（a）为六边形周期性排布形式，元胞常数为2a；图5（b）为单个桩土胞元，桩半径为R；图5（c）为六边形周期性排桩第一布里渊区和简约布里渊区。六边形排桩填充率f_{\mathrm{h}\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}}=\mathrm{\pi }{R}^{\mathrm{2}}/\left(\mathrm{2}\sqrt[]{\mathrm{3}}{a}^{\mathrm{2}}\right)。

对于六边形周期排桩，{\widehat{E}}_{\mathrm{s}}、{\rho }_{\mathrm{s}}取值如下：

（24）

（25）

采用人工弹簧法得到六边形周期排桩频散曲线，并与有限元结果对比，其中，元胞常数2a=2 m，填充率f_{\mathrm{h}\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}}=\mathrm{0.35}，桩半径r=0.621 m，结果如图6所示。由图6可以看出二者吻合良好。

2.3 效率对比分析

分析不同扫描波数时的计算时间以对比人工弹簧法和有限元法的计算效率，结果见表3。由表3可知，人工弹簧法的计算效率是有限元法的3倍以上。这是因为有限元法的动能和应变能矩阵中含有波矢k_x\mathrm{、}{k}_y，每次扫描波数时动能和应变能矩阵都需计算；而人工弹簧法只在边界弹性势能矩阵中含有波矢k_x\mathrm{、}{k}_y，动能和应变能矩阵不需要重复计算，减少了计算时间。

3 填充率与带隙关系

在周期性排桩设计时，需要考虑填充率对衰减域的影响。对比图4和6可知，当填充率为0.35时，六边形周期性排桩的衰减域宽度WAZ（width of attenuation zone）比正方形周期性排桩宽。故主要分析六边形排桩时填充率f_{\mathrm{h}\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}}与第一完全衰减域3个特征参量的关系，即与起始频率（lower bound frequency，LBF）、截止频率（upper bound frequency，UBF）和WAZ的关系。

3.1 实心桩填充率与衰减域的关系

图7为f_{\mathrm{h}\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}}与UBF、LBF、WAZ的关系。由图7可知，随着填充率增加，起始频率缓慢增大，截止频率比起始频率增大更多，衰减域宽度增大。根据布拉格散射型周期结构特性，当填充率逐渐增大、桩间净距逐渐减小，周期结构排桩能增大对弹性波的屏障效果。

3.2 管桩填充率与衰减域的关系

管桩不仅可以减少施工时的挤土效应，还可以减少混凝土用量。管桩形式如图8所示，管桩中间为填充土体。管桩填充率f_{\mathrm{p}\mathrm{i}\mathrm{p}\mathrm{e}}=\mathrm{\pi }{(R-r)}^{\mathrm{2}}/(\mathrm{2}\sqrt[]{\mathrm{3}}{a}^{\mathrm{2}})。分析f_{\mathrm{p}\mathrm{i}\mathrm{p}\mathrm{e}}与UBF、LBF、WAZ的关系，结果如图9所示。由图9可知：随着填充率的增加，LBF缓慢降低，UBF先增加随后略有下降；衰减域宽度先增大后略微下降，最宽衰减域在填充率为0.17附近。增加桩半径可以获得更宽的衰减域，但当半径较大时，桩之间净距较小，衰减域宽度反而会减小。

与图7对比可以发现，采用相同的填充率时，管桩LBF与实心桩LBF的起始频率基本一致，而带隙宽度比实心桩更大，管桩的隔振效果更经济高效。

4 填土特性与带隙关系

为了确定管桩内最优填充土体的密度和弹性模量，进一步分析其对衰减域的影响。典型管桩单元周期常数2a=2 m，填充率f_{\mathrm{p}\mathrm{i}\mathrm{p}\mathrm{e}}=\mathrm{0.15}，外径R=0.626 m，内径r=0.476 m，壁厚R-r=0.15 m，管桩和外侧土体特性见表1。常见土体弹性模量变化范围为1到100 MPa，密度为1 700~2 200 kg/m^3[30‒31]。

4.1 填充土体密度与衰减域关系

填充土体弹性模量为20 MPa，泊松比为0.35，密度为1 700~2 200 kg/m³时，填充土体密度与衰减域的关系如图10所示。由图10可知，随着填充土体密度增大，UBF基本没变化，LBF逐渐减小，衰减域宽度逐渐增大。这说明增加填充土体密度，有利于低频减振，并且能增加衰减域范围。

4.2 填充土体弹性模量与衰减域的关系

取填充土体密度为1 900 kg/m³，泊松比为0.35，弹性模量为1~100 MPa，填充土体密度与衰减域的关系如图11所示。

从图11可以看出：随着填充土体弹性模量增大，带隙宽度先增大后减小且在弹性模量E=16 MPa附近时最大；而UBF和LBF先增大后变化范围较小。

5 结 论

采用人工弹簧模型处理周期性边界和桩土密接边界，实现了边界条件与位移形函数的解耦。首先，对人工弹簧刚度和形函数截断项进行收敛性分析，并与有限元法对比了人工弹簧法的准确性和计算效率；最后，分析了填充率与填充材料对衰减域的影响。具体结论如下：

1）人工弹簧法具有较高的准确性，计算效率是有限元法的3倍以上。

2）利用人工弹簧模拟周期性边界和桩土密接边界，解耦边界条件与位移形函数，解决了传统能量法中形函数构造困难的问题。

3）分析填充率与衰减域的关系发现，填充率相同的情况下，六边形布置比正方形布置的衰减域宽度大，六边形布置的管桩比实心桩的衰减域宽度大。

4）分析填充材料特性与衰减域的关系可知：随着填充土体密度增大，衰减域宽度逐渐增大；随着填充土体弹性模量增大，带隙宽度先增大后减小且在弹性模量E=16 MPa附近时最大。

参考文献

[1] Pu Xingbo, Shi Zhifei.A novel method for identifying surface waves in periodic structures[J].Soil Dynamics and Ea-rthquake Engineering,2017,98:67‒71. doi:10.1016/j.soildyn.2017.04.011

[2] Xia Tangdai, Chen Chen, Sun Miaomiao.An improved met-hod for multiple scattering under SV incident waves and vibration isolation using rows of piles[J].Journal of Vibration and Shock,2011,30(4):86‒90.

[3] 夏唐代,陈晨,孙苗苗.SV波入射下多重散射问题的改进算法以及双排非连续刚性屏障的隔振研究[J].振动与冲击,2011,30(4):86‒90.

[4] Xiaoduo Ou, Bai Lu, Jiang Jie,et al.Calculation and analysis of force and deformation of double-row piles supporting structure in deep foundation pit in round gravel stratum[J].Highway Engineering,2022,47(3):15‒22.

[5] 欧孝夺,白 露,江杰,等.圆砾地层中深基坑双排桩支护结构受力变形计算分析[J].公路工程,2022,47(3):15-22.

[6] Richart F E, Hall J R, Woods R D.Vibrations of soils and fo-undations[M].Englewood Cliffs:Prentice‒Hall,1970:244‒ 262.

[7] Kattis S E, Polyzos D, Beskos D E.Modelling of pile wave barriers by effective trenches and their screening effectiveness[J].Soil Dynamics and Earthquake Engineering,1999,18(1):1‒10. doi:10.1016/s0267-7261(98)00032-3

[8] Tsai P H, Feng Zhengyi, Jen T L.Three-dimensional analysis of the screening effectiveness of hollow pile barriers for foundation-induced vertical vibration[J].Computers and Ge-otechnics,2008,35(3):489‒499. doi:10.1016/j.compgeo.2007.05.010

[9] Gao Guangyun, Yang Xianjian.Theory and application of r-ow pile vibration isolation[J].Journal of Building Structu-res,1997,18(4):58‒69.

[10] 高广运,杨先健.排桩隔振的理论与应用[J].建筑结构学报,1997,18(4):58‒69.

[11] Huang J, Shi Z.Vibration reduction for plane waves by using periodic pile barriers[C]//International Conference on Railway Engineering.Beijing:Beijing Key Laboratory of Tr-ack Engineering,China Railway Society,2012.

[12] Huang Jiankun, Shi Zhifei.Attenuation zones of periodic pi-le barriers and its application in vibration reduction for pl-ane waves[J].Journal of Sound and Vibration,2013,332(19):4423‒4439. doi:10.1016/j.jsv.2013.03.028

[13] Huang Jiankun, Shi Zhifei.Application of periodic theory to rows of piles for horizontal vibration attenuation[J].International Journal of Geomechanics,2013,13(2):132‒142. doi:10.1061/(asce)gm.1943-5622.0000193

[14] Meng Qingjuan.Research on the vibration reduction characteristics of periodic row piles in saturated soil[D].Beijing:Beijing Jiaotong University,2019.

[15] 孟庆娟.饱和土中周期性排桩的减振特性研究[D].北京:北京交通大学,2019.

[16] Li Jingru, Li Sheng.Multi-flexural wave band gaps of a ne-w periodic metamaterial plate[J].Journal of Vibration and Shock,2018,37(1):163‒171.

[17] 李静茹,黎胜.周期新型超材料板多阶弯曲波带隙研究[J].振动与冲击,2018,37(1):163‒171.

[18] Wu Liangyu, Chen L W.Acoustic band gaps of the woodpile sonic crystal with the simple cubic lattice[J].Journal of Physics D:Applied Physics,2011,44(4):045402. doi:10.1088/0022-3727/44/4/045402

[19] Zhu Xiaohui, Li Longqiu, Zhang Guangyu,et al.Investigation of acoustic manipulation by layered periodic composites[J].Journal of Mechanical Engineering,2017,53(6):10‒15.

[20] 朱晓辉,李隆球,张广玉,等.基于层状周期性结构的声波调控技术研究[J].机械工程学报,2017,53(6):10‒15.

[21] Huang Zigui, Chen Zhengyu.Acoustic waves in two-dim-ensional phononic crystals with reticular geometric stru-ctures[J].Journal of Vibration and Acoustics,2011,133(3):031011. doi:10.1115/1.4003201

[22] Hua Hongliang, Liao Zhenqiang, Zhang Xiangyan.An efficient dynamic modeling method of an axially moving cantilever beam and frequency response analysis[J].Chinese Journal of Theoretical and Applied Mechanics,2017,49(6):1390‒1398.

[23] 华洪良,廖振强,张相炎.轴向移动悬臂梁高效动力学建模及频率响应分析[J].力学学报,2017,49(6):1390‒1398.

[24] Chen Mingfei, Jin Guoyong, Zhang Yantao,et al.Three-di-mensional vibration analysis of beams with axial functionally graded materials and variable thickness[J].Composite Structures,2019,207:304‒322. doi:10.1016/j.compstruct.2018.09.029

[25] Hou Zhilin, Fu Xiujun, Liu Youyan.Singularity of the Blo-ch theorem in the fluid/solid phononic crystal[J].Physical Review B,2006,73(2):024304. doi:10.1103/physrevb.73.024304

[26] Bian Zuguang, Peng Wei, Song Jizhou.Thermal tuning of band structures in a one-dimensional phononic crystal[J].Journal of Applied Mechanics,2014,81(4):041008. doi:10.1115/1.4025058

[27] Cao Yongjun, Hou Zhilin, Liu Youyan.Convergence probl-em of plane-wave expansion method for phononic crystals[J].Physics Letters A,2004,327(2/3):247‒253. doi:10.1016/j.physleta.2004.05.030

[28] Hou Zhilin, Assouar B M.Modeling of Lamb wave propagation in plate with two-dimensional phononic crystal layer coated on uniform substrate using plane-wave-expansion method[J].Physics Letters A,2008,372(12):2091‒2097. doi:10.1016/j.physleta.2007.10.080

[29] Mei Jun, Liu Zhengyou, Shi Jing,et al.Theory for elastic wave scattering by a two-dimensional periodical array of cylinders:An ideal approach for band-structure calculations[J].Physical Review B,2003,67(24):245107. doi:10.1103/physrevb.67.245107

[30] Psarobas I E, Stefanou N, Modinos A.Scattering of elastic waves by periodic arrays of spherical bodies[J].Physical Review B,2000,62(1):278‒291. doi:10.1103/physrevb.62.278

[31] Guo Wenjie, Yang Zhou, Feng Qingsong,et al.A new meth-od for band gap analysis of periodic structures using virtual spring model and energy functional variational principle[J].Mechanical Systems and Signal Processing,2022,168:108634. doi:10.1016/j.ymssp.2021.108634

[32] 石志飞,程志宝,向宏军.周期结构理论及其在隔震减振中的应用[M].北京:科学出版社,2017.

[33] Lu Bin, Chen Yuehua, Feng Zhimin,et al.Bending vibration characteristics analysis of a rectangular plate with complex opening using Chebyshev-variational method[J].Journal of Vibration and Shock,2020,39(2):178‒187.

[34] 陆斌,陈跃华,冯志敏,等.基于Chebyshev变分法的复杂开口形状矩形薄板弯曲振动特性分析[J].振动与冲击,2020,39(2):178‒187.

[35] 温熙森,温激泓,郁殿龙,等.声子晶体[M].北京:国防工业出版社,2009.

[36] Du Jingtao, Li W L, Liu Zhigang,et al.Free vibration of two elastically coupled rectangular plates with uniform elastic boundary restraints[J].Journal of Sound and Vibration,2011,330(4):788‒804. doi:10.1016/j.jsv.2010.08.044

[37] Chen Yuehua, Jin Guoyong, Zhu Minggang,et al.Vibration behaviors of a box-type structure built up by plates and energy transmission through the structure[J].Journal of Sou-nd and Vibration,2012,331(4):849‒867. doi:10.1016/j.jsv.2011.10.002

[38] DAS B.Principles of Geotechnical Engineering SI version[M].Pws,2012.

[39] Briaud J L, Li Yanfeng, Rhee K.BCD:A soil modulus device for compaction control[J].Journal of Geotechnical and Geoenvironmental Engineering,2006,132(1):108‒115. doi:10.1061/(asce)1090-0241(2006)132:1(108)

基金资助

国家自然科学基金项目(12072112)

国家自然科学基金项目(52402425)

中国博士后科学基金项目(2021M700306)

江西省自然科学基金面上项目(20224BAB204069)

江西省自然科学基金面上项目(20242BAB214078)