Multiple solutions of the Navier-Stokes equations computing water flow in sand traps

Two cases are presented wherein the main flow pattern in sand traps changes considerably as a function of discretization scheme, grid resolution and turbulence model. Both cases involve channel flows directed into a desilting basin, where the main current changes from one part of the geometry to the other. The CFD computations are validated with field or laboratory measurements. The first case presented is one of the sand traps of Khimti hydro power plant in Nepal. According to the laboratory measurements, the recirculation zone for this case is close to the bed, with the main current following the water surface. This is reproduced by the numerical model when using a first-order upwind scheme. Using a second-order upwind scheme, the main current is close to the bed, and the recirculation is formed at the surface. The second case is one of the sand traps of Tonstad hydro power plant in Norway. CFD computations predict the main flow field to follow the right or the left sides or the centre of the expansion region, depending on the discretization scheme, grid resolution and turbulence model. Field measurements show that the main current follows the centre of the expansion zone.


Introduction
Since its early emergence in the field of aerodynamics (Hess & Smith, 1967), computational fluid dynamics (CFD) is currently employed to predict fluid flow characteristics in a diverse range of engineering applications. Real-life case studies using CFD have recently been performed in combustion engineering (Akbarian et al., 2018), thermal engineering (Ramezanizadeh, Nazari, Ahmadi, & wing, 2018) and also within numerous applications of water engineering. Numerical models to predict the amount of pollutants in rivers were developed by Jiang (2002, 2004) and used for the analyses of the Pearl River Estuary. In hydraulic engineering, CFD is considered an effective technique for computation of water and sediment flow in sand traps or desilting basins. The physical modeling alternative can be problematic with regards to the induced scale effects, particularly for scaling down the sediment particle size from prototype to model scale. The numerical model then has to compute both the flow pattern and the concentration of the sediments. This has been carried out successfully both for suspended particle movements (Olsen & Skoglund, 1994;Ruether, Singh, Olsen, & Atkinson, 2005) and for computation of bed elevation CONTACT Silje K. Almeland silje.k.almeland@ntnu.no changes (Esmaeili et al., 2017;Ruether & Olsen, 2006;Török, Baranya, & Rüther, 2017). Furthermore, Olsen and Kjellesvig (1999) computed bed elevation changes in a sand trap with satisfactory results. Even if successful CFD calculations have been achieved for sand traps, the geometry of the sand trap has an expansion zone that introduces a complex flow pattern which is challenging to reproduce numerically. This article presents two case studies where multiple flow fields are identified for the flow field downstream the expansion zone. The solution of the Navier-Stokes equations depends on several input parameters such as the overall geometry of the settling basin, its wall roughness, the discharge, to name a few. A number of solution algorithms are also available for the different terms in the Navier-Stokes equations. Diverse range of turbulence models exist, and it is possible to choose among several discretization schemes for the convective and transient terms. Furthermore, disparate variants of grid and cell configurations are available.
The grid resolution or cell density is often an important parameter in deciding the accuracy of the computed flow gradients. It is commonly acknowledged that variations in the numerical algorithms can give different accuracy of the results, meaning small variations in the solution appear. This article focuses on some rare instances where the main flow pattern of the solution changes significantly as a function of the choice of numerical algorithm. Such cases have also been observed in other fields of fluid mechanics. Durrani, Cook, and McGuirk (2015) computed thermally buoyancy driven flow in a ventilation system, where a box was filled with fluid and a thermal element was placed in the bottom of the box. Perforations in the box were opened/closed to model inflow and outflow of fluid. Three different steady-state flow fields were obtained, all of them were observed in experimental results. Kamenetskiy et al. (2014) computed flow over a wing of an air plane during stall, wherein both the Spalart-Allmaras and the k − ω turbulence models were used. A numerical technique termed implicit residual smoothing (IRS) (Jameson & Baker, 1983) was used to find different flow solutions. The k − ω model was found more robust to generate multiple solutions than the oneequation Spalart-Allmaras model. However, two solutions were also found using the k − ω model. Xu, Lin, and Si (2014) obtained multiple solutions for the Navier-Stokes equations when solved for an unsteady, laminar, incompressible flow in a porous expanding channel, maintaining constant the wall suction Reynolds number and the expansion ratio. Robinson (1976) found that three numerical solutions exist for laminar, incompressible, steady flow in a parallel plate porous channel with uniform suction at both walls, when the wall suction Reynolds number exceeded a certain value. Up to three different solutions were found for a single suction Reynolds number as a function of varying skin friction. Also, two solutions were found by an analytical approach. This verifies that the Navier-Stokes equations can have multiple solutions irrespective of the numerical method in use. Kantoush, Bollaert, and Schleiss (2008) performed a series of numerical and laboratory experiments on sediment deposition in a rectangular basin. An objective was to test the sensitivity of different flow and sediments parameters and different turbulence closure schemes. In the physical model experiments, deposition in the basin systematically developed along the left bank, although the inflow and outflow were positioned symmetrically along the centre line of the basin. Although asymmetric patterns were encountered most frequently, symmetrical behavioral patterns were also observed from time to time. This behavior was also identified in the results from the numerical model. The simulations generally produced an asymmetric flow pattern that easily switched side according to the assumptions made for the initial-and boundary conditions. The results were similarly sensitive to the choice of turbulence model.
In addition, performing three-dimensional simulations on a shallow basin, Esmaeili, Sumi, Kantoush, Haun, and Rüther (2016) found the symmetric behavior of the flow field to be sensitive to small disturbances in the boundary conditions. Viroulet et al. (2017) observed two different steady-state regimes for granular flow over a smooth two-dimensional bump in a small scale physical modeling study. Dependent on the initial number of particles placed in front of the bump, either the formation of a detached jet downstream or a shock upstream the bump was identified.
In the present study, multiple solutions were found for the flow field of desilting basins in two different hydro power plants. The numerical solution was investigated for different grids, discretization schemes, and turbulence models.

