Four free surface algorithms for the 3D Navier–Stokes equations

Three-dimensional (3D) numerical models solving the Navier–Stokes equations are today increasingly applied in water engineering projects. A particular problem for rivers and channels is determining the location of the free surface (Aguilar et al. 2013). The Navier–Stokes equations require special algorithms for computing how the surface changes in space and time (Erdbrink et al. 2014). The most classical hydraulic problem is to predict the coefficient of discharge for a spillway (Savage & Johnson 2001; Jacobsen & Olsen 2010; Li et al. 2010). This is also one of the most common investigations for commercial laboratory model studies. Another hydraulic engineering problem related to computation of the free water surface is reservoir flushing. Then the water level is drawn down to increase the water velocities and thereby erode the sediments in the reservoir (Haun & Olsen 2012). This problem requires an accurate and stable prediction method for the free surface. Predicting too low water levels will lead to a higher velocity and shear stress, giving too large erosion. Too high water levels reduce the computed erosion in the reservoir.

Laboratory studies are often used to test the accuracy of numerical models (Wilson et al. 2003). The test case can then be made fairly simple, ensuring that errors in the geometry modelling will not affect the results. A good test case for free surface algorithms is a weir, where the geometry is relatively simple compared to a natural river. Laboratory measurements of the water surface over the weir have been made in the current study, enabling a base for testing the different free surface algorithms.

The most common method for computing the free water surface is the volume of fluid (VOF) approach used in a fixed grid. This method is included in all commercial general purpose 3D computational fluid dynamics (CFD) programs (Temeepattanapongsa et al. 2013; Castillo et al. 2014; Liu & Yang 2014), and also in the popular OpenFOAM open source program. The VOF method has the advantage that very complex water surfaces can be computed. Multiple surfaces can be computed, for example, flow over a sharp-crested weir (Lv et al. 2011; Abd El-Hady Rady 2011). A disadvantage with the VOF method is that a fairly dense grid has to be used near the water surface to get an accurate prediction of its location (Rüther et al. 2010; Feurich & Olsen 2012). Also, a relatively short time step has to be used, respecting the Courant criteria. The required size of the time step may then be only fractions of a second. Longer time steps produce instabilities. In sedimentation engineering, it is often desirable to predict sediment transport in a reservoir over a very long time, for example, 50–100 years. The computational time using the VOF method would then be too long.

In most water engineering cases where a river is modelled, the water surface is fairly horizontal (Tritthart & Gutknecht 2007). If this assumption is made, it is possible to derive alternative numerical algorithms for the determination of the free water surface. An advantage with this simplification is the possibility to employ an adaptive grid that follows the water surface (Olsen & Kjellesvig 1998). A coarser grid can then give the same accuracy as a finer fixed grid. Also, alternative equations can be solved, giving different accuracies and stabilities for the computation of the free surface. The most interesting approach for simulations over a long time would be implicit methods, where long time steps can be used. Two implicit and two explicit algorithms for computation of the free surface are described and tested in the current paper.

An earlier study (Haun et al. 2011) investigated different algorithms for computing the water surface over a broad-crested weir for a 2D situation. The current paper also uses a broad-crested weir, but with a more complex 3D geometry. This enabled 3D effects to be investigated. The weir is V-shaped looking both in plan view and in a cross-section. The geometry is given in Figure 1, showing a plan view (Figure 1(a)), longitudinal profile (Figure 1(b)) and a cross-section (Figure 1(c)). A physical model of this weir was built at the laboratory of the Department of Hydraulic and Environmental Engineering at the Norwegian University of Science and Technology. The weir was built of plexiglass, and was 60 cm wide, 90 cm long and 20 cm high. A photograph of the weir is shown in Figure 2. The model was installed in a flume which was 60 cm wide, 1 m high and 15 m long. A pump lifted water into the flume at a rate of 21.6 L/s. The downstream water level was kept at 15 cm above the flume bed. The physical model with water is shown in Figure 3(a) from the side and Figure 3(b) from above. The water flowed over the weir where critical flow occurred. The Froude number was, however, so low that a large hydraulic jump did not form. Since the geometry of the weir was V-shaped also in plan view, the border between the subcritical and supercritical flow followed a curve in the horizontal plane, as seen in Figure 3(a) and 3(b). The water level upstream of the weir was measured to be 25.4 cm above the flume bed.

