1. 引言
吸力锚是一种被广泛应用于海洋工程的基础形式，具有造价低、施工简单、施工速度快以及可回收等优点，该基础也被称为吸力沉箱(Suction caisson)、吸力桶形基础(Suction bucket)或吸力桩(Suction pile)。在吸力锚基础沉贯安装过程中，首先依靠自身的重力，沉贯至海床一定深度后，内部形成密封的水体，然后将封闭水体排出锚外，使锚的内、外部产生足够大的压差即吸力，在吸力作用下沉贯至设计深度。然而，正是由于这种独特的沉贯安装方式，吸力锚在安装过程中存在不同程度的不确定性。例如：在安装过程中，锚内土体会向上隆起，使吸力锚无法沉贯至设计深度以及安装负压过大导致吸力锚基础周围土体强度失效等问题。吸力锚能够顺利安装并沉贯至指定的设计深度处，是保证其发挥最大抗拔承载作用的保证。
在吸力锚沉贯过程中，由于其独特的沉贯方式，贯入阻力会受到许多因素的影响，国内外学者对此进行了大量的研究[1]。朱兴运[2]通过对长径比为2的吸力锚沉贯过程进行分析，揭示了负压产生的吸力对贯入阻力的影响。Andersen等[3]对大量的现场测试结果和室内外实验结果进行总结分析，提出了吸力锚在砂土中沉贯时的贯入阻力计算公式；李大勇等[4]通过控制水位深度和抽吸速度条件进行沉贯实验对比，实验结果表明，吸力锚在中粗砂中沉贯时负压形成的渗流能够促使贯入阻力减小；国振等[5]在黏土中进行了吸力锚室内模型试验，对贯入阻力的变化趋势进行了分析，当沉贯方式由压力沉贯转变为吸力沉贯后，贯入阻力会随着沉贯的进行先逐渐减小，然后再逐渐增大，并超过压力沉贯阶段的贯入阻力。综上所述，在吸力锚沉贯过程中，贯入阻力受到多种因素的影响，想要准确地确定贯入阻力的大小还很困难，到目前为止，有关于贯入阻力的计算，还缺乏一套完善的理论体系。
与物理模型试验相比，数值模拟经济成本较低，同时可以获得通过试验很难获取的数据，因此许多学者采用数值模拟方法模拟吸力锚沉贯过程。韩智臣[6]使用ALE数值模拟方法建立吸力锚沉贯大变形二维有限元模型，分析土塞高度和贯入阻力的变化趋势，并将该方法计算的结果与实验和相关经验公式结果进行比较，验证了该数值模拟方法的可靠性；Zhou等[7]采用有限元方法分析了吸力锚在沉贯过程中土体流动机制对土塞效应的影响，同时也揭示了吸力锚壁端部形状效应对土塞效应产生的影响；Wang等[8]使用欧拉–拉格朗日耦合方法分析了吸力锚的几何形状和摩擦系数对贯入阻力和土塞隆起高度的影响。综上所述，众多学者开展了吸力锚负压沉贯过程的数值模拟研究，大多数学者只对土塞效应受到的影响因素进行分析，没有揭示土塞效应的演变机制，因而本文研究重点分析了吸力锚沉贯过程中土塞效应的演变机制。
在吸力锚成层土沉贯大变形数值模拟分析中，土体单元会受到吸力锚的挤压产生大变形，进而导致计算结果很难收敛，国内外学者使用了不同数值模拟方法来解决这一问题。王贺等[9]采用离心机模型试验和ALE流固耦合大变形有限元数值计算方法，揭示了吸力锚在砂土中的贯入阻力与其渗流场之间的关系以及各自的变化趋势。Zhou等[10]使用AFENA软件模拟吸力锚在正常固结黏土中分别采用压力式沉贯和吸力式沉贯过程，分析这两种方法的差异性；王胤等[11]采用基于离散单元法与计算流体动力学理论(CFD-DEM)流固耦合方法模拟吸力锚在砂土中的沉贯过程，通过实验验证了该方法的可靠性。Chen等[12]基于流体力学和多孔介质渗流理论，采用有限差分法模拟吸力锚的沉贯过程，研究负压吸力对土体应力变化的影响。由此可见，虽然有多种数值模拟方法可以被用于模拟吸力锚沉贯过程，但还有尚待进一步改进和探讨的方面。
在吸力锚成层土贯入大变形数值模拟分析中，受吸力锚挤压影响，会使土体网格产生较大的变形，进而导致计算结果收敛困难，因此采用有限元方法模拟吸力锚沉贯过程颇具有挑战性。欧拉–拉格朗日耦合方法综合了欧拉方法中材料和网格相互独立以及拉格朗日方法的高效性和准确性等优点，能够很好地解决此类问题。为此，采用欧拉–拉格朗日耦合方法，建立吸力锚成层土贯入大变形数值模型，分析土塞效应演变机制以及贯入阻力的变化趋势。
2. 欧拉–拉格朗日耦合方法
欧拉–拉格朗日耦合(Coupled-Eulerian-Lagrangian, CEL)方法综合了拉格朗日方法的高效性和准确性与欧拉方法中材料和网格相互独立等优点，可同时兼顾吸力锚的小变形和土体的大变形。在沉贯过程中，吸力锚结构为小变形体，采用拉格朗日网格，而土体则产生大变形，采用欧拉网格。该方法利用欧拉方法中背景网格固定，材料可以在背景网格中自由流动的机理模拟土体受到挤压产生大变形过程，同时利用罚接触方法模拟土体和吸力锚之间的接触关系，建立吸力锚–土体相互作用大变形模型，避免了数值模拟收敛困难等问题，因而该方法可以有效地模拟吸力锚沉贯过程的大变形问题。
2.1. 欧拉分析
欧拉分析方法的独特之处在于网格和材料是分开的，正是由于这一特点，才被广泛用于网格大变形数值模拟分析中。在该分析方法中网格结点固定，材料在固定的背景网格中自由流动，因此在每个增量步结束时需要确定欧拉域材料的边界位置。在计算过程中，一般采用EVF (Eulerian volume fraction)确定欧拉材料的运动，EVF表示欧拉网格中填充材料的程度，通过材料占据网格空间的比例确定，当网格中没有材料时EVF = 0，当网格中完全填充材料时EVF = 1。综上所述，建立欧拉网格的区域应该大于材料区域，为材料提供足够的移动和变形空间，因此在建立吸力锚成层土贯入大变形模型时，地基土上部应该预留一定的空隙空间，模拟土体受到挤压产生隆起现象。
2.2. 罚接触方法
罚接触方法是由Olovsson [13]提出的，主要用于模拟拉格朗日域与欧拉域界面之间的接触关系。在时间步长tn的初始时刻，已知上一时间步结束时的界面作用力，通过该作用力确定该时间步长初始时刻拉格朗日域和欧拉域中材料的运动。当该时间步结束时，利用欧拉材料节点和拉格朗日材料节点的相对位移计算节点罚力(界面上对应点之间的接触力)，采用式(1)计算罚力。
(1)
式中，Fp表示界面上对应点之间的接触力，被称为罚力；dp为罚位移，表示节点之间的相对位移；λp为罚刚度。
罚接触方法虽然能准确地确定界面节点上的作用力，但是不能够完全约束欧拉材料侵入拉格朗日区域。如果欧拉网格单元小于拉格朗日单元时，部分欧拉网格单元界面不与拉格朗日节点接触，欧拉网格界面将不会受到拉格朗日节点的反作用力，会导致部分欧拉域材料“渗漏”过拉格朗日的表面，侵入进拉格朗日区域内部。有很多方法能够解决这一问题，其中最简单直接的办法就是缩小拉格朗日网格的尺寸，当拉格朗日网格远小于欧拉网格时就不会出现“渗漏”问题。
2.3. 模型验证
周密等[14]为验证本研究有限元分析结果的准确性，将离心机试验观察到的土体流动结果和测得的贯入阻力与其研究结果相对比。模型验证分析中采用土体软化和率效应模型，选取数值包括参考剪切应变率γref = 0.10 s−1，速率参数µ = 0.10，相对延展性系数ξ95 = 10。结果显示土体的流动模式、地表隆起和回流入加劲肋之间的填充土的数量都有较高一致性。可见数值模拟方法具有极高的可靠性。
3. 数值模型及材料属性
在数值模型中，吸力锚沉贯模型的地基土的长度和宽度均为吸力锚桩结构直径的10倍，高度为18 m，其中地基土模型顶面一下0.4 m为空穴单元，无材料属性，用于体现吸力锚沉贯过程中土塞隆起高度的变化，其余部分为欧拉土体区域，土体区域包含淤泥质黏土、淤泥混砂、粉质黏土以及粉细砂，土层厚度分别为6.5 m、1.5 m、5.5 m及1.5 m，建立的模型结构如图1所示。考虑到该模型结构以及受力具有对称性，仅按1/2模型进行建模。约束吸力锚结构的水平位移和转角，对吸力锚结构整体施加竖直向下的位移，模拟其沉贯过程。在吸力锚内部土体表面施加均布荷载，等效模拟负压作用效果。采用罚接触方法模拟土体和吸力锚之间的接触关系，并通过接触对算法(contact pair)建立面面接触。采用扫略技术进行网格划分，吸力锚结构划分为1015个8节点拉格朗日单元(E3D8R)，土体划分为1,643,700个8节点欧拉单元(EC3D8R)。
Figure 1. Model rendering
图1. 模型效果图
吸力锚结构采用线弹性本构模型，弹性模量E = 2.06 GPa，泊松比v = 0.3。Mohr-Coulomb理想弹塑性本构模型。吸力锚桩结构的直径为3.0 m，桩身长为13 m，壁厚为50 mm。
数值模型中的地基土主要包含四种土质，自上而下依次是杂填土、淤泥质黏土、粉质黏土和粉细砂。各层土的力学特性参数如表1所示。
Table 1. Mechanical parameters of stratified soil
表1. 成层土力学特性参数
土体名称 |
重度γ/kN/m³ |
摩擦角φ/˚ |
弹性模量E/MPa |
粘聚力C/kPa |
泊松比v |
杂填土 |
17.80 |
12.00 |
7.41 |
14.00 |
0.38 |
淤泥质黏土 |
18.00 |
25.00 |
10.00 |
12.00 |
0.34 |
粉质黏土 |
19.00 |
32.00 |
12.00 |
17.00 |
0.35 |
粉细砂 |
18.00 |
35.00 |
48.00 |
8.00 |
0.31 |
4. 数值模拟结果
4.1. 土塞隆起高度
吸力锚沉贯过程中，土体受到挤压，锚内会产生土体向上隆起的现象。通过CEL方法建立上述工况中的数值模型，分析吸力锚在该工况下的土塞高度隆起的过程，土塞高度随沉贯深度变化如图2所示。根据土塞工期高度变化趋势，依然可将土塞效应划分为三个阶段，分别为隆起增长期、隆起过渡期以及隆起平缓期。吸力锚沉贯结束时，土塞隆起高度为0.56 m。根据土塞隆起高度随时间变化关系，在沉贯初期时，土体隆起高度随沉贯深度变化呈现出线性增长变化，该阶段被称为隆起增长期。当吸力锚继续沉贯时，土塞隆起高度随沉贯深度缓慢增长，该阶段被称为隆起过渡期。在沉贯后期时，土塞隆起高度随沉贯深度缓慢增长并趋于稳定时，该阶段被称为隆起平缓期。
Figure 2. Curve of soil uplift height variation with penetration depth
图2. 土体隆起高度随沉贯深度变化曲线
4.2. 土体应力
吸力锚沉贯过程中会对土体产生挤压，影响土体水平应力的变化。通过CEL方法模拟实际工况中吸力锚的沉贯过程，分析土体应力的变化。如图3和图4所示分别为沉贯深度h为2 m、6.2 m、6.8 m、7.7 m、8.3 m以及12 m时的土体水平应力云图和竖向应力云图。
如图3和图4所示，水平应力随沉贯深度变化呈现出先逐渐增大，之后缓慢减小的变化趋势，竖向应力随沉贯深度逐渐增大，并形成竖向应力拱，且竖向应力拱的分布范围逐渐减小。在沉贯过程中土压力逐渐增大，土体的水平应力和竖向应力逐渐增加，之后由于负压的作用，会对土体的水平应力产生影响，而对竖向应力影响较小。因此在沉贯后期，水平应力随沉贯深度逐渐减小。
(a) h = 2 m (b) h = 6.2 m
(c) h = 6.8 m (d) h = 7.7 m
(e) h = 8.3 m (f) h = 12 m
Figure 3. Horizontal stress contour (unit: Pa)
图3. 水平应力云图(单位：Pa)
(a) h = 2 m (b) h = 6.2 m
(c) h = 6.8 m (d) h = 7.7 m
(e) h = 8.3 m (f) h = 12 m
Figure 4. Vertical stress contours (unit: Pa)
图4. 竖向应力云图(单位：Pa)
4.3. 沉贯阻力
如图5所示为通过数值模拟方法计算得到的内、外壁及端部阻力随沉贯深度变化趋势。内、外壁及端部阻力随沉贯深度变化逐渐增大，锚外壁受到土体阻碍的作用最大，其次是锚内壁。同时，在沉贯过程中土体强度对沉贯阻力产生影响，当土体力学特性较强时，沉贯阻力会增大，内、外壁以及端部阻力皆受到土体强度的影响，且在同一层土中沉贯时变化趋势相同。
Figure 5. Variation trends of internal and external wall resistance and tip resistance with penetration depth
图5. 内外壁以及端部阻力随沉贯深度变化趋势
5. 结果与展望
本文采用欧拉–拉格朗日耦合方法，建立吸力锚成层土贯入大变形数值模型，以模拟吸力锚在自重和负压作用下的沉贯过程，进而得出以下主要结论：
1) 根据土塞隆起高度随贯入深度变化关系，对沉贯过程中吸力锚内部土体速度矢量场进行深入分析，将土塞效应划分为3个阶段：隆起增长期、隆起过渡期，隆起下滑期；
2) 在吸力锚沉贯过程中，水平应力形成的挤土效应对土塞效应的影响程度逐渐减小，竖向应力在锚端部附近的土体中形成向下凸起的应力拱，限制土体的隆起，土体中产生的应力应变主要集中在锚端部附近，随着贯入深度的增加逐渐增大；
3) 土塞效应能够使侧摩阻力增大，尤其对内壁侧摩阻力的影响程度最大，在沉贯初期，土塞效应不对外壁侧摩阻力产生影响；
4) 随着贯入深度的增加，端部阻力先逐渐增大，而后逐渐减小，土塞效应能够使端部阻力有所减小。