search
for
 About Bioline  All Journals  Testimonials  Membership  News


International Journal of Environment Science and Technology
Center for Environment and Energy Research and Studies (CEERS)
ISSN: 1735-1472 EISSN: 1735-2630
Vol. 6, Num. 1, 2009, pp. 1-12

International Journal of Enviornmental Science and Technology, Vol. 6, No. 1, Winter, 2009, pp. 1-12

Numerical simulations of ice accumulation under ice cover along a river bend

1J. Wang; 2*J. Sui; 1P. Chen

1 School of Civil Engineering, Hefei University of Technology, Hefei, Anhui Province, China
2 Environmental Engineering Program, University of Northern British Columbia, 3333 University Way, Prince George, BC, Canada
*Corresponding Author Email: sui@unbc.ca

Tel.: +250 960 6399; Fax: +250 960 5845

Received 22 May 2008; revised 30 October 2008; accepted 25 November 2008

Code Number: st09001

ABSTRACT

In this paper, the k-έ two-equation turbulence model has been used to simulate ice accumulation under ice cover along a river bend. A 2D depth-averaged numerical model has been developed in a nonorthogonal coordinate system with nonstaggered curvilinear grids. In this model, the contravariant velocity has been treated as an independent variable. To avoid the pressure oscillation in the nonstaggered grids, the momentum interpolation has been introduced to interpolate variables at the interface. The discretization equations have been solved by using pressure correction algorithms. An equation has been developed for describing the deformation of ice jam bottom. The thickness distribution of ice accumulation (ice jam) along the bend has been simulated. The developed model has been applied to the experimental studies under different conditions carried out at the Hefei University of Technology. Results indicate that all simulated thickness of ice accumulation agrees reasonably well with the measured thickness of ice accumulation in laboratory.

Key words: Nonorthogonal coordinate system, finite volume method, ice jam, two-equation turbulence model, transform of body-fitted coordinate system

Introduction

River ice represents an important hydrologic element in temperate and polar environments. An ice cover alters the hydraulics of an open channel by imposing an extra boundary to the flow, altering the water level and flow velocity compared to the ice free case (Shen and Wang, 1995). Ice alters the flow regime, incipient motion of bed material (Wang et al., 2008) and the transport of sediment by temporarily storing water through ice jams and results in differences in riverbed deformation as compared to that observed under open flow conditions (Shen and Wang, 1995; Sui et al, 2000; 2005). Under specific conditions such as continuous incoming of frazil ice from upstream due to low temperature, ice jam may be formed. Ice jam can result in significant increase in water level and thus ice flooding such as the most recent ice jam flooding of the Nechako river in winter 2007 in downtown Prince George, Canada. Without doubt, river ice attracts interest of many researchers and engineers in the world (Shen and Yapa, 1998; Wang, 2002; White, 2003; Hains and Zabilansky, 2004).

River bend is one important location which plays an important role during the formation of jams and is also one important location where jams such as the Hequ ice jams in the Yellow river often occurs (Sui et al., 2002) . Thickness distribution of ice jam along a river bend is one important research topic in river ice hydraulics which is the least understoodable subject. Due to two transverse circulation cells under ice covered condition along a river bend, ice jam thickness varies both in flow direction and the transverse direction (normal to flow direction) (Sui et al., 2008).

Research work in the field of river ice hydraulics can be classified in following three categories: Field observations, experimental study and numerical simulation. Based on long term field observations of ice jams at the Hequ reach of Yellow river and experimental study of ice accumulation along a bend flume, Sui et al. (2006; 2008) investigated the mechanisms of ice accumulation along a river bend. Urroz and Ettema (1992; 1994) described the mechanism of ice jam initiation, as well as the features of ice jams in a curved channel based on experimental studies.

On the other hand, some numerical simulations on river ice hydraulics have been carried out in the past 2 decades (Shen et al., 1995). However, most of the existing numerical simulations have been conducted based on 1D model (Zufelt and Ettema, 2000; Hu et al., 2002; Mao et al., 2003; She and Hicks, 2006). So far, only few researchers have carried out 2D numerical simulation on river ice hydraulics (Shen et al., 2004).