Numerical models
Two different CFD programs were utilised in the present study, SSIIM 1 and OpenFOAM. The SSIIM 1 program is a freeware, while OpenFOAM is an open source program. Both programs use a finite volume method to solve the Reynolds averaged Navier-Stokes equations for steady-and incompressible fluid flow. The system of equations that was solved, is given in Equation (1)-(2): where U [m/s] is the fluid velocity, ρ [kg/m 3 ] the fluid density, p the pressure, g [m/s 2 ] the gravitational acceleration, and μ [kg/m s 2 ] is the total dynamic viscosity. The SIMPLE method (Patankar & Spalding, 1972) was used to solve the pressure field. The programs used different types of three-dimensional grids. The SSIIM 1 program used a structured non-orthogonal grid, where geometry details were modeled by blocking out cells. The Open-FOAM program used orthogonal, unstructured grids, based on hexahedral cells. For both cases, sensitivity analyses for grid resolution and discretization scheme were performed. The turbulence was modeled by RANS. All simulations for the Khimti case used the standard k − model (Launder & Spalding, 1974) for turbulence modeling. For the Tonstad case, the standard k − model, the realizable k − model (Shih, Liou, Shabbir, Yang, & Zhu, 1995) and the RNG k − model (Yakhot, Orszag, Thangam, Gatski, & Speziale, 1992) were tested in a sensitivity analysis. Further, the near-wall behavior was modeled with wall functions. A brief description of the wall models used in these simulations is given in Section 2.1-2.2 of this article.

Discretization scheme
For both cases, the convective terms in the Navier-Stokes equations were discretized with two different numerical approaches, a first-and a second-order upwind scheme. The second-order scheme used for the Khimti case was not bounded, whilst for the Tonstad simulations, a bounded scheme was applied. Regarding the turbulence variables, the first-order upwind scheme was chosen due to its favorable stability properties.
The behavior of the flow field for the Khimti case seemed to be highly dependent on the discretization scheme in use. This is illustrated and presented in Section 3.1 of this article.

