The total flow control equations and the characteristics of unsteady gradually varied flow in open channels ＊

Shi-he Liu, Qian-yi Zhao, Qiu-shi Luo

1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

2. Yellow River Engineering Consulting Co., Ltd, Zhengzhou 450003, China

Abstract: Most unsteady channel flows in nature and practical engineering appear as gradually varied ones, and in the free surface,the deformation conforms to the long wave hypothesis. One-dimensional total flow models were usually used to for the numerical simulation of long-term and long-distance reaches to describe the water movements, however, the models lack a clear relationship between the three-dimensional flow field and the total flow field. Moreover, few studies of the variations of the roughness coefficient against the time in unsteady flows were conducted. The following results are obtained through the theoretical analysis and the numerical simulations in this paper. (1) One-dimensional total flow control equations of the unsteady gradually varied flow in open channels are obtained directly from the mathematical model of the viscous fluid motion, and can both reflect the influence of the turbulence and provide an explicit expression of the energy slope term. These equations establish a direct connection between the descriptions of the three-dimensional flow fields and the one-dimensional total flows. (2) Synchronous prototype observation data and planar two-dimensional numerical simulation results are used to extract the one-dimensional total flow information and discuss the total flow characteristics. (3)The orders of magnitude for terms in the total flow motion equation are compared, and the variation of the roughness coefficient against the time is analyzed.

Key words: Open channel unsteady flow, total flow characteristics, control equations

Introduction

The open channel flow[1-3]is defined as the flow that is driven by the gravity with a free surface. The flow within tidal rivers is a type of open channel flow,which is influenced by the runoff from upstream and the tides from downstream and with strong unsteady characteristics. Most flows within tidal rivers gradually vary, with the vertical and horizontal changes of the free surface conform to the long wave hypothesis.

The unsteady flow in open channels was extensively studied through experiments, observations of prototypes and numerical simulations. Through experiments and prototype observations, a substantial amount of data was accumulated with an enriched understanding of the unsteady flow in open channels.These include, for example, Nezu’s research[4-6]of the turbulence characteristics of the unsteady flow, as well as Graf and Song’s research[7-8]regarding the bed-shear stress and the distribution of the flow velocity in an unsteady flow. However, to date,systematic and complete data for the open channel unsteady flow are still few due to the complexity of its flow field characteristics, which change with time and space, and the limited means of measurement.Numerical calculations have increasingly evolved to become a primary method to study the open channel unsteady flow. Due to the effects of the turbulence closure mode and the free surface treatments, there was hardly any progress in the three-dimensional simulation[9-12]. The planar two-dimensional mathematical models were mainly adopted[13-14]based on the solution of the two-dimensional shallow water equations or the one-dimensional total flow models[15-18]based on the solution of the Saint-Venant equations for the numerical simulation of the unsteady flow in rivers. A one-dimensional total flow model comprises both the continuity equation and the motion equation.The continuity equation is obtained directly from the principle of mass conservation, while the motion equation, namely the Saint-Venant equation, is derived from either the theorem of momentum in the classical mechanics or the total flow energy equation[1-2]. The aforementioned motion equation,though still in use, cannot directly account for the influence of the turbulence and be used to obtain an explicit expression for the resistance. In addition, with the motion equation, the connection between the flow field description and the total flow description can not be established. The integral control equations for the homogeneous incompressible liquids in the open channel steady flow were constructed from the mathematical model of the viscous flow motion[19-20],but without much discussion on the differential control equations of the unsteady flow in the open channel. This paper investigates the open channel unsteady gradually varied flow within tidal rivers and has obtained the following results through the theoretical analysis and the numerical simulations. (1)The total flow control equations of the open channel unsteady gradually varied flow are obtained directly from the mathematical model of the viscous fluid flow.The equations reflect the influence of the turbulence and a connection between the field flow description and the total flow description is established directly.(2) The planar two-dimensional numerical simulations are conducted on the basis of the synchronous prototype observation data in a typical tidal river in order to extract the total flow information of the water level and the mean velocity in cross-sections. (3)Based on the total flow information, the order of magnitude of each term in the total flow motion equation is compared, and the variability of the roughness coefficient against the time is analyzed.

1. The existing control equations for open channel unsteady gradually varied flow

The existing control equations for the open channel unsteady gradually varied flow are deduced from the principle of mechanics. The continuity equation is based on the mass conservation principle,while the motion equation is based on the momentum theorem of the rigid body motion in the classical mechanics. They are discussed in the following subsections[2].

1.1 The continuity equation

