Next Article in Journal
Evaluating the Influence of Deficit Irrigation on Fruit Yield and Quality Indices of Tomatoes Grown in Sandy Loam and Silty Loam Soils
Previous Article in Journal
Drivers of Turbidity and Its Seasonal Variability at Herschel Island Qikiqtaruk (Western Canadian Arctic)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Study on the Shear Stress Characteristics of Open-Channel Flow over Rough Beds

1
State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China
2
School of Civil Engineering and Geomatics, Southwest Petroleum University, Chengdu 610500, China
3
State Key Laboratory of Urban Water Resource and Environment, Harbin Institute of Technology, Harbin 150090, China
4
Shenzhen International Graduate School, Tsinghua University, Shenzhen 518055, China
*
Author to whom correspondence should be addressed.
Water 2022, 14(11), 1752; https://doi.org/10.3390/w14111752
Submission received: 4 May 2022 / Revised: 24 May 2022 / Accepted: 26 May 2022 / Published: 30 May 2022
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
Bed shear stress is an important measure of benthic habitats since it is related to many ecological processes. In this study, we focused on the fluctuating characteristics of shear stress in rough-bed open-channel flows. The roughness element method was adopted to mimic natural rough beds and the Improved Delayed Detached Eddy Simulation (IDDES) model was used to obtain comprehensive information about shear stress near the rough bed. Three arrangement patterns of the roughness elements were simulated to compare their effects on flow structure and shear stresses. The arrangements of the roughness elements altered the Reynold stress and turbulent kinetic energy characteristics, due to the variance of blockage in lateral directions that led to flow detachment and changes in the flow directions. Quadrant analysis revealed the spatial variations of the instantaneous shear stress burst events at different locations in the wake. By using spectrum analysis, the accumulation of shear-stress energy from small to large vortex scales was estimated, which revealed that the instantaneous effect of the shear stress was significantly stronger than the effect of the time-averaged shear stress, especially on small scales. The results of this study suggest the significance of the fluctuation part of shear stress in further studies on ecological processes.

1. Introduction

The construction and operation of dams and hydropower plants dramatically change river hydrology and hydrodynamics, involving a series of hydrodynamic variables such as flow velocity, depth, and shear stress [1,2]. Among these variables, shear stress is an important measure of benthic habitats since it plays a critical role in many ecological processes including macroinvertebrate drifting [3] and bio-particle transportation [4]. An increase in shear stress destabilizes benthic macroinvertebrate habitats by moving sediments or directly forcing benthic macroinvertebrates to leave their habitats, which both lead to spatial variations in the benthic macroinvertebrate community [5,6,7]. Additionally, shear stress is widely regarded as the key indicator of bed material movements, which makes it significant in studies of river morphology.
In natural rivers, shear stress is a fluctuating signal and its related ecological processes may not only depend on its time-averaged value (usually regarded as the Reynolds stress) but also on its fluctuations. Some studies have highlighted the significance of the instantaneous fluctuation part of shear stress in turbulence [8,9,10]. However, the fluctuation part of shear stress has been seldom studied, compared to the time-averaged shear stress studied widely [5,11,12]. Known as ‘burst events’, shear stress fluctuations can be divided into four event types by quadrant analysis, among which ejections and sweeps are dominant. The fluctuation part of shear stress generates an additional force on the particles and thus reduces the shear stress threshold of dislodgement [9,10,13]. Blanckaert et al. revealed that the temporal fluctuation of flow velocity results in an increase in drag force by about 100% [9]. Turbulence fluctuation is related to the development of vortex structures, and the magnitude of fluctuation is significantly affected by vortex scales. It is indicated that the magnitude of shear stress may vary with the vortex scale, although the spectral analysis of turbulence is rarely adapted to shear stress.
The direct measurement of bed shear stress is a difficult task, and its spatial and temporal characteristics are usually analyzed based on near-bed hydrodynamic conditions. Flow field data can be acquired by direct measurements, such as Particle Image Velocimetry (PIV), Acoustic Doppler Velocity Profiler (ADVP), and intrusive Electromagnetic Current Meters (ECMs) [8,9,14,15]. However, these high-resolution measurements that satisfy the analysis of local shear stresses have a high cost and the results are also affected by the wall boundary or measurement instrument, and by now near-bed flow structures are still poorly investigated [16]. An alternative method is numerical simulation, by which the disadvantages in measurements could be overcome, and the comprehensive flow field data are available. Currently, in the ecohydraulic field, the most commonly used hydrodynamic models (e.g., PHABSIM, RHYHABISIM, and River-2D) only involve time- and space-averaged flow variables [17,18,19,20]. Morris et al. tried a three-dimensional Large Eddy Simulation (LES) in their study on Glossosoma biomass distribution, based on a high-resolution riverbed morphology measurement [21]. This provides insight into the application of three-dimensional high-resolution flow simulations in the ecohydraulic field, but this costly case study cannot be used widely. Thus, an adequately simplified method of rough channel beds is needed to balance computation costs and accuracy.
In numerical simulations, near-bed flows are commonly modeled by three methods, namely wall functions, wall with roughness elements, and wall applying measured rough bed morphology [21,22,23]. Wall functions over-simplify the rough bed and are unable to reflect the effect of large-size sediments on the local flow field. Additionally, as mentioned above, rough bed morphology measurements are costly. Roughness elements can be an appropriate tool to capture the flow characteristics of the roughness of a natural riverbed [16]. Among the many types of roughness elements, two-dimensional bars and three-dimensional cubes are most frequently used. Previous numerical studies have revealed the influence of the roughness element dimension, spacing ratios between the bar interval p and roughness element height d, cross-section shapes, and arrangement patterns upon turbulence characteristics such as turbulence intensities and drag coefficients [22,23,24,25,26]. In comparison, numerical studies and statistical analyses on shear stress characteristics, especially the fluctuation part of shear stress, are seldom found.
The main purpose of this study is to evaluate the shear stress characteristics in turbulent flow, especially the instantaneous characteristics, based on numerical results by using the roughness element method to mimic riverbeds. We first introduce the turbulent model applied in this study and test its validity, and then describe the roughness element arrangements of the conducted simulations. Following are the results of the numerical simulations analyzed by quadrant analysis and spectrum analysis. Discussions about the shear stress fluctuation events’ spatial variation and the accumulated effect of shear stress at different scales are also included in this section. Finally, the concluding remarks are presented.

2. Materials and Methods

2.1. Turbulence Model and Validation

2.1.1. Improved Delayed Detached Eddy Simulation Model