Wall function
The near-wall behavior was computed by wall functions. The near-wall region consists of three main parts, the laminar sub-layer, the buffer layer and the logarithmic layer. In the laminar sub-layer, the laminar law is valid (u + = y + ) 1 , and in the logarithmic layer, the logarithmic law (3) holds, Here κ ≈ 0.4, is the von Karman's constant and the additive constant B ≈ 5.0 − 5.4 (Schlichting, 1979). These equations have previously been developed and validated for smooth walls. Experiments for rough surfaces (Nikuradse & Nikuradse, 1933) have indicated that the behavior of the flow field follows the same logarithmic slope as for smooth walls. Nevertheless, for a rough wall, the curve of the logarithmic line in the u + vs y + diagram is shifted in the negative u + -direction by a magnitude of B. Then the logarithmic law for a rough wall is given as in Equation (4): where B is a function of the dimensionless sand-grain roughness k + s = u t k s /ν t and a roughness constant C s . y + and u + used in Equation (4) and (3) represent the dimensionless distance from the wall and the dimensionless velocity, respectively. These variables are defined in Equation (5): where u τ = τ wall ρ is the shear velocity, u is the velocity parallel to the wall and y is the wall normal distance (Versteeg, 2007).
All simulations for the Khimti case were run with a rough wall function approach. In SSIIM the near-wall behaviour for rough walls is based on Schlichting wall law (Schlichting, 1979), given in Equation (6): where k s is a roughness parameter. For the simulations on Tonstad sand trap, a rough wall function was adopted. In OpenFOAM, the choice of wall function is specified through the turbulent viscosity ν t . The rough wall function, implemented in an OpenFOAM case as nutkRoughWallFunction, utilises a roughness equation in the form of Equation (4) to account for the roughness effects (OpenFOAM the openfoam foundation, n.d.). This is equivalent to manipulating the parameter E in Equation (3).

The khimiti sand trap
The Khimti-I hydro power plant is located on the Khimti River in the Koshi basin of Eastern Nepal. This run-ofthe river hydro power project is equipped with a sand trap, as the sediment concentrations in the Khimti river can be substantial. The sand trap is concrete lined and operates under open channel flow conditions. It has two basins, where only one (the left) was investigated in the study. The length of the sand trap is 135 m, including the inlet section. Considering the fact that the inlet section is slightly skewed compared to the direction of the sand trap, a divide wall was constructed to establish a more uniform flow profile in the transverse direction. The maximum channel width is 12 m and the maximum depth is 7.7 m. The sand trap was investigated in a physical model study at the Tribhuvan University in Nepal, where velocity profiles of water flow within the sand trap were measured by anchored floats (Hydroconsult, 1997). The construction material of the physical model was plywood supported by steel frames. Transparent acrylic glass and extruded foam were used within the transition bends at the upstream end of the model. The physical model was built at a geometric scale ratio of 1:15 and the other hydraulic parameters were scaled according to Froude's law. The discharge into the prototype sand trap was 15.05 m 3 /s. The vertical profile of the water velocity was measured for different vertical lines at different locations downstream of the sand trap inlet. Flow anchors at three different water depths were used for these measurements. The floats were timed over a length of 1 m.
In the present study, the flow field of the sand trap was predicted by a numerical model, where the flowing water was treated as an incompressible fluid. The sand  trap was modeled in full size. The measurements from the physical model study were used to validate the results from the numerical model. A plan view of the computational domain is given in Figure 1. Vertical profiles of streamwise velocities were measured at a location 67.5 m downstream the inlet (see Figure 1). As illustrated in Figure 2, several vertical profiles were measured along this cross-section. These locations correspond to 1.5, 3, 6, 9 and 10.5 m from the upper wall in Figure 1.
The water discharge was specified at the sand trap inlet, and a zero gradient boundary condition was given for the outlet. The initial streamwise flow velocity (u x ) in the internal part of the sand trap was estimated according to the continuity equation. The remaining velocities (u y and u z ) were initially set to zero. The near-wall behaviour was calculated by wall functions (Equation (6)). A roughness height of k s = 0.0017 m was used for the concrete walls. This is similar to a Manning-Strickler value of 90 (Mayer-Peter & Mueller, 1948;Rijn, 1982), which is tabulated as a typical value for cement lined channels (Elger, Williams, Crowe, & Roberson, 2013).

Effect of discretization scheme
Longitudinal sections of the sand trap (side view), showing the velocity fields predicted by the different discretization schemes, are given in Figure 3. These figures show that the predicted flow field varied significantly depending on the discretization scheme in use. Figure 3(a) illustrates that the first-order upwind scheme predicts a jet at the water surface, and a recirculation zone towards the bed. On the other hand, the second-order upwind scheme computes the highest velocities close to the bed. Here a recirculation zone is found towards the free surface (see Figure 3(b)).
Measured and computed velocity values for the vertical lines (a), (b), and (e) in Figure 2 are compared in Figure 4. According to the measurements, the first-order upwind scheme gives a more accurate velocity profile than the second-order upwind scheme for all vertical lines (see Figure 4). The best correlation with the measurements are found for the vertical profile 3 m from the left edge of the sand trap (see Figure 4(b)). For the lines 1.5 and 9 m from the left edge of the sand trap, the correlation between the laboratory measurements and the first-order upwind simulations are not as good as for the profiles 3 m from the left edge. Nevertheless, the results from the calculations using the first-order scheme still appear to produce solutions closest to the measured values (see 4(a),4(c)).

Grid sensitivity
Three different grid configurations were used to carry out the computations. Identical boundary conditions were used for all grids, but they differed in cell density. The coarsest, middle and finest grid contained 38,000 cells, 303,000 cells and 1.2 million cells, respectively. The simulations using the different grid resolutions resulted in similar flow fields (see Figure 5). Nevertheless, the velocity measurements for the first-order upwind scheme tended to have a better correspondence with the measured data as the grid resolution was increased (see Figure 5(a)). Further grid refinement was problematic as SSIIM 1 was not parallelized with message passing interface (MPI) and therefore not able to fully utilise the advantages of high-performance computing (HPC) clusters.

Discussion
The fact that the first-order discretization scheme predicts a more correct velocity field than a second-order scheme is remarkable as most guidelines for CFD computations recommend the usage of a second-order scheme instead of a first-order scheme (Franke, Hellsten, Schlunzen, & Carissimo, 2011). Note that the second-order upwind scheme used in these simulations is not bounded. Therefore, it may produce unrealistic overshoots/undershoots in its extrapolation of a variable value. This provides a possible explanation for the findings in Section 3.1 of this article. Aryal and Olsen (2001) found results with similar overshooting problems for the second-order upwind scheme when doing simulation at the same sand trap.
Using the first-order upwind scheme, a recirculation zone was located close to the bed of the sand trap, whilst for the second-order upwind scheme, the recirculation zone was located close to the free surface. Regarding the sand trap efficiency, a recirculation zone close to the bottom could swirl up sediments from the bottom of the sand trap, and work against the sediment settling. This illustrates the importance of an accurate prediction of the flow pattern in the sand trap with respect to sand trap efficiency calculations. Similar observations with multiple stable flow configurations have been identified in physical model studies conducted for shallow water flows (Kantoush et al., 2008;Viroulet et al., 2017).

The Tonstad sand trap
The Tonstad hydro power plant is located in southwestern Norway. In terms of electricity production, it is the largest hydro power plant in Norway. Prior to being directed into the turbines for power generation, the water enters one of the three parallel sand traps. The flow within the sand traps is pressurised. The sand traps are rock blasted and unlined, resulting in large roughness elements with complex geometry. The flow in one of these sand traps was computed by Bråtveit and Olsen (2015) using the commercial program STAR-CCM+. Lateral profiles of horizontal velocities were measured in the sand trap using field Acoustic Doppler Current Profilers (ADCP) as described by Bråtveit and Olsen (2015). The numerical results obtained in connection to this study showed variability in the placement of the maximum velocity at these lines. The inlet jet tended to follow paths on different sides of the tunnel depending on the grid type and grid resolution.
The fact that multiple solutions can be produced by different grids for the numerical model is an important issue to consider when dealing with CFD. The intention of the present study was to investigate this phenomenon and to see which parameters influenced the results. From field measurements performed by Bråtveit and Olsen (2015), velocity measurements were available for a horizontal line at a cross-section just downstream the expansion zone (cf. Figure 6). The vertical position of   the ADCP is illustrated in Figure 7. These measurements were employed for validation of the numerical model. The tunnel section was 200 m long, had a maximum height of 10 m, and a maximum width of 15 m. A generated mesh of the model is illustrated in Figure 6. A stereo-lithography (STL) file, made from a scanned point cloud of the tunnel, was available and used for the mesh generation. The tunnel was modeled in full size. The mesh was constructed by using blockMesh and snappyHexMesh, which are the meshing tools incorporated in OpenFOAM. The different ADCP devices were placed at different locations along the tunnel. The simulations were run with the simpleFoam solver in OpenFOAM. This is the Open-FOAM implementation of the SIMPLE routine (Patankar & Spalding, 1972) and contains both the standard SIMPLE version and its consistent formulation, SIM-PLEC (OpenFOAM the openfoam foundation, n.d.). In this study, the consistent version SIMPLEC was utilised.

