肖芬芳， 张从合， 王慧， 叶亚峰， 张道林， 汪和廷，李波， 吴跃进， 刘斌美*
（1.中国科学院合肥物质科学研究院，合肥 230031；
2.中国科学技术大学，合肥 230036；
3.安徽荃银高科种业股份有限公司，合肥 230088）
授粉是杂交水稻制种的重要组成部分，我国主要采用人工授粉进行杂交，但近年来随着劳动力成本的大幅度增加，杂交水稻种子生产成本也大幅提升，严重影响杂交水稻生产效益，不利于我国种业持续健康发展[1-3]。不少研究者提出将父、母本分田种植，通过取粉和授粉机械分别完成父本田采粉、母本田授粉，其间通过低温贮存等技术维持花粉活力，可实现真正意义上的机械化制种，是我国杂交水稻规模化、效益化发展的必由之路[4-5]。
早在20世纪90年代，国内学者就采用采粉-贮存-喷粉模式对杂交水稻机械授粉装置进行研究。黄崇德等[6]和舒焜[7]发明了机械化采授粉技术，运用机械采粉、授粉以及相应的花粉贮藏技术完成授粉作业，将相关技术进行推广示范，获得了较好的的增产效果。但由于收集效率不高，收集后缺少合适的储存条件导致花粉活性较低，限制了该模式的进一步研究推广。近年来，随着花粉贮藏技术及雾化喷洒技术的发展，离体花粉可以在一定时间内保持活力及较高的萌发率，并通过特定技术均匀喷洒至母本完成授粉，为花粉贮藏和异地授粉提供了解决方案[8-11]。而在采粉方面，关于花粉高效收集装置研究较少，存在收集效率低等问题，本团队已研制了花粉机械收集装置[12]，目前需对装置进行优化，以提高收集效率，为采粉-贮藏-授粉模式中采粉环节提供花粉高效收集装置。
随着计算机仿真技术的发展，计算流体动力学（computational fluid dynamics，CFD）被广泛应用于谷物清选、壳仁分离、作物收集及气力输送等农业工程领域，用于指导机械设计及结构参数优化[13-16]。在气固两相流方面，为了更真实地模拟双流体运动状态、优化工作参数，当固体体积分数较小可在CFD中引入离散相（discrete phase model，DPM）模型建立连续-离散模型（computational fluid dynamics-discrete phase model，CFD-DPM）对双流体进行描述[17]。彭丽等[18]采用CFD-DPM方法结合响应面模型，对旋风分离器结构参数进行了仿真优化。敬志臣等[19]运用CFD-DPM对3种旋风分离器模型的小麦清选过程进行了仿真研究，根据清选效果确定最优模型并通过试验验证了优化结果的可靠性。本研究采用CFD-DPM对花粉收集装置关键部件气力输送系统进行仿真分析，在此基础上结合响应面法优化输送系统结构参数，以期解决杂交水稻父、母本分田种植下父本花粉机械化高效收集的问题。
1.1 整机结构与工作原理
气吸式花粉收集装置主要由气力输送系统、负压风机组件、过滤分离室、传送带组件、冷却干燥室、冷却气罐组件、传动系统、驾驶室及机架等部分组成，如图1所示。
图1 气吸式杂交水稻花粉收集装置结构Fig. 1 Structure of air-suction pollen collection device
作业时，高速运转的离心风机在装置内部产生较强的负压，将花粉吸入。在气流的运输下，花粉流经气力输送系统进入过滤分离室，经过分离、过滤，落至有一定速度的传送带，并在传送带的运送下进入冷却仓实现即时冷却干燥。
1.2 风机选型
负压风机使气力输送系统内产生大于花粉悬浮速度的空气速度将花粉吸落，同时完成花粉的运输。风机功率主要由风量和风压确定[16,12,20-21]。
1.2.1 风量确定 花粉能被顺利吸起时，吸粉口处风速应大于花粉悬浮速度，负压风机的风量需满足式（1）条件[21]。
式中，Q为风机风量，m3·h-1；
Smax为吸粉口端面面积，m2；
va为花粉的最大悬浮速度，m·s-1，取值3.5。
1.2.2 风压确定 风机的总压Pz主要包括动压损失Pd和静压损失Pj（摩擦压力Pm和局部压力Pj），计算公式如公式（2）。
式中，Pz为总压，Pa；
ρ为空气密度，取1.22 kg·m-2；
v为气流平均速度，取20 m·s-1；
γ为摩擦因数，取0.1；
R0为管路半径，为0.1 m；
ε1为局部阻力系数，取1.3；
ε2为90弯管压损系数，取1.3；
F为风管长度，m；
δ为漏风系数（实际值与理论值之比，约1.2～1.5），取1.4[20,22]。
1.3 气力输送系统仿真建模
气力输送系统在实际作业中完成作物花粉的收集及运输，其结构参数直接影响气流分布、颗粒运动，从而影响花粉收集效率及花粉活力，因此对气力输送系统进行仿真优化具有重要意义。
1.3.1 物理建模 简化后的气力输送系统结构如图2所示，由吸粉口和进粉管组成。主要结构参数包括进粉管内径（diameter of the tube，D）、进粉管水平长度（horizontal length of tube，A）、进粉管竖直高度（vertical height of tube，B）、吸粉口长度（length of the inlet，L）、吸粉口宽度（width of the inlet，W）、吸粉口高度（height of the inlet，H）、扩展域高度（height of the extended domain，H0）、弯曲半径（bending radius，R），8个结构参数受装置总体尺寸约束及预期工作效率等因素制约，取初值分别为150、350、500、700、200、50、10、100 mm[12]。
图2 杂交水稻花粉收集装置气力输送系统结构Fig. 2 Structure of pneumatic conveying system of hybrid rice pollen collection
1.3.2 网格划分 采用ICEM前处理软件作为网格划分工具对简化后的模型进行非结构性四面体网格划分。为避免因网格质量差影响仿真结果，对网格质量进行检查和优化，优化后的网格质量在0.4～1.0（最小值大于0.3则可视作网格较优[14]），网格数量为51 342，节点数60 921，划分结果如图3所示。
图3 气力输送系统数值模型计算域Fig. 3 Computational domain of numerical model of pneumatic conveying system
1.3.3 数值模拟 对花粉气力输送过程建立CFD-DPM模型计算，计算中将气相视为连续介质，求解欧拉坐标系下的流体控制方程；
把花粉颗粒视为离散介质，在拉格朗日坐标系下求解颗粒运动方程，采用迭代计算进行气固两相的双向耦合[18,23]。在众多CFD-DPM仿真工具中，ANSYS的Fluent因其优秀的网格再定义、强大的兼容性及计算能力而倍受青睐，因此，本研究使用Fluent软件进行流体仿真。
系统中气相可视为定常、恒温、不可压缩流体，考虑到实际工况下气力输送系统内气流流速较大，管内流动为湍流，故引用湍流模型对流体进行描述。相应的控制方程选取计算精度更高、实用广泛性更好的k-ε标准方程[24]。采用压力耦合方程组的半隐式算法和一阶离散格式对数值方程进行迭代计算，设置残差收敛值为1e-5，松驰因子调试优先选择较小值，根据模型是否收敛将松驰因子逐渐增加至默认值，其他项保持常规设置[25]。
离散相模型（DPM）可对花粉运动中所受流体力、固体力、场力进行受力平衡方程计算，通过求解颗粒运动方程[26]，得到气力输送系统内部固相的轨迹曲线。
1.3.4 参数设置 吸粉口端面部分为内外环境的分界面，气流运动较为复杂，压力、速度与流量均未知，因此在吸粉口端面处建立合适的扩展域，如图4所示，扩展域与外部区域直接相连，其压强即为1个标准大气压，以此确定气流入口边界条件为压力入口。气流出口边界条件根据1.2.2计算公式，设置为压力出口。具体仿真参数设置如表1所示。
表1 Fluent仿真参数设置Table 1 Fluent simulation parameter settings
图4 带入口扩展域的气力输送系统物理模型Fig. 4 Physical model of pneumatic conveying system with inlet expansion domain
1.4 单因素试验设计
气力输送系统吸粉均匀性及吸粉效率与系统各个结构参数息息相关。其中，进粉管内径（D）通过直接影响吸粉管内气流速度而影响系统内气流分布，吸粉口宽度（W）及吸粉口长度（L）对吸粉口吸粉效率及进入吸粉口内花粉颗粒速度有直接作用，相比之下本团队经预仿真试验结合前人研究结果发现吸粉口高度（H）、进粉管水平长度（A）、进粉管竖直高度（B）等对吸粉管内气流及花粉影响较小[17,27]，因此选择D、W、L作为本装置主要试验因素，研究其对系统内速度、压强分布均匀性及吸粉效率的影响。考虑到装置实际工况及自身尺寸约束，选取表2所示试验参数变化范围及水平。当D变化时，L、W分别固定为600、200 mm；
当L变化时，D、W分别固定为150、200 mm；
当W变化时，D、L分别固定为150、600 mm。其他参数保持初始值。
表2 主要因素变化范围及水平Table 2 Variation range and level of main factors
1.5 响应面优化试验设计
在单因素仿真试验的基础上，为进一步得到气力输送系统较优的结构参数组合，以较能反映系统内气流分布均匀性和气吸效率的速度变异系数[28]和全压差[17]作为仿真试验指标，采用BBC方法[29]进行仿真试验设计。仿真试验因素及水平如表3所示，试验方案与结果如表4所示。
表3 试验因素水平及编码Table 3 Levels and codes of experimental factors
表4 试验方案与结果Table 4 Experiment plan and results
1.6 优化结果仿真验证试验设计
为了获得气力输送系统最佳的结构参数组合，需要综合考虑进粉管内径（D）、吸粉口长度（L）、吸粉口宽度（W）对系统内气流速度（流速）变异系数和全压差的影响。在上述试验基础上，按照尽可能低的流速变异系数和全压差优化原则，利用Design-Expert10[30]进行优化。相应的约束条件如下。
2.1 初始结构仿真分析
对初始气力输送系统结构进行CFD-DPM仿真，仿真结果如图5和图6所示。由图5可知，系统内平均压强为-969.77 Pa，全压差为348.36 Pa，表明该结构下花粉收集输送过程能耗较大，局部损失和沿程损失较多；
大部分气流速度处于14.66～39.14 m·s-1，y=-550 mm（离底部扩展域下端面10 mm，为装置实际吸粉口）处平均速度为12.12 m·s-1，满足花粉收集时该处气流速度要求，但气流速度变异系数较大，为36.67%，说明系统内部速度分布不均匀性明显，加剧了花粉与花粉、花粉与壁面之间的碰撞从而影响了花粉活性。花粉运动轨迹与气流运动轨迹大体重合，大部分花粉速度大于悬浮速度（图6），说明花粉可快速通过气力输送系统进行装置主体。因此，为了降低气力输送系统内气流速度分布不均匀性，同时减少系统压力损耗从而提高吸粉效率，需对气力输送系统进行优化。
图5 气力输送系统中心截面气流场Fig. 5 Air field of central section
图6 气力输送系统内入射初速度（三维）为3.5 m·s-1时花粉运动轨迹Fig. 6 Pollen particles in pneumatic conveying system when the initial incident velocity is 3.5 m·s-1
2.2 单因素结果分析
通过单因素仿真试验得到图7所示不同D、L、W下系统内气流最大速度及平均速度变化情况，表5所示系统内气流速度均值、方差、速度变异系数以及全压差。在进粉管内径100～200 mm的变化范围内，随着内径的增加，气流最大速度、气流速度变异系数及全压差逐渐降低。这主要是由于随着D增加，进粉管与吸粉口连接处的截面积减小使得局部损失减少，因此全压差降低；
管内负压减小导致最大速度减小但气流均匀性增加，综合影响下，气流平均速度在直径为150～175 mm时增加，在150～200 mm区间内出现极值（为较优的D范围）。随着L、W的增加，气流最大速度逐渐降低，平均速度先增后减，全压差逐渐增加，这主要是由于随着L、W增加，吸粉口入口面积和气流量增加，加大了沿程损失（造成全压差增加）同时降低了吸粉口整体流速，但气流分布均匀性在一定范围内增加，因此平均速度先增后减，较优的L、W范围分别为500～700、150～250 mm。
表5 不同变量下系统速度变异系数与全压差Table 5 Variation coefficient and total pressure difference under different variables
图7 不同结构参数下气流速度分布Fig. 7 Airflow velocity distribution under different structural parameters
2.3 响应面结果分析
对响应面仿真试验结果进行统计学处理，得到如表6所示系统内气流速度变异系数、全压差的响应面方差分析结果。由表6可知,气流速度变异系数及全压差响应面方差分析结果模型P值均小于0.000 1，模型的失拟项均大于0.05，模型显著而失拟项不显著，说明模型成立且方程拟合较好。模型的一次项D、L、W的系数都是极显著，说明D、L及W对气流变异系数（气流均匀性）、全压差（气吸效率）均有显著影响，同时交互项DL、LW的交互作用对两指标影响显著。根据结果将不显著项去除后，得到以各因素编码值为自变量的系统内气流速度变异系数及全压差回归方程，分别为：
表6 响应面方差分析结果Table 6 Results of response surface analysis of variance
通过比较回归方程各因子系数的绝对值大小可知，对系统气流速度变异系数影响的主次关系依次为W2＞L2＞D＞D2＞L＞W＞DL＞LW；
对全压差影响的主次关系依次为D＞L＞W＞LW＞DL＞D2＞W2＞L2；
对两指标而言，一次项影响效果均为D＞L＞W。
进粉管内径（D）及吸粉口长度（L）对系统气流速度变异系数（CV）和全压差（Δp）的影响分别如图8所示，随着D的增大，CV和Δp逐渐减小；
随L的增加，Δp持续增大，CV则呈现先减小后增大的变化趋势。这是因为当D增加、L在一定范围内增加时，气流在系统内部流通面积增加且入口处的负压减小，通过吸粉口横截面的气体流速整体减小，系统内气流的均匀性增加，此外D变大使得进粉口拐角气流相对平缓，减小了此处的局部损失（避免了湍流漩涡现象的发生）;而当L过大，吸粉口边缘部分离系统中心较远，加大系统内沿程损失的同时极大降低了边界气流速度。两响应面中等高线图均为椭圆且更陡，说明DL的交互作用对CV及Δp影响显著，同时对Δp影响更大。由图8可知，当D为180～200 mm、L为550～650 mm时，CV有极小值；
当D为180～200 mm、L为500～600 mm时，Δp有极小值。
图8 回归模型响应曲面分析Fig. 8 Response surface analysis of regression model
吸粉口长度L和吸粉口宽度W对系统气流速度变异系数（CV）和全压差Δp的影响如图8所示，随着L和W的增大，CV先减小后增大,Δp逐渐增加。这是因为L和W的增加，加大吸粉口内部的沿程损失的同时极大增加了入口面积，使得单位时间进入的气流量变多，引起气流旋转速度增大，造成气流与内壁摩擦损失增加；
此外较大的L、W使吸粉口肩部距中心距离较远，部分区域负压较小速度较低，增加了系统内气流的不均匀性，因此CV随L、W增加呈现先减小后增大的变化趋势。图8显示等高线圆度较高，说明LW的交互作用尽管对CV有影响但影响优先级较低，图8中响应面等高线更密集，说明LW的交互作用对Δp有较显著的影响。由图8可知，当L为550～650 mm、W为190～210 mm时，CV存在极小值；
当L为500～600 mm、W为150～190 mm时，Δp存在极小值。
2.4 优化验证结果分析
响应面优化结果表明，当D为200 mm、L为564.40 mm、W为192.48 mm时，气流速度变异系数、全压差分别为15.28%，230.46 Pa，此时气流速度的均匀性较好，系统内局部压力损失较少，吸粉效率较高。
对优化后系统进行仿真，得到气流速度变异系数和全压差分别为16.03%、238.37 Pa，与预测结果的相对误差分别为4.91%、3.39%（均小于5%），误差较小（表7）；
较初始结构下两指标分别减少了56.29%、31.57%，说明吸粉效率及吸粉均匀性明显提升，具有较好的优化效果。将优化后结构参数代入式（1）、（2），可以选定离心风机型号为4-72-4A。
表7 仿真验证结果Table 7 Simulation verification results
制种机械化是我国杂交水稻产业发展的必然选择。一直以来，我国研究推广的机械化授粉模式主要有2种，一种是模拟人力授粉，在父母本分行种植的田块内采用机械碰撞或碰撞气吹（气吹）的方式让父本花粉落至母本，完成授粉[31-34]；
另一种是运用采粉-贮存-喷粉模式授粉，即父、母本分田种植，通过机械采粉、授粉以及与此相配套的花粉贮藏技术等，完成授粉作业[6-7]。前者操作简单，经优化后效率较人力式有所提升，但受双亲花期影响大,机械化移栽及田间管理难度大；
后者为完全意义上的机械化制种，不受父母本花期影响,可实现异时异地授粉，同时父、母本分田种植可在移栽收获和中期管理环节实现机械化操作（提高制种效率）。近年来随着农业劳动力成本的增加、相关技术的发展以及杂交水稻制种产业发展的需要，基于采粉-贮存-喷粉的全程机械化制种模式受到越来越多学者的关注[4-5]。而高效的采粉装置是后续贮粉以及授粉的基础，因此对设计的气吸式杂交水稻花粉收集装置进行气力输送系统优化，可为现阶段采粉-贮存-喷粉模式研究中的杂交水稻采粉环节提供花粉高效采集装置。
CFD-DPM是在欧拉-拉格朗日坐标系下对气固两相流建立连续-离散模型进行计算，适用于固相体积分数小于10%的情况[17]，不少研究者通过试验证实了运用CFD-DPM模拟两相流动，优化装置工作参数的可靠性[19,35]。花粉收集过程为典型气固两相流动，且花粉颗粒粒径小，在两相流中为稀相（体积分数小于10%），因此采用CFD-DPM模型对气力输送系统内气固两相进行仿真分析，结合响应面试验法对系统结构参数进行优化。气力输送系统结构参数中进粉管水平长度（A）、进粉管竖直高度（B）主要由车身尺寸和实际工况确定，扩展域高度（H0）结合徐云等[36]的研究结果固定为10 mm，本团队经预仿真试验得到吸粉口高度对系统内气流及花粉影响较小，结合整体尺寸约束及李炽[27]的研究结果将其固定为50 mm，综合考虑将进粉管内径（D）、吸粉口长度（L）、吸粉口宽度（W）作为试验因素，以较能反映系统内气流分布均匀性和气吸效率的流速变异系数和全压差作为优化指标进行仿真分析及优化试验。仿真结果表明，气吸效率和吸粉均匀性均随内径的增加而增加，吸粉均匀性随长度和宽度的增加先增加后减小，全压差随长度和宽度的增加而增加，这与彭宇源[17]和杨春朝等[37]的研究较接近。优化结果表明，花粉气力输送系统较优的结构参数组合为内径200.00 mm、长度564.40 mm、宽度192.48 mm，对优化后结构进行仿真，得到气流速度变异系数和全压差分别为16.03%、238.37 Pa，与预测结果的相对误差分别为4.91%、3.39%，较初始结构下两指标分别减少了56.29%、31.57%，说明吸粉效率及吸粉均匀性明显提升，具有较好的优化效果。
本研究通过仿真结合响应面法优化气力输送系统结构，一方面提高了花粉收集效率和吸粉均匀性，另一方面缩短设计优化周期同时降低了研制成本，为花粉收集装置气力输送系统快速优化设计提供了参考。但仅通过仿真验证了优化结果，没有进行样机试验，下一步将在前期设计以及仿真优化的基础上指导样机加工，完成样机田间试验证实优化结果的可靠性。此外，气力输送系统的外形、风机结构、颗粒进入速度等也是影响花粉收集效果的重要因素，还需要进一步研究其他因素对吸粉效果的影响规律；
在优化指标选择方面，考虑到花粉易失活，花粉收集后活性应是衡量吸粉效果的重要指标，由于本文优化的主要目标是高效收集，故只引入了吸粉均匀性这一间接指标研究各参数对花粉活性的影响，而未选择可定量描述的直接指标，后续研究可将花粉活性作为参数优化的重要指标。