Distribution regularity of downwash airflow under rotors of agricultural UAV for plant protection

In the plant protection spray operation of UAVs, the process of droplet from formation to sedimentation target is affected by airflow, easy to form uneven deposition. Accurately description of rotor downwash flow field, clarification of velocity vector distribution at different heights of the UAV rotor flow field, simulation of the flow field with high precision, which are the prerequisites for accurately analyzing the droplet deposition distribution in rotor downwash flow field. Based on CFD method, the detail of rotor flow field was numerically calculated. Taking LTH-100 single-rotor agricultural UAV as the research object, the three-dimensional solid model of UAV was established, the Reynolds average N-S equation was used as the control equation and the RNG κ-ε as the turbulence model to simulate the flow field of UAV in hover and lateral wind conditions, the wind velocity distribution at different altitudes of rotor downwash flow field was studied. The simulation results of the hover state showed that: In the flow field, the peak velocity appears in a circular distribution below the distal axis of the rotor. With the decrease of height, the peak velocity distribution area showed a tendency to expand gradually after small shrinkage; When the distance from the rotor was not more than 1.5 m, the downwash flow field presented an axisymmetric distribution based on the rotor axis, and the variation rate of velocity in the peak velocity was basically the same, turbulence in downwash flow field made the flow field more complex when the distance from rotor was larger than 2.0 m. On this basis, the optimal flight altitude of UAV is 1.5 m. Wind velocity test of the flow field was carried out on a rotor test bench, wind velocities at four altitudes of 0.5 m, 1.0 m, 1.5 m and 2.0 m were measured to verify the coincidence between the simulated and measured values. The test results showed that: the relative error between the measured and simulated values at four measurement heights were between 0.382-0.524, and the overall average relative errors was 0.430, which verified the confidence level of simulated values for measured values. When the lateral wind velocity was 3 m/s, 4 m/s and 5 m/s, the simulation results showed that: The distribution trend of airflow velocity at the same altitude in lateral-wind flow field with different wind speeds was similar; When the lateral wind speed was 5 m/s, the coupling field formed by the lateral wind and rotor airflow cannot reach the height of 2 m below the rotor. The results of this study can provide more accurate environmental conditions for theoretical analysis of droplet deposition regularity in the flow field, and also provide methodological guidance for the related research on rotor flow field of multi-rotor UAV.


