|
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 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:
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
© 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] |
|