Abstract
This study aimed to determine the split ratio, flow-field structure, and effect of different shaped channels to sudden pollution accidents in a generalized complex canal system of a wetland park, both experimentally and numerically. The three-dimensional instantaneous velocities at a typical section of each channel in the generalized model were measured experimentally using an acoustic Doppler velocimeter. The results showed that the split ratio calculation formula of three parallel channels could be derived under the condition of considering the frictional head and the local head losses. The water depth, velocities, and pollutant diffusion were widely influenced by changes in the cross-sectional shape and channel plane shape. The pollutants were trapped by stable vortices and transverse circulation due to shear force and secondary flow, thus delaying the diffusion of pollutants. The research results reported herein can help provide technical support for the normal operation of complex canal systems.
1. Introduction
In the planning and construction of cities and towns where water is available and recyclable, complex canals with multiple inlets (and outlets), including series (and parallel) branch channels and multivariate plane shapes, are important components of a town’s water ecosystem. The canal system can beautify the landscape and change its spatial structure. When a lake or reservoir is connected by a canal system, it can also regulate floods, material flow, information flow, and species flow. Therefore, it is of great significance to study the flow-field structure of a complex canal system and its influence on pollution accidents for the planning and operational management of water ecological towns.
At present, the flow characteristics and pollutant diffusion characteristics of natural rivers have been extensively studied by several scholars worldwide. In a study of curved rivers, the flow characteristics of single bend and continuous bend become the focus of research. On this basis, scholars also make a more comprehensive study on the flow of bend by changing the curvature, roughness, slope, and so on. The existence of centrifugal force in the bend leads to the formation of secondary flow, which also causes the local asymmetry in the variations of components of velocity. Secondary flow causes the lateral momentum transports toward the outer wall and disturbs the fractional velocity distribution of the longitudinal velocity component [1, 2]. At the top of the bend, the longitudinal average velocity distribution is close to a parabola, and the maximum point is below the water surface, and the overall average velocity is large on the top and small on the bottom [3]. Therefore, many scholars [4–7] have studied the movement of secondary flow, the change of water depth in transverse and longitudinal direction, and the distribution of bed shear force in the bend and found that there are multiple vortices in the secondary flow in the continuous bend, and there is a region of secondary flow direction reversal near the riverbed. The secondary flow also affects the redistribution of downstream velocity, and the transport of secondary flow also affects the redistribution of bed shear stress. In addition, Esfahani et al. [8] studied three curves with different curvatures and found that, contrary to the single sharp bend, multibend effects hinder the transfer of the kinetic energy in both directions in the entrance section of the strongly curved bend. Stoesser et al. [9] found that the RANS model can predict a primary as well as an outer bank secondary cell, and this outer-bank cell is found to be the residual of the main secondary cell of the previous bend and continues throughout the bending process. Blanckaert first studied three kinds of hydrodynamic processes at the sharp bend of the fixed bank [10], then added different roughness on the outer bank of the concave bank of the bend, and found that the increase in the roughness on the outer bank would lead to the considerable expansion and strengthening of the secondary flow on the outer bank [11]. At the same time, he also explored the location and reason of the separation of the flow on the convex bank of the open channel [12]. In addition to the existence of secondary flow in bends, branching or confluence rivers also exist. Therefore, Chen et al. [13] and Lyubimova et al. [14] found that while twin secondary circulations are found at channel confluence, their contribution to the mixing depends on their local positions with respect to the river streams, and the water depth drives the most rapid mixing and notably. Thomas et al. [15] evaluated the influence of secondary flow disturbance on bifurcation evolution. Jing et al. [16] used the RNG turbulence model to study that with the increase in unit width flow ratio, the recirculation zone in the splitter area increases. Balachandar et al. [17] and Singha et al. [18] studied the velocity distribution and turbulence characteristics in wide-and-narrow channel. Wang et al. [19] and Papanicolaou et al. [20] measured the three-dimensional instantaneous velocity of typical cross section of wide-and-narrow channel using an acoustic Doppler velocimeter (ADV) and analyzed the turbulence characteristics. Asnaashari et al. [21] used the turbulence model based on the VOF method to study the secondary flow distribution in a channel with contraction and diffusion sections. Aimed at the diffusion of pollutants in the confluence channel, Chen et al. [22] observed the pollutant diffusion in an asymmetric open channel at different intersection angles and pollutant concentrations and analyzed the distribution law of the pollutant concentration field. Hua et al. [23] found that the release position of pollutants was the main influencing factor in their transportation process when they discussed the mixing characteristics of pollutants under different intersection conditions and analyzed the plane flow field and mixing characteristics of pollutants from the bottom to the surface. Azevedo et al. [24] verified the mathematical model by water level and velocity in their experiment and analyzed the influence of different discharge conditions on the dynamics and pollutant diffusion of an estuary. Tang et al. [25] studied the influence of riverbed shape on pollutant concentration distribution based on the Reynolds equation and equation turbulence model.
Through the above research, we can know that previous studies have focused on flow movement characteristics and the pollutant diffusion law in a single bend, continuous bend, wide-and-narrow channel, or confluence channel. However, these studies were aimed at single channels, and the channel connections were simple. Conversely, research on canal systems with complex connection relationships and with many changes of section shape along the way—such as the complex canal system composed of gradually expanding and shrinking channels and continuous bends—has been scant.
The present study focused on the complex canal system of a generalized wetland park. The velocity field of the model canal system was measured using an ADV, and the flow-field structure and pollutant diffusion of the canal system were simulated using a three-dimensional turbulence model. The split ratio, flow movement characteristics, and effect of pollution accidents of each channel in the complex canal system are unique features of this paper.
2. Experimental Model and Numerical Method
2.1. Experimental Set Up and Measurement Method
The general scale of the generalized model was 5 × 3 m (length × width). The canal system consisted of two inlet channels in the upstream, an upstream pool, three parallel channels, an upstream pool, and downstream outlet, where the downstream outlet being connected to the outlet channel. The channels and pools were made of plexiglass plates with a surface roughness of 0.008. From left to right (see Figure 1), the three parallel channels—all with a rectangular cross-section—were the channel with a wide crest weir, a gradually expanding and shrinking channel, and plane “W-shaped” channel, respectively. The height of the three channels is 0.2 m. The channel with a wide crest weir was 0.2 m wide, and the wide crest weir was 0.3 m wide and 0.04 m high. The width of the gradually expanding and shrinking channel expanded from 0.1 to 0.2 m and then decreased back to 0.1 m, and the gradually expanding and gradually shrinking sections were 0.1 m long. The plane “W-shaped” channel was of 0.15 m width and consisted of five channels of lengths 0.5, 0.8, 0.9, 0.7, and 1.4 m, respectively. The channel connection angles were 120°, 120°, 110°, and 110°. Figure 1 shows a schematic view of the channel and the experimental setup.

