Analytical solution of non-uniform flow in compound channel

The present work investigates the analysis of depth-averaged velocity and boundary shear stress distribution in compound channels with non-uniform flow condition. A quasi two-dimensional model is proposed to assess the flow variables by accounting the physical processes that are specific to non-uniform flow. For analyzing the flow behavior, experimental data sets concerning compound channels with narrowing and enlarging floodplains of previous investigators are considered. The model accounts the influence of momentum transfer on the flow variables through additional shear stresses that are developed in non-uniform flow. Three types of effective stresses produced by molecular viscosity, turbulent and dispersion on the vertical planes are discussed. An analytical solution to the model is presented. Terms associated with the effective stresses are investigated relating them to the geometric and hydraulic parameters. The significance of lateral variation of energy slope is further investigated. For both homogenous and heterogeneous non-prismatic channels, the approach is examined to predict the flow variables with reasonable accuracy.


Introduction
Studies on the momentum exchange in compound river channels are focused for many decades especially in case where the overall channel width is constant along the length.These studies are commonly handled by simplification of the Navier-Strokes equation, to determine the transverse distribution of velocity and boundary shear-stress for steady uniform flow.The enhanced versions of these methods account for bed friction, lateral turbulent shear stress and secondary current term.However, most of compound rivers consist of non-prismatic floodplains.The flow variables have been numerically described by various researchers for this case.However, analytical studies on gradually varied flow are limited.For these channels, Bousmar and Zech [1] split up the secondary current term into two parts ; (1) dispersion term in uniform flow resulting from helical secondary currents which is similar to Γ in the Shiono and Knight Method [2]; (2) transverse convection term resulting due to the mass transfer because of non-prismaticity of the channel.Mc Gahey and Samuels [3] developed a general method by incorporating both these secondary-current terms for the straight and non-prismatic compound channels.Further, Bousmar and Zech [4] extended the lateral distribution method (LDM) to nonprismatic compound channels and explained the secondary-flow term in the converging compound channels.In many practical problems, accurate computations of the water elevation and flow variables in various sections of these channel is required.Numerical solution to the 1D Navier-Stokes equation are already investigated [3,4].However, the detail approximate solution has not been provided.So, an analytical solution is obtained for the depth-integrated turbulent form of Navier-Stokes equation for non-uniform flow to predict the depth averaged velocity and bed shear stress.The governing equation carries the energy slope in the gravitational force term and includes the dispersion and transverse convection in the secondary flow term.Independent calibrating coefficients are chosen for each subsection.The Reynolds shear stress term is modeled using the effective eddy viscosity concept that includes both horizontal coherent structures and three-dimensional bottom turbulence.To illustrate the implementation of proposed solution, laboratory data sets of compound channels of both converging and diverging floodplains are considered.

Methodology
In this study, two non-prismatic configurations are considered ; (1) converging geometry, where, the upstream geometry gets reduced towards downstream (Fig. 1a) and ( 2) diverging geometry, where, the upstream geometry gets enlarged along flow direction (Fig. 1b).(2) where,   =friction slope; the above equation is written by replacing the turbulent velocity fluctuations ' in x and ' in y direction by  −   and −   respectively [5], where,  and are the respective temporal mean velocity components in  and  direction,   and   are the depth averaged term of stream wise and lateral velocity.The last two terms are the dispersion terms.For steady and uniform flow, the 2D Saint-Venant momentum equation is The 1 st term in the left hand side is gravitational term; 2 nd term is the Reynolds shear stress, 3 rd term is the bed shear stress, 1 st term in right hand side is the dispersion term and 2 nd term is the convection term due to the molecular viscosity of secondary current.Here, the total secondary flow term is divided into two terms.Firstly, due to the helical secondary current known as dispersion and secondly due to the mass transfer known as transverse convection.However, for non-uniform flows, the variation of pressure and acceleration terms in the stream wise direction would be significant and should be accounted in this case.Bousmar and Zech [4] also stated that the stream wise turbulent friction    and dispersion terms   −      could be anticipated to remain small and the energy slope (  ) is no more equal to the bed slope ( ).Bed slope  ) together with the acceleration term For non-prismatic channels, it is essential to analyse the two profiles such as water surface line and energy gradient line.So, a brief idea about their estimations are described below.

