基于Porous模型的多旋翼植保无人机下洗气流分布研究

张 豪 祁力钧 吴亚垒 刘婠婠 程浈浈 MUSIUE

(中国农业大学工学院， 北京 100083)

0 引言

农业航空植保具有安全、高效和应急防治能力强等诸多优势[1-2]，可有效保障生态与粮食安全[3]，近年来得到了大力发展[4]。已有研究表明，采用植保无人机施药时，旋翼下洗气流是影响雾滴运动的关键因素，且对雾滴在作物冠层内部的穿透性有显著影响[5-6]。目前，果树冠层内部及周围旋翼下洗气流的时空分布尚不明确，有必要针对采用多旋翼植保无人机进行果树施药作业时旋翼下洗气流的时空分布展开研究。

围绕植保无人机授粉和施药技术，已有较多田间试验[7-8]和数值模拟[9-11]研究。ZHENG等[12]、王昌陵等[13]研究了植保无人机不同操作参数下雾滴的沉积分布规律，许童羽等[14]研究发现雾滴的沉积分布规律不仅与无人机作业参数有关，还受旋翼下洗气流和外界风场的影响。田间试验无法全面获取旋翼下洗气流的三维数据，结合计算流体力学方法，薛新宇团队对单旋翼[15]、六旋翼[16]植保无人机悬停状态的下洗气流进行了数值模拟与试验研究，但在数值模拟过程中都未考虑作物冠层对下洗气流的影响。另外，祁力钧等[17]、DELELE等[18]利用计算流体力学(Computational fluid dynamic，CFD)技术，针对果园风送式喷雾机无果树冠层的雾滴分布特性进行了研究。HONG等[19-20]、DUGA等[21-22]、SALCEDO等[23]则在数值模拟过程中引入Porous模型处理果树冠层，开展果园风送式喷雾机气流和雾滴分布规律数值模拟与试验验证研究，这些研究充分表明利用Porous模型处理果树冠层进行数值模拟研究的有效性。运用CFD技术，基于Porous模型的多旋翼植保无人机悬停施药果树冠层周围及内部下洗气流时空分布方面的研究未见报道。

为明确多旋翼植保无人机悬停施药时果树冠层周围及内部下洗气流时空分布规律，本文基于商用软件ANSYS Fluent 16.0，采用数值模拟与试验验证相结合的方法，针对六旋翼植保无人机，结合RANS方程、RNGk-ε湍流模型、Porous模型、滑移网格技术及SIMPLE算法，建立六旋翼植保无人机悬停施药下洗气流时空分布的三维CFD模型，并搭建数值模拟试验平台进行旋翼下洗气流速度测试试验。

1 材料与方法

1.1 数值模拟试验平台

试验平台(中国农业大学植保机械实验室)如图1所示。该试验平台包括可移动支架、悬挂装置和六旋翼植保无人机。所用六旋翼植保无人机定制于深圳市金铭睿电子有限公司，基本参数为：无人机轴距0.8 m，电机型号X4114 KV370，旋翼型号1555，最大载药量为5 kg。植保无人机主要包括机臂、中间盘、旋翼、动力系统和无线遥控器，加装有RC41智能测速仪，反馈旋翼转速。本文考虑载药量以及室内试验的安全性，旋翼转速取2 500 r/min，载药量可达2.5 kg。

图1 多旋翼植保无人机数值模拟试验平台Fig.1 Multi-rotor plant protection UAV numerical simulation experimental platform1.可移动支架 2.六旋翼植保无人机及悬挂装置 3.试验果树

1.2 多旋翼植保无人机气动布局

多旋翼植保无人机下洗气流时空分布与其旋翼系统的气动布局密切相关，本文研究对象六旋翼系统气动布局如图2所示，相邻旋翼夹角为60°。OiXiYiZi为固连在第i个旋翼上的旋转坐标系(i=1,2,…,6)，ObXbYbZb为植保无人机悬停状态下绝对坐标系。悬停状态下6个旋翼转速相同，相邻旋翼转向相反。

图2 六旋翼系统气动布局Fig.2 Aerodynamic layout of six-rotor system

1.3 数值模拟方法

1.3.1基本控制方程