For the open channel unsteady gradually varied flow, the variables Q, A, h and B, respectively,represent the discharge, the cross-sectional area, the depth and the width. The l, t terms, respectively,denote the longitudinal coordinates and the time.Because the change rate of the discharge with respect to the longitudinal distance is ∂Q / l∂, and the change rate of the depth with respect to the time is ∂h / t∂,the change of the water discharge in the time dt and the longitudinal distance dl is ( ∂Q / ∂ l )d ld t, and the corresponding change of the open channel storage volume is Bd l ( ∂ h/ ∂ t )= ( ∂ A/ ∂t )d ld t. The continuity equation can be obtained accordingly under the incompressible liquid conditions:

or

1.2 The motion equation

The open channel flow is driven by the gravity.U represents the section-averaged flow velocity, and the acceleration of a unit mass of fluid is ∂ U/ ∂t+U ( ∂ U/ ∂l). The external forces on the liquid in the process of motion include the gravity, the pressure and the wall resistance, which are expressed as.

Pressure: Including the resultant force of the pressure on the upper and lower sections, as well as the resultant force on the surrounding interface. Under a hydrostatic approximation,

The wall resistance: Estimated as the constant uniform flow, that is

where zb, χ, R, respectively, represent the elevation of the channel bottom, the wet perimeter and the hydraulic radius, and ib= -∂ zb/∂l is the open channel slope. From the theorem of momentum

Simplifying formula (2), we can obtain

Formula (3) represents the currently widely-used motion equation for the unsteady flow in an open channel, known as the Saint-Venant equation, whererepresents the water level, and Sf=represents the energy slope.

Formula (1) and formula (3) together constitute the existing differential hydraulic model for the unsteady gradually varied flow in open channels. The control equations are closed, but one has to supplement an empirical formula of the energy slope Sf.

Fig. 1 Sketch of control volume in the river

2. A new total flow control equation for unsteady gradually flow in open channels

Consider an open channel flow, as shown in Fig.1 (within a tidal river), with a control volume V at the time t. The control volume consists of two sections1S, S2, spaced at an interval of dl, the boundary of the river S31and the free surface S32have no inlet or outlet flows. S represents the surface of the control volume bounded by1S, S2,S31and S32. θ is the angle between x1axis (in the vertical direction) and the horizontal direction.The water in the river can be regarded as a homogeneous and incompressible fluid with a density of ρ and a dynamic viscosity coefficient of μ flowing within a gravity field. The control equation of the ensemble average is given as follows[19-20]:

The variablein in formula (6) is the unit vector along the outer normal direction of the water surface.The dynamic boundary conditions of the flow encompass the following. First, the time-averaged shear stress does not exist on the water-air surface,upon which the average normal stress is provided by the atmospheric pressure. Second, the flow in the river sees a small difference between the liquid flow velocity and the gas flow velocity on the water-air surface with a large radius of curvature. The corresponding dynamic boundary conditions are met.

The continuity equation and the motion equation of the one-dimensional total flow in the open channel under the aforementioned boundary conditions will be discussed in the following sections.

2.1 Total flow continuity equation

The width of the river is taken as B = ∂A / z∂,and thus

This formula can subsequently be modified into the following

The total flow continuity equation will then be obtained as follows

Formulas (7), (1b) are entirely consistent.

2.2 Total flow motion equation

(1) Integral of the first term on the left side of the equation

(2) Integral of the second term on the left side of the equation

Substituting formula (7) into the above formula,the following is obtained

Thus

(3) Integral of the third term on the left side of the equation

(4) Integral on the right side of the equation

Defining the energy slope as Sf, the integral of the right side of the equation can be represented as follows

Substituting formula (8) into the integral of formula (5) on a control volume V, and simplifying it with the continuity equation

Fig. 2 Map of the typical tidal river

Fig. 3(a) Comparison between calculated and measured values of water level at each measuring point

Fig. 3(b) Comparison between calculated and measured values of water level at each measuring point

The above new total flow motion equation contains one additional term, as compared to the Saint-Venant equation (the second term on the right side of formula (9)). This term is related to the changes of the flow cross the sectional area, and is derived from the change of momentum caused by the free surface deformation（the free surface deformation does not induce a mass change within the control volume). The term does not exist in formula (3)because a rigid body theorem of momentum is used to derive the motion equation. In addition, in contrast to formula (3), the first term (the energy slope item) on the right side of formula (9) can be calculated directly by using formula (8d). Thus, the relationship between the three-dimensional flow field description and the one-dimensional total flow description for the unsteady flow within open channels can be obtained.For an open channel steady flow, formula (9) can be simplified as follows

Formula (10) is completely consistent with the existing total flow motion equation for the open channel steady flow[1-2].

3. Total flow characteristics of unsteady flow in tidal rivers

