以二維有限元素形貌動力模式探討渠道強制效應對自由沙洲之影響

國立台灣大學生物環境系統工程學研究所 碩士論文

Graduate Institute of Bioenvironmental Systems Engineering

National Taiwan University Master Thesis

以二維有限元素形貌動力模式探討渠道強制效應 對自由沙洲之影響

Investigating the Influences of Channel Forcing Effect on Free Bars Using a 2D Finite-Element Morphodynamic

Model

邵允銓

Yun-Chuan Shao

指導教授：吳富春 Fu-Chun Wu 教授

中華民國九十七年六月 June 2008

b  

c 謝誌

研究所求學的日子裡衷心地感謝吳富春老師的諄諄教誨，引領我走入學術的 殿堂，開拓了我的視野與想法，不僅在學術研究上給與莫大的支持與指導，在生 活與待人處事的細節上亦受惠老師良多，研究所雖僅短短數年，但受吳老師教化 所學的是一生受用的。感謝卡艾瑋老師總是大方地的分享研究經驗與成果，讓我 在實驗的技術與執行上更加的靈活，希望能把論文與卡老師分享是驅使我以英文 撰寫原動力之一。感謝顏清連老師在數值模擬與物理意義之間討論給與我許多寶 貴的意見，使我的論文更為的完整。感謝楊德良老師對於數值方法上的建議與指 導讓我對數值領域有更深的興趣與探索。

生活在研究室的日子裡，要特別感謝陳嘉欣小姐替我們處理許多庶務，為我 們分憂解勞，讓我們無後顧之憂的專心在研究上。感謝葉子豪學長、陳臻琪學長 幫我解決了無數的大大小小的問題。謝謝潘泰斗、顧尚真、林鈺翔、陳昱辰等眾 多學弟一起分攤研究室的人力工作。感謝助理石武融、李依珊、吳子青一起為研 究室分憂解勞。感謝卡老師門下賴悅仁學長、柯文韜學長、粟群超、紀泰安、許 珮蓁同學大方分享研究設備與上課筆記，與你們互相學習交流讓我學到很多。謝 謝黃國文技師與水工所翁總管於流力實驗與水槽實驗上的幫忙，特別感謝莊詔傑 學長與何嵩晟學長分享許多人生的歷練與經驗，讓我於人生的路上有更明確的目 標。

感謝許兄中駿、謝兄清智、甘兄智仁、許兄詩涵、王兄翔生、范兄振國、塗 兄宗杰、塗氏俞慧、周兄宏泰還有郁潔老師，陪我一起互相打氣，打打鬧鬧的撐 過許多日子，還多了一個乾女兒范家甄，有你們真好！感協羅多倫台大分店的大 姐們總是讓我的早餐與晚餐 VIP 升級，讓我能有更充足的補給。

最後要特別感謝我的家人對我的包容與支持，讓可以全心全意的把心思放在 學校，謝謝你們的付出與犧牲。

d Contents

摘要 i

Abstrac ii

List of Tables iii

List of Figures iv

List of Symbols viii

Chap 1 Introduction

1.1 Statement of Problem 1-1

1.2 Literature Review 1-2

1.2.1 Study of free bars 1-2

1.2.2 Study of forced bars 1-5

1.2.3 Free bars affected by channel geometry 1-9

1.2.4 Morphodynamic model 1-11

1.3 Scope of Study 1-13

Chap 2 Mathematical Model

2.1 Governing Equations of Hydrodynamic Model 2-1 2.2 Governing Equations of Bed Evolution Model 2-3

2.2.1 Closure relations 2-3

2.3 Finite Element Method 2-5

2.3.1 Streamline upwind Petrov-Galerkin scheme 2-5 2.3.2 Applying CDG scheme to bed evolution model 2-9

2.4 Model Implementation 2-11

e Chap 3 Model Validation

3.1 Validation of Hydrodynamics Model 3-1

3.1.1 Channel with variable width 3-1

3.2 Validation of Bed Evolution Model 3-6

3.2.1 Forced bars – side bars 3-6

3.2.2 Forced bars – central bars 3-10 3.2.3 Free migrating alternate bars 3-13 3.2.4 Coexistence of free and forced bars 3-18

Chap 4 Forcing Effect of Width Variation on Free Bars

4.1 Numerical Experiments 4-1

4.2 Numerical Results 4-4

4.2.1 Coexistence of forced and free bars 4-4 4.2.2 Evolution of free bars in channels with variable width 4-6 4.2.3 Quantitative forcing effect on equilibrium stage 4-17

Chap 5 Conclusions

5.1 Conclusions 5-1

5.1.1 2D morphodynamic model 5-1

5.1.2 Influences of forcing effect on free bars 5-2

5.2 Suggestions 5-4

Reference R-1

i 摘 要

本研究利用二維有限元素河川形貌動力模式，探討自由沙洲受變寬渠強制效 應之影響。本研究所發展數值模式具有兩個特色：第一是使用流線曲率項來修正 底床拖移載受到水流二次流影響後的運動方向。第二是以流線上風演算法處理泥 沙連續方程式，使數值模式具有模擬強制沙洲(兩側沙洲、中央沙洲)與自由沙洲 共存之能力。數值模擬發現強制沙洲與自由沙洲共存狀態屬於疊加，因此可將自 由沙洲之演變取出討論。研究結果顯示自由沙洲之發展到達平衡狀態時，其波 高、波長與波速受到渠寬變化所影響，會產生穩定的周期性波動，其週期平均值 會隨變寬渠之振幅與波數增大而減小，本研究進一步將變寬渠之振幅與波數整合 為一強制因子，對自由沙洲之影響效應進行量化分析，證明自由沙洲會受渠寬變 化影響而被壓抑。