This experiment uses the precision 0.1 mm test needle to carry on the discharge control. The water depth of a typical section of each channel was measured using a steel ruler with an accuracy of 0.5 mm. If the water surface fluctuated, the water depth value was determined by the average value of five measurements. The velocity was measured using an ADV with a side probe and the highest frequency of 200HZ. The measurement range of ADV is 1 mm/s∼4 m/s, and the measurement accuracy is ±0.5%. In this experiment, the frequency of 100 Hz and sampling time of 10 s were used to measure the velocity.
2.2. Numerical Simulation Method
In this study, the software Flow-3D was employed to the numerical modelling, and the renormalization group RNG turbulence model based on the VOF method was adopted in this software as it deals with curved streamline flow well. The continuity, momentum, equations of the RNG turbulence model are as follows [26, 27]: The continuity equation is as follows: The momentum equation is as follows: The k equation is as follows: The ε equation is as follows:where is the coordinate component, is the velocity component in the i-direction, is the unit volume force, is the pressure intensity, is the density of water, is the dynamic viscosity of water, is the eddy viscosity, is the turbulent kinetic energy, is the turbulence energy dissipation rate, is the effect of the average strain rate on , is the ratio of the turbulence time scale to the mean flow time scale, S is the strain rate tensor, and is the generation term of turbulent kinetic energy, where The numerical model has its own model scalar; that is, scalar is used in standard models, e.g., Sediment Scour and Deposition, Defect Tracking, Air Entrainment, or Fluid residence time. Users can also customize scalars, which are not used in any standard model, and can be introduced as fluid markers (similar to Marker Particles). In this study, the user-defined method will be used to set the scalar for the simulation of sudden pollution accidents. At the same time, through the mass momentum source button, the position, shape, and direction of the defined source are set, and the solver automatically generates the defined scalar at the specified position so as to introduce the scalar into the calculation domain. Custom scalar is typically defined as concentration in terms of mass per fluid volume in cell. The transport equation is solved automatically. The transport equation is as follows:where D is the diffusion coefficient and is a source term. The longitudinal/transverse/vertical diffusion coefficient can be obtained by the following formula:where is the fluid density, is the dynamic fluid viscosity in three directions, which is obtained by turbulence model, and is the Schmidt number, and according to previous studies [28], the value of a in compound open channel is generally 0.5–0.9, and the value of this study is 0.77. In this study, the scalar concentration is defined as 2000 μg/L. The location of the source is (0.3, 1.1, 0.13), and the flow rate at the source is 0.02 L/s.
The VOF method was used to compute the free water surface. The volume ratio of gas and fluid in the calculation grid can be expressed by functions and , respectively. The solution equation is as follows [29]:where is the volume fraction of liquid, .
The simulation area included two upstream inlet channels, an upstream pool, a downstream pool, three parallel channels, and a downstream outlet. The inlet boundary condition was the flow rate, the outlet boundary condition was the free outflow, the free surface boundary is set to a pressure inlet condition assuming atmospheric pressure, and other boundary conditions are wall. Figure 2 shows the meshing of three different densities. The first is all general mesh (0.002 × 0.002 × 0.002). The second is some general grids (0.002 × 0.002 × 0.002) and some dense grids (0.001 × 0.001 × 0.001), and the dense grid was used for the location near the wide crest weir, the position of the gradually expanding and shrinking section, and the plane “W-shaped” channel bend. The third is all dense mesh (0.001 × 0.001 × 0.001). In the calculation condition of Case 2, the three meshes were, respectively, calculated. It can be seen from Figure 3 that the calculation results of different mesh densities are roughly the same. However, the calculation values of the all general grid and the other two meshes are slightly different in the gradually expanding and shrinking position and the bend position. Therefore, the dense grid should be selected in the gradually expanding and shrinking position and the bend position. In order to improve the calculation speed of the CPU, the second type mesh is used in this study for numerical simulation.