Calculation of water surface profile and energy slope
Due to non-prismaticity of channel geometry, the flow depth varies in longitudinal direction and consequencely, the flow variables vary along the flow direction although the discharge remains constant.It is found that, the water surface profiledepends upon the geometry and flow conditions of a compound channel.Investigating several diverging compound channel series, Das and Khatua [6] modeled the non-dimensional water surface profile as a function of width ratio (), relative flow depth   ), relative longitudinal distance   ) and diverging angle of floodplain ().Where,  is the ratio of the total flow depth in the main channel () to the bank full depth (),  is the ratio between the total compound channel width to the main channel bottom width,   is the ratio of floodplain flow depth ( − ) to the main channel flow depth (),   is the ratio of the distance () of a arbitrary reach to the total length () of the non-prismatic reach and  is the angle of the floodplain to the main channel.Utilizing multivariable regression technique, Das and Khatua [6] obtained a relationships for non-dimensional water surface profile ( ) as For determining the water surface profile in converging compound channels, the drop in water surface (ΔH) between two desired sections is given by Rezaei [7] as  (6) where,   and   ,   and   are the corresponding discharge and average flow velocities in main channel and floodplains.  and   are the change in sub-sectional discharge and flow velocity in main channel,  ′ = momentum transfer coefficient,  and  refer to a section and its forward section respectively.So the flow depth () of a desired section can be determined by adding  with the flow depth of upstream section.
For computation of energy slope (  ) in non-prismatic compound channels, Proust et al. [8] explored the role of an explicit modelling of mass conservation.The energy loss between two sections can be calculated from the energy conservation equations.Estimating the energy loss ∆ between two sections,   can be calculated as   ∆   −   .where,  is the distance and  and  are the total energies at those sections respectfully.Das and Khatua [6] also modeled the   in terms of non-dimensional parameters as

Secondary current and Reynolds stress modeling
For non-prismatic geometries, a complete secondary flow term is taken that includes dispersion and mass-transfer or convection terms.Ervine et al. [9] have taken the secondary current term proportional to the square of the velocity (  ) in view of both prismatic and non-prismatic flows.So, the assumption of Ervine et al. [9] is considered here for the convection term of the equation ( 4) which takes care of the lateral momentum exchange.
The assumption is that the depth averaged lateral velocity component   is a function of   i.e.,     .So the convection term is now rewritten as           .where,  is an empirical coefficient pertaining to the geometrical and roughness conditions of an open channel.The non-uniformity of the vertical profiles of stream wise and lateral velocities is manifested in terms of dispersion term.Shiono and Knight [2] assumed the dispersion term as a fraction of bed shear stress in prismatic compound channels.So, for non-prismatic compound channels, the dispersion term is considered according to Rezaei [7] as where,  is a dispersion coefficient.
In the shear layer regions, turbulent shear stress is predominant.Using the Boussinesq eddy viscosity model, the 2 nd term of the equation ( 4) is simplified by Shiono and Knight [2].They estimated   in terms of shear velocity ( * ) and flow depth () as     *   .where,  is the eddy viscosity coefficient due to the bottom friction.However, Van Prooijen et al. [10] stated that Reynolds shear stress is affected by coherent structure and bottom friction as well.So, in the present study, both of these contributions are taken in account.Firstly, Reynolds shear due to bottom turbulence is evaluated using the Boussinesq eddy viscosity approach.Previously, some mixing length approaches are proposed taking the mixing layer width as a length scale and the velocity difference between subsections as velocity scale.Here, Reynolds shear due to horizontal coherent structure is evaluated in terms of shear layer width () as length scale and   as velocity scale using the Prandtl's mixing length theory i.e.,        , where,  is an eddy viscosity coefficient due to coherent structure.Due to reduction in momentum exchange, the width of mixing layer reduces with increase in flow depth.As the turbulence at shear layer width varies with geometric and flow parameters, so, it is confirmed that the shear layer width depends upon these parameters.Devi and Khatua [11] An analytical solution to equation ( 8) is attempted for estimating   which is given by   4) upstream unit discharge  is known.Then, boundary shear stress   can be estimated from the predicted   by using the relationship     .Equations 8 and 10 depend upon five calibrating coefficients ,  ,  , and , which are related to local bed friction, eddy viscosity coefficients due to bottom friction and horizontal coherent structure, secondary flow coefficients due to convection and dispersion respectively.Out of these coefficients,  is calculated from Manning's  such as for main channel        and for floodplain        .where,   ,   and   ,   are the corresponding manning's roughness and hydraulic radius.The eddy viscosity coefficients due to the bottom friction  for both sub-sections are taken as suggested by Abril and Knight [12].For main channel,   is considered which is generally a constant for all channels and for floodplain, the expression     −   is taken.Significance of the rest three coefficients and their applicability is discussed below.