關鍵字：形貌動力模式、有限元素法、流線上風演算法、自由沙洲、強制效應。

ii Abstract

In this study a two-dimensional finite-element (2D FE) morphodynamic model is developed to investigate the forcing effects of periodic width variation on free bars.

Two features of the proposed model include: (1) The streamline curvature are used to correct the bed load direction effected by secondary flows. (2) The streamline upwind Petrov-Galerkin (SUPG) scheme is applied to solve the sediment continuity equation, which makes it possible to simulate the coexistence of free bars and forced bars (such as side bars and central bars). It is found that the coexistence of free and forced bars is a superposition of the two types of bedform. Thus the free-bar component can be extracted for our study. The results reveal that the bar height, wavelength and celerity of free bars affected by the effect of width variation lead to a periodic wavy pattern when the development of free bars reach the equilibrium state. The mean components of free bars in a cycle of channels with variable-width are inverse proportioned to the amplitude and wave number of width variation. We further derive a forcing factor by combining the amplitude with wave number of width variation and quantitatively prove that the free bars are suppressed by the forcing factor of channels with variable width.

Keywords: Morphodynamic model, finite-element method, streamline-upwind Patrov-Galerkin method, free bars, forcing effect.

iii

List of Tables

Table 4-1 Run number and simulation conditions used in B15 numerical experiments

iv

List of Figures

Figure 1-1 (a) Free alternate bars in Shi-hu Creek, Taiwan; (b) Central bars in Tai-ping Creek, Taiwan.

Figure 1-2 Sketch of alternate-bar structure. [Colombini et al., 1987]

Figure 1-3 A typical neutral curve for alternate bar formation. [Colombini et al., 1987]

Figure 1-4 Free bars wavelength and height as a function of bar celerity in Defina’s numerical experiments. [Defina, 2003]

Figure 1-5 Contour plots showing (a) central bars observed in experiment [Wu and Yeh, 2005] and results of (b) 2-D-Cs model [Wu and Yeh, 2005], (c) 2-D model [Bittner, 1994], (d) 2-D model [Repetto et al., 2002].

Figure 1-6 Contour plots showing (a) side bars observed in experiment [Bittner, 1994] and results of (b) 2-D-Cs model [Wu and Yeh, 2005], (c) 2-D model [Bittner, 1994], (d) 2-D model [Repetto et al., 2002].

Figure 1-7 Comparisons between observed values of beta and corresponding values of BetaC1 and BetaC2. [Wu and Yeh, 2005]

Figure 1-8 Unsteady migration of bars through meanders. Open squares indicate the position of leading edge of bars attached to the left bank; solid squares, the leading edge of bars attached to the right bank.

Figure 1-9 Comparison between measured values of (a) the amplitude and (b) dimensionless wavenumber of the leading Fourier component of bed topography associated with alternate bars in constant width and variable width experiments. [Lanzoni and Tubino, 2001]

v

Figure 3-1 Measured results of (a) bed deformation, and (b) flow depth in Run C1-11 [Bittner, 1994].

Figure 3-2 Averaged bed deformation of Run C1-11 [Bittner, 1994] used as the fixed-bed topography in the validation of the flow dynamics models.

Figure 3-3 (a) Computation domain of Run C1-11 [Bittner, 1994] used in the validation of the hydrodynamics models; (b) zoomed-in element mesh in a cycle of width variation, the wavelength of a cycle (π) is 1.6m.

Figure 3-4 Computed Run C1-11 [Bittner, 1994] results of flow depth by the VA models.

Figure 3-5 Comparison of measured and computed Run C1-11’s [Bittner, 1994]

results of flow depth at four specified sections of a cycle of width variation.

Figure 3-6 (a) Computation domain of Run C1-11 [Bittner, 1994] used in the validation of the bed evolution model; (b) zoomed-in element mesh in a cycle of width variation, the wavelength of cycle (π) is 1.6m.

Figure 3-7 Experimental results of Run C1-11 [Bittner, 1994] are compared with the linear solution and computed numerical bed topography of side bars.

Figure 3-8 Comparison of measured and computed Run C1-11’s [Bittner, 1994]

results of bed topography at four specified sections of a width-variation cycle.

Figure 3-9 (a) Computation domain of Run S6 [Wu and Yeh, 2005] used in the bed evolution model; (b) Zoomed-in element mesh in a cycle of width variation, the wavelength of a cycle is 3.351m.

Figure 3-10 Experimental result of Run S6 [Wu and Yeh, 2005] are compared with the linear solution and computed numerical bed topography of central bars.

Figure 3-11 Experimental result and computed lateral bed profiles of Run S6 [Wu and Yeh, 2005] at the wide and narrow sections.

vi

Figure 3-12 (a) Computation domain of Run P1505 [Lanzoni, 2000] used in bed evolution model; (b) zoomed-in element mesh.

Figure 3-13 Initial disturbance used to trigger the formation of alternate bars in the simulation of Run P1505 [Lanzoni, 2000].