The water level was measured from the top of the laboratory flume using a point gauge along two longitudinal profiles. One profile followed the centreline of the flume and another 3 cm from the side of the flume. The results from the laboratory measurements are shown in Figure 4. The two profiles are not identical, since the bed level was 5 cm deeper at the centreline than at the walls. The wall sections were also longer than the central section of the weir. The side of the weir reached further downstream than the central section of the weir. This situation is seen in the difference between the water levels on the sides and in the centreline in Figure 4. There were also other three-dimensional effects shown in the measurements, as lateral variation of the water surface at the top of the weir. Over the most upstream part of the weir, the water level was higher at the centre of the flume than at the sides. This reverted further downstream, before the longitudinal water level gradient became steep and ended in the small hydraulic jump.

Both computer programs used in the current study solved the Navier–Stokes equations in three dimensions using the k-epsilon turbulence model (Launder & Sharma 1974). A control volume method was used for the discretization of the equations. The pressure was computed based on the SIMPLE method (Patankar 1980) and wall laws were used to resolve the steep gradients close to the walls/bed (Schlichting 1979). The boundary conditions were a specified velocity and turbulence value at the inlet and zero gradients for all variables at the outlet. These numerical algorithms have been verified by a number of earlier studies where the results were compared with field or laboratory data (Fischer-Antze et al. 2008; Baranya et al. 2012; Ðordevic 2013).

The three-dimensional grid was based on a 2D depth-averaged grid of the spillway and the flume bed, as given in Figure 5(a). The grid had 150 cells in the longitudinal direction and 20 cells in the lateral direction. An initial horizontal water surface was specified above the bed and the 3D grid was made between these surfaces using vertical grid lines. This structured, non-orthogonal grid was made with the SSIIM 2 program and used for all the computations. The two-phase VOF method used this grid for the whole computation, and the other methods used this grid as a starting point before grid changes. The adaptive grid methods changed the grid after each time step, as the water level moved. This ensured that only the water fluid was modelled. The number of cells in the vertical direction changed during the computation, as seen in Figure 5(b). The adaptive grid method is described in detail by Olsen (2003).

VOF is the most commonly used method to predict the free surface with the solution of the 3D Navier–Stokes equations (Soleiman Beygi et al. 2014). A fixed grid with a horizontal upper boundary was used for the VOF computations in the current study. The OpenFOAM program was used to compute the two-phase flow of air and water. The most commonly used method to create grids with OpenFOAM is to place a structured, orthogonal grid around the geometry, which is defined by STL files. Then the OpenFOAM program snappyHexMesh is used to cut the orthogonal grid to its final shape. Another approach was used in the current study. A non-orthogonal grid was made (Figure 5(c)) that followed the bed of the flume and top of the weir. The grid was made by the SSIIM 2 program and exported to OpenFOAM in a blockMeshDict file. The interFoam solver was used to compute the two-phase flow of air and water. The top of the grid was then set to a level 29 cm above the bed level of the flume. The inflow region was set to cover the whole height of 29 cm. The outflow region initially only covered the lower 15 cm, corresponding to the measured downstream water level.

The VOF method requires a fairly short time step to give accurate and stable results. A time step of 0.001 s was used in the current study.

The OpenFOAM program produced a file called alpha1 for a number of time steps during the computation. The file contained the α values in each grid cell. The file was read back into SSIIM 2, where an interpolation algorithm changed the surface to correspond with α values equal to 0.5. The SSIIM 2 grid was then regenerated accordingly, and the resulting water surface written to a file. The SSIIM 2 program was thereby used as a post-processor for OpenFOAM. The resulting water elevations are shown in Figure 6, compared with the measured values.

A problem with the OpenFOAM computations was the downstream boundary condition. The downstream water level was kept at 15 cm in the physical model, and therefore the outlet area was initially set to be 15 cm above the bed in the numerical model. However, false diffusion caused the water level downstream of the weir to become higher than 15 cm. To solve the problem, the water level for the VOF computations was only allowed to flow out of the lower 10 cm of the downstream end of the grid. This caused a correct water level to form downstream of the weir during the first part of the simulations, as seen in Figure 6. However, when the simulations were done over a very long time, the downstream water level would rise. The water level over the weir and upstream of the weir was not affected by this problem.