(a)

(b)

(c)

(a)

(b)
2.3. Experiment Cases and Measuring Point Layout
Table 1 shows the flow rate divisions in three experimental cases. In the experiment, measuring points were arranged along the longitudinal center line of the channel with wide crest weir, gradually expanding and shrinking channel, and plane “W-shaped” channel. The number of measuring points was 19, 19, and 25, respectively, and the velocity and water depth of the section where the measuring point was located were measured.
3. Verification of Numerical Simulation
Figures 4 and 5 show the comparison between the measured and calculated values of water depth () and longitudinal velocity () along the center line of the three channels. It can be seen that the simulated water depth () and velocity () are in good agreement with the measured values. At the right-angle corner of the channel with a wide crest weir, the gradual expansion and contraction section of the gradually expanding and shrinking channel, and the bend part of the plane “W-shaped” channel, the flow pattern was complex, and the calculated values were slightly higher than the measured values. However, the relative errors of water depth and velocity were less than 5%.


4. Split Ratio of Parallel Channel
From the results of the numerical simulation and experimental measurements, it can be seen that the velocity from the distributary section to the confluence section of the three parallel channels was essentially the same—that is, the velocity head was also equal. Therefore, it can be said that the water level difference between the distributary section and the confluence section was mainly caused by the frictional head loss and the local head loss. The energy equation between the distributary section and the confluence section of the three parallel channels is as follows:where the corresponding section position of each variable is shown in Figure 6. is the average water level difference between the distributary section and confluence section. is the head loss of the channel with a wide crest weir. It should be emphasized that the average water level difference of the front and rear sections of the wide crest weir was taken as the calculation in this study. (i = 1–10) is the area of passage of the channel. (i = 1–10) is the frictional resistance coefficient. (i = 1–10) is the length of each channel. (i = 1–10) is the corresponding hydraulic radius of each channel. (j = 1–5) is the local resistance coefficient. are the flow rates of the three channels.

