Numerical Modeling of Local Scour at the Junction of Open Channels in Flow 3 D Numerical Model

At the junction of channels, the two corresponding flows of the main and submain channels are diverted from their main alignment and the form and the flow properties change at the junction. Changes in water level profile and depth of flow, velocity distribution, stagnation zone, constriction of public channel, energy loss and also formation of hydraulic jump are among the most important hydraulic variables in this location. For accurate recognition of hydraulic properties of flow and local scour at the junction of channels, physical models are made and constructed. Setting up a physical model requires many conditions and high costs which sometimes are not justifiable, hence appropriate numerical models could be proposed for such options. In this research using Flow3D numerical model, the numerical modelling of the flow has been performed in 3D form utilizing the available laboratory information which is calibrated and validated and accuracy of the numerical modelling, and the corresponding relative error are determined. The calibration and validation of the numerical model results demonstrate that the maximum relative error of the numerical model when simulating for maximum values of scour depth at the flow junction is equal to 8.2%. Also using the numerical model it was found that with passage of time in numerical model, from the start of scouring, the location of maximum scour is transferred towards the opposite wall of the sub main channel and is distanced from the junction position also the volume of sedimentation is increased and is translated toward the downstream main channel.


Introduction
At the channels junction, the two main and sub main channel flows are deflected from their main route and the form and properties of the flow change at the junction.Change in water surface profile and depth of the flow, velocity distribution, stagnation zone, constriction of the public channel, energy loss and also formation of hydraulic jump are among the most important hydraulic variables in this location.At the stagnation zone which is formed at the corner of junction upstream, the first encounter between two flows corresponding to the two main and sub main channels occurs and the flow velocity at this zone is nearly zero.The deflection zone is a zone in which the main channel flow is deflected from its rout and approaches the opposite wall of the junction.The sub main channel flow after the junction distances from the inner wall of the channel and creates a zone called the separation zone, in which the flow velocity is very low and sedimentation occurs at this zone.At this zone the maximum velocity or the compression zone of flow velocity is greatly increased and when the shear stress exceeds the critical stress, local scour occurs at this zone.Estimation of the local scour at the junction of rivers and channels is one of the most important problems in Hydraulic engineering.
Mosely, performed an experiment at the junction of two Y-shaped channels with erodible walls and demonstrated that erosion starts by creating sand mounds perpendicular to the flow direction at the sub main channel which in continuation meanders in downstream channel parallel to the flow direction and creates the erosion pit at the junction point and with increase in junction angle from 15 to 75 degrees, the depth of erosion pit is increased intensely [1].Best, studying the location of a junction with different junction angles and discharge ratios found that the maximum scour depth is a function of the encounter angle between two flows and demonstrated that with increase in junction angle, the deflection angle of the maximum scour depth increases from the junction upstream edge [2].Nazari in his laboratory studies of a 90-degree junction, utilizing an empirical relationship, demonstrated that by increase in mean diameter of the bed materials, the scour depth decreases and by increase in ratio of the main channel discharge to that of the sub main channel and decrease in the ratio of sub main channel width, the scour depth increases [3].

Materials and Methods
As numerical modeling of pure water flow field, bed scour and suspended sediments in channels and their reservoirs is of vital importance in presenting the design solutions, also experience of the researches concerning these kinds of simulations has shown that among the current software packages, the Flow3D software has a better capability in modeling of this type of hydraulic conditions, therefore in this study Flow3D (Ver.10.1) has been utilized to investigate the local scour at the junction of open channels.The numerical modeling in this software is performed based on the experimental studies.In continuation of this section, the geometry of the laboratory models is built by creation of identical numerical modeling conditions similar to the built hydraulic model and then the flow field hydraulic parameters are investigated and compared in these conditions.In this study according to the laboratory conditions, an intersecting channel with erodible bed is selected similar to the 3D simulation of the flow field laboratory studies and then in continuation the geometric properties of the channels and their 3D geometry are presented.The model which was investigated in this study was drawn in 3D form in SolidWork software and a schematic view of the intersecting channels with their erodible bed is shown in figure below in which the layout of weirs is also shown.

