The computational fluid dynamic modeling of downwash flow field for a six-rotor UAV

Abstract

The downwash flow field of the multi-rotor unmanned aerial vehicle (UAV), formed by propellers during operation, has a significant influence on the deposition, drift and distribution of droplets as well as the spray width of the UAV for plant protection. To study the general characteristics of the distribution of the downwash airflow and simulate the static wind field of multi-rotor UAVs in hovering state, a 3D full-size physical model of JF01-10 six-rotor plant protection UAV was constructed using SolidWorks. The entire flow field surrounding the UAV and the rotation flow fields around the six rotors were established in UG software. The physical model and flow fields were meshed using unstructured tetrahedral elements in ANSYS software. Finally, the downwash flow field of UAV was simulated. With an increased hovering height, the ground effect was reduced and the minimum current velocity increased initially and then decreased. In addition, the spatial proportion of the turbulence occupied decreased. Furthermore, the appropriate operational hovering height for the JF01-10 is considered to be 3 m. These results can be applied to six-rotor plant protection UAVs employed in pesticide spraying and spray width detection.

Keywords

CFD simulation / downwash flow field / numerical analysis / plant protection / six-rotor UAV

Introduction

In the application of spraying during plant protection, the downwash flow field of multi-rotor unmanned aerial vehicles (UAV) directly drives the spraying droplets to the crop canopy. The coverage and distribution of the spray can directly impact on the spray range as well as the deposition and drift of droplets, so the operational effectiveness of the UAV can be significantly affected^[1]. Therefore, studying the distribution and interaction of the downwash flow field over a crop canopy is an important prerequisite for using multi-rotor UAVs in plant protection.

Recently, some research has focused on the distribution and properties of the downwash flow field of UAVs used in plant protection. One of the main methods was to use a wind speed sensor array to collect wind speed data used for inversion modeling. Li et al.^[2,3] used a wind speed collection system to collect the UAV wind speeds on X, Y and Z axis directions, and analyzed the distribution of the wind field and the influence from operational parameters on the wind field. Another approach was to use the method of Spatial Quality-balance. Wang et al.^[4] conducted a practical test in the field. They collected the average distribution proportions of droplets in the top and bottom parts for upwind and downwind directions separately under different conditions. Although the proportions were utilized as analysis indicators, the sensor array was still applied for the study of wind field distribution. In real measurements, natural wind is one of the significant factors contributing to measurement errors, a problem that still needs to be resolved. Some studies used sensors to detect the distribution of droplets to verify the effect of downwash airflow. Zheng et al.^[5] used a laser scanner to measure the spray width of nozzles in a wind tunnel and establish a 3D outline model of the fog fields, offering an optional method to verify the entire distribution of downwash with high accuracy.

Numerical analysis and computational fluid dynamics (CFD) simulation has been widely used as a precise and intuitive simulation approach. Most of the studies contributed to the design and test of aerodynamic characteristics of fuselage and rotor. Lei^[6] from the University of the Chinese Academy of Sciences, numerically simulated and tested the fluid field of multi-rotor UAV. He analyzed the new dynamic characteristics of a six-rotor UAV incorporating two suggested non-planar dual-rotor units, and also qualitatively analyzed the influence from the aerodynamic interaction between multi-rotors on multi-rotor UAVs in a viscous low-Reynolds-number environment. Based on the analysis, the design was optimized. Wang et al.^[7] adopted the method of numerical simulation of fluid dynamic to calculate numerically the 3D flow field around the UAV and propeller engine. The overall temperature distribution of the UAV, the temperature and exhaust concentration of the propeller engine in the high temperature area were simulated. Qiao et al.^[8] applied moving meshes and transition model in numerical simulation. They analyzed the conditions when the propeller was in front and behind the wing separately. The results were compared with the situation in which there was only a propeller or wing. Some other related studies have been conducted, such as the topological optimized method of the wing by Oktay et al.^[9], the gradient optimized method for wings by James et al.^[10], and the design of the UAV vessel and intake pipe CFD approach used by Papadopoulos et al.^[11]. These works are valuable for UAV frame or rotor design, but few studies have focused on UAV operational conditions.