Let , , and where are the area ratio coefficient of the flow section. Therefore, Equation (9) to Equation (11) can be transformed into the following:
Let .
Then, Equation (12) to Equation (14) can be rewritten as follows:
Therefore, the calculation formula for the split ratio of parallel channels in a complex canal system is obtained as follows:
In the complex canal system, where the frictional head loss and the local head loss cannot be ignored, the size of the split ratio is related to the total head loss coefficient and the cross-sectional area. Table 2 shows a comparison of the split ratio between the three parallel channels, namely, the values calculated using Equations (17)–(19) and the numerical results. It can be seen that the split ratio of each channel solved using the two methods is essentially the same.
5. Flow-Field Structure and Its Influence on Pollutants Transport
5.1. Water Depth Variation
Figures 7(a), 7(b), and 7(c) show the water depth variation of the three parallel channels in Case 2. At the wide crest weir, the water surface profile rises and falls, which is mainly reflected in the water surface rising at the front of the weir, the water surface falling at the top of the weir, and the water surface rising behind the weir and at the right-angle corner after the weir. After passing through the weir, the water surface profile fluctuates. In the gradually expanding and shrinking channel, the water surface rises slightly before it enters the gradually expanding section. With an increase in the flow rate in this channel, the rising range decreases. When the water flow enters the gradually expanding section, the water depth decreases with an increase in the cross-sectional area. In the channel between the gradually expanding section and the gradually shrinking section, vortices appear on the water surface on the side wall, and the main stream in the middle is squeezed by the vortices on both sides. The water surface presents a “V” shape—high on both sides and low in the middle. In the gradually shrinking section, the water depth increases with a decrease in the cross-sectional area. In the plane “W-shaped” channel, the water depth of convex bank at the bend suddenly drops and the water depth of the convex bank is lower than that of the concave bank at each bend. The water level difference between the convex bank and concave bank decreases along the way.

(a)

(b)

(c)
5.2. Plane Velocity Distribution
In Case 2, the plane velocity distribution of the water level at z = 0.175 m is shown in Figure 8. The dead water zone is mainly distributed in the upstream and downstream pools, and there is reflux in the upstream pool, downstream pool, and downstream outlet.

