单桶和群桶基础有限元分析模型几何形状如图1．其中，D为桶径；L为桶长；t为桶壁厚；s为桶间距．单桶模型在模拟中是1个完整桶．为更清晰看到桶形基础与地基土尺寸，显示其一半进行标注，地基土深度为6D，地基土半径方向长度为5D，见图1（a）．地基土外圆周面只允许竖向变形，底面三向位移都约束，表面自由，如图1（b）．取群桶基础的1/4建立对称模型，桶形基础至对称面的距离为桶间距的1/2 ，即s/2，地基土深度为6D，地基土半径方向长度为5D+s/2．地基土两个竖直边界约束为沿该竖直面对称，具体为位移约束不允许沿该面法向发生位移，转动约束只允许沿该面法向转动，圆周面、底面和表面约束情况与单桶一致．单桶和群桶有限元模型具体边界尺寸通过边界影响分析确定．单桶和群桶模型均采用6面体8节点线性减缩积分单元（C3D8R）划分网格．桶壁及桶顶内侧与地基土均为绑定接触，桶壁外侧与地基土为摩擦接触，法向不可脱开.

软黏土地基某一深度处初始完整强度设置为s _u0=s _um+kz．其中， s _um为软黏土泥面强度； k为强度沿深度变化梯度值；z为软黏土深度．采用EINAV等^[17]提出的考虑应变软化和速率效应的修正Tresca屈服准则材料来模拟，软黏土地基某一深度处强度s_u可表达为

图1 数值模型几何形状（a）单桶；（b）群桶Fig. 1 Geometry of the numerical model of (a) single caisson and (b) caisson groups.

其中，δ_rem、ξ和ξ₉₅为软黏土应变软化参数，δ_rem为土体完全重塑后强度与初始强度的比值，即土体灵敏度S_t的倒数；ξ为土体累积塑性应变；ξ₉₅为土样强度降低到峰值强度5%时的累积剪应变，海洋软黏土的灵敏度一般为2～5，即δ_rem一般为0. 2～0. 5， ξ₉₅一般为10～50^[18]；μ为黏滞性系数，表示每对数周期土体强度的增加速率，取值范围一般为0. 05～0. 20^[19-20]；γ̇ref为参考剪切应变速率，一般取三轴试验中的典型应变速率，即3×10^-6 s^-1[21-22]；γ̇_max为土样在某时刻的最大剪切应变速率.

式（1）软黏土地基强度考虑应变软化和速率效应需要在ABAQUS中进行二次开发，此处通过场变量子程序USDFLD来实现，子程序相关变量解释可查看ABAQUS用户手册，并可通过建立一个单元体模型进行纯剪加载来验证子程序的准确性和可靠性．单元体模型尺寸为1m×1m×1m，底面固定， 4个侧面自由．利用位移控制方法在顶面施加水平向0. 01 m位移荷载，加载时间为1 s，固定增量为0. 001 s．单元体初始不排水抗剪强度s_u=1 kPa，弹性模量E=5 000s_u．为保证变形基本都为塑性变形，弹性模量选取比较大．根据应变软化和速率效应参数不同，考虑4种工况：理想土（工况1）、仅考虑应变软化（工况2）、仅考虑速率效应（工况3）、考虑应变软化和速率效应（工况4），如表1.

表1 子程序单元体验证工况Table 1 Unit cases of the subroutine verification

需要注意，实际现场软黏土灵敏度一般为2～5，ξ₉₅一般为10～50，单元体验证中为了使应变软化更加明显，假定灵敏度参数S_t=100，ξ₉₅=0. 05，二者并未根据实际常见值选取.

4种工况的数值模拟及式（1）计算结果如图2．其中，实线为数值模拟输出的不排水抗剪强度值，虚线为式（1）计算得到的不排水抗剪强度值；对于理想土，强度值为1. 00 kPa；只考虑速率效应，当剪切速率保持0. 01 m/s不变，根据式（1）计算出不排水抗剪强度为1. 35 kPa；若应变均为塑性剪应变，即0. 01，只考虑应变软化，最终强度为0. 55 kPa，然而由于存在弹性变形，实际上塑性剪应变为0. 009 94．依据此剪应变，由式（1）计算和数值模拟得到的最终强度均为0. 56 kPa，表明子程序准确有效.

