砂土中桩循环打入过程大变形数值模拟

黄中原; 杨仲轩; 郭宁

doi:10.11779/CJGE20211503

砂土中桩循环打入过程大变形数值模拟 English Version

黄中原^{1, 2,},
杨仲轩^{1, 2, ,},
郭宁^{1, 2}

1.
浙江大学岩土工程计算中心, 浙江杭州 310058
2.
浙江省城市地下空间开发工程技术研究中心, 浙江杭州 310058

基金项目:

国家自然科学国际合作重点项目 52020105003

国家杰出青年基金项目 51825803

详细信息

作者简介:
黄中原(1997—)，男，硕士研究生，主要从事打入桩数值模拟方面的研究工作。E-mail: 21912014@zju.edu.cn

通讯作者:
杨仲轩, E-mail: zxyang@zju.edu.cn

中图分类号: TU435
计量
- 文章访问数: 292
- HTML全文浏览量: 63
- PDF下载量: 132
出版历程
- 收稿日期: 2021-12-19
- 网络出版日期: 2023-02-23
- 刊出日期: 2023-01-31

Numerical simulation of large deformation of piles in sand during cyclic penetration

HUANG Zhongyuan^{1, 2,},
YANG Zhongxuan^{1, 2, ,},
GUO Ning^{1, 2}

1.
Computing Center for Geotechnical Engineering, College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China
2.
Engineering Research Center of Urban Underground Space Development of Zhejiang Province, Hangzhou 310058, China

摘要

摘要: 深入认识桩打入过程引起的桩周土体应力场变化对准确计算打入桩的竖向承载力十分重要。在有限元分析软件ABAQUS中采用大变形任意拉格朗日-欧拉（ALE）方法以及状态相关莫尔-库仑砂土本构模型，模拟了室内标准模型槽试验，得到了桩循环打入过程中桩周砂土全阶段的应力变化规律，并重点分析了径向应力，验证了桩侧砂土应力水平的h/R效应（h为测点与桩尖竖向距离，向上为正；R为桩半径），获得了应力与桩尖和桩中心线相对位置的分布规律。数值模拟结果与室内模型试验实测结果较一致，验证了数值模拟方法的可靠性；有效补充了试验中难以测量的近桩侧区域的应力分布数据，可为打入桩的承载力设计理论和设计方法研究提供依据。
- 砂土 /
- 打入桩 /
- 数值模拟 /
- 循环打入 /
- 应力场 /
- 大变形
Abstract: It is extremely important to understand the stress conditions around the pile shaft during pile installation for accurately predicting its bearing capacity. The large-deformation arbitrary Lagrangian-Eulerian (ALE) finite element method in ABAQUS with the state-dependent Mohr-Coulomb sand model is employed to simulate the highly instrumented calibration chamber tests. The stress regimes developed at different penetration stages are obtained with a focus placed on the radial stress. The h/R effects of the stress in the surrounding soil are validated, along with the stress distribution varying with the relative distance with respective to the pile tip and pile centerline. The numerical simulation results are in good agreement with the experimental ones, indicating the robustness of the numerical model. The stress distribution obtained from the numerical analysis fills the blank of the stress data near pile shaft in experiments and is considered to be useful for improving the design methods for driven piles.
- sand /
- driven pile /
- numerical simulation /
- cyclic penetration /
- stress field /
- large deformation

HTML全文

0. 引言

在波浪荷载作用下，海床内孔隙水压力增大以及垂直有效应力降低，进而导致海床土体失稳甚至液化^[1]。在人类探索海洋资源的近几十年内，海床的液化对人类在海洋中的作业安全产生了很大的影响^[2]。因此，针对波浪荷载对海床动力响应的研究具有重要的理论意义和工程价值。

自20世纪中叶起，学者们开始分别基于刚性非耦合模型^[3]、Biot流固耦合模型^[4]和多孔介质弹塑性模型^[5]对波浪荷载引起的海床动力响应问题进行了理论研究。其中基于可考虑饱和介质压缩性的Biot三维准静态固结理论^[6]和动力固结理论应用^[7]最为广泛。Yamamoto等^[8]基于Biot准静态固结理论得到了线性波作用下无限厚度海床内孔隙水压力、有效应力和土骨架位移的闭合形式解，为此后的研究奠定了基础。Madsen^[9]在此基础上，引入了Verruijt^[10]提出的考虑海床土体饱和度的储存方程，得到了线性波作用下无限厚度准饱和海床动力响应的解析解，发现海床的动力响应受到海床土体饱和度的影响显著。基于Madsen的研究，Okusa^[11]发现海床土体的参数对无限厚度海床动力响应的影响可以忽略不计。Mei等^[12]采用边界层近似法，即通过将海床分为上下两部分，求解无限厚度海床动力响应的方法求解上半区域海床，并得到下半区域海床的简化解，最终得到波浪荷载作用下有限厚度多孔弹性海床的动力响应的近似解。Hsu等^[13]则得到了波浪荷载作用下有限厚度海床动力响应的解析解，发现当海床厚度趋于无穷大时，有限厚度海床的动力响应与无限厚度海床的解析解结果基本一致。

Zienkiewicz等^[14]基于Biot准静态固结理论和Biot动力固结理论，整理总结得到了适应于饱和土体的准静态理论、部分动态理论和完整动态理论。Ulker等^[15]通过对准静态解、部分动态解和完整动态解的对比分析研究，发现在一定的海床和波浪条件下三者得到的结果基本一致。考虑到准静态理论的推导简单，便于用于实际工程分析，于是学者们在准静态理论的基础上进行了一系列深入的研究工作。海洋中波浪特性往往相当复杂^[16]，而且波浪非线性特性对海床动力响应的研究尚不够深入。Madsen^[9]虽然得到了无限厚度海床动力响应的解答，但是该研究没有考虑海床底部下卧基岩的情况以及波浪非线性的影响。王忠涛等^[17]则通过数值模拟的方法给出了波浪非线性对海床动力响应的影响。郭秀军等^[18]研究基于有限差分法研究了海床的分层对动力响应的影响。Zhou等^[19]和Zhang等^[20]考虑了海床分层的影响得到了二阶Stokes波作用下多层海床的半解析解。周晓智等^[21]研究了驻波作用下海床应力路径特性。但是前述学者的研究均没有考虑波浪参数对海床动力响应的影响，且没有得到二阶Stokes波作用下有限厚度海床动力响应的准静态解析解。

