移动射流切削黏土过程的ALE 数值方法

顾磊，倪雁，周凡，倪福生

（1.河海大学疏浚技术教育部工程研究中心，江苏 常州 213022；2.中交上海航道勘察设计研究院有限公司，上海 200120）

0 引言

射流在诸多黏土施工中得以广泛应用。港口航道的疏浚工程中，利用射流不仅可以提高挖掘产量，而且能够减少耙头的堵塞几率[1]。在海底开沟工程中，射流方式也因操作简单和移动部件少而得以采用[2]。

目前对移动射流切削黏土的研究还主要集中在实验领域。Shusen Zhang等实验获得了冲沟的截面形状，理论分析了该形状形成的原因[2]。唐立志针对施工需求通过实验进行施工参数比选，依据实验后的形状比较了不同参数下的射流破土效果[3]。作者也开展了一系列移动射流切削黏土实验研究[4]。但是实验仅能测得切削结果，由于喷嘴移动、作用时间短、流场混浊等因素，移动射流切削过程的信息获取非常困难。这不利于对移动射流切削黏土的机理性分析，导致对施工的理论指导不够，往往依据经验或试运行确定射流参数，必然造成不必要的能量消耗，加大施工成本。

数值模拟方法可以获知射流与黏土作用过程的详细信息。现有模拟射流与土作用的数值方法包括Fluent[5]、Flow-3d[6]、光滑粒子法[7]（Smoothed Particle Hydro-dynamics，简称SPH）和任意拉格朗日-欧拉法（Arbitrary Lagrange-Euler，简称 ALE）。实验可知，移动射流切削黏土过程中，黏土的破坏面变形较大，土粒间黏性是黏土抵抗变形的主要作用。根据对上述数值方法进行比较可知[8]，ALE方法对黏土的射流切削模拟更为适合。本文对移动射流切削黏土过程的ALE数值模拟方法进行探讨。

1 ALE方法理论基础

ALE算法兼具Lagrange方法与Euler方法二者的特长，既可像Lagrange网格一样准确反映边界的变形与运动，又可在求解过程中适当调整网格，使计算区域内部网格扭曲最小化，不致出现严重的畸变，从而解决黏土在射流切削过程中的大变形问题[9]。具体实现时通过ALE网格随变形网格的移动来实现，这种移动产生了相对速度，使得ALE算法的控制方程与Euler法有所不同：

质量守恒方程：

式中：ρ为密度；υi为速度张量；w为相对速度；σij为应力张量；bi为体积应力张量；E为内能；含wi的项即为由于坐标转换引起的额外项，反映了各物理量有Lagrange网格向Euler网格的映射。

射流与黏土的作用通过流固耦合算法模拟，本文选用Euler-Lagrange耦合算法，该方法在建模时可将射流与黏土的模型重叠在一起，再通过约束方法实现两者力学参数的传递。本文采用罚函数表征这一约束。通过判断黏土节点与射流网格间的相对位移d，当d贯穿射流界面网格时，即在两者之间产生界面力F：

式中：ki为刚度系数，反映了侵蚀距离对界面力的影响。这一界面力分别作用于射流和黏土上，即会引起其各自的速度、加速度以及变形或位移的变化，从而实现两者的耦合作用。

2 模型建立过程

考虑到射流的流体物质持续运动，同时在与黏土的作用面上具有大变形，射流区域采用ALE网格算法，而静止的黏土块则采用Lagrange网格算法。黏土块尺寸应考虑消除边界效应的影响。所建具体模型如图1所示。

图1ALE网格示意图Fig.1 ALE grid schematic diagram

2.1 移动射流源

首先必须实现移动射流源的模型，分两步骤实现，一是射流连续喷射，二是射流源移动。

射流连续喷射表现为喷嘴出口源源不断地出现带有垂直速度的流体，这通过设定流体物质速度即可实现。在整个计算期间设定喷射速度函数，并通过关键字*BOUNDARY_PRESCRIBED_MOTION_SET将该曲线与射流流体的速度值关联，即实现了连续喷射的射流源。

射流源移动同样如此，在计算时间周期内设定移动速度函数，并将该速度函数通过关键字*ALE_REFERENCE_SYSTEM_GROUP加载到射流的ALE物质组上。

2.2 黏土模型

软件中自带有土壤的材料模型，本文选用MAT_147材料。该材料模型中基本包含了所有黏土特性相关的基本参数，其中最为关键的参数包括：土粒比重Gs、含水率ω、密度ρ、剪切模量K、体积模量G、摩擦角φ和黏聚力C。为获得准确的模拟结果，依据SL-237—1999《土工测试规程》对上述参数进行了测定，其测定方法和具体数值见表1。

表1 黏土主要参数值及测定方法Table 1 Main parameters and determination methods of clay