以ZHU等^[22]中10 m桶径离心模型试验为例对模型有效性进行验证，此部分模拟中地基土强度参数、模型桶几何尺寸均按离心模型试验选取，弹性模量取E=300s_u，不排水条件体的应变为0，设置泊松比为0. 49，应变软化参数S_t和ξ₉₅分别取4. 2和25^[23]．黏滞性系数μ未给出，此处取0．模拟中上拔位移为0. 1倍桶径，由于率效应的存在，上拔速率设置为2 mm/s，与试验保持一致．

图2 模拟结果与计算结果对比Fig. 2 Comparison between simulation results and calculation results. Dashed lines are for calculation results. Solid lines are for simulation results. Black, green, red, and blue are for case 1, case 2, case 3, and case 4, respectively.

离心模型试验和数值模拟上拔荷载对比如图3．其中，上拔荷载V/（DLs _{u ，b} ）的量纲为1 ， V为上拔荷载，s_u，b为桶底处地基土不排水抗剪强度；w_c为桶形基础上拔距离．从图3可以看出，当考虑速率效应（即μ=0. 1）时，数值分析得到的上拔荷载峰值比试验结果大，不考虑速率效应时，数值计算与试验结果较为接近．这是因为与试验结果相比，数值计算中上拔荷载更早达到峰值，而试验中上拔荷载一直在增长，因此上拔至相同位置处，数值计算得到的上拔荷载要比试验结果高．图3中起始段数值计算和试验结果存在差别，是因为试验中桶基上拔，桶内土塞会迅速与地基土发生断裂，上拔荷载随之达到屈服段；数值分析中随着桶基上拔，需要桶底网格发生拉伸破坏，上拔荷载进入屈服，相比试验需要更大的上拔距离．因此，试验中桶形基础初始刚度明显大于数值分析．

图3 离心模型试验与数值模拟上拔荷载对比Fig. 3 Comparison of the uplift load between the centrifuge model test and the numerical model. Black lines are for the test results, in which the solid line represents the single caisson and the dashed line represents the caissons group, respectively. Color lines are for the simulation results, in which the red solid, blue solid and blue dashed represent single caisson with rate effect, without rate effect and caissons group without effect, respectively.

图4为应变软化对上拔荷载的影响．从图4可以看出，S_t越大或者ξ₉₅越小，上拔荷载也越小．根据式（1）可知，这是因为S_t越大或者ξ₉₅越小，相同塑性应变引起的软化效应也越明显．从图4还可见，上拔初始段上拔荷载基本重合，上拔至0. 02D左右上拔荷载才有所差异；S_t=10时，峰值荷载比未软化时（S_t=1）约低2%；定义上拔荷载最大值为峰值荷载，当ξ₉₅=5时，峰值荷载比ξ₉₅=50时约低5%，表明应变软化作用对上拔荷载影响较小．造成这一现象的原因有两点：上拔伊始，土体未到塑性屈服，上拔至0. 02D时，塑性应变仍然较小，此时应变软化最严重处强度降低4. 5%，见图5（a）；当上拔位移较大时，应变软化效应才明显，见图5 （b）．由于地基土塑性应变（应变软化）主要发生在桶端附近，影响范围很小，即使如图5（b）中软化效应最严重处强度下降20%，应变软化对上拔荷载影响亦不明显.

图4 应变软化对上拔荷载的影响Fig. 4 Influence of strain softening on the uplift load. Solid lines are for different sensitivities,in which black, red, green and blue represent the sensitivity of 1. 0, 2. 5, 4. 0, 10, respectively. Dashed lines are for different parameter ξ95, in which black, red, green and blue represent ξ95 of 5, 10, 25 and 50, respectively.

图5 应变软化使地基土强度降低（St=10） （a）wc/D=0. 02；（b）wc/D=0. 10 Fig. 5 Strength degradation of clay caused by strain softening (St=10) when (a) wc/D=0. 02, (b) wc/D=0. 10.