In the past, some researchers pointed out that the accumulation of frazil particles under ice cover should be similar to the motion of bed load on riverbed (Pariset and Hausser, 1961; Ettema 2002; Wang et al., 2008). Starosolszky (1970) and Sui et al. (2000) suggested that the theory of bed load movement could be borrowed to describe the movement of ice under ice cover. Atkinson and Waters (1978) used Meyer-Peter formula to describe the movement of ice particles. However, there has not been a 2D numerical model describing ice accumulation by applying the theory of bed load movement. In this research work, the k_e two-equation turbulence model has been used to simulate ice accumulation along a river bend. The theory of sediment transport has been borrowed in development of the numerical model. The incoming ice discharge has been considered in the model. By developing an equation for the deformation of ice jam bottom, the thickness distribution of ice jam along a river bend has been simulated. The developed model has been applied to the experiments carried out at the Hefei University of Technology from 2001 to 2005. Results indicate that all simulated results agreed reasonably well with the experimental results.

Materials and Methods

The conventional equation for turbulent flow

The k-ε model, which is also called as the two-equation model, is the most commonly used method among all turbulence models. The k-ε model has been derived by assuming that the flow is fully turbulent and the effects of molecular viscosity are negligible. The standard k-ε model is, therefore, valid only for fully turbulent flows. The k-ε model is a semiempirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε ). The model transport equation for k is derived from the exact equation while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart. In the Cartesian coordinate system, the continuity and momentum equations of a 2D turbulent flow can be described by the following conventional equation:

Where, H is water depth; φ is the studied dependable variable; u and v are the depth-averaged velocity component in x and y direction, respectively; the source term Sφ = Sφc + Sφpfp , Gφ is the diffusion coefficient and varies with φ.

If φ = 1, Eq. 1 will be the continuity equation. If φ = u and v, respectively, Eq. 1 will be the momentum equations in x and y direction, respectively. If φ = k and ε , respectively, Eq. 1 will be the equation for the depth-averaged turbulence kinetic energy and the equation for the depth-averaged dissipation rate of turbulence kinetic energy, respectively.

Computational domain for body-fitted coordinate system

In this study, the finite volume method has been used for numerical solution. The abnormal boundary has been transformed by using the elliptic differential equations.

It is very difficult to solve Eq. 2 directly. To solve these equations, the coordinate is usually inversely transformed (Tao, 2001; Wang et al., 2006). The coordinates (x, y) of the physical domain should be determined from the coordinates of the computational domain in a body-fitted coordinate system (ζ, η ) (here,ζrepresents abscissa and η ordinate) as the following equation:

Where,

Transformation of the equations in the computational domain

Since the computational domain has been transformed from the Cartesian coordinate system to a body-fitted coordinate system, Eq. 1 should also be converted. By transforming the unsteady state term, convective term, diffusion term and source term, Eq. 1 can be written as following:

with,

Where, U and V are contravariant velocity components in x and y direction, respectively. The relevant variables φ, Gφ and Sφ (ζ, η) are shown in table 1.

The source terms in Table 1 are expressed as follow:

Where, the composite roughness coefficient (nc) of the flow cross section under ice covered condition in Eq. 12 can be expressed as following (Beltaos, 1995):

Where, Pb and Pi are wetted perimeter of riverbed and ice cover/jam, respectively. Hb and Hi are water depth from riverbed mainly affected by riverbed and from ice bottom mainly affected by ice jam, respectively. nb and ni are Manning roughness coefficient of

riverbed and ice jam, respectively.

According to Launder and Spalding (1974), the turbulence constants in the model from Eqs. 5 to 11 have following default values:

These default values have been determined from experiments with air and water for fundamental turbulent shear flows, including homogeneous shear flows and decaying isotropic grid turbulence. They have been found to work fairly well for a wide range of wall-bounded and free shear flows (Launder and Spalding, 1974).

Discretion of the general equations

The discretion of the general equations can be achieved by using the finite volume method. Fig. 1 shows one control element in the computational domain. The nonorthogonal grid has been developed. To simplify the simulation with Dx = Dh = 1, the discretional form of the general equations can be expressed as following:

Where,

Here, symbol |…| represents the absolute value. The subscripts E, N, W and S represent the center of the adjacent control elements right, up, left and down to the center of the studied control element P, respectively (Fig. 1). The subscripts e, n, w and s represent the right, up, left and down boundaries of the studied control element (the control element is shown with oblique lines), respectively.

Interpolation method