FLOW 3D and Governing Equations
FLOW 3D is an appropriate model for complex fluids problems.This numerical model is widely used, particularly for unsteady 3-dimensional flows with free level and complex geometry.In this model, finite volume method is used in regular rectangular grid generation.Due to using finite volume method in a regular grid, the form of the employed discrete equations is similar to discrete equations in finite difference method.Accordingly, Flow3D enjoys first and second-order reliability methods which are explained in the following.Also, this software uses five turbulence models such as k-ε and RNG.In Flow3D, two methods have been simultaneously used for geometrical simulation.The first method is volume of fluid (VOF) which is used to show the behavior of fluid at free level.The second method is fractional area-volume obstacle representation (FAVOR) which is used to simulate solid levels and volumes such as geometrical boundaries.
Equations governing fluid flow are obtained from the law of conservation of mass and the law of conservation of momentum.These equations are in form of partial differential equations.In general, to obtain flow equations, three steps should be considered: selecting accurate base laws, applying laws by an appropriate model and adopting mathematical equations showing the above physical laws.The main equations to simulate 3-dimensional flow are three differential equations including continuity relations and movement size in x, y and z directions.
In Flow3D, to simulate transfer, erosion, deposition, and change of sediments establishment status is due to fluid flow.Sediment model of this numerical model uses two concentration fields including suspended sediments and bed sediments.Displacement of suspended sediments with fluid is due to local pressure's gradient changes.These suspended sediments may be created due to input flow containing suspended particles or due to bed erosion.Since bed erosions have been limited by neighboring particles and are not easily displaced, they can move only in case of changing into suspended load in shared level of bed and fluid.Suspended load can be changed into load when the velocity of depositing is higher than the velocity of bed erosion.A part of control volume occupied by solid particles of sediment is ( ) and the rest is defined from accumulated fluid of ( ) such that: Suspended load causes to the increase of real fluid viscosity.This increase continues until solid particles' volumetric element reaches to volumetric cohesion element.After that, increasing suspended load does not lead to the increase of viscosity but it causes that particles start to behave in solid manner.In this state, average viscosity of fluid is computed from the following relation: Where μ f indicates fluid's viscosity and μ* indicates average viscosity of sediment particles' critical element.Apparent density of is assumed as linear function of sediments volume where ρ s and ρ L indicate apparent density of sediment and fluid, respectively.
Drift refers to the deposition of sediment particles under the impact of floating forces affecting sediment particle.In model of sediment washing in Flow3D, sediment particles are assumed in spherical shape such that it is influenced by fluid's viscosity effect.Therefore, deposition coefficient is automatically computed according to the following relation.
Therefore, deposition velocity is computed through the following relation: Where in equation above, ( ̅ ) indicates mechanical potential of gradient or acceleration and is limited to 10 times more than particle's weight and causes to omit numerical fluctuations in pressure value.Near to fluid's free level, the value of, ( ̅ ) is replaced with acceleration (g).The coefficient of f L has been employed in above equation since sedimentation is possible only at the presence of solid particles (sediment).Therefore, if control volume is full of sediments, = and then,   = At the level of bed sediments, shearing stress is active and causes to erosion and displacement of sediment at bed sediment level.This erosion is a function of fluid's sear stress at surface, critical shearing stress and sediment and fluid's density.The parameter of critical shields indicates the minimum shear stress required to lift sediment particles from the shared surface of fluid and active bed.
Where θ crit indicates critical shields; ᴦ crit indicates the minimum shear stress necessary for bed length to lift sediment particles.The purpose of developing and explaining this model is to estimate and predict the magnitude of sedimental flow which has been worn out.To this end, the parameter of shear velocity is defined to measure removal power of flow.So, the velocity of removing sediments from bed can be presented through the following equation.
Where n s indicates normal vector of bed surface; α indicates dimensionless parameter indicating the probability of sediment particles removal from bed (usually 1).In static fluid, internal friction angle of sediment particles determine the minimum slope through which sediments' walls can be steady.Internal friction angle above sediments indicates steady wall slope in steep slopes such as clay.In low angels of walls, there is a strong intention to collapse and move forwards such as sand.
In downstream hole where sediments are compiled together and create a mass of sediments, sediments establishment status makes an angle with horizon surface which indicates internal friction angel.In the model, this angle is signified by .Natural establishment angle of sediments in various temporal and spatial conditions is computed through the following relation: Where n interface equals normal vector of surface and g indicates gravity acceleration.After scour or transferring sediments suspended at the surface, critical shear stress occurring in sloped surface for each surface is computed through the following relation: According to the above equations, when natural slope of sediments equals their internal friction ( = ), critical shear stress equals zero, indicating that bed surface undergoes erosion due to every kind of imposed shearing stress.Also, when ( ), we have (  ), indicating that higher internal friction angel of sediment particles leads to higher wall slope (φ) since the wall of scour or sediment washing hole undergoes erosion without shear stress (  = ).The movement of suspended sediments in system is expressed through convection-diffusion equation.
Where C s indicates sediments' concentration, Γ indicates diffusion coefficient and ω s indicates particles' collapse velocity which equals: Therefore, the above equation is entered into problem solving as following: The concentration of sediments suspended in shared surface of sediments' bed and water before starting sediment washing (t=0) equals: The above computational equations and algorithm in Flow3D are used to discharge bed sediments.