图6为率效应对上拔荷载的影响．由图6可见， μ或v越大，桶形基础上拔荷载越大．与μ=0相比，当μ=0. 05、0. 10和0. 20时，峰值上拔荷载分别增加约11%、21%和42%，荷载增加和μ增长呈线性关系．v每增加10倍，峰值荷载约增加8. 5%，上拔荷载与v的对数也存在线性关系．与应变软化影响不同，从图6可看出，对于不同μ和v，上拔荷载在上拔开始就出现差异，在整个上拔过程中速率效应对上拔荷载影响十分明显．这是由于地基土剪应变发生在整个上拔过程中，且和塑性剪应变主要集中在桶端附近不同，地基土剪应变范围桶底以下深度方向最大超过2D，桶周径向超过3D，如图7.

图6 率效应对上拔荷载的影响Fig. 6 Influence of rate effect on the uplift load. Solid lines are for different parameter μ,in which blue, green, red and black rep‐resent μ of 0, 0. 05, 0. 1 and 0. 2, respectively. Dashed lines are for different velocities,in which blue, green, red and black repre‐sent the velocities of 0. 2, 2. 0, 20 and 200, respectively.

桶形基础不排水上拔过程中，由于桶内负压吸力作用桶内土塞会随桶一起拔出，因此本模型中桶内壁与土的接触设置与MANA等^[11]一致，二者接触为绑定不允许脱开，桶外壁与地基土为摩擦接触.通过改变桶外壁接触粗糙度，进一步分析摩擦设置对桶形基础上拔响应的影响．图8为不同桶壁粗糙度对应的上拔荷载．其中，s_u，ave为不排水抗剪强度桶外壁深度方向平均值．由图8可见，显然桶壁粗糙度越大，桶外壁摩擦力越大，上拔荷载也越大．和光滑接触相比，桶壁粗糙度α=0. 2、0. 5和1. 0时，上拔荷载峰值分别增加7%、18%和36%，上拔荷载峰值和桶壁粗糙度呈线性关系.

图9为不同桶壁粗糙度归一化桶外壁摩擦力（A_o为桶外壁面积）．从图9可见，摩擦力变化基本与桶壁粗糙度变化线性增长关系．对于桶壁粗糙度α = 0. 2、 0. 5和1. 0的情况，摩擦力基本稳定在0. 2、0. 5和1. 0左右，此处处理后的摩擦力也就是桶壁粗糙度，因此数值结果是准确的.

图7 率效应使地基土强度增加（a）wc/D=0. 02；（b）wc/D=1. 00 Fig. 7 Soil strength enhancement caused by rate effect when (a) wc/D=0. 02, (b) wc/D=0. 10.

图8 桶外壁粗糙度对上拔荷载的影响Fig. 8 Influence of the caisson outer wall roughness on the uplift load. The black, red, green and blue lines are for parameter alpha of 0, 0. 2, 0. 5 and 1. 0, respectively.

图9中w_c/D=0. 04附近，存在数值下降后续又逐渐增加现象，这是因为随着桶基上拔，桶底网格发生拉伸畸形破坏，桶底端抗力基本不再增加，此时桶土间相对运动趋势减弱，桶土间摩擦力减小（w_c/D=0. 3～0. 5）．随着上拔的继续，摩擦力本应维持较为稳定，然而由于本研究考虑了速率效应和应变软化对土体强度影响，速率效应引起强度幅度大于应变软化强度降低幅度，故随着上拔继续，摩擦力会缓慢增加.

图9 归一化桶外壁摩擦力Fig. 9 Normalized caisson outer wall friction. The blue, green and red lines are for parameter alpha of 0. 2, 0. 5 and 1. 0, respectively.

图 10为上拔荷载减去桶外壁摩擦力后上拔荷载. 由图 10可见，不同桶壁粗糙度其图中荷载变化相当接近，只随桶壁粗糙度增大略有增加，上拔荷载增加基本与图9中摩擦力的增加一致，桶壁粗糙度变化引起的应变软化及率效应对上拔荷载影响可忽略.

图 10 减去桶外壁摩擦力后的上拔荷载Fig. 10 Uplift load minus the friction. The black line are for the single caisson. The wathet, blue, green, brown and red lines are for normalised caisson spacing of 0. 2, 0. 5, 1. 0, 1. 2 and 1. 5, respectively.