2.3 侵蚀算法

在射流切削黏土过程中，一旦黏土被冲破，即被水体带走，切削下的黏土颗粒对切削过程基本再无影响。MAT147材料模型中并不带有材料的侵蚀算法，若直接计算，则黏土在射流冲击下只会发生变形，直至网格发生畸变而计算无法进行为止。通过关键字*MAT_ADD_EROSION可以向黏土模型中引入侵蚀算法，该算法包含网格的失效准则，计算中网格一旦满足设定的失效准则，网格即会消失，从而较好地反映黏土在射流切削过程中的上述行为。

射流冲击在黏土上产生载荷，黏土在这一载荷下发生破坏，这一破坏与土体受压的破坏类似[10]。可通过三轴压缩实验测量其破坏特性。进行了2次无侧限压缩实验，实验结果重复性较好，其基本的主应力-应变曲线如图2所示。

图2 三轴压缩实验及结果Fig.2 Three axis compression experiments and results

图中，随着应变增大黏土主应力迅速增加，在较小应变时即达到最大主应力值，表明此时土样内部已出现破坏，将此时应变称作第一临界应变，记为εc1。图中显示，此时土体仍然具有抵抗变形的能力，其强度并未立即变为0，而是逐渐减小，这实际是破坏面逐渐扩展的过程，直到应变达到某一值后，主应力归于恒定，此后抵抗变形强度主要来源于破坏面的摩擦作用。将主应力达到恒定值时的应变称作第二临界应变，记为εc2。

旨在维持土体长期稳定的土力学研究中，普遍认为受到极限应力土体内部即出现破坏，历经长时间的演变即出现土体失稳，故此类研究一般将极限应力（图2中的A点应力）作为土体失效的临界条件。但黏土在射流切削下的破坏与此不同，黏土颗粒必须从土体中完全脱离才可视为破坏完成。因此，选择εc2作为土体失效的判据条件。

需要指出的是，在三轴压缩实验中，试样的压缩速率较慢，最大仅为4.14 mm/min，而在射流切削过程中其压缩速度极快，这会增强黏土的抗破坏强度[11]，模拟中并未考虑这一效应的影响。

3 模拟结果与分析

在LS-DYNA软件环境下，采用上述模型对移动射流切削黏土过程进行了模拟。模拟的具体参数为：射流喷射速度为33 m/s，移动速度为0.4 m/s，喷嘴直径为2 mm，喷嘴靶距3 cm。

俯视的切削沟槽形状如图3所示，与之对比的是相同工况下的实验结果，可以看到模拟与实验的冲坑在泥面处宽度都极窄，仅比射流直径稍大，这与Shusen Zhang的实验结果[2]一致。

图3 切削沟槽俯视图Fig.3 Overlook map of cutting groove

记沿射流移动方向过射流中心的铅垂面为纵截面，冲坑最深处沿宽度方向的铅垂面为横截面。则切削沟槽的纵截面与横截面形状分别如图4、图5所示。

图4 切削沟槽纵截面Fig.4 Longitudinal section of cutting groove

图5 切削沟槽横截面Fig.5 Cross section of cutting groove

图4中，纵截面形状包括3部分：初始斜面段、底部平整段和末端斜面段。在初始斜面段，切削深度逐渐增大，这是由于射流不仅对其直接冲击的土体产生作用，而且由于流体和土体的压力传递，会将应力传递至周围土体而使其发生变形，当射流冲击到这一变形的土体后，其失效条件更易达到，从而产生更大的深度。但是上述应力传递的范围有限，当射流移动超过该范围后，射流冲击区域的土体均已发生相同变形，深度即不再增加，出现了平整段。在冲刷末端模拟结果与实验稍有不同，模拟的斜面更为陡峭，这可能是由于实验时为了节约土样在末端覆盖钢板所引起的。

从图5可以看到，模拟与实验的切削沟槽侧面都基本垂直，底部相对平整。但与实验结果相比，模拟的沟槽宽度较小，而深度较大。这是由于模拟中并未考虑变形速率对土体强度的增强效用，将导致切削深度加大，而射流反冲水流的强度降低，对两侧的作用减小，故宽度也相应变窄。

综上所述，模拟获得的切削沟槽纵截面和横截面形状均与实验结果基本相同，表明所采用的计算方法能够有效模拟移动射流切削黏土的过程。

4 结语

利用ALE方法对移动射流切削黏土过程的模拟进行了探索，确定了准确模拟移动射流源的方法，确定了黏土模型中关键参数的实验获取手段，结果表明，通过侵蚀算法能够较好地描述黏土在射流冲击下的失效行为，本文方法可以有效模拟移动射流切削黏土过程。通过本方法模拟可以准确获取土体应力和射流流场的详细信息，从而为其机理性研究提供有效的数值手段。