Depth-averaged modeling of free surface flows in open channels with emerged and submerged vegetation

Applied Mathematical Modelling 37 (2013) 540–553

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Depth-averaged modeling of free surface flows in open channels with emerged and submerged vegetation Mingliang Zhang a,⇑, C.W. Li b, Yongming Shen c a

School of Ocean and Environment Engineering, Dalian Ocean University, Dalian, People’s Republic of China Department of Civil and Structure Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong c State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian, People’s Republic of China b

a r t i c l e

i n f o

Article history: Received 24 September 2009 Received in revised form 20 February 2012 Accepted 29 February 2012 Available online 9 March 2012 Keywords: Depth-averaged modeling Emerged and submerged vegetation Manning equation k–e Turbulent model SIMPLEC algorithm

a b s t r a c t The resistance induced by vegetation on the flow in a watercourse should be considered in projects of watercourse management and river restoration. Depth-averaged numerical model is an efficient tool to study this problem. In this study, a depth-averaged model using the finite volume method on a staggered curvilinear grid and the SIMPLEC algorithm for numerical solution is developed for simulating the hydrodynamics of free surface flows in watercourses with vegetation. For the model formulation the vegetation resistance is treated as a momentum sink and represented by a Manning type equation, and turbulence is parameterized by the k–e equations. An analytical equation is derived to represent the resistance induced by submerged vegetation by an equivalent Manning roughness coefficient. Numerical simulation is carried out for the flow in an open channel with a 180° bend, and the flow in a curved open channel partly covered by emerged vegetation, as well as the flow in a straight trapezoidal channel with submerged vegetation. The agreement between the computed results and the measured data is generally good, showing that the resistance due to emerged or submerged vegetation can be represented accurately by the Manning roughness equation. The computed results demonstrate that the depth-averaged modeling is a reasonable and efficient tool to study flows in watercourses with vegetations. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Vegetations, such as grasses, shrubs and mangroves, frequently grow in watercourses, floodplains and wetland waters. They increase the hydraulic resistance to water flow and reduce the shear stresses at the bed of the watercourse. The flow carrying capacity of the watercourse thus will be reduced and the potential of trapping and deposition of sediments will be enhanced. On the contrary, the increased flow resistance in a watercourse leads to an increase in the water depth and a reduction of the flow velocity, which are beneficial for the stabilization of riverbanks and the protection of floodplains. Furthermore, vegetations in rivers and coastal zones are important for the control balance of the ecosystem, as well as the solute transport process. In recent years unprecedented environmental concerns have prompted the need for a better understanding of the hydrodynamics and mixing in vegetated rivers, coastal zones and wetland waters. Knowledge on flow through rigid and flexible vegetation are primarily obtained from the studies of atmospheric flows over plant canopies [1,2], as well as the studies of flows in vegetated open channels [3–9]. One of the major objectives of these studies is to investigate the turbulence structure and the related transport processes in vegetated environments. Furthermore, other studies have been carried out to investigate the ⇑ Corresponding author. E-mail address: [email protected] (M. Zhang). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.02.049

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

541

