Sigma-coordinate numerical model for side-discharge into natural rivers

333

2009,21(3):333-340 DOI: 10.1016/S1001-6058(08)60153-5

SIGMA-COORDINATE NUMERICAL MODEL FOR SIDE-DISCHARGE INTO NATURAL RIVERS* LIU Zhao-wei, CHEN Yong-can, LI Ling State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China, E-mail: [email protected] ZHENG Jing-yun Wenzhou Oujiangkou Development and Construction Headquarters, Wenzhou 325027, China

(Received June 16, 2007, Revised October 22, 2007)

Abstract: Due to large topography slopes in natural rivers, pollutant concentration embodies a property of three-dimensional distribution when wastewater is discharged from effluents along the bank. With the sigma coordinate along the vertical dimension fitted to both the moving free surface and the bed topography, a three-dimensional numerical model was developed in the present work to address pollutant transport processes in the above-mentioned cases. To avoid the reduction in accuracy caused by spurious diffusion in the case of steep bottom slopes, a formula for horizontal diffusion in the sigma coordinate system was derived. A case study for the side discharge into a straight open-channel flow shows that numerical results are verified well by experimental data. Furthermore, the present model is also verified by the simulation of discharging wastewater from Fuling Phosphorus Factory effluent into the Three Gorges Reservoir and the agreement between the numerical simulation results and field observation data is satisfactory. The change of the mixing zone scope in the water surface versus the layers along the vertical dimension was also discussed in detail. The study shows that a more realistic calculation for pollutant discharge has been provided by the present model than by the depth-average model which predicts an unrealistically smaller mixing zone. Key words: sigma-coordinate, three-dimensional model, side-discharge, natural rivers

1. Introduction  In the prediction of mixing zones in a natural river, two-dimensional models have been commonly used with a simpler structure and higher computational efficiency. For example, with a depth-averaged finite element model, Jiang[1] simulated the concentration distribution in the region near the Taohuaxi effluent which discharges wastewater into the Three Gorges Reservoir, and Qi[2] investigated pollutant transport processes discharged from a effluent by means of hierarchy grids in a  * Project supported by the National Basic Research and Development Program of China (973 Program, Grant No. 2006CB403304). Biography: LIU Zhao-wei (1973), Male, Ph. D., Lecturer

two-dimensional domain. However, mixing zone predicted with the two-dimensional model is generally smaller than the actual size, because the pollutant concentration shows a three-dimensional distribution in the region adjacent to an effluent due to large topography slopes in the natural river. Consequently, three-dimensional models became the focus of some researchers. With the horizontal multi-layer three-dimensional model, Chen[3] predicted the mixing zones resulted by the Fuling Phosphate Fertilizer Factory effluent discharging wastewater into the Yangtze River. However, the water column in the shallow region near the effluent was represented by only one layer and thus the three-dimensional feature of concentration could not be predicted accurately. The sigma-coordinate model is another major form to solve the shallow water equations[4-10]. In this

334

approach, there could be more accurate numerical calculation along the vertical dimension in the region near the effluent because the sigma-coordinate system is a topography-fitted, appropriate for irregular bed bottoms. However, when bottom slopes are steep, encountered often in the natural river, the conventional definition of diffusion in the sigma transformation allows communication between the water in shallow areas and the water in adjacent deep regions, and induces spurious diffusion, which may result in considerable errors over a long integration period. In the following section, a three-dimensional numerical model with the sigma-coordinate system is presented to provide a realistic prediction of mixing zone formed by discharging from effluents into a natural river. To address the steep bed bottoms which are common in natural rivers, a formula for sigma-transformation is included to eliminate the spurious diffusion. A case study of the Fuling Phosphorus Factory Effluent discharging into the Three Gorges Reservoir is performed to validate the model and investigate the three-dimensional distribution of concentration in the regions near the effluent. 2. Mathematic model and numerical schemes The Navier-Stokes equations are a general model which can be used to model water flows in many applications. However, when considering a problem in which the horizontal scale is much larger than the vertical one, the three-dimensional shallow-water equations will suffice. 2.1 Governing equations in cartesian coordinates In some rivers or reservoirs, the vertical momentum is generally governed by the hydrostatic variation as the result of the fact that the vertical velocity component is so small that it can be neglected compared with the horizontal ones. Thus the hydrodynamic equations and the pollutant transport equation are reduced to the shallow-water equations as follows. Continuity equation:

wu j wx j

=0

Pollutant transport equation:

§ wc · ¨¨ H j ¸¸ © wx j ¹

(3)

Free surface equation:

w] w ] w ] =  ³ udz  ³ vdz wt wx  H wy  H

(4)

in which ui denote the velocity components in the x, y, z-directions in Cartesian coordinates respectively, Fi the Coriolis force components per unit mass in the x, y-directions, ] the free surface water level against the reference plane z = 0 , Q t the eddy viscosity coefficients, c the concentration of pollutant, H j the turbulence diffusion coefficients in the x, y, z-direction respectively and H the water depth below the reference plane z = 0 . 2.2 Transformed equations in sigma-coordinate system In natural river flows, the bottom and the water surface are usually not parallel to the reference plane ( z = 0 ). In order to cope with uneven bottom topographies and free surface water level in numerical applications, a transformation to so-called sigma-coordinate system is applied[11]. This transformation stretches the vertical direction, such that the transformed water depth is constant both in space and in time, as shown in Fig.1. z

]

x

H

(a) z-coordinate

V V 

(1)

V  x (b) V-coordinate

Horizontal momentum equations:

wui wu w] w + u j i + Fi + g = wt wx j wxi wx j

wu j c) wc w +uj = wt wx j wx j

§ wui · ¨¨Q t ¸¸ © wx j ¹

(2)

Fig.1 Transformation from z-coordinate coordinate system

system to

sigma-

The sigma-coordinate system is defined by,

335

x = x , y = y , V =

z+H  , t =t h

(5)

where h = h( x, y, t ) = H ( x, y , t ) + ] ( x, y, t ) is the water depth. The shallow water Eqs.(1)-(4) transformed from the z-coordinate system to the sigma-coordinate system are written as

wu wu wu Z wu w] + u + v + = g  fv + wt wx wy h wV wx Fx +

1 w § V wu · ¨Q t ¸ h 2 wV © wV ¹

 wV wu · § wu · H § wu f xx = Q tH ¨ ¸, ¸ =Qt ¨  + © wx ¹ © wx wx wV ¹ (6)

Fxx =

wv wv wv Z wv w] + u + v + = g + fu + wt wx wy h wV wy (7)

  1 w  hu 1 w  hv w] = ³ dV  ³ dV 0 0 wt wx wy

(8)

Z = ³

V

0

 º 1 w  hv ª 1 w  hu dV + ³ dV » 0 0 wx wy ¬ ¼

(9)

wc wc wc Z wc 1 w + u + v + = Fc + 2 < wt wx wy h wV h wV § V wc · ¨H ¸ © wV ¹

(12)

wu wu = wxwV wVwx

(13)

Substituting Eq.(13) into the last term of Eqs.(12) and performing some algebraic operation yield (10)

where the terms Fx and Fy represent the horizontal viscosity terms in the horizontal momentum equations, Fc denotes the horizontal diffusion term in the pollutant transport equation and Z is the transformed vertical velocity in the sigma-coordinate system, defined by,

Z := h

wV w § H wV wu · ¸ ¨Q t wx wV © wx wV ¹

The above expressions represent the full transformation of the horizontal viscosity term without any reduction. Considering that the second-order derivative of u is continuous, we have the following cross-derivative relation in the x - V space

 V w  hv w  hu dV  ³ dV + 0 wx wy

V «³

wf xx wf xx wV wf xx w § wu · = + = ¨Q tH ¸+ wx wx wx wV wx © wx ¹

w § H wV wu · wV w § H wu · ¨Q t ¸+ ¨Q t  ¸ + wx © wx wV ¹ wx wV © wx ¹

