The Airﬂow Field Characteristics of UAV Flight in a Greenhouse

: The downwash airﬂow ﬁeld of UAVs is insufﬁcient under the dual inﬂuence of greenhouse structure and crop occlusion, and the distribution characteristics of the ﬂight ﬂow ﬁeld of UAVs in greenhouses are unclear. In order to promote the application of UAVs in greenhouses, the ﬂow ﬁeld characteristics of UAVs in a greenhouse were studied herein. In a greenhouse containing tomato plants, a porous media model was used to simulate the obstacle effect of crops on the airﬂow. The multi-reference system model method was selected to solve the ﬂow ﬁeld of the UAV. Studies have shown that the airﬂow ﬁeld generated by UAV ﬂight in a greenhouse is mainly affected by the greenhouse structure. With the increase in UAV ﬂight height, the ground effect of the downwash ﬂow ﬁeld weakened, and the ﬂow ﬁeld spread downward and around. The area affected by the ﬂow ﬁeld of the crops became larger, while the development of the crop convection ﬁeld was less affected. The simulation was veriﬁed by experiments, and linear regression analysis was carried out between the experimental value and the simulation value. The experimental results were found to be in good agreement with the simulation results.


Introduction
Unmanned aerial vehicles (UAVs) offer a highly efficient and nondestructive method for spraying pesticides on crops and have been widely used in the agricultural industry, representing advanced agricultural technology [1][2][3]. Since agricultural UAVs are usually constructed of light materials with low stiffness, various airflows have a great impact on the stability, safety, and pesticide application of UAV flight [4][5][6], ultimately increasing the operation error. In order to improve the accuracy of UAV operation, it is necessary to study the characteristics of the UAV flow field, especially the downwash flow field.
In recent years, studies have shown that the airflow field of UAVs is generated downward from the rotor [7], as the maximum velocity field is maintained below the rotor [8], and its strength can be changed by adjusting the flight altitude [9]. Ni et al. [10] simulated the airflow field of different UAV flight altitudes using computational fluid dynamics (CFD) simulation software, which confirmed that the increase in UAV flight altitude will increase the area of the airflow field, but the velocity of the airflow field will decrease. At the same height, the farther away from the center of the flow field, the smaller the airflow velocity and velocity gradient. This is the conclusion when the UAV is still. Zhang et al. [11] successfully established a three-dimensional airflow field model of a six-rotor plant protection UAV using LBM (Lattice Boltzmann Method, LBM) and provided the UAV flight speed. The experiment showed that when the flight altitude was 3 m, the flight speed increased, and the airflow tilted obviously backward. At 1 m/s, due to the small velocity, the airflow diffused on the ground to produce a vortex ring and further moved backward. When the flight speed was above 2 m/s, the vortex ring disappeared, producing a Λ-shape transverse separation flow. When the velocity reached 4 m/s, the wake separated from the ground and formed a horseshoe vortex, which indicated that the velocity was no longer suitable for agricultural application. Li et al. [12] collected, analyzed, and processed experimental data from drone pollinators in the field and constructed a one-way two-dimensional wind field model of a drone rotor at the rice canopy using Fourier fitting, Gaussian fitting, and polynomial fitting. It was found that an important "steep wall" effect was observed when the wind speed under the drone rotor reached its maximum value, and the increased rate of forward wind speed was significantly higher than that of backward reduction. The whole wind field "steep wall" was symmetrical along the UAV flight direction. Wang et al. [13] analyzed from the perspective of the number of UAV rotors and measured the wind speed of rotor airflow perpendicular to the flight direction X, perpendicular to the ground direction Y, and parallel to the flight direction Z using a wind field sensor network measurement system under three types of UAV flight operation conditions. The experimental results showed that the wind power of a single rotor and four-rotor UAV in the downwash flow field perpendicular to the flight direction was the largest, while the wind power of the eight-rotor UAV in the vertical direction was the largest, which demonstrated the reason for the best application effect of an eight-rotor UAV. Guo et al. [14] combined compressible Reynolds-averaged Navier-Stokes (RANS) equations, the shear stress transport (SST) kω turbulence model, sliding grid technology and the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm to establish the computational fluid dynamics model of the downwash flow of a four-rotor agricultural UAV during hovering. The downwash flow field of the UAV was divided into four regions from time dimension analysis. It was found that the z-direction velocity of the downwash flow field differed with height, and the speed of stabilization also differed but stabilized in about 4 s. After the downwash flow field was stable, with the decrease in height, the influence of the rotor was weakened, and the airflow velocity right below the rotor was also reduced. However, due to atmospheric compression and airflow induction, the downwash velocity in the other three regions increased with decreasing height. In addition, at the same time and at the same height, the downwash velocity right below the rotor was greater than the other three regions. In another study, Li et al. [15] used an S-shaped pitot array and installed wind speed sensors to obtain wind field data. Vertical data, such as equivalent area, wind field width, and attenuation rate, were analyzed. It was found that when the plant distribution and density characteristics of rice were consistent, the lower the height, the smaller the wind field and the greater the decrease. That is, the dense leaves in the lower part significantly hindered the airflow, and the permeability of the wind field decreased with the decrease in height. Li et al. also analyzed the abnormal data in the range of 6.5 m/s (Pmedium) and 3.5 m/s (Pbottom) and found that this was caused by the vortex formed by the interaction between the wind field and the rice canopy. Finally, the ideal cone model was established based on the attenuation rate of the equivalent area of the rotor area and the equal height. Shi et al. [16] studied the wind-induced response of rice under the action of a multi-rotor UAV downwash flow field by using turbulence intensity and maximum wind speed as parameters. It was found that the downwash flow field was less affected by the external environment at low flight altitude, the pulsation characteristics were not obvious, and the turbulence intensity was not high. When the flight altitude was greater than 3 m, the turbulence intensity increased significantly with the increase in flight altitude. It is considered that the friction between the downwash flow field and the external low-speed air produces many different sizes of eddy currents during long-distance transmission. The eddy current system leads to a decrease in the downwash flow field velocity, an irregular change in the flow direction, and an increase in the turbulence intensity. With the increase of UAV flight height, the turbulence intensity at the rice canopy peaked and then gradually decreased to the final value consistent with the turbulence intensity of the ambient wind, but it was still much larger than that of the ambient wind conditions.
In summary, the research on UAV flow fields has mainly focused on the flow field under different flight conditions in outdoor environments and the blocking effect of crop convection field. There is still a lack of research on the characteristics of the UAV airflow field within a greenhouse structure. Therefore, in this study, a greenhouse planted with tomato was used as the research environment, and the CFD method was used to obtain the flow field simulation results of the UAV under different flight states in the greenhouse. The flow field distribution law of the UAV in the greenhouse was analyzed, and the accuracy of the simulation results was verified experimentally. These findings provide a theoretical reference for the application of UAVs in greenhouse plant protection.

