基于EDEM-Fluent耦合的气爆松土效果仿真分析

陈金楚，奚小波，殷慧子，张翼夫，张宝峰，金亦富，张瑞宏

（扬州大学机械工程学院，江苏 扬州 225127）

果园土壤合理有效深松有利于疏通土壤之间的间隙，降低水土流失，利于复合施肥作业，为肥料扩散提供空间，提高农作物产量[1-3]。传统果园松土施肥作业是机器或人工开沟松土、施肥、覆土[4]，易出现伤害作物根系、破坏土壤结构等情况[5]。气爆松土技术不破坏土层结构，通过输送气体对给定深度的土壤进行爆气疏松，并在土层中形成空隙，利于后续液肥扩散[6]。

近年来，基于土体气压劈裂原理，孔德刚等提出气压式土壤深松方法[7]；增子利一等通过向果树根部注入压缩空气实现快速松土和深层施肥[8]；张瑞宏等设计适用于液肥注施的气动深松施肥机[9-10]；奚小波等优化气爆松土参数[11-12]。目前气爆松土一般通过试验分析指标参数对土壤的大体扰动效果，对土壤颗粒的具体扰动效果缺乏理论研究。土壤颗粒运动的理论研究大多采用离散元方法[13]（Discrete element method,DEM），仅考虑土壤颗粒间产生的接触力、摩擦力，忽略实际土壤间粘聚力和塑性变形[14-15]。Walton等针对塑性材料创建迟滞弹性接触模型（Hysteretic spring contact model,HSCM）[16]，而后Janda等验证其准确性[17-18]。Ucgul等将迟滞弹性接触模型（HSCM）和法向粘聚模型（LCM）结合表征现实土壤颗粒之间的塑性变形和粘聚力[19-20]，使法向力随颗粒间迭代量呈理想线性关系：当外部施加力小于设定应力值时，执行弹性变形；当外力超过设定应力值时，发生塑性变形。考虑到土壤挤压变形属于非线性形变，Thankur等创建爱丁堡弹塑性附着力模型（Edinburgh elastoplastic adhesion model，EEPA）[21-22]，学者们根据土壤性质将土壤结构分为耕作层、犁底层和心土层，不断优化试验，验证EEPA模型合理性。

本文采用EDEM-Fluent耦合进行气爆松土仿真，研究不同爆气流孔结构及气爆参数对不同深度土壤扰动效果影响，拟合建立土壤扰动轨迹线方程，通过试验验证仿真结果，为气爆松土工艺参数制定及机具结构设计提供理论依据。

1 气爆松土总体结构参数设计

如图1所示，气爆松土工作部件主要由进气接头、气缸和钻杆3部分构成。气铲标准验收气压为机械行业标准JB/T 8412-2016[23]规定0.63 MPa，气铲尾部钻杆结构重量为2.43 kg。根据标准JB/T 8412-2016中气铲各类型技术规格，考虑到钻杆作业行程，选择冲击行程有优势、冲击频率低、运行平稳的C7型气铲。具体参数如下：机重7.4 kg，系统总重9.83 kg，机长445 mm，缸径28 mm，冲击频率13 Hz，冲击能17 J，耗气量16 L·s-1，冲击行程100.04 mm，气管内径13 mm，工作气压0.63 MPa。

图1 气爆松土机整体结构Fig.1 Overall structure of air burst soil loosening machine

考虑钻杆松土时顺畅程度和结构刚度，尽可能使钻杆外径趋于最小值，以减少作业时在轴向上产生的合阻力，优化后钻杆壁厚5 mm、内径15 mm、锥角2a取60°。在流孔部分，孔心距离钻杆平面l取20 mm，钻杆长度L为950 mm[12]。流体分布器孔口结构是影响液体分布器分布均匀性的因素之一[31]，本文在保持流孔体积不变情况下，根据截面形状变化设计3种爆气流孔结构杆，其流孔数量均为4，如图2所示。由图2可知，3种爆气流孔结构杆分别为：①外扩型，外口直径大于内口直径，且成60°凹角，横截面Φ1值为7 mm；②外缩型，外口直径大于内口直径，且成60°凹角，横截面Φ2值为7 mm；③直口型，外口直径等于内口直径，横截面Φ3值为9.9 mm。

图2 不同流孔的钻杆结构参数Fig.2 Drill pipe structural parameters of different flow holes

2 离散元土壤模型建立

2.1 土壤颗粒模型

土壤颗粒组成通常为单颗粒核状、双颗粒条状、多颗粒（≥3）团簇状等形式，如图3所示。

