- Email: [email protected]
Finite element current and sediment transport modelling
ContinentalShelfResearch,Vol. 13, No. 8/9,pp. 891-902, 1993.
0278-4343/93$6.00+ 0.00 © 1993Pergamon Press Ltd
Printedin GreatBritain.
Finite element current and sediment transport modelling* T. UTNESt (Received 25 September 1990; in revised form 26 May 1992; accepted 11 June 1992) Abstract--A three-dimensional circulation model is presented and applied to compute suspended sediment transport. The turbulence is modelled by the use of a two-equation (k, e) model. Examples are shown for some test cases where comparisons are being made with simplified analytical solutions. A realistic example is finally presented where suspended sediment transport is computed in a basin with river in- and outflow.
1. I N T R O D U C T I O N
SEDIMENTre-entrainment and transport in coastal and shallow waters is an important topic which is receiving increasing environmental concern. Among the main problems is the transport of pollutants bound to suspended sediments. In order to gain more knowledge about such processes, several numerical models have been developed in recent years trying to predict such flows. These models range from one-dimensional vertical profile formulations (e.g. DE VANTIERand NARAYANASWAMY, 1989; HAGATUNand EIDSVIK, 1986) via two-dimensional vertical or depth integrated models (e.g. VAN Rim, 1987 and refs therein; TEISSON and FRISCH, 1988; VEERAMACHANENIand HAYTER, 1988) to quasi- or full three-dimensional formulations (SHENG and BUTLER, 1982; VAN RIJN and MEIJER, 1988; O'CONNOR and NICHOLSON, 1988; EIDSVIKand UTNES, 1991; and others). Recent literature on presentations and applications of such models may be found, for example, in WANG (1989). The use of three-dimensional models is especially relevant in situations involving flows produced by wind action, flow separation or density currents. In such cases it is necessary to use some kind of three-dimensional description in order to compute the vertical distribution and the bottom stress of the flow. The present paper describes a three-dimensional formulation based upon the hydrostatic pressure assumption. The turbulence is represented by the use of a two-equation (k, e) model. The numerical formulation is a finite element method combined with a weak form of the characteristics method. This combination is applied in order to compute the advective part of the flow accurately. The paper presents some test examples together with a realistic computation of sediment transport in a basin.
*This paper was presented at JONSMOD '90, Birkenhead, U.K. tNHL, Kl~ebuveien 153, 7053 Trondheim, Norway. 891
892
T. UT~ES
2. THEORETICAL FORMULATION 2.1. Model equations By application of the hydrostatic pressure assumption, the horizontal momentum equations can be written as
aUhat
~- V - (UUh) = - - f
o_/Av au / a z ] -~- V .
X U h -- V(I) q- 0Z \
(AhVUh)
(1)
and the continuity equation
0u + av + a~ = 0. ax ay Oz
(2)
Here Uh = (U, V) is the horizontal part of the three-dimensional velocity vector u = (u, v, w), f i s the Coriolis parameter, Av and A h are the vertical and horizontal eddy viscosities, respectively. The hydrostatic potential ¢~ is defined by qb = Patm + g~ _ PO -- P dz P0 z Po
(3)
where Patm is the atmospheric pressure, ~ is the surface elevation, p is the fluid density and Po is a reference density. The surface elevation is found by solving for the vertically averaged form of the equations (1) and (2), which take the form -
a~f + V. (fi-fi) = - f x fi - gV(~ + b) at
Vpatm + (Xs - - 'Tb) + V. (AhVfi) po poll
-
-
-
0~ + V-(H'a) = 0. at
(4) (5)
In these equations fi denotes the vertically averaged velocity, H is the total water depth, % and % are the surface and bottom stresses, respectively, and b is the baroclinic pressure contribution defined by b=
Po - P dz' dz. -h
(6)
P0
As pointed out by MELLOR and BLUMBERG(1985), the horizontal diffusion will be of minor importance for models that resolve the dominant horizontal eddy scale or variable bottom topography. This is assumed to be the case here, and only the vertical eddy viscosity Av is therefore considered in the following. This quantitity is modelled by the use of a (k, e) turbulence model as follows: Av = C~ k2
(7)
ok + V.(uk) = a (Avak /
at
~ak
az/+ P+ G - ~
(8)
Sediment transport modelling
893
O.f_6.~_V . ( u E ) : L ( A v O E I . . ~ _ C I ~ ( p _ } _ C 3 G ) _ C 2 ~ Ok Oz \ ~e Ok/
(9)
where the conventional values of the constants are (Rom, 1980): (C~,, Ok, ae, C1, Ce, C3) = (0.09, 1.0, 1.3, 1.44, 1.92, 0.8). The shear production and buoyant destruction terms in (8) and (9) are defined by
\~Oz/
\Oz/ /
G = ag A---Y-~O__C_C (7c OZ
(11)
a - p~ - P0 P0
(12)
and the expansion coefficent is
The sediment concentration is computed from the transport equation
OC + V - ( u C ) = 0 (avOC] 0%C Ot ~z \-~ ~z / + --Oz
(13)
where ws is the vertical sediment fall velocity, measured positive downwards. The sediment velocity is thus assumed to be equal to the fluid velocity except for the vertical direction, where a constant slip velocity equal to the particle fall velocity (ws) is assumed. The Prandtl number in (13) is set equal to o¢ = 1. The density is given by the state equation p = po(1 + aC).
(14)
At the present stage the stratification effects in the turbulence model have been neglected in the computations. It is recognized, however, that such effects should be included for cases with high sediment concentration (see e.g. DE VANTIERand NARAVANASWAMY, 1989). 2.2. Boundary conditions The conditions applied at the different types of boundaries are as specified below. Free surface conditions: 0Uh =
ooAv~
~s
(15) _
0¢
OX
Oy'
OS -k- U 0¢ -+- V- -
W = Ot a v OC + WsC = 0 Oz
k = ~-C-~,
P ----Patm
(16) (17)
e = --'KA
(18)
Here r is the von Karman constant, and us. is the surface friction velocity defined by
894
T. UTNES
Us* =X/'rs/P0. The conditions for (k, e) are equilibrium boundary layer conditions which are to be applied at a small distance A below the surface. Bottom conditions: Bottomconditions: A Ouh
v -~z = %
(19) Oh + v
W=-- UOX C = Cb.
(20) (21)
The bottom stress in (19) is given by % = P0fblUb]Ub,where Ub is the near-bottom velocity and fb is the bottom friction coefficient. The bottom boundary conditions for (k, e) correspond to the boundary layer conditions applied at the surface. The condition for C in (21) is given as a local equilibrium bed concentration Cb. Following HAGATUNand EIDSVIK (1986) we assume C b = Cmaxfor S > $ 2 , C b = 0 for S <- S I , and a linear matching for $1 -< S -< S2, where $1 and $2 are given limits of the Shields number at the bottom. By introducing the sediment particle size D, the actual Shields number is defined by s =
I bl
(22)
(Pc - po)gD
Lateral conditions: At the lateral boundaries either ~ or O~/On is specified, while the conditions for u may be of similar form or "free", i.e. computed from the equations with known ~. For k, e and C the lateral conditions are usually specified at other lateral boundaries.
3. NUMERICAL FORMULATION In order to solve equations (1) and (2) a mode splitting procedure is introduced. The total velocity is split in U = U + Ur
(23)
where fi is the vertically averaged velocity and Ur is the residual part. The equations for Ur are formally obtained by substitution of equation (23) into the governing equations (1) and (2). The general solution procedure for each time step is as follows: The vertically averaged equations for (fi, ~) are solved first. These values are substituted into the residual equation system which can then be solved for Ur. The total velocity field is finally obtained from equation (23) and the continuity equation. A finite element formulation is applied to solve these equations. The domain is divided into three-dimensional trilineal elements, with a node configuration that results in perfectly vertical element sides. The projection of these elements in the horizontal plane is used to solve the vertically averaged problem. The numerical formulation is described in more detail elsewhere (UTNES, 1991), and only a short resume is given below. The vertically averaged equations are considered first. A weighted residual formulation in space and time of the momentum equation gives
Sedimenttransportmodelling
895
In
where fa represents all the terms in equation (4) which have not been expressed explicitly, is the space domain and (tn, t,,+l) is the time increment. After integration by parts, and by introducing the following restriction on the weight functions oa + (~. V)Q = 0 ot
in (t,, t,,+l)
(25)
it can be shown that equation (24) reduces to
The choice (26) implies that the weight functions are transported along the streamlines. At time t~+l these functions are chosen equal to the interpolation functions for g, and at t~ they are found by use of equation (26). The time integration can be performed by any suitable method. However, the boundary integral is computed at tn if the normal velocity un is not specified equal to zero. The continuity equation is expressed in ordinary weighted residual form as follows by introducing a two-level time discretization I q~"+Idlq=I
f~
q ~ n d f ~ - I qAt(OVH'fi~+l+(1-O)VH~n)dfl
f~
f~
(27)
where q represents space-dependent weight functions The solution strategy for these equations is as follows: A two-level time discretization and a time-splitting of the momentum equation is applied. The continuity equation is transformed to a Helmholtz equation for the elevation by performing some manipulation (see Ua~Es, 1991), and a mixed interpolation of the variables is used: biquadratic velocity and bilinear elevation. The three-dimensional residual velocity ur is considered next. The generalized weighted residual formulation of the horizontal part of the momentum equation (1) is given by ,n
\ Ot + V . (UUh) + fb df~ dt = 0
(28)
where fb denotes all terms which have not been written out explicitly. Application of the weak form of the characteristics method gives, in analogy with the vertically averaged momentum equation
/n÷l
[df~+l
By introducing the relation (24) and substituting for ~, this equation is used to solve for the horizontal part of the velocity Ur. A two-level time discretization is again applied, and the vertical diffusion is treated implicitly, while all other terms are either known or explicitly treated quantities. The vertical component of the velocity is found by use of the continuity equation. The interpolation of the velocity is biquadratic in the horizontal and linear in the vertical plane. The equations for the scalar variables (k, e, C) are solved in the same way as the
896
T. Ua'NES
ac
Vt ~
ap =
all= aX aC ~x
= - WsC
nst
i
C=C b Zo
Fig. 1. Domain and boundary conditions for computation of sediment concentration profile.
momentum equation. For example, the generalized weighted residual formulation for the sediment concentration equation becomes
+
I7
and the (k, e) equations are treated similarly. 4. COMPUTATIONAL RESULTS In the following two simplified test examples are presented together with a more complex computation of suspended sediment transport in a jet-like flow through a basin.
4.1. Sediment concentration profile This test is defined by a horizontally uniform pressure-driven flow over a flat bottom with specified constant sediment concentration, see Fig. 1. If we assume a parabolic vertical eddy viscosity of the form Av = Kz(1 - z]h)Ub,, where Ub* = ~ / ~ is the bottom friction velocity, the concentration profile is given by the analytical solution
C=
c {Zb(h - z)lW:( u ")
b!z( h _ Zb))
(31)
where z0 is a small distance from the bottom to the first grid point. This solution is shown in Fig. 2 together with model computations for Zb/h = 0.01 and Ws/Ub. = 0.4. In order to evaluate the numerical results, computations have been performed both with the parabolic eddy viscosity and with the (k, e) model. It is seen from Fig. 2 that the numerical result with parabolic viscosity is fairly close to the analytical solution. This is of course what should be expected if the numerical scheme is correct. The result from the (k, e) model also shows a general trend which is similar to the previous solution, but the sediment concentration is now somewhat lower above the bottom, except at the free surface. At this boundary the analytical solution gives zero concentration, while the turbulence model predicts a final value due to the free surface boundary conditions for the (k, e) model.
Sediment transport modelling
897
1.0 ---.......~=~,,~,,,,~. z
0.8
h
0.6
\\~. \X\ x ~ . ,
0.4
-,N
0.2
\\
\
0
I 10 -3
10 -4
L 10"2
I
10 -1
10 0
C/C max Fig. 2. Concentration profile. Zb/h = 0.01, Ws/U, = 0.4. - Analytical solution with parabolic eddy viscosity; - - . - - Numerical solution with parabolic eddy viscosity; - - - Numerical solution with (k, e) model.
Computations with other values of the input variables (bottom concentration, fall velocity etc.) show similar correspondence between the analytical and numerical results. The computations were performed on a grid with 18 points in the vertical direction, with fine resolution close to the bottom.
4.2. Horizontal concentration adjustment This problem is defined by a horizontally uniform flow over a fiat sediment bed, with zero upstream sediment concentration and a specified constant bed concentration downstream, see Fig. 3. Initially there is no suspended sediment, and the problem is to find the steady state concentration adjustment downstream. Again it is possible to derive the analytical solution if the eddy viscosity profile is assumed to be parabolic (cf. VAN RIJN, 1987). The analytical solution is shown in Fig. 4 for the input variables Zb/h = 0.05 and ws/ub* = 0.2, together with corresponding numerical computations have been performed with the parabolic eddy viscosity and with the (k, e) model. A different numerical solution with the actual parameters is also given by VANRIJN (1987) for the case with parabolic eddy viscosity. The computations were performed on a grid with 18 points in the vertical and 31 points in the horizontal direction, with grid stretching in both directions. Our numerical solution with parabolic eddy viscosity is seen to be close to the
a___Ru=
ax
c:o
[
_
I / / L ~ / / / / / / / / / /VT- 7 - / - - / ~ / - - / - - /
onSt
....
C=Cb
~
--
-
-
/ |
7- 7 - 7 - 7- 7- -7 - - / - - [ - - / - - 7 - T - /
rigid bed sediment bed Fig. 3. Domain and boundary conditions for computation of horizontal concentration adjustment.
898
T. UTNES
1.0 C / C m a x = 0 ., 0 5~ ~ . ~
0.8 Z
h 0.6 ~
0.4
•
~ ~-,-- -
-
~ ~.
-~ - ~ -
. . . .
0.2
0 0
~ 10
i 20
I
30
40
x/h
Fig. 4. Horizontal adjustment, zb/h = 0.05, Ws/U, = 0.2. Analytical solution with parabolic eddy viscosity; - - . - - Numerical solution with parabolic eddy viscosity; - - - Numerical solution with (k, e) model.
corresponding analytical solution (Fig. 4). This result also agrees well with VAN RIJN'S (1987) computations for the same case. The model result with the (k, e) turbulence model shows a similar evolution, but now with a general trend of somewhat lower concentration values above the bottom. This tendency with a certain reduction in the concentration values is similar to the previous result (cf. Fig. 2). In summary, these tests show that the numerical model is able to solve such simplified problems fairly accurately with a prescribed eddy viscosity. The inclusion of a twoequation turbulence model modifies these results quantitatively, but the general trend is the same.
4.3 Suspended sand transport in a jet-like basin flow Here we consider a jet-like river flow through a basin, and look for possible suspension and transport of bottom sediments. The geometry is a nearly circular basin with river inand outflow, as shown in Fig. 5. The depth is constant and equal to 10 m, and the uniform in/outflow velocity is 1 m s-1. Along the land boundary the velocity is set equal to zero, and a small horizontal diffusion Ah = 4 m 2 s -1 is included. For problems with significant non-linearities, such as the present one, a certain minimum horizontal diffusion must be included in order to avoid numerical instabilities. The in/outflow conditions for (k, e, C) are zero normal derivatives. The surface and bottom boundary conditions are as given by equations (10) and (11) with zero surface stress, and the sediment is specified by p = 2500 kg m -3, D = 2 x 10 -3 m, Cmax = 0.1, w s = 10 -4 m s -1. The limiting values of the Shields number, $1 and $2, have been chosen to 0.10 and 0.75, respectively. Figures 6 and 7 show the computed surface and near bottom velocity, respectively. It is seen that the jet-like inlet flow is passing through the basin, and recirculation eddies have developed at each side of the geometrical symmetry line. This characteristic jet-and-eddy structure of the flow is mainly a nonlinear effect (UTNES, 1991). The suspended sediment concentration near the bottom is shown in Fig. 8, and the corresponding values along the transversal midsection is given in Fig. 9. As expected, the main part of the suspended sediments are found in the jet-like domain of the current, where the velocity gradients are strong. This is
899
Sediment transport modelling
E
///1
l\\\
l
8 a 4500 Fig. 5.
I
I
- ~000
-500
I///
I
I
0
~;00
I
t000
I rn
4~00
Horizontal grid domain for river-basin flow. T h e nodes are marked on one of the elements.
especially the case in the vicinity of the in- and outflow areas. Generally, the vertical velocity variations are small, and the main advantage of using a three-dimensional model here is to compute more accurately the velocity gradients near the bottom. In many practical situations the addition of wind effects will be of interest, and in such cases the application of a three-dimensional model would be even more important (cf. EIDSVIKand UTr~ES, 1991). 5.
CONCLUSIONS
A three-dimensional hydrodynamic circulation model has been generalized to include turbulence and sediment models in order to compute sediment transport problems. The hydrodynamic part has previously been tested on some basic problems (UTNES, 1991), and the present paper describes some computations with sediment suspension and transport involved. Two simple tests have been performed, both compared with simplified analytical solutions. An example has also been shown for a somewhat more complicated problem, although in that case there are no quantitative results to compare with. However, certain effects may be judged qualitatively, and the result seems to be reasonable. It is recognized that for problems with a relatively high sediment concentration, the sediment and turbulence models should be generalized to include stratification effects. Such effects can be included in the present model.
900
T. UTNES
1 mls I
.I
E
I
~2000
-
1000
Fig. 6.
0
I000
m
2000
Computed velocity field at the free surface.
o o
1 rn/s
o-
Q O
o-
f'--'"~
t l t ~' . . . . .
)
c3 o
o-
\ \ "". " . .' . .1 } ltl liI ~\ ' ' •
",
,
,
s
l
~
. . .... . . / / \
~
.
t
t
.
o o o
#2000
1 -
0
t 000
F i g . 7.
Computed
t 000
m
velocity field near the bottom.
2000
901
Sediment transport modelling
o,I E
. " / I /F/I/
\\\'x~,'~ \
\\\\\\\'\
,."/1/1111 i
,; /71~DII~ lit/if/ ',, (
i _
'
if
/ II1
.lit
k/kXX\ \
I JJ)J)/~, '
\\
I
tilt\ \ \ I III I/I,,"
I
{000
,
I1111}/,,
~-\
--
t
°ltt\l\t
o
~tSO0
t tt
I IIIlll
",
Fig. 8.
III\\\\",,
U///I/
,,"
_
///
I
-SO0
I
0
I
SO(}
I
I000
ISO0
m
Computed suspended sediment concentration in the near-bottom domain. The iso-values are: (0.0025, 0.005, 0.01, 0.02, 0.03, 0.04, 0.05).
E
9'.0025
?-
O o
[
- 4000
Fig. 9.
I
-,gO0
I
I
0
,gO0
I
4000
I
03
4S00
Concentration isolines along the transversal m i d section.
902
T. Uan~ES
Acknowledgements--This work was partially supported by Konsesjonsavgiftsfondet (Fund of License Fees).
REFERENCES DE VANTIER B. A. and R. NARAYANASWAMY(1989) A suspended sediment flow model for high solids concentration using higher order turbulence closure. Advances in Water Resources, 12, 46-52. EIDSVIK K. J. and T. UTNES (1991) A model for sediment reentrainment and transport in shallow basin flows. Coastal Engineering, 15, 247-256. HA~ATUN K. and K. J. EIDSVIK(1986) Oscillating turbulence boundary layer with suspended sediments. Journal of Geophysical Research, 91, 13045-13055. MELLOR G. L. and A. F. BLUMnERG (1985) Modelling vertical and horizontal diffusivities with the sigma coordinate system. Monthly Weather Review, 113, 1379-1383. O'CONNOR B. A. and J. NICHOLSON(1988) A three-dimensional model of suspended particle sediment transport. Coastal Engineering, 12, 157-174. RODI W. (1980) Turbulence models and their application in hydraulics. IAHR State-of-the-Art Paper, Delft. SHENG Y. P. and H. L. BUTLER(1982)Modelling coastal currents and sediment transport. Proceedings of the 18th Coastal Engineering Conference, Cape Town, 14-19 November, pp. 1127-1148. TEISSON C. and D. FRISCH (1988) Numerical modelling of suspended sediment transport in the Loire estuary. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May-1 June, pp. 14-22. UTNES T. (1991) Finite element modelling of quasi-three-dimensional nearly horizontal flow. International Journal of Numerical Methods in Fluids, 12, 559-576. VAN RUN L. C. (1987) Mathematical modelling of morphological processes in the case of suspended sediment transport. Dr Thesis, Technical University Delft, Delft, 135 pp. and 2 appendices. VAN RIJN L. C. and K. MEIJER(1988) Three-dimensional mathematical modelling of suspended sediment transport in currents and waves. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May-1 June, pp. 89-99. VEERAMACHANENIR. and E. J. HAYTER(1988) Mathematical Modelling of sediment transport at tidal inlets. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May1 June, pp. 14-22. WANG S. S. Y. (1989) Preprints of International Symposium on Sediment Transport Modeling. New Orleans, LA, 14-18 August, 829 pp.
0278-4343/93$6.00+ 0.00 © 1993Pergamon Press Ltd
Printedin GreatBritain.
Finite element current and sediment transport modelling* T. UTNESt (Received 25 September 1990; in revised form 26 May 1992; accepted 11 June 1992) Abstract--A three-dimensional circulation model is presented and applied to compute suspended sediment transport. The turbulence is modelled by the use of a two-equation (k, e) model. Examples are shown for some test cases where comparisons are being made with simplified analytical solutions. A realistic example is finally presented where suspended sediment transport is computed in a basin with river in- and outflow.
1. I N T R O D U C T I O N
SEDIMENTre-entrainment and transport in coastal and shallow waters is an important topic which is receiving increasing environmental concern. Among the main problems is the transport of pollutants bound to suspended sediments. In order to gain more knowledge about such processes, several numerical models have been developed in recent years trying to predict such flows. These models range from one-dimensional vertical profile formulations (e.g. DE VANTIERand NARAYANASWAMY, 1989; HAGATUNand EIDSVIK, 1986) via two-dimensional vertical or depth integrated models (e.g. VAN Rim, 1987 and refs therein; TEISSON and FRISCH, 1988; VEERAMACHANENIand HAYTER, 1988) to quasi- or full three-dimensional formulations (SHENG and BUTLER, 1982; VAN RIJN and MEIJER, 1988; O'CONNOR and NICHOLSON, 1988; EIDSVIKand UTNES, 1991; and others). Recent literature on presentations and applications of such models may be found, for example, in WANG (1989). The use of three-dimensional models is especially relevant in situations involving flows produced by wind action, flow separation or density currents. In such cases it is necessary to use some kind of three-dimensional description in order to compute the vertical distribution and the bottom stress of the flow. The present paper describes a three-dimensional formulation based upon the hydrostatic pressure assumption. The turbulence is represented by the use of a two-equation (k, e) model. The numerical formulation is a finite element method combined with a weak form of the characteristics method. This combination is applied in order to compute the advective part of the flow accurately. The paper presents some test examples together with a realistic computation of sediment transport in a basin.
*This paper was presented at JONSMOD '90, Birkenhead, U.K. tNHL, Kl~ebuveien 153, 7053 Trondheim, Norway. 891
892
T. UT~ES
2. THEORETICAL FORMULATION 2.1. Model equations By application of the hydrostatic pressure assumption, the horizontal momentum equations can be written as
aUhat
~- V - (UUh) = - - f
o_/Av au / a z ] -~- V .
X U h -- V(I) q- 0Z \
(AhVUh)
(1)
and the continuity equation
0u + av + a~ = 0. ax ay Oz
(2)
Here Uh = (U, V) is the horizontal part of the three-dimensional velocity vector u = (u, v, w), f i s the Coriolis parameter, Av and A h are the vertical and horizontal eddy viscosities, respectively. The hydrostatic potential ¢~ is defined by qb = Patm + g~ _ PO -- P dz P0 z Po
(3)
where Patm is the atmospheric pressure, ~ is the surface elevation, p is the fluid density and Po is a reference density. The surface elevation is found by solving for the vertically averaged form of the equations (1) and (2), which take the form -
a~f + V. (fi-fi) = - f x fi - gV(~ + b) at
Vpatm + (Xs - - 'Tb) + V. (AhVfi) po poll
-
-
-
0~ + V-(H'a) = 0. at
(4) (5)
In these equations fi denotes the vertically averaged velocity, H is the total water depth, % and % are the surface and bottom stresses, respectively, and b is the baroclinic pressure contribution defined by b=
Po - P dz' dz. -h
(6)
P0
As pointed out by MELLOR and BLUMBERG(1985), the horizontal diffusion will be of minor importance for models that resolve the dominant horizontal eddy scale or variable bottom topography. This is assumed to be the case here, and only the vertical eddy viscosity Av is therefore considered in the following. This quantitity is modelled by the use of a (k, e) turbulence model as follows: Av = C~ k2
(7)
ok + V.(uk) = a (Avak /
at
~ak
az/+ P+ G - ~
(8)
Sediment transport modelling
893
O.f_6.~_V . ( u E ) : L ( A v O E I . . ~ _ C I ~ ( p _ } _ C 3 G ) _ C 2 ~ Ok Oz \ ~e Ok/
(9)
where the conventional values of the constants are (Rom, 1980): (C~,, Ok, ae, C1, Ce, C3) = (0.09, 1.0, 1.3, 1.44, 1.92, 0.8). The shear production and buoyant destruction terms in (8) and (9) are defined by
\~Oz/
\Oz/ /
G = ag A---Y-~O__C_C (7c OZ
(11)
a - p~ - P0 P0
(12)
and the expansion coefficent is
The sediment concentration is computed from the transport equation
OC + V - ( u C ) = 0 (avOC] 0%C Ot ~z \-~ ~z / + --Oz
(13)
where ws is the vertical sediment fall velocity, measured positive downwards. The sediment velocity is thus assumed to be equal to the fluid velocity except for the vertical direction, where a constant slip velocity equal to the particle fall velocity (ws) is assumed. The Prandtl number in (13) is set equal to o¢ = 1. The density is given by the state equation p = po(1 + aC).
(14)
At the present stage the stratification effects in the turbulence model have been neglected in the computations. It is recognized, however, that such effects should be included for cases with high sediment concentration (see e.g. DE VANTIERand NARAVANASWAMY, 1989). 2.2. Boundary conditions The conditions applied at the different types of boundaries are as specified below. Free surface conditions: 0Uh =
ooAv~
~s
(15) _
0¢
OX
Oy'
OS -k- U 0¢ -+- V- -
W = Ot a v OC + WsC = 0 Oz
k = ~-C-~,
P ----Patm
(16) (17)
e = --'KA
(18)
Here r is the von Karman constant, and us. is the surface friction velocity defined by
894
T. UTNES
Us* =X/'rs/P0. The conditions for (k, e) are equilibrium boundary layer conditions which are to be applied at a small distance A below the surface. Bottom conditions: Bottomconditions: A Ouh
v -~z = %
(19) Oh + v
W=-- UOX C = Cb.
(20) (21)
The bottom stress in (19) is given by % = P0fblUb]Ub,where Ub is the near-bottom velocity and fb is the bottom friction coefficient. The bottom boundary conditions for (k, e) correspond to the boundary layer conditions applied at the surface. The condition for C in (21) is given as a local equilibrium bed concentration Cb. Following HAGATUNand EIDSVIK (1986) we assume C b = Cmaxfor S > $ 2 , C b = 0 for S <- S I , and a linear matching for $1 -< S -< S2, where $1 and $2 are given limits of the Shields number at the bottom. By introducing the sediment particle size D, the actual Shields number is defined by s =
I bl
(22)
(Pc - po)gD
Lateral conditions: At the lateral boundaries either ~ or O~/On is specified, while the conditions for u may be of similar form or "free", i.e. computed from the equations with known ~. For k, e and C the lateral conditions are usually specified at other lateral boundaries.
3. NUMERICAL FORMULATION In order to solve equations (1) and (2) a mode splitting procedure is introduced. The total velocity is split in U = U + Ur
(23)
where fi is the vertically averaged velocity and Ur is the residual part. The equations for Ur are formally obtained by substitution of equation (23) into the governing equations (1) and (2). The general solution procedure for each time step is as follows: The vertically averaged equations for (fi, ~) are solved first. These values are substituted into the residual equation system which can then be solved for Ur. The total velocity field is finally obtained from equation (23) and the continuity equation. A finite element formulation is applied to solve these equations. The domain is divided into three-dimensional trilineal elements, with a node configuration that results in perfectly vertical element sides. The projection of these elements in the horizontal plane is used to solve the vertically averaged problem. The numerical formulation is described in more detail elsewhere (UTNES, 1991), and only a short resume is given below. The vertically averaged equations are considered first. A weighted residual formulation in space and time of the momentum equation gives
Sedimenttransportmodelling
895
In
where fa represents all the terms in equation (4) which have not been expressed explicitly, is the space domain and (tn, t,,+l) is the time increment. After integration by parts, and by introducing the following restriction on the weight functions oa + (~. V)Q = 0 ot
in (t,, t,,+l)
(25)
it can be shown that equation (24) reduces to
The choice (26) implies that the weight functions are transported along the streamlines. At time t~+l these functions are chosen equal to the interpolation functions for g, and at t~ they are found by use of equation (26). The time integration can be performed by any suitable method. However, the boundary integral is computed at tn if the normal velocity un is not specified equal to zero. The continuity equation is expressed in ordinary weighted residual form as follows by introducing a two-level time discretization I q~"+Idlq=I
f~
q ~ n d f ~ - I qAt(OVH'fi~+l+(1-O)VH~n)dfl
f~
f~
(27)
where q represents space-dependent weight functions The solution strategy for these equations is as follows: A two-level time discretization and a time-splitting of the momentum equation is applied. The continuity equation is transformed to a Helmholtz equation for the elevation by performing some manipulation (see Ua~Es, 1991), and a mixed interpolation of the variables is used: biquadratic velocity and bilinear elevation. The three-dimensional residual velocity ur is considered next. The generalized weighted residual formulation of the horizontal part of the momentum equation (1) is given by ,n
\ Ot + V . (UUh) + fb df~ dt = 0
(28)
where fb denotes all terms which have not been written out explicitly. Application of the weak form of the characteristics method gives, in analogy with the vertically averaged momentum equation
/n÷l
[df~+l
By introducing the relation (24) and substituting for ~, this equation is used to solve for the horizontal part of the velocity Ur. A two-level time discretization is again applied, and the vertical diffusion is treated implicitly, while all other terms are either known or explicitly treated quantities. The vertical component of the velocity is found by use of the continuity equation. The interpolation of the velocity is biquadratic in the horizontal and linear in the vertical plane. The equations for the scalar variables (k, e, C) are solved in the same way as the
896
T. Ua'NES
ac
Vt ~
ap =
all= aX aC ~x
= - WsC
nst
i
C=C b Zo
Fig. 1. Domain and boundary conditions for computation of sediment concentration profile.
momentum equation. For example, the generalized weighted residual formulation for the sediment concentration equation becomes
+
I7
and the (k, e) equations are treated similarly. 4. COMPUTATIONAL RESULTS In the following two simplified test examples are presented together with a more complex computation of suspended sediment transport in a jet-like flow through a basin.
4.1. Sediment concentration profile This test is defined by a horizontally uniform pressure-driven flow over a flat bottom with specified constant sediment concentration, see Fig. 1. If we assume a parabolic vertical eddy viscosity of the form Av = Kz(1 - z]h)Ub,, where Ub* = ~ / ~ is the bottom friction velocity, the concentration profile is given by the analytical solution
C=
c {Zb(h - z)lW:( u ")
b!z( h _ Zb))
(31)
where z0 is a small distance from the bottom to the first grid point. This solution is shown in Fig. 2 together with model computations for Zb/h = 0.01 and Ws/Ub. = 0.4. In order to evaluate the numerical results, computations have been performed both with the parabolic eddy viscosity and with the (k, e) model. It is seen from Fig. 2 that the numerical result with parabolic viscosity is fairly close to the analytical solution. This is of course what should be expected if the numerical scheme is correct. The result from the (k, e) model also shows a general trend which is similar to the previous solution, but the sediment concentration is now somewhat lower above the bottom, except at the free surface. At this boundary the analytical solution gives zero concentration, while the turbulence model predicts a final value due to the free surface boundary conditions for the (k, e) model.
Sediment transport modelling
897
1.0 ---.......~=~,,~,,,,~. z
0.8
h
0.6
\\~. \X\ x ~ . ,
0.4
-,N
0.2
\\
\
0
I 10 -3
10 -4
L 10"2
I
10 -1
10 0
C/C max Fig. 2. Concentration profile. Zb/h = 0.01, Ws/U, = 0.4. - Analytical solution with parabolic eddy viscosity; - - . - - Numerical solution with parabolic eddy viscosity; - - - Numerical solution with (k, e) model.
Computations with other values of the input variables (bottom concentration, fall velocity etc.) show similar correspondence between the analytical and numerical results. The computations were performed on a grid with 18 points in the vertical direction, with fine resolution close to the bottom.
4.2. Horizontal concentration adjustment This problem is defined by a horizontally uniform flow over a fiat sediment bed, with zero upstream sediment concentration and a specified constant bed concentration downstream, see Fig. 3. Initially there is no suspended sediment, and the problem is to find the steady state concentration adjustment downstream. Again it is possible to derive the analytical solution if the eddy viscosity profile is assumed to be parabolic (cf. VAN RIJN, 1987). The analytical solution is shown in Fig. 4 for the input variables Zb/h = 0.05 and ws/ub* = 0.2, together with corresponding numerical computations have been performed with the parabolic eddy viscosity and with the (k, e) model. A different numerical solution with the actual parameters is also given by VANRIJN (1987) for the case with parabolic eddy viscosity. The computations were performed on a grid with 18 points in the vertical and 31 points in the horizontal direction, with grid stretching in both directions. Our numerical solution with parabolic eddy viscosity is seen to be close to the
a___Ru=
ax
c:o
[
_
I / / L ~ / / / / / / / / / /VT- 7 - / - - / ~ / - - / - - /
onSt
....
C=Cb
~
--
-
-
/ |
7- 7 - 7 - 7- 7- -7 - - / - - [ - - / - - 7 - T - /
rigid bed sediment bed Fig. 3. Domain and boundary conditions for computation of horizontal concentration adjustment.
898
T. UTNES
1.0 C / C m a x = 0 ., 0 5~ ~ . ~
0.8 Z
h 0.6 ~
0.4
•
~ ~-,-- -
-
~ ~.
-~ - ~ -
. . . .
0.2
0 0
~ 10
i 20
I
30
40
x/h
Fig. 4. Horizontal adjustment, zb/h = 0.05, Ws/U, = 0.2. Analytical solution with parabolic eddy viscosity; - - . - - Numerical solution with parabolic eddy viscosity; - - - Numerical solution with (k, e) model.
corresponding analytical solution (Fig. 4). This result also agrees well with VAN RIJN'S (1987) computations for the same case. The model result with the (k, e) turbulence model shows a similar evolution, but now with a general trend of somewhat lower concentration values above the bottom. This tendency with a certain reduction in the concentration values is similar to the previous result (cf. Fig. 2). In summary, these tests show that the numerical model is able to solve such simplified problems fairly accurately with a prescribed eddy viscosity. The inclusion of a twoequation turbulence model modifies these results quantitatively, but the general trend is the same.
4.3 Suspended sand transport in a jet-like basin flow Here we consider a jet-like river flow through a basin, and look for possible suspension and transport of bottom sediments. The geometry is a nearly circular basin with river inand outflow, as shown in Fig. 5. The depth is constant and equal to 10 m, and the uniform in/outflow velocity is 1 m s-1. Along the land boundary the velocity is set equal to zero, and a small horizontal diffusion Ah = 4 m 2 s -1 is included. For problems with significant non-linearities, such as the present one, a certain minimum horizontal diffusion must be included in order to avoid numerical instabilities. The in/outflow conditions for (k, e, C) are zero normal derivatives. The surface and bottom boundary conditions are as given by equations (10) and (11) with zero surface stress, and the sediment is specified by p = 2500 kg m -3, D = 2 x 10 -3 m, Cmax = 0.1, w s = 10 -4 m s -1. The limiting values of the Shields number, $1 and $2, have been chosen to 0.10 and 0.75, respectively. Figures 6 and 7 show the computed surface and near bottom velocity, respectively. It is seen that the jet-like inlet flow is passing through the basin, and recirculation eddies have developed at each side of the geometrical symmetry line. This characteristic jet-and-eddy structure of the flow is mainly a nonlinear effect (UTNES, 1991). The suspended sediment concentration near the bottom is shown in Fig. 8, and the corresponding values along the transversal midsection is given in Fig. 9. As expected, the main part of the suspended sediments are found in the jet-like domain of the current, where the velocity gradients are strong. This is
899
Sediment transport modelling
E
///1
l\\\
l
8 a 4500 Fig. 5.
I
I
- ~000
-500
I///
I
I
0
~;00
I
t000
I rn
4~00
Horizontal grid domain for river-basin flow. T h e nodes are marked on one of the elements.
especially the case in the vicinity of the in- and outflow areas. Generally, the vertical velocity variations are small, and the main advantage of using a three-dimensional model here is to compute more accurately the velocity gradients near the bottom. In many practical situations the addition of wind effects will be of interest, and in such cases the application of a three-dimensional model would be even more important (cf. EIDSVIKand UTr~ES, 1991). 5.
CONCLUSIONS
A three-dimensional hydrodynamic circulation model has been generalized to include turbulence and sediment models in order to compute sediment transport problems. The hydrodynamic part has previously been tested on some basic problems (UTNES, 1991), and the present paper describes some computations with sediment suspension and transport involved. Two simple tests have been performed, both compared with simplified analytical solutions. An example has also been shown for a somewhat more complicated problem, although in that case there are no quantitative results to compare with. However, certain effects may be judged qualitatively, and the result seems to be reasonable. It is recognized that for problems with a relatively high sediment concentration, the sediment and turbulence models should be generalized to include stratification effects. Such effects can be included in the present model.
900
T. UTNES
1 mls I
.I
E
I
~2000
-
1000
Fig. 6.
0
I000
m
2000
Computed velocity field at the free surface.
o o
1 rn/s
o-
Q O
o-
f'--'"~
t l t ~' . . . . .
)
c3 o
o-
\ \ "". " . .' . .1 } ltl liI ~\ ' ' •
",
,
,
s
l
~
. . .... . . / / \
~
.
t
t
.
o o o
#2000
1 -
0
t 000
F i g . 7.
Computed
t 000
m
velocity field near the bottom.
2000
901
Sediment transport modelling
o,I E
. " / I /F/I/
\\\'x~,'~ \
\\\\\\\'\
,."/1/1111 i
,; /71~DII~ lit/if/ ',, (
i _
'
if
/ II1
.lit
k/kXX\ \
I JJ)J)/~, '
\\
I
tilt\ \ \ I III I/I,,"
I
{000
,
I1111}/,,
~-\
--
t
°ltt\l\t
o
~tSO0
t tt
I IIIlll
",
Fig. 8.
III\\\\",,
U///I/
,,"
_
///
I
-SO0
I
0
I
SO(}
I
I000
ISO0
m
Computed suspended sediment concentration in the near-bottom domain. The iso-values are: (0.0025, 0.005, 0.01, 0.02, 0.03, 0.04, 0.05).
E
9'.0025
?-
O o
[
- 4000
Fig. 9.
I
-,gO0
I
I
0
,gO0
I
4000
I
03
4S00
Concentration isolines along the transversal m i d section.
902
T. Uan~ES
Acknowledgements--This work was partially supported by Konsesjonsavgiftsfondet (Fund of License Fees).
REFERENCES DE VANTIER B. A. and R. NARAYANASWAMY(1989) A suspended sediment flow model for high solids concentration using higher order turbulence closure. Advances in Water Resources, 12, 46-52. EIDSVIK K. J. and T. UTNES (1991) A model for sediment reentrainment and transport in shallow basin flows. Coastal Engineering, 15, 247-256. HA~ATUN K. and K. J. EIDSVIK(1986) Oscillating turbulence boundary layer with suspended sediments. Journal of Geophysical Research, 91, 13045-13055. MELLOR G. L. and A. F. BLUMnERG (1985) Modelling vertical and horizontal diffusivities with the sigma coordinate system. Monthly Weather Review, 113, 1379-1383. O'CONNOR B. A. and J. NICHOLSON(1988) A three-dimensional model of suspended particle sediment transport. Coastal Engineering, 12, 157-174. RODI W. (1980) Turbulence models and their application in hydraulics. IAHR State-of-the-Art Paper, Delft. SHENG Y. P. and H. L. BUTLER(1982)Modelling coastal currents and sediment transport. Proceedings of the 18th Coastal Engineering Conference, Cape Town, 14-19 November, pp. 1127-1148. TEISSON C. and D. FRISCH (1988) Numerical modelling of suspended sediment transport in the Loire estuary. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May-1 June, pp. 14-22. UTNES T. (1991) Finite element modelling of quasi-three-dimensional nearly horizontal flow. International Journal of Numerical Methods in Fluids, 12, 559-576. VAN RUN L. C. (1987) Mathematical modelling of morphological processes in the case of suspended sediment transport. Dr Thesis, Technical University Delft, Delft, 135 pp. and 2 appendices. VAN RIJN L. C. and K. MEIJER(1988) Three-dimensional mathematical modelling of suspended sediment transport in currents and waves. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May-1 June, pp. 89-99. VEERAMACHANENIR. and E. J. HAYTER(1988) Mathematical Modelling of sediment transport at tidal inlets. IAHR Symposium on Mathematical Modelling of Sediment Transport in the Coastal Zone, Copenhagen, 30 May1 June, pp. 14-22. WANG S. S. Y. (1989) Preprints of International Symposium on Sediment Transport Modeling. New Orleans, LA, 14-18 August, 829 pp.