Other studies used CFD in some other areas, such as in an orchard spray assistant system, spray droplet distribution simulation and equipment optimization. For instance, Melese Endalew et al.^[12,13] used CFD to analyze the drift of pesticides in orchards for an air-assisted spraying system for different spray types. Duga et al.^[14] studied the verification and development of the drift and application of the CFD model. Dekeyser et al.^[15] evaluated an orchard spraying system based on CFD and laboratory experiments. Bartzanas et al.^[16] improved the crop production system using CFD. Delele et al.^[17] studied the characteristics of droplets and influence on the plant surface using a CFD model. Although these studies were not in the same area, they offered useful information for the study of UAV spraying.

Recently, research has been initiated to simulate the operation of plant protection UAV using CFD. Shi et al.^[18] used this approach to study the major movement form and speed change of the downwash flow field for a small unmanned helicopter in ultra-low-altitude flight; Zhang et al.^[19] simulated the deposition and drift of pesticide droplets under constrained conditions. An experiment in corresponding conditions was conducted to verify the simulated results, but the results need to be improved. Although numerical analysis already has been applied in UAV operation, there is limited CFD analysis related to the downwash flow field for multi-rotor plant protection UAVs.

Therefore, this paper used CFD numerical analysis to simulate the static wind fields of a six-rotor plant protection UAV at different hovering heights, to obtain the entire distribution of the downwash flow field at different hovering heights. By comparing the wind fields at different hovering heights, the entire distribution pattern of wind fields was established. These findings could be used for the pesticide spraying as well as the spray width detection of the six-rotor plant protection UAV.

Construction of the full-size physical model of the six-rotor plant protection UAV

Six-rotor plant protection UAV

A JF01-10 six-rotor plant protection UAV (Fig. 1; Table 1) was used as the subject for numerical simulation.

The main characteristics of JF01-10 UAV are shown in Table 1.

The relative positions between wheelbase and propellers are shown in Fig. 2. The wheelbase of the UAV was 1200 mm. The rotation diameter of each propeller was about 558 mm (22′′), and the angle between two adjacent arms is 60°. Six propellers are uniformly distributed on the circle, of which the central point of the frame was the center and the wheelbase is the diameter.

Construction of full-size physical model

The 3D full-size physical model of JF01-10 six-rotor plant protection UAV was constructed in SolidWorks 2013. The model mainly consisted of three parts: main frame, arms and propellers. Under the premise of ensuring all characteristics, the model construction was simplified as much as possible, and all parts were constructed as solid parts so as to reduce the number of meshes and facilitate the following calculation.

Construction of the main frame and arms of the UAV

Dimensions of the central disk, arms, spray rods and undercarriage of the UAV are given in Table 2.

The three feet of undercarriage were located on the circle; the center coincides with the central disk and the diameter is 80 mm. Also, the three feet lie on the apices of an equilateral triangle. The full-size physical model was constructed according to the dimensions provided in Table 2 and shown in Fig. 3.

Full-size physical model of propeller

According to the literature, to design and analyze the structure and dynamic properties of propellers, the exact propeller model should be used, while in the analysis of wind field and pressure distribution produced by the propeller, a simplified model could be used^[18,^20]. This study concentrated on the entire distribution of UAV downwash flow field. So the driving electric motor and the propeller were fused together to simplify the model. This not only meets the simulation requirement of entire flow field, but also fully considered the detailed influence from the motors on the flow field around the propeller.

The length of the propeller was about 558 mm (22′′), and its geometric pitch was about 140 mm (5.5′′), with a standard cross section, driving electric motor diameter of 117.4 mm, and height of 96 mm. The 3D model is shown in Fig. 4.

Full-size physical model of the UAV

The propellers, arms and the main frame of the UAV were assembled in SolidWorks 2013, and different colors were applied to different parts, after which the full-size physical model of the UAV was established and shown as Fig. 5.

Partition of the fluid field and meshes

Partition of the fluid field

As shown in Fig. 5, the six-rotor UAV consisted of 6 groups of propellers (rotational parts) and other supporting and auxiliary mechanical structures (fixed parts). The fluid field is the space where the fluid passes through. Therefore, in the CFD simulation, all these solid parts in Fig. 5 were removed from the fluid field. Based on the structure and dimensions of the UAV, the total fluid field was divided into the rotational domain and a fixed domain.

Given the UAV’s maximum diameter of D, considering the size of space to be simulated, the total wind field was set to be a cylindrical region with a diameter of 10 × D. According to the goal of simulation, the wind field was separated into upwind and downwind fields. The height of upwind field was set to 5 × D, while the height of the downwind field was set at 2, 3 and 5 m, and unlimited height.

The fluid field was partitioned in UG engineering software as shown in Fig. 6, which shows the total field partition, where the middle part was the UAV and the rotation domain. The section between the upwind and downwind fields was the boundary region. Figure 6b shows the rotation domains of six propellers. The partition of downwash field is shown in Fig. 6c.

Construction of the model meshes

Based on the previous partition of fluid field, the meshes were generated in ANSYS software. Considering the efficiency and the fact that the relative dimensions of the supporting part and accessories were small, an unstructured tetrahedral element was adopted. The meshes on the small surface were increased. The resultant meshes are shown in Fig. 7; the total number of nodes was 4248871, and total number elements was 21141012.

Numerical calculation

SST turbulence model

In the two-equation eddy viscous turbulence model, the k–e model can reliably simulate the full developed flow of turbulence far from the wall, while the k–w model is widely adopted in solving problems for boundary layers in different pressure gradients. To take advantage of the merits of both models, Menter first proposed the SST turbulence model^[19], and subsequently this combined model has been gradually and widely accepted in engineering.

The SST turbulence model kept the k–w model near the wall, while k–e model was applied far from the wall, after which the two models were combined together with a switching function F₁. So the model is expressed as follows:

(1){\displaystyle \frac{D\rho k}{Dt}}={\displaystyle \frac{\partial }{\partial xj}}\left[(\mu +\sigma k\mu t){\displaystyle \frac{\partial k}{\partial xj}}\right]+\tau {\displaystyle \frac{\partial ui}{\partial xj}}-{\beta }^{*}\rho \omega k\dots \dots

(2){\displaystyle \frac{D\rho \omega }{Dt}}={\displaystyle \frac{\partial }{\partial xj}}\left[(\mu +\sigma \omega \mu t){\displaystyle \frac{\partial \omega }{\partial xj}}\right]+{\displaystyle \frac{\gamma }{vt}}\tau ij{\displaystyle \frac{\partial ui}{\partial xj}}-\beta \rho {\omega }^{\mathrm{2}}+\mathrm{2}\rho (\mathrm{1}-F\mathrm{1})\sigma \omega \mathrm{2}{\displaystyle \frac{\mathrm{1}}{\omega }}{\displaystyle \frac{\partial k}{\partial xj}}{\displaystyle \frac{\partial \omega }{\partial xj}}\dots \dots

In Eq. (2), the turbulence shear stress is \tau ij=\mu t\left({\displaystyle \frac{\partial ui}{\partial xj}}+{\displaystyle \frac{\partial uj}{\partial xi}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\displaystyle \frac{\partial uk}{\partial xk}}\delta ij\right)-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\rho k\delta ij and eddy viscosity coefficient is vt={\displaystyle \frac{k}{\omega }}.

Equations (1) and (2) were combined by the switching function F₁. In other words, F₁ was 1 if it was nearest the wall and then activated the k–w model, while F₁ was approaching 0 when the layer was away from the wall, activating the k–e model. Given Ф₁, Ф₂ as the constants in k–w model and transformed k–e model, respectively, the constant Ф in the combined model is:

(3)\varphi =F\mathrm{1}\varphi \mathrm{1}+(\mathrm{1}-F\mathrm{1})\varphi \mathrm{2}\dots \dots

The switching function F₁ is defined as:

(4){{\displaystyle F}}_{\mathrm{1}}=\tanh ({{\displaystyle \arg }}_{\mathrm{1}}^{\mathrm{4}})\dots \dots

(5)\arg \mathrm{1}=\min \left[\max \left({\displaystyle \frac{\sqrt{k}}{\mathrm{0.09}\omega y}},{\displaystyle \frac{\mathrm{500}v}{y^{\mathrm{2}}\omega }}\right),{\displaystyle \frac{\mathrm{4}\rho \sigma \omega \mathrm{2}k}{CDk\omega y^{\mathrm{2}}}}\right]\dots \dots