1 w § V wv · ¨Q t ¸ h 2 wV © wV ¹

Fy +

2.3 Horizontal diffusion in sigma-coordinate system Mellor and Blumberg[12] noted that the full transformation formula for horizontal diffusion behaves poorly on sloping bottom and an approximation formula is needed. The formula proposed by Huang and Spaulding[13] is followed in this section. According to the rule of sigma transformation, the normal stress f xx and horizontal viscosity term Fxx in the x direction are expressed as

§ wV wV wV wV wV · = h¨ +u +v +w ¸ (11) wt wx wy wz ¹ © wt

Fxx =

w § H wu · w § H wV wu · ¨Q t ¸ + 2 ¨Q t ¸+ wx © wx ¹ wx © wx wV ¹ 2

§ wV · w § H wu · H wu 1 w < ¨ ¸ ¨Q t ¸ +Q t wV ¹ wV h wx © wx ¹ wV ©

wh · wV § wH V ¨ ¸+ wx ¹ wx © wx

§ wu wQ tH wu wQ tH · + ¨ ¸ wx wV ¹ © wV wx

The last two terms in the above expression can be neglected considering that the gradients of Q tH and

336

the second-order derivatives of H and h are relatively small. The final form of x-component of horizontal viscosity term in the x direction is given as

Fxx =

w § H wu · w wV wu [2Q tH + ¨Q t ¸+ wx © wx ¹ wV wx wx

Q tH (

wV 2 wu ) ] wx wV

(14)

Equation (14) is just the transformed form of horizontal viscosity used in the present sigma-coordinate model. Two terms are included in Eq.(14): one is the same as the conventional form, and another denotes the effect of steep bed bottoms. 2.4 Boundary conditions and turbulence closure model At the water boundary, the velocities and pollutant concentration are specified at the upstream boundary and the water elevation is prescribed at the downstream boundary. The stresses calculated by the mean of the wall function and the zero gradient condition of concentration are specified at the land boundary. Besides the impermeability conditions, wind stresses and friction forces are given at the free surface and bottom bed respectively. The horizontal eddy viscosity and the vertical one in the shallow water flow are distinguished. In this study, the vertical eddy viscosity is determined according to the formula proposed by Coleman[14] and the definition developed by Uittenbogaard et al.[15] is adopted for the horizontal eddy viscosity calculation, which is expressed as

Q tH = Q t2 D +Q tV

(15)

where Q t2D is the part due to “two-dimensional turbulence” computed by depth-averaged version of the k  H turbulence closure model proposed by Rastoqi and Rodi[16]. 3. Numerical schemes for governing equations The analytical solution of the governing Eqs.(6)-(10) can only be sought for very simple fluid flow problems, which however do not occur in natural rivers. A numerical approximation is therefore almost inevitable. In this work, a finite element scheme on a hierarchy grids and local time-step are chosen for solving the problem. Through layer-integration, the governing equations are uncoupled in the vertical dimension, so

that the three-dimensional problem can be reduced to a series of two-dimensional problem for simplicity. An approach of the operator splitting method combined with the Streamline Upwind Petro-Galerkin (SUPG) finite element methods on non-staggered grids is used to discretize the temporal domain and horizontal plane of each layer, and the trouble of numerical instability is overcame as a result. The detailed discussion of the numerical methods can be found in Ref.[17]. 4. Simulation of side discharges into open-channel flow In this section, results are presented for the model application to the side discharges into open-channel flow in situations studied experimentally by Mikhail et al.[18] for the verification of the present model. A typical flow situation was investigated for the case of a discharge at right-angles to the channel direction in this section. To compare with the experimental data of Mikhail et al.[16,18], the calculations were carried out for a discharge with width of b = 0.064m into a straight open-channel with width of B = 0.61m . A velocity of 0.1 m/s was specified at the inlet boundary of channel and the water depth at the outlet boundary was set to be 0.51 m. The discharge rates at the outfall were 0.05 m/s, 0.1 m/s, 0.2 m/s, 0.4 m/s and 0.8 m/s respectively. To ensure that boundary conditions have no significant influence on the predictions in the vicinity of the recirculation zone, it is necessary to move the position of the downstream boundary progressively further downstream until the size of the recirculation zone became unaffected by the downstream boundary location. The distance between the location of the discharge and the downstream boundary is 13.8 m, which is 22.6 times of the channel width. The grids used in calculations are shown in Fig.2, in which the mesh in the vicinity of the discharge were refined with the approach of the hierarchy grids for the accuracy of calculations and the size of the smallest element mesh reaches 0.008 m.