Establishment of the Simulation Model
The three-dimensional model of the UAV was established based on the E360-S2 UAV of CNOOC. The main parameters are shown in Table 1. The rotor 8045 was obtained by an EXAscan handheld three-dimensional laser scanner. In order to facilitate grid division and simulation convergence, the UAV frame was simplified. The UAV model is shown in Figure 1. consistent with the turbulence intensity of the ambient wind, but it was still much larger than that of the ambient wind conditions. In summary, the research on UAV flow fields has mainly focused on the flow field under different flight conditions in outdoor environments and the blocking effect of crop convection field. There is still a lack of research on the characteristics of the UAV airflow field within a greenhouse structure. Therefore, in this study, a greenhouse planted with tomato was used as the research environment, and the CFD method was used to obtain the flow field simulation results of the UAV under different flight states in the greenhouse. The flow field distribution law of the UAV in the greenhouse was analyzed, and the accuracy of the simulation results was verified experimentally. These findings provide a theoretical reference for the application of UAVs in greenhouse plant protection.

Establishment of the Simulation Model
The three-dimensional model of the UAV was established based on the E360-S2 UAV of CNOOC. The main parameters are shown in Table 1. The rotor 8045 was obtained by an EXAscan handheld three-dimensional laser scanner. In order to facilitate grid division and simulation convergence, the UAV frame was simplified. The UAV model is shown in Figure 1.  Taking the cuboid of 17 m × 0.8 m × 1.2 m as the crop model, the crop was treated by the porous medium model to simulate the blocking effect of crops on airflow. The porous media model is based on Darcy-Forchheimer law and expressed as where S is the momentum loss term; μ is the dynamic viscosity, N·s/m 2 ; α is permeability, 0.395 m 2 [17]; ρ is fluid density, 1.225 kg/m 3 ; and C2 is the inertia loss coefficient, 0.4 [17]. According to a Venlo greenhouse, a three-dimensional model was established. Its parameters are shown in Table 2. The established greenhouse model is shown in Figure 2. The diagram simplifies the structure of air exchange inside and outside the greenhouse and only depicts its geometric structure.  Taking the cuboid of 17 m × 0.8 m × 1.2 m as the crop model, the crop was treated by the porous medium model to simulate the blocking effect of crops on airflow. The porous media model is based on Darcy-Forchheimer law and expressed as where S is the momentum loss term; µ is the dynamic viscosity, N·s/m 2 ; α is permeability, 0.395 m 2 [17]; ρ is fluid density, 1.225 kg/m 3 ; and C 2 is the inertia loss coefficient, 0.4 [17]. According to a Venlo greenhouse, a three-dimensional model was established. Its parameters are shown in Table 2. The established greenhouse model is shown in Figure 2. The diagram simplifies the structure of air exchange inside and outside the greenhouse and only depicts its geometric structure.  The three-dimensional model was assembled to obtain the greenhouse-crop-UAV model, as shown in Figure 3. The greenhouse contained a tomato crop with an average height of approximately 2 m during the period of measurement and a leaf area index of about 1.5 m 2 leaf per m 2 ground. The height and leaf area of tomato were slightly lower than that of Molina-Aiz [18].

