Study on Fluid Movement Characteristics inside the Emitter Flow Path of Drip Irrigation System Using the Yellow River Water

: Vigorously developing e ﬃ cient water-saving agricultural technologies using the Yellow River Water is an important way to achieve sustainable use of water resources. In order to clarify the ﬂuid movement characteristics inside the ﬂow path of the emitter under complicated water quality conditions in a drip irrigation system using the Yellow River Water, the optimal simulation turbulence model for the ﬂow ﬁeld in the ﬂow path of the emitter was determined by comparing the macroscopic hydraulic characteristics with the microscopic ﬂuid motion characteristics of the ﬂuid in the emitter. On this basis, the two-phase ﬂow model was used to calculate and analyze the characteristics of water ﬂow movement and particle transport in the emitter. The results show that the RNG (Re- normalization group) k- ε turbulence model was the most suitable for the simulation of the ﬂow ﬁeld in the emitter, considering the macroscopic hydraulic performance and microscopic anti-clogging ability of the emitter synthetically, and both the comprehensive calculation accuracy and the calculation e ﬃ ciency. The pressure showed a step-like uniform decrease along the direction of water ﬂow. The ﬂuid ﬂow showed the regional movement characteristics of the mainstream and non-mainstream regions. The energy dissipation mainly occurred at the sudden change sites of the ﬂow path structure. The particle phase velocity was slightly lower than that of the water phase. The velocity at the near-wall surface was relatively lower than that at the center, and the velocity distribution along the depth direction of the ﬂow path was relatively uneven. The sediment was mainly deposited in the ﬁrst half of the ﬂow path. This study can provide a theoretical basis for solving the emitter clogging in the drip irrigation systems applying water from the Yellow River.


Introduction
Water is the basic material condition related to the survival and social development of mankind, is a limited and irreplaceable valuable resource, and is also an important guarantee to realize the sustainable development of economy and society. While the world's economy has rapidly grown and great achievements have been made in various fields of construction, great resource and environmental costs have been paid. Water pollution is one such cost, which makes the already severe shortage of freshwater resources more insufficient. Combined with the unreasonable use of water, the world is facing a water crisis. These water problems are threatening the survival and development of mankind. Therefore, only by insisting on sustainable development and saving water resources can we make the economy develop well and quickly. Since the 18th National Congress of the Communist Party of China, General Secretary Xi Jinping has conducted on-site inspections of various tasks in the Yellow River

Calibration Method of the Optimal Simulation Model by CFD
Large eddy simulation model (LES) and Reynolds-averaged Naiver-Stokes equations (RANS) (including standard k-ε model and Re-normalization group (RNG) k-ε model) were selected for comparative analysis due to the relatively widely application in non-direct numerical simulation. The accuracies of different turbulence models for overall calculation were verified by flow-pressure relationship curves from a macroscopic point of view. And the accuracies of different turbulence models for calculating the fine motion characteristics of the internal fluid in the structural units of the flow path were verified by the results obtained with the digital particle image velocimetry (DPIV) technique from a microscopic point of view [25].
The hydraulic performance curve of the emitter was obtained by the experiment according to the requirements of "technical specification and test method of agricultural irrigation equipment emitter (GB/T17187-1997)" [26]. The flow state index and flow coefficient were calculated by the least square method based on changing the flow-pressure relation formula into a linear relation curve by logarithm, and it was used as a regression equation [27].

Physical Model Establishment and Meshing
A high-precision CT (Computed Tomography) scanning technique was used to accurately test the geometric parameters of the physical prototype of the emitter, and the three-dimensional mapping software UGS NX10.0 was used to establish a three-dimensional physical model of the flow path of the emitter. Meshing is the basis of numerical simulation. The quality of the mesh has an important impact on the accuracy and speed of the numerical simulation. In this study, ANSYS ICEM was used to mesh. The physical model was meshed by a hexahedral mesh with good quality and simple data structure. The near-wall surface and the corner of the tooth tip were encrypted to more accurately simulate the flow characteristics in the dramatic changes of the fluid motion region.

Solving Method and Boundary Conditions
The finite volume method was adopted to separate the governing equations in the numerical simulation, and the second-order upwind was adopted to separate the parameters such as the Note: h-teeth height; θ-teeth angle; s-distance between tooth; w-path weight; d-path depth; l-path length (the length of path center line).