Fig.2 Computing grids for side discharge into open-channel flow

The contours of surface elevation for the case study with the discharge rate 0.2 m/s are plotted in Fig.3 and the velocity vectors in the same prediction situation are shown in Fig.5. Owing to the blockage

337

effect of the jet, the elevation for the channel water approaching the upstream of discharge rises in the near-bank region while it falls in the far-bank region. Consequently, the channel water flows around the jet and increases its velocity, which can be seen from the velocity vectors in Fig.4. The low pressure prevailing in the recirculation region which leads to the bending and reattachment of the jet can also be identified in Fig.3.

Fig.3 Predicted lines of constant surface level with the discharge rate 0.2 m/s (Unit: mm)

Fig.4 Predicted velocity vectors with the discharge rate 0.2 m/s

Fig.5 Variation of recirculation region with momentum flux ratio

Figure 5 presents the variation of recirculation region with the momentum flux ratio RM of the predicted non-dimensional length L / B and the width H i / B of the recirculation region, which are represented by the symbols of solid square. The

measurements of Mikhail et al. represented by the symbols of empty circle, are also included for comparison. In the figure, L denotes the length of the recirculation region, H i the width of the recirculation region, B the width of channel and R the momentum flux ratio. Although the predictions tend to underestimate somewhat the eddy length and overestimate the eddy width over the large range of momentum flux ratios, the results of predictions are consistent with those of experiments in the range of small momentum flux ratio. In view of some uncertainties both in experiments and predictions, the agreement between predictions and experiments can be regarded to be satisfactory and the present model is reliable for the predictions of side discharge into straight open-channel. 5. Simulation on Fuling Phosphate Fertilizer Effluent in Three Gorges Reservoir This section presents the model application to the river-side effluent of the Fuling Phosphate Fertilizer Factory in the Three Gorges Reservoir with the situations of field observation for the verifications of the model’s adaptability. Also the investigation of concentration distribution with different vertical resolution is performed. 5.1 Field observations In January, 1998, the field observation was carried out in the region near the Fuling Phosphate Fertilizer Factory, located on the right bank of the Yangtze River, by the Water Environment Monitoring Center for the Upper Reach of Yangtze River. Six sections were positioned for measurements along the river and the first section, labeled F1, was on 50 m upstream from the effluent of the factory while the others labeled F2, F3, F4, F5 and F6 on 20 m, 50 m, 100 m, 278 m downstream from the effluent respectively. With the flow rate of the Yangtze River of 4530 m3/s, the water surface evaluation of 137.40 m, and the COD concentration of 3.1 mg/L was measured on the section F1. The waste water at the flow rate of about 0.24 m3/s was discharge into the Yangtze River with the COD concentration of 45 mg/L and a pollutant mixing zone was formed at the downstream of the effluent[19]. 5.2 Results of flow simulation In computational domains, quadrilateral elements with refinement in the area near the effluent were employed in the horizontal plane and eight equally divided sigma layers were used in the vertical, so that the “hydrostatic consistency” condition was satisfied. With the situation of field observation mentioned above, the flow in the area near the effluent is predicted correctly and the predicted velocity vectors of surface layer are shown in Fig.6. It is observed that

338

the velocity in the mainstream area is higher than that in the region near the bank as the result of the influence of bottom topography of the river. With the refinement of hierarchy grids, a small backwater with the length of about 2 m is predicted just downstream of the effluent because of the wastewater discharge.