六旋翼植保无人机旋翼悬停流场采用RANS方程作为基本控制方程。对于单个旋翼，将旋转坐标系的坐标轴固连在旋翼上，其守恒积分控制方程具体可表示为[24-25]

(1)

其中

式中S——控制体表面积

n——控制体外法矢量

V——控制体体积

W——守恒变量

F——无粘通量

Fv——粘性通量

R——旋转通量

q——绝对速度

qw——牵连速度

Er——流体总内能

Hr——流体总焓

τxx、τyy、τzz、τxy、τxz、τyz——粘性通量Fv中粘性应力τ各分量

ix、iy、iz——控制体外法矢量n的分量

φx、φy、φz——粘性项

p——流体压强，Pa

ρ——流体密度，kg/m3

t——时间，s

Ω——旋翼旋转角速度矢量

果树冠层的阻力作用，会造成旋翼下洗气流的动量损失，本文在总结前人研究成果[26-30]后选择Porous模型处理果树冠层，用Porous域代替果树冠层区域。该模型通过在流体基本动量方程中增加动量损失源项来模拟Porous域对气流的阻力作用，冠层所造成的气流动量损失和湍流量变化可描述为[19]

(2)

Sk=CdLADβp|v|3-CdLADβd|v|k

(3)

(4)

式中 Δp——冠层压力损失，kg/(s2·m)

Cir——压力损失系数，m-1

v——风速，m/s

Δm——多孔介质厚度，m

k——湍流动能，J

ω——比耗散率，%

Sk——湍流动能k的源项，m2/s2

Sω——比耗散率ω的源项，s-1

Cd——冠层阻力系数，取0.25

LAD——叶面积密度，m-1

βp——平均流体动能由于作物阻力转化为湍流动能损失的比例系数，取1

βd——与作物相互作用的能量损失系数，取4

αp、αd——模型常数，取1.5、1.5

1.3.2旋翼模型与模拟果树结构

旋翼是植保无人机关键旋转部件，所用三维模型对数值模拟结果有较大影响，本文依托逆向工程技术，利用MCS五/四轴全自动三维扫描系统分别获取正、反旋翼表面点云数据，而后用Geomagic Studio软件对点云数据进行后处理，实现旋翼模型的逆向重建，得到精确的旋翼三维模型(图3)，用于数值模拟。

图3 旋翼逆向建模Fig.3 Rotor reverse modeling

本文数值模拟时参考文献[20]对试验果树进行简化处理，将模拟果树用球形冠层与圆柱树干的组合来代替；冠层压力损失系数与冠层叶面积密度关系式为Cir=2CdLAD，可用冠层压力损失系数表征冠层叶面积密度；一般叶面积密度冠层和高叶面积密度冠层的压力损失系数分别为0.8～6.9和1.2～13.4。综上，设置Ⅰ、Ⅱ、Ⅲ共3棵模拟果树，特征参数包括结构参数和压力损失系数，其中结构参数为：球形冠层直径1.2 m，圆柱树干直径0.15 m、高度0.7 m；Ⅰ、Ⅱ、Ⅲ号果树压力损失系数分别为2.0、6.0、10.0，依次表征冠层的稀疏、较密和稠密。结构参数由测量实际试验果树尺寸得到。

1.3.3计算区域与边界条件

总计算区域为直径4.5 m、高3.6 m的圆柱体(图4)，植保无人机、总计算区域上边界和下边界中心坐标分别为(0，0，0)、(0，0，1.2)、(0，0，2.4)，模拟果树垂直中轴线与Z轴重合。

为保证数值模拟的求解精度，将总计算区域划分为6个旋转域(设置转速为2 500 r/min)、冠层Porous域、树干固体域及空气外域，共9个计算区域，其中旋转域、冠层Porous域与空气外域之间的交界面采用Interface连接。进行无果树模拟时不包含冠层Porous域和树干固体域。无果树和有果树模拟总网格数分别为5 897 064和7 342 928；网格质量评价参数Skewness最大值分别为0.849 63和0.841 66，平均值分别为0.230 24和0.229 69，符合网格质量要求[31]。由于下洗气流场求解为瞬态计算，旋转域采用滑移网格技术处理，空气外域下边界为地面，采用无滑移壁面，其余边界为自由出口。

1.3.4计算方法