vegetation induced resistance on fluid flow [10–12]. A Manning type equation can be obtained for the estimation of the roughness for emerged and flexible vegetation in floodplain and vegetated zone of river by Stone [13]. Subsequently, some other empirical equations for evaluating the resistance due to emerged and submerged vegetation in terms of the vegetation parameters (e.g. diameter of stems, number density of vegetation, drag coefficient) have been developed by Fahi and Kouwen [14], Wilson and Horritt [15], Wu et al. [16] and Baptist et al. [17]. In the past several decades, numerical models have been extensively used in simulating and predicting velocity and solute transport fields in open channels, rivers and other water bodies, most models are developed based on the finite volume method, the finite element method, the finite difference method. The models developed for vegetated open-channel flows range from simplified one-dimensional models to fully three-dimensional models. Yoshida and Dittrich [18] presented a 1D steady-state flow simulation model with evaluating the roughness parameters due to the surface roughness of the main channel and floodplains, riparian forests on the floodplains. Nadaoka and Yagi [19] developed large eddy simulation (LES) models to simulate the hydrodynamic behavior of turbulent flow in open channels with domains of vegetation. David and Horritt [20] used TELEMAC-2D to simulate a flood event of the River Severn, UK. Vionnet et al. [21] calculated the values of the flow resistance and eddy viscosity coefficients in channel and floodplain areas using a 1D model based on the lateral distribution method and incorporated the results into a 2D numerical model. Wu and López [22,23] proposed a depth-averaged two-dimensional k–e model for computing the flow and sediment transport in vegetated open channels based on the finite volume method. With the development of computer technology, many 3D turbulent numerical models with vegetation have been developed such as mixing-length model, two-equation turbulence model, more accuracy LES for evaluating the resistance of the vegetation on water flow, the vegetation is considered by as internal source of resistant drag force and is added into the governing equations [24–26]. For the study of flows in watercourses with vegetation of different species and various degree of submergence, a depth-averaged models are enough efficient and of sufficiently computation accuracy. In calculating the vegetated channel or river, it is essential to express the depth-integrated vegetal resistance in terms of the vegetation properties such as height, number density and drag coefficient to obtain a Manning type equation for resistance. In this paper, a 2D depth-averaged model with k–e parameterization of turbulence is developed to compute the flows in open channels with vegetation and without vegetation. The mass conservative finite volume method is used to discretize the equations. Boundary-fitted grids are used to approximate closely the complex boundary geometries of watercourses. The resulting system of coupled algebraic equations of velocities and pressure is solved by the SIMPLEC algorithm. Vegetal resistance is represented by the Manning equation expression. A new analytical equation is derived to convert the drag resistance induced by submerged vegetation into an equivalent Manning roughness coefficient. The model is applied to simulate the flows in curved open channels with or without vegetation, and the flow in a trapezoidal channel with submerged vegetation. 2. Mathematical model 2.1. Basic equations The continuity equation, momentum equations and the k–e equations for turbulence closure on curvilinear coordinates are written as follows:

    @ 1 @ 1 @ 1 @ Hceff 1 @ Hceff ðHqV/Þ ¼ ðHq/Þ þ ðHqU/Þ þ ða/n  b/g Þ þ ðc/g  b/n Þ þ S/ ðn; gÞ; @t J @n J @g J @n J @g J J

ð1Þ

where / stands for 1, u, v, k and e, respectively; (U, V) are the contravariant velocities along the curvilinear coordinates (n, g); H is the water depth; q is the density of water; ceff is the effective viscosity coefficient; S/(n, g) is the source term of variable /, the expressions of ceff and S/ ðn; gÞ are given in Table 1.

         1 @z @z 1 @ Hceff 2 @u @u 1 @ Hceff @v @v þ þ  yn yg qgH yg  yn  xg y g yg yn xg J @n @g J @n @n @g J @n J J @g @n       1 @ Hceff @v @v 1 @ Hceff 2 @u @u þ  sbx ;  yn yg þ xn yg  xn yn yn J @g J @g @g @n J @n @g J

Su ¼ 

Table 1 Variables, diffusive coefficients and source terms of the general equations. Equation

u

ceff

Su

Continuity n -Momentum g-Momentum Turbulent kinetic energy

1 u v k

0

0 Su Sv H(Pk + Pkv  e)

Turbulent dissipation rate

e

leff leff leff rk leff re

H½ðke C z1 P k  C z1 qeÞ þ P ev 

ð2Þ

542

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

         1 @z @z 1 @ Hceff 2 @ v @v 1 @ Hceff @u @u þ þ qgH xn  xg  xg yg xg  xn xg xn y g J @g @n J @n J @n @g @n J @n @g J       H c 1 @ Hceff @u @u 1 @ @v @v eff þ þ  xn xg  xn yn y n xg x2n  sby ; J @g @n @g J @g J J @g @n

Sv ¼ 

U ¼ uyg  v xg ;

V ¼ v xn  uyn ;

ð4Þ

a ¼ x2g þ y2g ; b ¼ xn xg þ yn yg ; c ¼ x2n þ y2n ; J ¼ xn yg  yn xg ; Pk ¼ 2ceff

Pkv ¼

P ev ¼



@u yg @u yn  @n J @g J

2

ð5Þ

 2    2 @ v xg @ v xn @ v yg @ v yn @u xg @u xn þ 2ceff  þ ceff ; þ  þ  þ @n J @g J @n J @g J @n J @g J

ck qU 3 ; H ce qU 4 H2

ð6Þ

ð7Þ

ð8Þ

;

where 2

ceff ¼ qðm þ mt Þ; mt ¼ C l U ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cf ðu2 þ v 2 Þ;

k

e

ð9Þ

;

H ¼ z  zb ;

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

sbx ¼

qgn2 u u2 þ v 2

1 ck ¼ pffiffiffiffi ; cf

H1=3 ce ¼

ð10Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

;

sby ¼

3:6C e2 qffiffiffiffiffiffi Cl; c0:75 f

qgn2 v u2 þ v 2 H1=3

ð11Þ

;

cf ¼ 0:003;

ð12Þ

