Numerical Simulation and Experimental Verification of Rotor Airflow Field Based on Finite Volume Method and Lattice Boltzmann Method

Abstract

The primary focus of research in agricultural unmanned aerial vehicle (UAV) pesticide application technology is the investigation of droplet drift and deposition. The influence of the rotor airflow on droplets is particularly significant, making numerical simulations a crucial tool for airflow field analysis. Among existing numerical simulation methods, the Finite Volume Method (FVM) and the Lattice Boltzmann Method (LBM) are commonly used, but there is limited research that compares the two approaches. Therefore, this paper conducts numerical simulations of the rotor airflow of an agricultural UAV using Fluent, representing the FVM, and XFlow, representing the LBM. This research aims to reveal the distribution patterns of airflow field numerical simulations under different theoretical methods, validate them through practical experiments, and select the optimal method for simulating rotor airflow. The ultimate goal is to establish an effective airflow field model to enhance the precision of pesticide application by an agricultural UAV. The results indicate that the lift error calculated by XFlow in this paper is 2.57% smaller than that by Fluent. The wind field of Fluent entered the “stable state” earlier than XFlow, and the speed value of Fluent was smaller than that of XFlow. The difference between the two speed values became larger and larger as the distance from the rotor was longer. Compared with XFlow, Fluent changes more obviously in the core region, and the center region gradually disappears with the distance from the rotor. However, in the velocity field calculated by XFlow, there are still more turbulent flows outside the core region, indicating that the transient calculation method based on the LBM can better show the details of fluid flow than the steady-state calculation method based on the FVM. Through comparison with the actual test data, it is found that the relative error of the velocity value of XFlow at 0.2 m and 0.4 m is small, while that of Fluent at 0 and 0.2 m is small. It shows that XFlow simulation has higher accuracy for external turbulent flow, while Fluent simulation has higher accuracy for steady laminar flow. The research results provide data comparison and a basis for further analysis of the wind field model of the rotor wing of the plant protection UAV, and they lay a foundation for further research on precision application technology.

1. Introduction

The utilization of unmanned aerial vehicles (UAVs) for the protection of cultivated plants offers several advantages, including high operational efficiency, low operational cost, and the capacity to operate in a range of geographic areas and with a variety of crops. Therefore, this technology is a key enabler in the modernization of field management practices [1,2]. Since 2016, the scope of research into the use of UAVs for plant protection purposes has expanded rapidly, extending beyond rice to encompass a wider range of plant protection crops. For example, a study on the distribution of droplets between canopies under different operational parameters was conducted on a citrus orchard using a small UAV [3], a study on spraying droplets on maize with different growth periods was carried out using a JF01-10 plant protection UAV [4], a study on the prevention and control of wheat scab was conducted using the DJI T30 plant protection UAV [5], and a study on the relationship between payload and rotor speed was carried out using the TopXGun plant protection UAV [6]. These examples illustrate the extensive use of plant protection UAVs in modern precision agriculture.

As one of the most significant research areas in agriculture, the use of UAVs is a crucial aspect of this field. One of the most pivotal areas of research in this domain is the investigation of droplet drift and deposition. Additionally, understanding the influence of wind patterns on the dispersion of spray droplets is of paramount importance. Consequently, the analysis of the wind patterns associated with plant protection UAVs represents a fundamental element of research in this domain. The study of rotor wind fields typically encompasses three main aspects: field tests, wind tunnel tests, and numerical simulations. Field tests employ a range of materials, including water-sensitive paper, mylar sheets, Petri dishes, and polyethylene wires, to collect and observe moving droplets in a controlled environment. This approach indirectly demonstrates the rotor wind field by examining the deposition effect of the droplet and the performance of the droplet deposition [7,8,9]. Furthermore, researchers have employed a variety of measurement techniques, including wireless wind speed sensor networks, ultrasonic anemometer arrays, and tensiometers [9,10,11], to quantify the instantaneous wind field changes beneath the rotors of plant protection UAVs. However, the field test is susceptible to significant environmental influences, which introduces considerable uncertainty into the test results. Consequently, researchers have employed wind tunnel tests to observe the changes in the rotor wind field of plant protection UAVs, utilizing wind tunnels capable of generating stable wind conditions. In field conditions, the combination of optical equipment, including particle image velocimetry (PIV), a laser particle analyzer, and a phase Doppler particle analyzer (PDPA), is employed to measure flow field changes in rotor flow fields [12,13]. However, due to the limitation of the size of the test area, the wind tunnel can only be combined with some small plant protection UAVs to carry out relevant rotor wind field tests. In the case of the wind field of large plant protection UAVs, some researchers have also adopted the scaled-down model for experiments [14]. However, the scaled-down model must satisfy the conditions of similarity regarding fluid flow. It is also affected by the interference of the inner wall of wind tunnels and measurement errors. Furthermore, it is subject to significant limitations in terms of simulating the long-distance drift of aerial pesticide application. In comparison to field tests and wind tunnel tests, numerical simulation is distinguished by its low cost and short cycle time. It is capable of not only simulating complex flow field situations without the limitations and influences that are present in wind tunnel tests, such as the interference of the wind tunnel wall, but also of conveniently controlling and modifying various flow conditions in the fluid domain, thus allowing the study of the state of the plant protection UAV in a wind field environment.