Boundary and initial conditions
The volumetric flow rate was set constant at the inlet, with a magnitude of 75 m 3 /s. This value has been estimated by previous efficiency analyses of the turbines (Bråtveit & Olsen, 2015). The pressure flux was fixed at the inlet, whilst the turbulent variables k, , and ν t were set to initial values, in accordance with Equations (7)-(9) (Versteeg, 2007): The values for k and were based on the given discharge, and an assumption of a turbulent intensity (I) of 5%. The dissipation length scale, l, was estimated to be 10% of the width of the inlet, which was measured to 3.9 m. The turbulent model constant C μ was set to 0.09. The velocity at the outlet was governed by the pressure, which was set to the value of 6 · 10 5 Pa. This pressure was also the initial value for the pressure in the sand trap. k and were set to zero gradient at the outlet. At the walls, the velocity was set to zero, and a zero gradient boundary condition was applied for the pressure. The near-wall behavior of k, , and ν t was decided by wall functions, as described in Section 2.2. A rock-blasted tunnel has complex geometry with several roughness scales. When the grid is made finer, more roughness details are resolved. The outer boundary of the grid will therefore change according to the grid size. The sand grain roughness height k s should account for the roughness that will not be captured by the grid. A requirement stated in the CFD literature is that k s < x 2 (Blocken, Stathopoulos, & Carmeliet, 2007), where x in this case refers to the cell size of the cell adjacent to the surface. Based on the range of the near-wall grid sizes ( x = 0.03755 − 0.0625) used in the simulations, the sand grain roughness height, k s , was set to 0.01 m. The roughness constant, C s , was set to 0.5 (ref. Equation (4)). This is given as a default value in the literature (Blocken et al., 2007). The flow within the sand trap was considered incompressible. Initially, the tunnel was filled with stagnant water, with a pressure of 6 · 10 5 Pa.

Discretization scheme
The simulation outputs from the Khimti sand trap model predicted a total change in the flow pattern, strongly dependent on the order of the discretization scheme used. Based on this, a remarkable point for the analysis on Tonstad was to investigate whether the same would apply to this sand trap. The first-and secondorder upwind schemes were used to test this hypothesis. The effect of the discretization scheme was tested for different grid sizes. The results from these simulations are depicted in Figures 8 and 9. Figure 8 illustrates the results obtained from simulations on the grid with 19 million cells. From these (Figures 8 (a and c)), it is evident that the jet converges towards different sides of the tunnel depending on the discretization scheme, both for the standard k − -and the RNG k − models. The change of the discretization scheme appears to have limited effect on the results from the realizable k − model (Figure 8(b)). Figure 9 illustrates the corresponding results using a refined grid with 50 million cells. At this grid resolution, the choice of discretization scheme does not seem to have any significant influence on the flow pattern at the ADCP. For both grid resolutions, the results generated using the realizable k − model seem to demonstrate the highest correspondence with the field measurements. The magnitude of the maximum velocity is calculated with relatively good accuracy. Nevertheless, the location of this maximum velocity is relatively skewed compared to the field measurements. This applies to all computation results presented in Figure 9. The field measurements show a relatively symmetric flow pattern. Also, the recirculation zone found in the calculations is not identified in the field measurements.
The calculated flow field of the entire tunnel crosssection at the location of the ADCP is illustrated in Figure 10. These Figures 10(a,b) are from simulations on the 19 million cells grid using the RNG k − turbulence model. In Figure 10(a,b), recirculation zones are identified in the lower right and left corner, respectively. These zones are also recognized in Figure 11, which shows a top view of the flow situation. Apart from the fact that the jet follows trajectories on different sides of the tunnel, the flow situation seems similar in the two cases. A recirculation zone is formed at a lower side wall, and the maximum velocity is located in the upper part of the tunnel. Nevertheless, the results indicate that the length of the high-velocity jet is longer for the second-order scheme than for the first-order scheme ( Figure 11).

