岩土体贯入问题属几何非性线大变形问题，传统的数值模拟方法如拉格朗日（Lagrangian）法、欧拉（Eulerian）法已不适应，于是出现任意拉格朗日-欧拉（arbitrary Lagrangian-Eulerian，ALE）法与耦合欧拉-拉格朗日（coupled Eulerian-Lagrange，CEL）法［2］。关于ALE法和CEL法相较于传统方法的优越性引起研究界关注。何健等［3］研究CEL算法在土壤大变形仿真中的应用时发现，相较于Lagrangian法，CEL算法在计算桩打入土壤的过程中不会出现中断的情况。孙肖菲等［4］在研究管土相互作用时对比Lagrange、ALE、CEL方法的网格变形形态，发现Lagrange法模型网格畸变最大，且得到的荷载位移拟合曲线与实测值偏差过大，计算精度不如ALE、CEL法。陈榕等［5］用ALE法对扩底桩上拔承载力展开了研究，验证了ALE法的计算结果可靠性；王腾等［6］采用CEL法模拟黏土中海底管线竖向大变形贯入过程，研究其土体承载力变化规律。
本文从海上工装平台吸力桶贯入过程分析入手，基于ABAQUS软件，分别采用ALE和CEL方法建模分析，对比2种方法的建模工作量、计算精度、计算效率以及网格尺寸、贯入速率等因素的影响敏感性。
1 算法原理
1.1 ALE法
ALE描述中引入了一个独立于初始构形和现实构形的参考构形，基于任意参考构型坐标的随体导数方程［7］为：
[?fXi,t?t=?fxi,t?t+wi?fxi,t?xi] （1）
式（1）中：[Xi]是拉格朗日参考坐标；[xi]是欧拉参考坐标；[wi]是相对速度，[wi=v-u]（[v]是质点的物质速度，[u]是质点的网格速度）。
“桶体”贯入“土体”导致其变形后，土体上的网格也会变形甚至畸变，利用软件的网格重划分技术对变形“土体”进行重新划分，形成新的形状规则的网格，并由已知的参考构形求解出初始构形和现时构形，同时将初始网格计算的解答及状态变量如应力应变等传输到现时网格上，然后进行后续的增量步计算。重复以上过程，既能以高质量的规则网格完成计算，又能保证求解数据完整准确，从而解决大变形问题。
1.2 CEL法
CEL法中的质量、动量和能量欧拉形式的守恒方程可以写成通用的守恒形式［8］：
[???t+??Φ=S] （2）
式（2）中：[?]为场变量，[Φ]和[S]分别为通量函数和源项。
利用算子分裂方法求解，把式（2）分解为含有源项的拉格朗日步和含有对流项的欧拉步［9］：
[???t=S] （3）
[???t+??Φ=0] （4）
在拉格朗日步内，节点和材料暂时绑定在一起，变形一致，从而获得拉格朗日域的运动状况。在欧拉步内，材料变形停止，节点和材料解除绑定，变形的网格重新映射到初始空间位置，利用相邻网格的材料体积变化反映材料的流动，并将场内变量传输到材料上，通过求解方程（4）得到新的应力，并将在下一个时间步内继续循环拉格朗日步和欧拉步［2］，从而解决大变形问题。
2 计算模型
桶径5 m，桶厚0.03 m，桶长7 m。为消除边界效应，土体模型水平尺寸为5倍桶径（25 m），竖向尺寸为20 m（约3倍桶长）［10］。采用简化双层土以模拟实际工程复杂土体，土体弹性模量和抗剪强度随土体深度而变化［11］，选用Drucker-Prager模型作为本构模型，其它参数见表1。
表1 模型物理参数
Tab. 1 Physical parameters of model
[材料 高度 /
m 密度 /
（kg/m3） 泊松比 摩擦角 /
（°） 黏聚力 /
kPa 桶体 7 7 850 0.28 — — 黏土 4 1 630 0.40 10 28 砂土 16 1 800 0.33 30 2 ]
采用ALE法计算时，为了使计算结果更易收敛，建模需要预设贯入路径［12］、预设初始埋深［13］、修改贯入体端部形状［14］。CEL法中网格不会变形，土体流向顶部设置的高度为2 m的“空域”（只有网格没有材料的Void层）。2种土体模型的边界条件相同，侧面边界条件设置为水平向固定，竖直向自由；底部边界水平与竖直向都固定，上部设置为自由边界。2种方法所构建的初始模型如图1所示。
<G:\武汉工程大学\2024\第1期\陈 诚-1.tif>[Void层][ALE][CEL]
图1 ALE与CEL模型
Fig. 1 ALE and CEL models
3 相关技术问题
3.1 地应力平衡
软件通过荷载以外力形式施加给土体重力，一旦开始计算，土体中地应力产生的同时，土层会被压缩产生初始沉降位移。地应力平衡就是要使初始土体模型不存在初始位移，只存在初始应力［15］，与自然界土体性质吻合。
ALE采用输出数据库（output database，ODB）文件导入法，即通过静态通用分析步计算出土体的初始地应力和初始位移，形成ODB文件，再利用地应力平衡分析步消去初始位移。CEL法只能创建动力显式分析步，所以只有采用预定义场法，根据土壤性质计算出重力带来的土体应力，直接利用软件中的预定义场定义土体地应力场［16］。地应力平衡效果以土体具有初始地应力后土体各点竖向位移趋近于0的程度衡量。
图2为土体竖向位移图。如图2所示，由于土层性质复杂，侧向土压力系数只能通过经验确定，所以“预定义场”法中同一土层的竖向“归0”位移不一致，且平衡效果较“ODB导入法”差，但2种方法平衡后位移的数量级基本小于10-4 m，地应力平衡都达到了平衡标准。
3.2 贯入参数确定
3.2.1 ALE参数确定 若ALE中网格畸变太大，会出现虽有单元变形但没有应变能产生的情况（即沙漏现象），严重时会导致模型失效。根据ABAQUS软件帮助手册，沙漏现象严重程度用伪应变能来衡量，当伪应变能超过内能的2%后，仿真无效。ALE中伪应变能越小，模型精度越高。提高网格重划分频率和强度可以降低模型伪应变能［17］，频率可以直接选为1，强度不能设置成无限大，需要根据计算机能力和计算时间反复试算确定。
3.2.2 CEL参数确定 CEL中增量步要小于最小稳定增量步才能保证计算稳定，最小稳定增量步一般都是极小的值，所以时间消耗极大［18］。CEL计算时一般需要通过增大速度与质量因子节省时间，但这样可能会使准静态平衡状态卷入动态，求解难度大大增加。从能量角度考虑，模型的能量平衡方程是：
[E=EI+EV+EKE+EFD+EW] （5）
式（5）中：E是系统中的总能量，EI是内能，Ev是黏性耗散吸收的能量，EKE是动能，EFD是摩擦耗散吸收的能量，EW是外力所做的功。如果模拟是准静态的，那么EW几乎等于E，EV、EFD和EI一般都较小，推断出准静态过程中EKE很小，正确的准静态响应模型动能不会超过内能的5%。通过软件试算得到的模型动能，如表2所示。模型总内能为4.4×106 J，可根据动能占模型内能的比例，选择不同速度下最节约时间的质量放大因子。
表2 模型动能表
Tab. 2 Model kinetic energy table
[质量放大因子 动能 / J 0.20 m/s 0.10 m/s 0.02 m/s 100 1 726 840 431 560 4 315 64 1 105 177 276 198 2 761 49 846 151 211 464 2 114 36 621 662 155 361 1 553 25 431 710 107 890 1 078 16 276 294 69 049 690 9 155 415 38 840 388 4 69 073 17 262 172 ]
4 模拟结果与分析
4.1 网格尺寸对贯入阻力的影响
网格划分影响计算精度、时间、结果离散性。计算结果对于网格尺寸敏感性越低，越容易通过调整网格大小提高计算效率。
贯入体自身的厚度与其周围网格尺寸的关系对于贯入影响最显著，将土体模型分为2类区域（如图3所示，阴影区域为近桶区域，空白区域为远桶区域），近桶区域网格尺寸L与桶体厚度d的比值设为n，ALE中不同速率下n值对贯入深度的影响如图4所示。
<G:\武汉工程大学\2024\第1期\陈 诚-3.tif>[R 25 m][R 2.6 m][R 2.4 m]
图3 土体模型俯视图
Fig. 3 Top view diagram of soil model
由图4可知，n≥2.50时，桶体在不同速率下都可以贯入到预定深度7 m，n≤2.25时，同种网格尺寸下，贯入速率越小，贯入的深度越浅。这是由于网格尺寸相对过小或者贯入速率过小时，桶体贯入过程中会放大网格变形过程中的奇异性，从而导致计算无法收敛。要使桶体贯入到预定深度，需要在不同的贯入速率下选择合适的网格尺寸，不同的网格尺寸下选择合适的贯入速率。
<G:\武汉工程大学\2024\第1期\陈 诚-4.tif>[0 1 2 3 4
n][8
7
6
5
4
3
2
1][贯入深度 / m][0.01 m/s
0.02 m/s
0.05 m/s
0.07 m/s
0.10 m/s
0.20 m/s]
图4 n值对贯入深度的影响
Fig. 4 Effect of n-value on penetration depth
网格划分时，远桶区域网格尺寸/近桶区域网格尺寸可根据图4选择如下4种组合尺寸：0.6m/0.03m，0.8m/0.04m，1.0 m/0.05m和1.2m/0.06m，贯入速率选0.10 m/s，得到了如图5所示的不同网格尺寸下的贯入阻力-深度曲线。
<G:\武汉工程大学\2024\第1期\陈 诚-5-1.tif><G:\武汉工程大学\2024\第1期\陈 诚-5-2.tif>[8
7
6
5
4
3
2
1][贯入阻力 / MN][0 1 2 3 4 5 6 7 8
贯入深度 / m][1.2m/0.06m
1.0m/0.05m
0.8m/0.04m
0.6m/0.03m][（b）][（a）][0 1 2 3 4 5 6 7 8
贯入深度 / m][8
7
6
5
4
3
2
1][贯入阻力 / MN][1.2m/0.06m
1.0m/0.05m
0.8m/0.04m
0.6m/0.03m]
图5 不同网格尺寸下的贯入阻力-深度曲线：
（a）ALE，（b）CEL
Fig. 5 Penetration resistance-depth curves for different mesh sizes：（a）ALE，（b）CEL
由图5可以看出，2种方法在贯入深度小于4 m（0.8倍桶径）时，网格尺寸对计算数值大小和曲线离散性影响不大，这是由于初始贯入时2种方法都可以把模型的变形和接触问题处理好。当贯入深度大于0.8倍桶径时，随着网格更加精细，ALE法计算得到的阻力数值整体变大，但网格精细化到一定程度，结果基本趋近定值。使用ALE法时原始网格的精细化可以促使网格重划分更有效地应对畸变严重的土体区域（图6）。而CEL中网格不发生变形，初始网格大小对计算数值大小整体影响不大，得到的数值大小基本趋近。
<G:\武汉工程大学\2024\第1期\陈 诚-6.tif>
图6 ALE网格变形
Fig. 6 ALE mesh deformation
网格尺寸对于ALE法得到的曲线离散性影响不显著。但CEL法中网格尺寸越大，计算得到的曲线离散性越大。这是由于桶体对土体挤压产生的摩阻力受到界面处摩擦定律的限制，阻力达到极限后会突然下降，然后下一个土壤节点再次被压缩，阻力重新上升，导致了数值分析结果的发散。ALE中网格划分频率和强度足够高时，土壤节点不断变化，曲线离散性更小。
所以在利用ALE法计算时，要在计算不中断的前提下尽量使网格细化。CEL中网格尺寸过大时计算结果离散性较大，当需要某点的精确结果时，尽量细化网格计算。
4.2 贯入速率对贯入阻力的影响
桶体施工贯入过程中的速率选择不尽相同，实际工程中贯入速率会影响桶体贯入稳定性。针对ALE和CEL方法，分别取0.02、0.10、0.20 m/s的贯入速率，在远桶区域和近桶区域网格尺寸组合为1.0m/0.05m的模型中进行研究，得到了贯入阻力-深度曲线，如图7所示。
由图7可以看出，在ALE法中，吸力桶贯入土体中的贯入速率对计算结果与离散性影响不大。但是在CEL法中，通过拟合曲线发现在0.02、0.10、0.20 m/s的贯入速率下得到的决定系数R2分别等于0.975、0.950和0.945，速度越大，R2值越小，说明曲线离散程度越大，而且贯入速率越大得到的阻力值越大。这是由于贯入速率导致的惯性因素容易在CEL法中占据主导，从而影响计算结果稳定性。使用ALE法模拟吸力桶贯入时较稳定，可以适当提高贯入速率，以节约计算时间。
<G:\武汉工程大学\2024\第1期\陈 诚-7-2.tif><G:\武汉工程大学\2024\第1期\陈 诚-7-1.tif>[8
7
6
5
4
3
2
1][贯入阻力 / MN][0 1 2 3 4 5 6 7 8
贯入深度 / m][0.02 m/s
0.10 m/s
0.20 m/s][（b）][（a）][0 1 2 3 4 5 6 7 8
贯入深度 / m][8
7
6
5
4
3
2
1][贯入阻力 / MN][0.02 m/s
0.10 m/s
0.20 m/s]
图7 不同贯入速率下的贯入阻力-深度曲线：
（a）ALE，（b）CEL
Fig. 7 Penetration resistance-depth curves under different
penetration rates：（a）ALE，（b）CEL
4.3 与理论结果对比
将模拟结果与Houlsby等［19］的研究得到的吸力筒沉贯阻力理论值进行比较，网格尺寸选择最精细的0.6m/0.03m组合，贯入速率选择0.10 m/s，得到模拟结果与理论结果的比较，如图8所示。
<G:\武汉工程大学\2024\第1期\陈 诚-8.tif>[0 1 2 3 4 5 6 7 8
贯入深度 / m][8
7
6
5
4
3
2
1][贯入阻力 / MN][CEL法
ALE法
理论值]
图8 模拟与理论比较
Fig. 8 Comparison of values of simulation and theory
由于本模型为双层土，前4 m深度与后3 m深度的理论值曲线曲率不同，2种方法的模拟结果与理论值的曲线趋势接近。4 m深度以后，ALE得到的结果小于CEL得到的模拟结果与理论值，这是由于ALE法克服网格畸变能力有限，当桶体贯入一定深度时，模型中一部分土体的畸变无法通过网格重划分很好地解决，这部分土体与吸力桶壁脱离导致计算阻力值偏小。
5 结 论
为研究数值模拟在贯入问题中的应用，针对吸力桶贯入土体问题，建立ALE法与CEL法有限元计算模型，从实施过程、参数影响、计算效率与结果精度等方面，对比2种方法得到以下结论：
（1）ALE中使用的“ODB导入法”较之CEL使用的“预定义场法”地应力平衡效果更好。
（2）当吸力桶贯入深度小于0.8倍桶径时，网格尺寸对2种方法计算结果的影响不大。当贯入深度大于0.8倍桶径时，随着网格更加精细，ALE中得到的计算结果变大但逐渐趋近，CEL中计算结果无较大偏差但曲线离散性变小。使用2种方法计算时都应尽量使网格精细化。
（3）相较ALE法，采用CEL法计算得到的贯入阻力-深度曲线对贯入速率的变化更敏感，且阻力值随贯入速率的增大而增大。所以使用ALE法模拟时可以适当提高贯入速率以减少计算时间，对计算结果影响不大。
（4）选择合适的参数进行模拟，土体模型变形不大时，2种方法得到的阻力数值与理论值基本一致，但在ALE法中，当贯入体贯入到更深土体时网格畸变加大且土体向外挤压，造成部分土体与贯入体脱离，计算数值小于CEL法数值和理论值。