The typical numerical simulation method is based on the discretization of macroscopic continuous equations. The majority of commercial CFD software employs the traditional Finite Volume Method (FVM) to simulate fluid flow changes, which assumes continuity. The Navier–Stokes equations and other related control differential equations are then discretized, resulting in a set of algebraic equations. These are solved to directly obtain the physical parameters of each discrete point in the solution domain. The differential equations are transformed into a system of algebraic equations, which are then solved to directly obtain the velocity, pressure, and other related physical parameters of each discrete point in the solution domain. The FVM can be employed not only in structured grids but also in unstructured grids, which are widely utilized and can be utilized to address the flow problems of arbitrary geometrical models. Consequently, numerous scholars have employed the FVM to investigate the distribution characteristics of the rotor wind field of plant protection UAVs. In the initial stages of research, some scholars employed the FVM to examine the velocity field and time evolution during the jetting process of a large monoplane, the Thrush 510 G, which is widely utilized in China. The trajectory of the wake vortex and the circulation state were simulated by the realizable k-ε model in the RANS equations based on the FVM, which provides a reference basis for subsequent research on the wake velocity field and the tip vortex [11]. Several researchers have employed this turbulence model to investigate the interaction between the downwash airflow field and the droplet field of various rotorcrafts. Some scholars employed the SST k-ε model, regarded as a “universal” model in RANS equations for industrial applications, to conduct simulations and analyses of the downwash flow field of an agricultural UAV, and they validated the credibility of the simulation data through tests [10,15]. In addition to the discrepancy in turbulence model selection, some scholars have also conducted research on the evolution and distribution characteristics of the downwash flow of agricultural UAVs by integrating the turbulence model to simulate rotor rotation using the dynamic mesh technique [16]. An alternative approach is the relatively less complex slip mesh technique, which is used to simulate the rotor blade rotation to effectively establish the study of the downwash airflow of agricultural UAVs. The maximum error between the experimental results and the numerical results does not exceed 15% [15,17]. To simulate the spraying effect during plant protection UAV operation, many research scholars have carried out the simulation of spraying droplets from agricultural UAVs by combining the discrete phase model (DPM) in Fluent and simulating the spray nozzles by using the flat-fan atomizer model to study the motion characteristics of the droplets under the action of the rotor wind field [15,18]. This has been performed to verify the validity of its numerical simulation through experiments. The simulation of flight conditions for both single-rotor and multi-rotor aircraft indicates that the existence of a critical flight speed presents a challenge to the conduct of realistic validation tests [19]. In the field of plant protection, some researchers have employed the CFD method of the RNG k-ε turbulence model in conjunction with the canopy, defined as a porous medium, to examine the downwash airflow of an agricultural UAV in a fruit tree. To validate the numerical simulation, they have utilized a wireless measurement system to gather data for the numerical simulation. As computer arithmetic power continues to advance, scholars engaged in the study of droplet motion during aerial spraying by agricultural UAVs have begun employing the PISO algorithm to compute the motion behavior of fluids in transient states. This approach enables a more comprehensive understanding of droplet evolution and provides valuable insights for subsequent researchers [15,20,21].

Conversely, when considering the macroscopic flow behavior of a fluid system, the influence of microscopic particles within the flow is minimal. The Lattice Boltzmann Method (LBM) is an algorithmic approach based on the micro-particle model and the mesoscopic equations of motion. This method describes the motion of the fluid system at the microscopic level, regarding the fluid medium as many discrete particle micro-clusters that exist in the interlattice according to the rules governing their collision and migration. The macroscopic motion law of the fluid is then obtained through the statistics of the particle motion characteristics on the lattice points. This simplified kinematic model, which yields macroscopically averaged properties that obey the adopted macroscopic equations, not only circumvents the necessity for the full Boltzmann equation, but also the tracking of each particle during molecular dynamics simulations. The LBM provides a solution to the challenges encountered by conventional macroscopic techniques in modelling interfacial flow, as well as flow in porous media, boiling dynamics, and dendritic flow, which are significant and pertinent engineering applications. Consequently, several researchers and scholars have employed LBM-based studies for research on the variation in the rotor wind field of plant protection UAVs. By leveraging the transient nature of the LBM, researchers conducted large-vortex simulations of a full-scale agricultural unmanned helicopter and its generated droplet motion, as well as actual hovering tests. The results demonstrated that the experimental and simulated transient flow structures are comparable [22]. Additionally, some scholars conducted numerical simulations of the downwash wind field of a quadcopter UAV under varying operational conditions with large vortex simulations. They also employed the Lagrangian discrete phase particle tracking method to simulate the motion trajectories of droplets with differing particle sizes and validated the precision of the numerical simulation through wind tunnel tests [23]. Furthermore, to fully exploit the transient computational advantages of the LBM, researchers conducted a spatio-temporal distribution characterization of the airflow field of the six-rotor agricultural UAV. The results demonstrated that the flight speed and flight altitude exerted a significant influence on the distribution of the airflow field. Additionally, the predicted values in the vertical direction based on the average velocity attenuation model (Y-DAVA) exhibited a high degree of alignment with the experimental measurements [24].