Grid sensitivity
A basic requirement for a CFD calculation is to obtain a grid independent solution (Versteeg, 2007). As the calculation domain gets sufficiently large, and the flow attains a highly turbulent state, the cell density needed to obtain a grid independent solution might be computationally demanding. Nevertheless, the starting point in this study was to carry out simulations starting with a coarse grid, and refining the grid until a grid independent solution was obtained. The effects of grid refinement towards the wall were also examined. All simulations within this grid sensitivity analysis were calculated using the secondorder upwind discretization scheme for the convective term in Navier-Stokes equation.
The grid sensitivity analysis for the standard k − model is shown in Figures 12(a), 13 and 14. From these depictions, it is evident that the path of the jet at the ADCP is predicted at different locations depending upon the grid resolution. For the coarsest grid (3 million cells), the location of the jet is predicted to be in the upper centre of the tunnel. This location corresponds with the field measurements. Nevertheless, the magnitude of the velocity is considerably lower than in the field measurements. The maximum velocity calculated on the finer grids corresponds better with the field measurements (see Figure 12(a)). For the finest grids (40 and 50 million cells), the jet appears to converge to the left side of the tunnel. This is not in accordance to the field measurements. Increasing the number of cells from 40 to 50 million seemingly had limited effect on the result. Different from the results at the finer grid resolutions, a recirculation zone is found along the whole width of the tunnel floor at the cross-section of the ADCP for the coarsest grid resolution (see Figure 14). Unlike the results for the standard k − model, the realizable k − model predicts results that are relatively insensitive to the grid resolution (see Figure 12(b)). For all tested grid resolutions, the jet converged towards the left side of the tunnel. As the grid was refined, the path of the jet tended to be predicted further away from the middle of the tunnel and closer to the wall. This feature was also identified using the standard k − model.
In addition to the standard k − and the realizable k − model, the RNG k − turbulence model was tested. The corresponding results using this model are illustrated in Figure 12(c). Contradictory to what was observed for the other turbulence models, the jet tends to converge to the right side of the tunnel for the finer grid resolutions (40 and 50 million cells) for the RNG k − model. Nevertheless, at the coarsest grid resolution used for this turbulence model (19 million cells), the jet converges to the opposite side of the tunnel. This was another instance, where considerable variability in simulation outputs were observed, as a result of changing a numerical input parameter of the calculation.

Sensitivity of turbulence model
To investigate the effect of changing the turbulence model, several versions of the k − model were tested at different grids. Figure 15, 16 and 17 present the results from this analysis. At each grid resolution tested, one of the turbulence models predicts the jet to follow a path on the opposite side of the tunnel than the other two. At the coarsest grid resolution (19 million cells), the results using the standard k − model converges to the right side of the tunnel, while the jet shifts to the other side using the other turbulence models. For the finer grid resolutions, it is the RNG k − model that   predicts the jet to follow the right side of the tunnel, while the other models predict it to the left side. The situation is similar for the two finest grid resolutions (40 million cells and 50 million cells). These results can also be identified from the previous sections (Section 4.2.2 and Section 4.2.1).