Grid Division
ICEM-CFD was used to mesh the model. The greenhouse volume was over 1000 m 3 , and the maximum unit size was set to 0.3 m. The size of the UAV was small, and the grid division of the rotor had a great influence on the simulation results of the flow field. The grid of the UAV was encrypted, and the maximum unit size of the rotor was set to 0.0005 m. The grid model is shown in Figure 4. Number of rotor grids is shown in Table 3. The three-dimensional model was assembled to obtain the greenhouse-crop-UAV model, as shown in Figure 3. The greenhouse contained a tomato crop with an average height of approximately 2 m during the period of measurement and a leaf area index of about 1.5 m 2 leaf per m 2 ground. The height and leaf area of tomato were slightly lower than that of Molina-Aiz [18].  The three-dimensional model was assembled to obtain the greenhouse-crop-UAV model, as shown in Figure 3. The greenhouse contained a tomato crop with an average height of approximately 2 m during the period of measurement and a leaf area index of about 1.5 m 2 leaf per m 2 ground. The height and leaf area of tomato were slightly lower than that of Molina-Aiz [18].

Grid Division
ICEM-CFD was used to mesh the model. The greenhouse volume was over 1000 m 3 , and the maximum unit size was set to 0.3 m. The size of the UAV was small, and the grid division of the rotor had a great influence on the simulation results of the flow field. The grid of the UAV was encrypted, and the maximum unit size of the rotor was set to 0.0005 m. The grid model is shown in Figure 4. Number of rotor grids is shown in Table 3.

Grid Division
ICEM-CFD was used to mesh the model. The greenhouse volume was over 1000 m 3 , and the maximum unit size was set to 0.3 m. The size of the UAV was small, and the grid division of the rotor had a great influence on the simulation results of the flow field. The grid of the UAV was encrypted, and the maximum unit size of the rotor was set to 0.0005 m. The grid model is shown in Figure 4. Number of rotor grids is shown in Table 3.