The adaptive grid methods can be divided into two groups: one which includes gravity in the Navier–Stokes equations and one which omits this term directly. The latter methods take gravity into account when computing changes in the water surface.

The method using the continuity equation with gravity in the Navier–Stokes equations and an adaptive grid is in the following given the abbreviation CGA (continuity, gravity and adaptive grid). The algorithm was first developed by Olsen & Kjellesvig (1998), who used it successfully to compute the water flow over a spillway. A deviation of 0.5% was observed between the predicted and measured coefficient of discharge. Olsen (2013) also used the method to compute the movement of soil in a sand slide due to a water reservoir drawdown.

The time step is denoted Δt and A_z is the area of the cell projected in the horizontal plane. Since the formula is proportional to the time step, the method will require a short time step to be stable. A time step of 0.001 s was used in the current study, similar to what was used in the VOF method. The initial downstream water level was given as 0.3 m above the bed, and drawn down to 0.15 m after 2 s.

The results are given in Figure 7. The computed water surface profiles follow the measurements well over the weir. Upstream of the weir, the water level is slightly lower than the measurements. The water level downstream of the weir shows a dip at the hydraulic jump, a feature that was not observed in the physical model. Comparing the VOF and the CGA methods in Figures 6 and 7, the results are fairly similar, although the CGA approach gives slightly lower water levels upstream of the weir. The computational time was similar for the two methods.

The water density is denoted ρ, and g is the acceleration of gravity. The pressure, P, in each cell is computed from the SIMPLE method, when solving the Navier–Stokes equations. Equation (4) then gives the elevation difference between two cells as a function of the pressure difference. The equation can be used for each cell and solved implicitly. The method is abbreviated IPDA (implicit pressure differences with adaptive grid).

Equations (5)–(7) were solved implicitly using a Gauss–Seidel method. The boundary conditions were zero gradients at all sides, except for the downstream boundary. The initial downstream water level was equal to 0.3 m, the constant initial value for the whole grid. The downstream water level was drawn down to 0.15 m after the first time step and was kept fixed there during the computations.

Figure 8 shows good results for the water level behind the weir, and also reasonable results over the weir. Downstream of the weir, there is a dip in the water surface elevation where the hydraulic jump is located. The water level downstream of the jump is also too high. The reason is probably that the flow causing the energy loss over the jump is complex and therefore not computed correctly.

The method will be stable for long time steps if the weighting coefficients, a_i, are positive and the source term, S_p, in the equation is not proportional to the time step.

The shallow water equations and the water continuity equation have in principle three unknowns: the velocities in two horizontal directions and the water depth. The equations are usually solved in a coupled manner. The approach used in the current study decouples the equations, so that only one equation is solved for the water level. The other variables, the velocities, are taken from the solution of the Navier–Stokes equations.

The method was also used as boundary condition for the IPDA method.

Equation (9) was solved with a standard Gauss–Seidel solver. A time step of 100 s was used, with 100 inner iterations per time step for the implicit solution of the Navier–Stokes equations. Similar to the IPDA computations, the initial downstream water level was 0.3 m, and drawn down to 0.15 m after the first time step, where it remained fixed.

The result is shown in Figure 9. The water level upstream of the weir is computed fairly well and so is the water level in the steepest area upstream of the hydraulic jump. The maximum overprediction of the water level is 120%, at the acceleration zone on the side at the top of the weir. There are steep gradients in the water velocity at these locations, so the omission of the convective term could contribute to the discrepancy between measured and computed water levels. The largest velocity gradients are however at the hydraulic jump downstream of the weir. Downstream of the hydraulic jump, the water depth is overpredicted by more than 35%. The current grid is too coarse to resolve the flow field in this area properly, but also other factors complicate the computation of the hydraulic jump. The shallow water equation is based on a hydrostatic pressure assumption, which is not present in the hydraulic jump. There are also significant vertical velocity components upstream and in the hydraulic jump, pushing the water level down. A hydraulic jump may entrain a significant amount of air, complicating the physics. All these factors significantly reduce the accuracy of the solution of the shallow water equation for hydraulic jumps.

Note that the choice of the Manning's friction coefficient will affect the slope of the water surface. The energy losses over the weir are due to complex physics and the Manning's friction coefficient must take the losses into account. The value of the coefficient can therefore not be equal to what is usually given for a smooth wall, although the weir and side walls of the flume were made of plexiglass. The result in Figure 9 was obtained with a Manning's friction coefficient of n =0.1. Using a value of 0.2 gave a 3.2 cm higher water level behind the weir.

The OpenFOAM version used in the computations was made by the company BlueCFD. It was the single core version 2.0 for Windows. The computation was done on a 3.4 GHz PC with the Intel i7-4770 processor. The computational time was about 20 h until a steady situation was reached, as given in Figure 6. Note that most OpenFOAM versions can run on multiple cores, so on a four core PC, it is likely that the computational time could be reduced considerably.

The SSIIM computations were run on the same computer as the OpenFOAM program, but on two cores. The SSIIM program is parallelized with OpenMP, which gives good scaling for coarse grids. The computational speed was then about 40 time steps/minute, giving about 10 h of computational time until a steady solution emerged using the IDWA method. The IPDA method required much less computational time. Note that the adaptive grid will use roughly half the number of cells compared with the VOF method, reducing the computational time.

Figure 10 shows a time series of the water levels upstream of the weir for the three methods with an adaptive grid. The CGA method has an initial rise in the water levels that could come from an instability when the water hits the weir. As the time step is 0.001 s, the total simulated time in Figure 10 for the CGA method is only 20 s. The implicit methods give a faster convergence than the explicit method, with fewer time steps. The IPDA method converges much faster than the IDWA method. The reason is that the a_i coefficients in Equations (5) and (9) are higher in the downstream direction than the upstream direction for the IPDA method. This means that the information about the lowered downstream water surface travels faster upstream than in the IDWA method, which has the same a_i coefficients in the upstream and downstream directions. However, the IPDA method is more unstable. Note that smoothing and relaxation was required to get the IPDA to be stable, but this was not needed for the IDWA method in the current study.

Also note that since the IPDA and IDWA methods do not contain the time step, they will not give the correct solution for cases where the water level varies fairly rapidly, for example waves.

The computations were done over the weir, upstream of the hydraulic jump. The parameter D is dimensionless, scaled with the average water depth, H, at the weir. The number of measured water levels over the weir is denoted c. Table 1 gives the D values for the different algorithms. The VOF method gives the best result, and the IDWA method has the poorest accuracy. This corresponds well with the visual observations from Figures 6–9.

The VOF method used in, for example OpenFOAM, gives a very good estimate of the water surface over the weir. The main problem with the VOF and CGA method is the need for a relatively short time step. An adaptive grid method can use fewer cells with the same accuracy, but it has other drawbacks. The CGA method did not compute the water levels as accurately as the VOF method.

The implicit methods with adaptive grids (IPDA, IDWA) give reasonable results for the water level over the weir, although not as good as the VOF method. However, the implicit methods can use very long time steps, which is beneficial for simulating processes that take place over a long time. An important disadvantage with the IDWA method is that the resulting water level is a function of a friction coefficient, which requires calibration. All the other methods use a very low roughness for the skin friction and compute the other energy losses directly. Then calibration of a friction coefficient is not needed.

None of the free surface algorithms compute the water slope in the hydraulic jump correctly. This is probably not a serious problem for computation of reservoir flushing and other flows in most natural rivers. Large hydraulic jumps usually only occur downstream of hydraulic constructions. The coefficient of discharge for a weir is only based on the upstream water level, which all methods computed reasonably well.

The study also shows that a CFD program can be used to predict the water surface over a spillway on a relatively coarse grid and still produce high quality results. The computational times for all the tested algorithms were reasonable for engineering purposes on a desktop PC.

I want to thank Geir Tesaker for the construction work on the physical model of the weir and Heide Reppenhagen for assistance with the laboratory measurements.