Calibration Method of the Optimal Simulation Model by CFD
Large eddy simulation model (LES) and Reynolds-averaged Naiver-Stokes equations (RANS) (including standard k-ε model and Re-normalization group (RNG) k-ε model) were selected for comparative analysis due to the relatively widely application in non-direct numerical simulation. The accuracies of different turbulence models for overall calculation were verified by flow-pressure relationship curves from a macroscopic point of view. And the accuracies of different turbulence models for calculating the fine motion characteristics of the internal fluid in the structural units of the flow path were verified by the results obtained with the digital particle image velocimetry (DPIV) technique from a microscopic point of view [25].
The hydraulic performance curve of the emitter was obtained by the experiment according to the requirements of "technical specification and test method of agricultural irrigation equipment emitter (GB/T17187-1997)" [26]. The flow state index and flow coefficient were calculated by the least square method based on changing the flow-pressure relation formula into a linear relation curve by logarithm, and it was used as a regression equation [27].

Physical Model Establishment and Meshing
A high-precision CT (Computed Tomography) scanning technique was used to accurately test the geometric parameters of the physical prototype of the emitter, and the three-dimensional mapping software UGS NX10.0 was used to establish a three-dimensional physical model of the flow path of the emitter. Meshing is the basis of numerical simulation. The quality of the mesh has an important impact on the accuracy and speed of the numerical simulation. In this study, ANSYS ICEM was used to mesh. The physical model was meshed by a hexahedral mesh with good quality and simple data structure. The near-wall surface and the corner of the tooth tip were encrypted to more accurately simulate the flow characteristics in the dramatic changes of the fluid motion region.

Solving Method and Boundary Conditions
The finite volume method was adopted to separate the governing equations in the numerical simulation, and the second-order upwind was adopted to separate the parameters such as the convection terms. The coupling of speed and pressure was solved by the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm. For the turbulence model, the RNG k-ε two-equation turbulence model was selected after a comparative analysis of the flow change law with pressure and the velocity distribution inside the emitter measured by DPIV (Digital Particle Image Velocimetry) and simulated by the standard k-ε model, RNG k-ε model, and LES model. In the flow field calculation, the initial conditions are that the inlet of the emitter was set as a pressure inlet with 0.1 MPa and the outlet was set as a pressure outlet with 0 MPa. Except for the inlet and outlet of the calculation domain, all other fluid and solid surfaces were wall-type boundaries and set to non-slip boundaries. Standard wall functions were used to solve the near-wall flow. The multiphase flow model of standard Euler-Lagrange model was used to simulate particle motion. The sediment density was set to 2500 kg/m 3 , the volume fraction was 0.03, and the particle size was 100 µm.

Convergence Judgment and Postprocessing
The flow of the outlet and the residual value were used to judge the convergence. When the flow of the outlet was basically stable and the residual value was less than 10 −5 , the iteration was considered to converge. Tecplot 360 2011 was used for postprocessing the data.

Analysis Method of the Motion Characteristics of Water Flow and Particles
The pressure field, velocity field, and turbulence intensity distribution characteristics of the water flow in the emitter were mainly analyzed, and the particle motion of the full field, the velocity of the characteristic section, and the sand volume distribution were also analyzed.

Optimal Simulation Model Calibration of CFD
The flow change law with pressure and the velocity distribution inside the emitter are shown in Figures 2 and 3 measured by DPIV and simulated by the standard k-ε model, RNG k-ε model, and LES model. Compared with the measured values, the flow coefficient errors simulated by the standard k-ε model, RNG k-ε model, and LES model were 8.9%, 2.1%, and 1.4%, respectively; and the flow index errors were 3.8%, 0%, and 1.9%. In general, the calculation error of the Standard k-ε model was the greatest, and the calculation results of the LES model and the RNG k-ε model were relatively small. In Figure 3, the most accurate distribution position of the vortex area was simulated by the LES model, and the distribution position simulated by standard k-ε model had the farthest deviation from the position measured by DPIV. The characteristic flow velocity values simulated and measured in the characteristic area of the interior of the emitter (vortex center area, mainstream area, sharp angle disturbance area) are shown in Table 1. The velocity value errors between the values measured and calculated by the standard k-ε model, RNG k-ε model, and LES model in the center of the vortex were 12.50%, 8.33%, and 4.17%, respectively; in the mainstream area were 7.12%, 1.95%, and 1.30%; and the average speed errors at sharp angles area were 11.31%, 7.69% and 4.07%, respectively. The calculation error simulated by the LES model was the smallest, and the calculation error simulated by the Standard k-ε model was significantly larger than that simulated by the other two models.         Through macro and micro verification, the simulation accuracy by the LES model was the greatest, the one by the RNG k-ε model was in second place, and the one by the Standard k-ε model was the lowest. Although the simulation accuracy by the LES model was the highest, the LES model required Sustainability 2020, 12, 1319 6 of 12 more special mesh scale. The mesh must be encrypted enough to distinguish turbulent structures. The LES model was a non-steady-state calculation method. It took a long time to calculate the fluid motion and had high requirements for computer performance. Therefore, the LES model was not the main method for solving the complex turbulent flow field of the emitter. By the comprehensive analysis of the calculation accuracy and calculation efficiency, the RNG k-ε turbulence model was most suitable for the simulation of the flow field in the emitter.    Figure 6. The pressure showed a step-like uniform decrease along the flow direction, and the pressure sharply dropped at the sudden change of the flow path structure.      Figure 6. The pressure showed a step-like uniform decrease along the flow direction, and the pressure sharply dropped at the sudden change of the flow path structure.   Figure 6. The pressure showed a step-like uniform decrease along the flow direction, and the pressure sharply dropped at the sudden change of the flow path structure. At 0.1 MPa, the pressure distributions of the water flow along the flow direction at the longitudinal sections with the depth of 0.01 mm away from the cross-section of z = 0.5D of the emitter flow path are shown in Figure 6. The pressure showed a step-like uniform decrease along the flow direction, and the pressure sharply dropped at the sudden change of the flow path structure.  Figure 7. The fluid movement in the flow path was clearly divided into a high-speed mainstream region and a low-speed non-mainstream region. The velocity ranges in two regions were 1.5-2.6 and