Numerical Modeling of the Bed Local Scour in Flow3D
In Figure 2. below the built geometry in SolidWork software is added to Flow3D numerical model and utilizing the grid blocks, the computational space for numerical modeling of the local scour phenomenon at the intersection of open channels is generated.Existing boundary conditions should be introduced to the numerical model based on the accurate and feasible conditions.In Table 1.boundary conditions applied to the model are presented.

Location of boundary in numerical models
The amount of flow at the inlet channels Q Xmin

W Zmin
Symmetry boundary condition for the free surface and the boundary between the two blocks S Zmax

Figure 4. Boundary conditions applied to the numerical model
In all the performed models in this study, the 3D flow field is solved by the RNG turbulence model.The reason for utilizing this turbulence model is due to its features and advantages in comparison to models like k-ε.This model, due to possessing an extra term in ε equation, is improved for analysis of rapidly strained flows and flows over surfaces with large geometric changes.This model has a high capability in simulating the unsteady flows [4].Also according to the comparison made in studies of Taleb Biddokhti et al. (2012), between different turbulence models using Flow3D software for scouring, the RNG turbulence model is introduced as the best model of this software for scouring and on this basis the RNG turbulence model is utilized [5].Other needed parameters in the numerical model are also selected in accordance with the real conditions in this model, for example water is considered as a non-viscous, incompressible fluid and the air entrance occurs with a density of 1.2 kg/m 3 and the shear stress coefficient is 0.073.
For numerical analysis of the local scour at the junction of open channels, the results of laboratory studies performed by Borghei and Jabbari Sahebari (2010), are utilized [6].In these studies the main channel length is 12 m and that of the sub main channel is 2.5 m and the junction angle (θ) in this research is taken equal to 50 degrees.The width of both channels is 40 cm.Based on the laboratory studies in this research the longitudinal slope is ignored.In Figure 5. different parts of the model are shown.The erodible bed has a depth of 13cm and its length is 2m at upstream of the main and sub main channels, before the junction, and 2.5 downstream of the main channel after the junction.The dimensions of this zone which are depicted in the figure are chosen empirically.At the beginning of this zone in main and sub main channel, to create the laminar flow use has been made of flow smoother which is a gravel zone of 0.5 m length with a grading ranging from the course to fine grains.To eliminate the effect of non-uniformity and cohesion of the grains on the erosion pattern and sedimentation, the used bed materials have a uniform grading without cohesion.As is shown in Figure 5, the center of coordinates is taken at downstream corner of the junction and the positive direction of the x-axis is taken towards upstream flow and the positive direction of the Y-axis is taken towards inside of the main channel.