In the context of macroscopic hydrodynamic problems, the FVM is typically associated with high numerical accuracy and a wide range of applications, including the simulation of complex flow fields and irregular geometries. The numerical stability of the FVM is generally considered to be satisfactory in such cases. Numerical accuracy and order can be enhanced through the utilization of an appropriate numerical format and mesh refinement. However, when dealing with complex geometries and non-uniform media, specialized techniques are required to ensure numerical stability and minimize the error of the simulation. Additionally, the generation of non-structured meshes necessitates a certain investment of time and resources. In contrast, the LBM is a discretization method based on microscopic molecular motions. It treats the fluid system as a population of particles composed of molecules, simulates the motion of molecules at the lattice nodes, and simulates the macroscopic behavior of the fluid through collision and propagation steps. As a result, it performs well in dealing with complex boundaries and at the microscopic and mesoscopic scales. Furthermore, it can deal with the problems of multiphase flow and porous media flow with relative ease. In the case of high-speed flow or high Reynolds number, it is necessary to implement appropriate measures to enhance stability. Numerical accuracy is influenced by the velocity model, collision model, and numerical dissipation. Additionally, numerical error is introduced when dealing with the macroscopic scale. Nevertheless, no study has hitherto compared the two methods in relation to the rotor wind field of plant protection UAVs. Therefore, the present paper is dedicated to the study of the distribution changes in the rotor wind field of a plant protection UAV under the FVM and the LBM, with the view of elucidating the distribution of the rotor wind field under the two different methods. This study compares the changes in the physical quantities of the flow field of the rotor downwash wind field under different numerical simulation methods. It is hoped that this study will assist researchers in gaining a deeper understanding of the distribution characteristics of the rotor wind field under different methods, and in further elucidating the distribution characteristics and laws of the rotor wind field of the plant protection UAV under the application of different numerical simulation methods.

2. Materials and Methods

In the field of computational fluid dynamics, the main governing equations are the mass conservation equation (continuity equation), the momentum conservation equation, and the energy conservation equation, which correspond to the corresponding laws of mass conservation, momentum conservation, and energy conservation. Based on these three conservation laws, the Finite Volume Method (FVM) and the Lattice Boltzmann Method (LBM) have been developed gradually.

In the actual operational process, the rotor is situated at a considerable distance from the fuselage, and the influence of the fuselage on the evolutionary development of the downwash flow is minimal [17,24]. This paper focuses on the similarities and differences between the FVM and the LBM. To facilitate comparison, a single rotor of the plant protection UAV is selected for analysis. Given the pivotal role of the UAV rotor in generating the rotor flow field, it is of paramount importance to establish an accurate three-dimensional model for it. Considering the team’s prior research on the quadrotor plant protection UAV, this paper employs the blade of the 410S quadrotor plant protection UAV. The detailed parameters of the UAV are presented in

Table 1. Parameters of the quadrotor protection UAV.

Parameters	Technical Index
Rotor diameter/m	1.4
Rotor speed/rpm	0–3000
Load capacity/kg	10
Flight duration/min	8–10
Operating temperature/°C	−25–50

. Consequently, this paper utilizes the findings of the UAV HOBBYWING 3090 rotor model from a previous publication [12] to construct the 3D model of a single CW-type blade, as illustrated in Figure 1 below.

2.1. Simulation Model and Grid Independence Verification

2.1.1. Fluent Simulation Model and Grid Independence Verification

The FVM is based on the establishment of a discrete computational domain, which is defined as the computational space. This space is then divided into meshes, which are used to compute the required results. It is essential to ensure that the meshes within the computational domain do not overlap, as this could lead to inconsistencies in the results. Each discrete mesh constitutes a control body, with each mesh corresponding to a single control body. The control equations can be integrated over the aforementioned control bodies in order to obtain the corresponding discrete equations. Furthermore, the values of the unknown variables present in the equations are distributed over each control body, thus forming a variable field [25].

There are numerous software programs that facilitate the resolution of calculations through the FVM. Among these, Fluent is widely regarded as the most prominent, employing a discrete approach to locally identify approximate solutions. These solutions are then obtained through integration within the corresponding control volume. By assuming that the numerical solution in the entire calculation space is subject to a corresponding change rule, the weighted residual method, based on the sub-region method, allows for the selection of an appropriate region for integral processing. The crux of the matter is the discretization of the calculation region. This entails the reasonable and efficient calculation of the mesh, which in turn allows for the effective representation of the entire calculation region. The existence of the calculation mesh, which is a sub-region, is also a key consideration. Each sub-region is interconnected but does not overlap, with each sub-region representing the control body corresponding to the node, position, and so forth in this region. By integrating the control equations based on the fundamental principles of the FVM, the following Equation (1) can be derived:

$$\underset{V}{\int }{\displaystyle \frac{\partial \left(\rho \mathrm{\varnothing }\right)}{\partial t}}dV+\underset{V}{\int }\nabla ·\left(\rho U\mathrm{\varnothing }\right)dV=\underset{V}{\int }\nabla ·\left(\Gamma \nabla \mathrm{\varnothing }\right)dV+\underset{V}{\int }{S}_{\mathrm{\varnothing }}dV,$$

(1)

where the first term to the left of the equal sign denotes the transient term and the second term denotes the convective term, and the first term to the right of the equal sign denotes the diffusion term and the second term denotes the source term. The above Equation (1) can be rewritten by Gaussian dispersion theorem as follows:

$${\displaystyle \frac{\partial }{\partial t}}\left(\underset{V}{\int }\rho \mathrm{\varnothing }dV\right)+\underset{A}{\int }\mathrm{n}·\left(\rho U\mathrm{\varnothing }\right)dA=\underset{A}{\int }\mathrm{n}·\left(\Gamma \nabla \mathrm{\varnothing }\right)dA+\underset{V}{\int }{S}_{\mathrm{\varnothing }}dV,$$

(2)

where n denotes the unit vector normal to the surface of the control body being integrated.

Consequently, the FVM can be applied to a plethora of mesh types, including the most prevalent polyhedral meshes. Polyhedral meshes typically comprise a greater number of faces, which enables them to obtain more computational information in FVM interpolation calculations. This, in turn, facilitates superior convergence of the results.