Distribution of Turbulence Intensity
The turbulence intensity can be used as an indicator for characterizing the hydraulic and anticlogging performance of the emitter. The cloud map and contour map of the turbulent intensity distribution in the middle longitudinal section (z = 0.5D) of the flow path are shown in Figure 8. The turbulence intensity fields all showed the obvious intensity gradients. The turbulence intensity in the mainstream region was higher, and the turbulence intensity in the non-mainstream region was lower. The extremum of turbulence intensity appeared near the sudden change position of the flow path structure, which illustrates the fluid was the most disturbed and the fluid has the strongest turbulence ability, and it gradually decayed along the flow direction. This is consistent with the distribution results of the velocity field and pressure field inside the emitter measured by DPIV.

Distribution of Turbulence Intensity
The turbulence intensity can be used as an indicator for characterizing the hydraulic and anti-clogging performance of the emitter. The cloud map and contour map of the turbulent intensity distribution in the middle longitudinal section (z = 0.5D) of the flow path are shown in Figure 8. The turbulence intensity fields all showed the obvious intensity gradients. The turbulence intensity in the mainstream region was higher, and the turbulence intensity in the non-mainstream region was lower. The extremum of turbulence intensity appeared near the sudden change position of the flow path Sustainability 2020, 12, 1319 8 of 12 structure, which illustrates the fluid was the most disturbed and the fluid has the strongest turbulence ability, and it gradually decayed along the flow direction. This is consistent with the distribution results of the velocity field and pressure field inside the emitter measured by DPIV.
The turbulence intensity can be used as an indicator for characterizing the hydraulic and anticlogging performance of the emitter. The cloud map and contour map of the turbulent intensity distribution in the middle longitudinal section (z = 0.5D) of the flow path are shown in Figure 8. The turbulence intensity fields all showed the obvious intensity gradients. The turbulence intensity in the mainstream region was higher, and the turbulence intensity in the non-mainstream region was lower. The extremum of turbulence intensity appeared near the sudden change position of the flow path structure, which illustrates the fluid was the most disturbed and the fluid has the strongest turbulence ability, and it gradually decayed along the flow direction. This is consistent with the distribution results of the velocity field and pressure field inside the emitter measured by DPIV.

Distribution of the Velocity Field
The distributions of the velocity vector of water phase and sand phase in the intermediate structural unit (seventh and eighth structural unit) are shown in Figure 9. The velocities of water and particles were basically the same, and the velocity of the water phase was slightly greater than the velocity of the particle phase. This is because the two-phase flow model used in this calculation was the Eulerian model. In this model, the momentum equation was established for each phase in the calculation domain during the calculation, and the multiphase flow was regarded as a continuous medium. Therefore, the velocity distribution and values of the water phase and particle phase were relatively close.      Figure 10. The velocity distributions of the two phases were 246 relatively close, and the velocity of the particle phase was slightly lower than that of the water phase.