图3 土壤颗粒形状Fig.3 Soil particle shape

在EDEM软件中，基础球颗粒数目越多，仿真时间越长，为减少颗粒数量，保证土壤体积，选择相对较大粒径的颗粒，考虑到果园土壤常年不翻耕，土层结构紧实，故采用半径3.5 mm、形式单一团簇颗粒形状，模拟粗粒组中粒径较小的圆砾和角砾。犁底层颗粒间滚动摩擦系数为0.21，静摩擦系数为0.35[24]。颗粒仿真参数：密度2 600 kg·m-3，弹性模量1×106Pa，泊松比0.3，恢复系数0.6[25]。

2.2 土壤接触模型

为更好体现土壤力-形变关系，本文采取土壤模型为EEPA，区别于传统弹性接触模型，考虑土壤颗粒间产生非线性粘聚和塑性变形[26-27]，可有效描述黏附颗粒应力变化过程。颗粒间作用力为法向力，形变为压力卸载后颗粒间保留重叠量。土壤之间力关系如图4所示。

图4 力关系原理[28-29]Fig.4 Schematic of force relations[28-29]

当颗粒间法向力增加时，颗粒间重叠量δ相应增加。接触按照加载（Loading）曲线斜率加载；当施加力消失，重叠量相应降低，接触按照卸载/重新加载（Unloading/reloading）曲线的斜率卸载。颗粒在分离前继续施加力，按照k2加载。当法向力小于0时，且趋于峰值f max时，颗粒间粘结力随之消失。此后，法向力按照-kabh斜率恢复到颗粒间初始接触力f0。EEPA模型相对于HSCM+LCM模型结合优势体现在图4负半轴上。EEPA在EDEM软件中设置参数如下：黏附力强度0 N，接触黏附能50 J·m-2，接触塑性比0.7，非线性曲线幂指数1.5，黏附分支功率5，切向刚度系数0.28571[30]。

2.3 土壤和钻杆仿真模型

因气爆松土作业对象为果园林木土壤，果园田间管理机械和收获机械作业过程土壤压实，且常年不翻耕，造成土质板结，因此本文仿真可将其近似于土壤层结构中的犁底层，在EDEM中生成的土壤模型设置厚度为500 mm。通过土壤颗粒模型搭建形成虚拟土体，尺寸为1 200 mm（长）×800 mm（宽）×500 mm（高）。犁底层土壤团聚现象明显，所以尺寸分布为固定值1，颗粒工厂产生700 000个团簇体颗粒。生成的土壤模型如图5所示。

图5 土壤模型Fig.5 Soil model

采用Solidworks三维软件创建对应钻杆几何模型，材料选用304不锈钢，材料密度ρ为7850kg·m-3，剪切模量G为1×1010Pa，泊松比ν为0.25，钻杆和土壤之间恢复系数K为0.6。同犁底层之间静摩擦因素为0.639；滚动摩擦因素为0.13[16]。为更好满足软件间耦合时间步长匹配需要的整数倍关系，将时间步长设为0.00005 s，总体时间设为1 s，使EDEM中土壤颗粒生成后充分堆积接触，单元尺寸（Cell size）设置为3Rmin生成最终土壤颗粒模型。3类流孔钻杆与不同深度土壤层进行接触。

3 耦合仿真过程与结果分析

3.1 前处理过程

为简化模型结构，节约仿真时间，本文仅对钻杆部分进行相关处理分析。将钻杆三维模型导入Workbench软件中。由图2空腔作为钻杆内部流体域。上端面设为inlet，下端面气孔4个面设为outlet，钻杆圆柱内壁设为wall。自动划分网格，确保fluent流体域网格尺寸大于颗粒体积。

3.2 耦合计算过程

在Fluent软件中，选择压力出入口，根据实际工况设置出入口压力。初始排气压力值分别为0.4、0.6、0.8 MPa，容积流量为16.667 L·s-1。不计工作过程中磨损等各因素。根据空压机压降公式：

式中，Δp为压力降（Pa）；q为容积流量（L·s-1）；L为管子长度（m）；d为管道内径（mm）；p为排气压力（Pa）。最终计算得出口压力值分别为397 436.46、

598 290.97、798 718.23 Pa。

因气体与颗粒间作用仅发生在钻杆流体域管道外部，只需对Fluent产生气体与EDEM生成颗粒之间进行单向计算。耦合是基于不考虑体积分数影响的多相流耦合接口。在耦合时，接口以串行方式打开（并行时求解器进程设为0）。将动量交换源项的数学模型编译进edem_udf耦合文件，通过此文件加载连接EDEM和Fluent求解器。