植保无人机旋翼高速旋转形成下洗气流场，湍流模型选择适合于复杂剪切流动的RNGk-ε模型[32]，使用有限体积法将控制方程离散化。数值模拟时动量和压力为所关注变量，选择基于压力求解器的SIMPLE压力-速度耦合算法对控制方程进行数值求解。

1.4 旋翼下洗气流速度测试

为验证数值模拟旋翼下洗气流场的准确性，进行无果树下洗气流速度测试试验，另外为说明模拟果树与试验果树冠层内部气流速度分布规律的一致性，进行有果树下洗气流速度测试试验，旋翼转速稳定在(2 500±10) r/min(图5)。各测试点气流速度用固定于三角支架上的Kestrel 4500型风速仪依次进行测量，无自然风。风速仪设置每2.0 s保存一个数据，下洗气流稳定后开始测试，每个测试点的单点测试时间约20 s，导出数据，取所记录速度的平均值为该测试点气流速度。

图5 下洗气流速度测试试验Fig.5 Down-wash airflow speed test

无果树下洗气流速度测试点布置在各旋翼中心的正下方，分别位于旋翼正下方0.3、0.8、1.3、1.8、2.3 m处，共30个测试点，依次进行Z方向速度测试；图6给出了冠层内部气流速度测试点空间位置，将冠层分为上、中、下3层，距地面高度分别为1.6、1.3、1.0 m，共15个测试点，同样依次进行Z方向速度测试。

图6 冠层内部气流速度测试点示意图Fig.6 Schematic of measuring points for airflow speed distribution inside tree canopy

2 结果与讨论

2.1 数值模拟结果与分析

2.1.1旋翼下洗气流时变特征

CFD数值模拟的优势之一在于三维结果的可视化，图7分别给出了无果树和有果树时，不同时刻旋翼下洗气流总速度的分布。

无果树时，1.0 s已有速度较小的旋翼下洗气流接触地面但尚未沿地面向四周形成铺展；2.0 s时地面铺展已经形成，但与3.0 s和3.5 s比较，其铺展面积和速度都偏小，表明此时刻旋翼下洗气流已经发展至地面，且沿地面向四周的铺展正处在发展过程中；3.0 s和3.5 s旋翼下洗气流总速度分布情况基本一致，气流在计算区域内的地面铺展得到了充分发展，表明3.0 s时旋翼下洗气流场已趋于稳定。整个过程旋翼下洗气流近似呈“圆柱形”向下发展，到达地面后形成地面铺展，且与文献[16]研究结果一致，Z方向气流速度在旋翼下洗气流速度中占主体。

当存在果树时，旋翼下洗气流场与无果树时不同，果树冠层对旋翼下洗气流有明显的阻挡作用，气流发展至地面所需时间更长，地面铺展在计算区域内表现不明显，且随着冠层压力损失系数的增加该现象更为突出。Ⅰ、Ⅱ、Ⅲ号果树旋翼下洗气流在1.0 s时发展至冠层上半部，总速度分布情况基本相同，由于冠层的阻挡作用，此时刻在冠层上半部区域的周围均出现气流卷扬，且该卷扬在Ⅰ号果树中随气流的向下发展在冠层周围不再明显，逐渐发展成为近地面卷扬，在Ⅱ、Ⅲ号果树中随气流的向下发展该卷扬有增强趋势。2.0 s时旋翼下洗气流已穿透Ⅰ号果树冠层向地面发展，Ⅱ、Ⅲ号果树气流穿透深度较Ⅰ号果树小，这是由于冠层压力损失系数的增加造成的。3.0 s和3.5 s时，Ⅰ号果树旋翼下洗气流已发展至地面，冠层周围气流从冠层上半部区域开始呈“圆锥形”向下发展，以一倾斜角发展到地面形成小范围地面铺展，地面铺展末端出现近地面卷扬。由于Ⅱ、Ⅲ号果树的压力损失系数均较大，3.0 s和3.5 s时旋翼下洗气流总速度分布相似，冠层周围气流卷扬严重，这将增大雾滴的飘移，对果树施药不利，且计算区域内无明显地面铺展。

2.1.2旋翼下洗气流速度空间变化规律