Figure 3-14 Development and migration of alternate bar trains in the simulation of Run P1505 [Lanzoni, 2000].

Figure 3-15 Growth of bar height with time in the simulation of Run P1505 [Lanzoni, 2000]. The height of bar No.1 almost reaches a steady value.

Figure 3-16 Comparison between the calculated and measured Run P1505’s [Lanzoni, 2000] longitudinal profiles of bar height.

Figure 3-17 Experimental result and computed bed topography of coexistence of forced bar dominated case in Run F2 [Wu and Yeh, 2005].

Figure 3-18 Comparison of measured and computed results of bed topography at four specified sections of a width-variation cycle in the case of coexistence of forced bars dominated case in Run F2 [Wu and Yeh, 2005].

Figure 3-19 Experimental result and computed bed topography of coexistence of free bar dominated case in Run F7 [Wu and Yeh, 2005].

Figure 3-20 Comparison of measured and computed results of bed topography at four specified sections of a width-variation cycle in the case of coexistence of free bars dominated case in Run F7 [Wu and Yeh, 2005].

vii

Figure 4-1 Dimensionless bedforms in B15 series. Bedforms are side bars in most cases of B15 series, except in W02 series are central bars. The dimensions of bedforms are proportional to the amplitude of channels with variable-wdith.

Figure 4-2 A bed disturbance imposed at the upstream extended reach.

Figure 4-3 In B15A01W06 (a) Coexistence of free and forced bars; (b) Free bar components.

Figure 4-4 Definition sketch of bar height profile.

Figure 4-5 Free bar component in (a) A04 which the amplitude of channels is fixed and the wave number is altered and (b) W04 series at 8 hour. W00 and A00 represent the same straight channel run.

Figure 4-6 Bar height (BH), celerity (BC) and wavelength (BL) evolves with time in B15 series.

Figure 4-7 Acquire the equilibrium stage of target bars.

Figure 4-8 Bar height (BH), wavelength (BL) and celerity (BC) reach to the equilibrium state.

Figure 4-9 Ratio of equilibrium (a) bar height, (c) wavelength and (e) celerity vs.

Amp. factor, (b) bar height, (d) wavelength and (f) celerity vs. WN factor.

Figure 4-10 Ratio of equilibrium (a) bar height, (b) wavelength and (c) celerity vs.

forcing factor.

viii

List of Symbols

A    c   dimensionless amplitude of width variation A_e element area









B * channel half-width

B 0 mean half-width

b 0 channel wall perturbation

C f the bed friction coefficient

C s dimensionless local curvature of streamline D sm average grain size of sediment

d s grain size

D w flow depth (m) Fr Froude number









K s the roughness height N i the shape function of node i N_i a diagonal matrix of N _i

ix

ˆN_i the weighting function matrix of node i





^

^

P the pressure

Q b sediment transport discharge









Q n the fluid flux across the boundary segment









q 0 the unit discharge imposed in boundary segment

Rij the vertically averaged Reynolds stress

r an empirical coefficient

S 0 slope of channel S_t the source vector

U * the shear velocity













V h the vertically averaged eddy viscosities in the horizontal direction

x













z b bed elevation

 angle between the sediment trajectory and x-direction

β ratio of mean half-width to reference flow depth

 angle between the local bed shear stress and x-direction

 bed load intensity

 e the segment of the boundary element

 the solution of the governing equation

 the approximate solution of the element

 bed elevation

b dimensionless wave number of width variation

c wavelength of channels with variation width

p bed porosity

 Shields stress

c dimensionless critical shear stress

 density of fluid

s density of sediment

xi









 e an element domain

 an upwind coefficient

i the solution of node i

Chapter 1 Introduction

1.1 Statement of Problem

The evolution processes of bed configuration in the alluvial rivers are abundant and complicated. The bedform of natural rivers may be generated by the forcing effect of the channel geometry, e.g., width variation and curvature, and any disturbances that trigger free deformation on the riverbed. The former may yield the stationary forced bar formation and the latter may induce the free migrating alternate bars, as shown in Figure 1-1. In natural rivers, forced and free bars often coexist and interact with each other. The flowfield and morphology of a river are unsteady and dynamic, and thus affect river management issues such as channel training, flood defense, riverine habitat conservation/restoration, and transport of pollutants. In order to accurately predict the evolution of dynamic bedform in natural channels, a better understanding of the influences of these forcing effects on free migrating alternate bars is necessary.

The aims of this study are to develop a two-dimensional (2D) finite element (FE) morphodynamic model and systematically investigate these problems. The results of this study may provide useful information guiding more dynamic and comprehensive practices of river engineering.

Figure 1-1 (a) Free alternate bars in Shi-hu Creek, Taiwan; (b) Central bars in Tai-ping Creek, Taiwan.

(a)

(b)

1.2 Literature Review

In natural rivers bars are a kind of general large-scale bed topography which the bar height and wavelength are of the order of the flow depth and the channel width, respectively. According to the geometry, kinematic characteristics and emergence location in the channel, bars are grouped into free and forced bars. It is now well established that bars formation can be explained as the results of an instability mechanism which come from the channel geometry nonuniformities or the perturbation on the planimetric river bed configuration. Recent studies of free and forced bars are described as follows.

1.2.1 Study of free bars