flow prediction and the predicted results and field observation are compared as well. With solid circles and line representing the field observation data and predicted concentration respectively, the comparison on the sections of F3 and F5 are shown in Fig.8 and the agreement between the predicted results and observed concentrations can be regard to be satisfactory for the reason of the uncertainty both in the field observation and in the prediction.

Fig.6 Predicted velocity vectors

The predicted velocity results and field observation data are compared and the comparison on the sections F3 and F5 are plotted in Fig.7, in which the symbols of solid circle represent the field observation data and the solid line the predicted results. Predicted results are in agreement with field observation data, though there are some discrepancies between them, especially in the region near the bank. In view of the topography errors and field observation errors, the agreement between the predicted results and field observation data appears acceptable.

Fig.7 Comparison of predicted results and field observation data for velocity

5.3 Results of COD concentration simulation The distribution of COD concentration is predicted under the condition same as that used in the

Fig.8 Comparison of COD predicted results and field observation data for concentration

The predicted contours of COD concentration in both surface layer and bottom layer are shown in Fig.9. According to the Surface Water Standard of China, the pollutant mixing zone inside which the COD concentration is higher than that of Class II of the Standard is formed with the length of 74 m and the width of 20 m on the water surface as a result of the wastewater discharge. Since the water depth in the effluent channel is relatively small and the wastewater is discharge into to the river at the surface layer, the sizes of the pollutant mixing zone in surface layer are much greater than those in bottom layer with the length of 56 m and the width of 16 m. 5.4 Sensitivity study of vertical resolution In the case study above described, if the vertical resolution was performed with only one layer, the present sigma coordinate model was reduced to a depth-averaged model and the distribution of concentration was identical along vertical direction. To investigate the effect of vertical resolution on the predicted results, the simulations for 1, 2, and 4 layers were performed besides the simulation with 8 layers accomplished in last subsection. The predicted mixing zones at the water surface with different vertical resolutions are plotted in Fig.6 and the length of mixing zones versus layer number is shown in Fig.10.

339

unchanged when the count of layers is beyond required number. Additionally, it is shown that the length of mixing zone with 1 layer is 9% less than that with 8 layers, which illustrates that the presented model is more suitable in determination of the scope of mixing zone.

Fig.11 Change of the length of mixing zones versus layers in sensitivity study of vertical resolution

Fig.9 Predicted concentration contours of COD in the region near the effluent of Fuling Phosphate Fertilizer Factory

Fig.10 Predicted mixing zone with different vertical resolution

Predicted mixing zones for different vertical resolutions in Fig.10 vary with the change of layer in the region of downstream of the effluent, but keep unchanged in the region both upstream of the effluent and along transverse direction. The changing curves in Fig.11 indicate that the length of mixing zone increases with the increase of layers but the increasing speed will slow down. The length of mixing zone with 2 layers experiences an increase of 4.7% compared with that with 1 layer and the length with 8 layers increases only 1.9% compared with that with 4 layers. On this account, the predicted result will remain

6. Conclusions A three-dimensional numerical model with sigma-coordinate for vertical dimension is developed for the prediction of pollutant mixing zone as the result of the riverside effluent. To avoid spurious diffusion caused by the conventional definition of horizontal diffusion unfitted for steep bottom slopes, a sigma-transformation is included. By the use of layer-integration, the water column is divided into several layers of horizontal planes and an approach of the operator splitting method combined with the SUPG finite element method on non-staggered grid is used to discretize the temporal domain and horizontal plane of each layer for easy adaptation to complex boundary configurations. A case study for the side discharge into a straight channel flow is predicted correctly compared with experimentally studied situations. The good agreement between the prediction and experiment shows that the numerical model is suitable for the prediction of side discharge into a straight open-channel flow. Furthermore, the model is also applied to the simulation of discharge wastewater from the Fuling Phosphorus Factory effluent into the Three Gorges Reservoir. The good agreement between predicted results and field observed data shows that the numerical model is suitable for the prediction of effluent discharge into natural river flows. The sensitivity study of vertical resolution shows that the presented model can provide a realistic calculation for the pollutant discharge and the fact that the present predicted mixing zone is larger than that predicted by depth-average models is valuable for environmental evaluation.