为明确无果树和有果树时旋翼下洗气流速度的空间变化规律，选择3.5 s稳定流场模拟结果进行分析。图8分别给出了无果树和有果树时，3.5 s时刻各旋翼中心正下方的下洗气流Z方向速度沿Z轴的变化情况。

图8 各旋翼中心正下方下洗气流Z方向速度Fig.8 Down-wash airflow speed in Z-direction below center of each rotor

无果树时，旋翼中心正下方的下洗气流Z方向速度在0.2 m处最大，接近8.0 m/s，这是因为对于单个旋翼在近旋翼中心正下方存在一个小的“低速区”，旋翼中心外围下洗气流速度较大且气流有一个加速过程，随着外围气流的交汇叠加，“低速区”结束，在旋翼中心正下方0.2 m处气流速度达到最大。由于空气阻力作用，0.2～0.6 m区域内Z方向速度迅速由近8.0 m/s衰减至4.0 m/s，而后在0.6～1.7 m区域内出现一个“Z方向速度稳定区”，速度处在3.0～4.0 m/s范围内，速度衰减较小，表明该区域下洗气流呈稳定发展。1.7～2.4 m区域内由于下洗气流发展至地面，地面的阻挡作用造成Z方向速度逐渐衰减至0 m/s。

当存在果树时，在0～0.6 m区域内，下洗气流尚未受到冠层的影响，其Z方向速度变化规律与无果树时基本一致。在0.6～1.7 m区域内，由于存在冠层的阻挡作用，造成下洗气流的动量损失，不再出现“Z方向速度稳定区”，随着冠层压力损失系数的增加，旋翼下洗气流Z方向速度衰减加快。另外，Ⅰ、Ⅱ、Ⅲ号果树的冠层压力损失系数依次增大，其造成的旋翼下洗气流速度衰减同样依次增大，表现为Ⅰ号果树冠层下方Z方向速度大部分处在1～2 m/s范围内，而Ⅱ、Ⅲ号果树冠层下方Z方向速度较Ⅰ号果树小，均接近0 m/s，说明下洗气流在Ⅱ、Ⅲ号果树冠层中穿透性较Ⅰ号果树差。Ⅰ、Ⅱ、Ⅲ号果树各旋翼中心正下方Z方向速度曲线出现重合，这可能是由于Porous模型被假设为各向同性，且各旋翼在冠层下半部和冠层下方的Z方向速度较小造成的。

冠层的存在造成Z方向速度衰减加快，同时旋翼气流向四周产生扩散，图9给出了距地面高度1.6、1.3、1.0 m下，有Ⅰ、Ⅱ、Ⅲ号果树时，下洗气流X、Y方向速度分别沿X、Y轴的变化情况，其中速度值的正负分别表示该速度方向指向坐标轴的正方向和负方向。从图中可以看出，X、Y方向速度绝对值均呈对称分布；Ⅰ、Ⅱ、Ⅲ号果树冠层周围X方向和Y方向速度基本呈依次增大趋势，说明随着冠层压力损失系数的增大，冠层对旋翼气流的扩散作用增强，这将增加雾滴的横向飘移；同一果树、同一高度下，X方向速度较Y方向速度大，这可能是因为X轴处在文献[16]中描述的气流“引入”、“导出”区，受气流叠加的影响，而Y轴处在旋翼正下方不存在气流的叠加，但距地面越近，X、Y方向速度差异受气流“引入”、“导出”区气流叠加的影响越小；另外，距地面1.6 m处X、Y方向速度变化较快，气流向四周的扩散范围小，距地面1.3 m和1.0 m处X、Y方向速度变化相对缓慢，气流向四周的扩散范围增大，符合2.1.1节中冠层周围气流从冠层上半部区域开始呈“圆锥形”向下发展的特点。

图9 有果树不同高度下洗气流X和Y方向速度Fig.9 Down-wash airflow speed in X and Y directions at different heights with tree

2.1.3冠层内部下洗气流速度分布

按照有果树下洗气流速度测试试验，同样以距地面高度分别为1.6、1.3、1.0 m将模拟果树分为上、中、下3层，分别获取冠层内部各层下洗气流Z方向最大速度。