This study used the Improved Delayed Detached Eddy Simulation (IDDES) model to investigate the turbulent characteristics of near-bed flow, considering the balance between calculation cost and accuracy. The IDDES model is a Reynolds-averaged Navier–Stokes (RANS)/Large Eddy Simulation (LES) hybrid model based on the Detached Eddy Simulation (DES) model, whose feasibility in simulating turbulent fluctuations has been proved by previous research [27]. The model invokes either the RANS model or LES model in different flow regions according to the relative size of the turbulent length scale and grid scale. The regions in which the turbulent length scale exceeds the grid scale are solved by the LES model, and the regions near the solid boundaries and the locations in which the turbulent length scale is smaller than the maximum grid scale are solved by the RANS model. The governing equations of the Shear Stress Transport (SST) k ω IDDES model are shown in Equations (1) and (2):
ρ k t + · ( ρ U k ) = · [ ( μ + σ k μ t ) k ] + P k ρ k 3 / l I D D E S
ρ ω t + · ( ρ U ω ) = · [ ( μ + σ ω μ t ) ω ] + 2 ( 1 F 1 ) ρ σ ω 2 k ω ω + α ρ μ t P k β ρ ω 2
in which ρ is the fluid density, k is the turbulent kinetic energy, ω is the specific turbulent dissipation rate, U is the flow velocity vector, μ is the fluid viscosity, μ t is the turbulent eddy viscosity, Pk is the rate of production of k, and F1 is the blending function for α and β, α = α 1 · F 1 + α 2 · ( 1 F 1 ) , in which α 1 = 5 / 9 and α 2 = 0.44 , and β = β 1 · F 1 + β 2 · ( 1 F 1 ) , in which β 1 = 0.075 and β 2 = 0.828 . l I D D E S is the IDDES length scale. σk = 2. σω2 = 0.856 [28].
The original DES approach will incorrectly invoke the LES model when the grid is locally refined in multiple directions, especially in boundary regions. As a result, the eddy viscosity may be substantially underestimated from the RANS region to the LES region, triggering Grid Induced Separation. In the IDDES model, intricate blending and shielding functions are applied in characteristic length scales to protect the wall boundary layer from the incorrect invocation of the RANS model [28]. Equation (3) shows the length scale used in this model:
l I D D E S = f d ˇ ( 1 + f e ) l R A N S + ( 1 f d ˇ ) l L E S  
in which l R A N S , the RANS turbulence length scale, is calculated from the turbulence kinetic energy k and the specific rate of dissipation ω for the k ω model, i.e.,
l R A N S = k ˜ C μ ω
l L E S is the LES filter width, defined as:
l L E S = C D E S · min { C w max [ d w , Δ max ] ,   Δ max }
with C D E S = 0.65 , C w = 0.15   , and Δ max = max ( Δ x , Δ y , Δ z ) . f d ˇ is the empirical blending function that controls the selection of the characteristic length scale (and consequently the model applied in the grid) and f e is the elevating function that prevents an excessive reduction in the eddy viscosity in the vicinity of the RANS/LES interface. The details of f d ˇ and f e can be found in Gritskevich et al., 2012. [28].
When the near-wall mesh is fine enough, the region using RANS would be very limited; namely, the majority of the domain would be calculated by LES. In this study, we adopted the SST–IDDES model and the RANS region thickness was about 0.3 mm, which was thin enough to be neglected compared with the roughness height and the total flow depth.
The simulations were conducted by using Ansys Fluent 14.0, where the governing equations are numerically solved based on the finite volume method and PISO scheme.

2.1.2. Model Validation

The model was validated by predicting a pressure-driven flow in a channel with roughness elements and comparing the prediction against the experimental results found by Krogstad, et al. [29].
The plane channel in Krogstad’s experiment was 5 m long, 1.35 m wide, and 0.1 m high, with the bottom and top walls equipped with transverse square rods. The height of the rods was d = 1.7 mm and the distance between the front face of two adjacent rods was p = 8d. The cross-section average velocity was 0.1208 m/s and the Reynolds number based on the shear velocity u * was R e * = ρ u * h / μ = 600 , in which h was the half-height of the channel.
The bottom half of the plane channel with 15 rods was simulated in the present study, and the computational domain was Lx = 0.204 m long and Ly = 0.05 m high, with the roughness elements arranged on the bottom surface. The symmetry boundary condition was specified on the top boundary, and the periodic boundary condition was applied on the upstream and downstream boundaries. The width of the domain was Lz = 0.2 m and the periodic boundary condition was applied to the boundaries in the lateral direction.
The calculation domain and grid of the numerical model are shown in Figure 1. The computational domain was divided into 3 layers of different resolutions in the vertical direction in order to reduce the total grid amount while maintaining enough grid resolution in the near-wall region and a good aspect ratio in the domain. The layers were connected by the interface boundary condition. The total grid amount was 789,000. The grid size in the layer near the rough floor was 0.34 × 0.34 × 0.34 mm (Δy+ ≈ 4.1) while the grid size near the upper boundary was 1.36 × 3 × 3 mm. In the vertical direction, the cell height increased with a growth rate of 1.2. x, y, and z were the longitudinal, vertical, and lateral directions. u, v, and w were, respectively, the velocity components in these three directions, and u′, v′, and w′ were the fluctuation parts of the velocities. u′v′ represented the flow shear stress exerted by the velocity fluctuations, and its time-averaged value multiplied by the flow density ρ is the Reynolds shear stress (usually represented by ρ u v ¯ ).
Previous studies have concluded from the power spectrum that the major part of the water flow frequency ranges from 0–100 Hz. Considering the Fourier transform used in this research to analyze the turbulence characteristics, the sampling frequency should be no less than 200 Hz [30,31]. In this study, the time step and sampling interval were set to be 0.005 s, and the total sample size N = 6000 [30].

2.1.3. Validation Results

The normalized velocity, turbulence intensities, and Reynolds stress of the numerical results were compared with Krogstad et al.’s measurement for model validation [29]. Their profiles were averaged horizontally over the computational domain. Figure 2 shows the velocity profiles which were normalized by the shear velocity u * .
Based on:
u * 2 ( 1 y H ) = u v ¯ + ν u y
in open-channel turbulence, in which H represents the flow depth (and equals to h here), as u v ¯ ν u y in the outer region, u * can be calculated from the intercept of the linear part of u v ¯   ~   y at y = 0. The non-dimensional quantities are marked by the ‘+’ superscript in this paper. The velocity’s relative error of the IDDES model was calculated on each point on the profile, and their average was 5.0%. The normalized turbulence intensities and Reynolds shear stress of the IDDES model and Krogstad, et al.’s measurement are shown in Figure 3, where the turbulence intensities are represented by the root mean square error of the flow velocities and are marked by the ‘rmse’ subscript. The average relative errors of the turbulence intensity profiles in the x, y, and z directions were 9.5%, 9.0%, and 4.3%, respectively. As for the Reynolds stress profile, the absolute error was calculated, because the normalized Reynolds stress value approached zero on the top boundary. The absolute average error of the normalized Reynolds stress was 0.151, which was especially higher in the near-wall region where the numerical results exceeded the measurements. This was also observed in Krogstad, et al.’s direct numerical simulation (DNS) results, as the intrusive hot-wire method applied in the research caused errors in the vertical velocity measurements in the near-wall region, which tended to underestimate the measured Reynolds stress peak compared with the DNS. In other words, the real peak of the normalized shear stress should have been higher than the measured results in Krogstad, et al. [29]. For y/h > 0.2, the average error of the numerical model was reduced to 0.028, indicating that the agreement of the IDDES model was satisfactory.
It can also be seen from Figure 3 that the specification of the symmetry boundary on the upper boundary will decrease the vertical turbulent intensity and increase the lateral turbulent intensity. However, the influence is limited in the region of y/h > 0.9. In the region far from the upper boundary, especially in the region near the bottom floor (where we are mostly concerned), the influence of the symmetry boundary is neglectable. Therefore, it can be assumed that, in this study, the symmetry boundary could be used to represent the free water surface.
For a uniform flow, the usage of periodic boundary conditions on upstream and downstream boundaries is valid for time-averaged velocities but not for velocity fluctuations. To evaluate the influence of the periodic boundary condition, spatial correlation coefficients between the inlet and the downstream points on the centerline at y = 1 mm were calculated and are shown in Figure 4. It was observed that the use of the periodic boundary condition only affected a small range near the boundaries, indicating that in the near-wall region, the longitudinal correlation of the velocity fluctuation was weak, and that the periodic boundary condition had a minor effect on the near-wall flow field. This enables the application of the periodic boundary condition for the purpose of reducing the total grid amount and improving the efficiency of the simulation.