Accordingly, the rotor blade wind field is solved and calculated using Fluent software based on the FVM. The blade is meshed in ANSYS 2022, after which the computational domain of a single rotor blade is created. This is set up as a cylinder with a radius of 1.3 m and a height of 2 m, consisting of a stationary domain and a cylindrical rotating domain containing the rotor blade (ignoring the effect of the fuselage on the wind field). To describe the turbulence within the computational region, the SST k-ω turbulence model based on the FVM has been selected. This model is particularly suited to the treatment of the near-wall region. This approach accounts for the transport of turbulent shear stresses and predicts the separation problem at smooth walls by adding an eddy viscosity limiter, which significantly enhances the precision of the numerical simulations [10]. The upper surface of the cylinder represents the velocity inlet, the lower surface is the wall, and the cylindrical surface is the pressure outlet. These define the X-axis and Z-axis, which are perpendicular to each other and parallel to the ground. The Y-axis is upward and perpendicular to the ground. The gravitational acceleration along the Y-axis is in the negative direction, with a value of g = 9.81 m/s². This is illustrated in Figure 2. The results of previous studies indicate that the downwash airflow of the plant protection UAV is primarily concentrated directly below the rotor. Consequently, in this study, the mesh density of the region below the rotor is encrypted, and the mesh size is gradually increased along the main distribution region.

Mesh is a fundamental component of numerical simulation, and the efficacy of the mesh delineation directly impacts the precision and outcomes of the numerical simulation calculations. Consequently, to ascertain the reliability of the numerical simulation results, it is not sufficient to rely solely on the accuracy of the calculation method; it is also necessary to ensure that the model has been meshed with an appropriate level of quality and quantity. However, due to the inherent limitations of computer hardware performance configurations and other factors, the calculation speed tends to decline with an increase in the number of meshes, while the time required for numerical simulation computations also exhibits a notable increase. Consequently, the capacity to reduce the number of meshes employed in numerical simulation is a significant challenge that engineers must address and resolve to ensure the reliability of computational results. The utilization of mesh independence verification represents a prevalent methodology for addressing such issues. This entails the execution of numerical simulations based on varying numbers of meshes. By selecting an appropriate meshing scheme, it is possible to guarantee computational efficiency while simultaneously minimizing the impact of the number of meshes on the computational results [26]. To minimize the computational errors caused by the different numbers of meshes, and at the same time ensure the accuracy and efficiency of the computation, four mesh schemes, designated A, B, C, and D, were established in Fluent for the mesh independence test. The various mesh schemes are presented in

Table 2. Fluent mesh independence verification scheme.

Scheme	Number of Meshes
A	4,057,443
B	4,377,825
C	5,210,616
D	6,304,492

.

The four different mesh schemes in the table above were simulated, with the rotor speed set to 2000 rpm in the software and residuals of 1 × 10⁻⁵. Once the residuals are stabilized at each time step, convergence is deemed to have been reached, and the number of iterations is 9000. The server configuration employed is a 48-core Intel^® Xeon^® CPU 3.00 GHz (Intel Corporation, Santa Clara, CA, USA). Once convergence is achieved, the calculation is terminated, and the lift values under disparate scenarios are obtained by organizing the simulation data obtained. In Fluent, the equation used to calculate lift is as follows:

$$L={\int }_s\left(-\mathrm{p}\mathrm{n}+\mathsf{\tau }\mathrm{t}\right)·e_{L}dS,$$

(3)

where L—the lift, N; S—the surface area of the object, m²; p—the pressure on the surface, Pa; $\mathrm{\tau }$—the shear stress, Pa; n—the normal vector to the surface; t—the tangential vector to the surface; and $e_L$—the unit vector in the direction of the lift.

The results of the mesh independence test are obtained by analyzing and comparing the results shown in Figure 3 below.

As illustrated in Figure 3, the outcomes of numerical simulations utilizing varying numbers of meshes exhibit disparate effects. With an increase in the number of meshes, the lift value rises from 29.52 N to 30.14 N, subsequently declining to 28.87 N. In comparison to the manufacturer’s reference value of 33.97 N, there is a discrepancy of 5.10 N, representing an error of 15.01%. Conversely, there is also a discrepancy of 3.83 N, representing an error of 11.27% in the minimum value. The results demonstrate that an excessively dense mesh may result in numerical dissipation when Fluent calculates the rotor blade, which is a complex geometry. This can lead to a larger error in the calculation of the lift value. Following a comprehensive analysis of the time required for numerical simulation and the associated error situation, it was decided that the subsequent calculations in Fluent would utilize Scheme B as the experimental comparison. The maximum skewness of the mesh in the computational domain of Scheme B is 0.71, with a total of 4,377,825 meshes.

2.1.2. XFlow Simulation Model and Grid Independence Verification

The employment of dynamic mesh technology in the analysis of rotor characteristics during high-speed rotation in dynamic simulation, utilizing the FVM, frequently results in the consumption of a considerable amount of computational time, largely due to the necessity of mesh reconstruction. Furthermore, truncation errors associated with mesh quality are generated during the calculation, which affects the accuracy of the resulting data. This may even cause the calculation to terminate due to negative volume during the calculation process. Accordingly, the fluid dynamics analysis software XFlow (Next Limit Dynamics S.L., Madrid, Spain), which is based on the LBM, is employed for the comparison and analysis of the rotating motion of a single rotor. The LBM does not necessitate mesh partitioning for rotor simulation, thereby conferring advantages in the resolution of three-dimensional flow fields characterized by complex boundary conditions and non-steady moving objects. It is readily capable of addressing fluid-related issues at both the micro and macro levels [27].

