Pasternak地基中H(t)-V受荷桩水平响应有限杆单元法分析

Analysis of Lateral Response of H（t）-V Loaded Pile in Pasternak Foundation by Finite Beam Element Method

摘要

为探讨桩顶水平动荷载H（t）与竖向荷载V联合作用下桩基的水平响应，基于 Pasternak地基和Euler梁理论，建立了桩-土相互作用水平振动分析模型，采用改进的有限杆单元方法求解考虑P-Δ效应、土体剪切效应影响的综合刚度矩阵方程，结合桩土连续边界条件得到桩身内力与位移解答.通过与已有解析解、有限元解和模型试验的结果比较，验证了计算方法的合理性.最后，对影响桩身内力与位移的主要因素进行分析.结果表明：1）传统Winkler地基相较于Pasternak地基模型，忽略了地基土体的剪切效应，将夸大桩体结构的实际受力，使得计算得到的桩身水平位移和弯矩均大于Pasternak地基所得结果，且随着桩土弹模比E_p/E_s的降低，两种地基模型计算的桩身最大水平位移和弯矩的差异性呈现出增强趋势；2） 随着桩顶竖向荷载的增加，桩身水平位移和弯矩受P-Δ效应的影响显著.当桩顶竖向荷载特征参数λ由0增至2时，桩身最大水平位移和弯矩分别提高40.85%和78.57%；3） 相较无限长桩（L>20d_p），有限长桩的水平位移和弯矩的动力响应受桩身长径比L/d_p影响更大；桩身最大水平位移和弯矩随着水平简谐荷载幅值H₀的增加而增大，随着无量纲频率a₀的增大而减小.

Abstract

To investigate the lateral response of a pile under the combined action of lateral harmonic load H（t） and vertical load V at the top of the pile， based on the Pasternak foundation and Euler beam theory， a lateral vibration analysis model of pile-soil interaction is established. The improved finite beam element method is used to solve the comprehensive stiffness matrix equation considering the P-Δ effect and soil shear effect， and the corresponding analytical solution is obtained by combining the pile-soil continuous boundary conditions. The rationality of the proposed method is verified by comparing its predictions with the results of existing analytical solutions， finite element solutions， and model tests. Finally， the main factors affecting the internal forces and deformations of the pile are analyzed. The results show that： 1） Compared with the Pasternak foundation model， the traditional Winkler foundation ignores the shear effect of the foundation soil， which exaggerates the actual force of the pile so that the calculated lateral displacement and bending moment of the pile are larger than the results obtained by the Pasternak foundation. Additionally， as the pile-soil elastic modulus ratio （E_p/E_s） decreases， the difference in maximum lateral displacement and bending moment between the two foundation models becomes increasingly pronounced； 2） As the vertical load on the pile head increases， the lateral displacement and bending moment of the pile are significantly affected by the P-Δ effect. When the characteristic parameter λ of the vertical load on the pile head increases from 0 to 2， the maximum lateral displacement and bending moment of the pile increase by 40.85% and 78.57%， respectively； 3） Compared with the infinite pile （L>20d_p ）， the dynamic response of the lateral displacement and bending moment of the finite pile is more affected by the slenderness ratio L/d_p； The maximum lateral displacement and bending moment of the pile increase with the increase of the amplitude of the horizontal harmonic load H₀， and decrease with the increase of the dimensionless frequency a₀.

Graphical abstract

关键词

水平动荷载 / 轴向荷载 / 有限单元法 / Pasternak地基 / 动力响应

Key words

horizontal harmonic load / axial loads / finite element method / Pasternak foundation / dynamic response

桩基础作为海上风机、港口码头、桥梁和高速铁路等建设工程中所采用的重要基础型式，其在承受竖向力V的同时，还会承受来自波浪、风、地震等水平动荷载H（t）的作用^［1-3］.在轴、横向静-动荷载作用下，桩基结构的受力变形分析往往较为复杂^［4］，为此探讨复杂动力环境下桩基的水平振动响应，对揭示桩基的受力特性分析具有重要意义.

关于桩基水平振动响应的计算方法，目前常见的有分离变量法^［5］、传递矩阵法^［6］和有限单元法^［7］.其中胡安峰等^［5］和Mylonakis等^［6］基于动力Winkler地基模型，分别利用分离变量法和传递矩阵法计算得到水平简谐荷载作用下均质、层状地基中单桩的水平振动解析解，然而上述计算方法在复杂荷载和复杂边界条件下，数学推导涉及多项积分运算，求解过程较为繁琐，在工程中难以有效推广.有限杆单元法通过构造桩基的整体刚度方程，对复杂荷载和边界条件处理较为简便，计算结果可靠，在桩基础等杆系结构的受力分析中得到了广泛应用^［8-13］.利用有限杆单元法，赵明华等^［8-10］和黎学优等^［11-12］进行了桩顶静力荷载作用下桩基的受力分析；进一步地，蒯行成等^［13］基于动力Winkler地基模型，推导了单桩水平动力阻抗的计算方法；覃玉兰等^［7］基于动力Winkler地基模型求解了桩顶水平简谐荷载H（t）与扭矩T组合荷载作用下单桩的内力与位移.

上述研究普遍将桩周土体简化为Winkler地基模型，用单一参数表示土体的连续性，虽计算简便，但相较双参数Pasternak地基模型忽略了地基土体的剪切效应，计算结果存在一定的误差^［14］.因此，众多学者围绕Pasternak地基在桩基的动力响应中的应用进行了研究.例如，Bhaduri等^［15］和王珏等^［16-17］将地基等效为Pasternak地基，分别采用傅里叶变换和荷载传递矩阵法推导得到单桩的水平振动解析解. Qetin等^［18］探究了均质Pasternak地基弹性参数对桩基自由振动响应的影响.值得注意的是，上述研究大多仅考虑水平动荷载作用对桩基的动力响应.

而在实际工程中，桩基往往承受水平和竖向荷载的共同作用^［2］.对此，Halabe等^［19］求解了竖向静荷载作用下单桩的水平振动微分方程，认为竖向荷载对桩基的水平振动特性有显著的影响.吴鸣等^［4］和熊辉等^［20］分别建立桩-土动力相互作用模型，利用传递矩阵法得到了考虑竖向力影响的水平动力阻抗解析解.栾鲁宝等^［2］通过引入势函数分析了竖向荷载作用下大直径管桩的水平振动响应，得出了管桩水平振动响应频域解析解.可以看出，上述已有桩基水平振动响应研究中，较少有利用有限杆单元方法综合考虑地基剪切变形和桩顶竖向荷载影响的相关研究.

鉴于此，本文结合有限杆单元方法，同时考虑桩周地基土体的剪切效应和桩身P-Δ效应，构建基于Pasternak地基模型的H（t）-V受荷桩水平响应有限单元数值解.通过与已有解析解、有限元解和模型试验进行对比验证方法的正确性.最后，探讨土体剪切效应、桩顶竖向荷载、水平简谐荷载幅值、无量纲频率和桩身长径比等参数对桩身位移和弯矩的影响规律，可为工程实践提供参考.

1 计算模型与基本假设

1.1 基桩受力简化分析模型

单桩受力简化分析模型如图1所示，桩顶承受水平简谐荷载H（t）=H₀e^iωt，P₀为竖向静荷载.其中，H₀为激振力幅值，虚数\mathrm{i}=\sqrt[]{-\mathrm{1}}，ω为动荷载圆频率， t为时间，桩的长度和直径分别为L和d_p，z为地面以下计算结点深度，k_s、c_s分别为地基土体弹性系数与阻尼系数.

1.2 基本假设

为便于分析，基于上述计算模型的受力特征，引入如下基本假设：

1） 桩体为圆形等截面桩，视为均匀、各向同性的Euler-Bernoulli梁；

2） 桩周土体视为均质土体，简化为Pasternak地基模型以描述桩土相互作用，其对桩的作用等效为一系列沿桩身分布的阻尼器、弹簧和剪切单元.其中阻尼系数c_s和弹性系数k_s与土体性质有关；

3） 桩土体系振动满足小变形条件，处于线弹性范围内，桩-土完全接触，接触面保持位移协调，不发生相对滑移与脱离；

4） 忽略水平简谐荷载作用引起的桩身竖向位移.

2 改进的有限杆单元法

本节将在传统有限杆单元法的基础上，对单元刚度矩阵和结点荷载向量做出修正，建立可综合考虑P-Δ效应、土体剪切效应和桩侧摩阻力等因素的综合刚度矩阵方程，以期用于桩体水平振动响应的数值计算中.

2.1 基本平衡方程

将桩体离散为由有限个结点连接的杆单元，结点受力如图2所示.

深度z处杆单元长度为l，施加在i和j结点上的载荷包括竖向荷载{\overline{P}}_i和{\overline{P}}_j，剪力{\overline{H}}_i和{\overline{H}}_j，弯矩{\overline{M}}_i和{\overline{M}}_j.对应的结点变形分别为竖向位移{\overline{u}}_i和{\overline{u}}_j，横向位移{\overline{v}}_i和{\overline{v}}_j，转角{\overline{\theta }}_i和{\overline{\theta }}_j.根据杆单元虚位移原理，单元刚度方程可表示为^［8］：

\begin{array}{c}{\mathit{F}}^{\mathrm{e}}={\left\{\begin{array}{cccccc}{\overline{P}}_i& {\overline{H}}_i& {\overline{M}}_i& {\overline{P}}_j& {\overline{H}}_j& {\overline{M}}_j\end{array}\right\}}^{\mathrm{T}}=\\ {{\mathit{k}}_{\mathrm{0}}}^{\mathrm{e}}{\left\{\begin{array}{cccccc}{\overline{u}}_i& {\overline{v}}_i& {\overline{\theta }}_i& {\overline{u}}_j& {\overline{v}}_j& {\overline{\theta }}_j\end{array}\right\}}^{\mathrm{T}}={{\mathit{k}}_{\mathrm{0}}}^{\mathrm{e}}{\mathit{\delta }}^{\mathrm{e}}\end{array}（1）

式中： F^e为单元结点荷载向量； δ^e为杆端位移向量； k₀^e为单元刚度矩阵.

2.2 考虑P-Δ效应的单元附加刚度矩阵

在水平动荷载H（t）作用下，竖向荷载因桩身的水平偏移而产生附加弯矩，将进一步加剧桩身的挠曲变形，即桩身产生P-Δ效应.故本节将对上节传统杆单元刚度矩阵{{\mathit{k}}_{\mathrm{0}}}^{\mathrm{e}}进行修正，以考虑竖向荷载对桩身弯矩和水平位移的影响.根据杆单元虚位移原理，考虑P-Δ效应的单元刚度方程^［9］可表示为：

\begin{array}{c}{{\mathit{F}}_{\mathrm{N}}}^{\mathrm{e}}={{\mathit{k}}_{\mathrm{N}}}^{\mathrm{e}}{\left\{\begin{array}{cccccc}{\overline{u}}_i& {\overline{v}}_i& {\overline{\theta }}_i& {\overline{u}}_j& {\overline{v}}_j& {\overline{\theta }}_j\end{array}\right\}}^{\mathrm{T}}\\ ={\int }_{\mathrm{0}}^l[P_{\mathrm{N}}(\overline{z}){\mathit{N}}_{\mathrm{v}}^{\mathrm{T}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{v}}}{\mathrm{d}{\overline{z}}^{\mathrm{2}}}+P_{\mathrm{N}}(\overline{z})\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{v}}^{\mathrm{T}}}{\mathrm{d}{\overline{z}}^{\mathrm{2}}}{\mathit{N}}_{\mathrm{v}}]{\mathit{\delta }}^{\mathrm{e}}\mathrm{d}\overline{z}\end{array}（2）

式中： F_N^e为单元等效结点荷载向量； δ^e为杆端位移向量；P_{\mathrm{N}}(\overline{z})为距i单元上结点深度\overline{z}处的竖向荷载，可表示为P_{\mathrm{N}}(\overline{z})=P_{\mathrm{0}}+f_{\mathrm{s}}\cdot (h_i+\overline{z})，h_i 为i结点距桩顶的距离，f_s为桩侧摩阻力； N_v表示与杆单元水平位移对应的形函数矩阵，可表示为：

{\mathit{N}}_{\mathrm{v}}=\left[\begin{array}{cccccc}\mathrm{0}& N_{\mathrm{1}}& N_{\mathrm{2}}& \mathrm{0}& N_{\mathrm{3}}& N_{\mathrm{4}}\end{array}\right]（3）

其中：

N_{\mathrm{1}}=\mathrm{1}-\frac{\mathrm{3}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{2}}}+\frac{\mathrm{2}{\overline{z}}^{\mathrm{3}}}{l^{\mathrm{3}}},N_{\mathrm{2}}=\overline{z}-\frac{\mathrm{2}{\overline{z}}^{\mathrm{2}}}{l}+\frac{{\overline{z}}^{\mathrm{3}}}{l^{\mathrm{2}}}

N_{\mathrm{3}}=\frac{\mathrm{3}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{2}}}-\frac{\mathrm{2}{\overline{z}}^{\mathrm{3}}}{l^{\mathrm{3}}},N_{\mathrm{4}}=-\frac{{\overline{z}}^{\mathrm{2}}}{l}+\frac{{\overline{z}}^{\mathrm{3}}}{l^{\mathrm{2}}}（4）

k_N^e为单元刚度矩阵，可表示为：

{{\mathit{k}}_{\mathrm{N}}}^{\mathrm{e}}=\left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\frac{\mathrm{12}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{12}{h}_i+l)}{\mathrm{5}l}& -\frac{\mathrm{6}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{6}{h}_i+l)}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{12}{P}_{\mathrm{0}}+\mathrm{6}{f}_{\mathrm{s}}\cdot (\mathrm{2}{h}_i+l)}{\mathrm{5}l}& -\frac{P_{\mathrm{0}}+f_{\mathrm{s}}\cdot h_i}{\mathrm{5}}\\ \mathrm{0}& -\frac{\mathrm{6}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{6}{h}_i+l)}{\mathrm{5}}& -\frac{\mathrm{4}{P}_{\mathrm{0}}l+f_{\mathrm{s}}\cdot l(\mathrm{4}{h}_i+l)}{\mathrm{15}}& \mathrm{0}& \frac{P_{\mathrm{0}}+f_{\mathrm{s}}\cdot (h_i+l)}{\mathrm{5}}& \frac{\mathrm{2}{P}_{\mathrm{0}}l+f_{\mathrm{s}}\cdot l(\mathrm{2}{h}_i+l)}{\mathrm{30}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{12}{P}_{\mathrm{0}}+\mathrm{6}{f}_{\mathrm{s}}\cdot (\mathrm{2}{h}_i+l)}{\mathrm{5}l}& \frac{P_{\mathrm{0}}+f_{\mathrm{s}}\cdot (h_i+l)}{\mathrm{5}}& \mathrm{0}& -\frac{\mathrm{12}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{12}{h}_i+\mathrm{11}l)}{\mathrm{5}l}& \frac{\mathrm{6}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{6}{h}_i+\mathrm{5}l)}{\mathrm{5}}\\ \mathrm{0}& -\frac{P_{\mathrm{0}}+f_{\mathrm{s}}\cdot h_i}{\mathrm{5}}& \frac{\mathrm{2}{P}_{\mathrm{0}}l+f_{\mathrm{s}}\cdot l(\mathrm{2}{h}_i+l)}{\mathrm{30}}& \mathrm{0}& \frac{\mathrm{6}{P}_{\mathrm{0}}+f_{\mathrm{s}}\cdot (\mathrm{6}{h}_i+\mathrm{5}l)}{\mathrm{5}}& -\frac{\mathrm{4}{P}_{\mathrm{0}}l+f_{\mathrm{s}}\cdot l(\mathrm{4}{h}_i+\mathrm{3}l)}{\mathrm{15}}\end{array}\right]（5）

2.3 考虑双参数地基模型的单元附加刚度矩阵

Pasternak双参数弹性地基模型相较单参数Winkler地基模型，考虑了地基土体剪切的效应，计算结果更精确^［21］.基于Pasternak双参数弹性地基模型的桩身挠曲微分方程为^［22］：

E_{\mathrm{p}}{I}_{\mathrm{p}}\frac{{\mathrm{d}}^{\mathrm{4}}x}{\mathrm{d}{z}^{\mathrm{4}}}-G_{\mathrm{s}}{b}_{\mathrm{0}}\frac{{\mathrm{d}}^{\mathrm{2}}x}{\mathrm{d}{z}^{\mathrm{2}}}+Kx=\mathrm{0}（6）

式中：E_p、I_p、b₀分别为桩身弹性模量、截面惯性矩和计算宽度；x为桩身挠度；G_s为土体剪切模量；K为地基反力模量.

首先基于动力Winkler地基模型，桩周土产生沿桩身连续分布的水平抗力p（x），其与位移的关系式^［7］可表示为：

p(x)=k\cdot x(\mathrm{z})=(k_{\mathrm{s}}-{\omega }^{\mathrm{2}}{m}_{\mathrm{p}}+\mathrm{i}\omega c_{\mathrm{s}})\cdot x(\mathrm{z})（7）

式中：p（x）、m_p、ω分别为桩周土对桩体的水平抗力、单位桩长的质量和谐振频率.其中，k_s、c_s取值由下式^［23］确定：

\begin{array}{l}{k}_{\mathrm{s}}=\delta E_{\mathrm{s}}\\ c_{\mathrm{s}}=\mathrm{6}{a}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{4}}{\rho }_{\mathrm{s}}{V}_{\mathrm{s}}{d}_{\mathrm{p}}+\mathrm{2}\frac{k_{\mathrm{s}}{\beta }_{\mathrm{s}}}{\omega }\end{array}（8）

式中：δ为与土体参数、荷载类型、桩土相对刚度和桩顶约束条件有关的函数，与深度z无关.当桩顶自由时，δ≈1.5~2.5；桩顶约束转角时，δ≈1.0~1.2.本文考虑桩顶自由时k_{\mathrm{s}}=\mathrm{2.1}{E}_{\mathrm{s}}，桩顶约束转角时k_{\mathrm{s}}=\mathrm{1.2}{E}_{\mathrm{s}}.V_{\mathrm{s}}=\sqrt[]{E_{\mathrm{s}}/\mathrm{2}{\rho }_{\mathrm{s}}(\mathrm{1}+v_{\mathrm{s}})}，为土层的剪切波速，E_s、ρ_s、v_s和β_s分别为桩周土的弹性模量、密度、泊松比和阻尼比，a₀=ωd_p/V_s，为无量纲频率.

在此基础上，由式（6）知，当G_s=0时，Pasternak双参数地基模型即退化为Winkler地基模型.当无G_s实测值时，可由Tanahashi^［24］提出的经验公式确定：

G_{\mathrm{s}}=\frac{E_{\mathrm{s}}{t}_{\mathrm{s}}}{\mathrm{6}(\mathrm{1}+v_{\mathrm{s}})}（9）

式中：E_s、v_s同式（8）；t_s为剪切层厚度.

梁发云等^［25］基于Yao等^［26］的数值模拟结果，认为地基土剪切层厚度t_s近似等于水平荷载作用对桩侧土的影响范围S（S=11d_p，d_p为桩径），故本文分析中取t_s=11d_p，实际上剪切层厚度t_s与桩周土体的性质有关，本文对此不做深入研究.作用在杆单元上的荷载按照其作用位置不同，可分为结点荷载和非结点荷载两种.为建立各结点的平衡方程，需将分布作用在杆单元上的桩周土抗力（非结点荷载）按照虚位移原理等效为结点荷载 F^e，可表示为：