Figure 5. Applied initial conditions at the beginning of the channels[6]
The first step in a numerical model is its calibration, i.e. the effects of external factors should be reduced to minimum and the model conditions should resemble the real ones.To extract accurate and precise values for a numerical or real model, reaching the stable conditions is essential.In the studied model after investigating some models, the appropriate time for extracting results from the model was taken equal to 20 seconds.In Figure 6. the flow passage and intersection of the flows in the model could be observed.The flow becomes stable after 20 seconds in both channels and it moves outside from the exit boundary.According to the numerical model results, at the starting time of the numerical model computations (T=0.0Sec), the discharges are introduced into the model per considered laboratory data.By initiation of software computations (T=1 Sec) water flows over a portion of channel.At time (T=2 Sec), the flows encounter each other at the intersection of two channels.At time (T=5 Sec), the flow is transferred from the exit boundary to the outside and the effect of flow unsteadiness is observed in the main channel.Although at time (T=10 Sec), the flow appears to be steady but focusing on the difference between fluid fraction contours at consequent times it becomes clear that the flow from time (T=20 Sec) on becomes steady.In order to ensure stability and balance of the flow field in channels, the diagrams corresponding to variation in discharge passing through the exit boundaries are extracted from the numerical model with respect to time.The figure below indicates flow stability and uniformity after 15 seconds which verifies stability of the simulation, too.For investigating validity and calibration of the numerical model, variation in the ratio of the erosion depth to the downstream flow depth, (y d /h), in terms of variation in the sub main channel discharge to that of the main channel (Q bd ) is compared and assessed for both the numerical and laboratory cases assuming the flow depth value being constant.Comparison between the Flow3D numerical model results and the laboratory results shows that the maximum relative error for simulation of the maximum scour depth values at intersection of the flows is 8.2% which is acceptable with respect to the effective parameters used in the model and measurement of the laboratory data, hence the current numerical model is a calibrated and validated model.Also analysis of the above diagram reveals that within discharge (Q bd ) range of 0-0.30, variation in the y d /h ratio is more intense and the curve has a sharp slope and by increase in Q bd from 0.3 to 0.8, the curve slope becomes moderate.Considering the figure, it seems that in these experiments at higher discharge ratios, with increase in scour intensity, the erosion depth comes close to the flow depth and the y d /h ratio approaches unity.
The threshold of movement velocity (V c ) is the velocity at which the shear stress exerted on the bed particles equals the critical stress and particles begin to move.In this study, to find the threshold of movement velocity the relationships cited by Melville and Sutherland for the Shields diagram have been utilized [7].According to their definition the critical shear velocity ( *  ) for grains with mean diameter range (d) of 0.1-1mm is obtained using the following expression: *  = 115 + 125 1 4   (13) As the mean diameter of the particles in this modelling was 0.97 mm, the critical shear velocity is 0.033 m/s.By replacing the particles' diameter and the critical shear velocity values in formula presented by them, the threshold of movement velocity would be a function of the downstream flow depth according to the following expression: In Equation 14, y d is the downstream flow depth in meters and V c is the threshold of movement velocity in m/s.Concerning the moment at which local scour initiates at the channel junction, there are no available documented studies.Scour around bridge pier occurs for   5 but concerning the channels junction there is not a criterion at hand.As in this research the V c /V ratio has been selected as an important variable therefore some values should be chosen for this variable so that scour would occur for all the considered cases of other variables.Therefore for the junction angle of 50 degrees, the threshold to movement is investigated.In this research for four different sub main channel to main channel discharge ratios and the width ratio equal to unity, the initiation moment for the movement of particles is investigated.There is not an accurate criterion for selecting the moment of threshold to movement.In this study the initiation moment of the threshold to movement for generating the maximum erosion depth of 3 mm is selected 1800 seconds from the start of model execution.In Figure 9, variation of V/V c with respect to Q bd at the moment of the threshold to movement is drawn for four different downstream total discharge values (Q d ).In the executed models for investigating the threshold of the movement, by reaching the downstream discharge to the considered discharge, the bed particles began to move with decrease in depth and increase in flow velocity and secondary flow intensity.In these models a number of small holes are generated inside the main channel at the downstream edge of the junction, and the maximum depth of scour pit is taken as the criterion for specifying the threshold of movement.As is seen in Figure 9, with increase in ratio of the sub main channel to main channel discharge, intensities of vorticity flows and flow compression at the junction of channels are increased, and consequently erosion occurs at lower V/V c ratios.Also with decrease in the ratio of the sub main channel to main channel discharge, intensities of vorticity flows and flow compression are decreased and the flow condition approaches the condition of a flow without sub main channel.At a Q bd equal to zero, erosion takes place at the V/V c ratio equal to unity.It is interesting to note that contrary to the bridge pier, erosion begins at V/V c ratio less than 0.5 at the channels junction.

Conclusion
In this research the effects of various variables including the ratio of main channel to sub main channel discharge and the ratio of flow velocity to the threshold to movement velocity are investigated using the calibrated numerical model at a 50-degree junction.Calibration and validation of the numerical model results show that the maximum relative error for the numerical model in simulating the maximum scour depth values at intersection of flows is 8.2%.The threshold to movement occurs at lower V/V c ratios with increase in ratio of the sub main channel to main channel discharge and contrary to the bridge pier, the threshold to movement could also occur at V/V c ratios less than 0.5.By increase in ratio of the flow velocity to the threshold to movement velocity at downstream and consequently increase in the Froude number at downstream, the type two helical cells become larger and cause reduced type three cells and the erosion due to it is decreased in proximity to the main channel wall and opposite to the junction location.With increase in the ratio of the sub main channel to main channel discharge, the position of the maximum scour moves toward the junction and the erosion pit transfers toward the sub main channel.At higher discharge ratios with increase in scour amount, the erosion depth approaches the flow depth and the y d /h ratio approaches unity.

Figure 1 .
Figure 1.Simulation of 3D geometry of channels and erodible bed materials at the junction by SolidWork software

Figure 2 .Figure 3 .
Figure 2. 3D geometry at the intersection of two channels in Flow3D numerical model

Figure 6 .
Figure 6.Flow modelling in Flow3D numerical model and encounter of flows at the intersection of two channels

Figure 7 .
Figure 7. Variation in the discharge passing the outflow boundary with respect to time

Figure 8 .
Figure 8. Validation of the numerical model results with respect to the laboratory results

Figure 9 .
Figure 9. Curves depicting changes in V/V c with respect to Q bd at the threshold to movement moment