高压条件下石英砂中的静力触探试验及大变形模拟

孔壮壮，王栋，张民生，裴会敏; KONG Zhuangzhuang; WANG Dong; ZHANG Minsheng; PEI Huimin

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

高压条件下石英砂中的静力触探试验及大变形模拟 PDF

- ORCID：
孔壮壮
✉
- ORCID：
王栋
✉
- ORCID：
张民生
- ORCID：
裴会敏

中国海洋大学山东省海洋环境地质工程重点实验室，山东青岛 266100

中图分类号： TU411.93

最近更新：2023-07-09

DOI：10.11835/j.issn.2096-6717.2022.061

摘要

在深部土层中进行原位试验时，以往针对中低应力水平的数据解译公式不再适用。为探究砂土在高应力条件下的力学特性，利用自行研制的高压标定罐，改变围压水平，进行不同密实度石英砂中的静力触探试验，基于任意拉格朗日—欧拉格式的大变形有限元方法，模拟静力触探仪的完整贯入过程，有效避免了锥尖周围的网格扭曲。大变形模拟中采用修正摩尔—库伦模型描述砂土的内摩擦角和剪胀角随累积塑性剪应变的变化，并通过修正弯曲元试验，确定合理的土体剪切模量。标定罐试验和大变形有限元方法获得的锥尖阻力结果能很好地相互验证；在此基础上，建立高应力条件下静力触探锥尖阻力与砂土相对密实度之间的关系，与已有的低应力水平试验结果进行对比，证明建立的经验公式能通过锥尖阻力很好地预测砂土的相对密实度。

关键词

标定罐; 静力触探; 锥尖阻力; 深部土层; 力学特性; 砂土

深基础在陆上和海洋工程中的应用日益增多，例如，近海导管架平台的桩基埋深可达150 m以上，桩基底端的上覆有效应力水平约在1~2 MPa之间^[

1]；中国南海第二轮天然气水合物试采的储层埋深为208 m，水平开采套管承受的上覆有效应力约为1.5 MPa^{[参考文献 2

百度学术}2]。工程实践中常采用原位试验确定土层的物理和力学性质，但大多针对中等或低应力水平土层，已有的原位试验数据解译方法也很少考虑高应力水平的影响。

静力触探（CPT）试验是海洋工程勘察中最常用的原位试验方法，可通过量测的锥尖阻力经验性地预测砂层的相对密实度和内摩擦角^[

3-5]。为建立预测公式，需要进行一定数量的室内标定罐试验：在标定罐中制备均匀砂样，对砂样施加预定的上覆压力或围压，量测施加压力后砂样的相对密实度；将静力触探仪贯入砂样中，记录贯入阻力随深度的变化；改变砂土的相对密实度和压力水平，进行多个试验。Lunne等^{[参考文献 6

百度学术}6]总结了不同地区石英砂的标定罐试验和现场测试数据，建立起锥尖阻力q_c和砂土相对密实度D_r值及平均有效应力p′之间的关系，但该依托试验的p′值不超过500 kPa。Ahmadi等^{[参考文献 7

百度学术}7]开展了不同K₀条件下的标定罐试验，考虑不同K₀状态对锥尖阻力的影响，也提出过类似公式，但其研究中的p′值不超过200 kPa。目前公开报道的高应力条件下砂土标定罐试验和现场试验的数据极少。

标定罐试验砂样准备时间长、成本高，数值模拟可以作为有效辅助手段，以完善锥尖阻力与相对密实度之间的关系。然而，传统的小变形或有限变形框架内的有限元方法很难模拟CPT，原因是探头在土中的贯入导致周围土体网格的严重扭曲，使得计算无法进行。近年来，多种大变形数值方法已被用于追踪触探仪在砂土中的贯入过程，如任意拉格朗日—欧拉方法（Arbitrary Lagrangian Eulerian method，简称ALE）^[

8]、耦合欧拉—拉格朗日方法^{[参考文献 9

百度学术}9]和物质点法^{[参考文献 10

百度学术}10]。ALE方法结合了拉格朗日和欧拉步骤，允许网格独立于材料移动，保持剩余的网格拓扑不变。与传统的拉格朗日分析相同，满足平衡方程、边界条件、外部载荷和接触条件，然后重新定位节点，并将所有变量从旧网格映射到新网格。因此，当发生大变形时，可以保持高质量的网格。Kouretzis等^{[参考文献 11

百度学术}11]采用ALE法模拟砂土中静力触探试验，并与离心机试验结果相互验证，建立了相对密实度和锥尖阻力之间的关系，但其试验中p′值不超过340 kPa。

笔者采用高围压标定罐试验和ALE大变形有限元方法，进行砂土中静力触探贯入，建立高应力条件下锥尖阻力与石英砂相对密实度的关系。大变形模拟中采用修正摩尔—库伦模型（MMC）本构模型描述砂土内摩擦角和剪胀角随累积塑性剪应变的变化。

1 标定罐试验

采用某种青岛石英砂进行标定罐试验，砂的颗粒级配曲线如图1所示，平均粒径D₅₀=0.173 mm，不均匀系数C_u=1.70，最大孔隙比e_max=0.949，最小孔隙比e_min=0.490。

图1 青岛某石英砂颗粒级配曲线

Fig. 1 Particle size distribution curve of a silicaQingdao sand

早期标定罐侧壁是刚性的，当罐体直径不够大时，容易出现明显的边界效应。目前，标定罐都采用柔性双壁式标定罐，在刚性外壁和柔性内壁之间施加水头形成围压，从而有效消除刚性边界对贯入阻力的增强效应。笔者使用的高压标定罐（图2）也属于柔性双壁式：试样直径为600 mm、高度为750 mm，试样四周被柔性囊包裹，围压最大为5 MPa。采用砂雨法制备试样，罐口上部安装水平往复运动的砂斗，运动速度自动控制。砂斗底部开口宽度及其与砂面的相对高度为手动调节，试样的相对密实度取决于撒砂高度、开口宽度和运行速度^[

12]。

图2 中国海洋大学研制的高压标定罐

Fig. 2 High-pressure calibration chamber developed in Ocean University of China

砂中静力触探仪的贯入基本处于排水条件，Kluger等^[

13]进行了饱和砂和干砂中的CPT比较试验，认为当应力水平很低时，干砂样的锥尖阻力会略高于饱和砂样；但当平均有效应力p′增加到200 kPa时，干砂样和饱和砂样的锥尖阻力无明显差异。因此，采用干砂雨法制备不同相对密实度的试样：为了制备两种不同相对密实度的砂样，砂样下落高度分别采用0.5、1.3 m，每次撒砂的砂斗开口保持3 mm不变，运行速度控制在3 mm/s。在标定罐内部每隔0.1 m进行标记，撒砂厚度一旦达到标记线，即升高砂斗0.1 m，以保证砂样的下落高度基本不变。施加有效围压σ′₃分别为0.5、1.0、2.0 MPa，采用等向固结，即静止土压力系数K₀=1。提前在罐底布设4个体积为200 mm³的带底环刀，用于试验后确定砂样的相对密实度。

施加围压后砂样大致可分为两组，相对密实度分别约为46%和60%。共完成6个不同相对密实度和围压条件下的标定罐试验，对应的条件见表1。采用的CPT探头直径D=20 mm，锥角为60°，探头贯入速度为20 mm/s。

表1 标定罐中砂样相对密实度

Table 1 Relative density of sand sample in calibration chamber

σ′₃/MPa	D_r/%
0.5	45.2、59.8
1.0	46.2、58.1
2.0	46.6、62.1

2 有限元模拟

2.1 本构模型及参数标定

砂土在剪切时表现出的应变硬化或软化行为取决于砂土的相对密实度与应力水平。传统摩尔—库伦模型规定砂土内摩擦角和剪胀角为定值，不能描述砂土的软化和剪胀性。这里采用Hu等^[

14]提出的MMC模型，通过引进内摩擦角和剪胀角随累积塑性剪应变γ的变化来描述砂土的软化和剪胀性，其变化规律如图3所示，图中φ_i、φ_p和φ_cv分别为初始、峰值和临界内摩擦角。因为摩尔—库伦类模型无法描述剪缩，所以假定当γ < γ₁时，剪胀角为0。此后剪胀角迅速线性增大，在γ=γ₂时达到峰值剪胀角ψ_p；剪胀角保持峰值至γ=γ₃，然后线性减小，直至在临界应变γ₄处减小到0。当砂样表现为硬化和剪缩性时，MMC模型实质上退化为传统的摩尔—库伦模型：内摩擦角保持为φ_cv，剪胀角始终为0。对于一般石英砂，按照Hu等^{[参考文献 14

百度学术}14]和Zheng等^{[参考文献 15

百度学术}15]的建议，取φ_i=φ_cv、γ₁=1%、γ₂=1.2%、γ₃=5%和γ₄=15%。