Result and discussion
To verify the proposed model, homogenous and heterogeneous non-prismatic compound channels of previous investigators are taken.Converging channels of Bousmar and Zech [4] and Rezaei [7]; and diverging channels of Bousmar et al. [13] and Yonesi et al. [14] (both homogenous and heteregeneous channels) are considered.Table 1 shows summary of these data sets, where the codes such as /  / and /  / are used to explain the experimental data sets.Note that  and  stand for converging and diverging length in m respectively,   is the relative flow depth and  is the discharge in m 3 /s.According to the availability of the data sets, measured depth averaged velocity and boundary shear stress outputs at beginning, middle and end sections of the non-prismatic portions of the channels are taken.
Here, the energy loss (∆) between two sections is calculated and dividing the length with the ∆, the energy slope   is estimated.Utilizing the equation 5 and 6, the flow depth () of a desired section is determined.Then, model results of flow variables are compared with the measured values.Fig. 4(a) demonstrats the results of   and corresponding   at the middle section of the converging part of the Cv6/0.3/0.012channel.The   and   results for Rezaei [7] series are shown in Fig. 4 (b) and (c).It indicates that the narrowing floodplains compel the flow from floodplains to the main channel along the converging reach.So, maximum values of the flow variables occur at the end section ( * =6m).  results for diverging series of Bousmar et al. [13] are plotted in Fig. 4(d).Good agreements are obtained for both the main channel and floodplain regions.Fig. 4(e) to (h) present the results for homogenous and heterogeneous diverging data sets of Yonesi et al. [14].Bousmar and Zech (2004), (2) Rezaei ( 2006), (3)  In non-prismatic channel flows, a stronger mass transfer occurs between adjacent subsections due to its geometry and restrains the development of helical secondary current on the floodplains; so dispersion term is expected to be very less [4].Same observation is also found in the present study and a constant value of  equal to 0.0001 is taken for both sub-sections.Then, eddy viscosity coefficient due to horizontal coherent structure ( ) and the convection coefficient () are left to be analysed.The parameter () is explicitly linked to the non-prismatic channel geometry [4].In present study also, values of  and  are found to be sensitive to the model prediction and are varying in both subsections.This clearly confirms the significance of the transverse shear stress due to coherent structure and additional loss due to convection.These coefficients values for all experimental channels found from the proposed model are given in Table 1.For a channel, eddy viscosity coefficient in main channel (  ) is found to be varying at different longitudinal reach, but at the same time in the floodplain,   shows a single value for all the reaches.Lower values of   clarify the importance of Reynolds shear at the interface zone than near wall regions.Varying   values show that the turbulent stress at interface zone is varying along the flow direction.It can be seen that, the convection coefficient in main channel (  ) is higher than its value in floodplain (  ).It is caused due to the non-parallel streamlines in non-prismatic portions and the transverse velocity is of high order in main channel.In heterogeneous channel (Dv2/0.25/0.044), the   value is one order greater than the   value at the end position of the prismatic part.This confirms the occurrence of the higher turbulent shear stress at the interface region due to heterogeneous roughness.For data sets of Yonesi et al. [14] as shown in Fig. 4 (e), (f), (g) and (h), predictions are satisfactory for   but some deviations are found for   results in the main-channel zone.This deviation is pointed out due to the variation of local energy slope or two-dimensional energy slope across the section [4].Although    is very difficult to assss, however, it should be estimated for accurate estimation of velocity and boundary shear stres instead of one dimensional energy slope    , [4].As the velocity changes across a section, consequently the acceleration and pressure terms would vary with the velocity.So the energy slope   can not be constant for the whole section.Fig. 5 (a) and Fig. 5 (b) show the lateral variation of local energy slope    for homogenous and heterogeneous channels of [4].Values of    are estimated from the model results of   .It can be seen that,    is larger than the one-dimensional    in both the sub-sections.These values are maximum near the wall and at the interface as a consequence of large velocity gradient.So a better adjustment in prediction of   can be possible by including    .However, the practical application utilizing    can be difficult as it needs more field data sets for its modelling.As, the proposed solution provides good results for both velocity and bed shear for several non-prismatic series, so it is expected to predict both these flow variables for field cases.