Stationary vs. transient simulations
A working hypothesis was that the different solutions could probably result if the jet was in fact oscillating between the left and the right side of the tunnel. To investigate if this could be the case, the field data was further analysed. In Figure 18, development of the field data are plotted for a period with stable production. For the plotted period, the production was relatively stable at 300 MW. Each line on the plot is averaged over 10-min intervals. Figure 18(a) suggests that minor fluctuations in velocity peaks were observed. For the first 10 min, the velocity peak was located around 3.5 m from the tunnel wall. For the next 10 min, it was located more to the centre, and for the next period, the peak was located around 5.5 m from the tunnel wall. Further, it gradually transitions back for the next two time periods. Nevertheless, changes in position for the velocity peak were not as prominent as for the simulated data. Figure 18(b) plots the field data for another time period of stable production. From this figure, any prominent trend in position change of the jet is hard to identify. The small variations in the time averaged values might be caused by limited quality of the measurements (Bråtveit & Olsen, 2015). The variations between the time averaged values from the different time intervals are largest far away from the wall (see Figure 18(b)). Based on the findings of Bråtveit and Olsen (2015), the quality of the measurements are considered satisfactory only up to a distance of 5 m away from the ADCP apparatus. For the time period plotted in 18(b), the correlation between the time averaged measurements for the different time intervals, appears to be relatively high the first 5 m from the ADCP.