Free bars migrating downstream in a straight channel belong to an instability mechanism induced by the spontaneous perturbation on the planimetric erodible bed and are characterized by a sequence of steep consecutive diagonal fronts with deep pools at the lee face and gentler riffles along the stoss face. Figure 1-2 displays a top view of the free bars which are so called alternate bars due to the diagonal front across the transverse direction. Experimental observations support the perturbation on the river bed perturbs flow and trigger a series of bar formation downstream [Fujita and Muramoto, 1985; Garcia and Nino, 1993; Lanzoni, 2001].

  Figure 1-2 Sketch of alternate-bar structure. [Colombini et al., 1987]

In the early stage of studies of free bars, linear instability theory offered a convenient tool to investigate a selection of the most unstable wavelength promoting free bars to develop [Callander, 1969; Engelund and Skovgaard, 1973; Parker, 1976;

Fredsoe, 1978; Blondeaux and Seminara, 1985; Nelson and Smith, 1989; Lanzoni, 2000]. Figure 1-3 shows the representative result of the linear theory, the solid line represent the neutral state at which the growth rate of perturbation are eliminated and distinguish the stable and unstable conditions to judge the developing of free bars.

The shortcoming of linear theory appears when the width-to-depth ratio becomes significantly large because the nonlinear terms become more important. The nonlinear interaction between finite amplitude disturbances of different wavelength of free bars may not be ignored. As an overall trend, linear theory underpredicts bar wavelength and overpredicted the bar celerity.

Figure 1-3 A typical neutral curve for alternate bar formation. [Colombini et al., 1987]

The weakly nonlinear theory developed in the neighborhood of critical conditions has been applied to derive the finite amplitude equation for the marginally unstable bed forms. The Landau equation derived by Colombini et al. (1987) and Fukuoka (1989) have an ability to capture the long-term behavior of a single unstable wave. Schielen et al. (1993) based on the Landau equation to obtain the

Ginzburg-Landau equation that has the ability to capture the evolution of the envelope amplitude of the wave group. The weakly nonlinear theory improves the prediction of equilibrium bar height, but not wavelength and celerity, due to that the weakly nonlinear theory does not account for the modifying of wavelength of the perturbation when the free bars are developing.

The numerical approach applied to the fully nonlinear perturbation had been proposed in the literature [Nelson and Smith, 1989; Colombini and Tubino, 1991;

Defina, 2003; Bernini et al., 2006]. The perturbations are given artificially in time and space trigging bar formation in the numerical simulation. The process of bar evolution in the numerical simulation is similar to the experiment observation [Lanzoni, 2001].

A fully nonlinear numerical model has the ability to describe the nonlinear interaction between free bars in the evolution process in an infinite straight channel. The equilibrium bar height, wavelength and celerity are strictly related to one another regardless of the type of initial perturbation [Difina, 2003]. Figure 1-4 shows the main result of Defina’s numerical experiments, where the bar height and wavelength are inversing proportional to the celerity in the equilibrium.

Figure 1-4 Free bars wavelength and height as a function of bar celerity in Defina’s numerical experiments. [Defina, 2003]

1.2.2 Study of forced bars

The channel curvature and width variation are two types of channel geometry nonuniformities that generate the forced bars. The channel curvature induced point bars had been widely studied in the literature [Ikeda and Nishimura, 1985; Blondeaux and Seminara, 1985; Seminara and Tubino, 1989; Parker and Johannesson, 1989;

Whiting and Dietrich, 1993; Seminara et al., 2001].

The aim of our study is focused on the channels with variable width. According to the transverse bed deposition and scour at the wide section, the forced bar in the channels with variable width are grouped into central and side bars [Bittner, 1994;

Repetto et al., 2002; Wu and Yeh, 2005]. Figure 1-5(a) and Figure 1-6(a) display the experimental bed configurations of central and side bars [Bittner, 1994; Wu and Yeh, 2005]. In the numerical models the correction for the effect of secondary helical flow is necessary to simulate the forced bars. Figure 1-5 and Figure 1-6 display the model without helical flow would not predict the bed forms precisely. Repetto (2002) concluded that the wave number of width variation determines the bar type.

Figure 1-5 Contour plots showing (a) central bars observed in experiment [Wu and Yeh, 2005] and results of (b) 2-D-Cs model [Wu and Yeh, 2005], (c) 2-D model [Bittner, 1994], (d) 2-D model [Repetto et al., 2002]

Figure 1-6 Contour plots showing (a) side bars observed in experiment [Bittner, 1994] and results of (b) 2-D-Cs model [Wu and Yeh, 2005], (c) 2-D model [Bittner, 1994], (d) 2-D model [Repetto et al., 2002]

Wu and Yeh (2005) further concluded that the variable-width induced forced bars are a function of width-to-depth ratio, dimensionless wave number of width variation, dimensionless shear stress of reference uniform flow and dimensionless grain size.

The width-to-depth ratio and dimensionless wave number of width variation mainly determine the bed form development, however, dimensionless shear stress of reference uniform flow and grain size influence the bar height only. Figure 1-7 shows that the predictions of Wu and Yeh (2005) studies are in agreement with the experimental observations [Bittner, 1994; Wu and Yeh, 2005].

Figure 1-7 Comparisons between observed values of beta and corresponding values of BetaC1 and BetaC2. [Wu and Yeh, 2005]

1.2.3 Free bars affected by channel geometry

