Next Article in Journal
Informal Settlements and Flooding: Identifying Strengths and Weaknesses in Local Governance for Water Management
Previous Article in Journal
Adsorption Behaviors and Removal Efficiencies of Inorganic, Polymeric and Organic Phosphates from Aqueous Solution on Biochar Derived from Sewage Sludge of Chemically Enhanced Primary Treatment Process
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Difference of River Resistance Computation between the k ε Model and the Mixing Length Model

1
College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China
2
State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing 210098, China
3
National Engineering Research Center of Water Resources Efficient Utilization and Engineering Safety, Hohai University, Nanjing 210098, China
*
Author to whom correspondence should be addressed.
Water 2018, 10(7), 870; https://doi.org/10.3390/w10070870
Submission received: 22 May 2018 / Revised: 21 June 2018 / Accepted: 27 June 2018 / Published: 29 June 2018
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
River resistance characteristics, which can be reflected by the resistance factor, have an impact on flow and sediment transport. In the classical theory, Prandtl proposed the mixing length model for the simulation of the turbulence, and von Kármán established the logarithmic formula of the flow velocity distribution. Based on that, the expression of the resistance factor can be derived. With the development of the numerical technology, the k ε model has been widely applied in the channels computation. However, for the different closure ways of the Reynolds stress in turbulence equations, the outcomes of the k ε model and the Prandtl mixing length model are not exactly identical. In this paper, both qualitative and quantitative studies are carried out on the difference between these two models, with respect to the resistance factor. This difference is evaluated by the ratio of the resistance factor computed with the two models. The result shows that with the increment of the relative flow depth, the ratio first increases and then decreases. Moreover, it is also affected by the bed slope. Therefore, the difference should be taken into account when a comparison is made between the simulation results of the k ε model and the classical theory of river mechanics.

1. Introduction

In the (vertical) two-dimensional (2D) steady uniform open channel flow, as shown in Figure 1, the vertical distribution of velocity u can be expressed by the Prandtl-von Kármán logarithmic law, as follows:
u u * = 1 κ ln ( z k s ) + B s
where z is the vertical coordinate; u is the time-averaged velocity in the x direction; u * = τ b / ρ is the shear velocity, τ b is the bed sheer stress; κ is the von Kármán constant, κ = 0.4 ; k s is the equivalent bed-roughness and is discussed in Section 3; B s = ϕ B s ( R e * ) is the roughness function, R e * = u * k s / ν is the roughness Reynolds number; ν is the kinematic viscosity coefficient. B s can be determined by experimental curve [1,2] or calculated approximately by the following formula [1]:
B s = ( 2.5 ln R e * + 5.5 ) exp [ 0.0705 ( ln R e * ) 2.55 ] + 8.5 { 1 exp [ 0.0594 ( ln R e * ) 2.55 ] }
Remarkably, when R e * 70 , the flow is the rough turbulence, and thus B s 8.5 ; the value is independent of R e * . As a semi-empirical solution, Equation (1) is widely applied in the calculation of open channel flows.
The dimensionless Chézy resistance factor c is defined as c = u ¯ / u * [1], where u ¯ is the depth-averaged velocity. Considering u ¯ being equal to u at the relative height z / h = e 1 = 0.368 [1], where h is the depth, the Chézy resistance factor c is given by the following equation:
c = u ¯ u * = 1 κ ln ( 0.368 h k s ) + B s
which is also the classical formula for c .
For the two-dimensional open channel uniform flow, there exists the following equation:
u * = g h S
where S is bed slope along the channel. By combining Equations (3) with (4), u ¯ can be expressed as follows:
u ¯ = c u * = c g h S
Equation (5) is in the form of the Chézy formula, in which the dimensionless resistance factor c has a relationship with the usual Chézy coefficient given by the following:
C = c g
The deduction of Equation (1) relies on the closure of the Reynolds stress in turbulence equations using the Prandtl mixing length model [2] (the details on how Equation (1) is derived from the Prandtl mixing length model can be seen in [3,4]). Hence, Equation (3) is also based on the Prandtl mixing length model. However, with the development of numerical technology, the k ε turbulence model is more widely applied than the Prandtl mixing length model. In recent years, by employing the k ε turbulence model, Wu et al. [5] simulated the flow and sediment transport in a 180° channel bend with movable bed in the laboratory; Zhang et al. [6] developed the two-dimensional k ε turbulence model in curvilinear coordinates to simulate the behavior of turbulent flow in an open channel partially covered with vegetation; Thi et al. [7] investigated the turbulent flow around the local scour occurring around hydraulic structures with the k ε turbulence model. On the basis of the three-dimensional geometry model with a geological model, Wang et al. [8] presented a three-dimensional k ε turbulence model considering Bingham liquid character. Besides these, the k ε turbulence model is also applied in solving problems of other fields, such as air cooling technologies [9], automotive applications of on-board hydrogen storage [10], internal combustion engines scale-resolving simulation [11], pollutant removal [12], etc.
In general, the treatment of each layer in the k ε turbulence model is shown in Figure 2, in which u 1 is the velocity near the bed; z 1 is a selected near-bed height. This treatment assumes that u conforms to the Prandtl-von Kármán logarithmic law only within the range of z 1 . Above this range, however, the velocity is computed by solving the turbulence equation closed by the k ε model. On account of that, the resistance factor c corresponding to the k ε model is not exactly identical to the one calculated by Equation (3). Note that Equation (3) is indeed the outcomes corresponding to the Prandtl mixing length model applied to the whole range of the depth. Thus, the fundamental reason causing the divergence is that there are two turbulence models using different ways in the closure of the Reynolds stress in turbulence equations. In this paper, c represents the outcome of the k ε model and c represents that of the Prandtl mixing length model to distinguish them. In addition, the ratio c / c can be regarded as the quantitative parameter measuring the difference between the two models.
In this study, the hydraulic parameters of natural rivers are collected, from which c and c can be computed by using the above two models. On this basis, the variation pattern of c / c is studied.

