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

孔壮壮,王栋,张民生,裴会敏

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

深基础在陆上和海洋工程中的应用日益增多，例如，近海导管架平台的桩基埋深可达150 m 以上，桩基底端的上覆有效应力水平约在1～2 MPa 之间[1]；中国南海第二轮天然气水合物试采的储层埋深为208 m，水平开采套管承受的上覆有效应力约为1.5 MPa[2]。工程实践中常采用原位试验确定土层的物理和力学性质，但大多针对中等或低应力水平土层，已有的原位试验数据解译方法也很少考虑高应力水平的影响。

静力触探（CPT）试验是海洋工程勘察中最常用的原位试验方法，可通过量测的锥尖阻力经验性地预测砂层的相对密实度和内摩擦角[3-5]。为建立预测公式，需要进行一定数量的室内标定罐试验：在标定罐中制备均匀砂样，对砂样施加预定的上覆压力或围压，量测施加压力后砂样的相对密实度；将静力触探仪贯入砂样中，记录贯入阻力随深度的变化；改变砂土的相对密实度和压力水平，进行多个试验。Lunne 等[6]总结了不同地区石英砂的标定罐试验和现场测试数据，建立起锥尖阻力qc和砂土相对密实度Dr值及平均有效应力p′之间的关系，但该依托试验的p′值不超过500 kPa。Ahmadi 等[7]开展了不同K0条件下的标定罐试验，考虑不同K0状态对锥尖阻力的影响，也提出过类似公式，但其研究中的p′值不超过200 kPa。目前公开报道的高应力条件下砂土标定罐试验和现场试验的数据极少。

标定罐试验砂样准备时间长、成本高，数值模拟可以作为有效辅助手段，以完善锥尖阻力与相对密实度之间的关系。然而，传统的小变形或有限变形框架内的有限元方法很难模拟CPT，原因是探头在土中的贯入导致周围土体网格的严重扭曲，使得计算无法进行。近年来，多种大变形数值方法已被用于追踪触探仪在砂土中的贯入过程，如任意拉格朗日—欧拉方法（Arbitrary Lagrangian Eulerian method，简称ALE）[8]、耦合欧拉—拉格朗日方法[9]和物质点法[10]。ALE 方法结合了拉格朗日和欧拉步骤，允许网格独立于材料移动，保持剩余的网格拓扑不变。与传统的拉格朗日分析相同，满足平衡方程、边界条件、外部载荷和接触条件，然后重新定位节点，并将所有变量从旧网格映射到新网格。因此，当发生大变形时，可以保持高质量的网格。Kouretzis等[11]采用ALE 法模拟砂土中静力触探试验，并与离心机试验结果相互验证，建立了相对密实度和锥尖阻力之间的关系，但其试验中p′值不超过340 kPa。

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

1 标定罐试验

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

图1 青岛某石英砂颗粒级配曲线Fig.1 Particle size distribution curve of a silica Qingdao 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，以保证砂样的下落高度基本不变。施加有效围压σ′3分别为0.5、1.0、2.0 MPa，采用等向固结，即静止土压力系数K0=1。提前在罐底布设4 个体积为200 mm3的带底环刀，用于试验后确定砂样的相对密实度。

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

表1 标定罐中砂样相对密实度Table 1 Relative density of sand sample in calibration chamber

2 有限元模拟

2.1 本构模型及参数标定

砂土在剪切时表现出的应变硬化或软化行为取决于砂土的相对密实度与应力水平。传统摩尔—库伦模型规定砂土内摩擦角和剪胀角为定值，不能描述砂土的软化和剪胀性。这里采用Hu 等[14]提出的MMC 模型，通过引进内摩擦角和剪胀角随累积塑性剪应变γ的变化来描述砂土的软化和剪胀性，其变化规律如图3 所示，图中φi、φp和φcv分别为初始、峰值和临界内摩擦角。因为摩尔—库伦类模型无法描述剪缩，所以假定当γ<γ1时，剪胀角为0。此后剪胀角迅速线性增大，在γ=γ2时达到峰值剪胀角ψp；剪胀角保持峰值至γ=γ3，然后线性减小，直至在临界应变γ4处减小到0。当砂样表现为硬化和剪缩性时，MMC 模型实质上退化为传统的摩尔—库伦模型：内摩擦角保持为φcv，剪胀角始终为0。对于一般石英砂，按照Hu 等[14]和Zheng 等[15]的建议，取φi=φcv、γ1=1%、γ2=1.2%、γ3=5%和γ4=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量级）时的最大剪切模量Gmax。将弯曲元量测设备安装在三轴仪上，在三轴仪中分别制备初始孔隙比为0.694 和0.764 的砂样，逐级施加100～1 500 kPa 的有效围压。每一级的固结完成后，首先测定排水量，用于计算固结后的孔隙比e，然后进行弯曲元试验，获得微小应变水平的Gmax。Gmax是e和p′的函数[16]。

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

图4 青岛某石英砂的小应变剪切模量Fig.4 Small-strain shear modulus of a silica Qingdao sand

弯曲元试验确定Gmax时的应变约为10-6量级，而MMC 模型中剪切模量G的应变范围大致在10-2量级。Loukidis 等[17]提出一种由Gmax估计G的方法。

式中：T为退化因子。按照Papadimitriou 等[18]和Pei等[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

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

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

为了后续的大量变动参数分析，需要建立砂土峰值内摩擦角或峰值剪胀角与相对密实度、平均有效应力之间的关系。参考经典Bolton 公式[20]，拟合表2 的三轴试验结果，得到式（4）～式（6）。

式中：IR为剪胀系数。由于φcv为定值，一旦确定Dr和p′后，即可根据式（4）～式（6）计算φp和ψp。

2.2 有限元模型

采用通用有限元软件ABAQUS 中的大变形模块ALE 模拟CPT 试验。建立的ALE 模型如图5 所示。参考Fan 等[21]和Mahmoodzadeh 等[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 大变形方法避免网格扭曲，从而实现静力触探在砂土中的连续贯入。为了与标定罐试验对比，对土体施加K0=1的压力。

3 ALE 有限元模型的验证

3.1 试验和模拟结果对比

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

图6 有限元模拟结果和标定罐试验结果对比Fig.6 Comparison of finite element simulation results and calibration chamber test results

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

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

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

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

3）对于Dr约为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；选择典型相对密实度Dr=40%、60% 和80%。利用式（2）计算G，利用式（5）和式（6）计算φp和ψp。

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

图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]总结的相对密实度和锥尖阻力之间的关系为

式中：参数C1、C2与土体性质有关；Qtn为归一化锥尖阻 力；m值 取0.6。Ahmadi 等[7]考虑不同K0状态对锥尖阻力的影响，建议m值取0.5。

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

按照式（9），由6 个标定罐试验得到的稳定锥尖阻力预测Dr，所得结果如图10 所示。从图10 中可以看出，式（9）能较好地预测高应力条件下（p′=0.5～2.0 MPa）砂土的Dr，其误差在±10%以内；图10 同时包含了Lunne 等[6]、Pournaghiazar 等[23]和Bolton 等[24]的标定罐和离心机试验结果，其应力水平为p′=25～500 kPa，相对密实度Dr=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 范围内的标定罐试验和大量大变形模拟，提出了由锥尖阻力预测相对密实度的经验关系式。