表3给出了无果树和有Ⅰ、Ⅱ、Ⅲ号果树时，不同时刻冠层内部上、中、下各层下洗气流Z方向最大速度，并给出了3.5 s时刻Z方向最大速度的衰减比，该指标用于衡量旋翼下洗气流Z方向速度沿Z轴正方向的衰减幅度，若衰减比等于1说明速度无变化，大于1说明速度出现衰减，衰减比越大速度衰减幅度越大。其计算公式为

(5)

式中Ni——Z方向最大速度衰减比

v1——冠层内部上层或中层下洗气流Z方向最大速度，m/s

v2——冠层内部下层下洗气流Z方向最大速度，m/s

有果树时N1、N2分别指冠层上半部和下半部的Z方向最大速度衰减比；无果树时N1、N2分别指与有果树时冠层相同部位的Z方向最大速度衰减比。

从表3可以看出，各时刻均呈现出“上层-无果树”Z方向最大速度最大，“下层-Ⅲ号果树”Z方向最大速度最小的变化规律，表明随着冠层深度及冠层压力损失系数的增加，冠层内部Z方向最大速度均呈减小趋势。在同一层上，无果树、Ⅰ、Ⅱ、Ⅲ号果树Z方向最大速度随时间逐渐增大，但“3.5 s-中层-Ⅰ号果树”Z方向最大速度较“3.0 s-中层-Ⅰ号果树”的Z方向最大速度小，这可能是冠层内部气流受湍流的影响造成。3.0 s和3.5 s时刻，Ⅰ、Ⅱ、Ⅲ号果树中层和下层Z方向最大速度差距较小，可以说明3.5 s时刻下洗气流在冠层内部已得到充分发展。

计算Z方向最大速度衰减比发现，无果树和Ⅰ、Ⅱ、Ⅲ号果树Z方向最大速度衰减比依次增大(除Ⅲ号果树冠层下半部)，表明Z方向最大速度衰减幅度依次增大，出现Ⅲ号果树冠层下半部衰减比较Ⅱ号果树的小，说明冠层下半部Z方向最大速度衰减比与冠层压力损失系数之间无线性关系。无果树时冠层上半部和下半部衰减比分别为1.07和1.05，该区域处在“Z方向速度稳定区”。Ⅰ号果树冠层上半部、下半部衰减比分别为1.30和1.36，冠层下半部Z方向速度衰减幅度比冠层上半部的衰减幅度大。Ⅱ号果树和Ⅲ号果树冠层上半部Z方向最大速度衰减比分别为2.41和2.80，下半部分别为1.90和1.85，冠层上半部衰减比都比下半部的大，表明Ⅱ、Ⅲ号果树冠层上半部Z方向最大速度衰减幅度更大。

2.2 试验结果与分析

为保证旋翼下洗气流在冠层内部已充分发展，选择3.5 s时刻的数值模拟结果与实际试验结果进行对比分析。

2.2.1旋翼下洗气流场准确性验证与分析

无果树下洗气流速度测试试验结果和模拟结果见表4。沿Z轴正方向，下洗气流Z方向速度试验值与模拟值有相同变化趋势，试验结果表明，在旋翼正下方0.8、1.3、1.8 m处各测试点Z方向速度相差较小，衰减幅度小，同样存在“Z方向速度稳定区”。0.3、0.8、1.3、1.8 m处各测试点模拟值和试验值相对误差均在10%以内，2.3 m处各测试点模拟值和试验值相对误差不大于25%，造成此差异的原因可能是2.3 m处下洗气流已接近地面且向四周形成铺展，这增加了实际测量过程中误差增大的可能性。另外，2.3 m处各测试点Z方向速度试验值都比模拟值要小，主要原因是试验过程中的动量损失比模拟过程中的大，这对距旋翼较远处的下洗气流Z方向速度影响最大。虽然在近地面2.3 m处模拟值和试验值的相对误差较大，但考虑到果树冠层距此区域较远，冠层所在区域气流场模拟值和试验值相对误差较小，在10%以内，与文献[16]结果相比误差有所降低。其次，对试验值和模拟值进行线性回归分析(图10)，回归拟合方程为y=1.010 2x+0.006 3，总体拟合优度R2为0.984 6，试验值和模拟值显著相关。综上说明旋翼下洗气流场数值模拟结果准确。

表4 无果树下洗气流Z方向速度模拟值与试验值对比Tab.4 Comparison of test and simulation values of down-wash airflow speed in Z-direction without tree

