网箱养殖是我国海水鱼类养殖的主要方式，发展深水网箱养殖是我国海洋开发的重要战略和设施养殖的主要趋势^[1-2]。深水网箱鱼类养殖过程主要采用自动投饵系统进行投饵，在自动投饵系统中，采用管道气力输送颗粒饲料。颗粒饲料在输送过程中，常常会由于输送工艺参数选取不合理导致饲料颗粒堆积阻塞管道或饲料破损等问题^[3-7]。

目前对气固两相流的管道气力输送的研究主要针对粮食、橡胶等颗粒在直管弯管的运动特性模拟^[8-10]，祝先胜^[11]模拟了三维90°弯管结构内的气固两相流动过程，获得了不同工况下的颗粒运动轨迹、速度、碰撞作用力、碰撞频率及流体压力、速度分布等。王明旭等^[12]以小麦颗粒为研究对象，建立了以气流速度、小麦颗粒在输送过程中的体积百分含量为优化参数，流场内部压力为目标的小麦颗粒输送特性优化模型。卢洲等^[13]采用计算流体力学和离散单元法耦合模拟弯管内的柱状颗粒气力输送过程，对不同弯径比的90°弯管内柱状颗粒的运动状态、碰撞特性、破碎原因及相关的力学特性进行了研究。针对气力投饵管道输送的研究主要是通过现场实测比较输送工艺参数对投饵的影响。Oehme等^[14]研究了深水网箱鲑鱼养殖气力投饵系统中，不同气力输送速度和不同的颗粒粒径，颗粒饲料旋转撒料装置的抛洒面积分布状况。Halstensen等^[15]提出了一种鱼饲料颗粒速度的无损监测方法，给出了颗粒速度和空气流量的部分最小二乘验证结果。Turidsynnøve等^[16]利用3种不同物理品质三文鱼颗粒饲料，在不同风速及进料速率下，研究气力输送时的颗粒破碎情况，发现在高空速和低进料速率下，颗粒破碎率最高。国内外众多学者就气固两相流建立起了许多数学模型，并将这些模型应用于输送对象广泛的颗粒气力输送中，对气力输送下的颗粒运动特性进行了一些研究，而这些气力输送的数值模拟主要以二维气固两相流体模型为主，这与实际气力输送过程尚存在不少偏差^[17-22]。

本文基于气固两相流理论，利用计算流体力学和离散单元法进行耦合模拟，对气力输送下颗粒饲料在投饵系统管道初始阶段到稳定阶段运动过程中的流动状态与运动特性进行研究，并分析不同阶段饲料颗粒动能、管内的压降、物料所受的曳力。最后根据仿真分析结果，对设备结构改进提出了指导意见。

1.   数学模型

1.1   气固两相流基本方程

气体的流动规律用气固两相耦合的N-S方程描述。其中气相连续性方程：

气相动量方程：

固相连续性方程：

固相动量方程：

$ {\alpha _g}$为气相体积分数，$ {\alpha _p}$为固相体积分数，$ {\alpha _g}$+ $ {\alpha _p}$=1，

$ {u_{gi}}$、$ {u_{gj}}$、$ {u_{pi}}$、$ {u_{pj}}$分别为气相、固相在$ i{\text{、}}j$方向上的速度分量；$ {\rho _g}$、$ {\rho _p}$分别是气相、固相的密度，$ {g_i}$为重力加速度在$ i$方向上的分量；$ {\mu _e}$、$ {\mu _p}$为气相、固相在有效粘度；$ {\beta _i}$为气固两相交换系数在$ i$方向上的分量。

1.2   颗粒接触碰撞模型

采用Hertz接触理论，根据Mindlin切向接触力求解方法，建立Hertz-Mindlin无滑动接触模型。其中：

其中$ {R^*}$为等效粒子半径，$ {E^*}$为等效弹性模量，$ {m^*}$为等效质量，$ {G^*}$为等效剪切模量，$ \delta $为切向重叠量，$ v_n^{rel}$为相对速度的法向分量值，$ v_t^{rel}$为切向相对速度。

2.   投饵系统结构

自动投饵系统装备主要由罗茨风机、旋转供料器、料仓、水平管道、S形管道与输料管道等构成(图1)。罗茨风机产生的空气流进入管道，料仓的颗粒饲料通过旋转供料器排入输送管道，空气流混合颗粒饲料后，进入S形输料管，由出料管输送到输料管道。

图  1  自动投饵系统装备原理结构图

Fig. 1  Principle and structure diagram of automatic feeding system

下载: 全尺寸图片 幻灯片

3.   模型建立

3.1   离散元模型

实际饲料颗粒为膨化颗粒，形状可近似为高6 mm、底面直径7 mm的圆柱体，其尺寸见图2-a。为了更好地使离散元模型与真实饲料形状接近，在EDEM中利用particle工具创建颗粒模型，采用18个直径为3.5 mm的球形颗粒聚合成饲料离散元模型(图2-b、c)。物料属性及接触参数见表1。

图  2  多球面饲料离散元模型