In the channel with a wide crest weir, the wide crest weir reduces the cross section area of the channel, so the velocity increases at the top of the weir. The flow velocity increases near the right-angle bend, but the velocity in the corner outside the bend is almost zero. After the water flows through the right-angle bend, the velocity on the left is higher than that on the right because of the inertia effect of the main stream.
When the water flows into the gradually expanding and shrinking channel, the velocity increases sharply with the decrease in cross-sectional area. The water flow reaches the channel between the gradually expanding section and the gradually shrinking section, the velocity decreases gradually, and the main stream exists near the longitudinal center line of the channel. Therefore, the velocity near the center line is greater than that near the left and right banks. After the gradually shrinking section, the velocity increases once again.
In the plane “W-shaped” channel, the flow is affected by the plane shape of the bend. As a result of the imbalance between the lateral water surface gradient and the centrifugal force, the velocity of the convex bank increases and that of the concave bank decreases. At the bend, the mainstream is close to the convex bank. After the bend, the mainstream tends to the concave bank side. Owing to the inertial effect, the mainstream always approaches the concave bank side, and the flow velocity on the concave bank side is always greater than that on the convex bank side.
5.3. Velocity Distribution of a Typical Profile
Figure 9 shows the distribution of the flow field in the longitudinal section of the channel where the wide crested weir is located. In front of the wide crest weir, the velocity near the bottom and near the surface is very small, and the flow motion was irrotational [30]. The flow accelerated above the weir crest [31]. After the wide crest weir, the velocity in the upper part of the longitudinal profile is larger than that in the lower part of the profile, where the velocity distribution in the lower part being more complex, forming a vortex. In the three cases, the velocity increases with the increase in flow rate, the scale of the vortex also increases, and the vortex center is closer to the wide crested weir. The velocity distribution at the vortex increases gradually from the inside to the outside. In previous studies, it can be found that flow reattachment occurs behind the weir, and the reattachment length decreases with the increase in Reynolds number and finally reaches a fixed value [32, 33]. However, in this study, there is no flow reattachment behind the weir, which is mainly due to the lack of downstream distance. When the flow has not reached the reattachment distance, it is affected by the right-angle turning. Therefore, it can be concluded that the flow reattachment behind the weir is related to the downstream flow state.

(a)

(b)

(c)
Figure 10 shows the flow-field distribution in five typical cross-sections of the gradually expanding and shrinking channel. Under Case 2, in section CS1, there is transverse circulation on both sides. The circulation on the right is larger than that on the left, and the two circulation positions are close to the free water surface. This circulation trend brings water into the gradually expanding section, and the water begins to flow to both sides. In the CS2 section, owing to the influence of the circulation position and scale on both sides of the upstream, the diversion position location of the CS2 section is close to the left. The water flows into the channel between the gradually expanding section and the gradually shrinking section owing to the effect of the transverse shear force, and the velocity gradient of the water flow is generated at section CS3. Therefore, the shear force makes the flow form a small-scale circulation, the circulation position being located at the lower right corner of the CS3 section. The water flows into the gradually shrinking channel. With a decrease in the cross-sectional area, the flow on both sides is squeezed in the CS4 section. Due to the adjustment of the flow velocity gradient in the channel between the gradually expanding section and the gradually shrinking section, the flow velocity on the left and right sides is basically equal, so the extrusion center line is located near the channel center line [2]. Transverse circulation appears at the CS5 section. The circulation scale is large, and the circulation position is on the right side and close to the water surface.

(a)

(b)

(c)

(d)

(e)
Figure 11 shows the flow-field distribution of the three cross-sections (CS6, CS7, and CS8). At the bend, there is transverse circulation, and the circulation scale is large and close to the convex bank. After the bend, the flow always tends to the side of the concave bank in the transverse direction, and the flow field in the cross section becomes more stable with increasing distance along the way.

(a)

(b)

(c)
5.4. Pollution Incident Response
In Case 2, the diffusion of pollutants in sudden pollution accidents was simulated. The pollution event occurred at the left entrance channel of the complex canal system (see Figure 1). The discharge mode is a point source emission, the pollutant concentration was set at 2000 μg/L, and the pollutant discharge flow was constant at 0.02 L/s (0.37%). At time t0, the pollution incident began—that is, the pollutant discharged began. When the pollutant had discharged for 50 s, the discharge stopped (Figure 12(a)), and the complex canal system continued to provide a water supply for 100 s at a flow rate of 5.34 L/s.

(a)

(b)