Generally, linear interpolation method and harmonic interpolation method can be used to calculate physical variables on the surfaces of the control element. For different physical variables, different interpolation methods should be used. Generally, the linear interpolation method can be used to calculate water depth, concentration and flow velocity in the nonorthogonal grid. The harmonic interpolation method can be used to calculate the turbulence kinetic energy and the dissipation rate of turbulence kinetic energy. However, the linear interpolation method for calculating flow velocity in a nonorthogonal grid often results in an unreal pressure oscillation. To avoid this problem (Xia, 2002; Wu, 2004; Wu et al., 2005), the momentum interpolation method has been used in this research.

Here, symbol represents the linear means of variables.

As shown in Eqs. 21 and 22, by introducing a pressure gradient along the interpolation direction into the momentum interpolation, the unreal pressure oscillation can be avoided.

Results and Discussion

For pressurable fluid, the momentum equation can be directly solved by using the state equation describing relationship between pressure and mass density. However, for impressurable fluid, the unknown pressure field makes it very difficult to solve the momentum equation. In the fluid, the pressure field is usually indirectly determined by solving the continuity equation. The following four methods are usually used for determination of the pressure field:

  • Solving the simultaneous Eq. 1.
  • Solving the pressure poisson Eq. 2.
  • Using the pseudo compressibility method Eq. 3.
  • Applying pressure correction algorithm.

The pressure correction algorithm is one widely used method that has been used in this study. In the depth-averaged 2D model, the static pressure is usually used in the simulation. Thus, the pressure gradient in the equations can be replaced by the water surface slope.

As suggested by Wu et al. (2003; 2005) and Xia (2002), the equations can be solved using the trail and error method by assuming a water level. Then, the difference between the actual flow velocity and the calculated velocity can be determined. By restricting the angle of intersection of the grid from 45° to 135° in the x-h coordinate system and ignoring the cross derivative terms, the following equations were used for the contravariant velocity components at the center of the control element P ( Fig. 1):

Correspondingly, the contravariant velocity components on the interface of the control element have the following values of correction:

Corresponding to the velocity field in the Cartesian coordinate system, the calibration equation for water level is:

A'P Z'P = AE Z'E + AW Z'W + AN Z'N + AS Z'N + b (26)

Where,

Where, U* and V* are the estimated contravariant velocity components on the calculation interface and U ´ and V ´ are the corrected contravariant velocity components, respectively.

Equations describing ice jam bottom

Applying the above equations, the accumulation thickness of ice jam can be computed. The 2D equation describing the variation in ice jam bottom is as as follows:

Eq. 28 can be expressed in the body-fitted coordinate system (ζ, η) as following:

Where, Zice is the elevation of ice jam bottom; eu is the porosity of ice accumulation (usually 0.4); qix and qiy are ice discharge in the control element in x and y directions, respectively which can can be calculated as follow:

Distribution of ice accumulation along a bend

The laboratory study was conducted using a horizontal 36 m length (including one concrete inlet section and one concrete outlet section), 0.5 m width, 0.6 m depth recirculating flume at the Hefei University of Technology (Sui et al., 2008). As shown in Fig. 2, the flume section having transparent tempered glass wall is comprised of one 180° bend with a radius of curvature of 1.5 m (central line) that is connected by a 16.8 m length straight reach upstream and a 4.8 m long straight reach downstream, respectively. From the upstream entrance cross section of the bend (0 °cross section described as section C1 in Fig. 2) to the downstream exit cross section (180°, as section C7), the bend reach was divided into 6 sub reaches by 7 cross sections set 30° apart.

Model ice particle made of wax has dimensions of 5 mm × 5 mm × 5 mm. Styrofoam plates were used to model ice cover. A hopper located at the upstream straight reach of the flume was used to discharge model frazil particles into the flume at controlled rates. For each run, a constant rate of ice supply discharged from the ice hopper was maintained during the experiment.

The thickness of ice accumulation and water level gradually approach constant. When the rate of ice supply equals the output rate of ice particles collected at the tailgate, the thickness of ice accumulation and water level at any given cross sections do not change. Therefore, the ice accumulation process can be considered to reach an equilibrium phase and then the measurements of accumulation thickness were carried out.

Water level and jam thickness were measured by a photoelectrical meter (MZL-A). The MZL-A is a digital point gauge which was mounted on a traversing system to obtain the water level and thickness of ice accumulation. Water levels were recorded directly on a data logger. The uncertainty of the measurements with dimensions was less than about ± 1 mm.

Eight experiments have been carried out to examine the effects of flow rate, water depth and ice discharge on the thickness distribution of ice accumulation along the bend flume. Experimental results have been used here for calibrating the developed numerical model.

In this numerical simulation, the nonorthogonal grid having 91´16 grid notes has been developed. The elevations of the ice jam bottom along the bend flume have been simulated. The developed numerical model has been applied to all eight experiments. Results show that the simulated results are very close to the experimental results. Fig. 3 shows the simulated jam thickness compared to the experiment under experiment conditions of approaching flow velocity of Uo = 0.15 m/s, ice discharge of Qi = 0.098 L/s and approaching water depth of H = 0.136 m. Fig. 4 shows the simulated jam thickness compared to the experiment under experiment conditions of approaching flow velocity of Uo = 0.0673 m/s, ice discharge of Qi = 0.092 L/s, and approaching water depth of H = 0.167 m. Fig. 5 shows the simulated jam thickness compared to the experiment under experiment conditions of approaching flow velocity of Uo = 0.1497 m/s, ice discharge of Qi = 0.122 L/s and approaching water depth of H = 0.116 m. As Figs. 3, 4 and 5 indicate that the simulated results agree reasonably well with the experiment results, it also appears that the developed numerical model can be used to successfully simulate the thickness distribution of ice accumulation along a river bend.

Conclusion

In this study, the k-ε model, which is also called as the two-equation model, has been used for developing a turbulence model for flow under ice accumulation along a bend channel. The 2D depth-averaged numerical model has been developed in a nonorthogonal coordinate system with nonstaggered curvilinear grids having 9116 grid notes. The finite volume method has been used for numerical solution. The abnormal boundary has been transformed using the elliptic differential equations.

In this model, the contravariant velocity has been treated as an independent variable. To avoid the pressure oscillation in the nonstaggered grids, the momentum interpolation method has been introduced to interpolate the variables at the interface. The discretization equations have been solved using pressure correction algorithms. Equations have been developed describing the deformation of ice jam bottom.

To assess the developed numerical model for ice accumulation along the bend channel, eight experiments have been carried out to examine the effects of flow rate, water depth and ice discharge on the thickness distribution of ice accumulation along the bend flume at Hefei University of Technology. The thickness distribution of ice accumulation along this bend has been measured.

Using developed numerical model, the elevations of the ice jam bottom along the bend flume have been simulated for all eight experiments. Results indicate that all simulated results agree reasonably well with the experimental results. It appears that the developed numerical model can be used successfully to simulate the thickness distribution of ice accumulation along a river bend.

Acknowledgments

The authors would like to acknowledge that this research work is financially supported by the National Natural Science Foundation of China (10372028).