Fig. 2  Discrete element model of multi sphere feed

下载: 全尺寸图片 幻灯片

表  1  EDEM中物料属性及接触参数

Table  1  Material properties and contact parameters in EDEM

名称 name	参数 parameter	数值 volume
物料 material	泊松比	0.4
	剪切模量/Pa	1×108
	密度/kg·m–3	390
壳体 shell	泊松比	0.3
	剪切模量/Pa	7.0×1010
	密度/kg·m–3	7 800
物料-物料 materiel-materiel	恢复系数	0.2
	静摩擦系数	0.6
	动摩擦系数	0.01
物料-壳体 materiel-shell	恢复系数	0.15
	静摩擦系数	0.3
	动摩擦系数	0.01

下载: 导出CSV 

| 显示表格

为了提高仿真计算速度，节约计算机资源，对自动投饵颗粒饲料气力输送设备结构进行简化，只保留料仓、供料器、管道等部件。模型长2 919 mm，宽800 mm，高1 574 mm，简化后见图3。

图  3  自动投饵设备结构示意图

Fig. 3  Schematic diagram of auto-feeding equipment

下载: 全尺寸图片 幻灯片

3.2   计算流体模型

仿真过程中，忽略供料器与输送管道衔接部位的流场区域，只对输送管道进行网格划分。为了准确增加流场内网格节点和捕捉数据变化，采用Ansys Workbench中meshing对管道内部区域进行非结构四面体网格划分，共产生116 675个网格，节点数42 206个。导入fluent后经检查无负体积网格产生。

供料器转速为23 r·min^–1，EDEM中仿真时间步长设置为5e^–6 s，为瑞利时间步的35.3%，共进行4 s。仿真开始前利用颗粒静态生成法在料仓内产生了4万个多球面饲料颗粒。

在fluent中选择k-$ \varepsilon $湍流模型，入口边界速度为28 m·s^–1，湍流强度为5%，湍流黏度静速率10。速度压力耦合求解器选择SIMPLE，采用二阶迎风离散格式。

耦合算法选择欧拉-欧拉双向模型，fluent时间步长为5e^–4 s，为EDEM的100倍，共8 000个时间步长。

4.   仿真结果分析

饲料颗粒从供料腔进入输送管道后，由于主要受到气流的曳力作用开始做加速运动。饲料颗粒在气流入口段速度较小，基本维持在10 m·s^–1以下，随着运动距离的不断加大，颗粒也获得较高速度(图4)。图中管道末端饲料颗粒速度达到最大(26 m·s^–1)，并且出现悬浮现象，这与实际现象一致。该设备由于供料器是圆柱形壳体，落料口并没有直接与供料器壳体两端面连接，仿真过程中发现供料器靠近两端面处有颗粒滞留，导致无法充分供料。在后续结构改进时，要对供料器壳体做进一步优化，使其可以对颗粒运动有良好的导向性。

图  4  饲料管道内输送过程

Fig. 4  Transportation process of pellet feed pipeline

下载: 全尺寸图片 幻灯片

供料器接口处至管道末端位置水平距离为2 780 mm，对此段长度按695 mm等分成4段，对管道内颗粒平均动能进行统计。颗粒在气流推动下获得的动能沿运动方向不断增加(图5)。第1段区域内颗粒处于刚开始运动阶段，平均动能较小，基本为0.002 5 J；第2段区域颗粒速度维持在13 m·s^–1左右，处于加速阶段，平均动能为0.02 J，相比第一阶段获得运动较多；第3区域中管道有距离较短的弧形弯道，颗粒经过此段后速度继续增加，平均动能为0.025 J，同时该阶段较前2个阶段平均动能波动性增加；第4区域中，颗粒运动速度达到最大(26 m·s^–1)，平均动能接近0.03 J，颗粒运动趋于稳定状态。对比4个阶段可以发现，颗粒在第1、第2阶段处于初始运动状态，气流要对饲料颗粒做功，使其发生运动，所以这2个阶段较第3、第4阶段耗能更大，气流压降也更大。颗粒经过加速后进入第3、第4阶段，速度逐渐达到最大值，并维持稳定状态，但由于颗粒空隙率增加，动能波动性加强(图6)。

图  5  颗粒饲料动力统计区域

Fig. 5  Dynamic statistics area of pellet feed

下载: 全尺寸图片 幻灯片

图  6  不同阶段饲料颗粒动能

Fig. 6  Kinetic energy of feed particles at different stages

下载: 全尺寸图片 幻灯片