Experiments and prototype observations can not provide systematic, complete and highly precise threedimensional data of the instantaneous flow field characteristics for the unsteady flow, therefore, we adopt the following steps in order to obtain onedimensional total flow information: (1) based on the synchronous measured water level and the distribution of the mean velocity in the cross-section, twodimensional numerical simulations are employed to calculate the temporal and spatial variation characteristics of the water level, and in turns, the temporal and spatial variations of the cross-section area, (2) the continuity equation (7) is solved to obtain the changes of the mean velocity in the cross-section against the time. The following example is an illustration of the Nanjing reach of the Yangtze River, a typical tidal river.

3.1 Planar two-dimensional numerical simulation of flow motion in the Nanjing reach of the Yangtze River

A map of the Nanjing reach of the Yangtze River is shown in Fig. 2. The section of the river considered here is roughly 37 km long and 2 km wide and includes 4 separate beaches. We obtain the topography data for the region together with the prototype observation data for six measuring points in addition to the observation data of the instantaneous velocity at three measuring sections in this reach. The locations of the water level measurement sites and the flow velocity measuring sections are shown in Fig. 2. We conduct a two-dimensional numerical simulation of the flow in this reach utilizing the self-developed numerical simulation system RSS. In view of the characteristics of the dense beachfront regions in this reach, a non-structural triangle grid is used in the calculation. 31 548 grid nodes and 62 395 calculation units are arranged within the calculation area. The maximum grid spacing is 100 m with a minimum spacing of 50 m. Brief introductions to the calculations and the verification results of both the water level and the flow velocity are provided in the following sections.

3.1.1 Verification of the water level

The measurement process for determining the water level is consistent with the calculation process of the water at each of the six measuring points. The absolute error is less than 0.01 m. A comparison of the measured and the calculated water levels at each measurement point is given in Figs. 3(a), 3(b), from which it is clear that the calculated values are consistent with the measured values.

3.1.2 Verification of the flow velocity

The measured depth-averaged velocity fits well with the calculated velocity at the three measuring sections. The error of both values does not exceed 0.1m/s. Due to the space limitation, Fig. 4 provides a comparison of the measured and the calculated results of the flow velocity for the right branch of the Xinqianzhou section.

Fig. 4 Comparison between measured depth-averaged velocities and calculated velocities for the right branch of Xinqianzhou section

3.2 Total flow characteristics of flow within the tidal river

3.2.1 Extraction of the total flow information

Under the consideration that the planar twodimensional numerical simulation results can fundamentally reflect the flow characteristics within the Nanjing reach well, we extract the relevant onedimensional total flow information in formula (9)based on the calculated results for the water level and the flow velocity. The typical section is relatively straight and is located near the Cihu Estuary as shown in Fig. 2. We separately consider each calculation section from the upstream and downstream sections of the typical section to calculate the spatial derivatives.The distance between each calculating section and the typical section is Δl = 389 m. The water level Z,the cross-section area A and the section-averaged flow velocity U of the three sections mentioned previously are obtained.

In the physical hydrological measuring, the measurements of the water level are generally more accurate, while the measurements of the cross-section velocity show typically large deviations from the actual values due to the limitations of the measurement technique. We observe that the calculated water level values are in good agreement with the measured values derived from the verification results of the two-dimensional numerical simulation above. It is believed that the calculated water level is consistent with the real water level. Meanwhile, some error exists within the measured flow velocity. Since the real flow must satisfy the continuity equation,however, we utilize formula (7) jointly with the calculated value of the water level to correct the section-averaged flow velocity U. We believe that the modified flow velocity value is essentially equivalent to the actual value. After the modification,the one-dimensional total flow information of a typical cross-section as well as the upstream and downstream sections are obtained, and are shown in Fig. 5 (wherein the water level is represented with the vertical coordinate on the left, while the flow velocity is represented with the vertical coordinate on the right).

Fig. 5 The one-dimensional total flow information of a typical cross-section

3.2.2 Description of the total flow characteristics

Figure 5 reflects the one-dimensional total flow information of a typical cross-section within this tidal reach, as it is shown, which is discussed as follows:

(1) The flow within the typical reach manifests as an irregular half-day wave. The water level and the section-averaged flow velocity of a typical section have a periodic trend with a period of approximately 12 h.

(2) In one wave cycle, the variation of the water level and the section-averaged flow velocity against the time can be refined into two stages: in stage 1, the water level increases to the wave crest from the wave trough, while the section-averaged flow velocity decreases to the wave trough from the wave crest, in stage 2, the water level decreases to the wave trough from the wave crest, while the section-averaged flow velocity increases to the wave crest from the wave trough.

(3) The variations of the water level and the section-averaged flow velocity against the time are not synchronized in the same tidal wave cycle, but with a phase difference.

4. Comparison of the order of magnitude of each term in the motion equation