Solution Settings
The total calculation domain was within the whole greenhouse structure. In order to ensure the solution accuracy of the numerical simulation, the calculation domain was divided into four rotation domains (set the rotation speed of 6000 rpm), porous media area, cultivation groove solid area, and air area in the greenhouse. Rotation calculation domain is shown in Figure 5. The interface between rotation and air area was connected by interface.

Solution Settings
The total calculation domain was within the whole greenhouse structure. In order to ensure the solution accuracy of the numerical simulation, the calculation domain was divided into four rotation domains (set the rotation speed of 6000 rpm), porous media area, cultivation groove solid area, and air area in the greenhouse. Rotation calculation domain is shown in Figure 5. The interface between rotation and air area was connected by interface. Since the k-ω SST model can better predict the flow and swirl in the near-wall region in the simulation of the external flow field, the k-ω SST turbulence model was set up. The SIMPLE algorithm was used for the coupling of pressure and velocity, and the multi-reference system model method was used to solve the flow field of the UAV.
In order to simulate the actual operation situation of the UAV to study the influence of flight height, inter-row position, and crop conditions (whether there are any crops) in the greenhouse on the airflow field of a UAV in the greenhouse, an orthogonal experiment (Table 4) was designed to simulate the flow field. Flgiht position layout of UAV in greenhouse is shown in Figure 6. Since the k-ω SST model can better predict the flow and swirl in the near-wall region in the simulation of the external flow field, the k-ω SST turbulence model was set up. The SIMPLE algorithm was used for the coupling of pressure and velocity, and the multireference system model method was used to solve the flow field of the UAV.
In order to simulate the actual operation situation of the UAV to study the influence of flight height, inter-row position, and crop conditions (whether there are any crops) in the greenhouse on the airflow field of a UAV in the greenhouse, an orthogonal experiment (Table 4) was designed to simulate the flow field. Flgiht position layout of UAV in greenhouse is shown in Figure 6.

Effect of Flow Field on Crops
After the 1.5 m east-west development of the UAV downwash flow field, the wind speed tended towards zero. It can be seen from Figure 7 that the crop areas on both sides of the UAV were affected by the flow field of the UAV. The maximum wind speed in the crop canopy area was 2.8 m/s, which is not conducive to the image data acquisition of crops by UAVs.

Effect of Flow Field on Crops
After the 1.5 m east-west development of the UAV downwash flow field, the wind speed tended towards zero. It can be seen from Figure 7 that the crop areas on both sides of the UAV were affected by the flow field of the UAV. The maximum wind speed in the crop canopy area was 2.8 m/s, which is not conducive to the image data acquisition of crops by UAVs.

Effect of Flow Field on Crops
After the 1.5 m east-west development of the UAV downwash flow field, the wind speed tended towards zero. It can be seen from Figure 7 that the crop areas on both sides of the UAV were affected by the flow field of the UAV. The maximum wind speed in the crop canopy area was 2.8 m/s, which is not conducive to the image data acquisition of crops by UAVs.

Effects of Crops on the Flow Field
In this study, a greenhouse with tomatoes was used to obtain the wind speed on both sides of the tomato crop near the UAV. The results are shown in Figure 8. The two velocity curves of the crops on the far side of the UAV were approximately coincident, and the maximum velocity difference was 0.07 m/s. By contrast, the wind velocity curves of the crops near the UAV were quite different, particularly in the range of 0.35-1.1 m away from the ground, where the two curves were obviously not coincident. In the vicinity of the crops, the flow field near the UAV side was not fully developed, and the wind speed was slightly lower. The maximum velocity difference between the two was 0.36 m/s. In general, tomato plants have a small leaf area and low leaf density, which have no obvious effect on airflow blocking. maximum velocity difference was 0.07 m/s. By contrast, the wind velocity curves of the crops near the UAV were quite different, particularly in the range of 0.35-1.1 m away from the ground, where the two curves were obviously not coincident. In the vicinity of the crops, the flow field near the UAV side was not fully developed, and the wind speed was slightly lower. The maximum velocity difference between the two was 0.36 m/s. In general, tomato plants have a small leaf area and low leaf density, which have no obvious effect on airflow blocking.