湍流模型选择Standardk-ε模型，依据本文局部固相体积保持比例占总体积比例，本文选用Lagrangian法耦合。跨相参数交换中的阻力模型、热传递模型均采用默认设置，升力模型选择Saff⁃man Lift和Magnus Lift体现速度梯度以及颗粒自旋所产生的升力。

3.3 不同流孔结构的松土效果分析

通过Fluent设置的气压值、与EDEM生成土壤层不同接触深度及3种流孔结构钻杆在气爆压力0.4 MPa、松土深度100 mm条件下仿真。时间步长为0.0005 s，EDEM端的10倍，步数为2000，耦合总体时间为两者乘积值1 s，此时Residuals值开始收敛，适于本文分析。

图6a、b、c是3种钻杆出口气流速度云图，可见外扩型最大气流速度约为40 m·s-1，外缩型最大气流速度约为50 m·s-1，直口型最大气流速度约为60 m·s-1。由于重力作用，流孔出口速度云图分布呈向下分布趋势，外扩型出口速度分布效果最差，外缩型最好，直口型次之。结合速度值和分布均匀性情况，最终确定直孔型流孔对土壤颗粒扰动效果最佳。通过对3组数据生成的土壤颗粒从Y方向远离土壤模型中心70 mm、厚度150 mm的Slice处理，得出图6d、e、f 3种钻杆对土壤颗粒的扰动效果。本文通过选取土壤颗粒运动速度＞0.045 m·s-1部分作为土壤扰动活跃区域，通过EDEM中的速度调节及颜色划分并去除边界条件观察，可见，随出口气流速度增大，土壤颗粒扰动分布区域越大，颗粒轮廓线大体为抛物线，且抛物线顶点与气爆深度接近。从土壤扰动颗粒最外缘范围看，外扩型钻杆为（-370 mm，400 mm），外缩型为（-375 mm，400 mm），直口型为（-385 mm，400 mm），可进一步确定直口型为最佳钻杆流孔结构，其径向扰动区域范围达到785 mm。

图6 3种钻杆出口气流速度和土壤颗粒扰动效果Fig.6 Three kinds of drill pipe outlet airflow velocity and soil particle disturbance effect

3.4 不同气爆参数松土效果分析

上述仿真结果分析表明，直口型是符合作业能效最佳结构形式。为进一步研究气爆压力P与松土深度H对土壤扰动效果影响，对9组P-H条件下土壤扰动仿真分析，结果如图7所示。可看出，当H一致时，P越大，土壤颗粒扰动区域也越大；当P一致时，H不超过200 mm时，土壤颗粒扰动区域随H增大而增大，当H为300 mm时，随P变化，土壤颗粒扰动区域分布差异不明显，这是因为深度越大土壤本身重量越大，P在一定范围内对其扰动效果不明显。

图7 不同气爆压力P与松土深度H条件下的土壤扰动分布Fig.7 Distribution of soil disturbance under different air burst pressure Pand loose soil depth H

利用Solidworks软件进行截图采点，确定生成轨迹中代表性颗粒点坐标，再用Origin软件点状图分析拟合多项式方程，如图8所示。

图8 土壤颗粒速度分布的轨迹线Fig.8 Trajectory line of soil particle velocity distribution

拟合出的轨迹线方程为h=ar2+br+c型，由于理想条件下轨迹应为与h轴对称的抛物线，方程中b为0，同时为方便看出松土深度H最大扰动深度，可将方程简化为h=ar2-H+c´型，最终得出9组土壤颗粒扰动区域轨迹线方程：

从各组拟合曲线对应的a值和c´值来看，当作业深度为100mm时，土壤扰动深度与作业深度一致；随松土深度增大，颗粒扰动最大深度下移，大于松土作业深度。

图9为土壤表面土体扰动全域分布情况。其扰动范围类似椭圆，长轴Δx在785~900 mm、短轴Δy在560~620 mm范围内变化。当P一致时，H越大，Δx值变大，Δy值变小；当H一致时，P越大，产生同样变化趋势。在土壤表面扰动面积方面，通过椭圆面积公式计算得出，0.8 MPa-100 mm组合扰动面积最大，0.6 MPa-100 mm组合次之，0.4 MPa-100 mm组合最小。考虑到松土深度和实际作业时最佳效果，进一步比较300 mm深度时扰动范围面积大小，得出0.6 MPa-300 mm为最合适参数组合。土壤表面土体扰动范围值见表1。

表1 土壤表面的土体扰动范围值Table 1 Disturbance range value of the soil surface

图9 土壤表面的土体扰动全域分布Fig.9 Global distribution map of soil disturbanceon soil surface

9组仿真数据的抬升高度Δh与不同气爆压力-松土深度组合之间的关系见图10，可见，土壤模型整体抬升高度与深度呈反比，与气压值呈正比，随气压值增大，其整体抬高变化量趋于平缓。

图10 不同气爆压力P与松土深度H条件下土壤模型抬升高度Fig.10 Lifting height of thesoil model under the condition of different air burst pressure Pand loose soil depth H

4 试验验证

4.1 试验目的与条件

为进一步验证气爆松土实际效果并判断仿真结果准确性，制作标准试验土体，其土块质地为壤土，紧实度为3.28 MPa，含水量为16.6%，室外温度为18℃。试验装置如图11所示，其中空压机型号为罗翔W-1.0/8，满载气压为0.8 MPa，排气量1 m3·min-1，功率7.5 kW，气筒容量230L，转速2 900 r·min-1，气压可在0.4~0.8 MPa范围内调节。试验设备还有TJS-100型土壤紧实度测定仪（测量深度0~450 mm，精度±0.1%），TZS-I型土壤水分测定仪（含水率0~100%，相对误差≤3%）、钢板尺（量程50 mm，精度1 mm）、卷尺（量程5 m，精度1 mm）等各类用于测量标记的工具。

图11 气爆松土试验装置实物Fig.11 Picture of gas explosion subsoiling experimental equipment

4.2 试验方法

气爆松土试验时，气体在土壤层由内向外扩散，抬起土体，最终疏松的土体实现整体脱离，为更好观察气爆后土体扰动情况，处理土体剖面：距气爆中心点50 mm处标记，用铲子进行人工剖面。本试验参数选用具有代表性的气爆压力0.4 MPa、松土深度100 mm和气爆压力0.8 MPa、松土深度300 mm，测量土壤在横纵方向上的扰动范围，对比分析产生的裂隙采点标记坐标及仿真结果。

4.3 试验结果与分析

气爆后的土体内部产生裂隙，通过图12试验结果与仿真结果比较发现，0.4 MPa-100 mm组合的Δx试验值为730 mm，Δy试验值为540 mm，与仿真值误差均小于7.6%；0.8 MPa-300 mm组合的Δx试验值为820 mm，Δy试验值为515 mm，与仿真值误差均小于9.8%。

图12 气爆松土试验土体截面形貌Fig.12 Crosssection morphology of soil in air burst test

图13为气爆松土试验裂隙轨迹线的拟合曲线，可以发现，该拟合曲线亦是抛物线，与仿真结论一致，同时0.4 MPa-100 mm组合的轨迹线c´为-5.57 mm，对应仿真值1.678 mm；0.8 MPa-300 mm组合轨迹线的c´为-136.21 mm，对应仿真值-125.664 mm，二者相对误差较小，说明仿真结果合理性与准确性。试验与仿真结果存在误差，因为试验土体是根据果园土壤的紧实度通过人工堆积制作而成，与仿真模型存在一定偏差。

图13 气爆松土试验裂隙轨迹线的拟合曲线Fig.13 Fitting curveof crack track linein air burst soil loosening test

5 结论

a.本文采用EDEM-Fluent耦合作气爆松土仿真，研究不同气爆流孔结构及气爆参数对土壤扰动效果的影响，建立土壤扰动轨迹线方程，并通过试验验证仿真结果。

b.直孔型流孔对土壤颗粒扰动效果最佳，在气爆压力P为0.4 MPa、松土作业深度H为100 mm时最大气流速度达60 m·s-1，径向扰动范围达785 mm。

c.不同气爆参数对土壤扰动效果的仿真结果表明，当H一致时，土壤扰动区域随P增大而增大；当P一致时，H＜200 mm时，土壤扰动区域随着H增大而增大，而当H＞200 mm时，土壤扰动区域随H变化不明显。同时，H＜100 mm时，土壤扰动最大深度与H几乎一致；H＞100 mm时，土壤扰动最大深度大于H。土体抬升高度与H呈反比，与P呈正比，随着P持续增大，抬升高度趋于平缓。

d.气爆松土试验结果显示，气爆能使土体产生裂隙并扩散，其拟合曲线为抛物线，与仿真结论一致，试验土体扰动范围与仿真值的误差较小，仿真结果合理可靠。