XFlow employs the LBM with a uniform cubic cell computational domain. The LBM is a mesoscopic method, and the macroscopic Navier–Stokes equation group can be derived from the lattice Boltzmann equation through the Chapman–Enskog expansion method. Of the numerous models of the LBM, XFlow employs a three-dimensional lattice structure, as illustrated in Figure 4. This structure encompasses twenty-seven velocity vector directions (D3Q27), comprising one discrete velocity vector with zero velocity situated at the lattice center, six discrete velocity vectors extending from the center to the lattice face center, twelve discrete velocity vectors from the center to the lattice edge midpoint, and eight discrete velocity vectors from the center to the lattice vertex. Accordingly, in comparison with the conventional LBM, it exhibits a more sophisticated spatial discretization pattern [28].

In this method, the lattice Boltzmann equation is selected as the solution equation, and the lattice Boltzmann transport equation is discretized on the lattice as follows [29]:

$$f_i\left(x+e_i,t+dt\right)=f_i\left(x,t\right)+W_i\left(x,t\right),$$

(4)

$$W_i={\displaystyle \frac{1}{\tau }}\left(f_i-f_i^e\right),$$

(5)

where $f_i$—particle velocity vector distribution function in the equation; $e_i$—velocity of particles along the i-th direction; $dt$—time step; $f_i(x,t)$—velocity distribution function of particles in the i-th direction at time $t$ and point $x$; $W_i$—collision operator; $f_i^e$—single-particle equilibrium distribution function; $\tau$—dimensionless relaxation parameter.

After the simplification based on the collision operator, the equation is reduced to the control equation describing fluid flow—Navier–Stokes equation—which can also exhibit the fluid state at low Mach numbers.

The blade model employed is illustrated in Figure 1, which is identical to that utilized in the FVM. The same configuration of the server has been selected, and the transient computation and LES turbulence model have been chosen. As it is not permitted to define a computational domain of cylindrical shape in XFlow, the computational domain is set as a cube in XFlow, with a larger volume than that defined in Fluent. The upper surface of the computational domain represents the velocity inlet, the lower surface is the wall, and the remaining surface is the pressure outlet, as illustrated in Figure 5a.

Although the LBM does not have a traditional mesh structure, the velocity distribution function is discretized on a regular lattice structure. However, further investigation is required to ascertain the influence of mesh resolution on numerical results. It is possible that different mesh resolutions may result in changes to the accuracy of the numerical solutions, particularly in cases where strong gradients and rich flow details are present. To verify that the accuracy of the simulation is not affected by the number of meshes, the global spatial discretization resolution of the computational domain is set to 0.1 m, the resolution of the lattice around the rotor blades is set to 0.025 m, and the resolution of the trailing edge of the airflow is set to 0.0125 m. The results of the spatial discretization with an initial time of 0 s are presented in Figure 5b,c, which illustrate the different zones representing different levels of refinement. The octree structure is divided into 4 levels, from level 0 to level 3, according to the refinement scale. In this paper, four sets of meshes with different schemes are employed for calculation purposes. The rotor speed is set to 2000 rpm, and different time steps are set to obtain the number of meshes at the moment of 0.5 s, as detailed in

Table 3. XFlow mesh independence verification scheme.

Scheme	Time Step	The Number of Meshes at 0.5 s
E	0.001	524,625
F	0.001	4,626,401
G	0.0001	441,633
H	0.0001	2,507,522

.

By organizing the obtained simulation data, the lift values of the rotor under different schemes are obtained. In XFlow, the lift calculation equation is as follows:

$$L={\int }_s\left(-\mathrm{p}\mathrm{n}+\mathrm{\tau }\right)·e_{L}dS,$$

(6)

where L—the lift, N; S—the surface area of the object, m²; p—the pressure on the surface, Pa; $\mathsf{\tau }$—the shear stress, Pa; n—the normal vector to the surface; and $e_L$—the unit vector in the direction of the lift.

Since the lift value in XFlow is in a fluctuating state, the numerical values are averaged to obtain the mesh independence verification chart shown in Figure 6.

As illustrated in Figure 6, the numerical simulations conducted by XFlow yield disparate outcomes due to the selection of distinct mesh schemes. It is evident that the impact of the number of meshes on the rotor speed in XFlow is not particularly pronounced when the time step is set to 0.001. The highest number of meshes is observed in Scheme F at the moment of 0.5 s. However, the value of the rotor lift does not correspond to the manufacturer’s reference value of 33.97 N. Furthermore, the difference between the rotor lift values of Scheme E and G is significant when the number of meshes is closer. Furthermore, when the time step is 0.0001, the discrepancy between the rotor lift values of the two schemes is relatively minor. In the comparison of schemes G and H, it can be observed that the rotor lift values are near the manufacturer’s reference value despite a significant difference in the number of meshes, with a variation of 8.83% and 8.47%, respectively. The results indicate that the number of meshes is not the primary factor influencing the lift value in the dynamic simulation of the rotor in XFlow. Instead, the time step has a significant impact on the simulation calculation of the rotor. Therefore, it is more appropriate to select a time step of 0.0001. The limitation of the LBM itself causes the calculation to diverge when the lattice viscosity is close to 0.5. Therefore, the value of the characteristic length cannot be too small; that is, the number of meshes used for the simulation cannot be too small. To ensure the stability of the software and to minimize the time cost of the calculation, the H scheme is chosen for the subsequent calculations in XFlow as the experimental comparison in this paper.

2.2. Field Test

2.2.1. Test Environment