(6)CDk\omega =\max \left(\mathrm{2}\rho \sigma \omega \mathrm{2}{\displaystyle \frac{\mathrm{1}}{\omega }}{\displaystyle \frac{\partial k}{\partial xj}}{\displaystyle \frac{\partial \omega }{\partial xj}},{\mathrm{10}}^{-\mathrm{20}}\right)\dots \dots

In Eq. (5), y is the minimum distance from the present point to the wall.

To consider the effect of turbulent stress transfer, for the position where turbulent kinetic energy is generated more than dissipated, the eddy viscosity is limited and redefined:

(7)vt={\displaystyle \frac{a\mathrm{1}k}{\max (a\mathrm{1}\omega ,\left|\Omega \right|F\mathrm{2})}}\dots \dots

In Equation (7), \left|\mathbf{\Omega }\right| is the scalar constant for the eddy tensor, F₂ is:

(8){{\displaystyle F}}_{\mathrm{2}}=\tanh ({{\displaystyle \arg }}_{\mathrm{2}}^{\mathrm{2}})\dots \dots

(9){\arg }_{\mathrm{2}}=\max \left(\mathrm{2}{\displaystyle \frac{\sqrt{k}}{\mathrm{0.09}\omega y}},{\displaystyle \frac{\mathrm{500}v}{y^{\mathrm{2}}\omega }}\right)\dots \dots

Calculation environment and boundary condition

The working medium for the simulation was set as air, and the temperature was set to a constant 25°C. The heat transfer was ignored, and the relative flow pressure at inlet and outlet were set at zero. Since the simulation was for a six-rotor plant protection UAV in hover state, the entire calculation domain was set to 0 speed. Meanwhile, the outer wall was set by the non-slip boundary. The rotational velocity was set to 1500 r·min⁻¹, and the rotational axis and the direction of the rotating were set by real conditions. Furthermore, the boundary between the rotational and fixed domains was set to be the fixed-rotational boundary type, and the SST turbulence model was adopted to calculate 1000 steps, of which the convergence requirement was 10⁻⁴.

Analysis of the simulation results

The simulation results

The numerical analysis results of the downwash flow field for the multi-rotor UAV are shown in Fig. 8, where Fig. 8a–8d are the distributions of flow speeds for different hovering heights of 2, 3 and 5 m, and the relative boundless height, respectively.

The corresponding vector streamline diagrams of the downwash flow fields for different hovering heights are shown in Fig. 9, in which Fig. 9a–9d are the distributions of flow speeds for different hovering heights of 2, 3 and 5 m, and the relative boundless height, respectively.

The pressure distributions around the multi-rotor UAV at different heights, 2, 3 and 5 m, respectively, are shown in Fig. 10. Given the infinite height of the assumed relative boundless height, its result was not calculated or shown.

Analysis of the results

The distributions of speeds for the multi-rotor UAV were complex (Fig. 8). With the increase in height, the minimum flow speeds were 1.332 × 10⁻⁴, 2.017 × 10⁻⁴, 1.874 × 10⁻⁴ and 2.148 × 10⁻²⁶ m·s⁻¹ respective for hovering heights of 2, 3 and 5 m, and the relative boundless height. The trend of the minimum flow speed increased at first and then decreased till it approached 0. However, beside the boundless height, all the minimum flow speeds were of the order of 10⁻⁴ m·s⁻¹. Thus, the increase of the hovering height did not have a significant influence on the minimum flow speed.

With increased hovering height, the density of the streamlines reflected by the ground decreased (Fig. 8). Therefore, it could be inferred that with increased hovering height, the ground effect decreased and the flowing velocity became more concentrated until decreasing to 0 when the hovering height was boundless.

Regardless of the hovering height, the wind flow entered the region of rotors from above and exited downward (Fig. 9). However, as the hovering height increased, the ground effect decreased, and the fluid field was more uniformly distributed, so that the central part was more concentrated and moved vertically downward. Moreover, the turbulence in the fluid field was more stable, but the proportion of space occupied by the turbulence decreased compared to the whole downwash flow field of the UAV. When the hovering height was boundless, the entire flow field moved downward without turbulence.

The relative pressure (indicated by the color of the largest circular area in Fig. 10) around the multi-rotor UAV changed from 12.67 to 8.259, and 6.840 Pa, as the hovering height changed from 2 to 5 m, showing a decreasing trend. Finally, the relative pressure approached 0, when the height was boundless.

Discussion