图3 MMC模型中内摩擦角和剪胀角的变化

Fig. 3 Variation of friction and dilation angles of modified Mohr-Coulomb model

MMC模型包含的参数有：泊松比ν，取定值0.35；黏聚力c为0；剪切模量G；峰值内摩擦角φ_p；临界内摩擦角φ_cv和峰值剪胀角ψ_p。通过弯曲元试验推算G，利用三轴排水剪切试验获得φ_p、φ_cv和ψ_p。

2.1.1 弯曲元试验

弯曲元试验用于确定土体在微小应变（应变约为10^-6量级）时的最大剪切模量G_max。将弯曲元量测设备安装在三轴仪上，在三轴仪中分别制备初始孔隙比为0.694和0.764的砂样，逐级施加100~1 500 kPa的有效围压。每一级的固结完成后，首先测定排水量，用于计算固结后的孔隙比e，然后进行弯曲元试验，获得微小应变水平的G_max。G_max是e和p′的函数^[

16]。

图4为弯曲元试验得到的归一化最大剪切模量G_max/p_a随p′/p_a的变化，试验数据可拟合为

G_{m a x} = 900 \frac{{(2 - e)}^{2}}{1 + e} {(\frac{p'}{p_{a}})}^{0.6} \cdot p_{a}

（1）

图4 青岛某石英砂的小应变剪切模量

Fig. 4 Small-strain shear modulus of a silica Qingdao sand

弯曲元试验确定G_max时的应变约为10^-6量级，而MMC模型中剪切模量G的应变范围大致在10^-2量级。Loukidis等^[

17]提出一种由G_max估计G的方法。

G = G_max/T

（2）

式中：T为退化因子。按照Papadimitriou等^[

18]和Pei等^{[参考文献 19

百度学术}19]的建议，采用ALE方法试算CPT在标定罐中的贯入，根据试算反推，确定T = 5.5。

2.1.2 三轴排水试验

为了获取砂土的φ和ψ，进行6个三轴排水剪切试验，其相对密实度约为46%和60%，施加有效围压为0.5、1.0、2.0 MPa，剪切速率为0.5 %/min，具体方案和试验结果见表2。由于相对密实度约46%的砂土在围压1.0、2.0 MPa下表现为硬化和剪缩现象，故φ_p=φ_cv，ψ_p=0。

表2 三轴排水试验方案和结果

Table 2 Scheme and results of drained triaxial tests

σ′₃/MPa	D_r/%	φ_p/(°)	ψ_p/(°)
0.5	47.1	37.3	8.4
0.5	61.4	39.4	12.6
1.0	45.7	34.4	0.0
1.0	58.6	37.0	8.1
2.0	46.8	34.4	0.0
2.0	61.9	36.3	3.7

根据三轴排水剪切试验的应力—应变曲线峰值强度和临界强度可以分别求得砂土的φ_p和φ_cv。φ_cv基本不依赖应力水平与初始孔隙比，在6个试验中大致为定值34.4°。剪胀角利用式（3）根据体变—应变曲线得出。

s i n ψ = \frac{\frac{1}{k} \frac{d ε_{1}}{d ε_{3}} + 1}{\frac{1}{k} \frac{d ε_{1}}{d ε_{3}} - 1}

（3）

式中：dε₁和dε₃为主应变增量；三轴条件下k=2。

为了后续的大量变动参数分析，需要建立砂土峰值内摩擦角或峰值剪胀角与相对密实度、平均有效应力之间的关系。参考经典Bolton公式^[

20]，拟合表2的三轴试验结果，得到式（4）~式（6）。

I_R = D_r(10 - ln p′) - 1

（4）

φ_p - φ_cv = 3.67I_R

（5）

φ_p - φ_cv = 0.39ψ_p

（6）

式中：I_R为剪胀系数。由于φ_cv为定值，一旦确定D_r和p′后，即可根据式（4）~式（6）计算φ_p和ψ_p。

2.2 有限元模型

采用通用有限元软件ABAQUS中的大变形模块ALE模拟CPT试验。建立的ALE模型如图5所示。参考Fan等^[

21]和Mahmoodzadeh等^{[参考文献 22

百度学术}22]的成果，在锥尖下人为设置直径为D/20的刚性小管，此刚性小管与CPT一起移动，从而防止最左边的土体单元向对称轴的内部移动，有利于保证土体网格调整的有效性，同时不会影响计算结果。

图5 有限元轴对称模型

Fig. 5 Finite element axisymmetric model

按照二维轴对称计算，土体区域的宽和高分别是28D和42D。为了减少计算量，实施网格调整的ALE区域宽10D，根据前期试算，ALE区域中单元大小设置为D/8，这样能在保证计算精度的同时提高计算效率。ALE区域需要保持节点位置调整后网格的质量，以避免单元畸变，选用ABAQUS内置几何加强形式的网格平滑方法^[

11]。

采用硬接触算法模拟静力触探仪与土体之间的相互作用。假定二者交界面上的摩擦应力为0，即忽略贯入过程中探杆上的侧摩阻力，这是因为砂土中探杆上的侧摩阻力远小于锥尖阻力。

设置CPT贯入速度为20 mm/s进行匀速贯入，贯入深度为0~20D，采用ALE大变形方法避免网格扭曲，从而实现静力触探在砂土中的连续贯入。为了与标定罐试验对比，对土体施加K₀=1的压力。

3 ALE有限元模型的验证

3.1 试验和模拟结果对比

图6对比了有限元模拟和标定罐试验获得的锥尖阻力，其中H代表锥尖深度。从标定罐试验结果曲线可以看到，每条锥尖阻力曲线都分为3个阶段：1）H=0~1.5D时，q_c随深度几乎呈线性增加，增长斜率和围压与相对密实度呈正相关；2）H=1.5D~10D时，q_c的增长斜率变缓并逐渐趋于稳定；3）H>10D时，q_c略微减小或基本保持不变。q_c略微减小的原因可能是在试样上放置的顶盖使得浅层试样的密实度略高于深部。

（a） D_r=45%~47%

（b） D_r=58%~63%

图6 有限元模拟结果和标定罐试验结果对比

Fig. 6 Comparison of finite element simulation results and calibration chamber test results

高应力条件下锥尖阻力仍满足随相对密实度和围压增大而增大的现象。对比D_r约为46%的砂样，围压1.0 MPa或2.0 MPa时的锥尖阻力大约为0.5 MPa时的1.76倍和2.66倍；对比D_r约为60%的砂样，围压1.0 MPa或2.0 MPa时的锥尖阻力大约为0.5 MPa时的1.58倍和2.66倍。这表明对于不同相对密实度的砂样，高应力条件下的锥尖阻力仍依赖围压水平。如果围压相同，D_r约为60%砂样的锥尖阻力是D_r约为46%砂样的1.17~1.3倍。

整体来说，ALE方法能够较好地模拟高围压石英砂中的CPT贯入。从图6中可以看出：

1）对于不同相对密实度的砂样，当贯入深度H>10D时，有限元计算得到的q_c基本不再变化，存在一定程度的噪音波动，但在可接受范围内。

2）对于D_r约为46%的砂样，试验与有限元结果基本吻合。

3）对于D_r约为60%的砂样，H<10D时的试验结果比有限元结果约高15%，但二者给出的H>10D时的稳态锥尖阻力很接近。

3.2 剪切模量敏感性分析

为了进一步探讨获得的剪切模量G对锥尖阻力的影响，保持峰值摩擦角和剪胀角不变，人为变化G值。ALE得到的锥尖阻力和剪切模量的变化曲线如图7所示。随着G的增加，锥尖阻力逐渐增大，G=200 MPa时对应的锥尖阻力是G=50 MPa时的2.7倍，因此，数值模拟中不能直接采用利用弯曲元获得的低应变条件下的剪切模量，而应该按照式（2）计算。

图7 锥尖阻力随剪切模量变化曲线

Fig. 7 Changing curve of cone tip resistance with shear modulus

4 相对密实度和锥尖阻力的关系

为了建立高应力下砂土锥尖阻力和相对密实度之间的定量关系，进行变动参数分析：围压分别为0.5、0.8、1.0、1.3、1.6、2.0 MPa；选择典型相对密实度D_r = 40%、60%和80%。利用式（2）计算G，利用式（5）和式（6）计算φ_p和ψ_p。

图8给出了0.8、1.3 MPa围压下不同相对密实度的锥尖阻力曲线。对比围压为0.8、1.3 MPa的砂样，相对密实度80%和60%时的q_c大约是40%时的1.22倍和1.56倍左右。当D_r=40%、60%和80%时，围压1.3 MPa时，q_c约为0.8 MPa时的1.35倍。图9给出了CPT贯入过程中的土体位移矢量图，大变形条件下土体位移场可近似认为是破坏场。可以看出，锥面附近土体的位移大致沿锥面外法线方向，CPT探头下的土体处于三轴压状态，被径向挤出。

（a）围压为0.8 MPa

（b）围压为1.3 MPa

图8 不同围压下锥尖阻力曲线

Fig. 8 Cone tip resistance curves under different confining pressures

图9 CPT贯入过程土体位移矢量图

Fig. 9 Vector diagram of soil displacement during CPT penetration

Lunne等^[

6]总结的相对密实度和锥尖阻力之间的关系为

D_r=C₁ln Q_tn-C₂

（7）

Q_tn=（q_c/p_a）/（p′/p_a）^m

（8）

式中：参数C₁、C₂与土体性质有关；Q_tn为归一化锥尖阻力；m值取0.6。Ahmadi等^[

7]考虑不同K₀状态对锥尖阻力的影响，建议m值取0.5。

取贯入深度H=10D作为稳定时的锥尖阻力进行分析。总结18个工况的ALE模拟结果发现，在p′=0.5~2.0 MPa的高应力水平下，仍可以按照式（7）由锥尖阻力预测相对密实度，但式（7）的参数拟合为：C₁=0.58、C₂=1.91、m=0.64。即

D_r=0.58ln Q_tn-1.91

（9）

Q_tn=（q_c/p_a）/（p′/p_a）^0.64

（10）

按照式（9），由6个标定罐试验得到的稳定锥尖阻力预测D_r，所得结果如图10所示。从图10中可以看出，式（9）能较好地预测高应力条件下（p′=0.5~2.0 MPa）砂土的D_r，其误差在±10%以内；图10同时包含了Lunne等^[

6]、Pournaghiazar等^{[参考文献 23

百度学术}23]和Bolton等^{[参考文献 24

百度学术}24]的标定罐和离心机试验结果，其应力水平为p′=25~500 kPa，相对密实度D_r=20%~81%。可以看出，用于低应力水平时，对高密实度砂土的预测偏差相对较大，但其离散程度也能控制在±30%以内。

图10 砂土相对密实度预测值和实际值对比

Fig. 10 Comparison of the predicted and measured values of relative density of sand

5 结论

开展了高应力条件下的静力触探标定罐试验和大变形有限元模拟。大变形模拟中采用的本构模型考虑了砂土内摩擦角和剪胀角随累积塑性剪应变的变化。标定罐试验和大变形模拟能够很好地相互验证，主要结论如下：

1）土体剪切模量取值严重影响锥尖阻力，建议采用式（2）修正弯曲元试验结果。

2）静力触探的锥尖贯入到10倍的直径时，锥尖阻力达到稳定。锥尖阻力随相对密实度和围压的增大而增大。

3）结合有效应力水平在0.5~2.0 MPa范围内的标定罐试验和大量大变形模拟，提出了由锥尖阻力预测相对密实度的经验关系式。

参考文献

ELSAYED T, EL-SHAIB M, GBR K. Reliability of fixed offshore jacket platform against earthquake collapse [J]. Ships and Offshore Structures, 2016, 11(2): 167-181. [百度学术]

YE J L, QIN X W, XIE W W, et al. The second natural gas hydrate production test in the South China Sea [J]. China Geology, 2020, 3(2): 197-209. [百度学术]

ROURHE-NBURG L, JEFFERIES M G, CROOKS J H A, et al. The cone penetration test in sands: Part II, general inference of state [J]. Géotechnique, 1987, 37(3): 285-299. [百度学术]

沈小克, 蔡正银, 蔡国军. 原位测试技术与工程勘察应用[J]. 土木工程学报, 2016, 49(2): 98-120. [百度学术]

SHEN X K, CAI Z Y, CAI G J. Applications of in situ tests in site characterization and evaluation [J]. China Civil Engineering Journal, 2016, 49(2): 98-120. (in Chinese) [百度学术]

黄凯, 张明义, 白晓宇, 等. 基于单桥静力触探的静压桩沉桩阻力估算方法[J]. 土木与环境工程学报(中英文), 2019, 41(1): 55-61. [百度学术]

HUANG K, ZHANG M Y, BAI X Y, et al. Estimation of the penetration resistance of jacked pile based on the single bridge CPT [J]. Journal of Civil and Environmental Engineering, 2019, 41(1): 55-61. (in Chinese) [百度学术]

LUNNE T, POWELL J J M, ROBERTSON P K. Cone penetration testing in geotechnical practice [M]. Boca Raton: CRC Press, 2002. [百度学术]

AHMADI M M, KARAMBAKHSH P. K₀ determination of sand using CPT calibration chamber [C]//Proceedings of the 2nd International Symposium on cone penetration testing, California, 2010. [百度学术]

HUANG W, SHENG D, SLOAN S W, et al. Finite element analysis of cone penetration in cohesionless soil [J]. Computers and Geotechnics, 2004, 31(7): 517-528. [百度学术]

WANG D, BIENEN B, NAZEM M, et al. Large deformation finite element analyses in geotechnical engineering [J]. Computers and Geotechnics, 2015, 65: 104-114. [百度学术]

MARTINELLI M, GALAVI V. Investigation of the Material Point Method in the simulation of Cone Penetration Tests in dry sand [J]. Computers and Geotechnics, 2021, 130(7): 103923. [百度学术]

KOURETZIS G P, SHENG D, WANG D. Numerical simulation of cone penetration testing using a new critical state constitutive model for sand [J]. Computers and Geotechnics, 2014, 56: 50-60. [百度学术]

方浩, 段雪锋, 王永志, 等. 离心模型砂雨法制备技术研究与展望[J]. 世界地震工程, 2018, 34(4): 60-66. [百度学术]

FANG H, DUAN X F, WANG Y Z, et al. Research and prospect of preparation technique for centrifugal pluviation model [J]. World Earthquake Engineering, 2018, 34(4): 60-66. (in Chinese) [百度学术]

KLUGER M O, KREITER S, STÄHLER F T, et al. Cone penetration tests in dry and saturated Ticino sand [J]. Bulletin of Engineering Geology and the Environment, 2021, 80(5): 4079-4088. [百度学术]

HU P, WANG D, STANIER S A, et al. Assessing the punch-through hazard of a spudcan on sand overlying clay [J]. Géotechnique, 2015,65(11): 1-14. [百度学术]

ZHENG J B, HOSSAIN M S, WANG D. Numerical investigation of spudcan penetration in multi-layer deposits with an interbedded sand layer [J]. Géotechnique, 2017, 67(12): 1-17. [百度学术]

OZTOPRAK S, BOLTON M D. Stiffness of sands through a laboratory test database [J]. Géotechnique, 2013, 63(1): 54-70. [百度学术]

LOUKIDIS D, SALGADO R. Modeling sand response using two-surface plasticity [J]. Computers and Geotechnics, 2009, 36(1/2): 166-186. [百度学术]

PAPADIMITRIOU A G, BOUCKOVALAS G D, DAFALIAS Y F. Plasticity model for sand under small and large cyclic strains [J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001, 127(11): 973-983. [百度学术]

裴会敏, 王栋, 刘清秉. 钙质砂中静力触探试验的大变形有限元模拟[J].土木与环境工程学报(中英文), 2021, 43(1): 48-53. [百度学术]

PEI H M, WANG D, LIU Q B. Large deformation finite element analysis of cone penetration tests in calcareous sands [J]. Journal of Civil and Environmental Engineering, 2021, 43(1): 48-53. [百度学术]

BOLTON M D. The strength and dilatancy of sands [J]. Géotechnique, 1986, 36(1): 65-78. [百度学术]

FAN S, BIENEN B, RANDOLPH M F. Stability and efficiency studies in the numerical simulation of cone penetration in sand [J]. Géotechnique Letters, 2018, 8(1): 13-18. [百度学术]

MAHMOODZADEH H, RANDOLPH M F, WANG D. Numerical simulation of piezocone dissipation test in clays [J]. Géotechnique, 2014, 64(8): 657-666. [百度学术]

POURNAGHIAZAR M, RUSSELL A R, KHALILI N. The cone penetration test in unsaturated sands [J]. Géotechnique, 2013, 63(14): 1209-1220. [百度学术]

BOLTON M D, GUI M W, GARNIER J, et al. Centrifuge cone penetration tests in sand [J]. Géotechnique, 1999, 49(4): 543-552. [百度学术]