where (u, v) are the depth-averaged velocities in (x, y) directions, k is the turbulent kinetic energy, e is the rate of dissipation of turbulent kinetic energy, z is the water level, zb is the bed elevation, mt and m are the turbulent eddy viscosity and molecule viscosity coefficient respectively, Pk is the source of turbulent kinetic energy, g is the gravitational acceleration, U⁄ is the friction velocity, sbx and sby are the x- and y-components of the bed shear stress, n is the Manning’s roughness coefficient. a, b and c are the metric tensor components, J is the Jacobian matrix. Cl, C e1 ; C e2 ; rk and re are turbulent constants with values 0.09, 1.44, 1.92, 1.0 and 1.3 respectively. 2.2. Resistance due to emerged vegetation The effect of the vegetation on the flow can be simulated by an internal source of resistant force per unit fluid mass added into the momentum equations. Force equilibrium analysis can also be performed to obtain an equivalent Manning’s roughness coefficient for the vegetation resistance which is dependent on the flow depths, number density and the diameter of the vegetation elements. Fig. 1 shows the vegetation elements under submerged and emerged conditions in an open channel. For a group of emerged vegetation elements, the total shear stress, sv consists of the bed shear stress, sb and the vegetation resistance force per unit horizontal area mFd, where m is the number density (the number of vegetation elements per

Vegetation

Usub U

hv

H

(a) Emergent

H

Uveg

uu hv

uc

(b) Submerged

Fig. 1. Vegetation elements in open channel.

Vegetation

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

543

unit horizontal area) and Fd is the vegetation resistance force for each vegetation element. The equilibrium condition can be expressed by

sv ¼ sb þ mF d :

ð13Þ

Several investigators proposed a sheltering factor for the vegetation resistant force to account for the effect of velocity defect behind the vegetation elements [24], however the sheltering factor needs to be determined experimentally and is dependent on many factors, such as the size of the vegetation element, the number density, etc. For simplicity, the sheltering effect is generally absorbed into the drag coefficient of the vegetation element and the resistance force due to an emerged vegetation element is given by

1 qC D DHU 2v ; 2

Fd ¼

ð14Þ

where CD is the drag coefficient and Uv is the velocity against the vegetation element, D is the averaged diameter of the vegetation element. The number density is given by

1

m¼

2

dv

ð15Þ

;

where dv is the average distance between two adjacent vegetation elements. For emerged vegetation, Uv is the depth-averaged velocity U [13]. The local total shear stress can be related to the surface slope s by

sv ¼ qgHs:

ð16Þ

Applying the Manning equation to the bed shear stress,

qgn20 U 2

sb ¼

H1=3

ð17Þ

;

where n0 is the Manning coefficient for the bed roughness, the depth-averaged velocity can be obtained by substituting Eqs. (14), (16), and (17) into Eq. (13)

U¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hs n2b =H1=3 þ ðC D mDhv Þ=2g

ð18Þ

:

Applying the Manning equation to the total shear stress,

sv ¼

qgn2v U 2 H1=3

ð19Þ

;

where nm is the Manning coefficient for the total roughness, the following equation can be obtained

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C D mDH4=3 nv ¼ n20 þ : 2g

ð20Þ

For channels with densely distributed vegetation, the vegetation resistance will be dominant and n0 in Eq. (20) can be neglected. 2.3. Resistance due to submerged vegetation There are two distinctive flow regions in a flow with submerged vegetation: the region within the vegetation field in which the flow is subjected to the vegetation resistance, and the region above the vegetation field in which the flow is subjected to the shear resistance at the interface of the two regions. The average velocity of the flow within the vegetation field Uveg can be determined by the force equilibrium equation (13)

U v eg ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hs n2b =H1=3 þ ðC D mDhv Þ=2g

;

ð21Þ

where hv is the height of the vegetation elements. The vertical variation of the horizontal velocity over the vegetation field is due to the shear force at the interface of the vegetation field and the free water body. This vertical velocity profile has been studied extensively and is simulated by modified log-law profiles in which a zero-plane displacement parameter and a vegetation-dependent roughness parameter are introduced. Numerous expressions of similar forms have been proposed by Haber [27] and Baptist et al. [17], because the Haber’s expression can lead to a simpler expression of the equivalent Manning coefficient, the Haber’s expression is used in the present study,

544

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

uu ðzÞ ¼

u

j

ln

z

k

þ uc ;

ð22Þ

where uc is the ‘slip velocity’ to match the two profiles at z = k, j is the von Karman’s constant, and u⁄ is the shear velocity given by

u ¼

pffiffiffiffiffiffiffiffi gHs:

ð23Þ

The depth averaged velocity above the vegetation is given by

U sub ¼

1 hk

Z

h

uu dz ¼

u

k

j



   h h  1 þ uc : ln hk k

ð24Þ

Assume uc = aUveg, where a P 1. The discharge per unit width is