Increased hovering height led to a more uniform and symmetric downwash flow field. Therefore, the flight control of the multi-rotor UAV in both manual and auto modes became less difficult, due to the reduction in the influence from the ground effect. This result is consistent with previous research^[22,23]. Moreover, the pesticide can be sprayed more uniformly, potentially resulting in better crop protection outcomes, which is also consistent with results of earlier research^[4].

With increased hovering height, the pressure around the multi-rotor UAV decreased, leading to reduced resistance to the spray droplets, which can influence the penetration of spray droplets. Additionally, the deposition of the droplets decreased and the uniformity as well as the penetration of droplets decreased, as the spray width increased.

For practical spraying, the selection of the UAV operation parameters should guarantee the uniformity and penetration of droplets. Meanwhile, the turbulence caused by the ground will be important, so that droplets can be deposited on the underside of leaves. Moreover, the spray width should be a reasonable value, otherwise a small spray width will increase the spraying duration. Considering these factors and the results for the JF01-10 six-rotor plant protection UAV, the best hovering height is considered to be about 3 m. This can be validated by comparing experimental results in field studies.

It is difficult for the existing equipment to detect the whole distributions of the downwash flow field in real-world experiments. However, the miniature UAV model could be utilized to develop the airflow because of its consistent proportions. Laser doppler velocimetry (LDV) could also be used in experiments to detect the distributions of airflow.

Conclusions

This study simulated the downwash flow field of a multi-rotor UAV intended for practical spraying operations by establishing a full-scale physical model as well as by numerical analysis. The downwash flow fields at different hovering heights were simulated and analyzed.

It is concluded that the current speeds and streamline distributions of UAV are relatively complex. As hovering height increased, the minimum current velocity increased then decreased. Also, the pressure around the UAV and the ground effect decreased. As a result, the flight control is facilitated and the turbulence is more stable, but the space proportion occupied is reduced.

Also, as a result of increased hovering height, factors, such as spray range, deposition uniformity and penetration, limit each other, which could affect the spatial distribution of droplets and thus influence the effectiveness of the spraying. In practical operation, JF01-10 UAV using a hovering height of 3 m is likely to be appropriate.

References

[1] An J, Xiang W, Han Z W, Xiao K T, Wang Z F, Wang X H, Wu J B, Yan P Z, Li J, Chen Y, Li J, Li Y. Validation of the Institute of Atmospheric Physics emergency response model with the meteorological towers measurements and SF6 diffusion and pool fire experiments. Atmospheric Environment, 2013, 81: 60–67

[2] Li J Y, Zhou Z Y, Lan Y B, Hu L, Zang Y, Liu A M, Luo X W, Zhang T M. Distribution of canopy wind field produced by rotor unmanned aerial vehicle pollination operation. Transactions of the Chinese Society of Agricultural Engineering, 2015, 31(3): 77–86 (in Chinese)

[3] Li J Y, Zhou Z Y, Hu L, Zang Y, Xu S, Liu A M, Luo X W, Zhang T M. Optimization of operation parameters for supplementary pollination in hybrid rice breeding using round multi-axis multi-rotor electric unmanned helicopter. Transactions of the Chinese Society of Agricultural Engineering, 2014, 30(11): 1–9 (in Chinese)

[4] Wang C L, He X K, Wang X N, Jane B, Andreas H, Wang Z G, Pan H Y, He Z Z. Testing method of spatial pesticide spraying deposition quality balance for unmanned aerial vehicle. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(11): 54–61 (in Chinese)

[5] Zheng Y J, Yang S H, Lan Y B, Clint H, Zhao C J, Chen L P, Liu X X, Tan Y. A novel detection method of spray droplet distribution based on LIDARs. International Journal of Agricultural and Biological Engineering, 2017, 10(4): 54–65

[6] Lei Y. Aerodynamics of a Hex-rotor SUAV: Numerical Simulation and Experimental Study. Dissertation for the Doctoral Degree, Jilin: Changchun Institute of Optics, Fine Mechanics and Physics Chinese Academy of Sciences, 2013

[7] Wang W H, Li J L, Wang C, Lv Y M, Wang G F, Liu G M. Three dimension fluid field numerical simulation for airscrew unmanned aircraft vehicle. Infrared Technology, 2012, 34(5): 292–296