Kinoshita and Miwa (1974) first observed that the development of alternate bars is suppressed by channel curvature in their experiments. In particular, alternate bars do not develop when the channel sinuosity exceeds a threshold value. Turbino and Seminara (1990) used the perturbation expansion method to interpret this phenomenon theoretically with reference to a regular sequence of small-amplitude meanders. Their theory has the ability to determine the threshold value of channel curvature above which free bars are suppressed as a function of meander wavenumbe for given flow and sediment conditions. Whiting and Dietrich (1993) conducted a series of experiments to investigate the free bars migrating through channel bends and found that free bar migration was constrained by the wavelength of the meander channel. In their observations the migration of free bars were non-uniform and temporarily stalled when in phase with the curvature-induced topography. Figure 1-8 shows the experimental results regarding the unsteady migration of free bars through meanders.

Figure 1-8 Unsteady migration of bars through meanders. Open squares indicate the position of leading edge of bars attached to the left bank; solid squares, the leading edge of bars attached to the right bank.

Lanzoni and Tubino (2001) conducted a series of flume experiments to study the free bars development in channels with variable width. They compared the bar height and bar wavelength of free bars in a straight channel and channels with variable width, and found that both of them are suppressed when encountering the width variation, as shown in Figure 1-9. Tubino et al. (2000) studied analytically the suppression of alternate bars exerted by the channel width variation. A perturbation method with linear stability theory was used to analyze this issue. As a result, the suppression of alternate bars in channels with width variation was characterized by a correction factor of growth rate of free bars. However, it is not able to discuss the bar height, wavelength and celerity influenced by width variation. Numerical experiments are expected to capture the coexistence of free and forced bars and performed in this study to investigate the influences of forcing effect on free bars.

Figure 1-9 Comparison between measured values of (a) the amplitude and (b) dimensionless wavenumber of the leading Fourier component of bed topography associated with alternate bars in constant width and variable width experiments. [Lanzoni and Tubino, 2001]

1.2.4 Morphodynamic model

While 3D morphodynamic models can be used to study in detail the evolution of river morphology, 2D models are more efficient in practical applications. A number of researchers have devoted to the development of 2D morphodynamic models. Finite difference (FD) schemes have been commonly used to investigate the alluvial bend morphology and meandering channels. Koch and Folkstra (1981) applied a simplified 2D model to curved alluvial flumes of constant circular bends. Struiksma (1985) used a 2D FD model to reproduce the observed patterns of scour and deposition along a meandering reach of the Waal River. Struiksma et al. (1985) simulated the sour and deposition measured in the laboratory curved flumes. Shimizu and Itakura (1989) modeled the bed evolution in a sine-generated meandering channel. Kassem and Chaudhry (2002) applied a boundary fitted FD model to simulate some of the alluvial bend experiments of Koch and Folkstra (1981) and Struiksma et al. (1985).

The shortcoming of 2D morphodynamic models is that the momentum transport by secondary currents of 3D flow structures is neglected [Shimizu et al., 1990]. The streamline curvature is employed to reflect the effect of secondary flow for correcting the direction of sediment transport. The relation is expressed as follows [Engelund, 1974; Struiksma et al., 1985]:

tan _s ^w

c

a D

 _{ }^ r ^_

  (1-1)

where _s is the angle between the bed shear stress and depth-average flow direction, D is flow depth, w r is the local radius of curvature of the streamline, and _c a is a friction coefficient ranging between 5 and 12 [Engelund, 1974]. However, equation (1-1) tends to overestimate the effect of secondary flow in the case of strong curvature [Blanckaert and de Vriend, 2003; Blanckaert and Graf, 2004].

Vasquez (2005) incorporated the VAM model rather than the traditional VA

model into a 2D finite element (FE) morphodynamic model. VAM model assumes a vertical distribution of the velocity components across the flow depth, thus is by itself able to describe the secondary flows with the governing equations. Vasquez (2005) carried out numerical simulations on the scour and deposition in curved channels and meandering rivers by using VAM equations in the morphodynamic model. However, his model failed in the channels with variable width and the convection-dominated bed evolution because the streamline upwind Petrov-Galerkin (SUPG) scheme [Hicks and Steffler, 1992; Ghamry, 1999] was applied to the flow model but not to the bed evolution model [Vasquez, 2005]. The FE model with SUPG scheme has the ability to process the convective evolution of bed forms.

Nelson and Smith (1989) simulated the evolution of alternate bars using a standard FD scheme. The model reproduced the generation of free bars downstream an initial disturbance and the simulated results were similar to those observed by Fujita and Muramoto (1985). Defina (2003) used a 2D FE model to reproduce the experimental results of Lanzoni (2000), where migrating alternate bars developed in the straight channel from an initial flat bed. The initial disturbance was used in these numerical experiments to trigger the generation of free bars. Qualitatively speaking, the simulated results are in good agreement with the experimental results [Lanzoni, 2000] and the weakly nonlinear solutions [Colombini et al., 1987; Schielen et al., 1993]. However, different types of disturbance would lead to different bar characteristics. More recently, Bernini et al. (2006) simulated the generation of free bars in a straight, rectangular channel with both supercritical and subcritical uniform flows using the ADI scheme. The results were used to study the effect of gravity due to the transverse bed slope on the equilibrium geometric and kinematic bar characteristics. Although Defina (2003) and Bernini et al. (2006) successfully simulated the generation of free bars in straight, rectangular channels and

demonstrated qualitatively the agreement with experimental data and analytical solutions, none has incorporated the effect of secondary flow and simulated the case of variable-width channel.