本研究分析了桶形基础几何尺寸对上拔响应的影响，分别选取桶长径比L/D = 0. 5、 1. 0、 2. 0、3. 0和5. 0．图 11为桶形基础长径比对上拔荷载的影响．从图 11可以看出，桶形基础长径比越小，峰值上拔荷载越大，初始刚度也越大，这是由于不排水上拔地基土为整体破坏，桶底以下土体强度对上拔响应影响明显．对于L/D=0. 5和2. 0的情况，上拔至0. 1D时，地基土应力应变范围分别在桶底以下2. 5D～3D处的深度（图 12），此影响深度差异小于二者桶长差，因此L/D较大无量纲上拔荷载反而较小．此外，L/D越大桶底处不排水抗剪强度越大，模拟中越不容易达到屈服，单元网格较晚达到扭曲变形，达到峰值上拔荷载时上拔高度也越大.

图 11 长径比L/D对上拔荷载的影响Fig. 11 Influence of the parameter L/D on the uplift load. The red, green, blue, wathet and purple lines are for aspect ratio of 0. 5, 1. 0, 2. 0, 3. 0 and 5. 0, respectively.

本研究对3种长径比（L/D=0. 5、1. 0和2. 0）的群桶基础上拔进行模拟，群桶桶间距s/D=0. 2、0. 5、1. 0、1. 2和1. 5．桶间距对3种长径比群桶基础上拔荷载的影响可扫描文末右下角二维码查看图S1．其中，n为桶的数量．从图S1可见，群桶基础承载力要低于相同数量单桶基础，随着桶间距减小，上拔荷载更早达到峰值，且群桶基础承载力与桶间距变化正相关．此外，对于不同长径比桶基，当桶间距s/D>1时，群桶上拔荷载与单桶都差异较小，尤其在上拔前期( w_c/D<0.03)上拔荷载出现屈服前，单桶群桶上拔荷载几乎重合，表明当s/D>1时，群桶效应不明显.

定义群桶峰值荷载与相同数量单桶峰值荷载比值为极限承载力群桶效应系数，图 13为不同长径比群桶基础极限承载力群桶效应系数．从图 13可以看出，极限荷载群桶效应系数与桶间距正相关．当桶间距较小时，在群桶效应下桶形基础承载力降幅可达10%；而当桶间距超过1. 5D时，群桶效应可以忽略.

图 12 不同长径比地基土应力分布（a）L/D=0. 5；（b）L/D=2. 0 Fig. 12 Influence of the parameter L/D on the stress distribu‐tion of the clay (a) L/D=0. 5, (b) L/D=2. 0.

图 13 极限承载力群桶效应系数Fig. 13 Caisson group effect coefficients of the ultimate bearing capacity. Color lines are for this study, in which the black, red and green represent aspect ratio of 0. 5, 1. 0 and 2. 0, respec‐tively. Square symbolic lines are for literature 15, in which Solid and hollow represent aspect ratio 0. 25 and 1. 0, respectively. Circular symbolic lines are for literature 16, in which Solid and hollow represent aspect ratio 0. 25 and 1. 0, respectively.

与文献[15-16]研究结果对比可以看出，在桶间距s/D=1. 00时，群桶效应作用下桶形基础承载力均下降5%左右．本研究计算结果中，群桶效应随着桶间距增加而减弱，而文献[15-16]结果中桶间距对不同长径比群桶基础承载力影响变化规律不一致，这可能是由于土体本构选择引起的差异．此外本研究实现了对更小桶间距（s/D=0. 25、0. 50）群桶基础的分析模拟，可更全面反映桶间距对群桶承载力影响.

图 14为长径比为1的群桶基础上拔至0. 1D时不同桶间距桶周土体变形．从图 14可见，群桶效应与桶间距密切相关，当s/D=0. 50时，对称边界处土体最大变形约为0. 2 m，群桶效应明显；而当s/D=1. 5时，对称边界处土体变形很小，在0. 025~0. 050，此时群桶效应几乎不存在.

图 14 不同桶间距桶周土体变形(wc/D =0.1，L/D =1.0) （a）s/D=0. 5；（b）s/D=1. 5 Fig. 14 Soil deformation surrounding the caissons with (a) s/D=0. 5, (b) s/D=1. 5. Here, wc/D=0. 1, L/D=1. 0.