[8] Qiao Y H, Ma D L, Li Z. Unsteady numerical simulation of propeller/wing interaction. Journal of Aerospace Power, 2015, 30(6): 1366–1373

[9] Oktay E, Akay H U, Sehitoglu O T. Three-dimensional structural topology optimization of aerial vehicles under aerodynamic loads. Computers & Fluids, 2014, (92): 225–232

[10] James K A, Kennedy G J, Martins J R R A. Concurrent aerostructural topology optimization of a wing box. Computers & Structures, 2014, (134): 1–17

[11] Papadopoulos F, Valakos I, Nikolos I K. Design of an S-duct intake for UAV applications. Aircraft Engineering and Aerospace Technology, 2012, 84(6): 439–456

[12] Melese Endalew A, Debaer C, Rutten N, Vercammen J, Delele M A, Ramon H, Nicolaï B M, Verboven P. Modelling pesticide flow and deposition from air-assisted orchard spraying in orchards: a new integrated CFD approach. Agricultural and Forest Meteorology, 2010, 150(10): 1383–1392

[13] Melese Endalew A, Debaer C, Rutten N, Vercammen J, Delele M A, Ramon H, Nicolaï B M, Verboven P. A new integrated CFD modelling approach towards air-assisted orchard spraying—Part II: validation for different sprayer types. Computers and Electronics in Agriculture, 2010, 71(2): 137–147

[14] Duga A T, Delele M A, Ruysen K, Dekeyser D, Nuyttens D, Bylemans D, Nicolaï B M, Verboven P. Development and validation of a 3D CFD model of drift and its application to air-assisted orchard sprayers. Biosystems Engineering, 2017, 154: 62–75

[15] Dekeyser D, Duga A T, Verboven P, Endalew A M, Hendrickx N, Nuyttens D. Assessment of orchard sprayers using laboratory experiments and computational fluid dynamics modelling. Biosystems Engineering, 2013, 114(2): 157–169

[16] Bartzanas T, Kacira M, Zhu H, Karmakar S, Tamimi E, Katsoulas N, Lee Bok I, Kittas C. Computational fluid dynamics applications to improve crop production systems. Computers and Electronics in Agriculture, 2013, 93: 151–167

[17] Delele M A, Nuyttens D, Ambaw A, Lebeau F, Nicolai B M, Verboven P. Studying the impact characteristics of spray droplets on plant surfaces using a multiphase CFD model. Aspects of Applied Biology, 2016, (132): 291–298

[18] Shi Q. Numerical simulation for downwash fluid field of small-size unmanned helicopter hedgehopping. Journal of Drainage and Irrigation Machinery Engineering, 2015, 33(6): 521–525 (in Chinese)

[19] Zhang S C, Xue X Y, Qin W C, Sun Z, Ding S M, Zhou L X. Simulation and experimental verification of aerial spraying drift on N-3 unmanned spraying helicopter. Transactions of the Chinese Society of Agricultural Engineering, 2015, 31(3): 87–93 (in Chinese)

[20] Wu J Q. Quadrotor layout on cross lifting body integrated design and investigation of aerodynamic characteristic. Dissertation for the Master Degree. Jiangsu: Nanjing University of Aeronautics and Astronautics, 2009 (in Chinese)

[21] Menter F R. Zonal Two Equation κ–ω Model for aerodynamic flows. AIAA Paper 93–2906, 24th Fluid Dynamics Conference, Orlando, FL, 1993: 69

[22] Wu Z, Wu T W. Influence of downwash flow on helicopter airborne missile initial trajectory.Journal of Projectiles, Rockets, Missiles and Guidance, 2009, 29(1): 202–204 (in Chinese)

[23] Qin Y H, Zhu Q H, Shao S. Aerodynamic characteristics analysis for hovering coaxial rotors in ground effect.Journal of Nanjing University of Aeronautics & Astronauts, 2015, 47(2): 266–274

[24] Qin W C, Xue X Y, Zhou L X, Zhang S C, Sun Z, Kong W, Wang B K. Effects of spraying parameters of unmanned aerial vehicle on droplets deposition distribution of maize canopies. Transactions of the Chinese Society of Agricultural Engineering, 2014, 30(5): 50–56 (in Chinese)

RIGHTS & PERMISSIONS

The Author(s) 2018. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)