果园风送喷雾机出风口风场CFD建模与试验

doi:10.12133/j.smartag.2021.3.3.202106-SA007

果园风送喷雾机出风口风场CFD建模与试验

翟长远^,¹^,², 张燕妮¹^,³, 窦汉杰¹^,²^,³, 王秀¹^,², 陈立平^,¹^,²^,³

1.北京市农林科学院智能装备技术研究中心，北京 100097

2.国家农业智能装备工程技术研究中心，北京 100097

3.西北农林科技大学机械与电子工程学院，陕西杨凌 712100

CFD Modeling and Experiment of Airflow at the Air Outlet of Orchard Air-Assisted Sprayer

ZHAI Changyuan^,¹^,², ZHANG Yanni¹^,³, DOU Hanjie¹^,²^,³, WANG Xiu¹^,², CHEN Liping^,¹^,²^,³

1.Beijing Research Center of Intelligent Equipment for Agriculture, Beijing 100097, China

2.National Engineering Research Center of Intelligent Equipment for Agriculture, Beijing 100097, China

3.College of Mechanical and Electronic Engineering, Northwest A&F University, Yangling 712100, China

通讯作者: 陈立平（1973－），女，研究员，研究方向为农林业智能装备技术。电话：010-51503425。E-mail：chenlp@nercita.org.cn。

收稿日期: 2021-06-11 修回日期: 2021-07-08

基金资助:

国家自然科学基金项目 (31971775)
重庆市技术创新与应用发展专项 (cstc2019jscx-gksbX0089)

Received: 2021-06-11 Revised: 2021-07-08

作者简介 About authors

翟长远（1982－），男，研究员，研究方向为精准施药技术。E-mail：zhaicy@nercita.org.cn。 E-mail：zhaicy@nercita.org.cn

摘要

塔式果园风送喷雾机是目前普遍使用的果园喷雾机，塔式结构易产生旋转和不规则的垂直气流，导致喷雾机出风口两侧风场分布不对称，且分布规律不易预测。为探究适用塔式果园风送喷雾机出风口风场建模方法，本研究基于喷雾机出风口风场计算流体动力学（Computational Fluid Dynamics，CFD）建模方法，提出了用户定义函数（User-Defined Function，UDF）分段式三维风速入口边界条件设置方法，并研究了湍流模型和计算域尺寸对喷雾机CFD风场模拟结果的影响特性。采用Fluent软件，建立了三种CFD模型：模型一以11个区域的平均风速作为边界条件；模型二采用UDF分段式三维风速入口作为边界条件；为进一步研究计算域尺寸对风场模拟的影响，建立了小计算域尺寸的模型三。三种模型均采用基于雷诺时均控制（Reynolds-Averaged Navier-Stokes，RANS）方程的 $k - ε$ 湍流模型和 $k - ω$ 湍流模型进行风场计算。为了验证模型的可靠性，设计了空间风场立体测量系统，实现了精确快速空间风速测量。验证结果表明，Standard $k - ε$ 、Realizable $k - ε$ 、BSL $k - ω$ 和SST $k - ω$ 湍流模型较适合风场CFD建模，其中Standard $k - ε$ 湍流模型运算结果最优，模型决定系数R²为0.81。基于UDF分段式三维风速入口边界条件设置方法不仅提高了仿真结果的准确性（提高了5%），而且降低了计算的复杂性。在其他参数设置相同的情况下，大尺寸计算域模型的性能略优于小尺寸计算域。实际建模过程中建议根据计算机计算能力、仿真的实际要求设置合适的计算域尺寸。本研究结果可为喷雾机出风口风场CFD建模方法提供参考。

关键词： CFD ; 边界条件 ; UDF ; 湍流模型 ; 计算域

Abstract