The experimental site is located in the National Precision Agriculture Aviation Wind Tunnel Laboratory of South China Agricultural University (23°09′35.9′′ N 113°22′19.7′′ E). On the day of the experiment, the average indoor temperature was 27.2 °C, the average relative humidity was 67.1%, and the indoor air velocity was negligible.

2.2.2. Test Platform

The drone rotor used in the experiment was fixed on a plant protection UAV two-phase flow platform, which mainly includes a rotor platform, a spraying system, a control system, and a lifting device. The main body of the test platform is composed of aluminum alloy profiles, suspended under the gantry, with three adjustable attitudes ranging from −30° to 30°, which can support the use of multi-attitude drones such as quadcopters, hexacopters, and octocopters. Detailed information of the plant protection UAV platform is shown in

Table 4. Plant protection UAV platform.

Parameters	Technical Index
Platform height/m	3
Platform width/m	4
Operation duration time/min	0–15
Operation height/m	1–3
Motor speed/rpm	600–3000
Spray system pressure range/MPa	0–1.2
Maximum number of rotors supported	8

.

The rotor platform is specially customized and developed for this test platform. The motor speed of the motion platform is controlled in real time through the ground station software to achieve the corresponding wind field testing effect, while also monitoring the parameter information of the platform in real time and saving the test data, as shown in Figure 7.

2.2.3. Wind Speed Measurement

The wind speed data were collected and organized on site to verify the accuracy of the numerical simulation. As illustrated in Figure 8b, the corresponding rotor of the UAV was mounted on the plant protection UAV platform, and the rotor blades were aligned parallel to the ground. As illustrated in Figure 8a, the rotor was positioned at 2 m from the ground. A total of eight measurement planes, designated as H1-H8, were established with H1 situated at 0.3 m from the rotor and the remaining planes spaced at 0.2 m intervals. Each plane comprised five points, with the central point designated as P0 and the remaining points designated as P1, P2, P3, and P4. The P0 point was situated directly beneath the rotor, with the P1 and P3 points located at 0.2 m from the P0 point, and the P2 and P4 points situated at 0.4 m from the P0 point. As illustrated in Figure 8b, the anemometer was affixed to a bracket positioned parallel to the ground with the windward end oriented upwards. The position of the center point P0 was determined according to the plumb line. The anemometer employed was the Kestrel 5500 handheld weather meter, produced by Kestrel, Boothwyn, PA, USA. It was used to measure the vertical velocity (i.e., the velocity in the Z direction) of the downwash wind field with an accuracy of ±3%. The measurement range was 0.6 to 40 m/s, with a resolution of 0.1 m/s. The motor speed was regulated to 2000 rpm via the ground station of the plant protection UAV platform. Once the rotor speed had reached a stable state, the data pertaining to the measurement points were duly recorded. Each measurement point was recorded for a period of 30 s, with the anemometer set to sample every 1 s. This resulted in a set of data that was not steady over time. The data were obtained at the midpoint of the smoothing position, and the measurement data were repeated three times and averaged to ensure the accuracy of the resulting data set. This may serve to mitigate the impact of system error, which may be introduced because of the anemometer being susceptible to fluctuations in the surrounding environment.

3. Results

3.1. Distribution of Velocity Field at Different Times

Figure 9 illustrates the distribution of the rotor in the XOY plane at various points in time, employing the FVM and LBM in Fluent and XFlow, respectively. Figure 9(a1–d1) illustrate the distribution of rotor velocity at varying times in Fluent, while Figure 9(a2–d2) illustrate the distribution of rotor velocity at different times in XFlow. From Figure 9, it can be observed that when the rotor is situated at an elevation of 1.5 m above the ground, the velocity field is initially generated at 0.1 s, with much of the concentrated activity occurring around and below the rotor. This reaches the ground at 0.5 s. Both Fluent and XFlow demonstrate a symmetrical distribution about the vertical center of the rotor. However, there are still notable differences in the velocity fields generated by Fluent and XFlow. At 0.1 s, the velocity core generated by XFlow exhibits a distinctive shape in comparison to that generated by Fluent. The velocity core generated by XFlow is more concentrated, while that generated by Fluent diverges to both sides and is farther away from the rotor. Further investigation demonstrated that at 0.5 s, XFlow exhibited more pronounced upwash effects within the wind field in contact with the ground, whereas no upwash effects were observed in Fluent. After one second, the velocity field in Fluent reached a steady state, and the calculated results remained constant. In contrast, significant changes were observed in XFlow. The upwash wind field on the ground begins to extend and diverge in a manner that is characteristic of a turbulent flow, affecting the surrounding environment. This is because XFlow employs a transient calculation methodology, which results in a dynamic wind field. Further detailed insights into the wind field changes in quadrotors utilizing XFlow can be found in this research paper [12].

Figure 10 illustrates the variation in velocity values over time at two points situated 0.2 m from the center on the horizontal plane at heights of −0.5 m and −1 m, utilizing the Fluent and XFlow models. Figure 10a illustrates the velocity variation at a point on the horizontal plane at a height of −0.5 m. Figure 10b illustrates the velocity variation at a point on the horizontal plane at a height of −1 m. Figure 10 illustrates that at 0.3 s, the velocity below the rotor at a height of 0.5 m begins to approach a state of stability, while on the horizontal plane below the rotor at a height of 1 m, this occurs at 0.7 s. Figure 10a illustrates that the velocity at the point computed by XFlow stabilizes at approximately 14 m/s, while the velocity at the point computed by Fluent stabilizes at approximately 11 m/s. Figure 10b illustrates that the velocity at the point computed by XFlow stabilizes at approximately 12 m/s, while the velocity at the point computed by Fluent stabilizes at approximately 6 m/s. The velocity values obtained by XFlow are, overall, greater than those obtained by Fluent. Figure 10a illustrates that, despite the dynamic and time-varying nature of the velocity values computed by XFlow, the range of variation remains relatively limited and exhibits a tendency towards consistency. This suggests that the wind field in the −0.5 m plane below the rotor in XFlow has reached a state of equilibrium after 0.3 s. A comparison of the XFlow results in Figure 10b reveals that the discrepancy in velocity values within the wind field gradually increases. In conjunction with the wind field configuration in the XOY plane as illustrated in Figure 9, it can be postulated that the downwash wind field of the rotor reaches the ground at approximately 1 s, resulting in turbulence at the ground level and a notable alteration in the velocity values of the wind field near the ground. A further analysis of the Fluent results indicates that the velocity values in the −0.5 m plane tend to stabilize at 0.3 s, while those in the −1 m plane only begin to stabilize at 0.9 s, which is consistent with the XFlow results.

3.2. Distribution of Velocity Field at Different Altitude

Figure 11 illustrates the velocity distribution of a single rotor in the XOZ plane at varying heights, employing the FVM-based Fluent and LBM-based XFlow. In this representation, H represents the distance from the rotor center, with the negative Y-axis pointing downward below the rotor. From points A to F, the value of H increases, indicating that the distance from the rotor center also increases. Figure 11(a1–f1) illustrate the velocity distribution of a single rotor in Fluent, while Figure 11(a2–f2) illustrate the velocity distribution of a single rotor in XFlow. In general, the wind field velocity is concentrated below the rotor, with a region exhibiting lower wind speeds than those of the surrounding area in the central zone closer to the rotor. From (a1) to (f1), it can be observed that as the downwash airflow approaches the ground, the Y-component velocity in the velocity core region gradually decreases and disperses. Concurrently, given that the rotor utilized in the experiment is of the clockwise variety, the velocity core region also rotates in a clockwise direction, and the coverage area of the velocity core region gradually increases. As the distance from the rotor increases, the central area, characterized by lower wind speeds in the velocity core region, gradually diminishes. From (a2) to (f2), it can be observed that the downwash wind field calculated using the LBM is predominantly presented in a ring-shaped form, with a considerable degree of turbulence evident in the vicinity of the core region. Conversely, as the distance from the rotor increases, the central area represented by lower wind speeds in the core region also becomes more pronounced. Although there is also a decrease and dispersion of the Y-component velocity in the velocity core region as the downwash airflow approaches the ground, there are significant differences in the velocity core regions generated by the LBM and FVM. The reason is that Fluent uses a steady-state method, which averages the turbulent data through steady-state calculations and is insensitive to numerical dissipation and diffusion, which is not sufficient to capture the details of the flow and may result in distortion. On the other hand, the large eddy simulation method used in XFlow can directly simulate turbulence at smaller scales, and is more sensitive to numerical dissipation and diffusion, which can better capture the details of the vortex structure.

The distribution of the Y-component velocity has been demonstrated to exert a considerable influence on the efficacy of the spraying process, particularly in relation to droplet drift and deposition. It has the potential to suppress the drift of droplets in the air and facilitate the infiltration of droplets into plant leaves [5,16,30]. Figure 11 illustrates that the distance below the rotor is inversely proportional to the wind speed. This indicates that a considerable downwash velocity persists in the Y direction on the plane situated 1 m below the rotor.

Table 5. Maximum velocity values on each plane.

H (m)	V-Max (m/s)
	Fluent	XFlow
−0.1	14.08	14.29
−0.3	13.79	15.33
−0.5	12.32	15.37
−0.7	10.94	15.9
−0.9	10.4	14.8
−1.1	10	15.6

presents the maximum velocity values on each plane. The data in the table demonstrate that the maximum velocity declines gradually as the downwash airflow approaches the ground in Fluent. The maximum velocity calculated for the plane is 14.08 m/s, as illustrated in Figure 11(a1). As the descending airflow approaches the ground in XFlow, the maximum velocity remains relatively constant. The maximum value obtained from the planar calculation is 15.9 m/s, as illustrated in Figure 11(d2). Additionally, a comparison of the maximum wind speed values calculated in the two software programs reveals that the maximum wind speed value obtained from XFlow is consistently greater than that obtained from Fluent.

4. Discussion

The distribution data of the downwash wind field can be obtained through CFD simulation experiments and actual experiments. The relative error of the simulation results from Fluent and XFlow with respect to the experimental results is defined as follows:

$$e_F={\displaystyle \frac{\left|V_T-V_F\right|}{V_T}}\times 100\%,$$

$$e_X={\displaystyle \frac{\left|V_T-V_X\right|}{V_T}}\times 100\%,$$

where $V_T$ is the measured value in m/s, $V_F$ is the simulated value in Fluent in m/s, and $V_X$ is the simulated value in XFlow in m/s. $e_F$ is the relative error of the Fluent numerical simulation value compared to the experimental test value; $e_X$ is the relative error of the XFlow numerical simulation value compared to the experimental test value.

Table 6. Simulated values, measured values, and relative errors of marked points in Fluent, XFlow, and experimental tests.

The Horizontal Distance from the Origin (m)	Value	Measurement Height (m)
		0.3	0.5	0.7	0.9	1.1
0	$V_F$	3.71	6.67	7.80	7.79	6.85
	$V_X$	3.21	1.92	3.11	3.47	7.48
	$V_T$	4.51	7.08	8.53	8.71	8.95
	$e_F$	17.62	5.71	8.48	10.57	23.40
	$e_X$	28.77	72.85	63.57	60.10	16.32