Conclusions
Analytical solutions to the Navier Stokes equation for constant and variable flow depth domains are formulated and their applications to non-prismatic compound channels are analysed.The proposed approach requires an extra parameter determination i.e., the energy slope.By incorporating the energy slope and additional calibrating coefficients due to horizontal coherent structure, dispersion and convection, the depth averaged velocity and boundary shear stress profiles for homogenous and heterogeneous compound channels are well predicted.The eddy viscosity coefficient due to horizontal coherent structure and convection coefficient are found to be greatly influencing the model outputs.But, the effect of dispersion coefficient is very less.Individual calibrating coefficients for each subsection are assumed.Further, the significance of lateral variation of energy slope is discussed.

Fig. 1 .
Fig. 1.Top view of compound channels with (a) converging and (b) diverging geometries Fig. 2 shows half of a compound channel with four distinguished zones such as constant flow depth domain in main channel (Panel 1), varying depth domain in shear layer region (Panel 2), constant depth domain in shear layer region (Panel 3) and constant depth domain in flood plain (Panel 4).Due to momentum transfer between the subsections, stream wise velocity is dominant at the Panel 1 and 4, however, in the shear layer (Panel 2 and 3), both the stream wise and lateral velocity components are greatly affecting the flow variables.

Fig. 2 .
Fig. 2. Sketch of a half compound channel divided by four different zones The Navier-Stokes stream wise momentum equation is represented by and the pressure term   is known as the energy slope and represented as    − for steady and non-uniform flow, the equation 2 becomes , −,    and    −   .However, the equation (8) can not be applied for linear side slope domain (1 : ) as the flow depth () is a function of lateral dimension  (Fig.3).So, the varying depth of flow ∅ can be expressed as ∅  − − ′  .where,  ′ is the semi width of main channel and  is the side slope.

Fig 4 .
Fig 4. Depth averaged velocity and boundary shear stress distribution, measured and computed from proposed model for experimental data sets

Fig 5 .
Fig 5. Lateral variation of energy slope across half of the section for channels of Yonesi et al. (2013) modeled  as a function of sub-sectional velocity difference (  −   as      −    .So,  is expressed as a function of flow depth () such as    .where,   = coefficient of shear layer width.So, Reynolds stress term considering both contributions is      is the friction factor as  *    .So, using three expressions of the effective stresses, extended equation for estimating   for constant flow depth domain is written as

Table 1 .
Details of the experimental data sets