2. Computational Methods of Resistance Factors by Using Different Models

Either way, the basic steps of computing the resistance factor are as follows:
  • Solve the turbulence equations to obtain the vertical distribution of u ;
  • Average u over the depth to obtain u ¯ ; and
  • Calculate the resistance factor c or c by the ratio of u ¯ to u * .
Therefore, the core issue is how to solve the turbulence equations. The solving methods using the Prandtl mixing length model and the k ε model are given separately in the below.

2.1. Prandtl Mixing Length Model

As stated in the introduction, solving the turbulence equations closed by the Prandtl mixing length model actually gives the Prandtl-von Kármán logarithmic distribution of u (i.e., Equation (1) in Section 1), which is in the continuous form. After integrating this equation, it can be identified that u ¯ is at the relative height z / h = e 1 = 0.368 . For that, Equation (3) is used to compute c in this paper.

2.2. k ε Model

Generally, a continuous form of velocity distribution cannot be obtained by closing and solving the turbulence equation with the k ε model; thus, the discrete method is used. On the conditions of the two-dimensional uniform flow in open channels, the RANS (Reynolds-averaged Navier–Stokes) equation in the x direction [2] can be simplified into the following:
g S = z [ ν u z u w ]
where, u and w are the fluctuant velocity components in x and z directions, respectively;   denotes the time-averaged value of the content within the angle brackets, and thus u w actually denotes the Reynolds stress per unit mass.
According to the Boussinesq hypothesis, for the turbulent flow, ν t is defined as the eddy viscosity. Reynolds stress has a relationship with the time-averaged velocity gradient given by the following:
u w ¯ = ν t d u d z
Then, Equation (7) can be transferred to the following:
g S = z [ ( ν + ν t ) u z ]
Equation (9) is just the basic equation for the computation of c by using the k ε model. Then, the finite volume method is adopted to discretize this equation.
The integral of Equation (10) on the grid of the i-th layer, as shown in Figure 3, gives the following:
g S Δ z = [ ( ν + ν t ) u z ] i + 1 / 2 [ ( ν + ν t ) u z ] i 1 / 2
In which the differential terms can be further discretized with the central differencing scheme, as follows:
( u z ) i + 1 / 2 = Δ u i + 1 / 2 Δ z i + 1 / 2 = u i + 1 u i Δ z
( u z ) i 1 / 2 = Δ u i 1 / 2 Δ z i 1 / 2 = u i u i 1 Δ z
The substitution of Equations (11) and (12) into Equation (10) gives the following:
ν t i 1 / 2 + ν Δ z 2 u i 1 + ν t i 1 / 2 + 2 ν + ν t i + 1 / 2 Δ z 2 u i ν t i + 1 / 2 + ν Δ z 2 u i + 1 = g S
which is the general discrete form of Equation (9).
As shown in Figure 4, for the situation that i = 1 , Equation (10) is transferred to the following:
g S Δ z = ( ν + ν t 3 / 2 ) ( u z ) 3 / 2 ( ν + ν t 1 / 2 ) ( u z ) 1 / 2
The derivative term in the first term on the right-hand side can be discretized with the central differencing scheme as follows:
( u z ) 3 / 2 = Δ u 3 / 2 Δ z 3 / 2 = u 2 u 1 Δ z
According to the bed boundary condition, the second term on the right-hand side can be calculated as follows:
( ν + ν t 1 / 2 ) ( u z ) 1 / 2 = τ b / ρ = u * 2
Substituting Equations (15) and (16) into Equation (14), one obtains the following:
ν + ν t 3 / 2 Δ z 2 u 1 ν + ν t 3 / 2 Δ z 2 u 2 = g S u * 2 Δ z
As shown in Figure 5, for the situation that i = N , Equation (10) is transferred to the following:
g S Δ z = ( ν + ν t N + 1 / 2 ) ( u z ) N + 1 / 2 ( ν + ν t N - 1 / 2 ) ( u z ) N 1 / 2
The derivative term in the second term on the right-hand side can be discretized with the center difference scheme as follows:
( u z ) N 1 / 2 = Δ u N 1 / 2 Δ z N 1 / 2 = u N u N 1 Δ z
According to the free surface boundary condition, the first term on the right-hand side can be calculated as follows:
( ν + ν t N + 1 / 2 ) ( u z ) N + 1 / 2 = 0
Substituting Equations (19) and (20) into Equation (18), one obtains the following:
ν + ν t N 1 / 2 Δ z 2 u N 1 + ν + ν t N 1 / 2 Δ z 2 u N = g S
For the value of ν t , at the interface of the grid, the linear interpolation is resorted to the following:
ν t i 1 / 2 = 1 2 ( ν t i 1 + ν t i )
ν t i + 1 / 2 = 1 2 ( ν t i + ν t i + 1 )
whereas the value of ν t at the center of the grid is calculated by the k ε model, as follows:
ν t = c ν k 2 ε
z [ ν + ν t σ k k z ] = P z ε
z [ ν + ν t σ ε ε z ] = ( c ε 1 P z c ε 2 ε ) ε k
where k is the turbulent kinetic energy; ε is the dissipation rate; P z = ( ν + ν t ) (   u / z ) 2 , is the production term; the standard values of the model coefficients are employed [5,13]: c ν = 0.09, σ k = 1.0, σ ε = 1.3, c ε 1 = 1.44, c ε 2 = 1.92. Equations (25) and (26) are both differential equations, and they are discretized with the same scheme as Equation (9). In addition, for Equations (25) and (26), the bed boundary conditions [5] are as follows:
k 1 = u * 2 c ν ,   ε 1 = u * 3 κ z 1
and the free surface boundary conditions [5] are as follows:
( k z ) N + 1 / 2 = 0 ,   ε N = k N 3 / 2 0.43 h
Combing Equation (13) with Equations (17) and (21), the following tridiagonal linear system is finally derived, as follows:
[ β 1 γ 1 α 2   β 2 γ 2 α i β i γ i α N 1 β N 1 γ N 1 α N β N ] [ u 1 u 2 u i u N 1 u N ] = [ g S u * 2 / Δ z g S g S g S g S ]
where,
β 1 = ν + ν t 3 / 2 Δ z 2 ,   γ 1 = ν + ν t 3 / 2 Δ z 2 α 2 = ν + ν t 3 / 2 Δ z 2 ,   β 2 = ν t 3 / 2 + 2 ν + ν t 5 / 2 Δ z 2 ,   γ 2 = ν + ν t 5 / 2 Δ z 2 α i = ν + ν t i 1 / 2 Δ z 2 ,   β i = ν t i 1 / 2 + 2 ν + ν t i + 1 / 2 Δ z 2 ,   γ i = ν + ν t i + 1 / 2 Δ z 2 α N 1 = ν + ν t N 3 / 2 Δ z 2 ,   β N 1 = ν t N 3 / 2 + 2 ν + ν t N 1 / 2 Δ z 2 ,   γ N 1 = ν + ν t N 1 / 2 Δ z 2 α N = ν + ν t N 1 / 2 Δ z 2 ,   β N = ν + ν t N 1 / 2 Δ z 2
However, the determinant of the coefficient matrix of this system turns out to be zero and the rank is N 1 ( < N , non-full rank), which means that there exists infinitely many solutions to this system. Furthermore, this system is derived as follows.
By means of accumulative summation from the top to the bottom, Equations (29) can be transferred to the following:
[ γ 1 γ 2 γ i γ N 1 0 ] [ u 2 u 1 u 3 u 2 u i + 1 u i u N u N 1 0 ] = [ g S u * 2 / Δ z 2 g S u * 2 / Δ z i g S u * 2 / Δ z ( N 1 ) g S u * 2 / Δ z N g S u * 2 / Δ z ]
The last equation of Equations (30) is as follows:
N g S u * 2 / Δ z = 0
Note that N Δ z = h . Then, after transposition, Equation (31) can be transferred to
u * = g ( N Δ z ) S = g h S
which is the same equation as Equation (4).
From Equations (30), it follows that only the difference between the velocities of any two adjacent points can be solved. Thus, in order to get the absolute value of the velocity, another supplementary equation is needed. As stated in the introduction, it is generally assumed that the near-bed velocity within the range of z 1 ( z 1 0.2 h ) [4,14] conforms to the Prandtl-von Kármán logarithmic law, which means the following:
u 1 u * = 1 κ ln ( z 1 k s ) + B s
Then, the velocity u i (i = 1, 2, …, N) of any point can be obtained by successively solving Equation (33), and then Equations (30) from the top to the bottom. The depth-averaged velocity u ¯ can be evaluated as follows:
u ¯ = 1 N i = 1 N u i
Thus, the resistance factor c can be evaluated as c = u ¯ / u * .
The above calculation still needs iteration for γ in the coefficient matrix of Equations (30), involving ν t = c ν k 2 / ε , which is determined by the differential Equations (25) and (26).
The computation methods of c using the Prandtl mixing length model and the k ε model are given, as above. Figure 6 shows the typical computations of velocity distribution by using these two models. The parameters of the Rio Grande River [15] are used in both computations: h = 0.332 m, k s = 0.028 m, and S = 0.830 × 10−3. The result corresponding to Prandtl mixing length model, which is calculated by Equation (1), is shown in the continuous form; in contrast, the result corresponding to the k ε model, which is calculated by dividing the whole depth into ten layers, is shown in the discrete form. From Figure 6, an approximate coincidence can be detected between the two velocity distribution pattern computed by different models within the range of z 0.2 h ; however, the difference displays, to some extent, within the range of z > 0.2 h . Indeed, the Prandtl mixing length model and the k ε model treat the shear force between the flow layers in different ways. This is the main reason the resistance factors computed with these two models are not identical.

