考虑蠕变影响的堆载作用下被动排桩受力变形分析

Analysis on Deformation of a Passive Pile under Surcharge Load Conditions Considering Creep Effects

摘要

为分析堆载作用下被动桩受力变形的时间效应，首先引入分数阶Merchant模型描述土体的蠕变特性，然后根据对应性原理和Laplace变换推导Boussinesq黏弹性解，计算堆载引发水平附加应力，并通过Terzaghi土拱模型将附加应力传递到桩上，得到堆载作用在桩上的被动荷载.其次，将桩简化为黏弹性Pasternak地基上的欧拉梁，建立桩身挠曲微分方程，并采用有限差分法及Laplace逆变换求解，随后通过与已有试验结果的对比分析，验证了本文方法的可行性.最后，对分数阶Merchant模型参数（虎克体的弹性模量E_h、Kelvin体的弹性模量E_k及黏滞系数η、分数阶阶次α）、堆载-桩身水平距离以及堆载荷载对桩身水平位移的影响进行了分析，分析结果表明：E_h越大，桩身初始水平位移越小；E_k越大，桩身水平位移随时间的增量越小；η越大，桩身到达最终变形量的时间越长；α越大，Kelvin体的黏滞性越大；堆载-桩身水平距离越小、堆载荷载越大，桩身水平位移越大；且随时间的增加，桩身水平位移对E_h的敏感度降低，对堆载-桩身水平距离和堆载荷载的敏感度会增大.

Abstract

To analyze the time effect of passive pile deformation under the action of surcharge load， this paper first introduces the fractional-order Merchant model to describe the creep characteristics of the soil. Then， based on the correspondence principle and Laplace transform， the Boussinesq viscoelastic solution is derived to calculate the additional horizontal stress induced by the load. The Terzaghi arching model is used to transfer the additional stress to the pile， obtaining the passive load on the pile due to the action of the load. Next， the pile is simplified as an Euler beam on a viscoelastic Pasternak foundation， and the differential equation calculating the deflection of the pile is established. The finite difference method and Laplace inverse transform are used to solve the equation. Subsequently， by comparing the results with existing experimental data， the feasibility of the method proposed in this paper is validated. Finally， the effects of fractional-order Merchant model parameters （elastic modulus E_h of Hooke body， elastic modulus E_k of Kelvin body， viscosity coefficient η， and fractional-order order α）， surcharge load-pile horizontal distance， and surcharge load on pile horizontal displacement are analyzed. The analysis results show that the larger E_h is， the smaller the initial pile horizontal displacement is. The greater E_k is， the smaller the increment of pile horizontal displacement with time， and the greater η results in the longer time for the pile to reach the final deformation. The larger α is， the greater the viscosity of Kelvin body is. The smaller surcharge load-pile horizontal distance and larger load lead to larger horizontal displacement of piles. With increasing time， the sensitivity of pile horizontal displacement to E_h decreases， and the sensitivity of surcharge load-pile horizontal distance and surcharge load increases.

Graphical abstract

关键词

堆载 / 桩基础 / 分数阶Merchant模型 / 黏弹性地基 / 有限差分法

当修建的公路工程临近房屋建筑、桥梁等建（构）筑物时，公路路堤将对建筑物桩基产生堆载效应，即：路堤堆载引起土体侧移，对临近桩基产生水平附加应力，以致桩基被动受荷，甚至可能偏位，严重时引发工程事故^［1-3］.

目前，堆载作用下被动桩受力变形研究已取得了不少成果^［4-6］.当被动桩位于软土地区时，需考虑土体蠕变特性对桩基长期变形的影响，国内外学者对此开展了相关研究.陆建飞等^［7］基于Merchant模型，对层状两相介质固结和蠕变的水平受荷的单桩问题进行求解.邓会元^［8］假定临界土压力随着时间增加由主动土压力趋于静止土压力，提出考虑时间效应的拱荷载计算方法.Zhang等^［9］引入Boltzmann模型，得到基坑开挖引起的桩基水平变形时域解.江杰等^［10］基于三参量模型推导出Mindlin的时域解，提出基坑开挖引起的土体附加应力的计算方法.闵鹏等^［11］引入Burgers模型，对堆载引发邻近单桩水平变形进行了分析.

然而，上述软土蠕变特性的考虑大多基于整数阶微分，由于整数阶微分算子的局限，约束了土体应力或应变随时间变化的路径^［12］，针对这一问题，有学者引入分数阶导数，通过采用分数阶模型对土体蠕变曲线进行拟合，试验结果分析表明分数阶模型能以少量的参数准确描述软土蠕变曲线^［12-13］.

另外，对于被动受荷排桩，因桩、土刚度差异，堆载产生的水平向附加应力引起桩和桩间土的水平向相对位移，进而在桩土之间产生水平向剪应力迁移现象，即通过土拱效应将水平附加应力传至桩身.因此，堆载作用下被动桩尚需考虑土拱效应影响.对此，Ito等^［14］基于塑性变形理论将桩前土拱简化为梯形，提出了桩身被动荷载理论计算公式；竺明星等^［15］、张浩等^［16］将其改进并应用于堆载下既有桩基；李忠诚等^［17］将被动桩分为被动侧成拱和主动侧形成应变楔，分析了侧移土体成拱效应；彭文哲等^［18］引入应变楔模型，提出适用于斜坡地基桩前土抗力计算的桩前土楔模型.但由于土体自重应力的存在，随着深度的增加，桩侧土拱高度逐渐减小^［19］，而上述研究均未全面考虑深度对拱高以及时间效应对桩身拱荷载的影响.

因此，为解决临近堆载作用下既有桩基长期变形问题，本文引入分数阶Merchant模型以及Terzaghi土拱模型分别描述土的蠕变特性和桩前土拱效应，并结合有限差分法，考虑土体的成层性，对考虑蠕变效应的堆载作用下被动排桩受力变形进行研究.

1 堆载下桩-土相互作用模型

图1为路堤堆载下邻近桩基受力变形计算模型，桩长L、桩径D、中心距S₁（桩间净距S₂）的桩基置于成层土中（桩顶至桩端土层数为n），桩中心到堆载底部边缘的距离为a，堆载沿着车辆行驶方向视为无限长，尺寸如图1（b）所示.

将桩体视作竖直放置的欧拉-伯努利梁，桩后土体视作黏弹性Pasternak地基，堆载在桩上作用被动荷载p_p（z，t）并使其发生水平位移，挤压桩后土体，桩后土体随之产生水平地基反力q（z，t）.对如图2所示的桩身微段dz进行分析，根据微段水平力和弯矩的平衡，可得如下方程组：