因此本文基于准静态理论求解均质各向同性海床在二阶Stokes波作用下动力响应的解析解，通过与既有解进行对比验证，并结合数值计算分析了波浪和海床土体的参数对海床动力响应的影响规律。

1. 控制方程和边界条件

如图 1所示，本文考虑了下卧刚性不透水基岩的多孔弹性海床在二阶Stokes波作用下动力响应问题。其中二阶Stokes波的传播方向为x轴的正方向。与Li等^[1]的研究一致，取海床和海洋的接触面为笛卡尔坐标系的零点，z轴垂直向上时为正值，负值在海床和海洋的接触面以下。此外，图中H和L分别为二阶Stokes波的波高和波长；d为水深；h为海床的厚度。

图 1 波浪和多孔弹性海床之间的动力相互作用

Figure 1. Dynamic interaction between ocean waves and poroelastic seabed

下载: 全尺寸图片幻灯片

1.1 控制方程

基于Biot三维准静态固结理论^[6]与Verruijt储能方程^[10]，多孔弹性海床的控制方程可以表示为

${\sigma _{ij,j}} = 0 \text{，}$

(1)

$- {p_{,i}} = {{{\rho _{\text{f}}}g{{\dot w}_i}} \mathord{\left/ {\vphantom {{{\rho _{\text{f}}}g{{\dot w}_i}} K}} \right. } K} 。$

(2)

式中：下标i和j表示该场量沿x或z方向的分量；下标‘, i’和‘, j’表示该场量对x或z方向的一阶偏导；字母上的‘•’表示该场量对时间t的一阶导数； ${\sigma _{ij}}$ 为总应力张量；p为孔隙水压力； ${\rho _{\text{f}}}$ 为孔隙流体密度；g为重力加速度；K为渗透系数（水力传导系数）。

流体的质量守恒方程为

${\dot u_{i,i}} + {\dot w_{i,i}} = - n\beta \dot p = \partial {\rho _{\text{f}}}/\partial t 。$

(3)

式中：u为土骨架的位移；w为孔隙流体相对土骨架的位移；n为海床土体的孔隙率；β为孔隙流体压缩性系数。当海床土体中含有少量气体时，β可以通过下式计算：

$\beta = 1/{K_{\text{w}}} + (1 - {S_{\text{r}}})/{P_{{\text{w0}}}} 。$

(4)

式中： ${K_{\text{w}}}$ 为孔隙流体的实际体积弹性模量，一般取 ${K_{\text{w}}}$ =2×10⁹ Pa； ${S_{\text{r}}}$ 为海床土体的饱和度； ${P_{{\text{w0}}}}$ 为由水深d确定的绝对压力， ${P_{{\text{w0}}}}$ = ${\rho _{\text{f}}}$ gd。

将式（2）代入式（3），消去w项，可以将控制方程改写为

$K\left( {\frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{{\partial ^2}p}}{{\partial {z^2}}}} \right) - {\gamma _{\text{w}}}n\beta \frac{{\partial p}}{{\partial t}} = {\gamma _{\text{w}}}\frac{\partial }{{\partial t}}\left( {\frac{{\partial {u_x}}}{{\partial x}} + \frac{{\partial {u_z}}}{{\partial z}}} \right) 。$

(5)

式中： ${\gamma _{\text{w}}}$ 为孔隙流体的重度。

根据式（1），采用 ${\sigma^{\prime}_{xx}}$ ， ${\sigma{zz}^{\prime}_}$ 和τ_xz分别表示沿x方向和z方向的有效应力以及剪切应力，那么多孔弹性海床的整体平衡方程可以表示为

$\frac{{\partial {{\sigma^{\prime}}_{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{xz}}}}{{\partial z}} = - \frac{{\partial p}}{{\partial x}} \text{，}$

(6)

$\frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {{\sigma^{\prime}}_{zz}}}}{{\partial z}} = - \frac{{\partial p}}{{\partial z}} 。$

(7)

多孔弹性海床的本构关系可以表示为

$\left.\begin{array}{l}{{\sigma }^{\prime }}_{xx}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\partial {u}_{x}}{\partial x}+\frac{\mu }{1-\mu }\frac{\partial {u}_{z}}{\partial z}\right)\text{ }\text{，}\\ {{\sigma }^{\prime }}_{zz}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\mu }{(1-\mu )}\frac{\partial {u}_{x}}{\partial x}+\frac{\partial {u}_{z}}{\partial z}\right)\text{ }\text{，}\\ {\tau }_{xz}=-G\left(\frac{\partial {u}_{x}}{\partial z}+\frac{\partial {u}_{z}}{\partial x}\right)\text{ }。\end{array}\right\}$

(8)

式中：μ为泊松比；G为剪切模量。

1.2 边界条件

（1）海床表面的边界条件（z=0）

如图 1所示，海床表面的垂直有效应力和剪切应力趋于零，且海洋表面产生波浪运动引起的周期荷载：

${\sigma{zz}^{\prime}_} = {\tau _{xz}} = 0 \text{，} {p_{\text{f}}} = {p_{\text{b}}}(x,t) 。$

(9)

海床表面的波浪压力可以进一步表示为

${p_{\text{b}}}(x,t) = \sum\limits_{m = 1}^2 {{P_m}\cos } (mkx - m\omega t) \text{，}$

(10)

其中

${P_1} = \frac{{{\rho _{\text{f}}}gH}}{{2\cosh (kd)}} \text{，}$

(11)

${P_2} = \frac{3}{4}{\gamma _{\text{w}}}H\left( {\frac{{{\text{π }}H}}{L}} \right)\frac{1}{{\sinh (2kd)}}\left[ {\frac{1}{{{{\sinh }^2}(kd)}} - \frac{1}{3}} \right] 。$

(12)

式中： ${p_{\text{b}}}$ 为海床表面的周期荷载；m=1, 2，P₁和P₂分别为海床表面的一阶和二阶压力；k=2π/L为波数，由波长L确定；ω=2π/T为波的角频率，由周期T确定。

（2）海床底面的边界条件（z=-h）

如图 1所示，假定有限厚度海床的底面为刚性、不透水基岩，底部边界条件可以表示为

${u_x} = {u_z} = 0 \text{，} {{\partial p} \mathord{\left/ {\vphantom {{\partial p} {\partial z}}} \right. } {\partial z}} = 0 。$

(13)

2. 海床动力响应求解

由于作用在海床表面的波浪荷载呈现周期性变化，因此未知场量同样随着波数k与角频率ω的变化产生周期性变化，可以基于Madsen^[14]求解过程中使用的复变量对该问题进行求解。同时，考虑到二阶Stokes波由不同阶次的波浪组成，因此可以将所有场量改写成由不同阶次求和的复变量形式：

$\left\{ {\begin{array}{*{20}{c}} {{u_x}(x,z,t)} \\ {{u_z}(x,z,t)} \\ {p(x,z,t)} \end{array}} \right\} = \sum\limits_{m = 1}^2 {{Re} \left\{ {\begin{array}{*{20}{c}} {U_x^{\left( m \right)}(z)} \\ {U_z^{\left( m \right)}(z)} \\ {{P^{\left( m \right)}}(z)} \end{array}} \right\}} {{\text{e}}^{{\text{i}}m\left( {kx - \omega t} \right)}} 。$

(14)

式中：i=(-1)^1/2为虚数；Re代表对结果取实部； $U_x^{\left( m \right)}$ 、 $U_z^{\left( m \right)}$ 和 ${P^{\left( m \right)}}$ 为波浪荷载引起的海床动力响应的幅值。值得注意的是，由于将场量拓展到复数域进行求解，因此仅取最后得到结果的实部进行分析。

波浪作用在海床表面的荷载在复数域可以表示为

${p_{\text{b}}}(x,t) = \sum\limits_{m = 1}^2 {{P_m}{{\text{e}}^{{\text{i}}m\left( {kx - \omega t} \right)}}} 。$

(15)

通过将场量表示为复变量形式，式（5），（6）和（8）可以分别改写为

$\frac{{{\partial ^2}{P^{\left( m \right)}}}}{{\partial {z^2}}} - \left( {{m^2}{k^2} - \frac{{{\gamma _{\text{w}}}n\beta }}{K}{\text{i}}m\omega } \right){P^{\left( m \right)}}$

$\begin{array}{*{20}{c}} {}&{}&{}&{} \end{array} = - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}}}{K}\left( {{\text{i}}mkU_x^{\left( m \right)} + \frac{{\partial U_z^{\left( m \right)}}}{{\partial z}}} \right) \text{，}$

(16)

${P^{\left( m \right)}} = - \left( {{{\sigma^{\prime}}_{xx}} + \frac{1}{{{\text{i}}mk}}\frac{{\partial {\tau _{xz}}}}{{\partial z}}} \right) \text{，}$

(17)

$\left.\begin{array}{l}{{\sigma }^{\prime }}_{xx}=-\frac{2(1-\mu )G}{1-2\mu }\left(\text{i}mk{U}_{x}^{\left(m\right)}+\frac{\mu }{1-\mu }\frac{\partial {U}_{z}^{\left(m\right)}}{\partial z}\right)\text{ }\text{，}\\ {{\sigma }^{\prime }}_{zz}=-\frac{2(1-\mu )G}{1-2\mu }\left(\frac{\mu }{1-\mu }\text{i}mk{U}_{x}^{\left(m\right)}+\frac{\partial {U}_{z}^{\left(m\right)}}{\partial z}\right)\text{ }\text{，}\\ {\tau }_{xz}=-G\left(\frac{\partial {U}_{x}^{\left(m\right)}}{\partial z}+\text{i}mk{U}_{z}^{\left(m\right)}\right)\text{ }。\end{array}\right\}$

(18)

将式（17）代入式（7）和式（16）可以分别得到

$\frac{{\partial {{\sigma^{\prime}}_{zz}}}}{{\partial z}} + {\text{i}}mk{\tau _{xz}} = \frac{{\partial {{\sigma^{\prime}}_{xx}}}}{{\partial z}} + \frac{1}{{{\text{i}}mk}}\frac{{{\partial ^2}{\tau _{xz}}}}{{\partial {z^2}}} \text{，}$

(19)

$\frac{{{\partial ^2}{{\sigma^{\prime}}_{xx}}}}{{\partial {z^2}}} + \frac{1}{{{\text{i}}mk}}\frac{{{\partial ^3}{\tau _{xz}}}}{{\partial {z^3}}} - \left( {{m^2}{k^2} - \frac{{{\gamma _{\text{w}}}n\beta }}{K}{\text{i}}m\omega } \right)\left( {{{\sigma^{\prime}}_{xx}} + \frac{1}{{{\text{i}}mk}}\frac{{\partial {\tau _{xz}}}}{{\partial z}}} \right) \\ \;\;\;\;\;\;\;\;\;\; = \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}}}{K}\left( {{\text{i}}mkU_x^{(m)} + \frac{{\partial U_z^{(m)}}}{{\partial z}}} \right) 。$

(20)

将多孔弹性介质本构关系式（18）代入式（19），可以得到