3. Determination of Equivalent Bed-Roughness

In natural rivers, bed forms are induced by sediment transport. This means that the bed resistance to flow includes not only the grain resistance, but also the bed-form resistance—For some special cases, the resistance of other factors (e.g., the emergent vegetation) should also be taken into account [16], but they are not considered in this study. In this study, the van Rijn method [17] is employed to determine k s , which takes the influence of these two kinds of resistance into consideration.
The expression of k s proposed by van Rijn [17] is as follows:
k s = 3 D + 1.1 Δ [ 1 exp ( 25 δ ) ]   ( 0.01 δ 0.2 )
where D is the characteristic bed material diameter; Δ and δ are the bed-form height and the bed-form steepness, respectively, which are determined as follows:
Δ h = 0.110 ( h D ) 0.3 [ 1 exp ( 0.5 T ) ] ( 25 T )
δ = 0.015 ( h D ) 0.3 [ 1 exp ( 0.5 T ) ] ( 25 T )
where T is a transport stage parameter,
T = ( τ b ) f ( τ b ) c r ( τ b ) c r
where ( τ b ) f is the bed shear stress related to grains; ( τ b ) c r is the critical bed-shear stress for initiation of motion of grains, which can be determined by experimental curve [1,18] or calculated approximately by the following formula:
( τ b ) c r ( ρ s ρ ) g D = 0.13 Ξ 0.392 exp ( 0.015 Ξ 2 ) + 0.045 [ 1 exp ( 0.068 Ξ ) ]
where Ξ is the dimensionless material number, which is defined as the following:
Ξ = ( ρ s ρ ρ g D 3 ν 2 ) 1 / 3
where ρ s is the density of the bed material, ρ s = 2650 kg/m3; ρ is the density of water, ρ = 1000 kg/m3.
In the work of van Rijn [17], D was identified with D 90 ; in this study, however, as a result of the lack of information of D 90 reported in the references, D is thus identified from all the available data within a range of D 50 D 90 . Practices have indicated that although the exact value of c / c varies with the D selected, the pattern of c / c versus D does not change in a noticeable manner. This is further discussed and evidenced in Section 5.1.

4. Data, Results, and Analyses

4.1. Data Sources and Range

The data of natural rivers were collected from previous publications [15,19,20,21,22,23], including the discharge Q , the water depth h , the bed slope S , and the characteristic bed material diameter D . There are 213 groups data in total. The detailed information is shown in Table 1.

4.2. Analysis of Influential Factors

According to Section 2 and Section 3, c / c is influenced by the following seven dimensional parameters:
ρ , ρ s , g , ν , h , D , S
ρ , g and D are selected as the basic dimensionally independent variables that represent the dynamic, kinematic, and geometric characteristics, respectively. Then, the Buckingham π theorem is employed for the dimensional analysis. One determines for c / c that is influenced by the following four dimensionless parameters:
ρ s ρ , D g D ν , h D , S
within these parameters, ρ s / ρ is considered as constant (=2.65) in this paper. Thus, the influence of ρ s / ρ is not considered in the following analyses.
The product of D g D / ν , h / D and S gives the following:
D g D ν h D S = D g h S ν = u * D ν = R e D
where R e D is the grain size Reynolds number. According to Equation (35), the roughness Reynolds number R e * can be expressed as follows:
R e * = 3 R e D + 1.1 R e Δ ( 1 e 25 δ )
Here, R e Δ = u * Δ / ν , which is defined as the bed-form Reynolds number. According to Equation (44), R e D is included in R e * . Also, note that R e * only appears in the roughness function B s in Equations (3) and (33). For the data of natural rivers collected in this study, it can be proved that R e * 70 , meaning that the flow is an invariably rough turbulent, for which B s 8.5 ; as a consequence, in this study, B s is independent of R e * and thus R e D . Then, the influence of D g D / ν is neglected.
In the final analysis, c / c is a function solely depending on h / D and S , for which one can write the following:
c c = f ( h D , S )
Note that h / D is referred to as the relative water depth [1]. The function form of f will be determined in the following subsection.

4.3. Computational Results and Analyses

First, by clustering the data of S and h / D corresponding to all the collected natural rivers, the relation of S versus h / D is plotted in Figure 7 with log−log axes. From Figure 7, it can be seen that all the ( h / D , S ) points are approximately distributed around the following five horizontal lines: S = 10.050 × 10−3, 1.540 × 10−3, 0.806 × 10−3, 0.190 × 10−3, and 0.072 × 10−3. Note that despite a considerable scatter, the ( h / D , S ) points around S = 0.190 × 10−3 still follow this pattern. Thus, a classification will be made in the following contents for all the collected natural rivers on the basis of these five values of S .
Then, with the values of c / c corresponding to all the collected natural rivers computed, the variation pattern of c / c versus h / D is plotted in Figure 8 with lin−log axes. As shown in Figure 8, all the ( h / D , S ) points can be classified into five groups, according to different values of S , as specified above; each group is represented by a continuous curve. Although the ( h / D , S ) points around curve 4 ( S = 0.190 × 10−3) are scattered, the general pattern can still be followed. With the increment of h / D , only the ( h / D , c / c ) points around curve 3 show a trend of ascending followed by descending, which is considered to be ‘complete’. In contrast, the ( h / D , c / c ) points around curve 1 ( S = 10.050 × 10−3) only compose an ascending part, and so do those around curve 2 ( S = 1.540 × 10−3); the ( h / D , c / c ) points around curve 4 ( S = 0.190 × 10−3) only compose a descending part, and so do those around curve 5 ( S = 0.072 × 10−3). Despite insufficient evidence, it can be inferred that the ( h / D , c / c ) points around each of these five curves should have the same trend. It is due to the limited range of data that the ( h / D , c / c ) points around each of curves 1, 2, 4, and 5 cannot show a trend as ‘complete’ as those around curve 3.
The ‘complete’ curve 3 can be approximately expressed by fitting the ( h / D , c / c ) points around it as follows:
c c = [ 0.27 ln ( h D ) + 0.09 ] [ 1.75 ln ( h D ) 11.20 + 1 ]
It is assumed that curves 1, 2, 4, and 5 can be obtained by translating curve 3 on the ( h / D , c / c ) plane. Thus, two translation parameters Δ x and Δ y (see Figure 8) are introduced into Equation (46), which yields the following:
c c = { 0.27 [ ln ( h D ) + Δ x ] + 0.09 } { 1.75 [ ln ( h D ) + Δ x ] 11.20 + 1 } + Δ y
Apparently, both of Δ x and Δ y should be functions of S.
The values of Δ x and Δ y for curves 1, 2, 4, and 5 can be obtained by fitting Equation (47) to the corresponding ( h / D , c / c ) points. Taking also into account Δ x = 0 and Δ y = 0 for curve 3, the relations of Δ x versus S , and Δ y versus S are plotted in Figure 9a,b, respectively, with lin−log axes. It is indicated that both Δ x and Δ y are linear with the logarithm of S ( ln S ). With the premise that the fitting lines pass through the points corresponding to curve 3 in Figure 8, the fitting equations are constructed as follows:
Δ x = 1.11 ln S + 7.90
Δ y = 0.02 ln S + 0.15
The substitution of Equations (48) and (49) into Equation (47) yields the following:
c c = [ 0.27 ln ( h D ) + 0.30 ln S + 2.22 ] [ 1.75 ln ( h D ) + 1.11 ln S 3.30 + 1 ] + 0.02 ln S + 0.15 = ln [ 9.21 ( h D ) 0.27 S 0.30 ] { 1.75 ln [ 0.04 ( h D ) S 1.11 ] + 1 } + ln ( 1.16 S 0.02 ) = f ( h D , S )
This is the functional form sought for f . Note that there exists a vertical asymptote on the right side of each c / c h / D curve, as follows:
ln ( h D ) + 1.11 ln S 3.30 = 0
Thus, Equation (51) is only applicable to the following:
h D < 27.11 S 1.11
According to Equation (50), together with all the curves in Figure 8, the following aspects can be summarized:
(a)
For a specific S , c / c first increases and then decreases with the increment of h / D ; the variation pattern of c / c versus h / D is shown as a convex upward curve;
(b)
The peak value of the convex upward c / c h / D curve increases with the increment of S ; Figure 9b, along with Equation (49), indicates that the peak value is approximately linear with ln S ; and
(c)
The abscissa of the peak value decreases with the increment of S ; Figure 9a, along with Equation (48), indicates that the abscissa of the peak value is approximately linear with ln S as well.
Moreover, from the ( h / D , c / c ) points in Figure 8, it follows that the values of c / c are generally greater than unity. This means that for the collected natural rivers, the resistance factors computed with the k ε model are generally greater than those computed with the Prandtl mixing length model. However, it cannot be ignored that the circumstances where the values of c / c are smaller than unity can also be detected, as long as the corresponding values of h / D are sufficiently large (see Figure 8, the points on the lower right sides of curves 3 and 4).
In practice, if the data collected in a natural river include the water depth h , the bed slope S , and the characteristic bed material diameter D , all of which satisfy Equation (52), then the flow discharge Q corresponding to the k ε turbulence model can be simply calculated with the following steps, instead of by running a complex numerical model:
  • Calculate the resistance factor c corresponding to the Prandtl mixing length model by employing Equation (3) (see Section 2.1);
  • Covert the resistance factor from the one c corresponding to the Prandtl mixing length model into the one c corresponding to the k ε model by employing Equation (50);
  • Substitute c instead of c into Equation (5) to calculate the depth-averaged velocity u ¯ ; and
  • Calculate the flow discharge Q as Q = u ¯ B h .

4.4. Discharge Verification

The computed flow discharge Q and Q corresponding to the k ε model and the Prandtl mixing length model, respectively, are verified by comparing them with the field-measured flow discharge Q m . The results are shown in Figure 10. It can be seen that both the ( Q m ,   Q ) and ( Q m ,   Q ) sets of points are evenly distributed around the line y = x .
Furthermore, the frequency distributions of the computational errors of Q and Q m relative to Q m are shown in Figure 11. Both approximately conform to the normal distribution as follows:
Q   relative   error N ( μ = 21 . 65 % ,   σ 2 = 28 . 79 2 )
Q   relative   error N ( μ = 7.78 % ,   σ 2 = 31.11 2 )
where N is the normal distribution; μ is the mathematical expectation of normal distribution; and σ 2 is the variance.
From the values of μ in Equations (53) and (54), the mean relative errors of Q and Q are 21.65% and 7.78%, respectively. This means that the ( Q m ,   Q ) points in Figure 10b are closer to the line y = x than the ( Q m ,   Q ) points in Figure 10a. From the values σ2 of in Equations (53) and (54), the variance of the relative errors of Q and Q are 28.792 and 31.112, respectively. This means that the ( Q m ,   Q ) points in Figure 10b are more scattered than the ( Q m ,   Q ) points in Figure 10a. Therefore, each of two turbulence models has its own advantages and disadvantages, if considered from different aspects.
Considering the function f proposed in Section 4.3, which quantitatively describes the difference between two turbulence models in the computation of river resistance, Q and Q corresponding to the k ε model and the Prandtl mixing length model, respectively, seem interconvertible, as follows:
Q = Q f ( h D , S )   or   Q = Q / f ( h D , S )
As above, it is further indicated that the simulation results corresponding to the k ε model and the Prandtl mixing length model are not identical. When comparing them, the difference between the river resistance characteristics reflected by these two models should be taken into account.

5. Discussion

5.1. Influence of the Selection of the Characteristic Bed Material Diameter

As stated in Section 3, the calculation of the equivalent bed-roughness k s involves a certain characteristic bed material diameter D. In this study, the available data of D are within a range of D50D90. Indeed, with regards to Figure 8, even though the a different type of D is selected for each point, the same trend of the computed c / c h / D curves can still be obtained. The points are still distributed around the same one out of the five curves, although positioned at different positions.
Only a few references report various types of D for a certain river; (e.g., from Rio Grande River [15]) four groups of data corresponding to different river reaches of both D50 and D84 can be extracted. In Figure 8, the ( h / D , c / c ) points corresponding to these river reaches are computed with D50; these points, together with curve 3, which they are fitted to, are replotted in Figure 12. To show the influence of the selection of the type of D, the values of c / c corresponding to the same river reaches are further recomputed with D84, and the corresponding ( h / D , c / c ) points are plotted in Figure 12, too. It can be seen from Figure 12 that the ( h / D , c / c ) points computed with D84 are also distributed around curve 3. The pattern of c / c versus h / D does not change in a noticeable manner if a different type of D is selected.

5.2. Influence of the Coefficients in the k ε Model

As stated in Section 2.2, the k ε turbulence model involves five coefficients and their standard values are used, namely, c v = 0.09, σ k = 1.0, σ ε = 1.3, c ε 1 = 1.44, and c ε 2 = 1.92. However, according to Rodi [13], in a case where the production of the turbulence kinetic energy is much larger than its dissipation, the value of c v can be as small as about 0.05. Figure 13 displays a partial result of c / c versus h / D , with different values of c v : c v = 0.09 and c v = 0.05. From Figure 13, it follows that the computed values of c / c with c v = 0.05 are somewhat larger than those with c v = 0.09. Consequently, the result presented in this paper still depends on the values of c v in the k ε model. Nevertheless, it should also be noted that with either value of c v , c / c appears invariably as a function of h / D ; that is, it can never be constant.
For lack of knowledge on how the other coefficients change according to different cases, their influences are not included in this paper. How these coefficients, including c v , affect the resistance factor ratio c / c needs further studies.

6. Conclusions

With the data of natural rivers, the computational results corresponding to the k ε model and the Prandtl mixing length model are compared and analyzed. The following has been suggested:
(a)
The resistance factors corresponding to these two models are not identical. A new expression Equation (50), has been derived to quantify it. The difference of these two models, evaluated by the ratio c / c , first increases and then decreases with the increment of the relative water depth h / D . This variation pattern can be represented by a convex upward curve, of which the peak value increases and the abscissa of the peak value decreases with the increment of the bed slope S .
(b)
Both turbulence models can generate reasonable results with acceptable errors with respect to the flow discharges in natural rivers. Each model has its own advantages and disadvantages. Within the range of the collected natural rivers, the mean relative error of the computed flow discharges corresponding to the Prandtl mixing length model is less than that corresponding to the k ε model; in contrast, the error distribution corresponding to the former is more scattered than that corresponding to the latter. Furthermore, the flow discharges computed with these two models can be converted into each other by taking into account Equation (55).
(c)
The velocity distribution patterns corresponding to these two models within the range of z 0.2 h are approximately consistent, whereas the difference displays to some extent within the range of z > 0.2 h .
(d)
For river computation issues, the simulation results corresponding to the k ε model and the Prandtl mixing length model are not exactly identical, although both can well simulate the flow. Accounting for that, the difference between the river resistance characteristics reflected by these two models should be considered.

Author Contributions

W.D. carried out the research, including framing the scope of the work and writing the manuscript. M.D. computed the data with models and analyzed the results. H.Z. provided technical supports of models and review the manuscript.

Funding

This research was funded by the National Key R&D Program of China (2016YFC0402501), the NSFC (51479071), the 111 project (B12032, B17015), the Priority Academic Program Development of Jiangsu Higher Education Institutions (3014-SYS1401), and the College Students’ Innovation and Entrepreneurship Training Program Project of Hohai University (201710294016-Zingming He, 2017102941059-Xiaoting Wang).

Acknowledgments

The authors would like to acknowledge Ishwar Joshi and Chaoyang Dai for their help of reviewing the paper.

Conflicts of Interest

The authors declare no conflict of interests.

References

  1. Yalin, M.S.; da Silva, A.M.F. Fluvial Process; International Association of Hydraulic Engineering and Research (IAHR): Delft, The Netherlands, 2001. [Google Scholar]
  2. Zhang, Z.X.; Dong, Z.N. Viscous Fluid Mechanics, 2nd ed.; Tsinghua University Press: Beijing, China, 2011. (In Chinese) [Google Scholar]
  3. Yalin, M.S. Mechanics of Sediment Transport; Pergamon: Oxford, UK, 1972. [Google Scholar]
  4. Nezu, I.; Nakagawa, H. Turbulence in Open-Channel Flows; A.A. Balkema: Rotterdam, The Netherlands, 1993. [Google Scholar]
  5. Wu, W.; Rodi, W.; Wenka, T. 3D numerical modeling of flow and sediment transport in open channels. J. Hydraul. Eng. 2000, 126, 4–15. [Google Scholar] [CrossRef]
  6. Zhang, M.L.; Shen, Y.M.; Zhu, L.Y. Depth-averaged two-dimensional numerical simulation for curved open channels with vegetation. J. Hydraul. Eng. 2008, 39, 794–800. (In Chinese) [Google Scholar]
  7. Thi, H.T.N.; Jungkyu, A.; Sung, W.P. Numerical and physical investigation of the performance of turbulence modeling schemes around a scour hole downstream of a fixed bed protection. Water 2018, 10, 103. [Google Scholar] [CrossRef]
  8. Wang, X.L.; Wang, Q.S.; Zhou, Z.Y.; Ao, X.F. Three-dimensional turbulent numerical simulation of Bingham fluid in the goaf grouting of the south-to-north water transfer project. J. Hydraul. Eng. 2013, 44, 1295–1302. (In Chinese) [Google Scholar]
  9. Wibron, E.; Ljung, A.-L.; Lundström, T.S. Computational fluid dynamics modeling and validating experiments of airflow in a data center. Energies 2018, 11, 644. [Google Scholar] [CrossRef]
  10. Wookyung, K.; Volodymyr, S.; Dmitriy, M.; Vladimir, M. Simulations of blastwave and fireball occurring due to rupture of high-pressure hydrogen tank. Safety 2017, 3, 16. [Google Scholar] [CrossRef]
  11. Vesselin, K.K.; Luca, S.; Giacomo, F. A modified version of the RNG k–ε turbulence model for the scale-resolving simulation of internal combustion engines. Energies 2018, 10, 2116. [Google Scholar] [CrossRef]
  12. Liu, F.; Qian, H.; Zheng, X.; Zhang, L.; Liang, W. Numerical study on the urban ventilation in regulating microclimate and pollutant dispersion in urban street canyon: A case study of Nanjing new region, China. Atmosphere 2017, 8, 164. [Google Scholar] [CrossRef]
  13. Rodi, W. Turbulence Model and Their Application in Hydraulics—A State-of-the-Art Review, 3rd ed.; A.A. Balkema: Rotterdam, The Netherlands, 1993. [Google Scholar]
  14. Nezu, I.; Rodi, W. Open-channel flow measurements with a laser Doppler anemometer. J. Hydraul. Eng. 1986, 112, 335–355. [Google Scholar] [CrossRef]
  15. Toffaleti, F.B. A Procedure for Computation of the Total River Sand Discharge and Detailed Distribution, Bed to Surface; US Army Corps of Engineers: Washington, DC, USA, 1968.
  16. Meftah, M.B.; Serio, F.D.; Malcangio, D.; Mossa, M. Resistance and boundary shear in a partly obstructed channel flow. In River Flow 2016; Constantinescu, G., Garcia, M., Hanes, D., Eds.; Taylor & Francis Group: London, UK, 2016; pp. 795–801. [Google Scholar]
  17. Rijn, L.C.V. Sediment transport, part III: Bed forms and alluvial roughness. J. Hydraul. Eng. 1984, 110, 1733–1754. [Google Scholar] [CrossRef]
  18. Shao, X.J. Introduction to River Mechanics; Tsinghua University Press: Beijing, China, 2013. (In Chinese) [Google Scholar]
  19. Da Cunha, L.V. River Mondego, Portugal; 1969. In Brownlie, W.R. Compilation of Alluvial Channel Data: Laboratory and Field; California Institute of Technology: Pasadena, CA, USA, 1981. [Google Scholar]
  20. Einstein, H.A. Bed-Load Transportation in Mountain Creek; US Department of Agriculture, Soil Conservation Service: Washington, DC, USA, 1944. [Google Scholar]
  21. Leopold, L.B. Sediment Transport Data for Various U.S. Rivers; 1969. In Brownlie, W.R. Compilation of Alluvial Channel Data: Laboratory and Field; California Institute of Technology: Pasadena, CA, USA, 1981. [Google Scholar]
  22. Milhous, R.T. Sediment Transport in a Gravel-Bottomed Stream. Ph.D. Thesis, Oregon State University, Corvallis, OR, USA, 1973. [Google Scholar]
  23. Simons, D.B. Theory of Design of Stable Channels in Alluvial Materials. Ph.D. Thesis, Colorado State University, Fort Collins, CO, USA, 1957. [Google Scholar]
Figure 1. Sketch of two-dimensional (2D) open channel uniform flow.
Figure 1. Sketch of two-dimensional (2D) open channel uniform flow.
Water 10 00870 g001
Figure 2. Treatment of velocity in different layers by using the k ε turbulence model.
Figure 2. Treatment of velocity in different layers by using the k ε turbulence model.
Water 10 00870 g002
Figure 3. Grid discretization (i = 2, 3, …, N − 1).
Figure 3. Grid discretization (i = 2, 3, …, N − 1).
Water 10 00870 g003
Figure 4. Grid discretization for the near-bed layer (i = 1).
Figure 4. Grid discretization for the near-bed layer (i = 1).
Water 10 00870 g004
Figure 5. Grid discretization for the near-surface layer (i = N).
Figure 5. Grid discretization for the near-surface layer (i = N).
Water 10 00870 g005
Figure 6. Typical velocity distributions computed with two turbulence models (Rio Grande River [15]: h = 0.332 m, k s = 0.028 m, and S = 0.830 × 10−3).
Figure 6. Typical velocity distributions computed with two turbulence models (Rio Grande River [15]: h = 0.332 m, k s = 0.028 m, and S = 0.830 × 10−3).
Water 10 00870 g006
Figure 7. The relation of the bed slope S to the relative depth h/D.
Figure 7. The relation of the bed slope S to the relative depth h/D.
Water 10 00870 g007
Figure 8. The relation of the resistance factor ratio c/c’ to the relative depth h/D.
Figure 8. The relation of the resistance factor ratio c/c’ to the relative depth h/D.
Water 10 00870 g008
Figure 9. The relation of translation parameters (a) Δx and (b) Δy to the bed slope S.
Figure 9. The relation of translation parameters (a) Δx and (b) Δy to the bed slope S.
Water 10 00870 g009
Figure 10. Comparison of the computed discharge (a) Q and (b) Q′ with the measured discharge Qm.
Figure 10. Comparison of the computed discharge (a) Q and (b) Q′ with the measured discharge Qm.
Water 10 00870 g010
Figure 11. Frequency distribution of computed discharge relative error.
Figure 11. Frequency distribution of computed discharge relative error.
Water 10 00870 g011
Figure 12. The relation of the resistance factor ratio c/c′ to the relative depth h/D with different types of D.
Figure 12. The relation of the resistance factor ratio c/c′ to the relative depth h/D with different types of D.
Water 10 00870 g012
Figure 13. The relation of the resistance factor ratio c/c′ to the relative depth h/D with different values of cv.
Figure 13. The relation of the resistance factor ratio c/c′ to the relative depth h/D with different values of cv.
Water 10 00870 g013
Table 1. Data sources and range.
Table 1. Data sources and range.
No.Data SourcesQ (m3/s)h (m)S (×10−3)D (mm)Number of Groups
1da Cunha (1969)29.00∼159.990.55∼1.970.66∼0.952.244
2Einstein (1944)0.10∼0.450.07∼0.211.48∼1.650.929
3Leopold (1969)83.33∼499.300.96∼4.110.04∼0.350.14∼0.8155
4Milhous (1973)1.33∼3.400.37∼0.539.70∼10.808.20∼27.0014
5Simons (1957)1.22∼29.420.80∼2.590.06∼0.330.10∼0.7210
6Toffaleti (1968), Red River190.28∼1537.563.00∼7.380.07∼0.080.09∼0.2230
7Toffaleti (1968), Rio Grande River35.11∼282.310.33∼1.060.74∼0.890.21∼0.3631
Range0.10∼1537.560.07∼7.380.04∼10.800.10∼27.00213 (total)

Share and Cite

MDPI and ACS Style

Dai, W.; Ding, M.; Zhang, H. On the Difference of River Resistance Computation between the k ε Model and the Mixing Length Model. Water 2018, 10, 870. https://doi.org/10.3390/w10070870

AMA Style

Dai W, Ding M, Zhang H. On the Difference of River Resistance Computation between the k ε Model and the Mixing Length Model. Water. 2018; 10(7):870. https://doi.org/10.3390/w10070870

Chicago/Turabian Style

Dai, Wenhong, Mengjiao Ding, and Haitong Zhang. 2018. "On the Difference of River Resistance Computation between the k ε Model and the Mixing Length Model" Water 10, no. 7: 870. https://doi.org/10.3390/w10070870

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