\left\{\begin{array}{l}{p}_{\mathrm{p}}\left(z,t\right)D\mathrm{d}z+Q\left(z,t\right)=\\ q\left(z,t\right)D\mathrm{d}z+Q\left(z,t\right)+\mathrm{d}Q\left(z,t\right),\\ p_{\mathrm{p}}\left(z,t\right)D\frac{\mathrm{d}{z}^{\mathrm{2}}}{\mathrm{2}}+M\left(z,t\right)+Q\left(z,t\right)\mathrm{d}z=\\ q\left(z,t\right)D\frac{\mathrm{d}{z}^{\mathrm{2}}}{\mathrm{2}}+M\left(z,t\right)+\mathrm{d}M\left(z,t\right)\end{array}\right.（1）

式中：M（z,t）、Q（z,t）分别为桩身弯矩和剪力.

由于蠕变特性的存在，土体在受力变形过程中，应力-应变关系会随时间不断变化，导致桩身被动荷载p_p（z，t）与桩后水平地基反力q（z，t）均随时间变化.为模拟p_p（z，t）以及q（z，t）随时间的变化，本文引入分数阶Merchant模型，将传统整数阶微分模型的Newton黏壶替换为分数阶元件Abel黏壶（图3），其本构关系如下：

\epsilon ={\sigma }_{\mathrm{0}}\left[\frac{\mathrm{1}}{E_{\mathrm{h}}}+\frac{\mathrm{1}}{E_{\mathrm{k}}}\left(\mathrm{1}-E_{\alpha }\left(-\frac{t^{\alpha }}{\eta }\right)\right)\right]（2）

式中：E_h、E_k和η分别为虎克体的弹性模量、Kelvin体的弹性模量和黏滞系数；α为分数阶阶次；t为时间；E_α（）为Mittag-Leffler函数，表达式如下：

E_{\alpha }\left(t\right)={\displaystyle \sum _{n=\mathrm{0}}^{\mathrm{\infty }}}\frac{\mathrm{1}}{\Gamma \left(\alpha n+\mathrm{1}\right)}{\left(-\frac{t^{\alpha }}{\eta }\right)}^n（3）

此处分数阶导数采用目前应用较广的Caputo型分数阶导数^［20］：

D^{\alpha }\left[f\left(t\right)\right]=\frac{\mathrm{d}{f}^{\alpha }\left(t\right)}{\mathrm{d}{t}^{\alpha }}=

\frac{\mathrm{1}}{\Gamma \left(n-\alpha \right)}{\int }_{\alpha }^t\frac{\mathrm{1}}{{\left(t-\tau \right)}^{\alpha -n+\mathrm{1}}}\frac{{\mathrm{d}}^{n}f\left(\tau \right)}{\mathrm{d}{t}^n}\mathrm{d}\tau ,

\left(n-\mathrm{1}\le \alpha \le n,n\in {\mathbf{Z}}^{+}\right)（4）

式中：Г（x）为Gamma函数^［21］，表达式如下：

\Gamma \left(x\right)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{\mathrm{e}}^{-t}{t}^{x-\mathrm{1}}\mathrm{d}t（5）

由式（2）可得分数阶Merchant模型蠕变过程中任意时刻的应变与应力比值，即蠕变柔量J（t）为：

J\left(t\right)=\frac{\mathrm{1}}{E_{\mathrm{h}}}+\frac{\mathrm{1}}{E_{\mathrm{k}}}\left[\mathrm{1}-E_{\alpha }\left(-\frac{t^{\alpha }}{\eta }\right)\right]（6）

2 桩身被动荷载时域解

2.1 基于分数阶Merchant模型的Boussinesq时域解

对于堆载作用下水平附加应力的确定，目前最常用方法为Boussinesq弹性理论解，在弹性半无限体表面上作用一竖向集中荷载F，半无限体内任一点处引起的水平附加应力为：

\begin{array}{l}{\sigma }_{\mathrm{0}}=\frac{\mathrm{3}F}{\mathrm{2}\mathrm{\pi }}\left[\frac{x^{\mathrm{2}}z}{R^{\mathrm{5}}}+\frac{\mathrm{1}-\mathrm{2}{\nu }_{\mathrm{s}}}{\mathrm{3}}\left(\frac{R^{\mathrm{2}}-R^{\mathrm{2}}z-z^{\mathrm{2}}}{R^{\mathrm{3}}\left(R+z\right)}-\right.\right.\\ \left.\left.\frac{x^{\mathrm{2}}\left(\mathrm{2}R+z\right)}{R^{\mathrm{3}}{\left(R+z\right)}^{\mathrm{2}}}\right)\right]\end{array}（7）

式中：{\nu }_{\mathrm{s}}为地基土的泊松比；F为作用在坐标原点处的竖向集中荷载；R为原点至点M的距离，R= （x²+y²+z²）^1/2.

当考虑土体蠕变特性时，认为土体的弹性模量E_s与泊松比{\nu }_{\mathrm{s}}会随时间变化，采用张治国等^［22］基于弹性-黏弹性对应原理计算Laplace域内三维Merchant模型的弹性模量E（s）与泊松比\nu（s）的方法，可得三维分数阶Merchant模型的E（s）和\nu（s）为：

\left\{\begin{array}{l}{E}_{\mathrm{s}}(s)=\frac{\mathrm{9}K\left(E_{\mathrm{h}}{E}_{\mathrm{k}}+E_{\mathrm{h}}\lambda s^{\alpha }\right)}{\mathrm{3}K\left(E_{\mathrm{h}}+E_{\mathrm{k}}+\lambda s^{\alpha }\right)+E_{\mathrm{h}}{E}_{\mathrm{k}}+E_{\mathrm{h}}\lambda s^{\alpha }},\\ {\nu }_{\mathrm{s}}(s)=\frac{\mathrm{3}K\left(E_{\mathrm{h}}+E_{\mathrm{k}}+\lambda s^{\alpha }\right)-\mathrm{2}{E}_{\mathrm{h}}{E}_{\mathrm{k}}-\mathrm{2}{E}_{\mathrm{h}}\lambda s^{\alpha }}{\mathrm{6}K\left(E_{\mathrm{h}}+E_{\mathrm{k}}+\lambda s^{\alpha }\right)+\mathrm{2}{E}_{\mathrm{h}}{E}_{\mathrm{k}}+\mathrm{2}{E}_{\mathrm{h}}\lambda s^{\alpha }}\end{array}\right.（8）

式中：λ=E_k^1-αη_k^α；s为Laplace变换算子.

对于黏弹性地基，需引入赫维赛德函数单位阶跃函数H（t），并设竖向荷载F（t）=FH（t），对F（t）进行Laplace变换可得：

F\left(s\right)=\frac{F}{s}（9）

将式（8）、式（9）代入式（7），可得竖向集中荷载F作用下半无限体内任一点M（x，y，z）的水平附加应力在Laplace域的表达式：

\begin{array}{l}{\sigma }_{\mathrm{0}}\left(s\right)=\frac{\mathrm{3}F}{\mathrm{2}\mathrm{\pi }}\left[\frac{x^{\mathrm{2}}z}{R^{\mathrm{5}}s}+T_{\mathrm{1}}\left(\frac{R^{\mathrm{2}}-R^{\mathrm{2}}z-z^{\mathrm{2}}}{R^{\mathrm{3}}\left(R+z\right)}-\right.\right.\\ \left.\left.\frac{x^{\mathrm{2}}\left(\mathrm{2}R+z\right)}{R^{\mathrm{3}}{\left(R+z\right)}^{\mathrm{2}}}\right)\right]\end{array}（10）

其中，

T_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{3}s}-\frac{\mathrm{3}K\left(E_{\mathrm{h}}+E_{\mathrm{k}}+\lambda s^{\alpha }\right)-\mathrm{2}{E}_{\mathrm{h}}{E}_{\mathrm{k}}-\mathrm{2}{E}_{\mathrm{h}}\lambda s^{\alpha }}{\mathrm{9}Ks\left(E_{\mathrm{h}}+E_{\mathrm{k}}+\lambda s^{\alpha }\right)+\mathrm{3}{E}_{\mathrm{h}}{E}_{\mathrm{k}}s+\mathrm{3}{E}_{\mathrm{h}}\lambda s^{\alpha +\mathrm{1}}}（11）

对式（10）进行Laplace逆变换，可得水平附加应力时域解：

\begin{array}{l}{\sigma }_{\mathrm{0}}\left(t\right)=\frac{\mathrm{3}F}{\mathrm{2}\mathrm{\pi }}\left[\frac{x^{\mathrm{2}}z}{R^{\mathrm{5}}}+M_{\mathrm{1}}\left(\frac{R^{\mathrm{2}}-Rz-z^{\mathrm{2}}}{R^{\mathrm{3}}\left(R+z\right)}-\right.\right.\\ \left.\left.\frac{x^{\mathrm{2}}\left(\mathrm{2}R+z\right)}{R^{\mathrm{3}}{\left(R+z\right)}^{\mathrm{2}}}\right)\right]\end{array}（12）

式中：M₁为T₁的Laplace逆变换.

建立如图4所示的坐标系以及堆载-桩基计算模型，桩基位于堆载左侧，桩长为L，堆载宽度为b₁+b₂+b₃、长度为L_d，b₂×L_d部分作用的均布荷载为Q.将堆载分成三个区域，分别为b₁×L_d区域、b₂×L_d区域与b₃×L_d区域，由于堆载作用下桩前水平附加应力应为Boussinesq积分解的2倍^［16］，根据应力叠加原理，得梯形堆载产生的水平附加应力如式（13）所示，并将相同高度处的水平附加应力视为均布荷载.

\begin{array}{l}{\sigma }_{\mathrm{0}}\left(t\right)=\frac{\mathrm{6}F}{\mathrm{\pi }}\left[-\frac{x}{b_{\mathrm{1}}}{\int }_{\mathrm{0}}^{L_{\mathrm{d}}}{\int }_{-b_{\mathrm{1}}}^{\mathrm{0}}{M}_{\mathrm{2}}\mathrm{d}x\mathrm{d}y+{\int }_{\mathrm{0}}^{L_d}{\int }_{-b_{\mathrm{1}}-b_{\mathrm{2}}}^{-b_{\mathrm{1}}}{M}_{\mathrm{2}}\mathrm{d}x\mathrm{d}y\right.\\ \left.-\frac{\left(x+b_{\mathrm{1}}+b_{\mathrm{2}}+b_{\mathrm{3}}\right)}{b_{\mathrm{3}}}{\int }_{\mathrm{0}}^{L_{\mathrm{d}}}{\int }_{-b_{\mathrm{1}}-b_{\mathrm{2}}-b_{\mathrm{3}}}^{-b_{\mathrm{1}}-b_{\mathrm{2}}}{M}_{\mathrm{2}}\mathrm{d}x\mathrm{d}y\right]\end{array}（13）

其中，

M_{\mathrm{2}}=\frac{x^{\mathrm{2}}z}{R^{\mathrm{5}}}+M_{\mathrm{1}}\left[\frac{R^{\mathrm{2}}-Rz-z^{\mathrm{2}}}{R^{\mathrm{3}}\left(R+z\right)}-\frac{x^{\mathrm{2}}\left(\mathrm{2}R+z\right)}{R^{\mathrm{3}}{\left(R+z\right)}^{\mathrm{2}}}\right]（14）

2.2 考虑土拱效应的桩身被动荷载计算

当排桩周围存在大面积堆载时，堆载将产生水平附加应力并引起土体侧移，而土体的侧移又受到桩体的遮拦，在桩前产生土拱效应，通过水平向剪应力迁移，将土体所受到的荷载转移到桩上.本文采用Terzaghi土拱模型分析桩前土拱效应并计算桩身上的拱荷载p_p（图5）.作如下假定：①剪切面竖直；②剪切阻力τ_x 沿剪切平面呈线性分布，τ_x 可由式（15）计算：

\left\{\begin{array}{l}{\tau }_x={\tau }_{\mathrm{1}}+\frac{x}{h}\left({\tau }_{\mathrm{2}}-{\tau }_{\mathrm{1}}\right),\\ {\tau }_{\mathrm{1}}={\beta }_{\mathrm{1}}\left(K_{\mathrm{e}}{p}_{\mathrm{0}}\mathrm{t}\mathrm{a}\mathrm{n}\phi +c\right),\\ {\tau }_{\mathrm{2}}={\beta }_{\mathrm{2}}\left(K_{\mathrm{e}}{p}_{\mathrm{p}}\mathrm{t}\mathrm{a}\mathrm{n}\phi +c\right)\end{array}\right.（15）

式中：x为距土拱拱顶的距离；c、φ分别为土的黏聚力和内摩擦角；K_e=tan²（45°-φ/2）；h为拱高，即等侧移面高度，表达式为：h=ξS₂-ψz，其中ξ为拱高系数，根据不同情况进行取值，参考值见表1，ψ为拱高随深度衰减系数，主要受土体的自重应力影响.李登峰等^［23］对花岗岩残积土进行研究，提出拱高随桩深度线性减小的斜率可取0.24或0.29.当h降至0时不考虑土拱效应，即p_p=σ₀，当a≥h时，桩前土体形成完整土拱，当a<h时，土拱为不完全土拱，此时取h=a；β₁、β₂分别为柱顶和柱底侧向摩擦力发挥程度系数.

从桩柱中取微段dx进行分析，根据水平力的平衡可得：

\begin{array}{l}{p}_x\left(x,z,t\right)D+\mathrm{2}{\tau }_x\mathrm{d}x=\\ \left[\mathrm{d}{p}_x\left(x,z,t\right)+p_x(x,z,t)\right]D\end{array}（16）

将式（15）代入式（16）并结合以下边界条件：

\left\{\begin{array}{l}{p}_x(\mathrm{0},z,t)={\sigma }_{\mathrm{0}}(z,t),\\ p_x(h,z,t)=p_{\mathrm{p}}(z,t)\end{array}\right.（17）

可得作用在桩身上的拱荷载p_p：

p_{\mathrm{p}}\left(z,t\right)=\frac{D{\sigma }_{\mathrm{0}}\left(t\right)+h{\beta }_{\mathrm{1}}\left[c+{\sigma }_{\mathrm{0}}\left(t\right)\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\left(\mathrm{45}°-\frac{\phi }{\mathrm{2}}\right)\right]+h{\beta }_{\mathrm{2}}c}{D-{\beta }_{\mathrm{2}}h\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\left(\mathrm{45}°-\frac{\phi }{\mathrm{2}}\right)}（18）

3 桩身水平变形时域解

3.1 桩后地基反力计算

桩后土体因被动桩受荷挠曲变形产生的水平地基反力采用如图1所示的黏弹性Pasternak模型计算，即对Pasternak地基中弹性层进行改进，在需要考虑蠕变的土层，将离散弹簧替换为分数阶Merchant模型（图3），剪切层继续维持原模型.

地基反力q（z，t）与对应的位移w（z，t）的关系为：

q\left(z,t\right)=k_{\mathrm{s}}\left(z,t\right)w\left(z,t\right)-G_{\mathrm{s}}\left(z,t\right)\frac{{\mathrm{d}}^{\mathrm{2}}w\left(z,t\right)}{\mathrm{d}{z}^{\mathrm{2}}}（19）

式中：k_s（z，t）为地基的反力模量，kN/m²；G_s（z，t）为地基土剪切刚度，kN/m.

本文考虑地基土埋深影响和土体蠕变特性对 k_s（z，t）的影响^{［11，30］}.

k_{\mathrm{s}}\left(z,t\right)=\left\{\begin{array}{l}\frac{\mathrm{0.65}{E}_{\mathrm{s}}}{D\left(\mathrm{1}-{{\nu }_{\mathrm{s}}}^{\mathrm{2}}\right)}\sqrt[\mathrm{12}]{\frac{E_{\mathrm{s}}{D}^{\mathrm{4}}}{EI}}{z}^{\mathrm{0.5}},\mathrm{不考}\mathrm{虑蠕}\mathrm{变区}\mathrm{域}\\ \frac{\mathrm{1}}{J\left(t\right)}{z}^{\mathrm{0.5}},\mathrm{考虑}\mathrm{蠕变}\mathrm{区域}\end{array}\right.（20）

式中：E_s、{\nu }_{\mathrm{s}}分别为地基土的弹性模量和泊松比；z为地基土埋深.

G_s采取Tanahashi^［31］提出的公式并考虑土体蠕变特性：

G_{\mathrm{s}}\left(z,t\right)=\frac{E_{\mathrm{s}}\left(t\right)h_{\mathrm{t}}}{\mathrm{6}\left[\mathrm{1}+{\nu }_s\left(t\right)\right]}（21）

式中：E_s（t）、{\nu }_{\mathrm{s}}（t）分别为时域内分数阶Merchant模型的弹性模量和泊松比；h_t为地基土的剪切层厚度，h_t=11D^［5］.

3.2 桩身挠曲微分方程的建立

将桩视为欧拉-伯努利梁，桩身弯矩和剪力为：

Q\left(z,t\right)=\frac{\mathrm{d}M\left(z,t\right)}{\mathrm{d}z}=-EI\frac{{\mathrm{d}}^{\mathrm{3}}w}{\mathrm{d}{z}^{\mathrm{3}}}（22）

式中：EI代表梁的弯曲刚度.

将式（19）、式（22）代入式（1）可得：

EI\frac{{\mathrm{d}}^{\mathrm{4}}w\left(z,t\right)}{\mathrm{d}{z}^{\mathrm{4}}}=p_{\mathrm{p}}\left(z,t\right)D-\left[k_{\mathrm{s}}\left(z,t\right)w\left(z,t\right)-G_{\mathrm{s}}\left(z,t\right)\frac{{\mathrm{d}}^{\mathrm{2}}w\left(z,t\right)}{\mathrm{d}{z}^{\mathrm{2}}}\right]D（23）

对式（23）化简可得：

\frac{{\mathrm{d}}^{\mathrm{4}}w\left(z,t\right)}{\mathrm{d}{z}^{\mathrm{4}}}-\frac{G_{\mathrm{s}}\left(z,t\right)D}{EI}\frac{{\mathrm{d}}^{\mathrm{2}}w\left(z,t\right)}{\mathrm{d}{\mathrm{z}}^{\mathrm{2}}}+\frac{k_{\mathrm{s}}\left(z,t\right)D}{EI}w\left(z,t\right)=\frac{D}{EI}{p}_{\mathrm{p}}\left(z,t\right)（24）

式（24）是一个含时间变量t的非齐次四阶常微分方程，直接求解较困难.为求解w（z，t），将上式中的t进行Laplace变换：

\frac{{\mathrm{d}}^{\mathrm{4}}w\left(z,s\right)}{\mathrm{d}{z}^{\mathrm{4}}}-\frac{G_{\mathrm{s}}\left(z,s\right)D}{EI}\frac{{\mathrm{d}}^{\mathrm{2}}w\left(z,s\right)}{\mathrm{d}{z}^{\mathrm{2}}}+\frac{k_{\mathrm{s}}\left(z,s\right)D}{EI}w\left(z,s\right)=\frac{D}{EI}{p}_{\mathrm{p}}\left(z,s\right)（25）

式中：w（z，s）、p_p（z，s）分别为w（z，t）、p_p（z，t）的Laplace变换；k_s（z，s）、G_s（z，s）分别为k_s（z，t）、G_s（z，t）的Laplace变换.

3.3 桩身挠曲微分方程的求解

式（25）为高阶齐次微分方程，难以直接求解，且桩身置于成层地基中，土体参数沿深度变化.因此，采用有限差分法对此模型进行数值求解.该方法的基本过程如下：1）对被研究对象进行合理的离散化；2）采用有限差分公式重写微分方程中的导数项，得到一系列方程；3）根据边界条件求解代数方程.若基桩两端自由，在两端节点处，剪力和弯矩均为零.

图6给出了桩身离散化模型.将桩分为n个长度为β的单元，并在桩的两端分别附加2个虚单元，每个单元对应不同的地基参数k_s（s）、G_s（s）.

边界条件如下：

Q_{\mathrm{0}}=Q_n=\mathrm{0}（26）

M_{\mathrm{0}}=M_n=\mathrm{0}（27）

根据有限差分公式，式（25）可化成：

\frac{w_{i+\mathrm{2}}-\mathrm{4}{w}_{i+\mathrm{1}}+\mathrm{6}{w}_i-\mathrm{4}{w}_{i-\mathrm{1}}+w_{i-\mathrm{2}}}{{\beta }^{\mathrm{4}}}-

\frac{G_{\mathrm{s}}\left(z,s\right)D}{EI}\frac{w_{i+\mathrm{1}}-\mathrm{2}{w}_i+w_{i-\mathrm{1}}}{{\beta }^{\mathrm{2}}}+

\frac{k_{\mathrm{s}}\left(z,s\right)D}{EI}{w}_i=\frac{D}{EI}{p_{\mathrm{p},}}_i（28）

式中：p_p，i 表示i节点处的桩侧附加应力，i=0，1，…，n-1，n.

将式（28）用矩阵形式表示如下：

\mathit{w}={\left({\mathit{L}}_{\mathrm{1}}-{\mathit{L}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\mathit{P}

其中，

\mathit{w}={\left[w_{\mathrm{0}},w_{\mathrm{1}},w_{\mathrm{2}},\cdots ,w_{n-\mathrm{1}},w_n\right]}_{\mathrm{1}\times \left(n+\mathrm{1}\right)}^{\mathrm{T}}

{\mathit{L}}_{\mathrm{1}}={\left[\begin{array}{ccccccc}\mathrm{2}{A}_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& \mathrm{2}{A}_{\mathrm{1}}& & & & \\ -\mathrm{2}{A}_{\mathrm{1}}& \mathrm{5}{A}_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& A_{\mathrm{1}}& & & \\ A_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& \mathrm{6}{A}_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& A_{\mathrm{1}}& & \\ & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & A_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& \mathrm{6}{A}_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& A_{\mathrm{1}}\\ & & & A_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& \mathrm{5}{A}_{\mathrm{1}}& -\mathrm{2}{A}_{\mathrm{1}}\\ & & & & \mathrm{2}{A}_{\mathrm{1}}& -\mathrm{4}{A}_{\mathrm{1}}& \mathrm{2}{A}_{\mathrm{1}}\end{array}\right]}_{\left(n+\mathrm{1}\right)\times \left(n+\mathrm{1}\right)}

{\mathit{L}}_{\mathrm{2}}={\left[\begin{array}{ccccc}-C_{\mathrm{0}}& & & & \\ B_{\mathrm{1}}& -\mathrm{2}{B}_{\mathrm{1}}-C_{\mathrm{1}}& B_{\mathrm{1}}& & \\ & \ddots & \ddots & \ddots & \\ & & B_{n-\mathrm{1}}& -\mathrm{2}{B}_{n-\mathrm{1}}-C_{n-\mathrm{1}}& B_{n-\mathrm{1}}\\ & & & & -C_n\end{array}\right]}_{\left(n+\mathrm{1}\right)\times \left(n+\mathrm{1}\right)}

\mathit{P}=I{\left[{p_{\mathrm{p},}}_{\mathrm{0}},{p_{\mathrm{p},}}_{\mathrm{1}},{p_{\mathrm{p},}}_{\mathrm{2}},\cdots ,{p_{\mathrm{p},}}_{n-\mathrm{2}},{p_{\mathrm{p},}}_{n-\mathrm{1}},{p_{\mathrm{p},}}_n\right]}_{\mathrm{1}\times \left(n+\mathrm{1}\right)}^{\mathrm{T}}

A_{\mathrm{1}}=\frac{\mathrm{1}}{{\beta }^{\mathrm{4}}}，B_i=\frac{G_{\mathrm{s},i}\left(z,s\right)D}{EI{\beta }^{\mathrm{2}}}，C_i=\frac{k_{\mathrm{s},i}\left(z,s\right)D}{EI}，I=\frac{D}{EI}

解得w（z，s）后，进行Lalpace逆变换可得时域内桩基水平位移w（z，t）：

w\left(z,t\right)=L^{-\mathrm{1}}\left[w\left(z,s\right)\right]（29）

其他边界条件也可按相同方法得到，限于篇幅，未罗列.

4 算例分析

在宁波货运北站所处的软土地区，由于堆载作用，修建在该处的部分铁路桩基发生水平偏移.李双龙等^［32］为分析该地区软土蠕变特性对堆载作用下桩基水平偏移的影响，开展了一维软土固结-蠕变试验与数值试验.本节基于李双龙等的研究，对蠕变模型以及桩土相互作用模型进行验证.

4.1 蠕变模型验证

为对该地区软土蠕变特性进行理论描述，经初步对比，选取具有理论依据且各个参数具有明确物理意义的经典元件模型和分数阶蠕变模型进行拟合，这两类模型均为不同基本元件按串、并联的方式组合而成，通常，组合的元件越多，对土体蠕变的模拟效果越好，但也会增加待定参数的数量，为使用较少的元件达到较好的预测效果，重点分析Merchant模型、Burgers模型和分数阶Merchant模型.图7为各种蠕变模型的预测结果.

为定量分析不同蠕变模型对该地区软土蠕变特性的表征能力，表2以相关系数R²作为指标对各模型的预测效果进行了对比.

由表2可知，Merchant模型待定参数数量为3个，虽然模型简单、待定参数较少，但不能较好地预测大荷载长期作用下的土体蠕变.Burgers模型和分数阶Merchant模型待定参数数量均为4个，但分数阶Merchant模型的预测效果整体优于Burgers模型.

4.2 桩土相互作用模型验证

为验证本文桩土相互作用模型的有效性，对李双龙等^［32］的数值试验进行分析，试验布置如图8所示，堆载宽度B=10 m，堆载长度M=25 m，堆载到桩中心的距离为a=10 m，堆载大小为Q=70 kPa，邻近三桩承台的尺寸为7.4 m×1.4 m×2 m，桩中心距3 m，桩入土深度45 m，桩径1 m，桩弹性模量33.6×10⁶ kPa，地基土层情况及其物理力学参数见图8，试验其他信息详见文献［27］.

竺明星^［33］指出，长期荷载作用下边桩所受荷约为中间桩的67%，取3根桩的平均荷载作为桩上荷载，并考虑到自重应力会使拱高随着深度的增加逐渐变小，且重力越大，减小速率越大^{［19，34］}，由于文献［32］的软土重度在1.7 g/cm³左右，略小于李登峰等^［23］研究的花岗岩残积土重度1.9 g/cm³，故取0.24作为本文拱高随深度的减小斜率.根据文献［32］的数据进行反演得分数阶模型的参数见表3，采用本文方法计算的桩身水平位移与李双龙等的数值试验结果的对比见图9.

由图9可知，桩身变形本文解与数值解的基本规律一致，但本文按理论方法计算时未能考虑承台对桩身侧移的约束作用，在前期弹性变形以及桩间土初始蠕变阶段，桩身水平位移较小，承台约束效应影响并不明显，理论计算结果与考虑承台作用的数值分析较为接近.在后期蠕变稳定阶段，桩身水平位移增大，承台侧向约束作用增大，以致本文计算结果较数值分析结果偏大.由此可见，本文方法可较好地模拟单侧堆载作用下桩基的受力变形.

5 分数阶Merchant模型参数分析

本节在4.2节的基础上，对比分析分数阶Merchant模型中各个元件参数对被动桩长期受力变形的影响，计算分析时基本工况如下：堆载宽度B=10 m，堆载长度M=25 m，堆载到桩的距离a=10 m，堆载大小Q=70 kPa，邻近三桩的桩中心距为3 m，桩入土深度为45 m，桩径为1 m，桩弹性模量为33.6×10⁶ kPa.

5.1 虎克体弹性模量的影响分析

在基本工况其他参数不变的情况下，取土层参数E_k=600 kPa，η_k=16 kPa·d，α=0.4，K=10 000 kPa，E_h分别为1 000 kPa、2 000 kPa、4 000 kPa和8 000 kPa探究其对堆载下桩身水平位移的影响，t=0和t=100时E_h对桩身水平位移的影响如图10所示.由图10可知，E_h对桩身的瞬时变形和延时变形均有影响，E_h越小，桩身水平位移越大，桩身水平位移对E_h越敏感，但随着时间的增加，桩身水平位移对E_h的敏感度会降低.由此可知，E_h主要控制土的初始变形能力.

5.2 Kelvin体弹性模量的影响分析

在基本工况其他参数不变的情况下，取土层参数E_h=2 000 kPa，η_k=16 kPa·d，α=0.4，K=10 000 kPa，E_k分别为300 kPa、600 kPa、1 200 kPa和2 400 kPa探究其对堆载下桩身水平位移的影响，t=0和t=100时E_k对桩身水平位移的影响如图11所示.由图11可知，在堆载初期不同E_k下的桩身水平位移基本一致，随着时间的增加，桩身水平位移明显增大，且E_k越大，桩身水平位移随时间的增量越小.由此可知，E_k主要控制土的延时变形能力.

5.3 Kelvin体黏滞系数的影响分析

在基本工况其他参数不变的情况下，取土层参数E_h=2 000 kPa，E_k=600 kPa，α=0.4，K=10 000 kPa，η_k分别为8 kPa·d、16 kPa·d、32 kPa·d和64 kPa·d，得到桩顶位移随时间的变化如图12（a）所示；取土层参数E_h=8 600 kPa，E_k=6 200 kPa，α=0.6，K=6 800 kPa，η_k分别为8 000 kPa·d、16 000 kPa·d、32 000 kP·da和 64 000 kPa·d^［12］，得到桩顶位移随时间的变化如图 12（b）所示.由图12可知，η_k对桩身的瞬时变形和最终变形量几乎没有影响，但对到达最终变形量的时间有显著的影响，η_k越大，到达最终变形量的时间越长.由此可知，η_k主要控制土体的延时变形速率.

5.4 分数阶阶次的影响分析

在基本工况其他参数不变的情况下，取α分别为0.2、0.4、0.6、0.8和1.0，土层其他参数E_h=2 000 kPa，E_k=600 kPa，η_k=16 kPa·d，K=10 000 kPa，得到桩顶水平位移随时间的变化见图13（a）；取土层其他参数 E_h=8 600 kPa，E_k=6 200 kPa，η_k=160 000 kPa·d，K= 6 800 kPa，得到桩顶水平位移随时间的变化见图13（b）.由图13（a）可知，当η_k较小时，α对土体延时变形速率的影响显著，随着α的增大，土体延时变形速率迅速增大.由图13（b）可知，当η_k较大时，桩顶水平位移随时间的变化可分为两个部分：在第一部分，桩顶水平位移随α的增大而减小；在第二部分，桩顶水平位移随α的增大而增大，这是由于当α=0时，Abel黏壶退化为弹性体，当α=1时，Abel黏壶退化为牛顿黏壶，因此随着α的增大，Abel黏壶的黏滞性逐渐显著.

6 堆载参数分析

本节在4.2节的基础上，对堆载到桩身距离以及堆载大小对被动桩长期受力变形的影响进行分析，计算分析时土层参数如下：E_h=2 000 kPa，E_k= 600 kPa，η_k=16 kPa·d，α=0.4，K=10 000 kPa，其余参数同第5节基本工况.

6.1 堆载到桩身距离的影响分析

在基本工况其他参数不变的情况下，取堆载到桩身距离a分别为5 m、10 m、15 m、20 m，t=0和t=100时a对桩身水平位移的影响如图14所示.由图14可知，桩身水平位移随堆载到桩身距离a的增加呈非线性减小，且随时间的增加，不同a之间水平位移的差距逐渐增大.这是由于当堆载区域距桩较近时，桩身上部被动荷载急剧增大，荷载随时间的增量也同时增大.

6.2 堆载荷载的影响分析

在基本工况其他参数不变的情况下，取堆载荷载Q分别为50 kPa、100 kPa、150 kPa、200 kPa，t=0和t=100时Q对桩身水平位移的影响如图15所示.由图15可知，桩身水平位移随堆载荷载Q的增加而稳定增加，且随时间的增加，不同Q之间水平位移的差距逐渐增大.由此可知，当Q增加时，桩身被动荷载稳定增加，且桩身被动荷载越大，其随时间的增加也越大.

7 结 论

本文引入分数阶Merchant模型，结合Terzaghi土拱模型以及有限差分法，获得了考虑蠕变影响的堆载引起的邻近桩基水平位移时域解，并对分数阶Merchant模型以及堆载的各个参数影响进行了分析，得到以下结论：

1）分数阶Merchant模型中，E_h主要控制土的初始变形能力，E_h越大，桩身初始水平位移越小；E_k主要控制土的延时变形能力，E_k越大，桩身水平位移随时间的增量越小；η_k主要控制土体的延时变形速率，η_k越大，桩身到达最终变形量的时间越长；分数阶阶次α通过控制Abel黏壶黏弹性，进而控制整个Kelvin体的黏弹性.

2）堆载与邻近桩基的距离和堆载荷载均对桩身的水平位移有较大影响，堆载与邻近桩基的距离越近，桩身水平位移增长量越大；随着堆载荷载的增大，桩身水平位移稳定增加.

3）随时间的增加，桩身水平位移对E_h的敏感度会降低，但对堆载与邻近桩基距离和堆载荷载的敏感度会增大.

参考文献

[1] 赵明华，肖燕，陈昌富，等 ．考虑土体蠕变特性的桥台软基变形分析［J］．中国公路学报，2006，19（2）：56-61．

[2] ZHAO M H， XIAO Y， CHEN C F，et al ．Analysis of deformation on soft subsoil around bridge abutment considering soil creep property［J］．China Journal of Highway and Transport，2006， 19（2）：56-61．（in Chinese）

[3] 刘敦平，蒯行成，赵明华 ．软土运动作用下被动桩桩-土水平相互作用的三维有限元分析［J］．中国公路学报，2008，21（4）：18-24．

[4] LIU D P， KUAI X C， ZHAO M H ．3-D finite element analysis on pile-soil horizontal interaction of passive pile under soft soil movement action［J］．China Journal of Highway and Transport，2008，21（4）：18-24．（in Chinese）

[5] 蒯行成，刘敦平，赵明华 ．堆载作用下考虑界面接触行为的桩-土相互作用分析［J］．湖南大学学报（自然科学版），2006，33（6）：21-24．

[6] KUAI X C， LIU D P， ZHAO M H ．Analysis of pile-soil interaction under a surcharge load considering the contact behavior of interface［J］．Journal of Hunan University （Natural Sciences），2006，33（6）：21-24．（in Chinese）

[7] ZHANG H， SHI M， YANG L，et al ．A semianalytical solution for passively loaded piles adjacent to surcharge load［J］．Advances in Civil Engineering，2020，2020：1-19．

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

[9] 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）

[10] 竺明星， 卢红前， 王磊， 等. 深厚软土地区堆载作用下邻近基桩水平承载力研究［J］. 武汉大学学报（工学版）， 2017， 50（）：374-383．

[11] ZHU M X， LU H Q， WANG L，et al ．Research on lateral bearing capacity of pile foundation adjacentto surcharge loading on deep soft ground［J］．Engineering Journal of Wuhan University，2017，50（Sup.1）：374-383．（in Chinese）

[12] 陆建飞，王建华，沈为平 ．考虑固结和流变的层状地基中的水平单桩的理论分析［J］．岩石力学与工程学报，2001，20（3）：386-390．

[13] LU J F， WANG J H， SHEN W P ．Theoretical study of horizontally loaded single pile in layered soil considering consolidation and rheology of soil［J］．Chinese Journal of Rock Mechanics and Engineering，2001，20（3）：386-390．（in Chinese）

[14] 邓会元 ．滨海吹填围垦区堆载作用下桩基承载特性研究［D］．南京：东南大学，2021．

[15] DENG H Y ．Study on bearing characteristics of pile foundation under surcharge load in coastal reclamation area［D］．Nanjing：Southeast University，2021．（in Chinese）

[16] ZHANG Z G， HUANG M S， ZHANG C P，et al ．Time-domain analyses for pile deformation induced by adjacent excavation considering influences of viscoelastic mechanism［J］. Tunnelling and Underground Space Technology， 2019， 85（3）： 392-405.

[17] 江杰，张探，欧孝夺，等 ．软土地基基坑开挖对坑底桩受力与位移影响的时效分析［J］．科学技术与工程，2021，21（25）：10880-10886．

[18] JIANG J， ZHANG T， OU X D，et al ．Time effect analysis of foundation pit excavation on pile axial force and displacement in soft soil foundation［J］．Science Technology and Engineering，2021，21（25）：10880-10886．（in Chinese）

[19] 闵鹏，申玉生，林作忠，等 ．考虑时变效应的堆载引发邻近单桩水平变形分析［J］．中南大学学报（自然科学版），2022，53（8）：3129-3143．

[20] MIN P， SHEN Y S， LIN Z Z，et al ．Analysis of horizontal deformation of existing adjacent single pile induced by surface surcharge with time-dependent mechanism［J］．Journal of Central South University （Science and Technology），2022，53（8）：3129-3143．（in Chinese）

[21] 江杰，张探，欧孝夺，等 ．软土地基基坑开挖对临近桩变形影响的时效分析［J］．湖南大学学报（自然科学版），2022，49（11）：206-215．

[22] JIANG J， ZHANG T， OU X D，et al ．Time-effect analysis on influence of excavation on adjacent pile deformation in soft soil［J］．Journal of Hunan University （Natural Sciences），2022， 49（11）：206-215．（in Chinese）

[23] 艾智勇，叶梓坤，刘文杰 ．基于分数阶黏弹性地基模型的群桩时效行为研究［J］．岩土工程学报，2022，44（4）：749-754．

[24] AI Z Y， YE Z K， LIU W J ．Time-behavior of pile groups based on fractional derivative soil model［J］．Chinese Journal of Geotechnical Engineering，2022，44（4）：749-754．（in Chinese）

[25] ITO T， MATSUI T ．Methods to estimate lateral force acting on stabilizing piles［J］．Soils and Foundations，1975，15（4）：43-59．

[26] 竺明星， 龚维明， 何小元， 等 ．堆载作用下考虑土拱效应的被动桩变形内力半解析解［J］．岩土工程学报， 2013，35（11）：1997-2008．

[27] ZHU M X， GONG W M， HE X Y，et al ．Semi-analytical solution to deformation and internal force of passive piles under surcharge loads considering soil arching effect［J］．Chinese Journal of Geotechnical Engineering，2013，35（11）：1997-2008．（in Chinese）

[28] 张浩，石名磊，郭院成，等 ．不平衡堆载作用下邻近结构桩的侧向受力机制［J］．岩土工程学报，2016，38（12）：2226-2236．

[29] ZHANG H， SHI M L， GUO Y C，et al ．Lateral mechanical behaviors of structural piles adjacent to imbalanced surcharge loads［J］．Chinese Journal of Geotechnical Engineering，2016， 38（12）：2226-2236．（in Chinese）

[30] 李忠诚， 梁志荣 ．侧移土体成拱效应及被动桩计算模式分析［J］．岩土工程学报，2011，33（）：113-118．

[31] LI Z C， LIANG Z R. Analysis of arching effect of lateral soil and calculation mode of passive pile［J］. Chinese Journal of Geotechnical Engineering，2011，33（Sup.1）：113-118．（in Chinese）

[32] 彭文哲， 赵明华， 杨超炜， 等 ．斜坡地基桩前土抗力的应变楔模型修正［J］．中南大学学报（自然科学版），2020，51（7）：1936-1945．

[33] PENG W Z， ZHAO M H， YANG C W，et al ．Modification of strain wedge model for soil resistance in front of piles in sloping ground［J］．Journal of Central South University （Science and Technology），2020，51（7）：1936-1945．（in Chinese）

[34] 韩同春， 邱子义， 豆红强 ．基于颗粒离散元的抗滑桩土拱效 应分析［J］．中南大学学报（自然科学版）， 2016， 47（8）： 2715-2722．

[35] HAN T C， QIU Z Y， DOU H Q ．Soil arching effect between anti-slide piles based on YADE discrete element method［J］．Journal of Central South University （Science and Technology），2016， 47（8）：2715-2722．（in Chinese）

[36] 陈文，孙洪广，李西成，等 ．力学与工程问题的分数阶导数建模［M］．北京：科学出版社，2010．

[37] CHEN W， SUN H G， LI X C， et al. Fractional derivative modeling of mechanics and engineering problems［M］. Beijing： Science Press，2010．（in Chinese）

[38] B И 斯米尔诺夫. 高等数学教程［M］. 北京： 人民教育出版社， 1979.

[39] И SMIRNOV B ．Advanced mathematics course［M］．Beijing：People’s Education Press，1979．（in Chinese）

[40] 张治国，鲁明浩，宫剑飞 ．黏弹性地基中基坑开挖对邻近桩基变形影响的时域解［J］．岩土力学，2017，38（10）：3017-3028．

[41] ZHANG Z G， LU M H， GONG J F ．Time-domain solution for influence of excavation on adjacent pile deformation in visco-elastic foundation［J］．Rock and Soil Mechanics，2017，38（10）：3017-3028．（in Chinese）

[42] 李登峰，胡卸文，赵晓彦，等 ．花岗岩残积土边坡水平拱高竖向变化规律［J］．西南交通大学学报，2016，51（5）：1024-1032．

[43] LI D F， HU X W， ZHAO X Y，et al ．Variation of horizontal arch height of granite residual soil slope in vertical direction［J］．Journal of Southwest Jiaotong University，2016，51（5）：1024-1032．（in Chinese）

[44] Code of practice for Strengthened/reinforced soils and other fills：BS 8006-1：2010 ［S］．London： British Standards Institution，2010.

[45] TERZAGHI K ．Theoretical soil mechanics［M］. New York：John Wiley & Sons， Inc.， 1943．

[46] CARLSON B O ．Armerad Jord beräkningsprinciper för banker på pålar［J］．Terranova，Distr.SGI，Linköping，1987．（in Swedish）

[47] JENCK O， DIAS D， KASTNER R ．Two-dimensional physical and numerical modeling of a pile-supported earth platform over soft soil［J］．Journal of Geotechnical and Geoenvironmental Engineering， 2007， 133（3）： 295-305．

[48] CHEN R P， CHEN Y M， HAN J，et al ．A theoretical solution for pile-supported embankments on soft soils under one-dimensional compression［J］．Canadian Geotechnical Journal，2008，45（5）：611-623．

[49] RUI R， VAN TOL F， XIA X L，et al ．Evolution of soil arching； 2D DEM simulations［J］．Computers and Geotechnics，2016，73：199-209．

[50] 张红 ．水平受荷单桩分析方法研究［D］．哈尔滨：哈尔滨工业大学，2012．

[51] ZHANG H ．Research on analysis method of single pile under lateral load［D］. Harbin：Harbin Institute of Technology，2012．（in Chinese）

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

[53] 李双龙，魏丽敏，冯胜洋，等 ．基于扩展Koppejan模型的被动桩-软土时效性相互作用研究［J］．岩土力学，2022，43（9）：2602-2614．

[54] LI S L， WEI L M， FENG S Y，et al ．Time-dependent interactions between passive piles and soft soils based on the extended Koppejan model［J］．Rock and Soil Mechanics，2022，43（9）：2602-2614．（in Chinese）

[55] 竺明星 ．组合荷载作用下被动桩承载机理研究［D］．南京：东南大学，2016．

[56] ZHU M X ．Research on bearing mechanism of passive pile under combined loads［D］．Nanjing：Southeast University，2016．（in Chinese）

[57] 吕韶全， 孙狂飙， 王少锋， 等. 考虑土体自重应力影响的抗滑桩三维土拱效应［J］. 西安建筑科技大学学报（自然科学版），2021， 53（4）： 510-517．

[58] LÜ S Q， SUN K B， WANG S F，et al ．Three-dimensional soil arching effect considering the influence of soil gravity stress on anti-slide pile［J］．Journal of Xi’an University of Architecture & Technology （Natural Science Edition），2021，53（4）：510-517．（in Chinese）

基金资助

国家自然科学基金资助项目（52078205,52108317）; 湖南省自然科学基金资助项目（2020JJ3013）