$\frac{{{\partial ^2}U_z^{(m)}}}{{\partial {z^2}}} - {m^2}{k^2}U_z^{(m)} = \frac{1}{{{\text{i}}mk}}\left( {\frac{{{\partial ^3}U_x^{(m)}}}{{\partial {z^3}}} - {m^2}{k^2}\frac{{\partial U_x^{(m)}}}{{\partial z}}} \right) 。$

(21)

将式（7）代入式（20），可以得到

$\begin{array}{l} \left\{ {\frac{{{\partial ^3}U_z^{\left( m \right)}}}{{\partial {z^3}}} - {m^2}{k^2}\frac{{\partial U_z^{\left( m \right)}}}{{\partial z}} + } \right. \\ \left. {\frac{{{\text{i}}m\omega {\gamma _{\text{w}}}\left[ {n\beta + (1 - 2\mu )/G} \right]}}{K}\frac{{\partial U_z^{\left( m \right)}}}{{\partial z}}} \right\} \\ = - \frac{1}{{{\text{i}}mk}}\left\{ {(1 - 2\mu )\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} - \left[ {2(1 - \mu ){m^2}{k^2} + (1 - 2\mu ) \cdot } \right.} \right. \\ {\text{ }}\left. {\left( {{m^2}{k^2} - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}n\beta }}{K}} \right)} \right]\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} + \\ {m^2}{k^2}\left[ {2(1 - \mu ){m^2}{k^2} - {\text{i}}m\omega {\gamma _{\text{w}}} \cdot } \right. \\ \end{array}$

$\left. {\left. {\frac{{2n\beta (1 - \mu ) + (1 - 2\mu )/G}}{K}} \right]U_x^{\left( m \right)}} \right\} 。$

(22)

将式（21）代入式（22），可以得到

$\begin{array}{l} \frac{{\partial U_z^{\left( m \right)}}}{{\partial z}} = \frac{{2(1 - \mu )K}}{{{m^2}\omega k{\gamma _{\text{w}}}\left[ {n\beta + (1 - 2\mu )/G} \right]}} \times \\ \;\;\;\;\;\;\;\left\{ {\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} - \left[ {2{m^2}{k^2} - \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}n\beta }}{K}\frac{{(1 - 2\mu )}}{{2(1 - \mu )}}} \right]\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} + } \right. \end{array} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{m^2}{k^2}({m^2}{k^2} - {\kappa ^2})U_x^{\left( m \right)}\left. \begin{gathered} \\ \\ \end{gathered} \right\} 。$

(23)

式中：

${\kappa ^2} = \frac{{{\text{i}}m\omega {\gamma _{\text{w}}}\left( {n\beta + \frac{{(1 - 2\mu )}}{{2G(1 - \mu )}}} \right)}}{K} 。$

(24)

将式（23）代入式（22），可以得到