The tower-type sprayer produces swirling and irregular vertical airstream. The complex swirling results in airflow asymmetry between sides of the sprayer, and the vertical air velocity profile can be unpredictable when the rotational speed of the fan changes. The spray deposition is directly linked to the airflow pattern obtained from the sprayers. In order to study airflow field of this type of air-assisted sprayer, a CFD (Computational Fluid Dynamics) model for the tower-type sprayer was developed. A boundary condition setting method of UDF (User-Defined Function) sectional 3D air velocity was proposed. And the influences of turbulence models and the size of computational domain on CFD airflow simulation were studied. Using Fluent software, three different CFD models were established. The Model 1 took the average air velocity of 11 regions as the velocity inlet. The Model 2 used UDF segmented three-dimension air velocity line as the boundary condition. In order to further study the influence of the computational domain size on simulation, the Model 3 with a smaller computational domain was established. The turbulence model based on reynolds-averaged navier-stokes (RANS) control equation was used to calculate the airflow field in all models. In order to verify the reliability of the model, a three-dimensional measurement system of airflow field was designed, which was used for accurate and fast velocity measurement. The results showed that the Standard k-ε turbulence model, Realizable k-ε turbulence model, BSL k-w turbulence model, SST k-w turbulence model were suitable, and the Standard k-ε turbulence model was the best one. The CFD boundary condition setting method of UDF sectional three-dimension air velocity could improve the accuracy of simulation, and reduce the calculation complexity. With the same settings of other parameters, the performance of the CFD model with larger scale calculation domain was slightly better than that with smaller computational domain. The size of computational domain should be set to the appropriate extent, considering the calculation capacity and practical requirements of modelling. The research results could provide an important reference for CFD modeling of spray airflow field.

Keywords： CFD ; boundary condition ; UDF ; turbulence model ; computational domain

PDF (2094KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

翟长远, 张燕妮, 窦汉杰, 王秀, 陈立平. 果园风送喷雾机出风口风场CFD建模与试验[J]. 智慧农业（中英文）, 2021, 3(3): 70-81. doi:10.12133/j.smartag.2021.3.3.202106-SA007

ZHAI Changyuan, ZHANG Yanni, DOU Hanjie, WANG Xiu, CHEN Liping. CFD Modeling and Experiment of Airflow at the Air Outlet of Orchard Air-Assisted Sprayer[J]. Smart Agriculture, 2021, 3(3): 70-81. doi:10.12133/j.smartag.2021.3.3.202106-SA007

1 引言

果园的病虫害防治长期依靠化学农药，据统计，如停止使用化学农药，果园将损失66%~90%的果品产量^［1］。果园风送喷雾技术作为联合国粮食及农业组织（Food and Agriculture Organization of the United Nations，FAO）推荐的施药技术，具有药液贯穿能力强、叶片正反面沉积均匀的特性^［2,3］，在中国果园逐渐得到了应用推广。但目前果园施药存在大量农药飘移，不仅浪费农药、污染环境，更对果品安全生产带来巨大威胁。果园风送喷雾风力的有效调控是降低雾滴飘移、提高农药利用率的关键^［4-6］。Landers^［7］通过试验发现风力调节可提高30%的农药沉积、减少75%药液漂移。风力调控主要是通过改变喷雾机风箱进出风口面积、风箱角度、风机转速等方法以获取最佳的喷雾风场^［8-10］。

风送喷雾机在不同调控条件下的空间风场模型构建是进行风力精准调控的前提，空间风场的获得可以通过田间实际测量完成^［11-13］。Hołownicki等^［14］对关于风场的工作进行了梳理，研究者们为了确定风场参数对雾滴沉积的影响，进行了大量的风场试验与研究，但大空间多风场测量过程时间和经济成本巨大，再加上自然风方向和大小的骤变，田间试验收集的风场数据，对不同的环境条件具有不可重复性。近年来，研究者进行了大量的计算流体动力学（Computational Fluid Dynamics，CFD）风场建模方法研究，模拟不同类型喷雾机的外部风场，结果表明CFD风场模型可以准确模拟不同类型喷雾机的外部风场^［15-17］。对于特定类型的喷雾机，经过试验数据验证的风场模型将是一种便捷、可有效分析喷雾机外部风场的工具^［18-20］。

CFD风场建模主要包含了模型初始条件、控制方程与湍流模型的选择，以及设计合理的计算域3个方面的内容^［21］。Salcedo等^［22］采用简化的直线代替喷雾机出风口，采用基于雷诺时均控制（Reynold-Averaged Navier-Stokes，RANS）方程的湍流模型建立了CFD风场2D模型。Dekeyser等^［23］采用矩形截面和弯曲截面模拟横流式风机和轴流式风机出风口，尽可能靠近出风口测量出风口风速，采用基于RANS的湍流模型建立了CFD风场3D模型。Endalenw等^［24,25］测量距离出风口15 cm处三维风速作为模型的速度入口，采用平面和曲面组合的几何结构模拟不同类型的喷雾机出风口，采用基于RANS的湍流模型构建喷雾风场。Duga等^［26］采用矩形截面模拟出风口，以出风口风速轮廓作为计算的初始条件，使用非定常雷诺平均方程（Unsteady Reynold-Averaged Navier-Stokes，URANS）和 $k - ε$ 模型求解计算域气流。Garcia-Ramos构建风机的内部结构，以进入风机的总体风量代替出风口风速作为模型的入口边界条件，使用ANSYS-CFX软件，基于Spalart-Allmaras湍流模型（用于求解粘性湍流运动单方程模型）建立风场^［27］。Hong等^［28］以圆柱体代替出风口，以测量出风口风速值作为模型的入口风速值，采用基于RANS的SST湍流模型建立风场。Badules等^［29］进行出风口内部几何结构建模，并研究了不同湍流模型下仿真结果差异，表明准确的出风口内部结构建模显著提高了仿真结果的准确性，选择合适的湍流模型对风场建模至关重要。

研究者针对不同类型喷雾机进行了出风口外部CFD风场建模方法研究，发现塔式果园风送喷雾机将出风口风力与药量由垂直方向变为水平方向，可提高果树冠层顶部的药液沉积。但塔式结构内部导流装置布置不合理会导致出风口风速分布不均匀且存在湍流的现象。为探究适用塔式果园风送喷雾机出风口风场建模方法，本研究采用不同的边界条件和湍流模型，建立了喷雾机空间风场对比模型，通过搭建空间风场立体测量试验台，进行了风场测量和模型验证，以获得一种不同边界条件设置方法，以及计算域尺寸和湍流模型组合的最佳果园喷雾机外部风场CFD建模方法。

2 CFD风场建模

2.1　喷雾机结构

试验用果园风送喷雾机型号为3WG-1200A，出风口最大高度204 cm。喷雾机风箱外形为塔式，风箱高度154 cm，下边沿距地面高度为50 cm，出风口宽度8.3 cm，出风口处安装10个导流板，将出风口分为11个格段，导流板厚度为1 cm，如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 试验用喷雾机

Fig.1 Sprayer for test

2.2　控制方程与建模方法

喷雾风场模型采用CFD商用软件Fluent16.0建立。建模时，假设空气是恒温、不可压缩的粘性牛顿流体，空气参数为经验值（Fluent的默认值）。控制方程包括连续方程、动量方程和湍流方程^［30］。密度恒定的定常流动的连续方程如公式（1）所示。

\partial u / \partial x + \partial v / \partial y + \partial w / \partial z = 0

（1）

其中，u、v、w是速度矢量u在x、y、z方向的分量，m/s。

N-S动量守恒方程用于描述流体运动，但其对计算机内存和性能要求非常高，本研究采用基于RANS动量方程的湍流模型模拟喷雾机风场，RANS动量方程组如公式（2）所示。

\{\begin{matrix} d i v (ρ u u) = d i v (μ g r a d u) - \frac{\partial p}{\partial x} + [- \frac{\partial (ρ \bar{{u^{'}}^{2}})}{\partial x} - \frac{\partial (ρ \bar{u^{'} v^{'}})}{\partial y} - \frac{\partial (ρ \bar{u^{'} w^{'}})}{\partial z}] \\ d i v (ρ v u) = d i v (μ g r a d v) - \frac{\partial p}{\partial y} + [- \frac{\partial (ρ \bar{u^{'} v^{'}})}{\partial x} - \frac{\partial (ρ \bar{v^{' 2}})}{\partial y} - \frac{\partial (ρ \bar{v^{'} w^{'}})}{\partial z}] \\ d i v (ρ w u) = d i v (μ g r a d w) - \frac{\partial p}{\partial z} + [- \frac{\partial (ρ \bar{u^{'} w^{'}})}{\partial x} - \frac{\partial (ρ \bar{v^{'} w^{'}})}{\partial y} - \frac{\partial (ρ \bar{w^{' 2}})}{\partial z}] \end{matrix}

（2）

其中，u，v，w为三维风速的平均值，m/s； $u', v', w'$ 为三维风速的脉动值，m/s；ρ为空气密度，kg/m³；p为气流压力，Pa； $S_{u}$ 、 $S_{v}$ 、 $S_{w}$ 为对应源项， $- ρ \bar{u' v'}$ 为雷诺应力项（存在6个不同的雷诺应力项）。

$k - ε$ 湍流模型和 $k - ω$ 湍流模型是风场建模中较常用的湍流模型。 $k - ε$ 湍流模型的输运方程见公式（3）和公式（4）^［30］。

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ k u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{u_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + G_{k} + G_{b} - ρ ε - Y_{M} + S_{k}

（3）

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ ε u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{u_{t}}{σ_{k}}) \frac{\partial ε}{\partial x_{j}}] + C_{1 ε} \frac{ε}{k} (G_{k} + C_{3 ε} G_{b}) - C_{2 ε} ρ \frac{ε^{2}}{k} + S_{ε}

（4）

其中， $μ$ 为流体的动力粘度， $N \cdot s / m^{2}$ 。 $G_{k}$ 是由于平均速度梯度引起的湍动能k的产生项； $G_{b}$ 由于浮力引起的湍动能k的产生项； $C_{1 ε}$ 、 $C_{2 ε}$ 和 $C_{3 ε}$ 为经验常数； $σ_{k}$ 和 $σ_{ε}$ 分别是与湍动能k（公式（4））和耗散率 $ε$ （公式（6））对应的Prandtl数； $S_{k}$ 和 $S_{ε}$ 是用户定义的源项。

k = \frac{\bar{u_{i}^{'} u_{i}^{'}}}{2} = \frac{1}{2} (\bar{u^{' 2}} + \bar{v^{' 2}} + \bar{w^{' 2}})

（5）

其中， $u_{i}$ 为时均速度的张量符号，m/s；k为湍动能，m²/s²。

μ_{t} = ρ C_{μ} \frac{k^{2}}{ε}

（6）

其中， $C_{μ}$ = 0.09，为经验常数；ρ为空气密度，默认为1.293 kg/m³； $μ_{t}$ 为湍动粘度， $N \cdot s / m^{2}$ 。

对于定常流动， $\partial () / \partial t = 0$ 。不可压流体由于浮力引起的 $G_{b}$ =0，可压湍流中脉动扩张的贡献 $Y_{M}$ = 0，且在方程设置中用户无自定义的源项，本研究中k和 $ε$ 的输运方程可简化为：

\frac{\partial (ρ k u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{u_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + G_{k} - ρ ε

（7）

\frac{\partial (ρ ε u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{u_{t}}{σ_{k}}) \frac{\partial ε}{\partial x_{j}}] + G_{1 ε} \frac{ε}{k} G_{k} - C_{2 ε} ρ \frac{ε^{2}}{k}

（8）

定常流动的k-ω湍流模型的湍动能k和比耗散率ω的输运方程如下所示^［31］。

\frac{\partial (ρ k u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [\frac{u_{t}}{σ_{k}} \cdot \frac{\partial k}{\partial x_{j}}] + P_{k} - Y_{k}

（9）

\frac{\partial (ρ w u_{i})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [\frac{u_{t}}{σ_{w}} \cdot \frac{\partial w}{\partial x_{j}}] + P_{w} - Y_{w}

（10）

w = ε / k

（11）

其中， $P_{k}$ 和 $P_{w}$ 代表了k和w的产生项； $Y_{k}$ 和 $Y_{w}$ 代表了湍流有效耗散项，其他参数与湍流模型中的参数定义规则一致。

Standard $k - ε$ 模型用于强旋流或带有弯曲壁面的流动时会出现失真。RNG $k - ε$ 可以更好地处理高应变路及流线弯曲程度较大的流动。Realizable $k - ε$ 模型可用于旋转均匀剪切流，管道内流动，边界层流动，以及带有分离的流动^［30］。BLS $k - ω$ 湍流模型在内部区域采用 $k - ω$ 湍流模型，在自由剪切流中采用 $k - ε$ 湍流模型，避免了 $k - ω$ 模型的强自由流敏感性。考虑了主湍流剪应力的传递效应，在BSL $k - ω$ 湍流模型的基础上对湍动粘度进行修正，即SST $k - ω$ 湍流模型^［32］。为探究不同湍流模型对计算结果的影响，本研究计算了Standard $k - ε$ 、RNG $k - ε$ 、Realizable $k - ε$ 、SST $k - ω$ 、BSL $k - ω$ 和Standard $k - ω$ 六种湍流模型的风场，并进行了模型对比。

采用半隐式算法（Semi-Implicit Method for Pressure-Linked Equations，SIMPLE）进行流场数值计算，其特点是在交错网格的基础上求解压力场，从而求解动量方程，需要提供初始的压力场和速度场，通过不断地迭代逐渐逼近压力和速度的真值。离散格式为二阶迎风，在二阶迎风格式中实际上只是对流项采用了二阶迎风格式，而扩散项仍采用中心差分格式。模型中，喷雾机出风口域四周和空气域的地面均设置为壁面，壁面粗糙高度 $k_{s} = 0$ ，壁面粗糙度常数 $C_{s} = 0.5$ ，湍流强度I = 5%，湍流粘度比 $μ_{t} / μ$ = 10。

2.3　边界条件

为对比不同边界条件对模拟结果的影响，CFD风场建模分别设置两种边界条件。边界条件一设置方法：喷雾机出风口被10个导流板分割为11部分，对被导流板分割的出风口内风速平均，生成图2（a）所示散点，并以笛卡尔坐标输入作为模型一的风速入口。

边界条件二如图2（b）所示，采集的风速数据使用C语言写成分段式风速方程，再采用用户定义函数（User-Defined Function，UDF）函数和笛卡尔坐标法导入三维风速曲线，作为模型二和模型三风速入口。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 喷雾机出风口风速

Fig. 2 Air velocity at outlet of the sprayer

如图3所示，计算域分为空气域和喷雾机出风口，通过交界面连接空气域和喷雾机域，风力能够穿过交界面进入空气域；四周设置为无滑移壁面（图3（a））。空气域（图3（b））入口除去交界面的其他位置设置为压力入口，压力为0；侧面及顶面设为压力出口，压力也为0，地面设置为无滑移壁面。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 喷雾机域和空气域条件设置

注： W为宽度width， L为长度long， H为高度height

Fig. 3 Computational domains consisting of atmosphere part and sprayer part

2.4　计算域与网格划分

如图4所示，共建立三种对比模型。模型一（图4（a））对应边界条件一，计算域10 m（W）×10 m （H）×25 m（L）（其中W为宽度width，L为长度long，H为高度height）；模型二（图4（b）），对应边界条件二，计算域10 m（W）×10 m（H）×25 m（L）。风场试验在大田进行，即开放的空间。为了研究计算域尺寸对仿真结果影响特性，建立了模型三（4（c）），对应边界条件二，计算域为5 m（W）×7 m（H）×10 m（L）。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 关于不同边界条件研究的三种CFD对比模型

Fig. 4 Three comparison CFD models for different boundary conditions

采用ICEM CFD 16.0建立几何结构并划分网格。空气域为正六面体结构网格，喷雾机域为四面体非结构网格，网格数量与尺寸参数如表1所示。模型一由于其几何结构，需要在划分网格的过程中依据不同的出风口进行分块划分，除此之外，三种模型网格划分方法相同。以模型二为例，在喷雾机行进方向设置250个网格；在竖直方向，喷雾机出风口网格尺寸为0.04 m，并从0.04 m逐渐增大至0.5 m；在下风方向，喷雾机出风口网格尺寸为0.04 m，并从0.04 m逐渐增大到0.5 m。最后，细化喷雾机出风口附近网格，加密后的喷雾机出风口附近网格几何三维尺寸为0.02 m（W）×0.02 m（H）×0.02 m （L）。

表1 不同边界条件研究三种模型网格数量与尺寸

Table 1 Grid parameters of three models

网格数量与参数设置

模型一

模型二

模型三

四面体网格数量

四面体网格参数/m

4201

0.02

4474

0.02

8242

0.01

正六面体网格数量

3,820,375

3,434,833

1,114,074

正六面体网格参数/m

0.02

新窗口打开| 下载CSV

3 试验材料和方法

3.1　空间风场立体测量试验台设计

针对喷雾机空间风场区域大且复杂难以快速采集获得的问题，搭建了空间风场立体测量试验台，如图5所示。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 空间风场立体测量试验台

Fig. 5 Three dimensional measurement test-bed for airflow field

试验台由移动装置步进电机控制器、电机驱动器、导轨、滑块、同步带、支架和地轮、风速仪安装支架、风速仪，以及上位机三部分构成。移动装置中步进电机控制器（KH-01型，浙江热龙自动化有限公司）通过驱动器（HB860D型，深圳市鸿宝达有限公司）精确控制步进电机（HB80D型，深圳市鸿宝达有限公司）带动滑块沿导轨平稳移动，滑块移动速度和行程精确可控。风速仪支架与滑块刚性连接，带动固定在其上面的风速仪随滑块精确移动。

风速仪为热敏式风速仪（型号GM8903，深圳市聚茂源科技有限公司），测量范围为0~30 m/s，分辨率为0.001 m/s，采样频率为1 Hz，与上位机通过USB连接，并用Hot Wire Anemometer软件实时存储并显示风速值。

3.2　田间试验方法

试验在北京小汤山国家精准农业示范基地开展，喷雾机拖拉机头朝向正西，在试验过程中测量行进方向右侧的风场。

（1）喷雾机出风口风速测量。如图6（a）所示，喷雾机出风口按照几何形状分为竖直面、曲面和斜面三部分。喷雾机出风口测量原理和过程见图6（b），首先沿风箱外边缘用马克笔以1 cm测量间距标定，风速仪贴着标定刻度，在x、y、z方向，于喷雾机出风口左、中、右三个位置，间隔1 cm测量，共测量160组数据。将同一高度左、中、右3个风速值取平均值，得到160组喷雾机出风口三维风速。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 喷雾机出风口风速测量

Fig. 6 Air velocity measuring test of outlet of sprayer

（2）空间风场立体测量。如图7所示，在喷雾下风方向（Y方向），采样点与风箱出风口竖直部分的水平距离（S）分别为0.3、0.5、1.0、1.5、2和3 m（下文提到的水平距离S含义均与此处相同）。在竖直方向（Z方向），取10个采样高度，在0.23 至2.03 m之间采样间隔为0.2 m。试验过程中，在滑块去程进行数据测量，滑块位移为3 m，移动速度为0.02 m/s；回程不进行数据测量，滑块移动速度为0.3 m/s。在每个采样高度获取150个风速值，共获取60组数据，可以覆盖整个喷雾空间风场。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 空间风场测量采样分布图

注：彩色线条代表采样位置

Fig. 7 Distribution of sampling locations of airflow field measurement

喷雾机、试验台和基准台相对位置如图8所示，激光尺固定于基准台中心。首先，打开激光尺（型号为LM40，优利德科技股份有限公司），调整风速仪探头，使激光束穿过探头孔中心（孔尺寸为0.6 cm×0.8 cm），确保每一次测量起始点一致。以基准面水平度为基准水平，使用水平仪（EM5416-150型，乐清市伊莱科电气有限公司）调整试验台，确保试验过程中试验台的俯仰角一致。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 喷雾机、试验台和基准台相对位置图

Fig. 8 Relative position map of sprayer， test-bed and base station

超声波风速风向传感器安装在喷雾机侧后方，与喷雾机前进方向夹角为152°，距离基准面中心为14.8 m，安装高度2 m。风场测量试验共持续2 d，试验过程中自然风平稳，在测量过程中自然风有少量时间超过3 m/s，但整体的测量环境符合国际上风速喷雾机作业标准^［33］，可认为田间试验风场测量数据可靠，作为风场模型验证的基础数据。

4 仿真结果与验证

4.1　空间风场立体测量试验台设计

如表2所示，对三种CFD模型采用六类湍流模型，使用计算机（Intel i5，RAM，16GB，Win10）计算喷雾风场。

表2 不同湍流模型计算的3种CFD模型特性

Table 2 Characteristics of three CFD models for different turbulence models

特性	模型一	模型二	模型三
湍流模型与连续性残差	SST $k - ω$ ，2.43×10^-2	SST $k - ω$ ，2.70×10^-3	SST $k - ω$ ，4.05×10^-3
	BSL $k - ω$ ，2.47×10^-2	BSL $k - ω$ ，9.71×10^-2	BSL $k - ω$ ，2.31×10^-3
	Standard $k - ω$ ，5.30×10^-2	Standard $k - ω$ ，k发散	Standard $k - ε$ ，4.54×10^-2
	Standard $k - ε$ ，4.09×10^-3	Standard $k - ε$ ，1.91×10^-2	Standard $k - ε$ ，1.38×10^-2
	RNG $k - ε$ ，3.30×10^-1	RNG $k - ε$ ，1.88×10^-1	RNG $k - ε$ ，9.12×10^-2
	Realizable $k - ε$ ，3.99×10^-3	Realizable $k - ε$ ，1.88×10^-2	Realizable $k - ε$ ，2.76×10^-2
运算结果是否可接受	RNG $k - ε$ 湍流模型，连续方程残差偏大，运算结果不可接受；其余湍流模型运算结果均可接受	Standard $k - ω$ 湍流模型，k发散，模型不收敛；RNG $k - ε$ 湍流模型，连续方程残差不低于10^-2；其余湍流模型运算可接受	对于六种湍流模型，运算结果均可接受

新窗口打开| 下载CSV

模型网格数量可能会对模型精度产生影响，在网格建模同时进行网格无关性验证。对于模型一和模型二的RNG $k - ε$ 模型，残差曲线一直上下正弦波动，难以收敛到10^-2以下，连续性残差偏大。模型二Standard $k - ε$ 湍流模型的湍动能k计算发散。在Badules等研究中，喷雾机产生气流的复杂性使得连续性方程的残差难于小于10^-4，当连续性方程的残差稳定在3×10^-2时，视为可接受的计算结果，计算中止^［34］。本研究中SST $k - ω$ 、BSL $k - ω$ 、Standard $k - ε$ 和Realizable $k - ε$ 四种湍流模型计算风场的连续性残差（运算中最大值）均可稳定在10^-2，计算结果可接受。

风场中心面速度最大，具有代表性，选取长为5 m、高为3 m的平面，提取模型的中心面速度云图。如图9为表2得出的四种运算可接受的湍流模型计算的喷雾机外部风场中心面云图。模型一与模型二，或者模型一与模型三的风场云图分布轮廓存在明显差异，可见边界条件设置方法是模型风场分布的主要影响因素。对比模型二和模型三，可以发现计算域尺对风场有微小的影响。对于同一几何模型（以模型二为例），当湍流模型不同时，风场分布也有微小的差异。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 三种模型中心面风场分布云图

Fig. 9 Distribution of airflow field on the center of three modes

因此，边界条件设置方法是影响风场分布的主要因素。运算可接受的湍流模型和计算域的尺寸是影响风场分布的次要因素。

4.2　模型验证

图10为水平距离S=1.0 m，距地高度H=1.03 m的实测数据曲线，风场风速横向分布呈现不完全左右对称的山峰状，横坐标是风场横向距离，风速仪移动速度为0.02 m/s，采样总长度为3 m。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 当S=1 m，H=1.03 m的实测风速曲线

Fig. 10 Measurement air velocity curve，at S=1 m，H=1.03 m

考虑到实际风场存在湍动，选取实测风速曲线峰值风速，如表3所示，用于CFD模型验证，当H=1.03 m，S=2 m时，测量位置风速为1.33 m/s，同一高度处S=3 m时，风速为2.61 m/s。一般来说，距离出风口较近位置风速应大于较远位置风速，可能由于自然风风速风向变化，使得实际风场分布与理论风场分布不同，类似的点还有H=0.63 m，S=1.5 m和H=2.03 m，S=1.5 m。

表 3 喷雾风场不同高度和不同水平距离时风场宽度方向的实测最大风速 (m/s)

Table 3 Peak measurement air velocity of sampling points with different height and horizontal distances

水平距离/m	距地高度/ m
水平距离/m	0.23	0.43	0.63	0.83	1.03	1.23	1.43	1.63	1.83	2.03
0.3	1.099	5.962	14.73	22.03	17.65	21.07	11.78	15.7	11.36	7.208
0.5	1.744	2.489	5.901	18.67	10.68	19.61	12.09	10.67	10.26	8.847
1.0	1.151	5.362	4.866	8.034	11.13	10.5	12.72	9.014	6.39	8.544
1.5	1.439	1.901	1.413	4.302	5.784	5.401	7.678	8.974	6.578	4.173
2.0	2.126	1.604	1.932	2.554	1.33	3.336	5.845	6.593	2.717	6.337
3.0	1.901	1.547	1.59	1.81	2.61	1.182	3.749	1.964	2.264	2.303

新窗口打开| 下载CSV

在模型的中心面采用CFD Post提取试验采样点位置的风速点云，结合表3数据做出散点图和趋势线图，如图11所示。其中模型二对应Standard $k - ε$ 湍流模型，实测风速与仿真风速散点关于拟合直线聚合性最佳，模型的决定系数R²为0.81。空间风场试验在大田进行，图11（g）标记左下角的红色小圈内，模拟风速值均小于仿真风速值，可能的原因是在喷雾风场影响微弱的区域，自然风起主要作用，使得该区域测量风速大于模拟风速。考虑到田间复杂的自然环境，风场模型的决定系数R²仍可达到0.81，表明该CFD风场模型能够准确预测风场。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 模型一、模型二、模型三中心面仿真风速值与实测风速对比值

Fig. 11 Comparison of downwind air velocities between the CFD simulation values and field measurement values at the center surface of Model 1， Model 2 and Model 3

模型一几何建模过程复杂，且与模型二相比，其模型决定系数R²偏低，说明采用一个完整的出风口，通过UDF导入风速曲线是简单准确的建模方法。模型二与模型三相比较，模型二决定系数略高于模型三，可见大计算域的模拟风场更加符合实际风场，但计算域尺寸增加，模型的运算时长增加，且对计算机性能要求也大幅提高。大尺寸计算域的优势是便于后期进一步模拟动态的、整个果园内的风场。在风场建模中建议综合考虑运算能力和需要建模的果树、果园规模，设计风场计算域。

以模型二为例（图11（e）~11（h）），不同湍流模型对应的风场计算结果，其散点分布与模型决定系数存在差异。表明在CFD风场建模过程中，选择适合的湍流模型可以提高模型的准确性。

综上，以喷雾机出风口三维风速分段折线作为边界条件，选择较大的计算域和Standard $k - ε$ 湍流模型建模（图11（g）），模拟风场与实际风场最接近，是最佳的风场建模方法。

5 讨论

由喷雾机出风口风速云图（图9）可见，风速变化复杂，且存在湍流。如图2（b），喷雾机出风口实测风速曲线，在不到20 cm的高度差内风速变化超过20 m/s。Z向风速，在大于1 m（距地高度）时，其大小与Y向风速（垂直出风口方向）相近，所以Z向风速对风场的影响将不可忽略。X向风速关于轴线左右摆动，这是造成实际喷射气流左右摆动的主要原因。若采用一维风速作为模型的风速入口，其结果将不能反应实际中复杂湍流的喷雾机，会给模型带来严重的误差。

Hong等^［28］使用4个风速值作为速度入口边界条件却得到了较好的仿真结果，分析原因，其试验用喷雾机内部是由软管将气流分离，每一个出风口风速较稳定，而本研究使用的喷雾机出风口风速分布极为复杂，且存在湍流。国内主流风送喷雾机大都与本研究喷雾机相似，对于该类型的喷雾机，无法照搬Hong等的建模方法，仅考虑出风口一维风速获得的仿真结果不符合实际，应全面测量出风口三维风速，以相应的三维风速曲线作为模型的风速入口边界条件。

测量喷雾机出风口三维风速，借助UDF生成风速曲线，通过笛卡尔坐标导入三维风速，进行CFD风场模型建立，模拟结果显示，此边界条件设置方法不仅降低了建模复杂性，而且具有较高的准确性与工程价值。边界条件是决定模型准确性的关键，采用UDF出风口风速折线式边界条件设置方法，输入三维风速进行风场建模是一种简单准确的CFD建模方法。

6 结论

本研究通过研究风速入口边界条件设置方法、计算域和湍流模型的选取对CFD风场模型的影响，确定最优的CFD风场模型建模方法，并进行了田间试验验证，结论如下：

（1）提出了果园喷雾机风场CFD模拟UDF分段式风速入口边界条件三维风速设置方法。以Standard $k - ε$ 湍流模型为例，采用UDF分段式风速入口边界条件（模型二）相比采用离散点式三维风速入口边界条件（模型一），仿真风场准确性提高了5%，且采用UDF分段式风速入口边界条件设置方法，几何建模简单，模型运算经济。

（2）采用UDF分段式三维风速入口边界条件设置方法和离散点式三维风速入口边界条件设置方法，10 m （W）×10 m （H）×25 m （L）和5 m （W）×7 m （H）×10 m（L），两种计算域尺寸，建立了CFD对比模型。验证试验结果表明，采用Standard $k - ε$ 湍流模型，选择较大的计算域，仿真准确性最佳，模型决定系数R²为0.81。

（3）设计了风送喷雾机空间风场立体测量系统，以2 cm的间隔采集风速，获取10个高度、6个平面的三维立体风速，具有定位精确、测量快速且数据完整的特征，为风场仿真验证试验提供了良好的测量平台。