247
The velocity near the wall of each section was lower. The velocity distribution along the depth of the 248 Figure 9. The distribution of the velocity field of water phase and sand phase. Note: A-the area of the tooth tip dorsal water; B-the area of the tooth tip facing water; C-the area of the tooth heel dorsal water; D-the area of tooth heel facing water; the colors mean the velocity value, and arrows mean flow direction. Seven representative cross-sections were selected in the flow path structural unit (shown in Figure 9), and the contour maps of the velocity distribution of water and particulate matter at the different cross-sections are shown in Figure 10. The velocity distributions of the two phases were relatively close, and the velocity of the particle phase was slightly lower than that of the water phase. The velocity near the wall of each section was lower. The velocity distribution along the depth of the flow path was relatively uneven, with the maximum velocity and the largest change gradient at the tooth tip.
Distribution of the Sand Phase The particle volume distributions along the flow path depth direction and different cross-sections are shown in Figure 11. Under certain hydraulic condition and particulate matter conditions, the amount of particulate matter that the water stream can carry was limited. With the increase of the depth of the flow path, the volume fraction of the particle phase at the same position of the flow path showed an increasing trend, indicating that during the forward movement of the water flow, the particles were continuously deposited under the action of gravity. Sediment was mainly deposited in the first half of the flow path.
According to the particle phase volume fractions at different cross-sections of the flow path, the particle content distributions were very uneven along the depth direction of the flow path. Along the movement direction of water flow, the particulate matter contents in the area of the tooth tip facing water (B area) and the area of the tooth heel dorsal water (C area) were the highest, and the particulate matter contents were the lowest in the area of the tooth heel facing water (D area) and the area of the tooth tip dorsal water (A area). Therefore, the structure of the A and C areas needs to be optimized to make the sediment hard to deposit and easy to discharge outside the emitter.
tooth tip dorsal water (A area). Therefore, the structure of the A and C areas needs to be optimized to make the sediment hard to deposit and easy to discharge outside the emitter.

Discussion
In this paper, the RNG k-ε turbulence model was determined as the optimal simulation turbulence model of the flow field in the emitter flow path by comparing the macroscopic hydraulics and microfluid motion characteristics. Based on this model, the flow characteristics of fluid and particulate matter transport were calculated and analyzed. Through the analysis of the flow characteristics of the flow movement inside the emitter flow path, the pressure was reduced from the maximum to 0 after flowing through the whole flow path. This is due to the continuous energy dissipation effect of each structural unit of the flow path. When the fluid flowed through the whole path, its energy was eventually consumed, resulting in the formation of water droplets [28]. In addition, according to the distributions of the pressure field, velocity field, and turbulence intensity, it can be seen that there were significant changes of the pressure and velocity at the abrupt change of the flow path structure, and the fluid was greatly disturbed. This is because the abrupt change of the flow path structure (tooth tip) caused the fluid movement to be sharply adjusted, the flow velocity distribution to be reorganized, and the flow line to be denser. In this process, the relative motion between the viscous liquid particles was strengthened and the internal friction was increased, and then the energy loss was more serious. Therefore, the abrupt region of the flow path structure was the main area of the energy dissipation, and the effect of energy dissipation was the most significant [29].
By the analysis of the internal fluid flow and particle movement characteristics in the emitter flow path, it can also be found that the flow in the flow path presented the mainstream area and the non-mainstream area, and there were vortices. The formation of a low-speed vortex in the non-main zone is because of the following three factors: the fluid velocity differences between two zones at the junction of the mainstream zone and the non-main zone, the fluid viscosity, and the obstruction of the flow path wall. Due to the low velocity in the vortex zone, it is easy to deposit the clogging material, which is also the essential reason for the emitter clogging [30]. Therefore, it is necessary to improve the flow path structure of the emitter to make the vortex fully developed as far as possible, to increase the velocity of the vortex zone and its scouring effect on the inner wall of the emitter, and to make the clogging material not deposit and easily discharge with the water flow. This is also the most fundamental way to improve the anti-clogging performance [31].

Conclusions
The following conclusions can be got: (1) The RNG k-ε turbulence model was the most suitable for the simulation of the internal flow field of the emitter based on the macroscopic hydraulic performance and the microclogging resistance of the emitter, the comprehensive calculation accuracy, and the calculation efficiency.
(2) The pressure showed a step-like uniform decrease along the flow direction in the flow path of the emitter. The fluid flow presented regional motion characteristics-the mainstream region and non-mainstream region. The abrupt change sites of the flow path structure were the most important part of the energy dissipation of the emitter.
(3) The velocity of particle phase was slightly lower than that of the water phase. The near-wall velocity was lower than that in the center of the flow path, and the velocity along the depth of the flow path and the volume distribution of the sand phase were relatively uneven. The sediment was mainly deposited in the front half of the path.
The research in this paper only studies the water flow and particulate matter movement characteristics of a certain kind of emitters and the flow path structure under the condition of drip irrigation from the Yellow River. The most important thing is that the water quality of the Yellow River is complex, and the impact of complex water quality on the simulation results is not fully considered in the simulation process. Therefore, in future research we suggest: (1) Simulating the fluid motion characteristics of various emitter products and flow path structures under the condition of drip irrigation from the Yellow River.
(2) In determining the CFD simulation model, consider the impact of water quality on model selection.