\begin{array}{l}{{\mathit{F}}_{\mathrm{H}}}^{\mathrm{e}}=k{\int }_{\mathrm{0}}^l\left[{{\mathit{N}}_{\mathrm{v}}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{v}}{\mathit{\delta }}^{\mathrm{e}}\right]\mathrm{d}\overline{z}=\\ (k_{\mathrm{s}}-{\omega }^{\mathrm{2}}{m}_{\mathrm{p}}+\mathrm{i}\omega c_{\mathrm{s}}){\int }_{\mathrm{0}}^l\left[{{\mathit{N}}_{\mathrm{v}}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{v}}{\mathit{\delta }}^{\mathrm{e}}\right]\mathrm{d}\overline{z}=\\ \left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{13}kl}{\mathrm{35}}& \frac{\mathrm{11}k{l}^{\mathrm{2}}}{\mathrm{210}}& \mathrm{0}& \frac{\mathrm{9}kl}{\mathrm{70}}& -\frac{\mathrm{13}k{l}^{\mathrm{2}}}{\mathrm{420}}\\ \mathrm{0}& \frac{\mathrm{11}k{l}^{\mathrm{2}}}{\mathrm{210}}& \frac{k{l}^{\mathrm{3}}}{\mathrm{105}}& \mathrm{0}& \frac{\mathrm{13}k{l}^{\mathrm{2}}}{\mathrm{420}}& -\frac{k{l}^{\mathrm{3}}}{\mathrm{140}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{9}kl}{\mathrm{70}}& \frac{\mathrm{13}k{l}^{\mathrm{2}}}{\mathrm{420}}& \mathrm{0}& \frac{\mathrm{13}kl}{\mathrm{35}}& -\frac{\mathrm{11}k{l}^{\mathrm{2}}}{\mathrm{210}}\\ \mathrm{0}& -\frac{\mathrm{13}k{l}^{\mathrm{2}}}{\mathrm{420}}& -\frac{k{l}^{\mathrm{3}}}{\mathrm{140}}& \mathrm{0}& -\frac{\mathrm{11}k{l}^{\mathrm{2}}}{\mathrm{210}}& \frac{k{l}^{\mathrm{3}}}{\mathrm{105}}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_i\\ {\overline{v}}_i\\ {\overline{\theta }}_i\\ {\overline{u}}_j\\ {\overline{v}}_j\\ {\overline{\theta }}_j\end{array}\right\}=\\ {{\mathit{k}}_{\mathrm{1}}}^{\mathrm{e}}{\left\{\begin{array}{cccccc}{\overline{u}}_i& {\overline{v}}_i& {\overline{\theta }}_i& {\overline{u}}_j& {\overline{v}}_j& {\overline{\theta }}_j\end{array}\right\}}^{\mathrm{T}}\end{array}（10）

式中： k₁^e为考虑桩周土水平反力作用的单元附加刚度矩阵；k=(k_{\mathrm{s}}-{\omega }^{\mathrm{2}}{m}_{\mathrm{p}}+\mathrm{i}\omega c_{\mathrm{s}}).

\begin{array}{l}{{\mathit{F}}_{\mathrm{Q}}}^{\mathrm{e}}=G_{\mathrm{s}}{b}_{\mathrm{0}}{\int }_{\mathrm{0}}^l\left[{{\mathit{{\rm N}}}_{\mathrm{\phi }}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{\phi }}{\mathit{\delta }}^{\mathrm{e}}\right]\mathrm{d}\overline{z}=\\ \left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{6}{G}_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{5}l}& \frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& \mathrm{0}& -\frac{\mathrm{6}{G}_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{5}l}& \frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}\\ \mathrm{0}& \frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& \frac{\mathrm{2}{G}_{\mathrm{s}}{b}_{\mathrm{0}}l}{\mathrm{15}}& \mathrm{0}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}l}{\mathrm{30}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\frac{\mathrm{6}{G}_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{5}l}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& \mathrm{0}& \frac{\mathrm{6}{G}_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{5}l}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}\\ \mathrm{0}& \frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}l}{\mathrm{30}}& \mathrm{0}& -\frac{G_{\mathrm{s}}{b}_{\mathrm{0}}}{\mathrm{10}}& \frac{\mathrm{2}{G}_{\mathrm{s}}{b}_{\mathrm{0}}l}{\mathrm{15}}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_i\\ {\overline{v}}_i\\ {\overline{\theta }}_i\\ {\overline{u}}_j\\ {\overline{v}}_j\\ {\overline{\theta }}_j\end{array}\right\}=\\ {{\mathit{k}}_{\mathrm{2}}}^{\mathrm{e}}{\left\{\begin{array}{cccccc}{\overline{u}}_i& {\overline{v}}_i& {\overline{\theta }}_i& {\overline{u}}_j& {\overline{v}}_j& {\overline{\theta }}_j\end{array}\right\}}^{\mathrm{T}}\end{array}（11）

式中： k₂^e为考虑地基土等效弹簧间剪切力作用的单元附加刚度矩阵； N_φ表示与杆单元转角位移对应的形函数矩阵，可表示为

{\mathit{N}}_{\mathrm{\phi }}=\left[\begin{array}{cccccc}\mathrm{0}& N_{\mathrm{5}}& N_{\mathrm{6}}& \mathrm{0}& N_{\mathrm{7}}& N_{\mathrm{8}}\end{array}\right]（12）

其中，

N_{\mathrm{5}}=-\frac{\mathrm{6}\overline{z}}{l^{\mathrm{2}}}+\frac{\mathrm{6}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{3}}},N_{\mathrm{6}}=\mathrm{1}-\frac{\mathrm{4}\overline{z}}{l}+\frac{\mathrm{3}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{2}}}

N_{\mathrm{7}}=\frac{\mathrm{6}\overline{z}}{l^{\mathrm{2}}}-\frac{\mathrm{6}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{3}}},N_{\mathrm{8}}=-\frac{\mathrm{2}\overline{z}}{l}+\frac{\mathrm{3}{\overline{z}}^{\mathrm{2}}}{l^{\mathrm{2}}}（13）

2.4 考虑桩侧摩阻力的等效结点荷载向量

