(1)
Li Bei 1
and Tang
Shifang 2
1Tianjin Research Institute of Water Transport Engineering, Tanggu,Tianjin, 300456 China
Tel:+86 22 25707168-382, Fax:+86 22 25795125, E-mail:li3-pei2@sina.com
2Shanghai Municipal Government International Shipping Center Construction
Administration,
Port Planning Department, 220 Si Chuan Road(Central), Shanghai, 200002 China
Tel:+86 21 63292992, Fax:+86 21 63292992, E-mail:sifant @citiz.net
Abstract: In this paper a 3-dimensional numerical simulation of the tidal flow, salinity and sediment in the dredging channel in an estuary was set up by the method of coordinate transformation, thus solving the treatment of the complicated boundary. This approach had been used to study the regulation work on the Lingdingyang channel in the Pearl River Estuary and obtained satisfactory results. The forecast silt content and salinity distribution are in good agreement with the field measured results. The calculated back silting volume after dredging in the first phase channel work by the simulation is very close to the on-site measured one.
Keywords: 3-dimensional numerical simulation, salinity intrusion, tidal flow,
salinity, sediment
An estuary is a river reach where the river flows into the sea; the flowing downward river runoff and the up-bounding salt water meet here. After dredging a channel in the river mouth bar, the salt water in the channel will go up, meanwhile the distribution of velocity and silt content in the longitudinal and vertical directions are different from those of the upstream river reach. The back silting rate of the channel in the river mouth bar is closely interrelated to the salt-water intrusion.
The basic process of the salt-water intrusion in the estuary is a 3-dimensional unsteady physical phenomenon. The plane 2-dimensional numerical model can t satisfactorily describe the stratification of the salinity and the vertical circulating flow phenomenon in the estuary after dredging the channel. Therefore, to develop the 3-dimensional numerical simulation of tidal flow, salinity and sediment in the dredging channel in an estuary is of important scientific significance and practical value.
The set of basic equations for calculating the 3-dimensional flow, salinity and sediment in an estuary was formed by N-S equation, convection and diffusion equation of salinity and silt, state equation and river bed deformation equation.
The method of the coordinate transformation was developed to lengthen and shorten the X, Y, Z directional scales in this paper. The vertical conversion can conveniently simulate the sea area with uneven bottom topography. The horizontal conversion can accurately simulate the shape and location of the building, and the tortuous trend of the channel. This method can solve the problem of precision when dealing with a complicated boundary that other calculation methods can t does.
The equation of continuity
(1)
The momentum equation
(2)
(3)
(4)
The state equation
(5)
Convection and diffusion equation of salinity
(6)
Convection and diffusion equation of sediment
(7)
(8)
in which, x,y,z are respectively longitudinal, transverse and vertical coordinates; u, v, w are respectively the components of velocity in the direction x, y, z; g is the acceleration of gravity; Nx, Ny, Nz are respectively the turbulent viscosity coefficients in the direction x, y, z; Ax, Ay, Az are respectively the turbulent diffusion coefficients of salinity in the direction x, y, z; Dx, Dy, Dz are respectively the turbulent diffusion coefficients of sediment in the direction x, y, z; po is the density of fresh water; Ws is the setting velocity of silt in still water; s is the salinity; c is the silt content; t is time; f is Corialis parameter, f=2Ωsinθ, Ωis geostrophic angular velocity, θ is latitude; γc is the dry unit weight of silt; ΔH is the increment of river bed deformation; qx, qy are respectively the rate of transportation of bed load per unit width in the direction x, y; P is the water pressure
(9)
Where Po is
atmospheric pressure.
Subscript b in the equations indicates that the location of the value is near bottom.
The equations (1) to (9) are the basic equations for calculating the 3-dimensional flow, salinity and silt deposition.
The condition of determinate solution of 3-dimensional mathematical model are of two kinds: the boundary condition and the initial condition.
whereτsx, τsy are respectively the surface wind stress components. When the wind stress is not considered, then τsx=τsy =0.
is the net flux of surface salinity, when the evaporation of salt from the water
surface is not considered, then
=0.
Wf is the effective setting velocity of silt, Wf=W-Ws
(1) Water flow
u=v=w ( non-sliding condition)
(sliding condition)
(sliding condition)
where τbx, τby are respectively the components of bed shear in the direction x,y
(2) Salinity
(3) Silt
The exchange mechanism of the silt near the river bed is very complicated and there are many ways to describe the boundary condition,the following was used in this paper.
whereτb is the bed shear stress; τd, τe are respectively the critical setting shear and critical erosion shear; E is the erosion constant.
This kind of boundary condition of silt is widely used in river and coastal engineering.
(1) Water flow
There are two kinds of boundary condition of water flow at the open boundary. One is to give the water level process at the boundary; the other is to give the flow velocity process there.
(2) Salinity
There are also two kinds of boundary condition of salinity at the open boundary. One is to give the salinity process at the boundary, the other is to give the salinity process when flowing in, while to take the flux as zero when flowing out.
(3) Silt
The boundary condition of silt at the open boundary is similar to that of salinity.
There are also two kinds of lateral boundary condition: open and close boundaries.
There are two kinds of initial condition, one is to take the same value for a physical quantity in the whole calculation domain, and other is to give the instantaneous space distribution of every physical quantity.
The main idea of coordinate transformation in this paper is to lengthen and shorten the three coordinate directional scales of the calculation domain, turn the irregular 3-dimensional calculation domain into a unit cube. The calculation grid was set up in this cube and the basic equation set in original coordinate system was transferred to the new coordinate system to be carried out differential discreteness and calculation.
Transfer the original coordinate system (x, y,
z, t) to new coordinate system
, the relation between the two systems is defined as:
where L is the characteristic length, B is the characteristic width.
The new velocity component under the new
coordinate system
may be defined as follows:
The relation of the velocity components under two coordinate systems is:
When the surface wind stress and salinity flux are not considered
Non-sliding condition
Sliding condition
Where
Boundary condition at river bed for salinity and silt
Or
The silt boundary condition at the open boundary is similar to that of salinity.
For lateral open boundary
For lateral solid boundary
2.5 The setting up of calculation format
The discreteness of the basic equations is carried out in the transformed unit cube.
In driving the differential equation, uneven grid was used in all of the 3 directionsζ，λ，η. The grid space of horizontal directions(χ，λ)may be taken arbitrarily. If the change of topography is rather serious or it is near the engineering project, the space should be smaller, if it is far away from engineering project the space can be larger. It should be avoided changing space of adjacent grid suddenly. Gradually transit it so as not to lead to the unstability of calculation. If the water area to be calculated is relatively regular, the regular quadrilateral grid can be used.
In dispersing momentum equation the forward difference is used for time derivative, the explicit central difference is used for the space derivative inχ，λdirections,while the implicit central difference is used for the space derivative inηdirection.
As to the convection and diffusion equation of salinity and suspended sediment, the treatment of convection item is a difficult problem. An ideal difference format should have positive definition, accuracy and less numerical dissipation. In actual computation, for dispersing convection and diffusion equation, the explicit central difference is used for the horizontal convection item, the implicit central difference is used for the vertical convection item, the forward difference is used for time derivative, the explicit central difference is used for the horizontal diffusion item, the implicit central difference is used for the vertical diffusion item.
In dispersing the items in a basic equation which have an important influence on the dependent variable, in order to make the diagonal matrix of the tridiagonal matrix become dominant, a special treatment has been done to some items. That is, to take the average value at nΔt and (n+1)Δt as the calculation value.
Although the basic equation set N-S equation that describes the motion of fluid is still suitable for describing the laws of the average flow of turbulence flow, it is necessary to introduce an assumption into the turbulence item in it in order to close the equation set because the equation group is not closed when the N-S equation set studied and described the average quantity. In this way the turbulence item is related to the average quantity, i.e. to establish a turbulence model.
Because the variation range of the scale, which reflects the turbulence motion, is very great, the smallest scale of the turbulence motion is much smaller than that of the flow domain. Therefore, as to calculating in engineering application, the features of large and small vortex’s motion is not differentiated in turbulence model. In fact, the large vortex that absorbs energy from the main current is anisotropy and the small vortex consuming energy is isotropic. Because the structure of vortex is different, the available turbulence model has their own range of validity.
When solving the average N-S equation, the commonly used turbulence models are equation-0, equation-1 and equation-2, etc.
Between them, the K-εmodel of 2 equation has been most widely used in the past ten or more years. Many turbulence problems encountered in engineering have been solved quite well by using K-εmodel. But in the K-εmodel the turbulent viscous coefficient is assumed to be isotropic, which makes its application in complicated flow be limited.
Although the 0-equation model is the simplest turbulence model, it has satisfied the needs of precision of numerical simulation of tidal flow, salt and sediment Therefore the 0 equation of type Munk was used in this research. Its turbulent viscous coefficient is determined by the following equations:
where
where, k is the Karman
constant, u* is friction velocity, H is water depth, a1 is a constant that can
take a very small value in order to solve the diffusion problem produced when
calculating to the river bed, Ri is Richard parameter,
are coefficients.
In the case of fine grain silt, the turbulent diffusion coefficient of the silt is considered generally to correspond to the turbulent viscous coefficient of the water flow.
The research of regulation technology of Lingdingyang channel in the Pearl River Estuary is a key technical project of China s 8th Five-year Plan. In the near future the channel is planned dredging to depth of 11.5m (calculating from theoretical datum plan, similarly hereinafter). At a special future date it is planned to reach 12.5m. In order to provide accurate flow field, salinity field and silt content field for the study of the regulation of channel, a 3 dimensional numerical simulation of tidal flow, salinity and sediment was carried out using the mode provided in this paper.
Based on the good result of verifying, the 3-dimensional simulation calculation of two channel schemes was carried out. Scheme 1 is to dredge the channel to the depth of 11.5m, side slop 1:10, bottom width 160m. Scheme 2 is to dredge channel to the depth of 12.5m, side slop 1:10, bottom width 160m.
The following is the simulation results.
(1) Velocity After the channel has been deepened, water flows back to the channel, the flow direction deflects to the channel about 1°, the mean velocity on a vertical in the channel increases. For scheme 1, the velocity of flood tide current increases 1.2% to 4.8%, the velocity of ebb tide current increases 2.2% to 6.7%; For scheme 2, the velocity of flood tide current increases 1.1% to 3.2%, the velocity of ebb tide current increases 1.5% to 6.7%. To view the velocity variation on the vertical, the surface current velocity increases more, the bottom current velocity decreases.
(2) Salinity The mixing pattern of salt and fresh water is strong mixing pattern in the dry season, and the mild-mixing pattern in flood season. After dredging the channel, the seawater of high salinity up-bounds along the deepened channel. Due to the dredging depth is not too deep, the length of intrusion is limited. During the dry season the length of intru sion of the sea water with salinity of 2% moving upstream along the bottom is 960m for scheme 1, and 1280m for scheme 2. During the flood season this length of the seawater with salinity of 1% is 640m for scheme 1, and 960m for scheme 2. Fig. 1 is the sketch map of salinity distribution at the slack in tide flood during flood season.
(3) Suspended sediments During the dry season, the incoming marine sediments are limited. The sediment content in the water is very little, about 0.03kg/m3 and the variation on
Fig.1 The isograph of salinity (&) at the slack in tide flood during flood season
vertical is not large. During flood season, the sediment content in the water is rather large, about 0.13kg/m3 to 0.52kg/m3. It can be found from calculation that the time when the silt on the bottom was suspended and the sediment content reach the maximum is not at the moment when the flood or ebb current velocity is highest, but about 3 hours later after that, corresponding to the slack in tide flood or ebb. Moreover, the sediments content at the bottom layer varies rather obviously with tide flood or ebb. After dredging the channel, the sediment content does not vary more obviously than before, because the depth increase is limited and the situation of the incoming sediment is not changed. Fig. 2 is the sketch map of sediment content distribution at the slack in tide ebb during the dry season.
Fig.2 The distribution of suspended silt content along the longitudinal section at the slack in
Fig.3 The distribution curve of the deposition rate along the longitudinal direction of channel
(4) Sediment deposition in the channel The sediment deposition distribution in the channel in dry season, flood season and whole year was obtained and showing in fig. 3. From the result of numerical simulation, it is found that the amount of deposition in the dry season accounts for 35% of that of the whole year, in the flood season accounts for 65%. The average deposition rate of scheme 1 is 0.479m/a; its deposition amount is 3.24 million m3. The average deposition rate of scheme 2 is 0.5m/a, its deposition amount is 3.43 million m3. These numerical values are identical to the experiment results by physical model and calculation results by empirical equation. Moreover they are quite near the field measured back silting quantity after the engineering project was implemented.
The basic equation
set for calculating the 3-dimensional water flow, salinity and silt deposition
is composed of N-S equation set, convection and diffusion equation of salinity
and silt, state equation and river bed deformation equation.
In order to use the method of finite difference more conveniently the technique of coordinate transformation was developed to lengthen and shorten the x, y, z directional scale of the calculation domain. By this way, the vertical conversion can conveniently simulate the sea area with uneven bottom topography. The horizontal conversion can accurately simulate the shape and location of the building, and the tortuous trend of the channel.
Based on the relatively satisfactory results of the verifying calculation, the simulation calculation of two engineering schemes of dredging channel to 11.5m and 12.5m water depth was carried out. From the simulation results it can be found that the simulation format developed by this paper is feasible.
References
[1] Li Bei et al, Two-dimensional mathematical model of wave, tidal flow and silt deposition, Waterway and Port, 1993. (In Chinese).
[2] Li Bei et al, The two dimensional numerical model of wind-driven current and tidal flow, published by Tianjin Transport Engineering Research Institute, Ministry of Communication, 1994. (In Chinese).
[3] Liang Zaichao, Mechanics of Turbulence, Henan Publication House of Science and Technology, 1988. (In Chinese).