Introduction 
In recent years, agricultural UAV for plant protection has been widely concerned and valued by scholars because of its high operational efficiency and high flexibility [1][2][3][4] . When the UAV carried out plant protection operation, the rotor downwash airflow could increase the initial velocity of the pesticide droplet, which is helpful to improve the penetration of the droplet to the crop canopy, increase the sedimentary density on the back of the leaf blade, and improve the pesticide utilization rate [5][6][7] . There is great ability to coping with sudden pests and diseases by UAV, it can effectively solve the problem that difficult operation of ground machinery in the middle and late stage of crops.
The operation mode of UAV for plant protection usually adopts low capacity or ultra-low capacity spraying, deposition of droplet in target area is one of important indexes to evaluate the quality of spraying [8] .
Compared with ground application machinery, UAV spraying height is higher, droplets are easily affected by operating conditions and airflow. It is easy to produce pesticide droplet drift in the course of operation, causing phenomenon of repetition spray and leakage spray. Compared with fixed-wing aircraft, because of rotor rotation, the inflow velocity of blade varies with the azimuth angle, which leads to the uncertainty of droplet drift and deposition during spraying operation [9,10] . It can be concluded that the rotor downwash flow field characteristics are important factors affecting the drift and deposition of pesticide droplets. Comprehensive analysis of current research results [11][12][13][14][15][16] : UAV operation height, rotor downwash airflow and ambient airflow are the main factors causing unstable operating width, re-spray, leakage and drift. Therefore, it is necessary to explore the influence of rotor and ambient airflow on droplet deposition at different operating altitudes of UAV on the basis of existing research.
UAVs are divided into single rotor and multi-rotor. The rotor number, airfoil shape and size of each UAV are different, and the distribution of downwash airflow of its rotor is also greatly different [17] . At present, the research on UAV spraying was mostly focused on the measurement of pesticide distribution through specific experiments, and the influence of airflow on spraying was less explained from the perspective of principle. In this study, the difference between the simulated airflow field and the real airflow field was compared and analyzed from the principle point of view through the numerical simulation method to explore the accuracy of the rotor airflow, so as to provide a basis for further research on the influence of the rotor flow field on pesticide spraying. In order to ensure the accuracy of simulation and test as much as possible, the single-rotor UAV was selected as the research platform in this study. Compared with the interaction among the airflow fields of the multi rotor UAV, the downwash airflow field of the single rotor UAV had no influence by other airflow, with smaller error and higher accuracy in terms of simulation accuracy and actual test data acquisition [18][19][20][21][22][23] . Meanwhile, accurate simulation of rotor flow airfield of single-rotor UAV could provide method reference and basic data support for further study of rotor airflow field of multi-rotor UAV.
According to the aerodynamic principle of helicopter, rotor flow field is a kind of vortex formed by the diversion of incoming flow through the rotor surface when the rotor blade moves around the rotor axis. Airflow is accelerated downward by the action of the rotor, and the increase of airflow velocity is the induction velocity of that place. Therefore, the physical quantity that can best reflect the characteristics of the rotor flow field are vortex and induction velocity. The size and properties of the vortex are mainly designed for optimizing the rotor, however, in the field of aviation plant protection, the research on the flow field is focused on the air flow below the rotor [16] . Therefore, the induction velocity is taken as a parameter to consider the characteristics of the flow field, the flow field is described by analyzing the velocity distribution of flow field airflow.
Based on computational fluid dynamics (CFD) method, the rotor flow field of LTH-100 single-rotor UAV in hover and forward flight was numerically simulated. The velocity vector distribution trend at different altitudes in rotor downwash flow field was analyzed, the relationship model between the spatial position and velocity of the flow field was established, which provided theoretical support for further research on the deposition process of pesticide droplets in the rotor flow field.

Physical model establishment
The high speed rotation of rotor blades of single rotor UAV made the rotor flow field more complicated, it was necessary to make assumptions about the calculation model to improve the efficiency and accuracy of simulation calculation: (1) Assumed that air was incompressible, regardless the changes of air viscosity and temperature, and regardless the influence of blade selection on airflow [24,25] ; (2) Assumed that the top of the crop canopy was flat, ignored the different height of the plant, ignored the influence of the crop on the rotor flow field; (3) Ignored the interference of UAV tail to compute domain flow field.
The LTH-100 single-rotor UAV used in this study was developed and manufactured by Anyi Science and Technology Company (Liaoning). As shown in Figure 1, its maximum load was 25 kg, payload 20 kg, spray amplitude was 4-6 m, rotor diameter 2.1 m. The outer dimension of UAV is 2.05 m× 0.7 m×0.65 m. The forward direction of UAV was defined as X-axis positive direction, the forward direction of fuselage from right to left was defined as Y-axis positive direction, and the upward direction was Z-axis positive direction. The three-dimensional physical model of LTH-100 single-rotor UAV was established by using SolidWorks software. As shown in Figure 2, the rotor blade airfoil was NACA2415, the rotor radius was 1.05 m, the chord length was 0.07 m, the zero-lift angle of attack was -2°, the number of rotor blades was 2, and the speed was 1300 r/min.

Hover flow field establishment
In order to analyze the velocity distribution of the downwash flow field generated by high-speed rotation of UAV rotor, numerical simulation analysis was carried out with CFX (ANSYS Inc.) software in computational fluid dynamics (CFD) method.

Model computational domain
The calculation domain of UAV flow field was established by CFX gridding software, the flow field was divided into two parts, which were arranged in cylindrical shape respectively [20] , as shown in Figure 3a, one was the rotating domain surrounding UAV rotor; the other was the static domain outside the fuselage, as shown in Figure 3b [26]  In order to ensure the accuracy of the calculation results and reduced the simulation distortion caused by too small computational domain, it was necessary to ensure that the outer bounding static domain was large enough and the diameter D of the outer region of the computational domain was set to 10R [27] (R is rotor radius), the static domain radius was set to 11 m. The height of the rotor from the ground was 3.5 m, and the height of the external stationary region was 15 m. The computational domain of the flow field was 11 m×15 m (D 1 ×H 1 ) in the external stationary region and 2.2 m×0.08 m (D 2 ×H 2 ) in the rotor rotation region.