q ¼ hv U v eg þ ðH  hv ÞU sub ¼ HU v eg þ ðH  hv ÞðU sub  uc þ ða  1ÞU v eg Þ:

ð25Þ

Using the Manning equation

q¼

1 5=3 pffiffi H s: nr

ð26Þ

The representative Manning coefficient nr is then given by

8 91 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi >   > pffiffiffi    
þ 1=6 1 þ ða  1Þ 1   1 nr ¼ t ln : 1=3 > 2 h H j H > v H : nb þ C D mDhv H ; =2g

ð27Þ

CD is a function of the velocity, diameter of the vegetation element and fluid viscosity. For rigid cylinder, the experimental drag coefficient CD is dependent on the Reynolds number and is approximated by the following equations

C D ¼ ð103 =Re Þ0:25

for Re 6 103 ;

ð28Þ

C D ¼ Minð0:79 þ ½ð103 Re  2Þ=20:52 ; 1:15Þ for 103 6 Re 6 4  103 ; Re ¼

UD

m

ð29Þ ð30Þ

:

For flexible vegetation, the drag coefficient should be determined experimentally. 2.4. Numerical method The governing equations are discretized by the finite volume method on a curvilinear non-orthogonal coordinate system. The SIMPLEC algorithm is used to solve the resulting system of algebraic equations. In particular, the convection terms are discretized by a hybrid scheme, and the diffusion terms are discretized by the central difference scheme and the time derivative term is discretized by the first-order implicit scheme. To prevent the occurrence of inter-node oscillation, a partly staggered grid system is used. Special treatment of the discretized continuity equation has been done to obtain an improved flow field. The boundary conditions are specified as follow. At the inlet, the discharge or the velocity should be given. The turbulent kinetic energy and the dissipation rate of turbulent kinetic energy are then specified.

u ¼ u0 ;

v ¼ v0;

k ¼ 0:00375ðu2 þ v 2 Þ;

3

e ¼ 0:09k2 =ð0:05HÞ:

ð31Þ

At the outlet, the water level should be provided, assuming the fully developed condition is attained at the outlet, the normal gradients of the rest quantities are set to zero, namely,

@u @ v @k @ e ¼ ¼ ¼ ¼ 0: @n @n @n @n

ð32Þ

At the side wall boundary, the wall function is employed. The Manning equation is used to account for the bed roughness. 3.. Model applications 3.1. Curved open channel The model is first applied to the case of flow in an open channel with a 180° curved bend for which experimental data are available [28]. The channel is of rectangular cross section and of width 1.7 m. It has a 180° bend with 6 m long straight

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

545

sections before and after the bend. The discharge is 0.19 m3/s and the downstream water depth H = 0.18 m, giving a Froude number of 0.47. To replicate the channel numerically a curvilinear grid system of 101  26 points is used (Fig. 2). The computation stops when the convergence error is less than 105. Fig. 3 gives the depth-averaged velocity distribution at different sections of the channel. It can be seen that the agreement between the measurement and the computed results is generally good for different sections. Fig. 4 shows the comparison of the simulated and measured surface elevations at three sections: the center line section, the section at B/10 from the outer bank and section at B/10 from the inner bank. The longitudinal distance from the inlet is normalized by the length at the center of the channel and denoted by S; and the water surface depression is normalized by the downstream water depth and denoted by H: Generally there is a close agreement between the measured data and computed results at the three sections, although there are some discrepancies along the inner bank region. The water surface elevation field is shown in Fig. 5, it is noted that the water surface is higher at the outer bank due to the centrifugal force at the curved channel section. Overall, from these figures, one can see that this present model has capability of predicting complex flow in the open channel. 3.2. Curved open channel with emerged rigid vegetation

180

The model is then applied to simulate the flow in a curved open channel partly covered by rigid vegetation. The corresponding physical experiments have been carried out by Tominaga et al. [5], the channel cross section is rectangular and the channel is 17.2 m long, 0.9 m wide and 0.3 m deep. The channel has a 60 bend with 10.8 m and 3.6 m long straight sections before and after the bend. The radius of curvature of the center line of the bend is 2.7 m, as shown in Fig. 6. Vegetations are simulated by wooden sticks of 5 mm diameter and with vegetation density 400 per unit area. The vegetation zones are of width 225 mm and variable length and are placed at various locations at the inside and outside banks of the channel as shown in Fig. 7 and Table 2. The discharge is 0.04 m3/s with a constant downstream water depth H = 0.15 m. The cases of V-5, V-6 and V-7 are adopted to test interaction between flow and vegetation by using the numerical model, the case of R-5 is a reference case in which the channel has no vegetation. The measurement sections are located at h = 0°, 30°, 60° and at x1 = 135 cm. h is the angle measured clockwise from the y-axis, x and x1 are the streamwise coordinates measured from the bend inlet (h = 0°) and from the bend outlet (h = 60°), respectively. The computational mesh consist of 121  46 quadrilateral cells in longitudinal and transverse directions with a smallest cell area of 0.0009 m2, and the finer mesh is adopted at the contact region between the vegetation zone and curved channel section. The computational time step is 1 s, the computation stops until the steady state is reached, the convergence error bound is set at 105, the total simulating time is less than 1 min on a PC with AMD Athlon processor (2.4 GHz) for all cases. Fig. 8 gives the lateral distributions of the depth-averaged velocity with different grid points in transverse direction for the case of double-bank arrangement (V-7). The grid resolution (20, 46 and 60 grid points in transverse direction) effect on numerical results is presented in Fig. 8. From this figure, the grid with twenty grid points in transverse direction is very coarse and in turn, the simulated velocity profile is not accurate enough in these selected sections, especially near wall boundary. The calculated velocity profiles with 46 and 60 grid points almost match together, from the comparison, the 46 grid point used is enough to get the good computed accuracy and high computed efficiency for this case study. For case

0 12

90 60

30

0 Fig. 2. Grid system of the open channel.

546

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

Measured

Calculated

0.9

0°

30°

60°

90°

120°

180°

Velocity (m s-1)

0.6 0.3 0 0.9 0.6 0.3 0 0.9 0.6 0.3 0

0

0.2

0.4

0.6

0.8

0

0.2

0.4

0.6

0.8

1 y/B

Fig. 3. Depth-averaged velocity profiles at different sections.

Right bank

Center

Left bank

Cal

0.1

H

0.07 0.04 0.01 -0.02 0

0.2

0.4

0.6

0.8

1

S Fig. 4. Water surface elevation profiles along the channel.

10

0.1

0.18

893 92 3 28 3

8 0.18786 4

y (m)

6

4

6 44 92 1 0.1 641 9 0.1

2

0

0

2

4

6

8

10

12

x (m) Fig. 5. Computed water surface elevation field.

V-7 (Fig. 8), the relative maximum velocity is the largest at section 3 and section 4, the velocity is significantly decelerated and a high lateral shear is generated at the upstream vegetated region and the downstream outer-bank region due to vegetation effect on flow.

547

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

y 10.8m

R=2.25m

y1

0.9m

x

m

x1

3.6

Fig. 6. Sketch of the experimental flume.

V-5

V-6

V-7

Fig. 7. Plan view of the vegetation arrangements.

Table 2 Vegetation arrangements [5]. Case

Vegetation arrangements

R-5 V-5 V-6 V-7

Non-vegetation Inside Outside Inside and outside

x = 90 cm to h ¼ 15 / ¼ 45 to x1 = 90cm x = 90 cm to h ¼ 15 þ h ¼ 45 to x1 = 90 cm

Figs. 9 and 10 give the lateral distributions of the depth-averaged velocity in the case of short inner-bank arrangement (V-5) and short outer-bank arrangement (V-6) respectively, the simulation results and measurement data match well. For case V-5, the computed results show that the existence of vegetation causes the redistribution of the velocity field, the velocity in the vegetated domain is reduced by the presence of vegetation and is smaller than that of case R-5. On the contrary, the velocity in the non-vegetated zone is larger than that of case R-5. In the computation a large velocity gradient forms around the interface region between the non-vegetated and vegetated domains after the computation reaches a dynamic equilibrium state. For case V-6, because the measured sections 1 and 2 are far from the vegetated zone, the velocities at sections 1 and 2 are unaffected by the vegetal resistance (Fig. 10), a sharp lateral gradient of velocity appears at section 3 and maintains at the downstream straight region (section 4). The calculated turbulent kinetic energy and turbulent dissipation rate are presented in Fig. 11. It can be seen that the maximum turbulent kinetic energy and turbulent dissipation rate are generated at the interfacial region of the vegetated

548

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

Measurement

Section 1

grid20

2

grid20

2

grid46

grid46

grid60

1.6

grid60

1.2

u/u0

u/u0

1.6

0.8

1.2 0.8

0.4

0.4

0 0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

y (m) Measurement

Section 3

0.8

1

Section 4

Measurement grid20

2

grid46 grid60

1.6

0.6 y (m)

grid20

2

grid46 grid60

1.6

1.2

u/u0

u/u0

Section 2

Measurement

0.8

1.2 0.8

0.4

0.4

0

0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

y (m)

0.6

0.8

1

y (m)

Fig. 8. Velocity profiles with different grid points in transverse direction for case (V-7).

2

2

1.6

1.6

1.2

1.2

u/u 0

u/u 0

Section 1

0.8 0.4

Section 2

0.8 0.4

0 0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0

0.9

0

2

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.8

0.9

2 Section 3

Section 4

1.6

1.5

1.2

u/u 0

u/u 0

0.7

0.8

1 0.5

0.4 0

0 0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

Measurement

0.9

0 Non-vegetation

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

Model

Fig. 9. Lateral distributions of depth-averaged velocity in case V-5.

and non-vegetation zones, the major source of turbulence is due to the sharp velocity gradient occurred at the interfacial region. Fig. 12 displays the computed Reynolds stress, which clearly shows that the Reynolds stress is the maximum at the interfacial region. Figs. 13 and 14 show the contours of the computed water level field and the computed flow pattern for the case of short inner-bank arrangement (V-5), respectively. Fig. 15 and 16 show the contours of computed water level field and the computed flow pattern for the case of double-bank arrangement (V-7), respectively. From Figs. 13 and 15, it is found that the upstream water level of the vegetation zone is high due to the obstruction of the vegetation on the water flow. A local zone

549

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

2

2 1.6

1.2

u/u 0

u/u 0

Section 2

Section 1

1.6

0.8 0.4

1.2 0.8 0.4

0

0 0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

0

2

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.8

0.9

2 Section 3

Section 4

1.6

1.5

1.2

u/ u 0

u/u 0

0.7

0.8

1 0.5

0.4 0

0 0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

Measurement

0

0.1

0.2

Non-vegetation

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

Model

Fig. 10. Lateral distributions of depth-averaged velocity in V-6.

(a) 0.003

(b) 0.003

0.002

0.002

0.001

0.001

0 0

(c)

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

0 0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

0.8

0.9

(d) 0.004

0.003

0.003 0.002 0.002 0.001

0

0.001

0

0.1

0.2

0.3

0.4 0.5 y (m)

0.6

0.7

Turbulent kinetic energy

0.8

0.9

0

Turbulent dissipation rate

Fig. 11. Calculated turbulent kinetic energy and turbulent dissipation. (a) Section 1 of V-5 case, (b) section 3 of V-6 case, (c) section 1 of V-7 case and (d) section 3 of V-7 case.

of water level depression occurs downstream the vegetation zone, and the water level increases again further downstream. Due to the non-uniform velocity distributions over the water depth and the centrifugal force in the curved open channel section, the water level at the outer bank is higher than that at the inner bank. From Figs. 14 and 16, the velocity reduction effect by vegetation resistance is clearly observed. Additional numerical simulations have been carried out to evaluate the effect of vegetation density on water flow. The case of study is V-7 with three sub-cases of different values of vegetation density: m = 400/m2, 1111/m2, 12,345/m2. The computed center streamlines are shown in Fig. 17, they are deflected away from the vegetation zone towards the outer bank until they encounter the second vegetation zone downstream where they are deflected again toward the opposite bank. As expected, when the vegetation density increases, the friction increases and a larger deflection of streamline results.

550

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

0.149 0

0.1501

3

0.1503

Fig. 12. Calculated Reynolds stress u0 v 0 .

0.1480 0.1470

0.1458 0. 1 461

y (m)

2

1

62 0.14 0.1 46 1 0. 14 60

0

-1 0

2

4

6

x (m) Fig. 13. Contours of simulated water level in case (V-5).

3

y (m)

2

1

0

-1 0

2

4

6

x (m) Fig. 14. Simulated flow pattern in case (V-5).

3.3. Trapezoidal open channel with submerged rigid vegetation The last case of simulation is flow in a trapezoidal channel with submerged vegetation. A physical model study has been performed in a straight and trapezoid channel formed in a 1.52-m-wide rectangular Plexiglas flume [29]. The length of the test section in the flume is 4.9 m, the channel section consists of two banks, one is vertical and the another has a 2H:1V side slope, the largest width of the channel is 76.2 cm. The bank is covered by uniform sands with a median grain size d50 of 0.33 mm. To simulate willow posts, wood dowels with a diameter of 0.64 cm are placed on the bank with a staggered pattern, there are three rows of dowels and the above ground height of each dowel is 10.2 cm. The discharge is 0.0885 m3/s with

24 0.15

3

0. 15 25

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

0.1500 0.14 80 0.1 475

2

y (m)

551

465 0.1

1

0.14 59 0.1460

0

0.1461

-1 0

2

4

6

x (m) Fig. 15. Contours of simulated water level in case (V-7).

3

y (m)

2

1

0

-1 0

2

4

6

x (m) Fig. 16. Simulated flow pattern in case (V-7).