1.3 Scope of Study

The aim of this study is to conduct numerical experiments to investigate the influences of the forcing effect on free bars. Wu used a morphodynamic model composed of the hydrodynamic and bed evolution equations, both of them belong to the hyperbolic equation. The streamline upwind Petrov Galerkin (SUPG) scheme is applied to both the hydrodynamic and bed evolution equations to overcome the defect for which the traditional Galerkin scheme may fail in the hyperbolic equation. The details of the numerical models are described in Chapter 2. The vertically averaged (VA) model is used in the hydrodynamic model. The applicability of the VA model in the morphodynamic model is validated and their accuracies are examined in Chapter 3.

In the present study the VA model are chosen to conduct the numerical experiments due to its efficiency.

To investigate the influences of the forcing effect on free bars, a series of numerical experiments are conducted with different amplitudes and wave numbers of channel width variations. The ratios between the characteristics of the free bars developed in the variable-width and straight channels are used to describe the effect of channel width variations on free bars. These characteristics include the evolutions of bar height, wave length and celerity as the train of alternate bars passes through the periodic cycle of width variations. A forcing factor which quantifies the geometry of the channel with variable width is proposed to assess the influence of channel forcing effect on free bars. The numerical experiments and discussion are described in Chapter 4. Finally, overall conclusions are summarized in Chapter 5.

Chapter 2 Mathematical Model

In this chapter, the governing equations of hydrodynamics and bed evolution, along with their closure relations are described first. Then, the two-dimensional (2D) finite element (FE) morphodynamic model and the upwind scheme are presented.

2.1 Governing Equations of Hydrodynamic Model

The vertical average (VA) model is used in present study. It is derived from the fundamental full three-dimensional (3D) Reynolds equations. The full 3D Reynolds equations, including a continuity equation and three momentum equations, are given by

U V W 0

x y z

   

   (2-1a)

1 1 Rxx Rxy Rxz

U U U U P

U V W

t x y z  x  x y z

   

                  (2-1b)

1 1 Ryx Ryy Ryz

V V V V P

U V W

t x y z  y  x y z

  

 

                  (2-1c)

1 1 Rwx Rwy Rwz

W W W W P

U V W

t x y z  z  x y z

   

                  (2-1d)

where







 



density of fluid, and R is the Reynolds stress, defined as the stress in the j-direction _ij acting on a face whose normal is in the i-direction. The VA model is derived from integrating equation (2-1) over the flow depth with a constant velocity, as described below. The shallow water approximation is adopted here, which implies that the

vertical accelerations is negligible compared to the gravity. It is further assumed that the pressure is hydrostatically distributed and flow separation is ignored. [Tubino and Repetto et al, 2000]

The velocity distributions across the whole flow depth are given by

0

0

b w

b

b w

b

z D

z

w

z D

z

w

U Udz D V Vdz

D





 

 



 

 



(2-2a)

where





0 0

b w

b

z D w w w

z

D D D

Wdz U V

t x y

     

    (2-2b)

Integrating (2-1) over the flow depth vertically from the bottom to the flow surface, with the relations (2-2a, b) and an assumption that P is the hydrostatic pressure, leads to the following three equations [Molls and Chaudhry, 1995; Ghamry, 1999]:

0

w x Qy

D Q

t x y

   

   (2-3a)

2 2

2 0

xx x y xy

x x w w w b x

w

w w

Q Q gD D R Q Q D R z

t x D y D gD x



  

   

   

         

        (2-3b)

2 2

2 0

yx yy

y x y w y w w b y

w

w w

Q Q Q D R Q gD D R z

t x D y D gD y



  

 

     

         

        (2-3c)

where











 



2.2 Governing Equation of Bed Evolution Model

The governing equation of bed evolution is the sediment continuity equation, i.e., Exner equation. The processes involved in the Exner equation are bedform translation and diffusion [Lisle et al., 1997, 2001; Cui et al., 2003b], which describe different mechanisms of bed evolution. The Exner equation is given by





where









2.2.1 Closure relations

Closure relations for R^ij,









averaged Reynolds stresses R^ij are approximated with the Boussinesq model [Ghamry, 1999], i.e.,

2

2

xx h x

w

y

yy h

w

y x

xy yx h

w w

R V Q

x D

R V Q

y D

Q Q

R R V

x D y D







    

   



 

    

   



     

       

      



(2-5)

where V is the vertically averaged eddy viscosities in the horizontal directions. For _h simplicity, the case of bed-dominated turbulence is assumed, and the values of the

order of V_h 0.5U D_{* w} is used here [Ghamry, 1999], where

2 2

* 4

x y

U  

 

 

 

       is

the shear velocity,





2 2

0 0 0

2 2

0 0 0

x f

y f

C U U V

C V U V

 

 

  



  



(2-6)

where C is the bed friction coefficient, defined by _f

11 2

2.5log ^w

f

s

C D

K

  

   

 

  (2-7)

where K is the roughness height, defined as 2.5 times the grain size [Cui et al., _s 1996, 2003a]. The sediment transport rates in the horizontal plane are defined as





^

^

where  = angle between the sediment trajectory and x-direction, and is given by

sin sin r z^b

  y



  

 (2-9)

where r is an empirical coefficient reflecting the influence of transverse bed slope, and ranges between 0.3 and 1 [Talmon and Struiksma et al, 1995; Wu and Yeh, 2005];

 = Shields stress; = angle between the local bed shear stress and x-direction, expressed as