Based on the total flow information extracted from the above planar two-dimensional numerical simulation results, we compare the order of magnitude of each term in the motion equation (9). In formula (9),(ξ / g )( ∂ U / ∂t ) represents the time-varying inertia term, (α U / g )( ∂ U / ∂l ) denotes the position-varying inertial term, -∂ z/ ∂l is the water surface gradient,Sfrepresents the energy slope term, and η( U / gA) ( ∂ A / ∂ t) denotes the surface deformation term. We choose the following values for the parameters: ξ =1.0, α =1.0, η = 0.15. The relevant calculation results are shown in Fig. 6 (the dimensionless coordinates on the right hand side are to reflect the term’s values in the motion equation),which also provides the water level variations(presented with the vertical coordinate on the left),and are as follows:

(1) In the motion equation, the orders of magnitude of the river gradient term and the energy slope term are much larger than those of the other terms, which suggest that the gravity and the resistance are the controlling factors for the flow in the typical reach.

Fig. 6 Comparison of the orders of magnitude of each term

5. The variation of the roughness coefficient against the time in the motion equation

As mentioned previously, although the control equations of the differential model of the unsteady gradually varied flow in an open channels are closed,they still require the addition of the empirical formula for Sf. From formula (5d):, where

From a physical perspective, Sf1represents the energy slope term induced by a variation of the surface force (including the pressure, the shear stress and the Reynolds stress) on the inlet and outlet sections, while Sf2denotes the energy slope term related to the total flow mechanical energy loss. In the hydraulic calculations, a roughness coefficient is often used to reflect the change of the energy slope, namely,where U, R, respectively, represent the section-averaged flow velocity and the hydraulic radius of the cross-section A[1]. Figure 7 shows the variation of the roughness coefficient against the time for two wave cycles, as well as the water level against the time (presented with the vertical coordinate on the left), described as follows:

(1) In one wave cycle, the range of the variation for the roughness coefficient n is 0.015-0.035. The time-averaged roughness coefficient (dashed line in Fig. 6) is n = 0.0234 , the amplitude of the maximum value is ( nmax-n )/ n =38.6%and the amplitude of the minimum value is (nmin- n )/ n= -3 7. 6%.That is, the relative amplitude of the roughness coefficient does not exceed 40% of the time-averaged value for one wave cycle.

Fig. 7 Variations of Z, n against the time in different cycles

(2) In one wave cycle, the variations of the roughness coefficient and the water level against the time can be refined to four stages: in stage 1, the roughness coefficient increases to the maximum from the minimum, and the water level rises to the wave crest from an early rising period after the wave trough,in stage 2, the roughness coefficient decreases quickly to the time-averaged value from the maximum, and the water level decreases from the wave crest, in stage 3, the roughness coefficient sustains a value nearly equal to the time-averaged value, and the water level decreases to the wave trough; in stage 4, the roughness coefficient decreases quickly from the time-averaged value and the water level returns to the early rising period from the wave trough.

(3) In one wave cycle, the roughness coefficient changes around the time-averaged roughness coefficient (namely) in nearly 50% of the wave cycle.

6. Conclusions

The following progresses are made through theoretical analysis and numerical calculation:

(1) The total flow control equations of the unsteady gradually varied flow in open channels are obtained directly from the mathematical model for the viscous fluid motion, with the continuity equation consistent with the existing total flow continuity equation. The motion equation is generally consistent with the existing total flow motion equation, but since it is derived directly from the flow field control equations, it can reflect the changes of momentum in the control volume induced by the free surface deformation. The new control equations reflect the influence of the turbulent fluctuations, and provide an explicit expression of the roughness coefficient, and a relationship between the flow descriptions of the three-dimensional flow field and the one-dimensional total flow can be established.

(2) Numerical simulations are conducted for a typical tidal reach with a planar two-dimensional mathematical model based on the synchronous prototype observation data. The flow field information is calculated using the mathematical model, and the total flow information is extracted from the flow field calculation results. The results show that: (a) affected by the tide, the variations of the water level and the section-averaged flow velocity of the typical section in the tidal reach against the time are not synchronized in the same tidal wave cycle, but characterized by a phase lag, (b) considering the terms in the total flow motion equation, the magnitude of the surface deformation term is substantially smaller than that of the other terms. Moreover, the inertial terms are less significant than the water surface gradient and energy slope terms, indicating that the gravity and the resistance are the controlling factors for the water movement in this typical reach.

(3) The variation of the roughness coefficient of the unsteady gradually varied flow in the tidal reach against the time has the following features: (a) the roughness coefficient changes, in a single wave cycle,in a sequence of a rapid increase, a rapid decrease, a time-invariant plateau and a fast decrease, while it changes around the time-averaged roughness coefficient (the time-invariant or flat stage) during nearly 50%of the wave cycle, (b) the relative amplitude of the roughness coefficient oscillation does not exceed 40%of the time-averaged roughness coefficient in one wave cycle.