2.2. Simulation Domain Description and Rough Wall Configuration

In Krogstad et al.’s case, the half-channel height was h = 0.05 m and the roughness height d = 1.7 mm, which is much lower than the flow depth and roughness height caused by substrates in natural rivers, respectively. Therefore, the following simulations were based on Blanckaert et al.’s experiment to approximate a turbulent flow over a natural rough bed [9].
The open channel used in Blanckaert et al.’s experiments were 8.5 m long and 0.5 m wide. We adopted one test from their experiments, in which the flow depth was 0.2 m and the average flow velocity was 0.3 m/s. The bottom sediment characteristic sizes were D50 = 0.8 mm, Dm = 2.3 mm, and D90 = 5.7 mm [9]. In the following simulations, the computational domain was Lx = 0.24 m long, Ly = 0.2 m high, and Lz = 0.24 m wide. The boundary conditions were the same as those in the validation case. The three-layer layout of the mesh was also adopted, with the total grid amount ranging from 3,474,600 to 6,568,800 depending on specific cases. The grid size near the rough floor was about 0.5 mm × 0.5 mm × 0.5 mm, while the grid size near the water surface was 2 mm × 4 mm × 2 mm. The time step Δt and sample size were identical to the validation case. Sampling started after a simulation duration of 120 s when turbulence was fully developed in the domain.
The roughness elements were arranged on the bottom surface to mimic rough beds. Three scenarios with different roughness elements were investigated in the present study. They varied in terms of dimension, shape, spacing, and arrangements:
Case 1: two-dimensional bars with square cross-section shapes were perpendicularly arranged to the flow direction (Figure 5a).
Case 2: three-dimensional roughness elements of cubes were arranged in chessboard formation (Figure 5b).
Case 3: three-dimensional cubes were staggeringly arranged (Figure 5c). The spacing between the roughness elements in the longitudinal and lateral directions were larger than that of the chessboard formation, letting the fluid flow by their sides. The main difference between Case 2 and 3 was that in Case 2, no interstice that flow could go through existed between the cube elements connected to each other at the corners, while in Case 3, the spaces between the cubes were wide enough to let the flow go around the cubes.
We used these models to mimic the natural riverbeds where macroinvertebrates usually shelter themselves among substrates; therefore, the interstices and roughness element height were determined according to the benthic macroinvertebrate body size (1 mm to several centimeters). In the present study, the interstices were set to be 30 mm (except for the 3D-cu-s, in which the interstices were intentionally expanded). The roughness elements height was set to be d = 8 mm. Note that, by definition, the distance between the roughness elements p refer to the distance of a roughness element from its successive one in the x or z direction, and p varies according to the roughness element arrangements in different cases. Figure 5 labels the definition of p in each case. Specific magnitudes of roughness element size and spacing, as well as their arrangement patterns, are listed in Table 1. These cases were named by their dimension, roughness element shape, and arrangement patterns, and are also listed in Table 1.

2.3. Turbulence Analysis Methods

2.3.1. Statistical Variables for Near-Bed Flow

The main characteristic hydrodynamic variables include shear velocity u * , the roughness function Δ U + , and the roughness equivalent height ks. In this study, the shear velocity u * was calculated using Equation (6). The roughness function Δ U + was calculated according to the log profile of the time-averaged longitudinal flow velocity:
U + = 1 κ ln ( y + ) + B Δ U +
where the Karman constant κ = 0.41 , B = 5.5. On rough surfaces, the theoretical datum y = y0 should be considered for the application of the log profile of flow velocity. In the following part of this study, we redefine y+ = 0 to be acquired at the theoretical datum y = y0, giving:
U + = 1 κ ln ( y y 0 ν / u * ) + B Δ U +
and y + = y y 0 ν / u * . ks is derived by:
Δ U + = 1 κ ln ( k s + ) 3.5
In common practice, ks/d is used for analysis.

2.3.2. Quadrant Analysis Methodology

Quadrant analysis was first conducted by Wallace et al. to explore the unknown information on Reynolds shear stresses from velocity fluctuations [32]. The principle of this method is to classify shear stress u′v′ into four categories by instantaneous turbulent velocity fluctuations (u′, v′): Q1 (+u′, +v′), Q2 (−u′, +v′), Q3 (−u′, −v′), and Q4 (+u′, −v′), each of which represents an event of the turbulent burst phenomenon. In this way, shear stress characteristics including magnitude, contributions of each event, and velocity fluctuation components can be investigated by quadrants [9,33].
In this study, the velocity fluctuations were locally affected by the wakes behind the different roughness elements. At different locations of the wakes, the instantaneous shear stress −u′v′ was plotted by its u′ and v′ components on the u′-v′ plane to visualize the shear stress characteristics, such as the magnitude of shear stress, its velocity fluctuation components, and the quadrant it is in. The magnitude of the shear stress of each point is represented by the product of its coordinate values, so points with the same absolute magnitude of shear stress are located at the same | u v | = c o n s t hyperbolas.
Previous studies suggested a spatial variation in the time-averaged shear stress on the horizontal plane near the rough floor, and that its distribution depends on the pattern of the roughness element arrangements [16]. It is necessary to classify the locations around the roughness elements and to conduct an analysis by different locations. Considering the periodic arranging patterns of the roughness elements in the present study, the locations over the roughness elements and the locations behind the roughness elements were selected since they were, respectively, under the control of the main flow and the recirculation flow in the cavities. In the quadrant analysis of the instantaneous shear stress, we chose locations above each roughness element at y = 10 mm and 9 mm behind the front edge of the roughness element (type A) and locations 18 mm behind the back edge of each roughness element in the cavities at y = 6 mm (type B). This selection is shown in Figure 6 and would be applied in all cases for further comparison. The instantaneous shear stress in the sampled 6000 steps was analyzed.
In common quadrant analysis practices, longitudinal and vertical velocity fluctuations are considered as the studied flow is mainly longitudinal and lateral flow velocity could be neglected. In this study, for the type B locations, the roughness elements led to a strong three-dimensional flow in the roughness element layer, so the longitudinal and lateral velocity could be combined into a horizontal flow velocity.

2.3.3. Frequency Spectrum Analysis for Shear Stress

In general, spectrum analysis has been widely used to find the dominant frequency, focusing on the intensity of a vortex of a certain frequency (or certain length scale). Commonly analyzed variables in spectrum analysis include velocity, pressure, and turbulent kinetic energy [34,35].
Spectrum analysis was performed by using Fourier transform, by which a temporal signal A(t) can be transformed into a combination of trigonometric signals:
A ( t ) = a 0 2 + n = 1 [ a n cos ( 2 π n t T ) + b n sin ( 2 π n t T ) ]
In Equation (10), the constant term (the zeroth-degree harmonic term) represents the time average of the signal, and the trigonometric terms represent the signal fluctuation around the average, due to vortices of all sizes. This separates the effects of vortices of different frequencies, i.e., of different length scales.
Fourier transform is usually applied to the instantaneous velocity u ( t ) and the total energy of this Fourier series corresponds to the total energy of the kinetic energy of the flow. According to the Parseval theorem, this energy is the integral of the energy contained in the signal of each frequency and could be calculated from the amplitude of each harmonic term:
1 2 π π π | u ( t ) | 2 d t = a 0 2 4 + 1 2 n = 1 ( a n 2 + b n 2 )
In the Fourier series of instantaneous velocity, a 0 / 2 = u ¯ , n = 1 ( a n 2 + b n 2 ) = u 2 ¯ , therefore the integral should be u ¯ 2 + u 2 ¯ / 2 (the symbols with an overbar represent the time-averaged values). This integral is the sum of the effects of all sizes of vortices.
In this study, we followed this idea by analogizing the fluctuating part of shear stress to the fluctuating velocity and accumulating the effects of vortices from low frequency to high frequency on shear stress. In this way, we discovered the differences between the same turbulent shear stress acting upon particles such as benthic macroinvertebrates of different sizes. Similarly, Fourier transform was applied to the instantaneous shear stress signal and its energy was analyzed. The accumulation of the shear stress signal energy from low frequency to high frequency could be written as a function of frequency:
F ( f ) = { a 0 2 4 ( f = 0 ) a 0 2 4 + 1 2 n = 1 f ( a n 2 + b n 2 ) ( f > 0 )
where f is the frequency of the shear stress signal, a 0 2 / 4 = ( u v ¯ ) 2 , and n = 1 f ( a n 2 + b n 2 ) = ( u v u v ¯ ) 2 ¯ when f→∞. Compared to the frequency, the vortex length scale is a more practical indicator when relating hydrodynamic characteristics to specific objects in the flow. It is natural to transfer the frequency value into a vortex length scale value by adopting the local time-averaged flow velocity uc as the characteristic velocity:
L v = u c / f
An instantaneous shear stress signal of type B was selected for spectrum analysis to reflect the shear stress characteristics in the roughness element layer.

3. Results and Discussion

3.1. Characteristics of Time-Averaged Hydrodynamic Parameters

Figure 7 displays the predicted and measured velocity profiles, normalized by the shear velocity u * (calculated by Equation (6) and listed in Table 2). These profiles have approximately equal slopes and their different intercepts separate them from each other. The roughness functions Δ U + were obtained by subtracting the intercepts from 5.5 and are listed in the third row of Table 2. Then, the roughness equivalent heights ks/d were derived by using Equation (9) and are listed in the fourth row of Table 2.
With the same flow rate, the roughness element arrangements had a great impact on the flow resistances. According to the values of ks/d and Δ U + , 2D-sq generated the largest flow resistance, followed by 3D-cu-c and 3D-cu-s, which corresponds to Volino et al.’s conclusion [23]. These cases are listed in the same descending order by the shear velocity u * . Among the studied cases, the flow resistance in 3D-cu-s was the closest to the experiment of Blanckaert et al., while the 2D-sq and 3D-cu-c profiles in Figure 7 were approximately coincident. It is not surprising that 3D-cu-s had the smallest flow resistance since its roughness elements were arranged in the sparsest way. In the region y + < 100 , these velocity profiles deviated from the logarithmic law. Compared to the logarithmic law, the velocity-profile slopes for 2D-sq and 3D-cu-c in the y + < 100 region increased while that for 3D-cu-s decreased, indicating different vortex structures among these cases. For the case of 3D-cu-s, the horizontal vortices generated by the flow around the roughness elements were responsible for the nearly uniform velocity profile in the region y + < 100 .
Figure 8 displays the normalized Reynolds shear stress profiles and Figure 9 shows the normalized turbulent kinetic energy profiles, both of which were averaged horizontally over the computational domain. A sharp increase in the Reynolds shear stress in 2D-sq and 3D-cu-c was observed in the near-wall region, as the vertical recirculation flow in the cavities between the roughness elements produced a strong shear. A jump in the Reynolds stress profile was found at the top face of the roughness elements, above which the form drag due to the roughness element vanishment. Above the roughness element layer, the three-dimensional cases were roughly 20~30% lower than the two-dimensional case in the y+ = 800~2600 range, which was noted by Volino et al., indicating the influence of the different roughness dimensions. There was a local peak for the Reynolds shear stress above the roughness elements, respectively, at y+ = 237, y+ = 86, and y+ = 142 in 2D-sq, 3D-cu-c, and 3D-cu-s. The increasing and decreasing trends in the turbulent kinetic energy between the three cases were relatively similar, while their turbulent kinetic energy was basically dependent on the value of u * . However, the peak in 2D-sq was greater than in the two three-dimensional cases, indicating the difference in the flow structures among the three cases.
The influence of the rough wall to the time-averaged shear stress u v ¯ could be partially reflected by its components u′ and v′, while urmse, vrmse, and wrmse are evaluations of the magnitude of the velocity fluctuations. Figure 10 displays the urmse, vrmse, and wrmse profiles, normalized by u * . For the open-channel rough turbulence, urmse+ > wrmse+ > vrmse+ was observed [22,36]. The roughness elements exerted a greater impact on the velocity fluctuation in the near-wall region. Below y+ = 100 in the roughness element layer, 2D-sq and 3D-cu-c had a similar urmse+ and vrmse+, while 2D-sq and 3D-cu-s had a similar wrmse+. This similarity is related to the overlapping of the normalized Reynolds stress in this region for 2D-sq and 3D-cu-c. In the outer region, the urmse+ differences between the cases were similar to that in the Reynolds stress profiles, with the three-dimensional cases around 12~22% lower while the vrmse+ values were almost identical. This indicates that the differences in the Reynold stress might have mainly resulted from the differences in the longitudinal velocity fluctuation characteristics.
Volino et al. suggested that the complete obstruction over the channel width in two-dimensional cases forces flow detachment and acceleration over the bars and leads to a larger ks/d and vertical range affected by the rough wall [23,37]. The complete blockage in the lateral direction forces the flow to go over the roughness bars, offsetting the position of the maximum time-averaged shear stress away from the wall [23]. This phenomenon was also observed in this study, as 2D-sq had the largest u * , ks/d, and urmse+. Comparatively, a horizontal flow around the roughness elements was allowed for the typical three-dimensional 3D-cu-s case, where uniform distributions of the horizontal velocity (both time-averaged values and their fluctuations) were found in the element roughness layer. With regard to 3D-cu-c, the top surface drag of the roughness element forced its upstream flow to its two side cavities by going over the corner connections, although the chessboard-arranged roughness elements let no interstices between them. As a result, the flow in the roughness element layer was similar to 2D-sq, while the flow over the roughness element was similar to 3D-cu-s.
The overlapping of the urmse+, vrmse+, and wrmse+ profiles in the roughness element layer was another consequence of the roughness elements’ blockage effect in each direction. In 2D-sq and 3D-cu-c, the arrangements of the roughness elements cut off the continuous longitudinal flow in the roughness element layer, which limited the longitudinal velocity fluctuations. In contrast, the arrangement of the roughness elements in 3D-cu-s allowed fewer limitations to both the longitudinal and lateral flow, resulting in a higher urmse+ than 2D-sq and 3D-cu-c and a similar wrmse+ as 2D-sq in the roughness element layer, where flow in the lateral direction was also less limited.

3.2. Quadrant Analysis for Shear Stress

The scattered points on the u′-v′ plane representing the type A and B locations each formed an elliptic cloud, and based on the number of the points that formed the cloud, the joint probability density function (PDF) of u′ and v′ can describe the probability distribution of fluctuations in each quadrant [38]. Here, 2D-sq is taken as an example and the joint PDF contours were calculated and are displayed in Figure 11. The major axes of the elliptic clouds and hyperbolas of the average shear stress are also plotted. On the u′-v′ plane, the points of shear stress were mainly distributed in Q2 (36.38% in type A and 32.31% in type B) and Q4 (24.79% in type A and 36.05% in type B), with a stronger magnitude than the Q1 and Q3 fluctuations, indicating the dominance of the ejection and sweep event. The fluctuations in these two quadrants had positive contributions to the shear stress –u′v′, and for all three cases, the Q2 and Q4 fluctuations dominated and their total proportion in the type B points was higher than in the type A points, resulting in a higher average shear stress magnitude in the type B points. Table 3 provides the proportion of each quadrant in all three cases. In 2D-sq, the Q2 proportion was higher in the type A points, which is typical at the boundary layer edge from the front edge of the bars [38]. In the 3D cases, the Q4 proportions were higher in the type A points, as the arrangements of the staggeringly distributed cubes allowed velocity exchange between the type A points and the points beside them in the lateral direction with a higher-speed flow, which contributed to the Q4 fluctuations. For the type A and type B points, the two-dimensional case and three-dimensional cases had opposite top-dominating fluctuation events, which may explain the different profiles of the Reynolds shear stress between the two-dimensional case and three-dimensional cases. For the type A points in 2D-sq, the top-dominating event of the Q2 transported the larger shear stress upward and resulted in a smaller slope above the local peak point (Figure 8).
The fluctuations in Q2 and Q4 represent the exchange of energy in vertical directions and the accompanying mass transportation, such as the dislodgement of sediment and benthic macroinvertebrates. Previous research has revealed that Q2 fluctuations contribute to sediment dislodgement while Q4 fluctuations are related to mass transportations toward the bed. Q2 fluctuations contribute to particle dislodgement from the bed, which is partly related to the positive u’ and thus the greater longitudinal velocity, and Q4 fluctuations maintain their suspension in the water column [8,9]. The spatial variation of the Q2 and Q4 events corresponds with previous research on microflow regimes around stream boulders, and is related to benthic macroinvertebrate distributions [39]. While shear stress is commonly regarded as a driving force of particle dislodgement from riverbeds, it may also benefit the transportation of particulate organic matter and the exchange of dissolved gases, which is another factor that relates to benthic macroinvertebrate distributions. Further analysis by decomposing shear stress components can build links to Stokes equations, converting shear stress components into lift force [9]. This builds a bridge between the flow shear stress and the actual force exerted on the particles, which allows an intensive understanding of sediment and benthic macroinvertebrate movements from a mechanical perspective.

3.3. Spectrum Analysis for Shear Stress

Figure 12 displays the power spectrum of the shear stress temporal signal. The major part of the shear stress energy is contained in low-frequency fluctuations, and a −5/3-law can be observed in the power spectrum, which is similar to turbulent kinetic energy and scalar transportation [34,40,41,42]. The frequency range of the −5/3-law is roughly from 8 Hz to 30 Hz. Taking the local time-averaged velocity as the characteristic velocity, the corresponding length scale range is about from 1 mm to 7 mm. The largest length scale has the same order of magnitude as the roughness element size and the smallest length scale has the same order of magnitude as the Kolmogorov length scale.
The power function F(f) added from low frequency to high frequency represents the accumulation of the shear stress energy from large scale vortices to small-scale vortices, and its maximum is the total energy of the shear stress temporal signal. In Figure 13, the frequency f is replaced with characteristic length scale   L v , and the curves represent the accumulation of the shear stress signal energy from the largest scale to the smallest scale. The integral increases with the decrease in length scale, and it reaches the maximum at around 1 mm, corresponding to a frequency of around 30 Hz where the power spectral density is about two orders of magnitude smaller than its maximum.
The total energy of the shear stress temporal signal was about 261.0~375.8% of the energy of time-averaged shear stress. According to the deduction of Equation (12), the theoretical value of this ratio should be 1 + ( u v u v ¯ ) 2 ¯ / 2 ( u v ¯ ) 2 , which is decided by the mean and the temporal fluctuation magnitude of the shear stress. The stronger the fluctuation (relative to the time-averaged shear stress), the larger the ratio. The power function revealed that the shear stress energy was mainly contained in fluctuations of 0~30 Hz, which corresponds to the energy cascade. For 3D-cu-s, the sparse distribution of the roughness elements resulted in low resistance to the flow, reflected by the lowest ks/d and u * . The rough wall in this three-dimensional case affected a relatively smaller range in the vertical direction (Figure 10) [23]; therefore, it had a weaker influence on the large-scale fluctuations, leading to the smallest normalized shear stress energy in 3D-cu-s.
In turbulent flows, large-size vortices contain the majority of the turbulent kinetic energy and are dominated by force from inertial range scales, exerting a pressure gradient force upon the particles in the fluid [43]. Small-size vortices dissipate the turbulent kinetic energy through viscosity, and the viscosity force is proportional to the torque upon the particles [44]. The available literature pointed out that for a finite-sized particle, as the Particle Reynolds number increases along with the particle size, the effect of small-size vortices decreases [43]. Experiments by Qureshi et al. revealed that the transport and rotation of particles are forced only by turbulent pressure fluctuations at scales larger than the particle [44].
By extending this argument to shear stress fluctuations, the accumulation of shear stress energy should be conducted from the largest length scale to the size of the object. The total effect of the flow fluctuations exerted upon the object decreases with the increase in its length scale. At the range of the benthic macroinvertebrate size (mm~cm), smaller benthic macroinvertebrates (~2 mm) are usually affected by all scales of vortices; therefore, small-scale flow fluctuations should not be neglected. For benthic macroinvertebrates and other particles of larger sizes (>2 mm), the effect of small vortices weakens with a scale increase; therefore, the effect from shear stress fluctuations is reduced. For larger objects in the flow, the effect of the shear stress fluctuation is even weaker, and finally the total effect of the shear stress is reduced to merely the effect of the time-averaged Reynolds stress which, together with the time-averaged flow velocity, has the main contribution.
When taking the difference of particle mass into consideration, fluctuations in turbulent flows may exert a stronger influence on small-size particle movement than expected. A previous experiment estimated that instantaneous velocity could be 40% higher than its time average, resulting in an increase of about 100% in lift force [9]. In this study, we evaluated the total energy of the shear stress signal, which could be 2~4 times the time-averaged part. This evaluation was based on the current flow field simulation, while the diversity of the turbulent flows and the complexity of fluid–solid coupled dynamics may result in different phenomena.

4. Conclusions

Both the mean value and the fluctuations of turbulent flow shear stress play critical roles in river ecological processes, but the fluctuation of shear stress has been seldom studied because of difficulties in direct measurements. In this study, we used numerical simulations to study shear stress characteristics over a natural rough riverbed, including the time-averaged and fluctuating characteristics. The roughness element method was used to mimic a natural rough bed. The arrangements of the roughness elements resulted in different degrees of blockage over the channel width by changing the flow pattern in the roughness element layer. Quadrant analysis revealed that ejections and sweep events were stronger behind the roughness element, which indicates a strong momentum exchange. This corresponds to the benthic macroinvertebrate distribution around boulders in the river. Spectrum analysis provides a method to quantitatively estimate the accumulated effect of instantaneous shear stress energy by vortex scales. With a decrease in the length scale, particles endure a larger influence from the fluctuating part of turbulence, and its magnitude could be about 2~4 times the time-averaged part of the flow. At scales of benthic macroinvertebrates, which is a key indicator of aquatic ecosystems, the effect of the instantaneous part of shear stress is not negligible and should be considered in ecological models.
The above fluctuation characteristics of shear stress are related to the dislodgement of sediment and benthic macroinvertebrates on the riverbed. The results and methods of this study will further support future intensive studies on river morphology studies and ecological processes, such as benthic macroinvertebrate drifting and bio-particle transport from a mechanical perspective.

Author Contributions

Data curation, J.W. and Z.L.; methodology, J.W. and Z.L.; supervision, Y.C.; writing—original draft, J.W.; writing—review and editing, Z.L. and M.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Key Research and Development Plan Project of China: Analysis and Evaluation of the Safety Risk of Diversion Tunnel Structures (2019YFB1310504), the National Natural Science Foundation of China (NO.51979142, No. 52009060), the Open Research Fund Program of State Key Laboratory of Hydroscience and Engineering (sklhse-2019-B-06), the China Postdoctoral Science Foundation (2021T140384), and the Open Project of State Key Laboratory of Urban Water Resource and Environment, Harbin Institute of Technology (No. QA202132).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, H.; Chen, Y.; Liu, Z.; Zhu, D. Effects of the “Run-of-River” Hydro Scheme on Macroinvertebrate Communities and Habitat Conditions in a Mountain River of Northeastern China. Water 2016, 8, 31. [Google Scholar] [CrossRef] [Green Version]
  2. Wen, J.Q.; Wang, H.R.; Chen, Y.C.; Liu, Z.W. Longitudinal Distribution of Benthic Macroinvertebrates Affected by a Hydropower Plant Cascade in the Mudan River. Huan Jing Ke Xue = Huanjing Kexue 2020, 41, 3266–3274. [Google Scholar] [PubMed]
  3. Hayes, J.W.; Goodwin, E.O.; Shearer, K.A.; Hicks, D.M. Relationship between background invertebrate drift concentration and flow over natural flow recession and prediction with a drift transport model. Can. J. Fish. Aquat.Sci. 2019, 76, 871–885. [Google Scholar] [CrossRef]
  4. Steuer, J.J.; Newton, T.J.; Zigler, S.J. Use of complex hydraulic variables to predict the distribution and density of unionids in a side channel of the Upper Mississippi River. Hydrobiologia 2008, 610, 67–82. [Google Scholar] [CrossRef]
  5. Gibbins, C.; Vericat, D.; Batalla, R.J. When is stream invertebrate drift catastrophic? The role of hydraulics and sediment transport in initiating drift during flood events. Freshwater Biol. 2007, 52, 2369–2384. [Google Scholar] [CrossRef]
  6. Greimel, F.; Schuelting, L.; Graf, W.; Bondar-Kunze, E.; Auer, S.; Zeiringer, B.; Hauer, C. Hydropeaking Impacts and Mitigation. In Riverine Ecosystem Management: Science for Governing Towards a Sustainable Future; Schmutz, S., Sendzimir, J., Eds.; Springer: Cham, Switzerland, 2018; Volume 8, pp. 91–110. [Google Scholar]
  7. Gibbins, C.; Vericat, D.; Batalla, R.J.; Gomez, C.M. Shaking and moving: Low rates of sediment transport trigger mass drift of stream invertebrates. Can. J. Fish. Aquat.Sci. 2007, 64, 1–5. [Google Scholar] [CrossRef]
  8. Salim, S.; Pattiaratchi, C. Sediment resuspension due to near-bed turbulent coherent structures in the nearshore. Cont. Shelf Res. 2020, 194, 104048. [Google Scholar] [CrossRef]
  9. Blanckaert, K.; Garcia, X.F.; Ricardo, A.M.; Chen, Q.; Pusch, M.T. The role of turbulence in the hydraulic environment of benthic invertebrates. Ecohydrology 2013, 6, 700–712. [Google Scholar] [CrossRef] [Green Version]
  10. Yang, Y.; Gao, S.; Wang, Y.P.; Jia, J.; Xiong, J.; Zhou, L. Revisiting the problem of sediment motion threshold. Cont. Shelf Res. 2019, 187, 103960. [Google Scholar] [CrossRef]
  11. Malakauskas, D.M.; Willson, S.J.; Wilzbach, M.A.; Som, N.A. Flow variation and substrate type affect dislodgement of the freshwater polychaete, Manayunkia speciosa. Freshw. Sci. 2013, 32, 862–873. [Google Scholar] [CrossRef]
  12. Schnauder, I.; Aberle, J.; Rudnick, S.; Garcia, X.F. Incipient motion and drift of benthic invertebrates in boundary shear layers. In International Conference on Fluvial Hydraulics; Bundesanstalt für Wasserbau: Karlsruhe, Germany, 2010; pp. 1453–1461. [Google Scholar]
  13. Sotiropoulos, F.; Yang, X. Immersed boundary methods for simulating fluid-structure interaction. Prog. Aerosp. Sci. 2014, 65, 1–21. [Google Scholar] [CrossRef]
  14. Mignot, E.; Barthelemy, E.; Hurther, D. Double-averaging analysis and local flow characterization of near-bed turbulence in gravel-bed channel flows. J. Fluid Mech. 2009, 618, 279–303. [Google Scholar] [CrossRef]
  15. Wilkes, M.A.; Maddock, I.; Visser, F.; Acreman, M.C. Incorporating Hydrodynamics into Ecohydraulics: The Role of Turbulence in the Swimming Performance and Habitat Selection of Stream-Dwelling Fish. In Ecohydraulics: An Integrated Approach; Maddock, I., Harby, A., Kemp, P., Wood, P., Eds.; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2013; pp. 9–30. [Google Scholar]
  16. Li, J.; Li, S.S. Near-bed velocity and shear stress of open-channel flow over surface roughness. Environ. Fluid Mech. 2020, 20, 293–320. [Google Scholar] [CrossRef]
  17. Kim, S.K.; Choi, S.-U. Prediction of suitable feeding habitat for fishes in a stream using physical habitat simulations. Ecol. Modell. 2018, 385, 65–77. [Google Scholar] [CrossRef]
  18. Orth, D.J.; Maughan, O.E. Evaluation of the Incremental Methodology for Recommending Instream Flows for Fishes. Trans. Am. Fish. Soc. 1982, 111, 413–445. [Google Scholar] [CrossRef]
  19. Jowett, I.G.; Davey, A.J.H. A Comparison of Composite Habitat Suitability Indices and Generalized Additive Models of Invertebrate Abundance and Fish Presence–Habitat Availability. Trans. Am. Fish. Soc. 2007, 136, 428–444. [Google Scholar] [CrossRef]
  20. Holmquist, J.G.; Waddle, T.J. Predicted macroinvertebrate response to water diversion from a montane stream using two-dimensional hydrodynamic models and zero flow approximation. Ecol. Indic. 2013, 28, 115–124. [Google Scholar] [CrossRef] [Green Version]
  21. Morris, M.; Mohammadi, M.H.; Day, S.; Hondzo, M.; Sotiropoulos, F. Prediction of Glossosoma biomass spatial distribution in Valley Creek by field measurements and a three-dimensional turbulent open-channel flow model. Water Resour. Res. 2015, 51, 1457–1471. [Google Scholar] [CrossRef]
  22. Orlandi, P.; Leonardi, S.; Antonia, R.A. Turbulent channel flow with either transverse or longitudinal roughness elements on one wall. J. Fluid Mech. 2006, 561, 279–305. [Google Scholar] [CrossRef]
  23. Volino, R.J.; Schultz, M.P.; Flack, K.A. Turbulence structure in a boundary layer with two-dimensional roughness. J. Fluid Mech. 2009, 635, 75–101. [Google Scholar] [CrossRef] [Green Version]
  24. Li, B.; Jiang, A.; Wu, J. Research on wind resistance characteristics of transverse rod-roughened surface. Acta Aerodyn. Sin. 2016, 34, 517–523. [Google Scholar]
  25. Ahn, J.; Lee, J.H.; Sung, H.J. Statistics of the turbulent boundary layers over 3D cube-roughened walls. Int. J. Heat Fluid Flow 2013, 44, 394–402. [Google Scholar] [CrossRef]
  26. Cheng, H.; Castro, I.P. Near wall flow over urban-like roughness. Bound. Layer Meteorol. 2002, 104, 229–259. [Google Scholar] [CrossRef]
  27. Li, Z.; Liu, Z.; Wang, H.; Chen, Y.; Ling, L.; Wang, Z.; Zhang, D. Investigation of Aerator Flow Pressure Fluctuation Using Detached Eddy Simulation with VOF Method. J. Hydraul. Eng. 2022, 148, 04021052. [Google Scholar] [CrossRef]
  28. Gritskevich, M.S.; Garbaruk, A.V.; Schütze, J.; Menter, F.R.J.F. Turbulence; Combustion, Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model. Flow Turbul. Combust. 2012, 88, 431–449. [Google Scholar] [CrossRef]
  29. Krogstad, P.A.; Andersson, H.I.; Bakken, O.M.; Ashrafian, A. An experimental and numerical study of channel flow with rough walls. J. Fluid Mech. 2005, 530, 327–352. [Google Scholar] [CrossRef]
  30. Wang, M. Progress in data processing, engineering application and mechanism research of water flow fluctuating pressure. Adv. Sci. Technol. Water Resour. 1990, 10, 28–43. [Google Scholar]
  31. Singh, A.; Porte-Agel, F.; Foufoula-Georgiou, E. On the influence of gravel bed dynamics on velocity power spectra. Water Resour. Res. 2010, 46. [Google Scholar] [CrossRef] [Green Version]
  32. Wallace, J.M.; Eckelmann, H.; Brodkey, R.S. The wall region in turbulent shear flow. J. Fluid Mech. 1972, 54, 39–48. [Google Scholar] [CrossRef]
  33. Guan, D.; Agarwal, P.; Chiew, Y.-M. Quadrant Analysis of Turbulence in a Rectangular Cavity with Large Aspect Ratios. J. Hydraul. Eng. 2018, 144, 04018035. [Google Scholar] [CrossRef]
  34. Das, R. Respond of Bedforms to Velocity Power Spectra of Acoustic-Doppler Velocimetry Data in Rough Mobile Beds. Water Resour. 2020, 47, 835–845. [Google Scholar]
  35. Nie, M. Fluctuant characteristics of two-phase flow behind a bottom aerator. Sci. China Ser. E Technol. Sci. 2001, 44, 291–297. [Google Scholar] [CrossRef]
  36. Grass, A.J. Structural Features of Turbulent Flow over Smooth and Rough Boundaries. J. Fluid Mech. 1971, 50, 233–255. [Google Scholar] [CrossRef]
  37. Tang, Z.; Wu, Y.; Jia, Y.; Jiang, N. PIV Measurements of a Turbulent Boundary Layer Perturbed by a Wall-Mounted Transverse Circular Cylinder Element. Flow Turbul. Combust. 2018, 100, 365–389. [Google Scholar] [CrossRef]
  38. Nolan, K.P.; Walsh, E.J.; McEligot, D.M. Quadrant analysis of a transitional boundary layer subject to free-stream turbulence. J. Fluid Mech. 2010, 658, 310–335. [Google Scholar] [CrossRef]
  39. Bouckaert, F.W.; Davis, J. Microflow regimes and the distribution of macroinvertebrates around stream boulders. Freshwater Biol. 1998, 40, 77–86. [Google Scholar] [CrossRef]
  40. Kolmogorov, A.N. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J. Fluid Mech. 1961, 13, 82–85. [Google Scholar] [CrossRef] [Green Version]
  41. Shraiman, B.I.; Siggia, E.D. Scalar turbulence. Nature 2000, 405, 639. [Google Scholar] [CrossRef]
  42. Warhaft, Z. Passive scalars in turbulent flows. Annu. Rev. Fluid Mech. 2000, 32, 203–240. [Google Scholar] [CrossRef]
  43. Qureshi, N.M.; Bourgoin, M.; Baudet, C.; Cartellier, A.; Gagne, Y. Turbulent transport of material particles: An experimental study of finite size effects. Phys. Rev. Lett. 2007, 99, 184502. [Google Scholar] [CrossRef] [Green Version]
  44. Zimmermann, R.; Gasteuil, Y.; Bourgoin, M.; Volk, R.; Pumir, A.; Pinton, J.F. Rotational intermittency and turbulence induced lift experienced by large particles in a turbulent flow. Phys. Rev. Lett. 2011, 106, 154501. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Domain and mesh of the validation simulation.
Figure 1. Domain and mesh of the validation simulation.
Water 14 01752 g001
Figure 2. Normalized time-averaged longitudinal velocity profile of the IDDES calculation results and Krogstad et al.’s measurements [29].
Figure 2. Normalized time-averaged longitudinal velocity profile of the IDDES calculation results and Krogstad et al.’s measurements [29].
Water 14 01752 g002
Figure 3. Normalized turbulent intensities and Reynold stress profiles of IDDES calculation results and Krogstad et al.’s measurements [29].
Figure 3. Normalized turbulent intensities and Reynold stress profiles of IDDES calculation results and Krogstad et al.’s measurements [29].
Water 14 01752 g003
Figure 4. Space correlation along longitudinal direction at y = 1 mm in the prediction model.
Figure 4. Space correlation along longitudinal direction at y = 1 mm in the prediction model.
Water 14 01752 g004
Figure 5. Arrangement of roughness elements on the channel bed: (a) 2D square bars, (b) 3D cube chessboard arrangements, and (c) 3D cube stagger arrangements. Dark squares in (b,c) indicate the locations of roughness elements. Flow direction is along the x direction.
Figure 5. Arrangement of roughness elements on the channel bed: (a) 2D square bars, (b) 3D cube chessboard arrangements, and (c) 3D cube stagger arrangements. Dark squares in (b,c) indicate the locations of roughness elements. Flow direction is along the x direction.
Water 14 01752 g005
Figure 6. Sketch of the locations selected for quadrant analysis (taking 3D-cu-s as an example).
Figure 6. Sketch of the locations selected for quadrant analysis (taking 3D-cu-s as an example).
Water 14 01752 g006
Figure 7. Normalized velocity profiles for predicted results and experiment data [9].
Figure 7. Normalized velocity profiles for predicted results and experiment data [9].
Water 14 01752 g007
Figure 8. Normalized Reynolds stress profiles averaged over roughness elements.
Figure 8. Normalized Reynolds stress profiles averaged over roughness elements.
Water 14 01752 g008
Figure 9. Normalized turbulence kinetic energy profiles averaged over roughness elements.
Figure 9. Normalized turbulence kinetic energy profiles averaged over roughness elements.
Water 14 01752 g009
Figure 10. Normalized (a) urmse, (b) vrmse, and (c) wrmse profiles.
Figure 10. Normalized (a) urmse, (b) vrmse, and (c) wrmse profiles.
Water 14 01752 g010aWater 14 01752 g010b
Figure 11. Quadrant analysis of velocity fluctuations for 2D-sq, (a) type A, (b) type B. The contour represents the joint probability density function of velocity fluctuations u′ and v′. Hyperbolas represent the average shear stress isoline of each type of location, respectively. The dashed line represents the major axis of the elliptic cloud formed by the scattered points.
Figure 11. Quadrant analysis of velocity fluctuations for 2D-sq, (a) type A, (b) type B. The contour represents the joint probability density function of velocity fluctuations u′ and v′. Hyperbolas represent the average shear stress isoline of each type of location, respectively. The dashed line represents the major axis of the elliptic cloud formed by the scattered points.
Water 14 01752 g011
Figure 12. Power spectrum of shear stress at type B points.
Figure 12. Power spectrum of shear stress at type B points.
Water 14 01752 g012
Figure 13. Normalized energy of instantaneous shear stress signal of each case on characteristic length of vortex (from large scale to small scale) at type B locations. Normalization is performed by dividing ( u v ¯ ) 2 , the energy of time-averaged shear stress.
Figure 13. Normalized energy of instantaneous shear stress signal of each case on characteristic length of vortex (from large scale to small scale) at type B locations. Normalization is performed by dividing ( u v ¯ ) 2 , the energy of time-averaged shear stress.
Water 14 01752 g013
Table 1. Configuration of roughness elements.
Table 1. Configuration of roughness elements.
Case No.ArrangementCase NameWidth of Roughness Element b/mmRoughness Height d/mmDistance between Roughness
Elements p/mm
12D square bars2D-sq30860
23D cube chessboard3D-cu-c308 p x = 60 ,   p z = 60
33D cube stagger3D-cu-s308 p x = 120 ,   p z = 120
Table 2. Statistics of different flow resistance parameters.
Table 2. Statistics of different flow resistance parameters.
Case2D-sq3D-cu-c3D-cu-s
u * /m·s−10.02430.02410.0211
Δ U + 11.61911.4249.740
ks/d2.5512.3741.362
Table 3. Proportions of shear stress fluctuations in each quadrant.
Table 3. Proportions of shear stress fluctuations in each quadrant.
Case and LocationsQ1Q2Q3Q4Q2 + Q4
2D-sq A15.38%36.38%23.44%24.79%61.17%
2D-sq B13.68%32.31%17.96%36.05%68.36%
3D-cu-c A17.82%29.14%19.44%33.60%62.74%
3D-cu-c B16.35%34.12%17.82%31.71%65.83%
3D-cu-s A16.49%30.50%17.97%35.04%65.54%
3D-cu-s B11.28%36.71%15.80%36.21%72.93%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wen, J.; Chen, Y.; Liu, Z.; Li, M. Numerical Study on the Shear Stress Characteristics of Open-Channel Flow over Rough Beds. Water 2022, 14, 1752. https://doi.org/10.3390/w14111752

AMA Style

Wen J, Chen Y, Liu Z, Li M. Numerical Study on the Shear Stress Characteristics of Open-Channel Flow over Rough Beds. Water. 2022; 14(11):1752. https://doi.org/10.3390/w14111752

Chicago/Turabian Style

Wen, Jiaqi, Yongcan Chen, Zhaowei Liu, and Manjie Li. 2022. "Numerical Study on the Shear Stress Characteristics of Open-Channel Flow over Rough Beds" Water 14, no. 11: 1752. https://doi.org/10.3390/w14111752

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop