Mathematical modelling of morphological processes in the case of suspended sediment transport

in the case of

Suspended Sediment Transport

Leo C. van Rijn

( \

^

TR diss

1553

in the case of

Suspended Sediment Transport

Proefschrift

Ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, prof. dr. J.M. Dirken, In het openbaar te verdedigen ten overstaan van een commissie aangewezen door het College van Dekanen op dinsdag 23 juni 1987 te 16.00 uur.

door

Leonardus Cornells van Rijn

geboren te Amsterdam

civiel Ingenieur

1987

Grafische verzorging Waterloopkundig Laboratorium Delft

TR diss

1553

Dit proefschrift is goedgekeurd door de promotor prof. dr. ir. M. de Vries

behorende bij het proefschrift:

Mathematical Modelling of Morphological Processes

in the case of Suspended Sediment Transport

L. C. van Rijn

Juni 1987

SAMENVATTING

ABSTRACT

ACKNOWLEDGEMENTS

page

1. INTRODUCTION 1

2. MASS-BALANCE AND MOMENTUM-BALANCE EQUATIONS FOR FLUID-SEDIMENT

FLOWS 5 2.1 Mass-balance equations 5

2.2 Momentum balance equations 8

3. EQUILIBRIUM SEDIMENT TRANSPORT 11

3.1 Introduction . 11 3.2 Sediment transport modes 11 3.3 Bed load transport 12 3.4 Suspended load transport 14 3.5 Bed forms and effective bed roughness 16

4. TWO-DIMENSIONAL MATHEMATICAL MODELLING OF SUSPENDED SEDIMENT

TRANSPORT 18 4.1 Introduction 18 4.2 Simplified mass-balance equation for suspended sediment

transport 18 4.3 Flow velocity field 20

4.3.1 K-EPSILON model 21 4.3.2 PROFILE-model 23 4.3.2.1 Longitudinal flow velocity 23

4.3.2.2 Vertical flow velocity 27 4.3.2.3 Bed-shear velocity 27 4.3.2.4 Computation procedure 28 4.3.2.5 Calibration of PROFILE model 28 4.3.2.6 Verification of PROFILE model 31 4.4 Fluid and sediment mixing coefficients 32

4.4.1 Introduction 32 4.4.2 K-EPSILON model 33 4.4.3 PROFILE model 33

CONTENTS (continued)

4.4.3.1 Vertical distribution of fluid mixing coefficient 33 4.4.3.2 Longitudinal distribution of fluid mixing coefficient 35

4.4.3.3 Calibration of coefficients 35 4.5 Boundary conditions for sediment concentration model 38

4.5.1 Computation domain 38 4.5.2 Inflow boundary 39 4.5.3 Outflow boundary 39 4.5.4 Water surface boundary 39 4.5.5 Sediment bed boundary 39 4.6 Bed level changes 40 4.6.1 Mass-balance equation 40 4.6.2 Suspended load transport 41 4.6.3 Bed-load transport 41 4.7 Numerical solution method and accuracy for sediment

concentrations 42 4.7.1 Solution method 42 4.7.2 Numerical accuracy 43 4.8 Numerical solution method and accuracy for bed level changes... 48

4.8.1 Solution method 48 4.8.2 Numerical accuracy 49 4.9 Sensitivity analysis: Influence of errors in boundary

conditions and basic physical parameters 50

4.9.1 Introduction 50 4.9.2 Influence of errors in boundary conditions 51

4.9.3 Influence of errors in basic physical parameters.... 55 4.9.4 Combined effect of errors in basic physical parameters 64

4.9.5 Conclusions of sensitivity analysis 64 4.10 Verification of SUTRENCH-2D model 64

4.10.1 Introduction 64 4.10.2 Sediment concentrations and transport in a horizontally

uniform flow over a perforated bottom in a flume 65

4.10.3 Migration of a trench in a flume 66 4.10.4 Sedimentation in a trench In the Western Scheldt, The

Netherlands 73 4.11 Special applications 77

4.11.1 Adjustment length of suspended sediment transport 77 4.11.2 Trapping efficiency of dredged trenches and channels 80

5. THREE-DIMENSIONAL MATHEMATICAL MODELLING OF SUSPENDED

SEDIMENT TRANSPORT 86

5.1 Introduction 86 5.2 Basic equations and simplifications 86

5.2.1 Basic equations for suspended sediment transport 86

5.2.2 Basic equations for fluid flow 87 5.3 Fluid and sediment mixing coefficients 91

5.3.1 General 91 5.3.2 Horizontal mixing 92

5.3.3 Vertical mixing 94 5.4 Boundary conditions for sediment concentration model 95

5.4.1 Computation domain 95 5.4.2 Inflow boundary 96 5.4.3 Outflow boundary 96 5.4.4 Bank boundary 96 5.4.5 Water surface boundary 96

5.4.6 Sediment bed boundary 97 5.4.7 Equilibrium concentration profiles and transport 99

5.5 Bed level changes 100 5.5.1 Mass-balance equations 100 5.5.2 Suspended load transport 101 5.5.3 Bed load transport 101 5.6 Numerical solution method and accuracy for sediment

concentrations 101 5.6.1 Solution method 101 5.6.2 Numerical accuracy 102 5.7 Numerical solution method for bed level changes 105

5.8 Test case; partially closed channel 106

5.8.1 Introduction '• 106

5.8.2 Hydraulic conditions 106 5.8.3 Boundary conditions and input parameters 106

5.8.4 Computation results 110

6. SUMMARY, CONCLUSIONS, DISCUSSION AND RECOMMENDATIONS 123

6.1 Summary 123 6.2 Conclusions 123 6.3 Discussion and recommendations 125

CONTENTS (continued)

APPENDIX A EQUILIBRIUM SEDIMENT TRANSPORT

1. Bed load transport

2. Suspended load transport

3. Bed forms and alluvial roughness

APPENDIX B BED-BOUNDARY CONDITIONS FOR SUSPENDED SEDIMENT CONCENTRATIONS

1. Introduction B.l

2. Schematization of sediment transport process B.l

3. Sediment exchange between bed-load and suspended load B.3

3.1 General description B.3 3.2 Functions for bed concentration and upward diffusive

sediment flux _ B.4

3.2.1 Stochastic function for bed concentration B.5 3.2.2 Stochastic function for upward sediment flux B.ll

4. Location of boundary level B.12

4.1 Flat bed B.12

4.2 Bed forms B.13

REFERENCES APPENDIX B

LIST OF SYMBOLS

Wiskundige Modellering van Morfologische Processen In het geval van gesuspendeerd Sedimenttransport

Deze studie geeft een beschrijving en toepassing van een twee-dimensionaal vertikaal en een drie-dimensionaal wiskundig model, waarmee morfologische pro­ cessen in het geval van suspensietransport kunnen worden gesimuleerd. De basisvergelijking van beide modellen is de massa-balans vergelijking van de sedimentdeeltjes in suspensie, waarmee het horizontale en het vertikale trans­ port door konvektie, diffusie en de zwaartekracht wordt beschreven. Uit deze vergelijking kunnen op numerieke wijze de sedimentkoncentraties worden bere­ kend, als de stroomsnelheden, de diffuslekoëfficlënten en de valsnelheid van het sediment bekend zijn, en er tevens geschikte randvoorwaarden zijn opge­ legd.

In het twee-dimensionale morfologische model (SUTRENCH-2D) worden de stroom­ snelheden en diffusiekoëfficlënten berekend door middel van een PROFIELEN-mo-del, dat is gebaseerd op de toepassing van flexibele profielen in vertikale richting. Het snelheidsprofiel bestaat uit een logarithmisch profiel om de stroming nabij de bodem weer te geven en een machtsprofiel om de invloed van drukgradiënten in de stroming ten gevolge van steile bodemhellingen weer te geven. De verdeling van de vertikale diffusie-koëfficiënt over de diepte wordt beschreven door een parabolisch profiel in de onderste helft en een konstante waarde in de bovenste helft van de vertikaal. Een nadeel van de PROFIELEN-methode is de noodzaak tot ijking van een aantal koëfficlënten met behulp van meetgegevens of berekeningsresultaten van meer geavanceerde wiskundige model­ len. Bij de ijking van het stroomsnelheidsmodel is gebruik gemaakt van goot-proeven betreffende snelheidsmetingen in geulen van verschillende afmetingen en dwars op de stroom gesitueerd. De diffuslekoëfficiëntverdelingen zijn ge­ ijkt met behulp van berekende waarden afkomstig van het twee-dimensionale ver­ tikale K-EPSILON model, het meest geavanceerde wiskundige model dat thans be­ schikbaar is.

In het drie-dimensionale morfologische model (SUTRENCH-3D) worden de stroom­ snelheden berekend met behulp van een diepte-gemiddeld stromingsmodel in kom-binatie met logarithmische snelheidsprofielen om een quasi-driedimensionaal snelheidsveld te verkrijgen. De diffuslekoëfficiëntverdeling in dit model be­ staat eveneens uit een parabolisch-konstante verdeling in vertikale richting. De bodemrandvoorwaarde wordt in beide modellen weergegeven door het opleggen van een evenwichtskoncentratie of een vertikale evenwlchtssediment flux aan de bodem als funktie van de lokale bodemschuifspanning. De randvoorwaarde ter plaatse van de instroom, uitstroom, en oeverranden bestaat uit het opleggen van het lokale evenwichtskoncentratieprofiel. Met behulp van het 2D en 3D mo­ del is er uitgebreid onderzoek uitgevoerd naar de invloed van de belangrijkste fysische parameters en randvoorwaarden op de morfologische processen (gevoe­ ligheidsanalyse) . De resultaten van het 2D model zijn geverifieerd met behulp van gootproeven en veldgegevens betreffende het sediment transport over een poreuze bodem in een goot en de aanzanding van geulen dwars op de stroom. De resultaten van het 3D model zijn nog niet experimenteel geverifieerd, maar de berekeningen voor een proefgeval betreffende de morfologie in een gedeeltelijk afgesloten waterloop lijken realistische resultaten te tonen.

ABSTRACT

Mathematical Modelling of Morphological Processes In the Case of Suspended Sediment Transport

A two-dimensional vertical and a three-dimensional mathematical model have been developed to study the morphological processes in case of suspended sedi­ ment transport. The basic equation of these models is the convection-diffusion equation for the suspended sediment particles. To solve this equation numeri­ cally, the flow velocities, the sediment mixing coefficients and the particle fall velocity must be known, while also appropriate boundary conditions must be prescribed.

In the two-dimensional morphological model (SUTRENCH-2D) the flow velocities and sediment mixing coefficients are computed by a PROFILE model, which is based on the application of flexible profiles in vertical direction. Flume data have been used to calibrate the PROFILE model for the velocities. Compu­ tation results of the sophisticated K-EPSILON model have been used to cali­ brate the PROFILE model for the sediment mixing coefficients.

In the three-dimensional morphological model (SUTRENCH-3D) the flow velocities are computed by a two-dimensional horizontal depth-averaged flow model in com­ bination with logarithmic velocity profiles. The bed boundary condition in the 2D and 3D model are represented by prescribing the equilibrium bed concentra­ tion or the equilibrium upward sediment flux at the bed as a function of hy­ draulic variables. Equilibrium concentration profiles are prescribed at the inflow, outflow, and bank boundaries.

Applying the 2D and 3D model, the influence of the basic physical parameters and the boundary conditions on the morphological processes has been Investi­ gated in detail.

The SUTRENCH-2D model has been verified extensively using flume and field data. The SUTRENCH-3D model has not yet been verified, but computations for a

test case seem to show realistic results.

ACKNOWLEDGEMENTS

The writer gratefully acknowledges the financial support for this study by the basic research programme of Delft Hydraulics and by the applied research pro­ gramme (TOW) of the "Rijkswaterstaat" of the Delft Ministry of Transport and Public Works.

Dr. C.B. Vreugdenhil and Dr. K. Meijer of Delft Hydraulics are acknowledged for the numerical modelling work.

Messrs H. van Zanten and C. ten Napel are acknowledged for the programming work.

A fundamental phenomenon of non-equillbrlum sediment transport Is the conti­ nuous adjustment of the sediment transport to the sediment transport capacity, especially in environments with a predominant suspended load transport. The morphological processes in such conditions can be simulated sufficiently accurately by one-dimensional or two-dimensional horizontal morphological mo­ dels based on an equilibrium sediment transport formula when the horizontal grid size of the model is much larger than the adjustment length scale of the sediment transport process. In that case the local sediment transport can be assumed to adjust almost instantaneously to the local flow velocity. General­ ly, such an approach is applied for short-term and long-term (decades) predic­ tions at large scales (10 to 1000 km) in rivers (De Vries, 1965, 1971), estu­ aries and seas (Boer et al, 1984). When the horizontal grid size is much smal­ ler than the adjustment length of the sediment transport process, the adjust­ ment process itself should be modeled. Various approaches are possible: (1) depth-integrated models, (2) two-dimensional vertical models and (3) three-di­ mensional models.

In this study the main objective is the two-dimensional vertical and three-di­ mensional mathematical modelling of suspended sediment transport (particles larger than about 50 p.ra) in non-stratified free-surface flows, while also some basic phenomena of the equilibrium sediment transport process such as the par­ ticle movement near the bed, the bed concentrations and the vertical transfer of fluid momentum and sediment mass, are considered.

Equilibrium sediment transport

The transport of sediment particles by a flow of water can be in the form of bed load and suspended load. Although in natural conditions there is no sharp division between the bed load and the suspended load layer, it is necessary to define two distinct layers for mathematical representation. Bed load transport can be defined as the transport of particles by rolling, sliding and saltating processes. The trajectories and velocities of bed load particles were studied in detail by Kallnske (1947), Einstein (1950) and Bagnold (1956). Later the writer (Van Rijn, 1984a) studied the motion of bed load particles by solving

the equation of motions for a saltating particle.

A general equation for equilibrium suspended sediment concentration profiles was introduced by Rouse (1937). This equation can be used to compute the equi­ librium sediment concentrations as a function of height above the bed when the particle fall velocity, the sediment mixing coefficient and the reference

con-

-2-centratlon are known. Einstein (1950), Engelund and Fredsde (1976) and later the writer (Van Rijn, 1984b) proposed expressions for the reference concentra­ tion as a function of sediment properties and flow conditions. The experiments of Coleman (1970) provided a better understanding of the vertical distribution of the sediment mixing coefficient. The experimental data of Coleman were used by the writer to study the relationship between the transfer of sediment mass and fluid momentum. The writer also investigated the influence of the sediment particles on the vertical transfer of turbulence energy (damping effects). The work of the writer related to the equilibrium sediment transport is summarized in Chapter 3; a detailed description is presented in Appendices Al, A2 and A3.

Two-dimensional vertical mathematical modelling of suspended sediment trans­ port

Early attempts of two-dimensional vertical mathematical modelling for non-uni­ form conditions have been presented by Kalinske (1940) and Dobbins (1943). Later a more general mathematical approach was presented by O'Connor (1971). Smith and O'Connor (1977) presented a two-dimensional vertical model based on the laterally-integrated momentum and continuity equations for the fluid and sediment phases, and an "one-equation" turbulence closure model to represent the fluid shear stresses and diffusion coefficients. Celik and Rodl (1984) presented a similar model as that of Smith and O'Connor. However, the former applied a "two-equation" turbulence closure model (K-Epsilon model). Bechteler and Schrimpf (1984) presented a relatively simple two-dimensional model ne­ glecting vertical convection and horizontal diffusion. The fluid field was represented by logarithmic velocity profiles. In the Netherlands two-dimen­ sional modelling of suspended sediment was initiated by Kerssens (1974). Kerssens et al (1977 and 1979) presented a two-dimensional vertical model for gradually varying flows neglecting vertical convection and horizontal dif­ fusion. Logarithmic velocity profiles were used to represent the fluid velo­ city field. The vertical sediment mixing coefficients were represented by a parabolic-constant distribution. A concentration-type boundary condition was applied at the bed, assuming instantaneous adjustment to equilibrium condi­ tions close to the bed. A six-point implicit finite-difference method was used to solve the basic convection-diffusion equation. To obtain a better description of the velocity profiles in horizontally non-uniform flow, experi­ ments concerning the flow field in trenches (perpendicular to the flow direc­ tion) were performed. This resulted in a semi-empirical representation of the velocity profiles in the two-dimensional vertical mathematical model (Van Rijn, 1980c).

From 1981 onwards the mathematical model was further improved by applying a finite-element solution method. The flow width was introduced to extend the applicability of the model to gradually varying flows in transverse direction. A (PROFILE-) method based on the application of flexible profiles (shape

functions) was Introduced to obtain a better description of the velocity profiles. This method was also applied to describe the spatial variation of the sediment mixing coefficients. New expressions for the bed boundary con­ dition (bed concentration or upward sediment flux) were proposed and implemen­ ted in the model. The developments of the two-dimensional vertical model (SUTRENCH-2D) since 1981 are presented in Chapter 4 (see also Van Rijn, 1986b).

Three-dimensional mathematical modelling of suspended sediment transport

Until recently, most applications of suspended sediment models have been restricted to transport processes in the two-dimensional vertical plane, mainly for economic reasons (computer cost). For computations at wide horizon­

tal scales (estuaries, coastal zones, seas), usually, another approach has been followed applying depth-averaged flow models in combination with a simple equilibrium sediment transport formula (Boer et al, 1984; Coeffë and Péchon, 1982; Lepetlt and Hauguel, 1978). Such an approach is only valid when the horizontal length scale of the adjustment process of the local suspended sediment transport to the local equilibrium transport is smaller than the maximum allowable horizontal grid size of the applied mathematical model. If this latter requirement is not satisfied, it is of essential importance to represent the horizontal and vertical adjustment of the sediment transport process in the model. Basically, two types of modelling can be used: (1) the depth-integrated approach as introduced by Galappatti and Vreugdenhil (1985) and (2) the three-dimensional approach as applied by Sheng and Butler (1982), Miller (1983), Wang and Adeff (1986) and others. The application of the depth-integrated approach is limited to situations where the difference between the local true suspended sediment transport and the local equilibrium transport is relatively small (Wang and Ribberink, 1986), otherwise the results are not accurate. An advantage of the depth-integrated approach is the relatively low computer cost compared with that of three-dimensional models which are be­ coming increasingly popular because of impressive advancement of computer technology (memory size and speed). A good example of these latter develop­ ments is the full unsteady three-dimensional model for fluid velocities and sediment concentrations of Wang and Adeff (1986). Relationships proposed by the writer (Van Rijn, 1984b) were applied by them to model the vertical dis­ tribution of the sediment mixing coefficient and the sediment entrainment at the bed.

-4-In Chapter 5 the development of a three-dimensional mathematical model for suspended sediment transport in gradually varying non-stratified free-surface flows is presented.

As a first approach the velocity field has been represented as simple as pos­ sible by applying a depth-averaged flow model in combination with logarithmic velocity profiles. This restricts the applicability of the model to (thin) shear layer flows over a very gradually varying bottom. At a later stage fully three-dimensional flow will be studied.

2. MASS-BALANCE AND MOMENTUM-BALANCE EQUATIONS FOR FLUID-SEDIMENT FLOWS

2.1 Mass-balance equations

Applied to Instantaneous variables, the mass-balance equations of a unit volume are given by:

fluid : |

o

- {p(l-c)} + 1 ^ - {p(l-c) u

} = 0 (2.1)

In which:

c = local sediment volume concentration (-)

Uf = local fluid velocity (m/s)

us = local sediment velocity (m/s)

p = fluid density (kg/m3)

p = sediment density (kg/m3)

x = spatial coordinate (m)

t = time (s)

The effect of turbulence can be introduced by applying the Reynold's-procedure in which the variables are represented as the sum of a time-averaged (overbar) component and a fluctuation (prime) component, as follows:

u = Ü + u' (2.3)

c = c + c' (2.4)

Substitution of Eqs. (2.3) and (2.4) in Eqs. (2.1) and (2.2), averaging over time and assuming the fluid and sediment density to be constant, results in:

fluid : ^- (1-i

'ka-o ^ N v - V

1 =0 (2

-

5)

sediment: ^ (c) + — (c u ^ + c ' u ^ ) = 0 (2.6)

The sediment velocity is assumed to be equal to the fluid velocity with excep­ tion of the vertical direction where a constant slip velocity equal to par­ ticle fall velocity (wg) ±s assumed. Thus:

u = üf . - w 6. (6 = 6 = o, 6 = 1) <2-7)

-6-The eddy-viscosity concept is applied to represent the turbulence-induced transport components. In this study the coefficients expressing the transfer of fluid momentum and sediment mass are shortly called fluid and sediment mixing coefficients.

Applying the eddy-viscosity concept, the turbulence-induced components are:

PTTT - -c. SS_ f,1 f ox

s.1 s ox i in which:

E = fluid mixing coefficient e = sediment mixing coefficient

Substitution of Eqs. (2.7), (2.8) and (2.9) in Eqs. (2.5) and (2.6) results in: fluid : |- (1-c) + I — {(1-c) Zc , + ec | M = 0 (2.10) ( 2 . 8 ) ( 2 . 9 ) ( m2/ s ) ( m2/ s ) : ^ ( l - c )+f ^ { ( l - O uf > 1 + £ f^ - } = 0 sedlment:

IT

' )

h; t« <«

,i -

s

i> -

s ! y ■ ° <

-

In the longitudinal-vertical (x-z) plane the mass-balance equations read, as follows:

Ü H Ü :§ r ( l - c ) + f e { ( l - c ) u + ef( g ) }+|r{ ( l - c ) w + E f ff} = 0 (2.12)

s e d l m e n t : It ( O + k Cc Ü " ts |f) + |j- {C ( W ^S) - .e, |f} = 0 (2.13) In which:

z = vertical coordinate (m) u = fluid velocity in x-direction (m/s)

w = fluid velocity in z-direction (m/s)

For a steady and uniform flow (öu/ot = 0, öc/öt = 0, ou/8x = 0, öc/ox = 0) i the equations reduce to:

fluid : (1-c) i + E f ^ « 0 (2.14)

sediment: ê (w-w ) - e ~r- = 0 (2-15)

Elimination of the vertical fluid velocity (w) yields:

( 1_ c ) c ws + {E s + c (Ef - es» Jf- = 0 (2.16)

Equation (2.16), first presented by Hallbron (1949) and Hunt (1954), repre­ sents the sediment concentration profile for equilibrium conditions. By the introduction of the mass-balance equation for the fluid the vertical return flow due to the fluid displaced by the falling particles has been taken into account. Assuming that the fluid and sediment mixing coefficients are approxi­ mately equal for fine sediment particles (e = e ) , it follows that:

(1-S) c w + e ^- = 0 (2.17) s s dz

Equation (2.17) can also be expressed, as:

c w ,, + e ;r- = 0 (2.18) s.eff s dz

in which:

w ,, = (l-c) w = effective particle fall velocity (m/s)

Experimental research by Richardson and Zaki (1954) has shown that the fall velocity is not only affected by the return flow due to the displaced fluid but also by additional effects such as: particle collisions, particle-induced turbulence and modified drag coefficients. The overall effect can be represen­ ted by:

w

s,eff = d - V

s <

'

in which:

ws = particle fall velocity in a clear, still fluid (m/s)

a = coefficient (= 4 to 5 for particles in the range of 50 to 500 ^m) (-)

The influence of the sediment particles on the turbulence characteristics resulting in a damping of the turbulence and hence a reduction of the effec­ tive mixing coefficient should also be represented.

For small concentrations (l-c = 1) Equation (2.18) reduces to the following expression, first presented by Rouse (1937):

c w + e ^ - = 0 (2-20) s s dz

2.2 Momentum-balance equations

To represent the modification of the fluid mixing coefficients and hence the fluid velocity profiles due to the presence of the suspended sediment par­ ticles, the momentum-balance for the fluid-sediment mixture is considered. In the present analysis it is assumed that the fluid and sediment velocities are equal in the horizontal directions, while the usual assumption of a constant slip velocity equal to the particle fall velocity is applied for the vertical direction (Rouse 1937, Hunt 1954).

Applied to instantaneous variables, the momentum-balance for the fluid-sedi­ ment mixture is given by:

3_ at

0 0 0 V

(p u .) + » — (p u . u .) = - ~z.— (p ) + » — (t ,,) + p s 6, (2.21)

vl^m n,i' öx v hm m,i m,y dx VHm'' dx v m,ij' ^m 6 1 v ' in which:

p = p(l-c) + p c = density of fluid-sediment mixture ' (kg/m3) m s

u . = u , . - w 6. = local velocity of fluid-sediment mixture (m/s) m,i f,l s i = local pressure (N/nr) Pm „.v

X i j

g

= local viscous shear stress (N/m2)

= acceleration of gravity (m/s2)

The viscous shear stress due to the shearing of the intergranular fluid and the shearing induced by the sediment particle interactions is represented by:

du

(2.22)

= dynamic viscosity coefficient of the fluid as

modified by the presence of the particles (kg/sm) = dynamic viscosity coefficient of the fluid phase (kg/sm)

= coefficients (-)

Usually, the viscous shear stress is neglected in the momentum-balance equation for a clear flow. For a two-phase flow over a movable bed, however, the sediment concentrations may be relatively large near the bed. Consequent­ ly, the viscous shear stresses near the bed may become important because of the increased viscosity.

The effects of turbulence are introduced by applying the Reynolds'-procedure, as follows: V lm i n *m

^

u m < " u i + u'm i <2'2 3> m,i m,i m,i

P = P + P' (2-24>

'm 'm r m

p = p + p' (2.25) Km 'm H m

Substitution of Eqs. (2.22), (2.23), (2.24) and (2.25) in Eq. (2.21) and averaging over time, yields:

ïït (Pmum,i) + öT" (Pmum,iun,,j> = ~ ox- <Pm> + dx~ ^ . i j + \ i ] ' + Pm 8 6i (2.26) in which:

-v

t , . = H O u ,/ox.) = viscous shear stress (N/nr) m,ij m m,i i

-t

T = - p U1 U1 — U pfUf +

m,ij m m,i m,j m,i m m,j

- u p' u' = effective turbulence-related m, j m m,i

shear stress ( i * j ) (N/m2)

pm = p (1-c) + ps c = mean density of mixture kg/m3)

M-m = ^o (!■ + <*lc) 2 = dynamic viscosity of mixture (kg/sm)

Pm = (Ps~P) c' = density fluctuation of mixture (kg/m3)

To derive Equation (2.26), the following terms have tentatively been neglec­ ted:

time dependent fluctuations, &(p'u' , ) / ö t m m,i

fluctuation terms due to viscous shear stress

tripple fluctuation terms of the turbulent shear s t r e s s , p ' u ' v P ^ m m.i m,j turbulence-induced pressure, p u ' u

r m,i m,i

Further study of Eq. (2.26) applying scale analysis is necessary to determine the relative importance of the various components of Eq. (2.26) including those which have tentatively been neglected.

The turbulence-induced shear stress (t ) is represented by (i # j):

- - ö um l

. . - (p - pc + p c) e, .. ' - (p - p) e u ., a — (2-27) m,ij vf^ ^ 's f 5x V Ms v' s m,i ox

-t , - , u m,i , 5c

j j In the x-z plane (with velocities u and w) Equation (2.26) reads, as follows:

-10-2 - (p ü) + | - (p ü2) + | - (p ü w) = - T ^ + f- ( TV + TC ) (2.28) öt v fm ' öx V Km ' öz V Km ' öx öz v m,xz m.xz' Ir- (P w) + | - (p ü w) + f- (p w2) = - 7T-S+ | - ( TV + ÏC ) + p g (2.29) ot m ' öx m ' öz V Km ' öz öx m,zx m,zx Kn ° in which: •c = t = H (-5— + -r—) (2.30) m,xz m,zx m öz öx T_{m,xz m,zx}C =xt = - (p-pc + p c)(ef)(|ii + Ist) - (p -p)(e0)(ü Is- + w Is-) (2.31) v v 's f' öz öxy v rs f/\"-s/v oz öx'

Equations (2.12), (2.13), (2.28) and (2.29) define a system of 4 equations with 5 unknown variables, being c, u, w, p and E . By relating the mixing co­ efficients to local flow variables applying a simple mixing length model or a more sophisticated K-Epsilon model (including terms for the solid phase), the system can be closed and solved using appropriate boundary conditions. The K-Epsilon model which is based on the application of transport equations for the kinetic energy (K) of the turbulence and its dissipation rate (Epsilon), of­ fers the most promising representation because the damping of the turbulence energy and hence the reduction of the mixing coefficients can be represented in a straightforward way. This latter approach was followed by DeVantler and Larock (1983).

In the present study Eqs. (2.28) and (2.29) are not considered further, be­ cause only conditions with relatively small sediment concentrations are inves­ tigated.

3. EQUILIBRIUM SEDIMENT TRANSPORT

3.1 Introduction

This Chapter presents a brief summary of a series of papers published in the Journal of Hydraulic Engineering of the American Society of Civil Engineers (ASCE):

Part I, Bed Load Transport Part II, Suspended Load Transport

Part III, Bed Forms and Alluvial Roughness

In the following sections the expressions relevant for the physical descrip­ tion of the 2D and 3D mathematical models (see Chapter 4 and 5) are presented. The original papers are presented in Appendix A.

3.2 Sediment transport modes

The transport of sediment particles by a flow of water can be in the form of bed-load and suspended load, depending on the size of the bed material par­ ticles and the flow conditions. The suspended load may also contain some wash load, which is generally defined as that portion of the suspended load which is governed by the upstream supply rate and not by the composition and pro­ perties of the bed material. Although in natural conditions there is no sharp division between the bed-load transport and suspended load transport, it is necessary to define a layer with bed-load transport for mathematical represen­ tation. Usually, three modes of particle motion are distinguished: (1) rolling and sliding motion or both; (2) saltation motion; and (3) suspended particle motion. When the value of the bed-shear velocity just exceeds the critical value for initiation of motion, the particles will be rolling and sliding or both, in continuous contact with the bed. For increasing values of the bed-shear velocity, the particles will be moving along the bed by more or less regular jumps, which are called saltations. When the value of the bed-shear velocity exceeds the fall velocity of the particles, the sediment particles can be lifted to a level at which the upward turbulent forces will be compa­ rable with or of higher order than the submerged weight of the particles and as a result the particles may go in suspension.

Usually, the transport of particles by rolling, sliding and saltatlng is called bed-load transport. For example, Bagnold (1956) defines the bed-load transport as that in which the successive contacts of the particles with the bed are strictly limited by the effect of gravity, while the suspended load

-12-transport is defined as that in which the excess weight of the particles is supported wholly by a random succession of upward impulses imparted by turbu­ lent eddies. Einstein (1950), however, has a somewhat different approach. Einstein defines the bed-load transport as the transport of sediment particles in a thin layer of 2 particle diameters thick just above the bed by sliding, rolling and sometimes by making jumps with a longitudinal distance of a few particle diameters. The bed load layer is considered as a layer in which the mixing due to the turbulence is so small that it cannot influence the sediment particles, and therefore suspension of particles is impossible in the bed-load layer. Further, Einstein assumes that the average distance travelled by any bed-load particle is a constant distance of 100 particle diameters, indepen­ dent of the flow conditions, transport rate and the bed conposition. In the view of Einstein, the saltating particles belong to the suspension mode of transport, because the jump lengths of saltating particles are considerably larger than a few grain diameters. This approach of Einstein is followed by Engelund and FredsiSe (1976).

In the present study, the approach of Bagnold is followed, which means that the motion of the bed-load particles is assumed to be dominated by gravity forces, while the effect of turbulence oiv the overall trajectory is supposed to be of minor importance. This latter statement means that the trajectory of a saltating particle can be somewhat wavy at certain locations due to the hig­ hest (turbulent) fluid velocities of the spectrum. However, the overall dimen­ sions of the trajectory are typically those of a saltating particle.

It is assumed that the (theoretical) maximum saltation height for a given flow condition can be determined from the equations of motion (excluding turbulent fluid forces) for a bed-load particle. As the start of a particle trajectory will be caused by random fluid forces, the height and length of a trajectory are random variables. However, by estimating the most unfavorable initial con­ ditions, a maximum saltation height can be computed. If for given flow condi­ tions there are sediment particles with a jump height larger than the (compu­ ted) maximum saltation height present in the flow, then these particles are assumed to be transported as suspended load. All particles with a jump height smaller than the maximum saltation height are transported as bed load.

3.3 Bed load transport

In this study the bed-load transport is defined as the transport of particles by rolling and saltating along the bed surface. The transport rate (st>) is de­ fined as the product of the particle velocity (ub), the saltation height (£.)

8b " Ub 6b % ( 3 a> Based on numerical solution of the equations of motion for a saltating par­

ticle and the analysis of experimental data, new relationships have been pro­ posed for the saltation height (6, ) , the particle velocity (ub) and the bed-load concentration (cb), being:

0.7 0.5 6b= 0-3 d5 0D* T ( 3 > 2 ) 0.5 0.6 ub = 1.5 (A g d5 Q) T (3.3) c = 0.117 D"1 T (3.4) b * in which:

ho

D = d (Ag/v2) = particle diameter (-)

K )

- (u* „ )

^3 d9 0

transport stage parameter (-)

(u.

V

*,cr

u* = (g0,5/C')ü = effective bed-shear velocity related

to the grains. (m/s)

uA = critical bed-shear velocity according

to Shields (m/s)

A = (p - p ) / p = r e l a t i v e density (-)

C' = 18 log(v~5—) = Chézy-coefficient r e l a t e d to grains ( m ° '5/ s )

depth-averaged flow velocity (m/s)

d , d = particle diameters of bed materical (m)

v = kinematic viscosity coefficient (m2/s)

g = acceleration of gravity (m/s2)

Substitution of Eqs.(3.2), (3.3) and (3.4) in Eq.(3.1) yields: 0.5 1.5 _0.3 2.1

-14-in which:

s^ = bed load transport rate (volume) (m2/s)

Equation (3.5) is valid for particles in the range of 200 to 2000 um. About 50% of the calibration data are in the range of 200 to 500 \im.

3.4 Suspended load transport

The equilibrium suspended load transport is described in terms of an equili­ brium sediment concentration profile (ce) and a logarithmic velocity (u), which after multiplication and integration yields:

h

s = ƒ u c dz (3.6) s,e e

* a in which:

s„ „ = equilibrium suspended load transport (volume) (m/s)

h = water depth (m) a = reference level (= upper limit of bed load layer) (m)

z = height above bed (m) The general expression for the equilibrium concentration profile is (see Eq.

(2.18)):

c w

+

<Ë. = o (3.7)

s, m s dz

in which:

ws m = effective fall velocity of a particle in a fluid-sediment (m/s) mixture

e = sediment mixing coefficient (m2/s)

The fall velocity of a particle in a fluid-sediment mixture shows considerable reduction in case of high sediment concentrations (hindered settling). This latter effect can be represented by:

w = (l-c)a w (3.8)

s,m s v '

in which:

ws = particle fall velocity in still clear water (m/s)

The fluid mixing coefficient, usually, is described by a parabolic distribu­ tion over the depth:

Ef = n" ( 1 ~ h} K u*h ( 3-9 )

In the present study a parabolic-constant distribution is used, mainly because it may give a better description of the concentration profile near the water surface (see also Section 4.4.3.1):

£f - « ï ï '1" ^ £f,max for ^ < 0.5 (3.10a)

ef " ef,max= °-2 5 *»*h for ^ > 0.5 (3.10b)

The sediniönt mixing coefficient is related to the fluid mixing coefficient by:

e = p * e. (3.11) s r

The (3-factor describes the difference in the diffusion of a discrete sediment particle and the diffusion of a fluid particle and is assumed to be constant over the depth. Based on the analysis of the data of Coleman (1970), the fol­ lowing expression has been proposed:

P = 1 + 2(w /u )2 for 0.1 < w /ut < 1 (3.12)

s * s *

The ^-factor expresses the influence of the sediment particles on the turbu­ lence structure of the fluid (damping effects) and is assumed to be dependent on the local sediment concentration. Based on theoretical and experimental work, the following function has been proposed:

0.8 o.i»

* = 1 + (c/cQ) - 2(c/cQ) (3.13)

in which:

(-)

For concentrations smaller than 2500 mg/1 (c=0«001) the ^-factor is in the range of 0.9 to 1 and may, therefore, be neglected.

Applying a parabolic-constant es~ distribution and assuming small sediment

concentrations (ws m = wg and $ = 1 ) , the equilibrium sediment concentration

1 6

-= [(-) ( r ^ ) ]

for f < 0.5

Z = w /(pKUj,) = suspension parameter (-) s *

c = equilibrium concentration at reference level (-) a,e

The sediment concentration profile can be computed when the reference concen­ tration is known. The following expression has been proposed for this latter variable:

d T1-5

c = 0.015 — - 5 - 3 (3.15)

a,e a D ° -J

Equation (3.15) specifies a dimensionless volume concentration; multiplying by p yields a value in kg/i • Equation (3.15) is valid for particles in the range of 100 to 500 \xm.

The reference level (z = a) is proposed to be equal to half the bed form height.

3.5 Bed forms and effective bed roughness

The morpological behaviour of an alluvial channel bed is rather complicated. The fundamental difficulty is that the channel bed characteristics (bed forms) and hence the hydraulic roughness depend on the flow conditions (flow velocity and depth) and sediment transport rate. The flow conditions are in turn, strongly dependent on the channel bed configuration and its hydraulic rough­ ness. For predictions it is of essential importance to know the relationships between the bed form characteristics, the hydraulic roughness and the flow conditions. Based on analysis of flume and field data, the following expres­ sions for the bed form dimensions (dunes) have been proposed:

j p - 0.11(d50/h)°'3(l - e "°'3T)(25 - T) (3.16)

\fa = 7-3 h (3.17)

in which:

A, = dune height (m) \ = dune length (m)

The effective hydraulic roughness height (ks) of the bed forms has been ex­

pressed, as:

k . - 3 do n+ l . l Ah( l - e b b) (3-18)

4. TWO-DIMENSIONAL MATHEMATICAL MODELLING OF SUSPENDED SEDIMENT TRANSPORT

4.1 Introduction

An overview of two-dimensional vertical mathematical models for suspended se­ diment transport is presented in Chapter 1.

This Chapter presents a detailed description of a two-dimensional vertical model (SUTRËNCH-2D) for suspended sediment particles (larger than about 50 urn) in non-stratified free surface flows developed by the writer (see also Van Rijn, 1986b).

The basic equations and the boundary conditions are presented. The applied numerical solution method is described briefly. Considerable attention is paid to the numerical inaccuracies of the solution method. The results of a de­ tailed sensitivity analysis concerning the influence of the basic physical pa­ rameters on the morphological processes are given. Three test cases to verify the mathematical model are considered. Finally, two special applications of the model are presented.

4.2 Simplified mass-balance equation for suspended sediment transport

Basically, the computation of the time-averaged sediment transport in the two-dimensional vertical plane requires the simultaneous solution of the mass-ba­ lance equations for the fluid and the sediment phases (see Eqs. (2.10) and (2.11)). Generally, the flow conditions in rivers and estuaries are such that the sediment concentrations are considerable smaller than about 5 kg/nH. po r

these conditions the influence of the sediment particles on the fluid phase is negligible small (Appendix A 2 ) . In that case it is sufficiently accurate to solve the mass-balance equation for the sediment phase only (Eq. (2.11)), which reads as (overbars indicating time-averaged values are omitted):

oc a . „ a ,

, a , ac,

"57

aï<

) + ay <

) + ai" (»-w

>

- ax" <

s,x ax"> +

" ^

s , y a ? - a r

, z ^ >

in which:

c = time-averaged sediment concentration (-) u,v,w = time-averaged fluid velocities in x,y,z directions (m/s)

ws = particle fall velocity (m/s)

e = sediment mixing coefficient (m2/s)

s

-19-A d e f i n i t i o n sketch i s given in Fig. 4 . 1 . surface z=z,

-»■ x

Fig. 4.1, Definition sketch

steady flow

Assuming steady state conditions and neglecting the horizontal diffusive transport components which are an order of magnitude smaller than the other transport components, as shown by Kerssens et al (1979), Eq. (4.1) reduces to:

5— (uc) + T — (vc) + T — (w-w )c ox K ' oy v ' öz s

T~ ( E T~)_{oz s,z oz} = 0 (4.2)

To represent the suspended sediment transport in gradually converging or di­ verging channels (or stream tubes), Eq. (4.2) is integrated in the y-direc-tion. The' flow is assumed to have a rectangular cross section (flow width b is a function of x only). Integration yields: y2 a y2 a Y2 y2 a ö c ƒ S T <U C> dy + / a ^v c ) d y + ' aT (w-ws)c d v " / ^ <E= ^ d v or 3_ 9x öz s öz dy dy ,

( ƒ (uc) dy) + (uc)1~ - (uc)2 -j^-+ (vc)2 - (vc)x+

y2 y2

+ f- ( ƒ (w-w )c dy) - £ - ( ƒ (e ^ ) dy = 0 oz y s öz y s öz

(4.3)

(4.4)

Assuming the fluid velocities, the mixing coefficients and the sediment con­ centrations to be approximately constant in tranverse direction, it follows that:

| j (buc) + | j - {b(w-ws)c} - |- (b e8 g ) + ((uCl) A _ ( v c > i } + dy,

+ {(uc)2 -^f-- (vc)2) = 0 (4.5)

Applying the kinematic boundary condition (zero transport normal to the boun­ dary), the latter two terms of Eq. (4.5) are zero, resulting in:

£ — (buc) + ~ (w-wjc - f- (eQ £7) = 0 (4.6)

b öx oz s' oz * s oz'

Equation (4.6) which is the basic equation of the SUTRENCH-2D model, can be solved numerically when the fluid velocities (u,w), the sediment mixing coef­ ficient (e ) , the fall velocity (wg) and the geometrical and physical boundary conditions are known.

The particle fall velocity is assume^ to be a constant, which implies that the proposed model is restricted to relatively narrow graded bed materials.

non-steady flow

Equation (4.6) can also be used to represent the sediment transport process in tidal flow by schematizing the tidal cycle to a number of quasi-steady flow periods.

Kerssens et al (1979) have shown that the unsteady concentration term (öc/öt) is relatively small for sediment particles larger than about 150 um transpor­ ted in a channel with a depth smaller than about 10 m. For these conditions the time lag between the integrated sediment transport and the depth-averaged flow velocity is relatively small and the sediment transport can be considered to be in a quasi-steady state during the tidal cycle.

4.3 Flow velocity field

Various models can be applied to describe the velocity field, depending on the complexity of the flow (bed-level gradients). In case of complicated flow con­ ditions including those with flow separation and flow reversal a refined math­ ematical approach is of essential importance.

The most universal and widely used mathematical model for these type of flows is the K-EPSILON model which is based on the equations of continuity and mo­ tion and two additional transport equations for the turbulence kinetic energy (K) and its dissipation rate (E) to describe the turbulent fluid shear stress (Rodi, 1980; Alfrink and Van Rijn, 1983). A disadvantage of the application of the K-EPSILON model Is the relatively large computation time needed to solve the complete set of equations. Consequently, the K-EPSILON model is not yet an attractive model for long-term morphological computations.

-21-Therefore, a simple PROFILE model was developed to compute the velocity field. This PROFILE model is based on the application of relatively simple profiles to describe the vertical distribution of the velocities. A disadvantage of the PROFILE model is the need for empirical data to calibrate the applied coeffi­ cients.

In this study the sophisticated K-EPSILON model as well as the simple PROFILE model are used to compute the velocity field. Flume experiments concerning the flow in a trench are used to compare the results of both models.

4.3.1 K-EPSILON model

Considering the flow in a channel with a constant width (b) and a plane water surface, the equations of continuity and motion for steady state conditions in the vertical plane read (overbars Indicating time-averaged values are omit­ ted): continuity o u Öw ,, „. T~ + TT = 0 (4-7) OX OZ v ' motion

t ( *

2

) + h <

uw

>

+

i h (p-e

R

,x> -

i I T

^

,XZ>

- ° <

4

-

8

>

| _

| _

i | _

.

_ i | _

^

u = time-averaged fluid velocity in x-direction (m/s) w = time-averaged fluid velocity in z-direction (m/s)

p = time-averaged static fluid pressure (N/m2)

pR = (Reynolds') turbulent fluid pressure (N/m2)

i = (Reynolds') turbulent fluid shear stress (N/m2) R

Basically, the Reynolds stresses represent a viscous part and a fluctuating (turbulence) part as follows:

ou PR.X " 2 p v ^ - p u ' u ' (4.10) P 2 p v — - p ^P~ (4-n> R,z r dz ÖU ÖW.

, =

=

< £ U £ > - P ^ (4-12)

Usually, the terms representing the turbulent fluctuations are approximated in a way analogous to the viscous stresses (Boussinesq hypothesis, see Rodi, 1980): ou p = - p W = 2 p E f - i p K (4.13) t,x f ox 3 p = - p -W^ Tr= 2 p e I"- - | p K (4.14) t,Z I OZ 1 ,ou öw. T = T = - p u'w' « p E j f + f ) (4.15) t.xz t,zx f 8z 9x' in which:

Pt = turbulence fluid pressure (N/m2)

T = turbulence fluid shear stress (N/m2)

ef = fluid mixing (or eddy viscosity) coefficient (m2/s)

K = turbulence kinetic energy (m2/s2)

Neglecting the viscous stresses, the equations of continuity and motion can be solved numerically when the fluid mixing coefficient is specified ("turbulence closure"). The turbulence closure implies the specification of a set of equa­ tions which relate the fluid mixing coefficient to additional variables of the main flow.

The most universal and widely used closure for complicated flows is the K-EPSILON model (see Rodi, 1980), which is a two-equation closure based on

transport equations for the turbulence kinetic energy (K) and its dissipation rate (E). The K and E variables are defined as:

K = \ ["(^y + c o

+ lyl

] <

-

)

, öu' 2 ou' 2 öw' 2 öw' 2 r, ^^

E - - v [ ( — ) + ( — - ) + ( — - ) + ( - — ) ( 4 . 1 7 ) öx Öz dx oz

The variables K and E are related to the eddy viscosity (Ef) by:

2

e

f = «v r-

in which:

c = turbulence constant (-)

The transport equations for the turbulence energy (K) and the dissipation rate (E) read:

-23-MH!Q.

I W _ I t ^ ö u _ ^ ^ | u _ a ff |K ^.

lf Ö K )

o

öx T <5z p tfx p 1>Ï öx '•a 3x' dz U oz-1 + b u ^ -l y^ K K

S(uE) Q(wE) E. , Pt,x du^ Tt,xz ^u öx öz IE K *• p öx p öz'

Öx V öx; öz V Öz' + C2E K U l*-^;

E E

in which C , C , o and o are universal constants. IE 2E K E

Assuming equal production and dissipation of turbulence energy (K) and neglec­ ting the convection of E, the transport equation for E reduces to (Rodi, 1980):

C1F = C2E " — ( 4 , 2 1 )

IE 2E c 0 > 5

E u

Launder (1972) has proposed a "standard" set of constants, which have been ob­ tained by computer calibration for plane jets and mixing layers:

C = 0.09, C = 1.44, C = 1.92, o = 1, a = 1.3, < = 0.435 |i IE 2E K E

In the present study the following values were used:

C = 0.09, C1 E = 1.6, C2 E = 1.92, aK = 1, ^ = 1.3, < = 0.35

which yield somewhat better results for hydraulic rough flow conditions (Alfrink and Van Rijn, 1983)

Equations (4.7), (4.8), (4.9), (4.18), (4.19) and (4.20) represent a set of six equations with six unknowns (u, w, p, e,, K and E) which can be solved nu­ merically applying an appropriate set of boundary conditions.

Computation results for the flow in a trench are shown in Fig. 4.4. More de­ tailed results have been presented by Alfrink and Van Rijn (1983).

4.3.2 PROFILE model

4.3.2.1 Longitudinal flow velocity

Coles (1965) showed that the velocity profiles in a non-uniform flow can be described by a linear combination of a logarithmic profile representing the law of the wall and a perturbation profile representing the influence of pres­ sure gradients.

In the present study a similar approach has been used, as follows:

u = A, u, In (z/z ) + A. u, F

1 h o Z n (4.22)

in which:

u = flow velocity at height z above bed u, = flow velocity at water surface

h A = dimensionless variable A = dimensionless variable 2 F = perturbation profile (m/s) (m/s) (-) (-) (-)

The perturbation profile (F) is represented by:

o t 2t

F = 2 Ti - n (4.23)

in which:

T) = (z-zo)/(h-zQ) (-)

Figure 4.2 shows the perturbation profile for various values of t. 1.0 0.5 t : 4 t =

2y

Fig. 4.2, Perturbation profiles

The A2-variable can be related to the A.-variable by applying the boundary condition, u = u^ for z = h resulting in:

-25-A = 1 - -25-Ax ln(h/zQ) (4.24)

Combining Equations (4.22), (4.23) and (4.24) yields:

u = A: uh ln(z/zo) + ",, (l " Ax ln(h/zQ)) (2 nt- n2 t) (4.25) The flow velocity profile, as described by Equation (4.25), is completely de­

fined when the unknown variables A^, t and uj, are specified. Therefore, three additional equations must be specified, which are:

• equation of continuity, • equation for the t-parameter,

• equation for the surface flow velocity ( uh) . continuity

The width-integrated discharge (Q) can be represented by: h

Q = b ƒ u dz (4.26) z

o

Substitution of Equation (4.25) in Equation (4.26) and integration yields:

Q = A1 (-1 + ln(h/zQ)) b h uh + (l - Ax ln(h/zQ)) ( 3 t + 1— ) b h uh (4.27)

2t2+3t+l

t-garameter

Analysis of flow velocity profiles measured in a trench perpendicular to the flow direction (Van Rljn, 1980a), showed that the mlddepth (z = 0.5h) velocity at each particular location is approximately equal to the equilibrium mlddepth velocity according to a logarithmic profile at that location.

Thus:

umiddepth=umiddepth, equilibrium (4.28)

The mlddepth velocity is described by Equation (4.25) resulting in:

um = Ax uh ln(0.5h/zo) + uh (l - Al ln(h/Z())) (2(0.5)' - (0.5)2t) (4.29) The mlddepth velocity for an equilibrium flow is described by a simple loga­

ln(0.5h/z )

° " (4.30) m > e -1 + ln(h/zo)

in which:

u = Q/bh = crossection-averaged flow velocity

Combining Eqs. (4.28), (4.29), (4.30) and (4.27), it follows that

t -

<

y i 3 t + i

o [2t + 3t + 1][2(0.5) -(0.5)

= 0.16t2 - 0.29t + 1.02 (4.31)

surface velocity

The surface velocity is described by a first order differential equation which yields an exponential adjustment of the surface velocity to the equilibrium surface velocity (u^ e) , as follows:

du, u, u, u,

n h . e h n n ...

ÏT =

i T T ^ ~

2 TT "

3 F"

-

in which:

u^ e = surface velocity for equilibrium flow (Eq.(4.34)), (m/s)

h = water depth (m) b = flow width (m) a ,<x ,a, = empirical coefficients (see section A.3.2.5) (-)

To demonstrate that the PROFILE method is capable of representing a wide range of velocity profiles (including flow reversal), some examples are shown in Fig. 4.3 for given values of Q,b,h and uh.

For uniform flow (9u/dx=0) the velocity profile can be expressed, as:

u

" InTWO

^l\->

-

o

Substitution of Eq. (4.33) in Eq. (4.26) and integration yields: ln(h/zn)

uh e » ( 4 , 3 4 )

* -1 + ln(h/zQ)

Equation (4.34) yields similar values as the more commonly used expression: g°-5ln(h/zo) o

u 2_5 (4.35)

-27-ln which:

C = 1 8 log(12h/kg) = Chézy-coefficient

K = constant of Von Karman

1.0 'S 0.5

/

!

'/ j

/

Si

I

ir

\

i

/

f

^

Fig. 4.3, Velocity profiles according to PROFILE model.

4.3.2.2 Vertical flow velocity

Applying the (width-integrated) equation of continuity 1 d(bu) 8w _ „

(m°'5/s) (-)

b öx dz (4.36)

the vertical flow velocity (w) can be computed as:

z,+z zh+ z

w - - ƒ |2-dz' - i ^ / u dz' (*-37)

ox b bx

z.+z z,+z b o b o

Substitution of Equation (4.25) and integration yields a (complicated) analy­ tical expression for the vertical flow velocity.

4.3.2.3 Bed-shear velocity

The bed-shear velocity (u*) is determined from the flow velocity (ufe) computed at a height z = 0.05 h above the bed assuming a logarithmic profile in the near-bed layer (z < 0.05 h) as follows:

u* " l'ri(O.Ö5ri7^y

For uniform flow it follows that:

(4.38)

u*,e " -1 + ln(h/z ) u

4.3.2.4 Computation procedure

(4.39)

Equation (4.32) can be solved numerically by using a Runge-Kutta method. The surface velocity at the inlet (u ) must be known. The complete set of Equa­ tions (4.27), (4.31) and (4.32) is now defined and can be solved to determine the A\, t and uj, variables. Using Equation (4.25), the velocity profiles can be computed at each location. The input data for the PROFILE model are: dis­ charge (Q), width (b) and depth (h) along the traject, effective bed roughness (ks), constant of Von Karman (K) and the surface velocity (u^ ) at the inlet.

4.3.2.5 Calibration of PROFILE model

The a and o coefficients have been determined by fitting computed and mea­ sured velocity profiles in a flow of constant width (b). The calibration data consisted of velocity profiles measured in trenches perpendicular to the flow direction. The bottom was covered with gravel material. In all, the data of six flume experiments have been used (Van Rijn, 1980a). The basic data are given in Table 4.1. test T6 T7 T8 T13 T14 T16 uo (m/s) 0.4 0.4 0.395 0.395 0.4 0.39 ho (m) 0.206 0.204 0.204 0.206 0.317 0.2 Ll (m) 0.4 0.4 0.4 0.8 0.8 1.6 L2 (m) 3.2 3.2 1.0 2.4 2.4 0.8 hi (m) 0.2 0.2 0.2 0.2 0.2 0.2 k (m) 0.006 0.011 0.006 0.006 0.006 0.011 -»K-

-»K-Table 4.1, Basic data of calibration tests

-29-The calibration procedure consisted of the computation of the velocity pro­ files in the trench for a given set of a and a coefficients. The computed velocity profiles were compared with the available measured velocity profiles.

If necessary, the a and a coefficients were adjusted and the procedure was repeated until the overall agreement between computed and measured velocity profiles was satisfactory. The procedure was carried out for all cases (Table 4.1). The a and a -coefficients were found to be dependent on the local bot­ tom slope (dh/dx). The following expressions yield the best results for all available data:

a = 0.28 + 0.11 tanh[6(dh/dx) - 0.15] (4.40)

a = 0.235 + 0.065 tanh[17(dh/dx - 0.035)] (4.41)

The a -coefficient represents the adjustment of the surface velocity to varia­ tions in the transverse direction. Since experimental data were not available, the a -coefficient could not be calibrated. Therefore, an expression is ap­ plied which yields a gradual adjustment of the surface velocity, as follows:

a = 0.1 tanh[10(db/dx)] (4.42)

This restricts the application of the model to situations with db/dx < 0.1. For a straight uniform flow (dh/dx = 0, db/dx = 0) it follows that a = a = 0.2 and a, = 0 resulting in

-(a1/h)x

u, - u, = (u, - u, ) e (4.43) h,x h,e h,o »,e

in which:

u = surface velocity at distance x (m/s)

uh 0 = surface velocity at inlet x=0 (m/s)

u = equilibrium surface velocity (m/s) n, e

<*! = 0.2 (-)

Assuming u, = 2 u, , the 90%-adjustment length is about 15h and the 99%-ad-h ,o 99%-ad-h,e

justment length is about 30h which are commonly observed values for the ad­ justment of the velocity profiles in a straight channel.

Figure 4.4 shows the calibration results of the PROFILE model for the T6-ex-periment. Velocity profiles computed by the sophiscated K-EPSILON model (Alfrink and Van Rijn, 1983) are also shown.

vertical velocity longitudinal velocity * 0.20 Ï 0- 10 ~ 0.00 '-0-10 -0-20 -0.30 9 10

/

A

\ h 1 _J

COMPUTED L0NG1TU0INAL FLOW VEL0C1TT COMPUTED LONGITUDINAL FLOW VEL0C1TT COMPUTED VERTICAL FLDW VELOCITY EQUIVALENT ROUGHNESS - 0.02 M VON KARMANN CONSTANT - 0.35 DISCHARGE - O.OB M2/S MEASURED OF K-E MODEL OF PROFILE MODEL OF PROFILE MODEL -01STANCE.XIM)

Fig. 4.4, Measured and computed velocity profiles in a trench (T6-experiment)

Flow separation and reversal phenomena can be observed in the deceleration zone of the trench (locations 2,3 and 4 ) . These effects are represented by the K-EPSILON model and the PROFILE model, but on a somewhat smaller scale. The calibrated PROFILE model yields the largest reversed flow velocities (location 3 ) . Some three-dimensional effects due to the relatively small flume width can be observed in the deceleration zone of the trench, as indicated by the unit width did occur (in centre line of flume) which Is not constant in the dece­ leration zone. At locations 6,7 and 8 both models produce similar results. In the acceleration zone of the trench (locations 9,10 and 11) the PROFILE model shows somewhat better results in the near-bed region.

Figure 4.5 shows computed and "measured" bed-shear velocities. The "measured" values have been derived from the flow velocities measured in the lowest mea­ suring point near the bottom assuming a logarithmic velocity profile near the bottom and a known bottom roughness (k ) . The computed bed-shear velocities are obtained in a similar way using the velocities computed in the lowest grid points.

3 1 -■ / " a 0.04

i™

- •

V

p;

<

S>"

T^>

ED

computed K-E model

£?■■'computed PROFILE model

x measured

F i g . 4 . 5 , M e a s u r e d and c o m p u t e d b e d - s h e a r v e l o c i t i e s ( T 6 - e x p e r i m e n t )

4.3.2.6 V e r i f i c a t i o n of P R O F I L E m o d e l

CJf.'ïjr.'.ïïy*.

F i g u r e 4.6 shows c o m p u t e d a n d m e a s u r e d v e l o c i t y p r o f i l e s for a n o t h e r e x p e r i ­ m e n t , n o t used i n the c a l i b r a t i o n p r o c e d u r e . T h e r e f o r e , these r e s u l t s c a n be .'.'caJH-.tai.e': c o n s i d e r e d a s an i n d e p e n d e n t v e r i f i c a t i o n of the P R O F I L E m o d e l . T h e e x p e r i m e n t SI13 \<i b3 c o n s i s t e d o f flow o v e r a d a m f o l l o w e d b y a t r e n c h . The w i d t h i s c o n s t a n t . T h e &;.T . a l a r d e p t h - a v e r a g e d v e l o c i t y at the i n l e t (x = 0 ) is 0.185 m / s . T h e e f f e c t i v e bed HO' js.iol.'

r o u g h n e s s is k = 0.01 m. T h e b o t t o m s l o p e s of the d a m a n d the t r e n c h are 1 : 2 . riB.-i -iJili

A s can b e o b s e r v e d , the a g r e e m e n t b e t w e e n c o m p u t e d and m e a s u r e d v e l o c i t i e s is

•.'.'TW :iJ!J

rather good at the locations 2 , 3 , 4 , 9 , 1 0 and 11 showing maximum deviations less -3.'iSb 31 t h a n a p p r o x i m a t e l y 2 0 % . 't'JOP' A t the l o c a t i o n s 5 , 6, 7 , 1 2 a n d 13 the d e v i a t i o n s b e t w e e n c o m p u t e d a n d m e a ­ sured v e l o c i t i e s at s p e c i f i c p o i n t s a r e c o n s i d e r a b l y l a r g e r , u p to 1 0 0 % . A s r e g a r d s v e r i f i c a t i o n for f i e l d c o n d i t i o n s , o n l y s o m e v e l o c i t y m e a s u r e m e n t s in -1 t h e m i d d l e zone of a p i p e l i n e t r e n c h i n a w i d e e s t u a r y i n t h e N e t h e r l a n d s w e r e a v a i l a b l e . C o m p a r i s o n of c o m p u t e d a n d m e a s u r e d p r o f i l e s for t h i s s i t u a t i o n s h o w s r e a s o n a b l e a g r e e m e n t ( V a n R i j n , 1 9 8 5 b ) . M o r e v e r i f i c a t i o n f o r f i e l d c o n -f>. d i t i o n s i s n e c e s s a r y .

© © © ©

T~7T

© ©

®

®

'

r

r'

*

•

®

•:

_!

1

0

*

©

V .

i

/•

'

©

J.

ƒ

/

/

computed longitudinal tlow velocities [PROFILE - model)

Fig. 4.6, Computed and measured velocity profiles

4.4 Fluid and sediment mixing coefficients

4.4.1 Introduction

The eddy viscosity concept is applied to represent the transfer of fluid mo­ mentum and sediment mass. In the present study the coefficients expressing the transfer of fluid momentum and sediment mass are shortly called fluid and se­ diment mixing coefficients.

By definition the fluid mixing coefficient reads as (Hinze, 1975):

E

f

= Ï Ï

? 7

in which:

(4.44)

(m/s) u' = velocity fluctuation in i direction

lj = mixing length scale in j direction normal to 1 direction (m)

Generally, the fluid mixing coefficient is a tensor of the second order.

The sediment mixing coefficient is related to the fluid mixing coefficient, as follows (see Eq. (3.11):

-33-e = p * -33-e, (4.45) s f

in which:

p = proportionality factor related to the difference in the transfer of

fluid momentum and sediment mass (-) $ = turbulence damping factor (-)

Two approaches have been used to compute the spatial distribution of the fluid mixing coefficient: (1) the two-dimensional vertical K-EPSILON model and (2) the one-dimensional PROFILE model.

4.4.2 K-EPSILON model

A detailed description of the K-EPSILON model is given in section 4.3.1. The computed fluid mixing coefficients for one test case are shown in Fig. 4.8.

4.4.3 PROFILE model

4.4.3.1 Vertical distribution of fluid mixing coefficient

Assuming a linear shear stress and a logarithmic velocity distribution over the depth for equilibrium flow (du/ox = 0 ) , the vertical distribution of the fluid mixing coefficient is represented by a parabolic function, as follows:

e , - e, - e, (l - %r) (4.46)

f f,max f,max *■ h ; in which:

e,_ = 0.25 K u.h = maximum value of fluid mixing coefficient (m /s) f.max *

Based on the experimental data of Coleman (1970), Kerssens (1974, 1977) intro­ duced the parabolic-constant sediment mixing coefficient distribution; para­ bolic in the lower half of the depth and a constant value in the upper half of

the depth. This approach has been adopted by van Rijn (see Appendix A2) to mo­ del the fluid mixing coefficient, as follows (see Fig. 4.7):

£f " £f .max " ef ,max ^ " W f o r ÏÏ < °'5 ( A"4 7 a )

e. - ef for ^ > 0.5 <*-«b)

f f.max h

Equations (4.46) and (4.47) yield different values for the fluid mixing coef­ ficient in the upper half of the flow depth.