桩基础受外荷载作用后，将相对桩侧土体发生一定程度的偏转与沉降.如图3所示，桩表面单元m处的竖向变形量为z_{\mathrm{m}}=d_{\mathrm{p}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \cdot \theta (\overline{z})/\mathrm{2}，单元面积为\mathrm{\delta }{A}_{\mathrm{m}}=z_{\mathrm{m}}\mathrm{\delta }\alpha \cdot d_{\mathrm{p}}/\mathrm{2}，其中径向夹角为α，水平面偏转角为\theta (\overline{z}).基于虚位移原理，桩侧摩阻力所做的外力虚功δW_τ的关系式^［9］为：

\begin{array}{l}\mathrm{\delta }{W}_{\mathrm{\tau }}={\int }_{\mathrm{0}}^l\tau {d_p}^{\mathrm{2}}\delta \theta (\overline{z})d\overline{z}=\tau {d_{\mathrm{p}}}^{\mathrm{2}}{\int }_{\mathrm{0}}^l{\mathit{N}}_{\phi }\delta {\mathit{\delta }}^{e}d\overline{z}=\\ \tau {d_{\mathrm{p}}}^{\mathrm{2}}{\int }_{\mathrm{0}}^l\frac{\mathrm{d}{\mathit{N}}_{\mathrm{v}}}{\mathrm{d}\overline{z}}\delta {\mathit{\delta }}^{e}d\overline{z}={\mathit{F}}_{\mathrm{\tau }}^{\mathrm{e}}\mathrm{\delta }{\mathit{\delta }}^{\mathrm{e}}\end{array}（14）

式中，{\mathit{F}}_{\mathrm{\tau }}^{\mathrm{e}}为桩侧摩阻力的等效结点荷载向量， N_v、 N_φ意义同上.可得，

{\mathit{F}}_{\mathrm{\tau }}^{\mathrm{e}}={\left[\begin{array}{cccccc}\mathrm{0}& -\tau {d_{\mathrm{p}}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \tau {d_{\mathrm{p}}}^{\mathrm{2}}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}（15）

综上所述，根据式（10）和式（11）可确定考虑 土体剪切效应的单元刚度矩阵 k_p^e，再结合式（1）、 式（5）、和式（15）可得到考虑桩身单元刚度、P-Δ效应、土体剪切效应和桩侧摩阻力等因素的单元综合刚度方程为：

{\mathit{F}}^{\mathrm{e}}={\mathit{K}}^{\mathrm{e}}{\mathit{\delta }}^{\mathrm{e}}=({{\mathit{k}}_{\mathrm{0}}}^{\mathrm{e}}+{{\mathit{k}}_{\mathrm{N}}}^{\mathrm{e}}+{{\mathit{k}}_{\mathrm{p}}}^{\mathrm{e}}){\mathit{\delta }}^{\mathrm{e}}（16）

3 有限杆单元法求解内力、位移步骤

基于改进的有限杆单元方法，建立考虑上述因素的综合刚度矩阵，编制MATLAB程序求解受荷桩基的内力与位移，主要包括如下计算步骤：

1） 单元离散化：将桩基离散为若干具有较高计算精度的微单元和结点.其中结点位置可取集中力作用点、土层分界处、桩截面面积突变点、支承点等.

2） 整体特性分析：结合式（1）、（5）、（10）和

式（11），按照后处理法将单元刚度矩阵组装为整体刚度矩阵，并建立整体刚度方程.

3） 结点位移求解：引入桩基的位移边界条件，求解整体刚度方程，得到桩身任意结点的位移向量 δ^e.

4） 杆端力求解：将结点的位移向量代入单元刚度方程式（1）中，求得杆单元的杆端力，进而确定桩身内力与水平变形量，绘制桩身水平位移与内力图.

4 算例验证

4.1 算例1

为验证本文解的合理性，将本文与文献［27］考虑竖向荷载影响的水平振动响应室内模型试验解、有限元解进行对比分析.计算参数如下：桩长L= 0.9 m，桩径d_p=0.06 m，泊松比v_p=0.3；桩周土的密度ρ_s=1 900 kg/m³，弹性模量E_s=8.5 MPa，泊松比v_s=0.3，阻尼比β_s=0.05；第一组桩顶组合荷载为f =3 Hz，水平动荷载幅值H₀=500 N，竖向荷载P₀=0 N；第二组桩顶组合荷载为f =2 Hz，水平动荷载幅值H₀=150.2 N，竖向荷载P₀=150 N；边界条件为桩顶、桩底自由.图4和图5分别为两组桩在桩顶受不同组合荷载作用一个周期内桩身弯矩随深度z的变化曲线.

由图4及图5可知，本文方法所得的桩身弯矩计算值与文献［27］的有限元值较为接近，两组桩身弯矩最大值均发生在距桩顶1/3深度处，桩身弯矩随深度的变化趋势一致.

4.2 算例2

为进一步验证本文方法的正确性，将本文解与文献［7］的水平振动响应解进行对比.令P₀=0，即退化到水平简谐荷载H（t）作用下的水平振动响应解，并将本文退化解与已有案例解进行对比.相关计算参数为：桩长L=4.5 m，桩径d_p=0.3 m，桩身弹性模量E_p= 20 GPa，泊松比v_p=0.2；桩周土的密度ρ_s=2 000 kg/m³，弹性模量E_s=4 MPa，泊松比v_s=0.4，阻尼比β_s=0.05；边界条件为桩顶自由，桩底嵌固.图6和图7给出了水平动荷载幅值H₀=100 kN，无量纲频率a₀=0.5作用一个周期内桩身位移、弯矩的最大值随深度z的变化曲线.

由图6、图7可知，当水平动荷载幅值H₀=100 kN时，桩顶最大位移出现在桩顶位置，最大值为 14.12 mm；最大桩身弯矩发生在距桩顶1/3深度处，最大值为42.96 kN·m，与文献［7］计算结果吻合较好.文献［7］基于Winkler地基模型，忽略地基土体的剪切效应，使得桩体结构的受力计算结果较实际偏大.结合图6、图7的对比验证结果，表明利用本文计算方法得到的桩身内力和位移与文献结果基本一致，验证了提出方法的准确性.

4.3 算例3

文献［28］基于Winkler地基理论，得到了层状地基中单桩水平振动解析解，为进一步验证本文解的适用性，本文Pasternak地基可退化为Winkler地基模型（G_s→0）.计算参数如下：桩径d_p=0.3 m，桩长径比L/d_p=10，15，泊松比v_p=0.25，桩体密度ρ_p=2 500 kg/m³，弹性模量E_p=20 GPa；土层参数n=2，两层土的密度、阻尼比和泊松比取相同值，分别为ρ_s=2 000 kg/m³、β_s=0.05和v_s=0.4，桩土刚度比E_p/E_s2=5 000，上下两层土的弹性模量比E_s1/E_s2=0.5；边界条件为桩顶约束转角，桩底固定.图8、图9为水平动荷载幅值H₀=100 kN，无量纲频率a₀=0.5，竖向荷载P₀=0，作用一个循环周期内桩身水平位移、弯矩的最大值随z/d_p的变化曲线.为便于分析，对深度z处桩身水平位移和弯矩的最大值x_max（z）和M_max（z）进行无量纲化处理，其中\overline{x}=E_{\mathrm{p}}\cdot d_{\mathrm{p}}\cdot x_{\mathrm{m}\mathrm{a}\mathrm{x}}(z)/(\mathrm{500}{H}_{\mathrm{0}})和\overline{M}=M_{\mathrm{m}\mathrm{a}\mathrm{x}}(z)/(\mathrm{2}{H}_{\mathrm{0}}\cdot d_{\mathrm{p}}).

由图8、图9可以看出，不同桩身长径比条件下，本文得到的桩身无量纲水平位移和弯矩与文献［28］已有解析解结果吻合较好，进一步验证本文方法的适用性.

结合上述算例的对比验证结果，可见当单元数划分合适，本文基于综合刚度矩阵的改进有限杆单元法对于复杂边界、荷载和层状地基等工况下得到的桩身内力和位移解与已有文献的解析解、有限元解和试验结果基本一致，计算精度较高，验证了提出方法的准确性和适应性.

5 参数分析

本节将侧重分析土体剪切效应、桩顶竖向荷载P₀、水平简谐荷载幅值H₀、无量纲频率a₀和桩身长径比L/d_p等参数对桩身水平位移、弯矩的影响.相关计算参数选取为：L/d_p=20，桩径d_p=0.3 m，桩身弹性模量E_p=31.5 GPa，泊松比v_p=0.2；桩周土的密度ρ_s= 2 000 kg/m³，弹性模量E_s=15 MPa，泊松比v_s=0.3，阻尼比β_s=0.05，边界条件为桩顶自由、桩底固定.水平动荷载幅值H₀=60 kN，无量纲频率a₀=0.5，竖向荷载P₀=200 kN.为便于对比分析，对深度z处桩身水平位移和弯矩的最大值x_max（z）和M_max（z）进行无量纲化处理，其中x(z)=x_{\mathrm{m}\mathrm{a}\mathrm{x}}(z)/x_{\mathrm{m}\mathrm{a}\mathrm{x}}和M(z)=M_{\mathrm{m}\mathrm{a}\mathrm{x}}(z)/M_{\mathrm{m}\mathrm{a}\mathrm{x}}.

5.1 土体剪切效应的影响

为探讨土体剪切效应对桩身内力和变形的影响规律，图10和图11采用Pasternak地基模型与传统Winkler地基模型（G_s→0）计算了不同桩土弹模比E_p/E_s和土体剪切模量G_s下单桩的无量纲水平位移与弯矩的变化情况，其中固定E_s=15 MPa，E_p/E_s分别取 1 500、2 500和3 500^［7］；G_s =0、1.0G_s和5.0G_s.由图10、图11可以看出：1）地基土体的剪切效应对约束桩身的内力和变形具有一定的贡献.当G_s =0和1.0 G_s时，采用Winkler地基计算得到的桩身水平位移和弯矩大于Pasternak地基所得结果，说明未引入土体剪切模量G_s将忽略地基土体的剪切效应，使得桩体结构的受力计算结果较实际偏大，因此，土体剪切效应对桩身水平振动响应的影响不容忽视.2） 土体剪切模量G_s对桩身水平位移与弯矩具有削弱作用.当E_p/E_s=2 500，G_s由0增至5G_s时，桩身最大水平位移由0.856减至0.510，减小40.42%，最大弯矩由0.960减至0.343，减小64.27%.3） Pasternak地基和Winkler地基计算所得的桩身最大水平位移和弯矩的差异性随着桩土弹模比E_p/E_s的降低呈现出增强趋势.当土体剪切模量G_s =0和1.0G_s，E_p/E_s由3 500降至1 500时，桩顶水平位移逐渐增加，变化幅度由9.79%增至13.33%；桩身最大弯矩逐渐降低，变化幅度由20.42%增至34.92%.故当桩土弹模比较小时，土体剪切效应对桩身水平振动响应的影响不容忽视.

5.2 桩顶竖向荷载P₀的影响

梁发云等^［29］定义桩顶荷载H₀ = πG_sd_p²和P₀ = λH₀，探究不同竖向荷载对桩身内力和变形的影响.图12和图13给出了竖向荷载特征参数λ取0、0.5、1和2时桩身无量纲水平位移、弯矩随深度的变化.由图12、图13可知：1） 随着桩顶竖向荷载特征参数λ的逐渐增加，桩身无量纲水平位移和弯矩均显著增大.当λ由0增至2时，桩顶水平位移由0.71增至1，增加了40.85%；桩身最大弯矩值由0.56增至1，提高了78.57%.2） 竖向荷载特征参数λ的变化在较浅深度范围内（z=0~0.4L）对桩身无量纲位移影响显著，而随着深度的增加，这种影响逐渐减缓，下半段桩身位移曲线基本重合.桩身无量纲弯矩在距桩顶L/3深度位置出现反转拐点，而桩身最大位移和最大弯矩发生位置受竖向荷载特征参数λ的影响较小.上述计算实例进一步表明，桩顶竖向荷载和水平荷载引起的P-Δ效应会影响桩身最大水平位移和弯矩，可能威胁结构的安全，在工程设计计算中应引起重视.

5.3 水平简谐荷载幅值H₀的影响

赵明华等^［30］结合水平受荷桩模型试验，认为桩顶水平位移x=0.15d_p时的荷载近似等于桩基的水平极限承载力H_u.结合算例，求得其水平极限承载力为H_u=1 100 kN.为探究不同水平简谐荷载幅值H₀对桩身内力和变形的影响，图14和图15给出了水平简谐荷载幅值H₀分别为0.2H_u、0.4H_u、0.8H_u和1.0H_u时的桩身无量纲水平位移、弯矩随深度的变化，其中P₀=100 kN，其他计算参数同上.由图14、图15可知：1） 随着桩顶水平简谐荷载幅值H₀的逐渐增加，桩身无量纲水平位移和弯矩均显著增大.当H₀由0.2H₀增至1.0H₀时，桩身最大无量纲位移和弯矩值均由0.2增至1，提高了4倍.2） 桩身位移随深度的增加呈现减小趋势，且桩身最大位移出现在桩顶；桩身弯矩随深度的增加呈现出先增加后减小趋势，最大弯矩发生位置在距桩顶2.8d_p处.而桩身最大位移和最大弯矩发生位置受水平简谐荷载幅值H₀的影响较小.

5.4 水平简谐荷载无量纲频率a₀的影响

图16、图17给出了桩身长径比L/d_p为10和30，无量纲激振频率a₀=0.01、0.5和1.0时^［28］桩身无量纲水平位移、弯矩随深度的变化.由图16和图17可知：1）在不同桩身长径比的情况下，随着a₀降低，桩身最大水平位移、弯矩均逐渐增大.其中，当L/d_p=10和30，a₀由1.0降至0.01时，桩身最大水平位移分别提高105.89%和141.67%，桩身最大弯矩分别增加15.07%和22.45%.表明桩身水平位移和弯矩受无量纲频率a₀的影响较大，这是由于频率越高，桩的变形还未来得及发生，荷载已反向施加^［28］，从而造成位移量和弯矩减小，因此低频振动环境对桩身内力和位移的影响更大.2） 不同长径比下的桩身最大无量纲位移、弯矩发生的位置与无量纲频率a₀基本无关，最大位移均位于桩顶，最大弯矩均在距地面z=3d_p深度处.

5.5 长径比L/d_p的影响

图18和图19分析了桩身长径比L/d_p分别取5、10、20和40四种情况下^［7］，桩身无量纲水平位移与弯矩随深度的变化曲线，其中纵坐标归一化深度取z与最大桩长L_max的比值（L_max取40d_p）.从图18、图19中可以看出：1） 当桩身长径比较小时（L/d_p=5），桩身水平位移随深度增加而减小，而桩身弯矩随深度增加而增大，且桩底弯矩最大，桩呈刚性短桩特性.2） 随着长径比的增大（L/d_p=10，20，40），桩身出现挠曲变形，水平位移呈曲线，桩表现为弹性长桩特性，桩身弯矩随深度先增加后减小，最大弯矩均出现在z/L_max=0.07~0.08处，且最大弯矩位置不随L/d_p改变而变化.3） 当z/L_max≥0.5（L_p=20d_p），桩身的位移和弯矩均接近于0，故可将z=20d_p=L_p视为桩的有效长度，若再增加桩长（L_p=20d_p，40d_p），桩身水平位移、弯矩曲线已基本重合，不再随桩长的增加而变化.

6 结 论

本文基于Pasternak地基和Euler梁理论，建立了桩-土相互作用水平振动简化分析模型，考虑了竖向荷载和水平简谐荷载引起的P-Δ效应、土体剪切效应的影响，采用改进的有限杆单元方法求解综合刚度矩阵方程，结合桩土连续边界条件得到了桩身内力和位移解答.并对土体剪切效应、桩顶竖向荷载P₀、水平简谐荷载幅值H₀、无量纲频率a₀和桩身长径比L/d_p进行参数分析，得到如下结论：

1） 本文所采用的有限杆单元法通过构造桩基的整体刚度方程可综合考虑地基土体剪切效应和竖向荷载变化的影响，所得结果均与已有文献解进行对比验证.相较传统求解方法可以更简便地考虑复杂荷载、地层和边界条件的影响，计算效率高，可为进一步分析复杂工况下桩基的水平振动响应分析和设计提供参考.

2） 采用Pasternak地基模型计算得到的桩身水平位移和弯矩均小于Winkler地基所得结果，说明未引入土体剪切模量G_s将忽略地基土体的剪切效应，使得桩体结构的受力计算结果较实际偏大.当G_s由0增至5G_s时，桩身最大水平位移和弯矩分别减小40.42%和64.27%.且随着桩土弹模比E_p/E_s的降低，两种地基模型计算的桩身最大水平位移和弯矩的差异性呈现出增强趋势.

3） 随着桩顶竖向荷载的增加，桩身水平位移和弯矩受P-Δ效应的影响显著.随着桩顶竖向荷载特征参数λ的增大，桩身最大水平位移和弯矩均显著增加，其中当λ由0增至2时，桩身最大无量纲水平位移由0.71增至1，增加40.85%；桩身最大无量纲弯矩由0.56增至1，提高78.57%.因此，在工程设计计算中应重视P-Δ效应对桩身内力和变形的影响，避免威胁结构的安全.

4） 整体上，桩身长径比L/d_p对有限长桩的位移、弯矩的动力响应影响较大；当桩长L>20d_p后（无限长桩），影响逐渐减小.且桩基础的水平动力响应随着水平简谐荷载幅值H₀的增加而增加，随着无量纲频率a₀的降低而增加.其中H₀由0.2H₀增至1.0H₀时，桩身最大无量纲水平位移和弯矩均提高了4倍；当L/d_p=30，a₀由1.0降至0.01时，桩身最大无量纲水平位移和弯矩分别提高了141.67%和22.45%.故在工程设计中应注重水平简谐荷载高幅值、低频率振动环境对桩身变形和内力的影响.

参考文献

[1] 江杰， 付臣志， 王顺苇， 等 ．考虑实际分布形式的水平受荷桩桩周土抗力分析方法［J］．工程力学， 2021， 38（11）： 199-211．

[2] JIANG J， FU C Z， WANG S W，et al ．Analytical method of soil resistance around laterally loaded piles considering its actual distribution［J］．Engineering Mechanics， 2021， 38（11）： 199-211．（in Chinese）

[3] 栾鲁宝， 丁选明， 刘汉龙， 等 ．考虑竖向荷载影响的大直径管桩水平振动响应解析解［J］．岩土工程学报， 2016， 38（10）：1859-1868．

[4] LUAN L B， DING X M， LIU H L，et al ．Analytical solution of lateral dynamic response of a large diameter pipe pile considering influence of axial load［J］．Chinese Journal of Geotechnical Engineering， 2016， 38（10）： 1859-1868．（in Chinese）

[5] 任宇， 朱斌，陈仁朋，等 ．大型风电机组群桩基础受荷特性及长期累积沉降控制［J］．土木工程学报， 2010， 43（3）：75-80．

[6] REN Y， ZHU B， CHEN R P， et al ．Performance and control of long-term settlement of pile group foundations for large wind turbines［J］．China Civil Engineering Journal，2010， 43（3）： 75-80．（in Chinese）

[7] 吴鸣， 赵明华， 邹银生 ．层状地基中轴、横向荷载下桥梁基桩的水平动力反应［J］．土木工程学报， 2005， 38（12）： 88-93．

[8] WU M， ZHAO M H， ZOU Y S ．Horizontal dynamic response of piles under simultaneous axial and lateral loads in layered soils［J］．China Civil Engineering Journal，2005， 38（12）： 88-93．（in Chinese）

[9] 胡安峰， 谢康和， 应宏伟，等 ．粘弹性地基中考虑桩体剪切变形的单桩水平振动解析理论［J］．岩石力学与工程学报， 2004， 23（9）： 1515-1520．

[10] HU A F， XIE K H， YING H W，et al ．Analytical theory of lateral vibration of single pile in visco-elastic subgrade considering shear deformation［J］．Chinese Journal of Rock Mechanics and Engineering， 2004， 23（9）： 1515-1520．（in Chinese）

[11] MYLONAKIS G， GAZETAS G ．Lateral vibration and internal forces of grouped piles in layered soil［J］．Journal of Geotechnical and Geoenvironmental Engineering，1999，125（1）：16-25．

[12] 覃玉兰，邹新军，曹雄 ．均质砂土中水平简谐荷载与扭矩联合受荷单桩内力、位移分析［J］．岩土力学，2020，41（1）：147-156．

[13] QIN Y L， ZOU X J， CAO X ．Internal forces and deformations of a single pile in uniform sand under combined action of horizontal harmonic load and torque［J］．Rock and Soil Mechanics， 2020， 41（1）：147-156．（in Chinese）

[14] 杨超炜， 赵明华， 陈耀浩， 等 ．高陡横坡段桩柱式桥梁双桩基础受力分析［J］．湖南大学学报（自然科学版）， 2018， 45（3）：129-135．

[15] YANG C W， ZHAO M H， CHEN Y H，et al ．Force analysis of bridge double-piles foundation in high and steep cross slopes［J］．Journal of Hunan University （Natural Sciences）， 2018， 45（3）：129-135．（in Chinese）

[16] 赵明华， 陈耀浩， 杨超炜， 等 ．基于有限杆单元法的陡坡基桩非线性分析［J］．岩土力学， 2018， 39（8）： 3020-3028．

[17] ZHAO M H， CHEN Y H， YANG C W， et al ．Nonlinear analysis of pile in steep slope based on finite bar element method［J］．Rock and Soil Mechanics， 2018， 39（8）： 3020-3028．（in Chinese）

[18] PENG W Z， ZHAO M H， ZHAO H，et al ．A two-pile foundation model in sloping ground by finite beam element method［J］．Computers and Geotechnics，2020，122：103503．

[19] 赵海鹏， 黎学优， 万建宏， 等 ．基于高性能有限单元法的多层地基水平受荷桩受力分析［J］．岩土力学， 2021， 42（7）： 1995-2003．

[20] ZHAO H P， LI X Y， WAN J H，et al ．Analysis of laterally-loaded piles embedded in multi-layered soils using efficient finite-element method［J］．Rock and Soil Mechanics，2021，42（7）：1995-2003．（in Chinese）

[21] LI X Y， ZHAO H P， LIU S W，et al ．Innovative Timoshenko soil-pile integrated element for lager diameter laterally-loaded piles considering soil-pile interactions［J］．Computers and Geotechnics，2022，152：105020．

[22] 蒯行成，沈蒲生 ．层状介质中群桩水平动力阻抗的简化计算方法［J］．振动工程学报，1998，11（3）：258-264．

[23] KUAI X C， SHEN P S ．Simplified method for calculating horizontal dynamic impedances of pile groups in layered media［J］. Journal of Vibration Engineering， 1998， 11（3）： 258-264．（in Chinese）

[24] PASTERNAK P L. Fundamentals of a new method of analyzing structures on an elastic foundation by means of two foundation moduli［J］. Proceedings of the Gosudarstvennoe Izdatelstro Liberaturi Po Stroitelstvui Arkhitekture， 1954： 4-7.

[25] BHADURI A， CHOUDHURY D. Steady-state response of flexible combined pile-raft foundation under dynamic loading［J］. Soil Dynamics and Earthquake Engineering， 2021， 145： 106664．

[26] 王珏， 周叮， 刘伟庆， 等 ．考虑成层土剪切效应的相邻桩基动力相互作用［J］．振动工程学报， 2013， 26（5）： 732-742．

[27] WANG J， ZHOU D， LIU W Q，et al ．Dynamic interaction between adjacent piles considering soil shear effect in layered foundation［J］．Journal of Vibration Engineering，2013， 26（5）： 732-742．（in Chinese）

[28] WANG J， LO S H， ZHOU D ．Effect of a forced harmonic vibration pile to its adjacent pile in layered elastic soil with double-shear model［J］．Soil Dynamics and Earthquake Engineering， 2014， 67： 54-65．

[29] QETIN D， SIMSEK M ．Free vibration of an axially functionally graded pile with pinned ends embedded in Winkler-Pasternak elastic medium［J］．Structural Engineering and Mechanics，2011，40（4）： 583-594．

[30] HALABE U B， JAIN S K ．Lateral free vibration of a single pile with or without an axial load［J］．Journal of Sound and Vibration，1996，195（3）：531-544．

[31] 熊辉， 杨丰 ．文克尔地基模型下液化土桩基水平振动响应分析［J］．岩土力学， 2020， 41（1）： 103-110．

[32] XIONG H， YANG F ．Horizontal vibration response analysis of pile foundation in liquefied soil under Winkler foundation model［J］．Rock and Soil Mechanics，2020，41（1）：103-110． （in Chinese）

[33] FENG Z H， COOK R D ．Beam elements on two-parameter elastic foundations［J］．Journal of Engineering Mechanics，1983， 109（6）：1390-1402．

[34] 张玲， 岳梢， 赵明华， 等 ．基于改进Pasternak地基模型的桩柱式桥墩受力变形分析［J］．岩土工程学报， 2022， 44（10）： 1817-1826．

[35] ZHANG L， YUE S， ZHAO M H，et al ．Behaviors of pile-column piers based on modified Pasternak foundation model［J］．Chinese Journal of Geotechnical Engineering，2022， 44（10）： 1817-1826．（in Chinese）

[36] GAZETAS G， DOBRY R ．Horizontal response of piles in layered soils［J］．Journal of Geotechnical Engineering， 1984， 110（1）：20-40．

[37] TANAHASHI H ．Formulas for an infinitely long Bernoulli-Euler beam on the Pasternak model［J］．Soils and Foundations，2004，44（5）： 109-118．

[38] 梁发云， 李彦初， 黄茂松 ．基于Pasternak双参数地基模型水平桩简化分析方法［J］．岩土工程学报， 2013， 35（）：300-304．

[39] LIANG F Y， LI Y C， HUANG M S ． Simplified method for laterally loaded piles based on Pasternak double-parameter spring model for foundations［J］．Chinese Journal of Geotechnical Engineering，2013，35（Sup.1）：300-304．（in Chinese）

[40] YAO W J， YIN W X， CHEN J，et al ．Numerical simulation of a super-long pile group under both vertical and lateral loads［J］．Advances Structural Engineering，2010，13（6）：1139-1151．

[41] 吕冰 ．水平振动荷载作用下支盘桩动力特性研究［D］．洛阳：河南科技大学，2015．

[42] LÜ B ．Study on dynamic characteristics of branch pile under horizontal vibration load［D］．Luoyang：Henan University of Science and Technology，2015．（in Chinese）

[43] 胡安峰 ．桩基水平振动理论与性状研究［D］．杭州： 浙江大学， 2002．

[44] HU A F ．Study on horizontal vibration theory and characteristics of pile foundation［D］．Hangzhou：Zhejiang University，2002． （in Chinese）

[45] 梁发云， 陈海兵， 宋 著．轴向荷载对横向受荷桩承载性状影响的参数分析［J］．土木建筑与环境工程， 2011， 33（6）： 25-30．

[46] LIANG F Y， CHEN H B， SONG Z ．Parameter analysis of effects of axial load on the lateral response of piles［J］．Journal of Civil， Architectural & Environmental Engineering， 2011， 33（6）： 25-30．（in Chinese）

[47] 赵明华， 李微哲， 杨明辉，等 ．成层地基中倾斜偏心荷载下基桩位移特性室内模型试验研究［J］．土木工程学报，2006， 39（12）：95-99．

[48] ZHAO M H， LI W Z， YANG M H，et al ．A model test study on displacement of piles under inclined and eccentric loads in layered soils［J］．China Civil Engineering Journal，2006，39（12）：95-99．（in Chinese）

基金资助

国家自然科学基金资助项目(52068004)

National Natural Science