Discussion Tonstad sand trap
The observed results demonstrate that there exist multiple solutions of the Navier-Stokes equations for the flow field within this sand trap. This corresponds to findings by Bråtveit and Olsen (2015). All simulations appear to converge to steady solutions, but different computational setups predict different trajectories for the jet at the expansion zone of the tunnel. Using identical boundary conditions, the jet converges to solutions on either sides of the tunnel, and in the centre, depending on the discretization scheme, grid resolution, and the turbulence  Kantoush et al. (2008). In a numerical experiment of sedimentation settling, the main flow tended to follow one of the sides of the basin. By applying small changes to the assumptions in the initial-and boundary conditions, the main flow switch to follow the other side of the basin. This behavior was further confirmed by laboratory experiments.
The role of the sand trap in a hydro power tunnel is to reduce the size of the sediments that are transported with the flow. The main purpose of this is to avoid severe damage to the turbine blades. This is done by decreasing the velocity of the water flow, which in turn allows for sediments to settle. To reduce the velocity of the inlet flow, many sand traps are equipped with an expansion zone. Increasing the width of the tunnel will reduce the flow velocity, which increases the amount of sediment settling. In addition, the expansion will introduce turbulence and recirculation zones and create a complicated flow pattern. Recirculation zones and turbulence will lead to swirling of the sediments, and thereby work against the process of sediment settling. This is very unfavorable for effective hydraulic performance of the sand trap (Brox, 2016). In addition to the magnitude of the velocity, important factors influencing the trap efficiency will then be the length and extension of the high-velocity inlet jet, as well as the placement of the recirculation zone. In this study, the path of the inlet jet and the placement of the recirculation zone were changed dependent of the set up of discretization scheme, grid resolution and turbulence model. In many sand traps, a series of tranquilizing racks are placed in the expansion zone to enhance the flow conditions with regards to sediment settling. The purpose of these racks is to homogenise the flow and thus reduce the turbulent velocity fluctuations (Paschmann, Fernandes, Vetsch, & Boes, 2017) in turbulent regions of the expansion zone. Nevertheless, construction of such racks has not been common in Norwegian hydro power plants and was not present in the expansion section of Tonstad sand trap. For this sand trap, a high-velocity jet downstream of the expansion zone was identified both from the field measurements and in the results from the numerical simulations.
Changing the discretization scheme from first-to second-order upwind moved the jet from the left-to the right-hand side of the tunnel, when calculated on the grid of 19 million cells. These results were found both for the standard k − model and the RNG k − model. Apart from the fact that the jet followed trajectories on either sides of the tunnel, the flow patterns were similar for the different discretization schemes tested. For both schemes, a high-velocity jet followed one side of the tunnel, and a recirculation zone was found at the opposite side. Whether the jet follows the left or right side of the tunnel is of little importance for the sediment settling, as long as the remaining flow features are similar. The length of the high-velocity jet seemed to be larger using the more accurate second-order discretization scheme. This could affect the efficiency of the sediment settling (Nøvik, Dudhraj, Olsen, Bishwakarma, & Lia, 2014).
For the finer grid resolutions, no significant differences were detected between the different discretization schemes. More equal results for first-and secondorder discretization schemes are expected as the grid resolution increases (Versteeg, 2007). This is due to the fact that numerical diffusion, which might be an issue for first-order discretization schemes, is reduced (Versteeg, 2007). In other words, numerical diffusion will introduce inaccuracy to the calculations, but should in principle not lead to results that are totally different as seen from this study.
Changing the grid from relatively coarse (3 million cells) to finer (19, 40 and 50 million cells) made the jet follow the middle, right and left side of the tunnel, when using the standard k − model for turbulence modeling. Similar to the sensitivity analysis on discretization scheme, the flow pattern was comparable, but on different sides of the tunnel. This means that with regard to sand trap efficiency, this should not be of significant importance. Nevertheless, for the coarsest grid resolution, the jet was predicted to follow the upper centre of the tunnel and a recirculation zone was spanning the whole width of the ground section of the tunnel. This might affect the efficiency of the sand trap.
All converged solutions considered, it is evident that in the majority of the cases, the jet follows the left wall of the tunnel. This is not according to what is seen from the field data, where the jet seems to follow the middle of the tunnel. Nevertheless, similar behavior was identified by (Kantoush et al., 2008) in a physical model. For the majority of the experiments, the flow was aligned along one of the sides of the basin, and occasionally it changed to the other side.
To investigate whether the multiple solutions found in the simulations could be due to some transient features of the flow, the field measurements were analysed to consider whether the jet was in fact oscillating between the left and right side of the tunnel. Any prominent oscillations could not be proved to exist based on the available field measurements. The fact that velocity measurements are only available for one horizontal line at three crosssections of the tunnel, and also that questions are raised concerning the validity of the measurements data for the right side of the tunnel (Bråtveit and Olsen (2015), limits the ability to draw strict conclusions regarding the behaviour of the flow and the accuracy of the flow field calculations.
To be able to draw further conclusions on this, field measurements covering more features of the flow would be required. Moreover, transient simulations using detached eddy simulation (DES) or large eddy simulation (LES) could be performed. In these models, the motions of the turbulent eddies are captured, which gives an improved capability to capture transient features within the flow. Nevertheless, the use of LES in flows with Reynolds number in the range of 10 6 , as in this case, has shown to be challenging (Catalano, Wang, Iaccarino, & Moin, 2003). The use of LES might have greater applicability replicating physical model studies, where the Reynolds number will be smaller.
Based on the analysis of the field data and the magnitude of the Reynolds number, simulations based on LES were not performed in the current study. Durrani et al. (2015), Kamenetskiy et al. (2014) and Robinson (1976) reported multiple solutions when solving the steady Navier-Stokes equations. The current study introduces new examples where the solutions of the Navier-Stokes equations solved by CFD are not unique. The findings from the mentioned studies suggest that interpretation of results from CFD computations should be treated with caution. The results should be validated by physical measurements, and sensitivity analyses for grid resolution, turbulence model and discretization scheme should be carried out.

Conclusion
The results from CFD simulations done at Khimtiand Tonstad hydro power plants, failed to produce a unique solution for a flow field within a sand trap, using a RANS model. The main reason for this is that the expansion zone, which is the crucial part of the sand trap, creates a complicated flow pattern where a high turbulent zone is created. None of the sand traps included in this study were equipped with tranquilizing racks, which could have had a beneficial impact on the flow pattern by homogenising the flow.
According to the current study findings, changes in input parameters as discretization scheme, grid resolution and turbulence model could have significant impacts on the computed flow pattern within sand traps. If CFD computations are performed with the aim of deciding the settling efficiency, discrepancies could arise due to changes in the numerical parameters investigated in this study. Nevertheless, major part of the sediment settling takes place downstream the expansion zone. Here, the fluid velocity is reduced, and a more uniform flow field exists.
The study of the Khimti sand trap suggests that the boundedness of the discretization scheme might be of higher importance than its order, with regards to accurately predicting the flow pattern within a sand trap.
To further validate the results from these case studies, transient analyses using DES could be performed. Field measurements covering more features of the flow would be required to be able to further validate the CFD calculations.