(1)
Junqiang Xia1 and
Guangqian Wang2
1 Doctor Student, Key Laboratory for Water and Sediment Sciences of Education Ministry, Tsinghua University, Beijing, 100084,China
2 Prof. Key Laboratory For Water and Sediment Sciences of Education Ministry, Tsinghua University, Beijing, 10084, China
Abstract:
The method of fractional steps along the spatial direction is applied to split
vertical plane 2-D suspended sediment transport equation into two 1-D
advection-diffusion equations. Then two different computation schemes are
proposed to solve them. Through comparison and analysis, this paper suggests
that: as the advection process is weak in transport, the computation scheme (I)
is proper through using a Crank-Nicholson (C-N) implicit difference scheme both
in the x and y directions; as the advection process is preponderant in
transport, the computation scheme (II) is very reasonable through using an
exponential explicit difference scheme in the x direction and the C-N scheme in
the y direction, respectively. Experimental results are used to verify two kinds
of computation schemes, and good agreements between the calculated and measured
results are obtained.
Keywords: vertical plane 2-D, suspended sediment transport, fractional steps method, advection-diffusion equation
In natural rivers, the movement of suspended sediment affects such practical problems as aggradation or degradation of channel beds, diversion works, and flushing sediment projects. In order to settle these engineering problems, it is very important to know how the suspended sediment moves in flow. In some regions with deep water depths and wide water surface widths, we usually apply vertical plane 2-D transport equation to describe the movement process of suspended sediment.
This transport equation is usually solved by the alternating direction implicit method. The fractional steps method has been used widely in practical engineering problems recently owning to its flexible forms and high numerical stability. In this paper, the fractional steps method along the spatial direction is used to split vertical plane 2-D suspended sediment transport equation into two 1-D advection-diffusion equations. Then two different kinds of computation schemes are employed to solve these equations. Through comparison and analysis, this paper suggests as following that: as the advection process is weak in transport, the computation scheme (I) is proper through using a C-N implicit difference scheme both in the x and y directions; as the advection process is preponderant in transport, the computation scheme (II) is very reasonable through using an exponential explicit difference scheme in the x direction and a C-N scheme in the y direction respectively. Experimental results are used to verify two kinds of computation schemes, and good agreements between the calculated and measured results are obtained.
The basic equation for vertical 2-D suspended sediment transport is described (Zhang, 1993) as:
(1)
where S= suspended sediment
concentration; u= local flow velocity in the longitudinal (x) direction; v=local
flow velocity in the vertical (y) direction;
=sediment diffusion coefficient in the x direction;
= sediment diffusion coefficient in the y direction;
=the settling velocity of suspended sediment particle; and t=time. Through
rearrangement Eq. (1) can be also expressed as
(2)
where
and
can be defined as
，
(3)
Because of complicated initial and
boundary conditions, it is very difficult in finding the analytical solution of
Eq.(2), and numerical solution methods should be employed for computing the
suspended sediment concentration. The influence of sediment diffusion
coefficients on the calculated results is rather obvious. With
regard to how to calculate the sediment diffusion coefficient
in the vertical direction, we usually assume that
is a constant along the y direction, which is equal to 0.067
. Many researchers proposed their views on how to determine
(Coleman,1970;Van Rijn,1988; Motohiko,1992). Based on the results of some
experiments, we concluded that the value of
should be related to not only the water depth y, but also the turbulence
strength of flow. In this study, the vertical sediment diffusion coefficient
proposed by Tan (1998) is given as:
(4)
in which
=correction factor, which presents the diffusion difference between sediment
particles and flow mass point;
=shear velocity; h=water depth;
=relative depth; y=the distance from the bed.
=parameters which are 3.0 and 0.8, respectively. With regard to the longitudinal
diffusion coefficient
, according to the magnitude order analysis accomplished by Graf (1997), the
following relationship is used:
=50.0
(5)
We suppose that the suspended sediment concentration in the whole computation domain is zero in the initial state and that there is a steady concentration resource of suspended sediment at the inlet boundary subsequently. The vertical distribution form of suspended sediment at the inlet suggested by Van Rijn (1988) is given:
(6)
where the Rouse suspension number Z can be defined as:
(7)
in which k is
Karman constant with the value of 0.4;
is the near-bed sediment concentration; and a= 0.05h is the reference level from
the bed.
We usually assume that at the water surface, there is an equilibrium state between the upward turbulence diffusing of suspended sediment and the downward falling of sediment particles. Namely, no sediment particles penetrate the water surface, so the boundary condition is defined as:
at
y=h
(8)
It
is very complicated to determine the near-bed condition for suspended sediment
transport under the non-equilibrium condition. Based on some experimental data
and numerical results computed from different bed conditions, Zhou (1997)
proposed that it is reasoned that the Dirichlet boundary condition is correct
for the process of erosion; while the Neumann boundary condition holds for
deposition. However, for general bed aggradation or degradation problems, the
Neumann condition used as the near-bed boundary condition is expressed:
at
y=0.05h
(9)
where
is the near-bed concentration in the equilibrium state. This treatment method on
the near-bed boundary is based on such the assumption that the suspended
sediment near the bed can be supply from the bed material simultaneously, and
the near-bed sediment concentration is equal to the local saturated
concentration.
The outlet cross section should be located in the position, which is enough away from the inlet cross section. On this assumption, we hypothesize that its boundary condition is defined as
(10)
Generally speaking, to get an analytical solution, some additional assumptions should be made, for example, the longitudinal diffusion, compared with the vertical diffusion is negligible, and so on. Only for a few simplified cases can the analytical solution be obtained (Zhang,1980; and Wei,1982). For many complicated cases, we have to fall back on numerical simulation methods to solve Eq.(2).
The fractional steps method is such a way that split multi-dimensional equations into some simple 1-D equations, and then use all kinds of effective computation schemes to solve them. In this study, the fractional steps method along the spatial direction is applied to split Eq.(2) into two 1-D advection-diffusion equations:
The fractional steps equations Eq.(11)
and its corresponding difference equations are equivalent to the original
equations Eq.(2) in terms of mathematics. Detailed principles of the fractional
steps method can be referred to Cao(1994). According to characteristics of
fractional steps equations, we can use all kinds of effective numerical
computation schemes to solve them. Usually the whole time step is divided into
two half steps. In the first half time step, Scanning the x direction using Eq.(11-1),
the intermediate value
can be obtained; in the second half one, Scanning the y direction using Eq.(11-2),
the value
at the time level (n+1) can be obtained.
This paper suggests that for the
vertical plane 2-D suspended sediment transport equation, the advection movement
is preponderant in the longitudinal transport of suspended sediment; and the
diffusion movement is dominant in the vertical transport. As the advection term
(
and
) is very week, the Crank-Nicholson implicit difference scheme is used to solve
Eq.(11-1) and Eq. (11-2) respectively and obtain reasonable results. As the
advection term is very strong, especially
, we should avoid applying this scheme, because an unphysical oscillation will
occurs if we use C-N scheme. Under the circumstances, so we employ the
combination of explicit and implicit schemes. Owning to the high numerical
stability for the exponential scheme proposed by Lu (1988), which can apply to
the situation with dominant advection action, this scheme will be used in the
longitudinal direction to eliminate such spurious oscillations. And in the
vertical direction, because the value
is relatively small, the C-N implicit scheme is still used to compute the Eq.(11-2)
in order to take larger time increment.
As the advection term is very weak, we can use C-N scheme to compute the x and y direction fractional steps equations. In the first half time step, we can calculate the central finite-difference approximation for the space and time derivatives, then by substituting these approximations into Eq.(11-1), the corresponding difference equation in the longitudinal direction can be yield:
(12-1)
Rearrangement of Eq.(12-1) can obtain the following form:
(12-2)
where
are coefficients which relate to
the values in the known time level. Substituting the corresponding
discretization scheme of initial and boundary conditions Eq.(6) and (10), during
the course of solving the Eq.(12-2),
each step in the algorithm leads to a system of equations that have tridiagonal
matrix which are comprised of coefficients
. In this case, we can use the tridiagonal matrix algorithm (TDMA) technique to
solve these equations and yield the intermediate value
.
In the second half time step, using the same treatment as the first half time step, the finite difference scheme in the vertical direction can be yield:
(13-1)
rearrangement of Eq.(13-2) can also obtain the following form:
(13-2)
where
are coefficients which relate to
the intermediate value. Similarly, substituting the corresponding discretization
of boundary conditions Eq.(8) and (9), and adopting the same solution method as
Eq.(12-2), we can obtain the value
in the next time level (n+1).
Obviously, the difference schemes Eq (12-1) and (13-1) have second-order accurate both in space and in time. Through using the analysis method of freezing up the coefficients, we can rigorously prove that the computation scheme (I) is unconditional stable.
As the advection term is very strong, we can use the C-N scheme to discretize Eq.(11-1) and use the exponential scheme to discretize Eq.(11-2). In the first half time step, the difference scheme in the x direction can be expressed as:
(14)
in which
is the fixing factor in the longitudinal direction. This difference equation is
explicit, which has first-order time accurate and second-order space accurate.
As the coefficient of advection term is larger than the longitudinal sediment
diffusion coefficient, namely
, using the exponential scheme can produce very high numerical stability and can
eliminate numerical oscillations. In the vertical direction, the difference
equation is the same as Eq.(13-1). We can use the C-N scheme in the y direction
mainly due to weak advection movement in the vertical direction, and due to
taking larger time increment.
Using the analysis method of “freezing up” the coefficients, the computation scheme (II) is conditionally stable if
(15)
During the course of computations, we
always hope to take larger time increment, but selecting of time increment must
be restricted to stability condition and computation accurate. Although the
computation scheme (I) is unconditional stable, to satisfy the domination of
diagonal elements in the coefficient matrix as we use the solver of TDMA, the
time increment has to be limited in order to meet the requirement. The
computation scheme (II) is conditionally stable one, and the time increment is
approximately determined owning to the variance of coefficients. Through the
preliminary analysis, the time increments for the computation scheme (I) and
(II) are
respectively, and they must meet
the following conditions in order to satisfy the numerical stability
requirement.
;
(16)
where
is the maximum value of
.
At present, experiments data for the
transport of suspended sediment in the unsteady state is insufficient and
unavailable. But we can obtain the data of depth-averaged concentrations in the
steady state. So in this paper, during the course of computing, assuming that
the bed topography is unchanged, and that the computation time consumed is
enough long, and that the value of
approach to zero, we can verify the
validity of these schemes using the data in steady state.
We use the foresaid two computation schemes to simulate the transport process of suspended sediment in a settling basin. This data provided by Wei (1982) are as following: the water depth is 1.53m; the depth-averaged velocity is 0.12m/s; the settling velocity of suspended sediment particle is 0.0176cm/s; the Manning coefficient of this basin is 0.014; the depth-averaged sediment concentration is 3.012 kg/m3 in the inlet boundary; and the Rouse suspension number is 0.09. At the same time, we also assume that the vertical distribution of flow velocity in the settling basin is subject to logarithmic velocity formula. The computation domain comprises 800 by 29 nodes. In order to analyze and compare these two computation schemes under the same conditions, identical time and space increment are applied to these schemes with △x=4.0m , △y=0.05m and △t=12.0 s.
Fig. 1 Variation process of concentration with time from scheme (I) and (II)
Fig.1 shows the variation process of suspended concentration with time at the downstream distance x=1000.0m. As can be seen from this figure, we find that computation results of scheme (I) and scheme (II) are very approximate. The reason for obtaining similar computation results is that the advection term is very weak, and the results of scheme (I) may not produce an unphysical oscillation phenomenon.
Fig. 2 Comprison between calculated and observed results
As the value of the time derivative of
the suspended concentration
at all the nodes approaches to zero
during the computation, we can deem that the computation results have already
reached the solution in the steady state, and then we can use the experiment
data in the steady state to verify the numerical simulation results. It can be
seen from Fig.2 that the results from these two computation schemes are both
very approximate to the observed depth-averaged suspended sediment
concentration. It is evident that very low longitudinal velocity in the settling
basin and very weak advection process leads to non-oscillatory computation
results of the scheme (I).
Fig. 3 Vertical distribution of suspended sediment concentration
Furthermore, this paper provides the numerical results of vertical distribution of concentration for the scheme (II), as shown in Fig.3. At the inlet boundary, the vertical distribution of suspended sediment concentration is very non-uniform, which the concentration near the water surface is relatively low, and the one near the bed is relatively high. Through the advection-diffusion process for a long distance, at the downstream location x=1000m, the vertical distribution of suspended sediment tends to be uniform, which approaches to the distribution form in the equilibrium state.
In order to verify the validity and high numerical stability of the computation scheme (II) as the advection term is very strong, we suppose another numerical example. Realistic hypothetical data are used in computations because real data of sufficient detail are unavailable. All channel and flow data are: the manual channel length is 50km; the water depth is 2.0m; the depth-averaged velocity is 2.0m/s; the Manning coefficient is 0.01. We also assume that the suspended sediment concentration is zero in the initial state, then a steady concentration source is input at the inlet boundary, with the median diameter of suspended sediment of 0.05mm,and the corresponding settling velocity of 0.156 cm/s. Similarly, we also assume that during the computation, the bed topography is unchanged. Moreover, the time increment and the space step in the schemes (I) and (II) are identical in order to analyze and The suspended sediment concentrations for the downstream cross sections gradually increase and approach to some steady value in the end through the processes of advection-diffusion. Fig.4 shows the numerical results of the schemes (I) and (II) at the downstream location x=2000m. Seen from the figure, we know that the results of the scheme (I) produce an unphysical oscillation phenomenon, and that of the scheme (II) can not take on oscillations and maintain the numerical stability, and corresponds to actual transport laws of suspended sediment. The reason is that as the advection term is very strong, the scheme (II) in the longitudinal adopts the exponential scheme with relatively high stability, and eliminates oscillations in the computation results. The scheme (I) still adopts the C-N implicit scheme in the x direction as the advection term is strong, which lead to appearance of oscillatory results.
Fig. 4 Variation process of concentration with time from scheme (I) and (II)
In terms of the time increment, the available maximum time increment for the scheme (I) can reach twice time than that for the scheme (II); in terms of the calculating speed, the differences between consumed computation times of these two computation schemes are very little because the scheme (II) employs explicit scheme in the x direction, and need not solve the triangular equations, which can save a lot of CPU time.
In terms of the computation accuracy, the scheme (I) can reach second-order accurate both in time and in space and can take larger time step with good stability as the advection term is very weak, so in this case, we should employ the scheme (I). As the advection term is very strong, in order to avoid the appearance of unphysical oscillations, we should apply the computation scheme (II), which can eliminate oscillations and take larger time increment at the same time
In this paper, based on applying the fractional
steps methods to the vertical 2-D transport equation for suspended sediment, two
kinds of different computation schemes are proposed. Moreover advantages and
disadvantaged for these schemes are analyzed. Through analysis, this paper
suggests that: as the advection term is weak, it is reasonable to apply the
computation scheme (I) to solve Eq.(2); as the advection term is strong, it is
proper to apply the computation scheme (II) to solve Eq.(2). In the end,
experimental results are used to verify two kinds of the computation schemes,
and good agreements between the calculated and measured results are obtained.
Acknowledgements
This project was supported by the National Natural Science Foundation of China and the Ministry of Water Resource (Grant, No.59890200).
References
Cao zhude and Wang yunhong(ed), Hydrodynamics simulation of sediment, Tiangjing University Press, 1994, P122-136 (in Chinese).
Coleman N.L.,Flume studies of the sediment transfer coefficient, Water resource research, Vol.6,No.3, 1970, P801-809.
Graf W.H. and Altinakar M.S.(ed), Hydraulics of fluviology, translated by Zhao wenqian and Wan zhaohui, Chengdu University of Science and Technonlogy Press, Chengdu, P417-474 (in Chinese).
Motohiko, Vertical distribution of suspended sediment in uniform open-channel flow, Journal of Hydraulic Engineering , ASCE,Vol.118,No.6, 1992,P936-941.
Lu jinfu and Gu lizhen (ed), Difference method for partial differential equations, High Education Press,1988, P200-215 (in Chinese).
Tan guangming, Xia junqiang and Zheng bangmin, Study on the turbulence diffusion coefficient for the sediment-laden flow, Journal of Hydraulic, No.12, 1998, P23-27 (in Chinese).
Van Rijn L.C., Mathematical modeling of suspended sediment in non-uniform flows, Journal of Hydraulic Engineering ,ASCE, Vol.112, No.6,1986,P443-455.
Wei zhilin, Settling process of sediment in two-dimensional steady uniform flow, Journal of Wuhan University of Hydraulic and Electric Engineering, No.4, 1982, P35-47 (in Chinese).
Zhang qisun, Diffusion process of sedimentation in open channel flow and its application, Journal of sediment research, 1980,P37-51 (in Chinese).
Zhou jianjun, bed conditions of non-equilibrium transport of suspended sediment, International Journal of sediment research, Vol.12, No.3, 1997, P 241-247.
Zhang, ruijin and Xie jianheng (ed), sediment research in China (systematic selections),China water and power press,Beijing,1993, P43-66.