饲料颗粒具有不同的速度时在管道内分布的形态也不一样。颗粒从供料器进入输送管道后，由于喂入量较大，颗粒加速缓慢，颗粒与颗粒之间空隙率低，颗粒成整体移动方式前进(图7-a)。在颗粒整体前进一段距离后，由于存在颗粒间碰撞作用力、颗粒与壁面碰撞作用力、场内漩涡等较复杂情况，导致颗粒间空隙率增加，单个颗粒速度增加明显，团聚移动方式减弱(图7-b)。在弯管中，饲料颗粒由于离心力作用紧贴外侧管道壁面前进，在弯管末端颗粒被抛向水平管道上表面，此时与壁面产生冲击力作用，颗粒回弹到水平管道下表面，此过程中可以看到颗粒速度有所下降，产生了能量损失(图7-c)。对管道结构改进时，要减小弯管曲率，以降低弯管与直管连接处颗粒碰撞的能量损失。在图7-d中，颗粒运动速度达到最大，同时出现悬浮状态，颗粒位置布满整个管道截面，颗粒间空隙率也达到最大。整个过程中并未见到颗粒堆积现象发生，管道内排料顺畅。

图  7  饲料颗粒不同运动段形态分布

Fig. 7  Morphological distribution of feed particles at different motion stages

下载: 全尺寸图片 幻灯片

4.1   管内压降

饲料颗粒沿着管道前进所需的能量是由气流压降来实现的。在仿真过程中，根据能量守恒定律，一定质量的流体在入口截面所具有的能量等于颗粒运动的能量、沿程压力损失、局部损失、出口截面能量之和。仿真采用简化的管道模型，长度较短，并且没有不同段的辅助衔接装置。因此，忽略气体流动过程的局部压力损失。将饲料颗粒从供料器下落至管道内的位置作为起点，运动方向为正向，每隔695 mm做一管道总压云图(图8)。随着饲料颗粒运动距离增加，管道内总压下降明显。饲料颗粒在管道内刚开始发生运动时，气流总压为4 600 Pa。颗粒发生水平距离2 780 mm时速度达到26 m·s^–1，管道总压3 614 Pa，压降986 Pa。这个过程饲料颗粒运动速度从0增加到26 m·s^–1，压降能量大部分用于转化为颗粒运动的动能及颗粒群提升的重力势能。

图  8  颗粒运动不同距离处管道截面总压云

Fig. 8  Total pressure cloud of pipe cross section at different distance of particle movement

下载: 全尺寸图片 幻灯片

4.2   物料所受曳力

物料所受耦合力有曳力、马格努斯旋转提升力、赛夫曼剪切提升力、压力梯度力等，在仿真中后3项作用力微小，因此计算过程中将其忽略。饲料颗粒所受曳力公式为：

式中$ {C_D}$为阻力系数；$ {\mu _p}$为球形颗粒在静止流体中的运动速度(m·s^–1)；$ {d_p}$为球形颗粒直径(m)；$ \overline \rho $为流体密度(kg·m^–3)。

将非球形饲料颗粒等效为球形，$ {d_p}$=7 mm，Re=2.02×10⁶对应的$ {C_D}$近似为0.2，$ {\mu _p}$=60 m·s^–1，代入公式(12)可得$ {F_D}$=1.79×10^–2 N，则单个饲料颗粒运动1 m气流作功为1.79×10^–2 J。

4.3   空载时管道压降

在此工况下，空气为不可压牛顿流，因此，由雷诺数计算公式$ {\mathop{\rm Re}\nolimits} = \displaystyle\frac{{\rho vd}}{\mu }$可得Re=2.02×10⁶。所以，管道内空载时压降计算公式为：

式中$ \Delta p$为直管压降(Mpa)；$ {f_{sw}}$为流体摩阻系数；$ {\rho _w}$为流体密度 (kg·m^–3)；$ {v_w}$为流体流经管道的速度(m·s^–1)；$ {L_y}$为直管长度(m)；D为直管直径(m)。

此工况下Re > 2 100为湍流运动，摩阻系数公式为：

式中$ {e_p}$为直管的表面粗糙度(m)。

由以上公式代入数据计算可得1 m管道内$ \Delta p$=24 Pa。

5.   结论

1) 利用CFD-DEM耦合法对饲料颗粒在气力输送设备管道中初始阶段的运动过程进行了分析，得到颗粒从初始状态到运动稳定阶段的颗粒位置分布情况，以及不同阶段颗粒的平均动能大小。整个初始阶段过程中并未见到颗粒堆积现象发生。

2) 饲料颗粒运动速度从0 m·s^–1增加到26 m·s^–1时，颗粒饲料管道入口到颗粒运动稳定状态时刻的压降为986 Pa，压降能量大部分用于转化为颗粒运动的动能及颗粒群提升的重力势能。

3) 根据仿真分析结果，对设备结构改进提出了意见。对供料器壳体做进一步优化，使其可以对颗粒运动有良好的导向性。减小弯管曲率，以降低弯管与直管连接处颗粒碰撞能量损失。

4) 经现场测试，投饵机水平管道出口风速为24.03 m·s^–1，与仿真模拟结果比较接近，说明采用CFD-EDM两相流耦合技术进行饲料颗粒气力输送运动特性的模拟仿真是可行的。

5) 后续研究中，将着重讨论饲料尺寸、管道尺寸、管道弯曲程度、风速、风压等工艺参数对颗粒饲料运动特性的影响规律，获得管道尺寸、输送参数和饵料特性的最优化方案。