0.2	$V_F$	13.02	11.20	9.49	7.76	6.04
	$V_X$	13.69	13.75	15.00	12.19	12.44
	$V_T$	10.72	10.52	10.24	10.16	10.15
	$e_F$	21.44	6.47	7.32	23.70	40.42
	$e_X$	27.71	30.75	46.48	19.92	22.58
0.4	$V_F$	0.57	2.27	4.09	3.26	3.22
	$V_X$	0.79	1.43	1.16	1.70	3.16
	$V_T$	0.75	1.23	1.76	1.41	2.16
	$e_F$	24.36	85.00	131.70	131.15	49.08
	$e_X$	4.84	16.95	34.15	20.27	46.41

presents the simulated, measured, and relative error values for the marked points in Fluent, XFlow, and experimental tests. From the data presented in the table, it can be observed that the calculated values of points on different horizontal planes at varying heights have been determined. Furthermore, the relative errors are summarized and plotted in a line chart, as illustrated in Figure 12. Figure 12 illustrates that the relative error of XFlow’s simulated values is relatively minor at 0.2 m and 0.4 m, with simulated values below 50%. However, the relative error is relatively significant at the 0 m position, and the relative error on the horizontal plane at heights of 0.5 m, 0.7 m, and 0.9 m all exceed 50%. In contrast, the simulated values produced by Fluent exhibit relatively minor relative errors at 0 m and 0.2 m, with many points displaying relative errors of approximately 20%. This suggests that the data are of a higher degree of accuracy. However, the relative error at 0.4 m, particularly on the horizontal planes at heights of 0.5 m, 0.7 m, and 0.9 m, is notably high, with relative errors exceeding 80%. While there may be instances where the simulated values of test points exceed the measured values, this is predominantly attributable to the loss of downwash airflow generated by the rotor in the actual test process, which is not accounted for in the simulation calculation.

The discovery of such elevated relative errors in the data prompted the researchers to conduct a more thorough investigation of the wind field data. Nine vertical line segments on the vertical plane, ranging from 0 to −1.0 m, were selected for analysis, as illustrated in Figure 13. The velocity values on each vertical line segment were then calculated and plotted in Figure 14. Figure 14a,b depict the velocity value curves of Fluent and XFlow, respectively, on the vertical lines. From Figure 14a,b, it can be observed that the velocity values of the corresponding line segments on either side of line 0 are nearly identical, indicating that the wind field beneath the rotor is predominantly symmetrical. Furthermore, it can be observed that the velocity values are typically higher at line segments 0.2 and −0.2 than at other line segments. This suggests that the maximum velocity under the rotor is situated on both sides of line segments 0.2 and −0.2. Nevertheless, a notable discrepancy is observed in the velocity values of Fluent and XFlow at line segment 0. In Fluent, the velocity values exhibit a gradual increase with distance from the rotor, whereas in XFlow, this trend is less evident. A further analysis of the data indicates that in Fluent, the velocity values at a vertical height of −0.5 m demonstrate a decreasing trend within the range of line segments −0.3 to 0.3. This suggests that the simulation values of the flow field in Fluent have reached a steady state and are no longer being calculated. In contrast, the simulation values in XFlow are perpetually undergoing a state of dynamic change, and they even exhibit an increase in velocity values as they approach the ground. This phenomenon is likely attributable to the ground effect previously discussed.

To further analyze the accuracy of the simulated wind field data, the researchers in this study select the average value of each line segment as a parameter for more in-depth analysis. The measured values are also averaged accordingly, and the results are plotted in a line chart as shown in Figure 15. From Figure 15, it can be seen that from line segment −0.4 to line segment 0.4, the trends of the three are basically the same, showing an “M”-shaped trend, but there are differences in the corresponding average value positions. The simulated values calculated using XFlow are generally less close to the measured values than those calculated using Fluent. The largest difference between XFlow and the measured values is at line segment 0, with a difference of 4.33 m/s, while the smallest difference is at line segment −0.4, with a difference of 0.01 m/s. At line segments −0.2 and 0.2, the maximum differences are 3.96 m/s and 3.62 m/s, respectively.

5. Conclusions

This study investigated the differences in the numerical simulation of a single rotor wind field using the FVM-based Fluent software and the LBM-based XFlow software. The Y-axis velocity of the wind field was tested to verify the simulation results. The main findings of this study can be summarized as follows:

(1).

Through grid independence verification, the minimum error of lift calculated by Fluent is 11.27%, while the minimum error of lift calculated by XFlow is 8.47. It can be seen that the calculation accuracy of the lift by Fluent is not as good as that of XFlow, indicating that the setting of model parameters in Fluent needs to be further improved. In XFlow, the influence of time step on rotor simulation is very important.

(2).

After the rotor wind field contacts the ground, XFlow produces a very obvious roll-up phenomenon, while the steady-state calculation results of Fluent do not show a very obvious performance. As the distance from the rotor increases, the center region that exists in Fluent gradually disappears, while the center region that exists in XFlow still exists in a more obvious form. Additionally, there is more turbulence around XFlow outside the core region. This indicates that the turbulence SST k-ω model in Fluent, by averaging the physical quantities, cannot capture the details of the wind field.

(3).

Compared with the test results, it is found that the average result of numerical simulation using Fluent is more likely to approximate the test value than that using XFlow. However, the wind field of the UAV is an unsteady system, and the steady-state method will lose the details of the numerical simulation stroke field. Next, it is suggested that relevant researchers carry out transient numerical simulation studies based on the FVM to compare the similarities and differences between FVM and LBM numerical simulation in transient states.