基于CFD-DEM的排肥用波纹管结构优化设计与试验

张一帆，何瑞银，段庆飞，徐 勇

(南京农业大学 工学院，江苏 南京210031)

肥料是作物的“粮食”，合理使用肥料可以改良土壤肥力、提高单位面积产量[1-3]。然而随着肥料使用量的增大与不科学施肥方式的增多，国内施肥行业产生了肥料用量大、利用率低、面源污染严重等问题[4-6]。传统粗放的施肥方式已不适合新时代绿色农业的发展要求，急需新型施肥技术提高肥料利用率，减少使用量[7-9]。气力输肥系统由于排肥均匀性好，输送效率高近年来在施肥领域得到了广泛应用[10-11]。波纹管具有良好的理化特性，许多专家学者在气力输肥系统中使用波纹管作为关键输送部件并开展了相关研究。杨庆璐等[12]研究了波纹管直径对气固两相流动性的影响，通过颗粒碰撞次数与颗粒分布图判断均匀性。张晓辉等[13]以气流场内速度分布均匀性作为评价指标，通过对比直管与褶皱管内部速度场分布判断输送性能。戴亿政等[14]通过空气流的单相仿真分析了波纹管结构对内部气体速度场的影响。李中华等[15]研究了波纹管内部单相压力场，得出了波纹管结构导致气流压力平缓降低。常金丽[16]设计了一种褶皱管，用来增压提升输送均匀性。以上学者大都只针对波纹管某一结构参数变量进行定性研究，鲜有研究不同结构参数变量之间交互作用对输送均匀性的影响。在判断输送均匀性时大多依据非标准化的流场分布图，没有参数化的指标进行均匀性判断。

本文旨在通过CFD-DEM耦合仿真分析技术，在单因素试验的基础上研究波纹管幅宽、波纹间距和波长之间交互作用对输送均匀性的影响，通过仿真试验明确波纹管结构参数与输送均匀性的关系。构建以排肥均匀性变异系数为响应值的二次响应面模型，并利用响应面分析法获得波纹管的最佳参数组合以提高气力排肥均匀性。最后通过台架试验对仿真模型与最优参数组合进行验证，确定仿真模型的可靠性。

1 基于CFD-DEM的耦合仿真模型建立

1.1 几何模型建立

参考戴亿政等[14]的研究，波纹管几何结构如图1-a所示，图中d表示波纹管内径，A表示幅宽，L表示波长，S表示波纹间距，H表示总长。三维建模效果如图1-b所示。

a，波纹管几何结构图；b，三维建模效果图。

1.2 Fluent参数设置

1.2.1 数学模型

本次研究模拟工作介质为不可压缩牛顿流体空气，忽略流体流动时黏性耗散作用，假设波纹管壁面为刚性[17]。连续性方程、动量方程、能量方程分别表示如下[18]：

连续性方程：

∇u=0。

(1)

动量方程：

(2)

能量方程：

(3)

式中：u为速度矢量，m·s-1；t为时间，s；T为温度，K；ρ为密度，kg·m-3；f为单位质量上的质量力，N·kg-1；p为压强，Pa；μ为动力黏度，Pa·s；cp为定压比热容，J·kg-1·K-1；κ为导热系数，W·m-1·K-1。

1.2.2 网格划分与模拟方法

利用前处理软件Meshing进行网格划分，由于波纹管本身曲率变化较大，对整体进行网格划分时采用Tet/Hybrid非结构性网格[19]。划分完成后通过软件自带的Statistic功能进行质量检验[20]。划分后网格数量233 634，节点数量80 911，平均偏度(skewness)为0.25，正交质量(orthogonal quality)为0.823，说明网格质量较好可以用来模拟，网格划分效果如图2所示。考虑到实际使用情况，将进出口适当延长以消除仿真过程中物理边界对进出口边界影响。

图2 网格划分结果示意图

模拟方法根据流动状态雷诺数判断，雷诺数计算公式如下：

(4)

式中：Re表示雷诺数；ρ表示密度,kg·m-3；ν表示速度,m·s-1；L表示特征长度,m；μ表示动力黏度,Pa·s-1。

试验模拟工作介质为标准条件下空气，密度为ρ=1.29 kg·m-3，动力黏度μ=1.79×10-5Pa·s，特征长度L=0.06 m。参考文献[21]，模拟风速ν=25 m·s-1，将空气物性参数与风速带入公式(4)可得波纹管入口处雷诺数Re=108 136>4 000，因此本次模拟需要使用湍流模型。湍流模型采用标准K-epsilon模型，算法选择Simple二阶迎风差分格式离散。入口边界条件选择速度入口，出口边界条件选择压力出口，各项变量的收敛残差均为10-6。

1.3 EDEM参数设置

参考文献[22]，仿真试验所用肥料颗粒为尿素，泊松比0.51，弹性模量8.9×107Pa，剪切模量3.56×107Pa，密度1 337 kg·m-3，粒径3.77 mm 。在EDEM中接触模型均选择Hertz-Mindlin (No Slip)碰撞接触模型，颗粒与颗粒之间碰撞恢复系数为0.28，静摩系数0.36，滚动摩擦系数0.15，颗粒与设备之间碰撞恢复系数为0.36，静摩系数0.41，滚动摩擦系数0.04[23]。参考高观保[24]、杨庆璐等[12]的研究，初步设定本研究中起始施肥速率350 g·s-1。

1.4 耦合模型参数设置

采用基于Lagrangian法的多相流模型的耦合算法，设置动量亚松弛、体积亚松弛因子为0.3[25]。由于气流速度相对于模型尺度较大，为捕捉到足够数据点，设置Fluent时间步长0.000 1 s，每100步记录1次，Edem时间步长5×10-6s，记录间隔0.01 s。

2 单因素仿真试验

2.1 排肥均匀性指标建立

利用虚拟网格划分技术，在波纹管靠近出口处设置网格组，该网格组以管道中心为基准均分为6个扇形网格，在仿真时长内对通过各个网格的颗粒数量进行累加计数，并由EDEM计算得出通过该网格颗粒的总质量，网格划分与计数效果如图3所示。

图3 划分网格箱体并对各个箱体内颗粒累加计数

再通过公式(5)、公式(6)、公式(7)计算排肥均匀性变异系数：

(5)

(6)

(7)

式中，aj为各网格累计颗粒质量平均值，单位g；ax为第x个网格内累计颗粒质量，单位g；Sj为各网格累计颗粒质量标准差；Vj为排肥均匀性变异系数；nj为测定网格个数。

2.2 单因素试验设计与分析

参考杨庆璐等[12]、戴亿政等[14]的研究，模拟波纹管内径固定为60 mm，总长500 mm，以波长、波纹间距、幅宽3个最能反映波纹构造的因素作为试验因素，不同因素水平选取采用等距排列，每组试验重复3次取平均值，因素水平表如表1所示。

表1 单因素试验因素水平表

2.2.1 幅宽对排肥均匀性的影响

如图4所示，排肥均匀性变异系数随幅宽增加先减小后增加，在幅宽A<6 mm时变化幅度较大，幅宽A>6 mm后变化幅度较小，变异系数总体变化幅度>44%，推测幅宽处于6～10 mm，有最小值点。综上，选取幅宽A为6、8、10 mm作为响应面试验的3个水平。

图4 不同幅宽下排肥均匀性变异系数变化

2.2.2 波长对排肥均匀性的影响

如图5所示，排肥均匀性变异系数随波长增加先减小后增加，波长L<10 mm时变化幅度较大，波长L>10 mm后变化幅度较小，变异系数总体变化幅度>25%，推测波长处于10～20 mm，有最小值点。综上，选取波长L分别为10、15、20 mm作为响应面试验的3个水平。

图5 不同波长下排肥均匀性变异系数变化

2.2.3 波纹间距对排肥均匀性的影响

如图6所示，排肥均匀性变异系数随波纹间距增加先平缓减少再缓慢增加，变异系数总体变化幅度<5%。综上，选取波纹间距S分别为0、10、20 mm作为响应面试验的3个水平。

图6 不同波纹间距下排肥均匀性变异系数变化

3 基于BDD法的响应面优化试验

3.1 响应面模型建立

为进一步研究不同波纹管结构参数之间交互关系对输送均匀性能的影响，进行BDD法响应面分析。选取幅宽、波纹间距、波长作为影响因素，排肥均匀性变异系数作为响应值。参考上一节分析结论，实验组因素编码及水平如表2所示。利用Design Expert 10软件计算仿真数据，得到响应面试验结果如表3所示。

表2 BDD法实验组因素编码及水平

表3 实验组BDD法响应面试验设计与结果

3.2 响应面结果与模型分析

3.2.1 残差分析

BDD法残差分析结果如图7。图7-a中数据点基本处于一条直线上，呈现线性分布，且置信度在95%以内，说明残差符合正态分布。图7-b中所有残差都围绕0随机分布，说明拟合效果较好。图7-c中实际值与预测值的点基本贴近于拟合曲线，说明两者拟合程度较高。

a，内学生化残差正态分布概率；b，外学生化残差与预测值；c，预测值与实际值。

3.2.2 方差分析与显著性检验

利用Design Expert软件对试验数据进行拟合分析，响应面的回归方程为：

Y=8.20-1.53A+0.93B-1.48C-1.47AB-1.7AC-1.83BC+2.08A2+1.80B2+3.39C2

(8)

表4 方差分析结果

3.2.3 基于3D响应面的最优解分析

为了进一步考察不同因素之间的交互关系以及对均匀性的影响，并寻求最优参数组合。根据试验结果建立3D响应曲面如图8所示。

由图8-a可知，当幅宽一定时，变异系数随着波长与波纹间距的增大均显示出先减小后增大的趋势。最优解范围对应的波长14～20 mm，波纹间距5～15 mm。由图8-b可知，当波长一定时，变异系数随着波纹间距与幅宽的增大均显示出先减小后增大的趋势。最优解范围对应的波纹间距5～15 mm，幅宽7～9 mm。由图8-c可知，当波纹间距一定时，变异系数随着幅宽与波长的增大均显示出先减小后增大的趋势。最优解范围对应的幅宽7～9 mm，波长14～20 mm。可以看出，3个响应曲面都呈现非线性关系，并且曲面向下弯曲，说明有最小值解。

a，波长与波纹间距对变异系数的交互作用；b，波纹间距与幅宽对变异系数的交互作用； c.波长与幅宽对变异系数的交互作用。

为了获得最优参数组合，根据响应面分析初步获得的参数范围，利用Design Expert软件进行分析。最终得到拟合的结构参数最优解，波长18.103 mm，波纹间距12.158 mm，幅宽8.863 mm，预测的变异系数7.5%。

4 台架验证试验

4.1 试验装置与材料

参考文献[26]制作试验台架，台架如图9所示主要包括风机、风机调速器、排肥器、排肥控制盒、供料喷射器、波纹管、分配器以及其余管路部件。其中排肥控制盒可以调节排肥器转速以控制肥料排量，风机调速器可以调节风机转速以控制风速。试验材料选择尿素颗粒与史丹利复合肥。波纹管结构参数参考上文2.2小节表1制作，最优解波纹管结构参数实际加工精度为波长18.08 mm，波纹间距12.19 mm，幅宽8.84 mm。

1，排肥管；2，分配器；3，波纹管；4，弯管；5，供料喷射器；6，排肥控制盒；7，风机调速器；8，排肥电机；9，风机；10，排肥器；11，肥箱。

4.2 试验方法

试验方法参考文献[27]与GB/T 20346.1—2006《施肥机械试验方法》，试验条件为施肥速率350 g·s-1，入口风速25 m·s-1。台架试验主要的评价指标为排肥均匀性变异系数，具体计算公式如下。

(9)

(10)

(11)

式中，ah为各行排肥量平均值，单位g；ax为第x行排肥量，单位g；Sh为各行排肥量标准差；Vh为排肥均匀性变异系数；nh为测定行数。

4.3 试验设计与结果分析

4.3.1 波纹管结构参数单因素验证试验

为验证仿真模型精确性对尿素颗粒仿真值与实际值进行对比试验分析，单因素试验用波纹管结构参数参考2.2节表1制作。尿素颗粒台架试验结果与仿真试验结果共同绘制如图10所示，由图10可知，仿真预测值与台架实际值趋势基本一致，数值基本重合，说明仿真模型精确可靠。

a， 不同幅宽仿真值与实际值对照；b，不同波长仿真值与实际值对照；c，不同波纹间距仿真值与实际值对照。

为验证仿真模型普适性，对尿素颗粒与史丹利复合肥颗粒实际值进行对比试验分析。史丹利复合肥颗粒与尿素颗粒台架试验结果如图11所示。由图11可知，各因素水平下复合肥颗粒变化趋势与尿素颗粒基本一致，说明仿真模型基本具有一定普适性。

a，不同幅宽尿素与复合肥对照；b，不同波长尿素与复合肥对照；c，不同波纹间距尿素与复合肥对照。

4.3.2 最优参数优化结果验证试验

为验证仿真试验优化结果，针对最优参数组合波纹管进行台架验证试验，共重复5次。最优解结构参数实际加工精度为波长18.08 mm，波纹间距12.19 mm，幅宽8.84 mm。

试验结果如表5所示。由表5结果可知，最优参数组合波纹管排肥均匀性变异系数台架试验结果与仿真试验预测值相对误差的均值为5.84%，台架试验结果与仿真试验结果基本一致，说明仿真模型的优化结果可靠，可以利用仿真模型进行优化设计。

表5 最优参数台架试验结果

5 结论

(1)建立了一种基于CFD-DEM 气固耦合仿真模型，构建了以排肥均匀性变异系数为响应值的二次响应面模型，并进行了单因素试验与响应面试验。仿真试验结果表明，排肥均匀性变异系数随波长、幅宽、波纹间距的增加先减小后增加，幅宽变化造成的影响最为显著；波纹管结构参数对排肥均匀性的影响显著性顺序为幅宽>波长>波纹间距。当波长18.103 mm、波纹间距12.158 mm、幅宽8.863 mm时输送均匀性最优，仿真试验预测的排肥均匀性变异系数为7.5%。

(2)为检验仿真模型精确性与普适性，搭建气力排肥系统试验台对尿素颗粒与复合肥颗粒分别进行台架试验。单因素台架试验中尿素颗粒仿真预测值与台架实际值趋势基本一致，数值基本重合，说明仿真模型精确可靠。复合肥颗粒变化趋势与尿素颗粒基本一致，说明仿真模型基本具有一定普适性。

(3)为验证仿真试验优化结果，针对最优参数组合波纹管进行台架验证试验。最优参数组合波纹管排肥均匀性变异系数台架试验结果与仿真试验预测值相对误差的均值为5.84%，台架试验结果与仿真试验结果基本一致，说明基于仿真模型的优化可靠，本研究可为波纹管的使用与优化提供参考。