图10 无果树下洗气流Z方向速度试验值与模拟值回归分析Fig.10 Regression analysis of Z-direction down-wash airflow speed test values and simulation values without tree

2.2.2冠层内部气流速度分布规律一致性验证与分析

将有果树下洗气流速度测试试验的15个测试点数据和对应数值模拟数据按照测试位置上、中、下层分为3组，获取冠层内部各层Z方向气流速度范围，并计算各层速度分布变异系数，该指标用于衡量各层Z方向速度分布的均匀性，计算公式为[33]

(6)

其中

(7)

式中Cv——Z方向速度分布变异系数，%

σ——各层Z方向速度标准差，m/s

Vi——各测试点Z方向速度测量值，m/s

n——各层Z方向速度测试点个数，取5

图11给出了Ⅰ、Ⅱ、Ⅲ号果树和试验果树冠层内部各层Z方向气流速度范围及速度分布变异系数。对比分析可以发现，试验果树与模拟果树冠层内部气流速度分布规律具有很好的一致性。

从速度范围看，各果树上、中、下层Z方向速度范围依次减小。试验果树与Ⅲ号果树上、中、下层最大速度及上层最小速度较为接近，但中层和下层的最小速度二者存在较大差异，这可能是因为在有果树下洗气流速度测试试验过程中由于树叶扰动，风速仪被瞬间遮挡，造成风速仪在该时刻所记录的数值较小，使得测试点平均速度值偏低，而果树上层的气流速度较大，树叶扰动对上层测量结果影响较小。从速度分布变异系数来看，各果树上、中、下层Z方向速度分布变异系数依次减小，说明随着冠层深度的增加，Z方向速度分布趋于均匀。另外，综合速度范围和速度分布变异系数，试验果树与Ⅲ号果树情况更为接近，可以估计试验果树压力损失系数在10左右。

图11 模拟果树和试验果树各层气流速度范围及变异系数Fig.11 Airflow speed range and coefficient of variation for each layer of simulated tree and experimental tree

3 结论

(1)建立了轴距0.8 m、旋翼转速2 500 r/min的六旋翼植保无人机无果树和有果树的旋翼下洗气流时空分布三维CFD模型，对比分析了无果树和Ⅰ、Ⅱ、Ⅲ号果树旋翼下洗气流的时变特征。无果树时，旋翼下洗气流近似呈“圆柱形”向下发展，2.0 s时已到达地面且形成地面铺展；有果树时，冠层周围气流从冠层上半部区域开始呈“圆锥形”向下发展，以一倾斜角发展到地面形成小范围地面铺展；Ⅱ、Ⅲ号果树冠层周围气流卷扬严重。

(2)研究了3.5 s时无果树和有果树时旋翼下洗气流速度空间变化规律。无果树时，在旋翼正下方0.6～1.7 m区域内出现速度范围为3.0～4.0 m/s的“Z方向速度稳定区”；有果树时，冠层对旋翼下洗气流有明显的阻挡作用，不再出现“Z方向速度稳定区”；旋翼中心正下方Z方向速度最大接近8 m/s，随着冠层压力损失系数的增大，Z方向速度衰减加快，冠层对旋翼气流的扩散作用增强。

(3)讨论了冠层内部下洗气流速度分布。3.5 s时冠层内部气流已发展充分；不同时刻冠层内部Z方向最大速度，随着冠层深度及冠层压力损失系数的增加而减小；除Ⅲ号果树冠层下半部，无果树和Ⅰ、Ⅱ、Ⅲ号果树冠层内部Z方向最大速度衰减幅度依次增大；冠层下半部Z方向最大速度衰减比与冠层压力损失系数之间无线性关系。

(4)开展了无果树和有果树的下洗气流Z方向速度测试试验：旋翼正下方0.3、0.8、1.3、1.8 m处和近地面2.3 m处试验值与模拟值相对误差分别在10%以内和不大于25%，总体拟合优度为0.984 6，验证了旋翼下洗气流场数值模拟结果的准确性；试验果树冠层内部Z方向速度随着冠层深度的增加逐渐减小，且速度分布趋于均匀，与模拟果树冠层内部气流速度分布规律具有很好的一致性。