$\frac{{{\partial ^6}U_x^{\left( m \right)}}}{{\partial {z^6}}} - {m^2}{k^2}\left( {3 - \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)\frac{{{\partial ^4}U_x^{\left( m \right)}}}{{\partial {z^4}}} +$

${m^4}{k^4}\left( {3 - 2\frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)\frac{{{\partial ^2}U_x^{\left( m \right)}}}{{\partial {z^2}}} - {m^6}{k^6}\left( {1 - \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right)U_x^{\left( m \right)} = 0 。$

(25)

式（25）为六阶偏微分方程，该方程的解可以假设为

$U_x^{\left( m \right)} = A{{\text{e}}^{\gamma mkz}} 。$

(26)

将所假设的解式（26）代入式（25）中，可以得到

$({\gamma ^2} - 1)({\gamma ^2} - 1)\left( {{\gamma ^2} - 1 + \frac{{{\kappa ^2}}}{{{m^2}{k^2}}}} \right) = 0 。$

(27)

由式（27），可以得到 $U_x^{\left( m \right)}$ 的解析解为

${U}_{x}^{\left(m\right)}=({A}_{1}+{A}_{2}z){\text{e}}^{mkz}+({A}_{3}+{A}_{4}z){\text{e}}^{-mkz}+{A}_{5}{\text{e}}^{{k}_{*}mkz}+{A}_{6}{\text{e}}^{-{k}_{*}mkz}。$

(28)

式中： ${k_*}$ =[1-κ²/(mk)²]。

考虑到后续计算中矩阵求逆的正指数溢出问题，通过将式（28）重构为

$U_x^{\left( m \right)} = ({A_1} + {A_2}z){{\text{e}}^{mkz}} + {A_5}{{\text{e}}^{{k_*}mkz}} + ({A_3} + {A_4}z){{\text{e}}^{ - mk\left( {z + h} \right)}} +$

${A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}} 。$

(29)

将式（29）代入式（23），并对z求积分，可以得到 $U_z^{\left( m \right)}$ 的解析解为

$U_z^{\left( m \right)} = - {\text{i}}\left\{ {\left( {{A_1} + {A_2}z - \frac{{1 + 2\delta }}{{mk}}{A_2}} \right){{\text{e}}^{mkz}}{\text{ + }}{k_*}{A_5}{{\text{e}}^{{k_*}mkz}} - } \right.$

$\left. {\left( {{A_3} + {A_4}z + \frac{{1 + 2\delta }}{{mk}}{A_4}} \right){{\text{e}}^{ - mk\left( {z + h} \right)}} - {k_*}{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} 。$

(30)

式中：

$\delta = \frac{{(1 - 2\mu )n\beta }}{{n\beta + {{(1 - 2\mu )} \mathord{\left/ {\vphantom {{(1 - 2\mu )} G}} \right. } G}}} 。$

(31)

同理可得孔隙水压力和有效应力的解析解为

${P^{\left( m \right)}} = 2G{\text{i}}\left\{ {\frac{{(\delta + 2\mu - 1)}}{{1 - 2\mu }}{A_2}{{\text{e}}^{mkz}} + \left[ {\frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }}} \right]mk{A_5}{{\text{e}}^{{k_*}mkz}} + } \right. \\ \left. {\frac{{1 - 2\mu - \delta }}{{1 - 2\mu }}{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(32)

$\sigma_{x x}^{\prime(m)} = - 2G{\text{i}}\left\{ {mk{A_1}{{\text{e}}^{mkz}} + \frac{{(1 - 2\mu )mkz + 2\mu \delta }}{{1 - 2\mu }}{A_2}{{\text{e}}^{mkz}} + } \right.$

${\text{ }}\frac{{1 - \mu - \mu k_*^2}}{{1 - 2\mu }}mk{A_5}{{\text{e}}^{{k_*}mkz}} + mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} +$

${\text{ }}\frac{{(1 - 2\mu )mkz - 2\mu \delta }}{{1 - 2\mu }}{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} +$

$\left. {\frac{{1 - \mu - \mu k_*^2}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(33)

$\begin{gathered} \sigma_{z z}^{\prime(m)} = 2G{\text{i}}\left\{ {mk{A_1}{{\text{e}}^{mkz}} - \left[ { - mkz + \frac{{2(1 - \mu )}}{{1 - 2\mu }}\delta } \right]{A_2}{{\text{e}}^{mkz}} - } \right. \\ {\text{ }}\frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }}mk{A_5}{{\text{e}}^{{k_*}mkz}} + mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \\ {\text{ }}\left[ {mkz + \frac{{2\delta (1 - \mu )}}{{1 - 2\mu }}} \right]{A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} - \\ \end{gathered}$

$\left. {\frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} \text{，}$

(34)

$\begin{gathered} \tau _{xz}^{\left( m \right)} = 2G\left\{ { - mk{A_1}{{\text{e}}^{mkz}} - (mkz - \delta ){A_2}{{\text{e}}^{mkz}} - {k_*}mk{A_5}{{\text{e}}^{{k_*}mkz}}} \right. + \\ {\text{ }}mk{A_3}{{\text{e}}^{ - mk\left( {z + h} \right)}} + (mkz + \delta ){A_4}{{\text{e}}^{ - mk\left( {z + h} \right)}} + \end{gathered} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \left. {{k_*}mk{A_6}{{\text{e}}^{ - {k_*}mk\left( {z + h} \right)}}} \right\} 。$

(35)

式中：未知参数A₁~A₆可以代入边界条件（式（9），（13））确定，如下式所示：

$\left[\begin{array}{c}{A}_{1}\\ {A}_{2}\\ {A}_{3}\\ {A}_{4}\\ {A}_{5}\\ {A}_{6}\end{array}\right]={\left[\begin{array}{cccccc}{B}_{1}& -{B}_{2}& {B}_{1}{B}_{5}& {B}_{2}{B}_{5}& -{B}_{3}{B}_{1}& -{B}_{3}{B}_{1}{B}_{10}\\ -{B}_{1}& \delta & {B}_{1}{B}_{5}& \delta {B}_{5}& -{B}_{4}& {B}_{4}{B}_{10}\\ {B}_{5}& -{B}_{5}h& 1& -h& {B}_{10}& 1\\ {B}_{5}& -{B}_{5}h-{B}_{7}{B}_{5}& -1& h-{B}_{7}& {k}_{*}{B}_{10}& -{k}_{*}\\ 0& {B}_{8}{B}_{5}& 0& {B}_{8}& {B}_{9}{B}_{4}{B}_{10}& -{B}_{9}{B}_{4}\\ 0& 2G\text{i}{B}_{8}& 0& -2G\text{i}{B}_{8}{B}_{5}& 2G\text{i}{B}_{9}{B}_{1}& 2G\text{i}{B}_{9}{B}_{1}{B}_{10}\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ {p}_{b}\end{array}\right]。$

(36)

式中：

$\begin{array}{l} {B_1} = mk;{B_2} = \frac{{2\delta (1 - \mu )}}{{1 - 2\mu }};{B_3} = \frac{{\mu - (1 - \mu )k_*^2}}{{1 - 2\mu }};{B_4} = {k_*}mk; \\ {B_5} = {{\text{e}}^{ - mkh}};{B_6} = {{\text{e}}^{mkh}};{B_7} = \frac{{1 + 2\delta }}{{mk}};{B_8} = \frac{{\delta + 2\mu - 1}}{{1 - 2\mu }}; \end{array} \\ \;\;\;\; {B_9} = \frac{{(1 - \mu )(1 - k_*^2)}}{{1 - 2\mu }};{B_{10}} = {{\text{e}}^{ - {k_*}mkh}};{B_{11}} = {{\text{e}}^{{k_*}mkh}} 。$

(37)

3. 二阶Stokes波作用下海床动力响应分析

3.1 与既有解的对比

为了对比验证本文得到的二阶Stokes波作用下海床动力响应的解析解，本文分别给出了线性波和二阶Stokes波作用下的对比结果，对比结果中实线代表本文解析解，点代表既有解。首先通过选取与Hsu等^[13]一致的参数：T=15 s，d=70 m，L=324 m，h=25 m，G=10⁷ Pa，n=0.3，μ=1/3以及K=10^-2 m/s，将本文解析解与Hsu等^[13]得到的线性波作用下的解进行对比。由图 2可知本文解析解得到的孔隙水压力和有效应力曲线与Hsu等的曲线完全一致，从而可以证明本文解析解的正确与可靠。

图 2 线性波作用下本文解析解与Hsu等的解之间的对比

Figure 2. Comparison of seabed response between present solution and solution by Hsu et al. under linear waves

下载: 全尺寸图片幻灯片

进一步地，将本文解析解与Zhang等^[20]二阶Stokes波作用下的解进行对比如图 3所示。对比中选用的参数与Zhang等一致：T=12 s，H=8 m，d=0.125L，h=24 m，S_r=0.975，n=0.35，μ=0.4，K=10^-4 m/s以及G=5×10⁶ Pa。对比结果表明本文解析解与Zhang等^[20]的解有良好的匹配性，进一步表明了本文解的正确性和可靠性。

图 3 二阶斯托克斯波作用下本文解析解与Zhang等的解之间的对比

Figure 3. Comparison of seabed response between present solution and solution by Zhang et al. under second-order Stokes waves

下载: 全尺寸图片幻灯片

3.2 波浪和海床参数对动力响应的影响

由二阶Stokes波作用在海床表面周期压力的表达式（10）~（12）可知，波浪的周期T、波高H以及水深d对二阶Stokes波的影响较大；而在对海床动力响应的求解过程中可知，海床土体的剪切模量G和渗透系数K同样会对海床的动力响应产生影响。为了探究波浪参数和海床土体参数对海床动力响应的影响，给出了如表 1所示的波浪和海床土体参数并进行算例计算和分析。在下列算例计算和分析中，没有提及的参数的取值均来源于表 1。

表 1 多孔弹性海床和波浪的基本参数

Table 1. Basic poroelastic properties and wave conditions

波浪和海床土参数	数值
波的周期T	10 s
波高H	5 m
水深d	0.125L^*
海床厚度h	20 m
孔隙流体的实际体积弹性模量K_w	2×10⁹ Pa
孔隙流体密度ρ_f	1000 kg/m³
孔隙率n	0.35
泊松比μ	0.33
渗透系数K	10^-4 m/s
剪切模量G	1×10⁷ Pa
海床土饱和度S_r	0.975
^*表中L代表波长。

下载: 导出CSV

| 显示表格

（1）波浪参数的影响

首先对波浪参数对海床动力响应的影响进行研究。图 4是二阶Stokes波作用下波浪周期T对海床动力响应的影响。其中波浪周期T的取值分别为8，9，10，11，12 s。从图 4中可以发现不同的波浪周期T对海床动力响应随深度方向z/h的变化趋势是一致的，还可以从图中发现：在海床相同深度z/h时，最大孔隙水压力|p|随着波浪周期T的增加而增加；最大水平有效应力|σ′_xx|在较浅区域随着波浪周期T的增加而减小，在较深区域则呈现相反的趋势；最大垂直有效应力|σ′_zz|和最大剪切应力|τ_xz|在较浅区域随着波浪周期T的增加而减小，在较深区域则存在转折点。同时还可以看出在相对水深d/L不变的情况下，波浪周期T越小，海床的动力响应受到二阶Stokes波的二阶项的影响越明显。

图 4 二阶Stokes波作用下波浪周期T对海床动力响应的影响

Figure 4. Influences of wave period on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 5是二阶Stokes波作用下水深d对海床动力响应的影响。Le Méhauté^[22]认为当水深d小于0.1L时，Stokes波不能严格满足自由水面的边界条件，故水深d的取值分别为：0.1L，0.12L，0.125L，0.13L，0.14L。从图 5可以得知水深d对海床的动力响应有明显的影响，同时水深d还对二阶Stokes波的二阶项引起的海床动力响应的大小有相当显著的影响。与波浪周期T的影响类似，海床动力响应随深度方向z/h的变化趋势是一致的。最大孔隙水压力|p|、最大水平有效应力|σ′_xx|、最大垂直有效应力|σ′_zz|和最大剪切应力|τ_xz|随着水深d的增加与波浪周期T增加时均表现出相似的规律。当波浪周期T不变时，随着相对水深d/L较小时（例如d/L=0.1），海床的动力响应受二阶Stokes波二阶项的影响有显著的影响。

图 5 二阶Stokes波作用下水深d对海床动力响应的影响

Figure 5. Influences of water depth on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 6是二阶Stokes波作用下波高H对海床动力响应的影响。其中波高H的取值分别为3，4，5，6，7 m。由图 6可以得知波高H对海床动力响应的影响较小。在海床的相同深度z/h处，最大孔隙水压力|p|、最大垂直有效应|σ′_zz|和最大剪切应力|τ_xz|均随着波高H的增大而增大。同时，波高H对二阶Stokes波二阶项引起的海床动力响应基本没有影响。

图 6 二阶Stokes波作用下波高H对海床动力响应的影响

Figure 6. Influences of seabed thickness on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

（2）海床土体参数对海床动力响应的影响

图 7是二阶Stokes波作用下渗透系数K对海床动力响应的影响。其中渗透系数K的取值分别为10^-2，10^-3，10^-4，10^-5，10^-6 m/s。由图 7可得：渗透系数K对海床的动力响应有显著的影响。渗透系数K对海床动力响应的影响体现在沿深度方向z/h的变化率上，随着渗透系数K的减小在深度较浅处的最大孔隙水压力|p|减小速率越快，最大垂直有效应力|σ′_zz|增大速率越快，而在较深区域这些场量的变化速率与渗透系数K的相关性减小。当渗透系数K较大时，最大水平有效应力|σ′_xx|和最大剪切应力|τ_xz|受到的影响较大，渗透系数K较小时对最大剪切应力|τ_xz|基本没有影响，且对最大水平有效应力|σ′_xx|的影响较小。

图 7 二阶Stokes波作用下渗透系数K对海床动力响应的影响

Figure 7. Influences of permeability on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

图 8是二阶Stokes波作用下剪切模量G对海床动力响应的影响。其中剪切模量G的取值分别为1×10⁶，5×10⁶，1×10⁷，5×10⁷，1×10⁸ Pa。从图 8中可以得出：剪切模量G对海床的动力响应有显著的影响。在海床的相同深度z/h处，最大孔隙水压力|p|和最大垂直有效应力|σ′_zz|随着剪切模量G的增大而分别减小和增大。而最大剪切应力|τ_xz|在深度z/h大约为0.5处存在一个转折点，在转折点以上的相同深度z/h处随着剪切模量G的增大而增大，在转折点以下则呈现相反的变化趋势，此现象可能是由于竖向和水平向有效应力状态改变引起的，即在z/h=0.5以下部分，水平有效应力出现了快速的增大而竖直有效应力出现了减小，在竖向偏应力状态改变的情况下，剪切应力状态呈现了相反的变化趋势。

图 8 二阶Stokes波作用下剪切模量G对海床动力响应的影响

Figure 8. Influences of shear modulus on dynamic response of seabed under second-order Stokes waves

下载: 全尺寸图片幻灯片

4. 结论

本文研究引入了二阶Stokes波理论模拟波浪荷载，考虑了波浪非线性的影响，最终建立了二阶Stokes波作用下有限厚度海床动力响应的解析解。并将本文得到的解析解与线性波既有解进行了对比，验证了本文解的正确性和可靠性。最后分析了波浪和海床土体参数对海床动力响应的影响。通过本文的研究主要得出以下3点结论。

（1）在相对水深d/L确定，波浪周期T较小时，以及波浪周期T确定，相对水深较小时，海床的动力响应受到二阶Stokes波的二阶项影响均增加。

（2）波浪周期T、水深d和波高H均会对海床的动力响应产生影响，但波高H对二阶Stokes波的二阶项引起的海床动力响应基本没有影响。

（3）海床土体的渗透系数K和剪切模量G均会对海床的动力响应有显著的影响。渗透系数K对海床动力响应的影响表现为沿深度方向的变化速率，而剪切模量G对海床动力响应体现为幅值大小的影响。

图 1 模型槽试验布置图（改自Jardine等^[9]）

Figure 1. Schematic diagram of calibration chamber tests (after Jardine et al^[9])

下载: 全尺寸图片幻灯片

图 2 锥尖阻力q_c结果对比

Figure 2. Comparison of q_c

下载: 全尺寸图片幻灯片

图 3 桩打入完成后砂土平均有效正应力等值线图

Figure 3. Contours of mean normal stress of sand after pile penetration

下载: 全尺寸图片幻灯片

图 4 桩打入完成后砂土体积应变等值线图

Figure 4. Contours of volumetric strain of sand after pile penetration

下载: 全尺寸图片幻灯片

图 5 桩打入完成后砂土状态参数等值线图

Figure 5. Contours of state parameter of sand after pile penetration

下载: 全尺寸图片幻灯片

图 6 桩打入过程中距桩轴r/R = 2处归一化径向应力 ${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$ 结果对比

Comparison of ${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$ at r/R = 2 at driving stage

下载: 全尺寸图片幻灯片

图 7 桩打入过程中距桩轴r/R = 2处归一化径向稳态应力 ${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$ 结果对比

Comparison of ${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$ at r/R=2 at pause stage

下载: 全尺寸图片幻灯片

图 8 桩打入过程中距桩轴r/R = 2和3处归一化径向应力随h/R变化结果对比

Figure 8. Comparison of normalized radial stresses at h/R= 2 or 3

下载: 全尺寸图片幻灯片

图 9 桩打入完成后桩尖上方区域归一化径向稳态应力随r/R变化结果对比

Figure 9. Comparison of normalized radial stresses at pause stage with r/R after pile penetration

下载: 全尺寸图片幻灯片

表 1 状态相关莫尔-库仑模型参数

Table 1 Sand parameters adopted for state-dependent Mohr- Coulomb model

弹性参数临界状态参数状态相关参数

E = 60 MPa φ_c = 30° n_d = 1.70

$\nu$ = 0.3 e_Γ = 0.9 n_b = 2.30

λ_c = 0.019

ξ = 0.23

下载: 导出CSV

参考文献(29)

[1]	LEHANE B M, JARDINE R J, BOND A J, et al. Mechanisms of shaft friction in sand from instrumented pile tests[J]. Journal of Geotechnical Engineering, 1993, 119(1): 19-35. doi: 10.1061/(ASCE)0733-9410(1993)119:1(19)
[2]	CHOW F C M. Investigations into Displacement Pile Behaviour for Offshore Foundations[D]. London, UK: University of London, 1997.
[3]	张明义, 邓安福. 预制桩静力贯入层状地基的试验研究[J]. 岩土工程学报, 2000, 22(4): 490-492. doi: 10.3321/j.issn:1000-4548.2000.04.021 ZHANG Mingyi, DENG Anfu. Experimental study on static penetration of precast pile into layered foundation[J]. Chinese Journal of Geotechnical Engineering, 2000, 22(4): 490-492. (in Chinese) doi: 10.3321/j.issn:1000-4548.2000.04.021
[4]	GAVIN K G, LEHANE B M. The shaft capacity of pipe piles in sand[J]. Canadian Geotechnical Journal, 2003, 40(1): 36-45. doi: 10.1139/t02-093
[5]	王浩, 周健, 邓志辉. 砂土中桩端阻力随位移发挥的内在机理研究[J]. 岩土工程学报, 2006, 28(5): 587-593. doi: 10.3321/j.issn:1000-4548.2006.05.008 WANG Hao, ZHOU Jian, DENG Zhihui. Mobilization of toe resistance of piles with local displacement in sand[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(5): 587-593. (in Chinese) doi: 10.3321/j.issn:1000-4548.2006.05.008
[6]	周健, 邓益兵, 叶建忠, 等. 砂土中静压桩沉桩过程试验研究与颗粒流模拟[J]. 岩土工程学报, 2009, 31(4): 501-507. doi: 10.3321/j.issn:1000-4548.2009.04.002 ZHOU Jian, DENG Yibing, YE Jianzhong, et al. Experimental and numerical analysis of jacked piles during installation in sand[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(4): 501-507. (in Chinese) doi: 10.3321/j.issn:1000-4548.2009.04.002
[7]	YANG Z X, JARDINE R J, ZHU B T, et al. Sand grain crushing and interface shearing during displacement pile installation in sand[J]. Géotechnique, 2010, 60(6): 469-482. doi: 10.1680/geot.2010.60.6.469
[8]	JARDINE R J, ZHU B T, FORAY P, et al. Measurement of stresses around closed-ended displacement piles in sand[J]. Géotechnique, 2013, 63(1): 1-17. doi: 10.1680/geot.9.P.137
[9]	JARDINE R J, ZHU B T, FORAY P, et al. Interpretation of stress measurements made around closed-ended displacement piles in sand[J]. Géotechnique, 2013, 63(8): 613-627. doi: 10.1680/geot.9.P.138
[10]	HENKE S, QIU G, GRABE J. A Coupled Eulerian- Lagrangian Approach to Solve Geotechnical Problems Involving Large Deformations[M]// Numerical Methods in Geotechnical Engineering. Boca Raton: CRC Press, 2010: 233-238.
[11]	王腾, 薛浩, 吴瑞. 黏土中静压管桩土塞机制研究[J]. 岩土力学, 2018, 39(12): 4335-4341, 4350. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201812006.htm WANG Teng, XUE Hao, WU Rui. Mechanism of soil plug for jacked pipe pile in clay[J]. Rock and Soil Mechanics, 2018, 39(12): 4335-4341, 4350. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201812006.htm
[12]	ZHANG C, NGUYEN G D, EINAV I. The end-bearing capacity of piles penetrating into crushable soils[J]. Géotechnique, 2013, 63(5): 341-354. doi: 10.1680/geot.11.P.117
[13]	ZHANG C, YANG Z X, NGUYEN G D, et al. Theoretical breakage mechanics and experimental assessment of stresses surrounding piles penetrating into dense silica sand[J]. Géotechnique Letters, 2014, 4(1): 11-16. doi: 10.1680/geolett.13.00075
[14]	YANG Z X, GAO Y Y, JARDINE R J, et al. Large deformation finite-element simulation of displacement-pile installation experiments in sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2020, 146(6): 04020044. doi: 10.1061/(ASCE)GT.1943-5606.0002271
[15]	周健, 陈小亮, 周凯敏, 等. 静压开口管桩沉桩过程模型试验及数值模拟[J]. 岩石力学与工程学报, 2010, 29(增刊2): 3839-3846. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2010S2057.htm ZHOU Jian, CHEN Xiaoliang, ZHOU Kaimin, et al. Model test and numerical simulation of driving process of open-ended jacked pipe piles[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(S2): 3839-3846. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2010S2057.htm
[16]	詹永祥, 姚海林, 董启朋, 等. 砂土中开口管桩沉桩过程的颗粒流模拟研究[J]. 岩土力学, 2013, 34(1): 283-289. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201301041.htm ZHAN Yongxiang, YAO Hailin, DONG Qipeng, et al. Study of process of open-ended pipe pile driven into sand soil by particle flow simulation[J]. Rock and Soil Mechanics, 2013, 34(1): 283-289. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201301041.htm
[17]	WANG J F, ZHAO B D. Discrete-continuum analysis of monotonic pile penetration in crushable sands[J]. Canadian Geotechnical Journal, 2014, 51(10): 1095-1110. doi: 10.1139/cgj-2013-0263
[18]	LIU S, WANG J F. Depth-independent cone penetration mechanism by a discrete element method (DEM)-based stress normalization approach[J]. Canadian Geotechnical Journal, 2016, 53(5): 871-883. doi: 10.1139/cgj-2015-0188
[19]	LORENZO R, DA CUNHA R P, CORDÃO NETO M P, et al. Numerical simulation of installation of jacked piles in sand using material point method[J]. Canadian Geotechnical Journal, 2018, 55(1): 131-146. doi: 10.1139/cgj-2016-0455
[20]	CIANTIA M, O'SULLIVAN C, JARDINE R J. Pile penetration in crushable soils: Insights from micromechanical modelling[C]//Proceedings of 17th European Conf on Soil Mechanics and Geotechnical Engineering (ECSMGE-2019). Alexandria, 2019.
[21]	BEEN K, JEFFERIES M G. A state parameter for sands[J]. Géotechnique, 1985, 35(2): 99-112. doi: 10.1680/geot.1985.35.2.99
[22]	LI X S, WANG Y. Linear representation of steady-state line for sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(12): 1215-1217. doi: 10.1061/(ASCE)1090-0241(1998)124:12(1215)
[23]	LI X S, DAFALIAS Y F. Dilatancy for cohesionless soils[J]. Géotechnique, 2000, 50(4): 449-460. doi: 10.1680/geot.2000.50.4.449
[24]	AGHAKOUCHAK A. Advanced Laboratory Studies to Explore the Axial Cyclic Behaviour of Driven Piles[D]. London: Imperial College London, 2015.
[25]	DANO C, HICHER P Y, TAILLIEZ S. Engineering properties of grouted sands[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(3): 328-338. doi: 10.1061/(ASCE)1090-0241(2004)130:3(328)
[26]	GAUDIN C, SCHNAID F, GARNIER J. Sand characterization by combined centrifuge and laboratory tests[J]. International Journal of Physical Modelling in Geotechnics, 2005, 5(1): 42-56. doi: 10.1680/ijpmg.2005.050104
[27]	ANDRIA-NTOANINA J, CANOU J. DUPLA J C. Caractérisation Mécanique Du Sable De Fontainebleau NE34 À L'Appareil Triaxial Sous Cisaillement Monotone[R]. Paris: Laboratoire Navier-Géotechnique (CERMES, ENPC/LCPC), 2010. (ANDRIA-NTOANINA J, CANOU J. DUPLA J C. Mechanical Characterization of the Fontainebleau NE34 Sand in Monotonic Triaxial Shear Tests[R]. Paris: Navier Laboration-Geotechnical Engineering, 2010. (in French))
[28]	ALTUHAFI F N, JARDINE R J, GEORGIANNOU V N, et al. Effects of particle breakage and stress reversal on the behaviour of sand around displacement piles[J]. Géotechnique, 2018, 68(6): 546-555. doi: 10.1680/jgeot.17.P.117
[29]	ARSHAD M I. Experimental Study of the Displacements Caused by Cone Penetration in Sand[D]. West Lafayette: Purdue University, 2014.