340

References [1]

[2] [3]

[4] [5]

[6]

[7]

[8]

[9]

JIANG Chun-bo, LI Kai and LI Ping et al. Finite element simulation of pollutant mixing zones in streams entering the Three Gorges Reservoir[J]. Journal of Tsinghua University (Science and Technology), 2004, 44(6): 808-811(in Chinese). QI Xue-chun. The application of adaptive grid algorithm in 2-D flow and water quality[D]. Beijing: Tsinghua University, 2002(in Chinese). CHEN Y., LIU Z. and SHEN M. The numerical simulation of pollutant mixing zone from riverside discharge outlet in Three Gorges Reservoir[C]. Process of the 6th Congress of ICHD. West Leederville, Australia, 2004, 401-408. MASK A. C, PRELLER R. H. A numerical simulation of East Asian Seas in March 2002: Effect of vertical grid choice[J]. Ocean Modelling, 2007, 18(2): 81-96. CHAO Xiao-bo, JIA Ya-fei and SHIELDS D. et al. Numerical modeling of water quality and sediment related processes[J]. Ecological Modelling, 2007, 201(3-4): 385-397. JIANG Y. W., WAI Onyx W. H. Drying-wetting approach for 3D finite element sigma-coordinate model for estuaries with large tidal flats[J]. Advances in Water Resource, 2005, 28: 779-792. JIN Kang-ren, JI Zhen-gang and JAMES R. T. Three-dimensional water quality and SAV modeling of large shallow lake[J]. Journal of Great Lake Research, 2007, 33(1): 28-45. WANG Chao, WANG Yan-ying and WANG Pei-fang. Water quality modeling and pollution control for the eastern route of South to North Water Transfer Project in China[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(3): 253-261. WU Wei, YAN Zhong-min. Evaluation of threedimensional numerical model for saline intrusion and purging in sewage outfalls[J]. Journal of Hydrodynamics, Ser. B, 2007, 19(1): 48-53.

[10] YANG Li-ling, WANG Yun-hong and ZHU Zhi-xia et al. Three-dimensional numerical model for winding tidal river with branches[J]. Journal of Hydrodynamics, Ser. B, 2007, 19(2): 249-254. [11] PHILLIPS N. G. A Coordinate system having some special advantages for numerical forecasting[J]. Journal of Meteorology, 1957, 14: 184-185. [12] MELLOR G., BLUMBERG A. Modeling vertical and horizontal diffusivities with the sigma coordinate system[J]. Monthly Weather Review, 1985, 113: 1379-1383. [13] HUANG Wen-rui, SPAULDING M. Modeling horizontal diffusion with sigma coordinate system[J]. Journal of Hydraulic Engineering, 1996, 122(6): 349-352. [14] COLEMAN N. L. Flume studies of the sediment transfer coefficient[J]. Water Resources Research, 1970, 6(3): 801-809. [15] UITTENBOGAARD R. E., Van VOSSEN B. Subgridscale model for quasi-2D turbulence in shallow water[C]. International Symposium on Shallow Flows. Delft, The Netherlands: Delft University of Technology, 2003, 572-582. [16] RASTOQI ASHOK K., RODI W. Two-dimensional mathematical model for the dispersion of heat in rivers[J]. SAE Special Publications, 1975, 4: 146-153. [17] LIU Zhao-wei. Studies on characteristics and carrying capacity of the water environment near the riverside in Three Gorges Reservoir[D]. Beijing: Tsinghua University, 2005(in Chinese). [18] MIKHAIL R., CHU V. H. and SAVAGE S. B. The reattachment of a two-dimensional turbulent jet in a confined cross flow[C]. Proceeding of 16th IAHA Congress. San Paulo, Brazil, 1975, 3: 414-419. [19] LI Chuang. Computation and analysis of mixing zone of pollutant for side discharge into natural river using FEM[D]. Master Thesis, Beijing: Tsinghua University, 1998(in Chinese).