(c)
Figures 12(b) and 10(c) show the variation of pollutant concentration with time at the inlet section (WCS1, WCS2, and WCS3) and outlet section (WCS4, WCS5, and WCS6) of the three parallel channels. According to Figure 6, the average velocity in the channel with wide crest weir and gradually expanding and shrinking channel was higher than that in the plane “W-shaped” channel, and most of the water flow in the left inlet channel where the pollution event occurred flowed to the channel with wide crest weir and gradually expanding and shrinking channel. Consequently, after the pollution event, it was evident that the amount of pollutants transported by the channel with wide crest weir and gradually expanding and shrinking channel was far greater than that of the plane “W-shaped” channel. The split ratio of the channel with a wide crest weir was greater than that of the gradually expanding and shrinking channel, as shown in Table 2, so that the channel with a wide crest weir could transport more pollutants than the gradually expanding and shrinking channel. After the end of the pollution event, when the pollutants had been discharged from each channel, the pollutants in the gradually expanding and shrinking channel discharged the fastest. Therefore, the split ratio, average velocity, and the location of sudden pollution accidents are important factors that affect the amount of pollutants transported in each channel.
Every table must have a descriptive title and, if numerical measurements are given, the units should be included in the column heading. Vertical rules should not be used (see Table 1). Tables should be cited consecutively in the text.
Figure 13(a) shows the concentration distribution of pollutants in six typical sections of the canal system at 30 s after the pollution event. The WCS7 section is located behind the wide crest weir. Due to the existence of large-scale vortices (i.e., dead water zone) [34], the diffusion is delayed and pollutants are trapped in the vortex zone. Transverse circulations are observed at the WCS9 and WCS10 sections (Figures 11(c) and 9(e)). The transverse circulation makes the pollutants fully mix with the water flow, so the diffusion of pollutants in WCS9 and WCS10 was also delayed. The WCS11 section and the WCS12 section are located in the plane “W-shaped” channel farthest from the pollution event, and the water flow of the channel comes from the right inlet channel (Figure 8), so the pollutant concentration here is at a minimum. However, it is also evident that the pollutant concentration at the bend (WCS11 section) is higher than that behind the bend (WCS12 section) because of the influence of transverse circulation at the bend. Therefore, it can be found that the dead water zone caused by vortex and circulation will delay the diffusion of pollutants, which is consistent with the results of Ouro [35].

(a)

(b)
Figure 13(b) shows the concentration distribution of pollutants in six typical sections of the canal system 40 s after the end of the pollution event. At this point, the pollutants originally retained behind the broad crest weir are carried out by the flow and gradually spread downstream. Therefore, the pollutant concentration of the WCS8 section is higher than that in the WCS7 section. Originally, the concentration of pollutants in the WCS9 section and WCS10 section is very large, and the concentrations of pollutants in the WCS11 and WCS12 sections are the smallest. However, the pollutant concentrations of these sections are very close at this point because the velocities of the WCS9 and WCS10 sections are higher than those of sections WCS11 and WCS12 (Figure 8). Therefore, it is evident that the pollutant emission rate is positively related to the velocity. At this point, the pollutant concentration of the WCS11 section and WCS12 section is greater than that at t0+30 s because there is stable transverse circulation in several cross-sections before the WCS11 and WCS12 sections (Figure 11), which retains pollutants and keeps them in the channel for a longer period of time. When the pollution event stops, the remaining pollutants gradually spread downstream.
6. Conclusions
(1)In a complex canal system where the frictional head loss and the local head loss cannot be ignored, the split ratio of three parallel channels can be calculated using the energy equation and continuous equation of constant total flow. The sum of head loss coefficients of each channel has an important influence on the split ratio of this type of the complex canal system.(2)The plane shape and flow-field structure of the complex canal system affect the pollutant diffusion. The stable transverse circulation on the cross section and the stable vortex on the vertical section can delay the diffusion of pollutants and increase the concentration of pollutants in the vortex and circulation position.(3)The calculated results are in good agreement with the experimental data, which indicates that the RNG model in Flow-3D software can be used to simulate the complex canal system. At the same time, the relationship between the flow-field structure and the pollutant diffusion can be used to control the pollutants in the complex canal system and provide some theoretical support for the river treatment after the sudden pollution accident.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China, grant no. 51679157, the National Key Research and Development Program, grant no. 2016YFC0401705, and the Sichuan Science and Technology Program, grant no. 2019JDTD0007.