0

0 0

sin V _{w s}

U V aD C

  

 (2-10)

where a is the helical flow coefficient; C is the dimensionless local curvature of _s streamline [Wu and Yeh, 2005], defined by

 

 

0 0

2 3/ 2

0 0

1 /

s

V C x U

V U

  

   



 (2-11)

The angle  accounts for the deviation of the zero-average helical flow from the

depth-averaged flow, driven by the bedform or channel variation.

The Meyer-Peter and Muller formula is suitable for evaluating the sediment transport rate of well-sorted grains [Wu and Yeh, 2005], which is given by

 

8  _c

   (2-12)

where  is the bed load intensity, and _c 0.04 is the dimensionless critical shear stress. The sediment transport discharge Q could be evaluated by _b





b sm

Q   gD



   (2-13)

where _s and D are the density and average grain size of sediment, respectively. _sm

2.3 Finite Element Method

2.3.1 Streamline upwind Petrov-Galerkin scheme

The term ‘upwind’ originates from the manner in which the discretization is applied depending on the direction of wave propagation. This wave can be the characteristic waves of the conservation laws or disturbance wave [Giraldo, 1995].

Godunov (1959) introduced the idea that the information from the exact local solution to the Euler equations could be included in the discretization for the purposes of computational studies. He applied this method to the finite volume (FV) method.

Brooks and Hughes (1982) first applied this concept to the finite element (FE) method to solve the convection-dominated flow.

Dendy (1974) and Wahlbin (1974) derived the dissipative Gelerkin scheme to solve the first-order hyperbolic equations. Katopodes (1984) applied the dissipative Galerkin scheme to the non-conservation form of de St. Venant equations and simulated the 1D and 2D hydraulic jumps. The dissipative Galerkin scheme, however,

only considered the progressive characteristic velocity of the hydrodynamic equations in the upwinding term. The characteristic velocities of the hydrodynamic equations, in fact, include the progressive and reprogressive directions. Hughes and Mallet (1986) examined the application of the Petrov-Galerkin method to the symmetric systems of hyperbolic equations. Based on their work, Hicks and Steffler (1992) developed the characteristic dissipative Galerkin (CDG) scheme which used each of the characteristic velocities in determination of the upwinding matrix. Both of the above studies examined their models using the 1D hydrodynamic equations. Recently, Ghamry (1999) succeeded in applying the CDG scheme to solve the 2D vertically averaged and moment (VAM) equations.

The FE method employed in this study is the streamline upwind Petrov-Galerkin scheme [Brooks and Hughes, 1982]. The phrase ‘streamline upwind’ implies that the direction of advection is incorporated into the discretization, thus it is more suitable for hyperbolic equations than the traditional Bubnov-Galerkin scheme [Brooks and Hughes, 1982; Giraldo, 1995]. The advantage of using the FE method is that boundary conditions can be easily imposed on the discretized domain of natural environments with complicated geometry. Throughout this study, triangular elements are used for discretization of the computational domain.

The governing equations of the VA model has a general form that can be expressed as

 

   

x

t x y t

 

 ^{ } 

    

  

f f

S (2-14)

where  is the solution of the governing equations; f_x

 

 

 

 





triangular elements, then the approximate solution of an element is given by

 ³

1

e e e

i i i

 N 

  (2-15)

where N is the shape function of node i, the Lagrange interpolation function is _i^e used in this study; _i is the solution of node i. Substituting (2-15) into (2-14), multiplying the resulting equation with a specified weighting function and integrating it over an element domain  leads to the general form of the equation for _e streamline upwind Petrov-Galerkin scheme [Hicks and Steffler, 1992]:



 

 

 

ˆ 0

e

e e e

y e

e x

i t d e

t x y

 

 



   

      

     

f f

N S for i1, 2,3 (2-16a)

where N^ˆ_i^e is the weighting function matrix of node i in an element, which is defined as :

ˆ 0 0 00

0 0

e

e e

e e e i e e i e i

i i i i x y

e i

N

N x y

x y

N

 

     

 

         

N N

N N W W W (2-16b)

where N is a diagonal matrix of _i^e N ; _i^e  is an upwind coefficient ranging from 0.25 to 0.75, in this study a value of 0.5 is used (Note that (2-16) would reduce to the traditional Bubnov-Galerkin scheme if 0); the element sizes

2 Ae

x y

  ,

with A_e element area [Ghamry, 1999]; W and _x^e W_y^e are the upwind matrices in the x- and y-direction, which characterize the advection mechanism and are the key terms in (2-16). To specify W and _x^e W_y^e, a characteristic dissipative Galerkin (CDG) scheme is adopted here [Hicks and Steffler, 1992; Ghamry, 1999], which gives

 

 

2 2 1

2 2 1

e

x x x y

e

y y x y





  



  



W A A A

W A A A

(2-17)

where A and _x A are the advection matrices in the x- and y-direction, respectively, _y which are derived from (2-14) and expressed as follows:

 

x y t

t x y

   

      

 A  A  S (2-18)

The details regarding the applications of the CDG scheme to the flow dynamics and bed evolution models are described below. The general forms of A and _x A are _y defined as

     

     

     

c c c

mx mx mx

my my my

x x x

w x y

x x x

x w x y

x x x

w x y

f f f