Dashdotdotm=400 Solidm=1111 Dashdotm=12345

Fig. 17. Simulated streamline with different vegetation density for case (V-7).

a flow depth 0.381 m. The computational mesh consisted of 210  70 grid node locations, because the domain is simple, a uniform mesh is used here, with a size of 0.05 m  0.02 m for each rectangle cell. The computational time step is 1 s, the computation stops until the steady state is reached.

552

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

u (m s-1 )

0.4

Data Model Baptist equation

0.3 0.2 0.1 0

0

0.3

0.6

0.9

1.2

1.5

y (m) Fig. 18. Predicted and observed depth-averaged velocity in test 7.

1.5

1.2

Dowels zone

y (m)

0.9

0.6

0.3

0

5

5.5

6

6.5

7

x (m) Fig. 19. Simulated flow pattern in test 7.

0.4 Data Model Baptist equation

u (m s-1)

0.3 0.2 0.1 0

0

0.3

0.6

0.9

1.2

1.5

y (m) Fig. 20. Predicted and observed depth-averaged velocity in test 5.

The computed velocities are shown in Figs. 18 and 19. It can be seen that the present model predicts the velocities quite accurately as compared with the measured data. The results also show that the velocities decrease in the vicinity of the dowels where the head loss is the maximum. The equation for the equivalent Manning roughness due to submerged vegetation derived in this study, Eq. (27), is slightly more accurate than that due to Baptist et al. [17] for this case. One more case of simulation is performed to further validate the model. In this case five rows of dowels were placed in a staggered pattern. The two rows of dowels closest to the toe of the slope are 5.1 cm high and the remaining rows of dowels are 10.2 cm high. The discharge, inlet velocity and water depth are 0.0885 m3/s, 0.23 m/s and 0.381 m, respectively. The computed transverse distribution of the streamwise depth-averaged velocity is shown in Fig. 20. Again the computed results are in good agreement with the measured data. Compared to the case above, the velocity reduction is more significant at the dowel-placing region due to the increase in the dowel placing area. 4. Conclusions A 2D depth-averaged k–e model for flows in vegetated watercourses with complex geometries has been developed. The mass conservative finite volume method with curvilinear coordinates is used to discretize the equations. The SIMPLEC

M. Zhang et al. / Applied Mathematical Modelling 37 (2013) 540–553

553

algorithm is used to solve the resulting coupled algebraic equations of velocities and pressure (water surface elevation). Hydraulic resistance due to submerged or emerged vegetation is represented by a Manning type expression. The model is applied to simulate the flows in curved open channels with or without vegetation, and the flow in a trapezoidal channel with submerged vegetation. The model predicts correctly the velocity field and the water level distribution for the flow in an open channel with curved bend. It reproduces accurately the vegetal resistant effect on the velocity profiles at various sections of a 60° bended channel, and predicts the occurrence of the peak turbulence kinetic energy and Reynolds stress at the interface of the vegetated and non-vegetated region. The newly derived expression of the equivalent Manning roughness estimates accurately the resistance due to submerged vegetation in the cases of flow in a trapezoidal channel. The model has been proven to be an efficient and accurate tool to study flows in watercourses with vegetation. Acknowledgments This study is funded by Special Funds for Talent Projects of Dalian Ocean University No. 017211, Special Funds for Postdoctoral Innovative Projects of Liaoning Province No. 2011921018, the Research Grant Council of the Hong Kong Special Administrative Region No. 5221/06E and by the National Basic Research (973) Program of China No. 2006CB403302. The authors are thankful the comments from anonymous reviewers, which significantly improved the paper. References [1] G.G. Katul, L. Mahrt, D. Poggi, C. Sanz, One and two equation models for canopy turbulence, Bound. Lay. Meteor. 113 (2004) 811–819. [2] D. Poggi, A. Porporato, L. Ridolfi, J.D. Albertson, J.D. Katul, The effect of vegetation density on canopy sub-layer turbulence, Bound. Lay. Meteor. 111 (2004) 565–587. [3] D. Naot, I. Nezu, H. Nakagawa, Hydrodynamic behavior of partly vegetated open channels, ASCE J. Hydraul. Eng. 122 (11) (1996) 625–633. [4] H.M. Nepf, Drag, turbulence, and diffusion in flow through emergent vegetation, Water Resour. Res. 35 (2) (1999) 479–489. [5] A. Tominaga, M. Nagao, I. Nezu, Flow structures and momentum transport processes in curved open channels with vegetation, in: 28th IAHR Congress, 1999, Graz, Austria. [6] H.M. Nepf, E.R. Vivoni, Flow structure in depth-limited, vegetated flow, J. Geophys. Res. 105 (C12) (2000) 28547–28557. [7] C.I. Thornton, S.R. Abt, C.E. Morris, J.C. Fischenich, Calculating shear stress at channel overbank interfaces in straight channels with vegetated floodplains, ASCE J. Hydraul. Eng. 26 (12) (2001) 929–936. [8] I. Nezu, O. Kouki, Turbulent structures in partly vegetated open-channel flows with LDA and PIV measurements, J. Hydraul. Res. 39 (6) (2001) 629–642. [9] S.J. Bennett, Taner Pirim, B.D. Barkdoll, Using simulated emergent vegetation to alter stream flow direction within a straight experimental channel, Geomorphology 44 (2002) 115–126. [10] N. Kouwen, T.E. Unny, H.M. Hill, Flow retardance in vegetated channels, J Irrigat. Drain Div. 95 (2) (1969) 329–342. [11] S. Petryk, I.I.IG. Bosmajian, Analysis of flow through vegetation, J. Hydrol. Div. 101 (7) (1975) 871–884. [12] V. Kutija, H.T.M. Hong, A numerical model for assessing the additional resistance to flow introduced by flexible vegetation, J. Hydraul. Res. 34 (1) (1996) 99–114. [13] A.M. Stone, H.T. Shen, Hydraulic resistance of flow in channels with cylindrical roughness, ASCE J. Hydraul. Eng. 128 (5) (2002) 500–506. [14] M.M. Fahi, N. Kouwen, Nongrid, nonsubmerged vegetative roughness on Floodplains, ASCE J. Hydraul. Eng. 123 (1) (1997) 51–57. [15] C.A.M.E. Wilson, M.S. Horritt, Measuring the flow resistance of submerged grass, Hydrol. Proc. 16 (13) (2002) 2589–2598. [16] F.C. Wu, H.W. Shen, Y.J. Chou, Variation of roughness coefficients for unsubmerged and submerged vegetation, ASCE J. Hydraul. Eng. 125 (9) (1999) 934–942. [17] M.J. Baptist, V. Babovic, J.R. Uthurburu, M. Keijzer, R.E. Uittenbogaard, A. Mynett, A. Verwey, On inducing equations for vegetation resistance, J. Hydraul. Res. 45 (4) (2007) 435–450. [18] H. Yoshida, A. Dittrich, 1D unsteady-state flow simulation of a section of the upper Rhine, J. Hydrol. 269 (2002) 79–88. [19] K. Nadaoka, H. Yagi, Shallow-water turbulence modelling and horizontal large-eddy computation of river flow, ASCE J. Hydraul. Eng. 124 (1998) 493– 500. [20] C.M. David, M.S. Horritt, Floodplain friction parameterization in two-dimensional river flood models using vegetation heights derived from airborne scanning laser altimetry, Hydrol. Proc. 17 (2003) 1711–1732. [21] C.A. Vionnet, P.A. Tassi, M.J.P. Vide, Estimates of flow resistance and eddy viscosity coefficients for 2D modelling on vegetated floodplains, Hydrol. Proc. 18 (2004) 2907–2926. [22] W. Wu, F.D. Shields Jr., S.J. Bennett, S.S.Y. Wang, A depth averaged two-dimensional model for flow, sediment transport, and bed topography in curved channels with riparian vegetation, Water Resour. Res. 41 (2005) 1–15. [23] F. López, M.H. García, Mean flow and turbulence structure of open-channel flow through non-emergent vegetation, ASCE J. Hydraul. Eng. 127 (5) (2001) 392–402. [24] X.H. Su, C.W. Li, Large eddy simulation of free surface turbulent flow in partly vegetated open channels, Int. J. Numer. Mech. Fluids 39 (2002) 919–937. [25] M.L. Zhang, C.W. Li, Y.M. Shen, A 3D non-linear k–e turbulent model for prediction of flow and mass transport in channel with vegetation, Appl. Math. Modell. 34 (2010) 1021–1031. [26] B.N. Rajani, A. Kandasamy, Sekhar Majumdar, Numerical simulation of laminar flow past a circular cylinder, Appl. Math. Model. 33 (2009) 1228–1247. [27] B. Haber, Über den Erosionsbeginn bei der Überströmung von flexiblen Rauheitselementen. Mitteilungen des Leichtweiss Institutes für Wasserbau der Technischen Universität Braunschweig, Deft 74, 1982. [28] Y. Jian, J.A. McCorquodale, Depth-averaged hydrodynamic model in curvilinear collocated grid, ASCE J. Hydraul. Eng. 123 (5) (1997) 380–388. [29] V.W. Gregory, Flow through trapezoidal and rectangular channels with rigid cylinders, ASCE J. Hydraul. Eng. 133 (5) (2007) 521–533.