Control equation
The Renault average N-S equation was used as control equation: where, W is unknown conserved variable vector; F c is convective flux; F v is viscous flux; Ω is an arbitrary element in the flow field; For various Newtonian fluids, according to Stokes viscous coefficient hypothesis, velocity deviation can be calculated to obtain stress components: where, τ xx is the shear stress in xx direction, N; V  is divergence of vector V; μ=μ 1 +μ t is the viscosity coefficient; μ 1 , μ t are respectively the laminar and turbulent viscous coefficients.

Grid partitioning
The unstructured grids were used to divide the outflow field of rotor and fuselage. The contact surface between the rotation field and the stationary domain was set as the interface, encrypt grids at the computational domain boundary [28,29] , the total number of meshes are 6.3 million, and the computational domain grid is shown in Figure 4.

Turbulence model
RNG k-ε model was used as the turbulence model, rotation effect was added to the model, higher accuracy could be obtained in the calculation of rotating flow [30,31] .
Its continuity equation: The momentum equation: where, ρ is fluid density, kg/m 3 ; B is volume force synthesis; μ eff is effective viscosity, Pa· s; p′ is corrected pressure, Pa. Its expression was as follows: where, μ t is turbulent viscosity, Pa: The turbulent kinetic energy equation: where, products of viscous force and buoyancy; η is viscosity coefficient of fluid; k is turbulent kinetic energy, J; ε is turbulent dissipation rate; U is velocity vector; C μRNG is constants of RNG k-ε model, its value is 0.085; β RNG is constants of RNG k-ε model, its value is 0.012; σεRNG is constants of RNG k-ε model, its value is 0.7179. The SIMPLE algorithm was used to couple of pressure and velocity, the convection phase and diffusion phase in the control equation were separated by second order upwind format, and the convergence residuals was set to 0.0001.

Coupled flow field establishment
The deposition distribution of droplets in plant protection UAV was affected not only by the downwash flow field of the rotor, but also by the disturbance of the lateral wind [32] . It is necessary to accurately describe the coupling effect between the lateral wind and the rotor flow field to study the movement regularity of droplets in the flow field. Therefore, numerical simulation of UAV hover flow field under lateral wind conditions was carried out.

Computational domain and grid
The computational domain of the flow field was divided into two parts as shown in Figure 5. The computational domain included the rotational domain following the rotor rotation and the flow domain outside the rotational domain including the fuselage. The rotational domain was cylindrical in shape, and the rotating region is 2.2 m×0.08 m (diameter× height) in size. The flow domain set to cuboid structure, 12 m×12 m×8 m in size. Set the flow domain inlet as the left plane of the UAV forward direction, for the speed inlet, to simulate the effect of cross wind on the flow field, the range of velocity was set to 3-5 m/s. Unstructured grids were used to partite computational domain, local grids encryption was performed on the interface between two computational domains and the interface between the flow domain and the fuselage. Rotor flow field space contains a total number of meshes of 6.8 million, mesh unit mass of 0.822.

Turbulence model
The SST k-ω model was used as the turbulent model for numerical calculation, the most significant difference was that the propagation of turbulent shear stress was considered in turbulent viscosity compared with the standard k-ω model.

Test contents
Taking the rotor axis of UAV as the center and the heading direction as the starting direction, a direction was set clockwise every 90°, from D 1 to D 4 , totaling four measuring directions [33] , as shown in Figure 6. In each measurement direction, starting from the center of the rotor, a marker point was made per 30 cm, with a total of 10 measuring points from 30 cm to 300 cm. UAV  Figure 7 shows the distribution of velocity vectors at different altitudes from 1.0 to 2.5 m. Comparing the streamline diagrams at four heights, it can be concluded that: (1) The annular high velocity region enlarges its diameter and decreases its velocity with the increase of height; (2) With the increase of the distance from the spinning rotor axis, the trend of velocity distribution increases first and then sharply decreases outward. (3) With the increase of the distance from rotor axis, the distribution of velocity streamline becomes more chaotic, and the range of flow field by the rotor is larger.

Overall velocity distribution of hover flow field
In order to obtain the flow field characteristics of rotor downwash flow field of single rotor plant protection UAV, the hovering state of UAV was simulated and analyzed. The simulated cloud images of UAV in hover are selected for analysis. As shown in Figure 8, the X-Y plane velocity clouds at 1.0 m, 2.0 m and 3.0 m below the rotor were selected respectively. Figure 8 shows that: (1) the maximum wind speed position is located below the far axis end of the rotor, and the total velocity of the flow field shows the trend to increase first from the near Axis end to the far axis end and then reduce; (2) speed gradually reduced from the upper plane to lower plane, the trend of distribution range gradually increased; (3) Because of the blocking effect of the fuselage, the velocity distribution at the lower position of the fuselage is disordered, because of the close distance to the ground at 3.0 m, the velocity trend distribution is not obvious due to the effect of ground effect. Open Access at https://www.ijabe.org Vol. 14 No. In order to clarify the velocity distribution of the rotor flow field at different altitudes, the X-Z plane at Y=0 and the Y-Z plane at X=0 were selected as reference planes, and six speed detection lines L were set on the two reference planes respectively. The detection lines were horizontal with the ground and were 0.5 m, 1.0 m, 1.5 m, 2.0 m, 2.5 m and 3.0 m below the UAV rotor, respectively, as shown in Figure 9. Figure 10 shows that the total velocity distribution of each height on the X-Z plane: (1) The peak velocity appears directly below the radius of the 9/10 rotor; (2) On the tail side of UAV, the peak velocity decreases with height and the position expands outward in turn. It shows that the head has a larger volume and irregular shape than the tail, which has a greater impact on the airflow.  Figure 11 shows the velocity distribution in the X, Y and Z directions at six detection lines altitudes on the X-Z plane. The change rate of the upward velocity of X positive direction is significantly weaker than that of X negative direction, and the air disturbance effect of UAV fuselage on X direction is obviously stronger than that of the tail, which is consistent with the results obtained above. On the inside of the peak velocity, the velocity mutation is small, while at the lateral peak speed, the velocity mutation is large, which indicates that the downward airflow is mainly concentrated directly below the rotor. Figure 12 shows the total velocity distribution at the position of 6 detection lines on the Y-Z plane. As is shown from the graph, the velocity distribution at 0.5 m, 1.0 m and 1.5 m below the rotor is basically symmetrically distributed with Y=0, and the peak velocity appears just below the far-axis end of the rotor blade. Figure 13 shows the velocity distribution in the X, Y and Z directions at six detection lines altitudes on the Y-Z plane. When the height is not more than 1.5 m, the velocity on both sides of the rotor axis is symmetrically distributed, the peak velocity appears below the rotor far axis end, and decreases with the increase of height. When the height is not more than 1.5 m, the velocity is symmetrically distributed with the rotor axis, the velocity increases with the increase of height. When the height is greater than 2.0 m, the velocity distribution is irregular.

Variance analysis of measurement data
Based on the analysis of the data collected in the wind speed measurement experiments of 16 groups in 4 directions, the wind speed scalar values of X, Y and Z directions of each velocity point are obtained, the average value is sampled 15 times per measurement position, 5 s per sampling time is recorded after the value is stable, and the data are processed according to Equation (14). The processing results are shown in Figure 14. As shown in Figure 14, in the direction of D1 to D4, the standard deviation of the measured points in X, Y and Z directions are all less than 0.8 m, the fluctuation of the collected data relative to the average value is small, which indicates that the error of the collected data is low. The maximum standard deviation appears in both X direction and Y direction at 0.5 m measurement height, which indicates that the data collected in X direction and Y direction fluctuate greatly. The minimum standard deviation appears at the 2.0 m height, which indicates that the data collected at the height of 2.0 m have little fluctuation from the average value, and the measured value is of good accuracy.

Wind speed distribution at different heights of flow field
When the flight altitude of UAV is 3.5 m, the measured altitudes are 0.5 m, 1.0 m, 1.5 m and 2.0 m, respectively. Airflow velocity in X, Y and Z direction of each collected data in four directions from D1 to D4 is shown in Figures 15-17.
As shown in Figure 15, the airflow velocity of each collected data in the X direction shows a decreasing trend with the increase of the measurement height, i.e., the airflow velocity in the X direction is smaller as it is closer to the UAV rotor plane. As a whole, the X-direction partial velocity increases first and then decreases. The maximum velocity appears below the far-axis end of the rotor, which is in good agreement with the position of the peak velocity of the simulated flow field.
In Y direction as shown in Figure 16 that with the increase of the height of measurement, the integral velocity of each acquisition data in Y direction decreases, and the peak velocity decreases with the increase of height. In Z direction as shown in Figure 17, the peak velocity of each height is mostly concentrated around the radial 100 mm of the UAV rotor axis, which is very close to the lead hammer position at the far axis end of the rotor. The larger position of wind speed is the horizontal distance between 50 mm and 150 mm from the rotor axis of UAV, which shows that the peak distribution of vertical downward air flow in the rotor flow field is in good agreement with the position obtained by simulation.

Error analysis of measurement value and simulated value
In order to judge the uniformity of the measured and numerical simulations intuitively, the data measured by the actual test were compared with the simulated data of the corresponding position. Because some of the conditional factors set in the simulation process were ideal, the relative error equation of the simulated data to the test data was as follows, taking the test measured value as the standard value and the relative error as the evaluation index.
|| cm where, q c is the test measurements, m/s; q m is simulation value, m/s. The actual measured and simulated values of the airflow velocity in the X, Y, and Z directions of the D1-D4 direction are measured by the relative error statistical table calculated by Equation (15), which is limited to the length and lists the simulated values and test measurements of each measuring point under the measurement height 0.5 m, as shown in Table 1.
From Table 1, it can be concluded that the maximum relative error between the actual measured value and the simulated value is 3.761 at 0.5 m measurement height, and the minimum value is 0.004. There are 7 locations with relative errors exceeding 1, which account for a small proportion of 120 locations. In the four directions of D1-D4, the average relative error of wind speed in X direction is 0.382, the average relative error in Y direction is 0.524, and 0.385 in Z direction, and the overall average of relative error is 0.430.
For other measuring heights, the relative error of wind speed measured by measuring points in X, Y and Z directions are all less than 1, the overall average relative error at each height increases slightly with the increase of measuring heights. The overall average value of relative error at each height is as follows: 0.430 at 0.5 m, 0.499 at 1 m, 0.535 at 1.5 m and 0.569 at 2 m, which indicates that the confidence level of the simulated data relative to the actual measurement data is good.

Analysis of simulation results of coupled flow field
The X-Z plane at Y=0 and the Y-Z plane at X=0 were selected as reference planes. Four speed measuring lines were set on each reference plane.
Each speed measurement line was horizontal, located below the UAV rotor and the distance were 0.5 m, 1.0 m, 1.5 m and 2.0 m. The X, Y and Z direction velocity values were taken for each speed measuring line, as shown in Figure 18. Figures 19-21 show the airflow velocity distribution of 4 measuring lines in X, Y and Z directions on X-Z plane under different lateral wind conditions. As shown in Figure 19, on the same measuring height, the velocity extremum decreases with the increase of lateral wind speed. The position of extreme velocity at 0.5 m below the rotor is similar with other different lateral winds condition. In Figure 20, with the increase of lateral wind speed, the fluctuation of speed at the same measuring height decreases, and the velocity extremum decrease. In Figure 21, with the increase of lateral wind speed, the extreme velocity of each height in negative direction of X decreases, and the extreme velocity in positive direction of X increases. Figures 22-24 show the airflow velocity distribution of 4 measuring lines in X, Y and Z directions on Y-Z plane under different lateral wind conditions. Figure 22 shows that, with the increase of cross wind speed, the air velocity at each height decreases as a whole, and the extreme velocity in the negative direction of Y shifts to the left side of the figure. As shown in Figure

Discussion
The above studies show that the peak velocity in the flow field appears in a circular distribution below the radius of 9/10 rotor, with the decrease of height, the velocity peak distribution area presents the trend of gradual expansion after a small shrinkage. Velocity is relatively stable in flow field below the fuselage of the near-axis end of the rotor, the airflow is relatively stable in the distance no more than 2 m below the fuselage. This reveals the distribution trend of airflow velocity and the position of peak velocity in the wash flow field under rotor. Xu et al. [34] using SST k-ω turbulence model to calculate the near-ground flow field of UAV, the results showed that the velocity of the rotor increased first and then decreased, and the maximum velocity appeared near y/R=0.8. Xue et al. [33] , Yang et al. [36] , Zheng et al. [37] , Wu et al. [38] and Zhao et al. [39] simulated the downwash airflow distribution in the helicopter flow field, the general trend of wind speed change in the same radial direction was determined as follows: wind speed first changed from small to large, then from large to small, the maximum wind speed was below the 3/4 rotor radius. Compared with the previous research results, the conclusions of each study are consistent in the trend of velocity distribution of the downwash field under the rotor, while the peak velocity position is slightly different. The difference of the result may be due to the different selection of rotor parameters and turbulence model, the rotor radius and rotational speed had a direct influence on the eddy current of rotor wingtip, an important characteristic of the eddy current of UAV rotor wingtip is shrinkage distortion [40,41] . When the wingtip vortex left the rotor, it shrinked rapidly after being induced by the wake, and the peak velocity appeared near the wingtip vortex. Turbulence model was the embodiment of mathematical method in fluid mechanics, simulation accuracies of different turbulence models in different research conditions would also be different.
At present, there is no relevant standard to judge the accuracy of the simulation results. The relative error method was used to judge the coincidence between the simulation results and the actual measured values by scholar. Zhang et al. [42] and Xiahou et al. [43] using relative error method to measure the confidence of simulation value and the test value, the maximum overall relative error did not exceed 0.70, which can verify the credibility of simulation data. Existing research shows that the relative error method can validate the simulation credibility effectively.
The influence of side wind convection field is very complicated, Xu et al. [44] and Shi et al. [45] studied the flow field of UAV in the upwind, and analyzed the trend of flow field distribution. They describe the flow field mainly through the vertical downward airflow in the flow field. Compared with other people's research, this study describes the analysis mode of flow field characteristics in the direction of X, Y and Z axes in the three-dimensional coordinate system of flow field, which can describe the convection field more accurately, and can also provide a more accurate calculation model of airflow velocity for the later study of the sedimentary regularity and drift motion of droplets in the flow field.

Conclusions
(1) Based on CFD method, the rotor downwash flow field of the single-rotor UAV was studied, the downwash velocity along the rotor axis increased first and then decreased, airflow speed reached its peak near the rotor end. The peak velocity of downwash airflow was located below far-axis end of the rotor at different height.
(2) When the distance from the rotor was less than 1.5 m, the downwash flow field presented an axisymmetric distribution based on the rotor axis. The variation rate of the peak velocity was basically the same. The turbulence in the downwash flow field made the flow field more complex when the distance from the rotor was more than 2.0 m. The optimal operating height of the UAV was determined to be 1.5 m.
(3) Comparison and fitting analysis of wind field test results with simulation results, the average relative errors of wind speed in X, Y, Z direction were 0.382, 0.524, and 0.385 respectively, and the overall average of relative error is 0.430, less than 0.75, which indicates that the confidence level of simulation data relative to actual measurement data is good.
(4) The influence of lateral wind on the trend of velocity distribution was similar in most cases, with the increase of lateral wind. The velocity distribution of flow field was shifted to the downwind direction. When the lateral wind speed was 5 m/s, the height of 2.0 m below the rotor was almost unaffected, the rotor flow field cannot reach this planar position under this lateral wind condition.