References

  • Atkinson, C. H.; Waters, T. W., (1978). Ice regime at churchill falls, labrador - a comparison of design expectations with actual performance., Proc. IAHR symposium on ice problems, Lulea, Sweden, 2, 163-186.
  • Beltaos, S. (1995). River ice jams., water resources publication, LLC. Colorado, USA.
  • Ettema, R., (2002). Review of alluvial-channel responses to river ice., J. Cold Reg. Eng., 16 (4), 191-217.
  • Hains, D.; Zabilansky, L., (2004). Laboratory test of scour under ice: Data and preliminary results., US Army CREEL research report, tr-04-09, 77-118.
  • Hu, H.; Liu, C.; Chen, C.; Li, G.; Yang, K. A.; Liu, Z., (2002). Simulation of ice jam in river channels., Water Resour. Hydropower Eng., 33 (10), 40-47.
  • Mao, Z.; Wu, J.; Zhang, L., (2003). Numerical simulation of river ice jam., Adv. water Sci., 14 (6), 700-705.
  • Launder, B. E.; Spalding, D. B., (1974). The numerical computation of turbulent flow., Comput. Method. Appl. M., 3 (2), 269-289
  • Pariset, E., Hausser, R., (1961). Formation and evolution of ice covers on rivers., T. Eng. Inst. Canada, 5 (1), 41-49.
  • She, Y.; Hicks, F., (2006). Modeling ice jam release waves with consideration for ice effects., Cold Reg. Sci. Tech., 45 (3), 137-147.
  • Shen, H. T.; Liu, L.; Tuthill, A., (2004). Modeling ice passage at navigation locks., J. Cold Reg. Eng., 18 (3) 89-109.
  • Shen, H. T.; Wang, D., (1995). Under cover transport and accumulation of frazil granules., J. Hydraul. Eng-ASCE, 121 (2) 184-195.
  • Shen, H. T.; Wang, D. S.; Wasantha, L. A. M., (1995). Numerical simulation of river ice processes., J. Cold Reg. Eng., 9 (3), 107-118.
  • Shen, H. T.; Yapa, P. D., (1998). Flow resistance of river ice cover., J. Hydraul. Eng-ASCE, 112 (2), 173-183.
  • Sui, J.; Hicks, F.; Menounos, B., (2006). Observations of riverbed scour under a developing hanging ice dam., Can. J. Civil Eng., 33 (2), 214-218.
  • Sui, J.; Karney, B.; Fang, D., (2005). Variation in water level under ice-jammed condition _ Field investigation and experimental study., Nord. Hydrol., 36 (1), 65-84.
  • Sui, J.; Karney, B.; Sun, Z.; Wang, D., (2002). Field investigation of frazil jam evolution _ A case study., J. Hydraul. Eng-ASCE, 128 (8), 781-787.
  • Sui, J.; Wang, J.; Balachandar, R.; Sun, Z.; Wang, D., (2008). Accumulation of frazil ice along a river bend., Can. J. Civil Eng., 35 (2), 158-169.
  • Sui, J.; Wang, D.; Karney, B., (2000). Suspended Sediment concentration and deformation of riverbed in a frazil jammed river reach., Can. J. Civil Eng., 27 (6), 1120-1129.
  • Starosolszky, Ö., (1970). Ice in hydraulic engineering., report No.70-1, Division of hydraulic engineering, University of Trondheim, Norway.
  • Tao, W., (2001). Numerical heat transfer, 2nd. (Ed.), Press of Xian Jiatong University, China, 353-357.
  • Urroz, G.; Ettema, R., (1992). Bend ice jams: Laboratory observations., Can. J. Civil Eng., 19 (5), 855-864.
  • Urroz, G.; Ettema, R., (1994). Small-scale experiments on ice-jam initiation in a curved channel., Can. J. Civil Eng., 21 (5), 719-727.
  • Wang, J., (2002). A study on ice jams in balance transporting discharge. J. Hydroelectronic Eng, 26 (1), 61-67.
  • Wang, J.; Sui, J.; Karney, B., (2008). Incipient motion of non-cohesive sediment under ice cover _ an experimental study., J. Hydrol, (in press).
  • Wang, J.; Zhang, C.; Ni, J., (2006), Methods of generating two-dimensional body-fitted grids, transactions of Hefei University of Technology, 12, 1549-1551.
  • White, K. D., (2003). Review of prediction method for breakup ice jams., Can. J. Civil Eng., 30 (1), 89-100.
  • Wu, X., (2004)., Turbulent model for distribution and transport of contaminants in water in curvilinear coordinate system., Ph.D. thesis.
  • Wu, X.; Shen, Y.; Zheng, Y., (2003). 2-D flow SIMPLEC algorithm in non-orthogonal curvilinear coordinates., Chinese J. Hydraul. Eng., 34 (2), 25-30.
  • Wu, X.; Shen, Y.; Zheng, Y., (2005). Algebraic-stress turbulent model for 2-D plane flow in curvilinear coordinates., Chinese J. Hydraul. Eng., 36 (4), 383-390.
  • Xia, Y., (2002). Research on application of 3-D hydrodynamic, sediment transport model with non-staggered curvilinear grid for tidal rivers., Ph.D. thesis, Hohai University, Nanjing, China.
  • Zufelt, J. E.; Ettema, R., (2000). Fully coupled model of ice-jam dynamics, J. Cold Reg. Eng., 14 (1), 24-41.

© Copyright 2008 IRSEN, CEERS, IAU


The following images related to this document are available:

Photo images

[st09001f1.jpg] [st09001f5.jpg] [st09001t1.jpg] [st09001f3.jpg] [st09001f4.jpg] [st09001f2.jpg]
Home Faq Resources Email Bioline
© Bioline International, 1989 - 2024, Site last up-dated on 01-Sep-2022.
Site created and maintained by the Reference Center on Environmental Information, CRIA, Brazil
System hosted by the Google Cloud Platform, GCP, Brazil