Effect of Greenhouse Structure on Flow Field
The UAV was limited by the height of the greenhouse, and so the flight height was generally less than 3 m. There was an obvious ground effect, and a near-surface regional flow field developed on the wall on both the north and south sides of the greenhouse. Occluded by the cultivation trough, the flow field extended 1.5 m in the east-west direction (Y-Z plan)and then decreased to 0 m/s, and the development area was far less than the north-south (X-Y plane) flow field, as shown in Figure 9. The cultivation tank was the main factor restricting the development of the UAV flow field.

Effect of Greenhouse Structure on Flow Field
The UAV was limited by the height of the greenhouse, and so the flight height was generally less than 3 m. There was an obvious ground effect, and a near-surface regional flow field developed on the wall on both the north and south sides of the greenhouse. Occluded by the cultivation trough, the flow field extended 1.5 m in the east-west direction (Y-Z plan)and then decreased to 0 m/s, and the development area was far less than the north-south (X-Y plane) flow field, as shown in Figure 9. The cultivation tank was the main factor restricting the development of the UAV flow field. maximum velocity difference was 0.07 m/s. By contrast, the wind velocity curves of the crops near the UAV were quite different, particularly in the range of 0.35-1.1 m away from the ground, where the two curves were obviously not coincident. In the vicinity of the crops, the flow field near the UAV side was not fully developed, and the wind speed was slightly lower. The maximum velocity difference between the two was 0.36 m/s. In general, tomato plants have a small leaf area and low leaf density, which have no obvious effect on airflow blocking.

Effect of Greenhouse Structure on Flow Field
The UAV was limited by the height of the greenhouse, and so the flight height was generally less than 3 m. There was an obvious ground effect, and a near-surface regional flow field developed on the wall on both the north and south sides of the greenhouse. Occluded by the cultivation trough, the flow field extended 1.5 m in the east-west direction (Y-Z plan)and then decreased to 0 m/s, and the development area was far less than the north-south (X-Y plane) flow field, as shown in Figure 9. The cultivation tank was the main factor restricting the development of the UAV flow field.   Figure 10 shows the flow field distribution of the X-Y plane and the crop near the UAV side when the UAV is in the middle of flight. From top to bottom, the flight height increased continuously. The downwash flow field of the UAV developed more fully downward and around, and the ground effect weakened, while the impact on crops increased. As shown in Figure 11, the affected area of the crops doubled with the increase in flight height by 40 cm. Figure 10 shows the flow field distribution of the X-Y plane and the crop near the UAV side when the UAV is in the middle of flight. From top to bottom, the flight height increased continuously. The downwash flow field of the UAV developed more fully downward and around, and the ground effect weakened, while the impact on crops increased. As shown in Figure 11, the affected area of the crops doubled with the increase in flight height by 40 cm. According to the 3D velocity volume data, all the velocity volume details were displayed on the 2D image at the same time, which was called velocity volume rendering. The velocity volume rendering was used to display the comprehensive distribution of the airflow field in the two-dimensional image, and through the control of the opacity, the situation of the velocity isosurface was reflected. Figure 11 shows a speed volume rendering of four locations in the middle of the row, in the quarter of the row, at the end of the row, and in the aisle. It can be seen that the flow field mainly developed between the rows when the UAV hovered. When the UAV hover position gradually moved to the northern According to the 3D velocity volume data, all the velocity volume details were displayed on the 2D image at the same time, which was called velocity volume rendering. The velocity volume rendering was used to display the comprehensive distribution of the airflow field in the two-dimensional image, and through the control of the opacity, the situation of the velocity isosurface was reflected. Figure 11 shows a speed volume rendering of four locations in the middle of the row, in the quarter of the row, at the end of the row, and in the aisle. It can be seen that the flow field mainly developed between the rows when the UAV hovered. When the UAV hover position gradually moved to the northern side of the greenhouse, the symmetry of the flow field distribution in the greenhouse disappeared, and the wind speed in the northern side of the greenhouse increased, resulting in a larger wind speed on the side of the UAV flight altitude, which affected the flight state of the UAV and caused a decrease in flight stability.

Effect of Flight State on Flow Field
Agriculture 2021, 11, x FOR PEER REVIEW 9 of 12 side of the greenhouse, the symmetry of the flow field distribution in the greenhouse disappeared, and the wind speed in the northern side of the greenhouse increased, resulting in a larger wind speed on the side of the UAV flight altitude, which affected the flight state of the UAV and caused a decrease in flight stability.

Test Validation
The UAV was suspended at 1.5 m height and had a midline position. The wind speed at 0.3 m, 0.6 m, and 0.9 m height was measured by a KIMO VT110 anemometer. Greenhouse test is shown in Figure 12. The test and simulation results are shown in Table 5.

Test Validation
The UAV was suspended at 1.5 m height and had a midline position. The wind speed at 0.3 m, 0.6 m, and 0.9 m height was measured by a KIMO VT110 anemometer. Greenhouse test is shown in Figure 12. The test and simulation results are shown in Table 5.
side of the greenhouse, the symmetry of the flow field distribution in the greenhouse disappeared, and the wind speed in the northern side of the greenhouse increased, resulting in a larger wind speed on the side of the UAV flight altitude, which affected the flight state of the UAV and caused a decrease in flight stability.

Test Validation
The UAV was suspended at 1.5 m height and had a midline position. The wind speed at 0.3 m, 0.6 m, and 0.9 m height was measured by a KIMO VT110 anemometer. Greenhouse test is shown in Figure 12. The test and simulation results are shown in Table 5.  Linear regression analysis was performed on the experimental and simulated values ( Figure 13). The regression fitting equation was y = x − 0.08139, the R 2 was 0.97486, and the significance level was 0.05. In summary, the numerical simulation results of the UAV greenhouse airflow field were accurate.  Linear regression analysis was performed on the experimental and simulated values ( Figure 13). The regression fitting equation was y = x − 0.08139, the R 2 was 0.97486, and the significance level was 0.05. In summary, the numerical simulation results of the UAV greenhouse airflow field were accurate.

Conclusions
During the flight of the four-rotor UAV in the greenhouse, the flow field was mainly hindered by the cultivation groove, and the flow field was not fully developed. The flow field was mainly concentrated at the bottom of the crop and the junction of the cultivation groove, and the flow field only affected the adjacent crops. Tomatoes have a small leaf area and low leaf density and thus had little effect on the flight flow field of the UAV in the greenhouse.

Conclusions
During the flight of the four-rotor UAV in the greenhouse, the flow field was mainly hindered by the cultivation groove, and the flow field was not fully developed. The flow field was mainly concentrated at the bottom of the crop and the junction of the cultivation groove, and the flow field only affected the adjacent crops. Tomatoes have a small leaf area and low leaf density and thus had little effect on the flight flow field of the UAV in the greenhouse.
When the UAV flight height increased, the flow field developed more fully downward and around, and the affected area of the crops became larger. For every 40-cm increase in flight height, the affected area doubled. When the UAV moved to one side between rows, the flow field distribution in the greenhouse was asymmetric, and the flight stability decreased.
Through the experimental verification, linear regression analysis was carried out on the measured value and the simulated value, and the goodness of fit (R2) was 0.97486, which verified the accuracy of the simulation results of the airflow field of the UAV in the greenhouse.

Limitations
This study did not consider the coupling effect of window ventilation and UAV downwash flow field. According to Espinoza et al.'s [19] research, window ventilation could significantly change the microclimate in greenhouse. Mostly, the air flow of the window ventilation is much larger than that of the UAV downwash flow field, for the whole greenhouse. When the window ventilation and the UAV downwash flow field worked together, the window ventilation factor may occupy the dominant position. However, due to the large local wind speed of the UAV downwash flow field, the downwash flow field may occupy the dominant position near the UAV.
It is necessary to further study and analyze the change law of greenhouse microclimate under the competitive coupling of window ventilation and UAV downwash flow field.