D Q Q

f f f

D Q Q

f f f

D Q Q

  

  

  

    

    

 

 

  

 

 

    

 

  

 

    

 

 

A (2-19a)

     

     

     

c c c

mx mx mx

my my my

y y y

w x y

y y y

y w x y

y y y

w x y

f f f

D Q Q

f f f

D Q Q

f f f

D Q Q

  

  

  

    

    

 

 

  

 

 

    

 

  

 

    

 

 

A (2-19b)

where the subscripts c, mx, my refer to the continuity equation, momentum equations in the x- and y-direction, respectively. The specific forms of A and _x A can be _y derived from (2-3), and are given by

2 0 0

0 0 0 0

0 1 0 2 0

x w

gD U U

U V V U

 

 

  

  

  

 

A (2-20a)

0 0 0 0

2

0 0

0 0 1

0 2

y

w

U V V U

gD V V

 

 

  

  

  

 

A (2-20b)

The eigenvalues of (2-20) are the characteristic celerity of the system of equations (2-3) [Hicks and Steffler, 1992], the reason why it is called ‘characteristic dissipative Gelerkin’ scheme. The eigenvalues of the advection matrices of the VA model, U , ₀

0 w

U  gD , and U₀ gD_w , are the characteristic celerity of (2-3). To calculate the inverse of A_x²A_y² in (2-17), the numerical procedure based on the Cayley-Hamilton theorem is employed [Hoger and Carlson, 1984].

2.3.2 Applying CDG scheme to bed evolution model

While the CDG scheme has been successfully applied to flow dynamics models, to date applying the CDG scheme to solve the sediment continuity equation has never been carried out due to the complexities involved in the formulation [Vasquez, 2005].

In the river morphology model of Vasques (2005) only the hydrodynamic model applied the CDG scheme and the bed evolution model was treated with a hybrid method. As a results, the point bars in the bend were successfully simulated but the free bars and the force bars in the channels with variable width were not reproduced.

It is demonstrated that applying the CDG scheme to the bed evolution may overcome such defects. Because the morphodynamic model is unable to capture the advection mechanism of bed evolution without incorporating an upwind scheme, some simplifications must be made in the formulation. As such, equations (2-3) and (2-4) are modified as

0 0 0

U V

x y

 

 

  (2-21a)

 

0 0 0

0 0 ^{w b x} 0

w

D z

U U U

U V g

t x y x D





 

       

    (2-21b)

 

0 0 0

0 0 ^{w b y} 0

w

D z

V V V

U V g

t x y y D





 

       

    (2-21c)





t x y

 ^ ^  ^

   

   (2-21d)

Since Q is a function of bed shear stress _b  , the gradients of Q in the x- and _b y-direction may be rewritten as

0 0

0 0

2

b b b

f

Q Q Q U V

C U V

x x x x

 

 

          (2-22a)

0 0

0 0

2

b b b

f

Q Q Q U V

C U V

y y y y

 

 

          (2-22b)

Multiplying equations (2-21b) and (2-21c) by U and ₀ V , respectively, and then ₀ summing up the resulting equations leads to

0 0 0 0

0 0 U 0 V 0 0 U 0 V 0 z^b 0 z^b _t' 0

U U V V U V gU gV S

x x y y x y

     

        

         

    (2-23)

Substituting equation (2-23) into equation (2-22) gives

0 0

1

0 0

2

b b b b

f t

Q Q U z V z

gC S

x  U x U y



         (2-24a)

0 0

2

0 0

2

b b b b

f t

Q Q U z V z

gC S

y  V x V y



         (2-24b)

where S and _t1 S are the source terms, whose values do not affect the final results. _t2 Substituting equation (2-24) into equation (2-21d) gives

0 0

0 0

cos sin

cos Q^b, sin Q^b 2 _f z^b, z^b

gC U V

x y U V x y

 

 ^  ^  ^{ }

     

        

     (2-25)

Equation (2-25) implies that flow dynamics is the main driving force for the advection of bed deformations, while the contributions of bed slope and gravity are categorized

as the diffusion mechanism. Substituting equations (2-24) and (2-25) into (2-21d) would result in





t x y

 ^{    }

         (2-26)

where C' is a constant; S is the source term. Equation (2-25) is similar to (2-18), _t3 however, the upwind matrices become two upwind components, which are given by





^{ }

0 0 0 0

, , cos ,sin

b b

z x z y

U V

W W

U V U V

 

 

   

(2-27)

2.4 Model implementation

To implement the morphodynamic model, the flow is assumed to be quasi-steady.

The time derivative terms in the governing equations of the VA model can be neglected, implying that the flow adapts to the altering bed topography immediately.

Such a morphodynamic model is essentially a decoupled one. The flow, sediment, and bed topography computations are executed iteratively with the following procedure:

Step 1: Solve the flow dynamics model with the given bed topography.

Step 2: Evaluate the sediment transport rates with the calculated bed shear stress.

Step 3: Compute the bed topography at the next time using the bed evolution model.

Step 4: Go back to step 1 with the updated bed topography.

Implementations of the flow dynamics and bed evolution models are described below.

Implementing the flow dynamics model with the CDG scheme is described here.

The weak form of the finite element equation and boundary conditions are introduced first. The Newton-Raphson algorithm is then adopted to solve the resulting system of equations.

Equation (2-3) may be rewritten in the form of (2-16a) as the following: