A higher-efficient non-hydrostatic finite volume model for strong three-dimensional free surface flows and sediment transport

In order to accurately simulate strong three-dimensional (3-D) free surface flows and sediment transport, the fully 3-D non-hydrostatic pressure models are developed based on the incompressible Navier-Stokes equations and convection–diffusion equation of sediment concentration with the mixing triangle and quadrilateral grids. The governing equations are discretized with the unstructured finite volume method in order to provide conservation properties of mass and momentum, and flexibility with practical application. It is shown that it is first-order accurate on nonuniform plane two-dimensional (2-D) grids and second-order accurate on uniform plane grids. A third-order approximation of the vertical velocity at the top-layer is applied. In such a way, free surface zero stress boundary condition is satisfied maturely, and very few vertical layers are needed to give an accurate solution even for complex discontinuous flow and short wave simulation. The model is applied to four examples to simulate strong 3-D free surface flows and sediment transport where non-hydrostatic pressures have a considerable effect on the velocity field. The newly developed model is verified against analytical solutions with an excellent agreement.


Introduction
An accurate prediction of water flow and transport is important to hydraulic and coastal engineering. In order to simulate complex free surface flow, hydrodynamic models based upon the three-dimensional (3-D) Navier-Stokes equations seem a promising approach with the development of computation technology. In recent years, more and more experts and scholars have devoted themselves to the research of the mathematical model of three-dimensional hydrodynamics.
Based on the full Navier-Stokes equations of incompressible flow, the SLMPLE series algorithm (Stansby and Zhou, 1998;Jia et al., 2010;Hu et al., 2011) based on the computational heat transfer is widely applied to solve the Navier-Stokes equations of incompressible flow. In the computation, the computation time of the Possion pressure equation often accounts for more than 70% of the total computing time. It is because the incompressible flow pressure wave transmission speed is infinite, the instantaneous pressure perturbations spread to the whole flow field of the elliptic equation properties. Therefore, the pressure correc-tion iteration convergence is relatively slow for meeting the incompressible condition of pressure.
A series of the 3-D non-hydrostatic pressure mathematical models are structured based on the idea of the explicit projection method (Chorin, 1968;Li and Fleming, 2001;Lin and Li, 2002), implicit pressure method (Namin et al., 2001;Yuan and Wu, 2004) and semi implicit fractional step (Casulli, 1999;Casulli and Zanolli, 2002;Chen, 2003). The 3-D mathematical model is mostly developed by using structured or unstructured grid in the horizontal direction and vertical layer in the vertical direction. The total computation time increases exponentially with vertical layers. Therefore, how to reduce the number of vertical layers by improving the calculation efficiency is the research focusing on their model under the premise of reducing the accuracy of the model.
For simulating three-dimensional non-hydrostatic flows, the computational cost is a critical issue. Many efforts have been made to save the computational cost by reducing the number of vertical layers. For example, Stelling (2005, 2008) proposed the Keller-box method to replace the staggered grid in the vertical direction, which enables the pressure to be located at the cell faces rather than the cell centers. The zero pressure boundary condition at the free surface can be exactly assigned to zero without any approximation. The coefficient matrix of the Possion pressure equation is not symmetrically positive definite or diagonally dominant, which will cost a considerable computational time. Yuan and Wu (2004) proposed an integral method, different from the Keller-box scheme, to obtain a non-hydrostatic pressure condition at the free surface cell under a staggered grid framework. The order of the Possion pressure equation coefficient matrix is four times as much as the other way, which also needs a lot of computation time. Badiei et al. (2008) and Ai et al. (2011) presented non-hydrostatic finite volume models based on a vertical boundary fitted grid system. Their satisfactory results demonstrated the capability of simulating strong 3-D free surface flows with a small number of vertical layers. However, their models are based on structured grids in the horizontal projection of 3-D physical domain. Therefore, continuous efforts to develop an efficient non-hydrostatic model are highly desired.
In curved open channels, the secondary currents in the transverse direction cannot be neglected (Zeng et al., 2008), as the secondary flows are combined with the primary flow to produce the well known spiral motion, which are turbulent and of strongly three-dimensional behavior. They play an important role in the sediment transport, bed shear stress distribution and even the development of the riverbed profiles. Prediction and knowledge of the three-dimensional flow patterns as well as the mass transport in the channel bends are important and challenging in river engineering.
Over the years many papers about 3-D numerical models of flow in curved open channels have been published (Wu et al., 2000;Lesser et al., 2004;Rüther and Olsen, 2005;Khosronejad et al., 2007;Zeng et al., 2008). The three-dimensional models represent much more comprehensive flow field information, including detailed vertical structure of the flow field. Furthermore, it also accounts for historical or non-local effects and can accurately capture the intensity of the secondary currents directly in curved bends, which prevents it from being corrected by extra semi-empirical sub-models.
In the present study, a finite-volume formula is employed for the 3-D Navier-Stokes equations and convection-diffusion equation of the sediment concentration on the mixing triangle and quadrilateral grids. An Eulerian scheme for advection diffusion terms is employed in this model, with the attractive property of being conservative. By incorporating the integral method of the top-layer pressure treatment, the developed model first calculates the intermediate velocities for the advection and diffusion terms. The intermediate velocities are updated by a non-hydrostatic step of which its Possion pressure equation is obtained by employing a divergende-free velocity field. A bi-conjugated gradi-ent method with a preconditioning procedure is used to solve the pressure correction equation. Finally, the model is applied to four examples including linear 3-D standing short wave in deep water, long wave resonance in a parabolic basin, the U-shape bend flume test and the S-shape bend flume test to demonstrate the capability, efficiency and accuracy of the developed model using very few vertical layers.

Governing equations
A non-hydrostatic free surface flow is governed by the continuity equation and three-dimensional incompressible Navier-Stokes equations based on the conservation of mass and momentum, which can be expressed in the following forms by splitting the pressure into hydrostatic and non-hydrostatic ones, i.e. p=g(η-z)+q.
where u, v, and w are the velocity components in the x-, yand z-directions, respectively; ν is the turbulent eddy viscosity; g is the gravitational acceleration; η is the free surface elevation; q is the non-hydrostatic pressure component. The continuity Eq. (1) is integrated from the bottom z=-h(x, y) to the surface z=η(x, y, t), which yields the free surface equation To calculate the suspended load, the convection-diffusion equation of the sediment concentration is solved to find the sediment concentration distribution throughout the suspended-load layer.
where S is the sediment concentration; w s is the fall velocity of the sediment particles; K h and K v are the diffusion coefficient in the horizontal and vertical direction, respect-ively, K h is constant, and K v =v/σ c , where σ c is the turbulent Schmidt number. With knowing the velocity profiles and the concentration of suspended sediment load, the suspended sediment transport rates can be integrated over the suspended layer The sediment continuity equation, integrated over the water depth, is used to calculate the bed-elevation changes where p′ is the porosity of the bend; z b is the bed elevation; q s is the suspended load rate.

Boundary condition
At the free surface, for the case of a wind stress on the surface, the shear stress at the free surface is given by τ w xz τ w yz where and are the components of the wind stress in the x-and y-directions, respectively.
At the free surface, the vertical sediment flux is zero, and it can be expressed as: At the bottom, the no-slip condition is applied together with a zero normal velocity component to the bottom. At the bottom boundary of the suspended layer, the net sediment flux equals the difference between deposition rate D b and entrainment rate E b at the interface, and hence the condition applied is where S a is the equilibrium sediment concentration, and it can be determined by employing the empirical formula of van Rijn (1987) where d 50 is the mean diameter of sediment particle; d * is the dimensionless particle diameter; T * is the dimensionless shear stress; and δ is the thickness of the bed-load layer.
At the solid impermeable boundaries, the slip boundary conditions are used by setting all normal components to the vertical solid equal to zero.
At the inflow, the velocity components are specified by use of the standard formulae for open channel flow, and the sediment concentration is set to zero.
At the outflow, η is defined and zero normal gradients are specified for all the other variables.

Spatical discretization and definition of variables
Owing to the flexibility of the triangle mesh and minimization of the quadrilateral mesh, the complex boundary can be accommodated and the amount of meshes can be efficaciously reduced compared with the triangle mesh. In this model, the horizontal projection of the 3-D physical domain bounded by the moving free surface, z=η, and the bottom, z=-h is discretized by means of the orthogonal unstructured triangular and quadrilateral grids. The 3-D space discretization consists of 3-D prisms which usually has five or six faces, depending on whether the 2-D grids are the triangles or quadrilaterals. A finite volume method is used to discretize the governing equations and boundary conditions on the horizontally orthogonal unstructured triangle and quadrilateral grids. The vertical spacing is plan-wise constant at any z-level. The velocities and free-surface elevations are defined at staggered locations in the horizontal domain. In the vertical direction, the velocities are assigned at the layer faces, and the pressures and scalar transport variables are assigned at the center of the layer. Details about the grid arrangement and definition of variables can be found in Lv et al. (2010). The space discretization and definition of variables are shown in Fig. 1.

Numerical algorithms
A semi-implicit fractional step scheme is adopted to solve the Navier-Stokes equations in two major steps. In the first step, the intermediate velocities are achieved by solving the momentum equations which contain the advection terms, diffusion terms and the old step pressure terms. The advection and diffusion terms in the momentum equations are discretized explicitly with an Euler scheme. In the second step, the new step pressure terms are computed by the pressure correction equation, which is obtained by the combination of the discretized continuity equation and momentum equation.

Discretization of the horizontal momentum equations
Since the horizontal momentum Eqs. (2) and (3) are invariant under solid rotation of the x-and y-axis around the zaxis, by neglecting the pressure term at the new time step, the horizontal velocity component on each vertical face of a prism (see Fig. 1) can be discretized as follows: where denotes the horizontal normal velocity component to the j-th side of the grid at vertical level k and time step n. θ=0.5 is an implicit factor. F is an explicit operator which accounts for the contributions from the discretization of the advection and diffusion terms , and are the horizontal and vertical advection and horizontal diffusion terms.
By using Perot's scheme (Perot, 2000), the form of the horizontal advection terms can be given by In Eq. (15), and are discretized advection terms, which are defined at the centers of the two neighboring cell of the j-th side. For defined on cell i, it can be obtained by integrating the horizontal advection term over the finite volume V i where A i denotes the area of each cell i. N j is a sign function associated with the orientation of the normal velocity defined on the j-th side of cell i. Specifically, N j =1 if a positive velocity on the j-th side corresponds to the outflow, and N j =-1 if a positive velocity on the j-th side corresponds to the inflow to the i-th water column. .
C V (U j,k ) The form of vertical advection terms can be obtained by The component of the vertical advection at the Voronoi point in cell i in the direction of the normal vector is n j given by Approximating this velocity vector with a linear interpolation at the upper and lower faces yields Horizontal diffusion term is discretized with the same conservative discretization as for the advection, such that where the face-normal gradient is given by ∂u ) and the eddy viscosity is approximated with the linear inter-

Discretization of the vertical momentum equations
In analogy with Eq. (13), the vertical momentum Eq. (4) can be discretized by the finite volume method as follows: where H(w) is an explicit operator where the explicit terms , and account for the contributions from the discretization of the horizontal advection, vertical advection and horizontal diffusion terms, respectively. After a finite volume discretization, they have a similar form as Eqs. (15), (18) and (22).

Discretization of the free surface equation
The new intermediate velocity field must satisfy each grid cell i of the finite volume representation of the face surface Eq. (5). Thus, the free-surface Eq. (5) is discretized by a semi-implicit finite volume method and expressed in the matrix notation giving where B denotes coefficient matrix of the horizontal and vertical advection and horizontal diffusion terms, and G denotes vectors of at the j-th side of the cell i. This linear system is symmetric and positive-definite, which can be solved efficiently with the preconditioned conjugate gradient method. With the free surface determined from Eq. (28), the predicted velocity field can be obtained from Eqs. (13) and (26).

Non-hydrostatic pressure correction
In the second step, continuity Eq. (1) together with the momentum equations excluding advection and diffusion terms are solved as follows: A finite volume method is used to discretize the continuity Eq. (1) as: By substituting the expressions for the new step velocities from Eqs. (31)-(33) into continuity Eq. (34), the non-hydrostatic pressure correction equation excludes half top layer.
In some non-hydrostatic models, non-hydrostatic pressure is defined at the two horizontal faces of a prism and the zero pressure boundary conditions at the surface cannot be given in these models. So, a large number of vertical layers are necessary to simulate strong 3-D free surface flows.

0.25∆z i,nzt
To minimize the number of vertical layers and subsequently the computational cost, an integral method of the vertical momentum equation at the top-layer is applied. Lv (2014) employed a third-order approximation of the vertical velocity at the top layer located at a distance of from the surface. This model obtained satisfactory results by employing a small number of vertical layers. In this model, a third-order approximation is adopted. Details of the discretized non-hydrostatic pressure correction equations have been described in Lv (2014).
Non-hydrostatic pressure correction equation for the top layer together with the non-hydrostatic pressure correction equations of the lower layers which are obtained by substituting Eqs. (31)-(33) into the continuity Eq. (6) to form a tri-diagonal block matrix system which can be written as where A is a sparse coefficient matrix, P is a vector of the non-hydrostatic pressure correction, and B is a known vector related to the intermediate velocities. Specifically, if the 2D grid is built from triangles, the sparse matrix A contains eight non-zero diagonals at the bottom layer and twelve non-zero diagonals at other layers. If the 2D grid is built from quadrilaterals, it contains ten non-zero diagonals at the bottom layer and fifteen non-zero diagonals at other layers. Thus, Eq. (36) can be solved by a direct matrix solver. The details of this solving method have been described in Ahmadi et al. (2007).
Once the pressure terms are obtained, the new velocity field is updated by substitution in Eqs. (31)-(33).

Wet-dry boundary tracking
Similar to the work of Yamazaki et al. (2009) andAi (2012), the following algorithm is applied for the treatment of wet-dry fronts. According to the wet or dry category of cells and sides, the governing Eqs. (1)-(4) and (5) will be solved as follows. The continuity Eq. (1) and the momentum Eqs. (2)-(4) are solved by a semi-implicit fractional method mentioned above in wet cells and sides, respectively. In this procedure, the non-hydrostatic pressure terms are set to be zero at the wet-dry boundary and dry cells, so the Poisson equation is only solved for the wet ones. Finally, the free surface Eq. (5) is solved only for the wet cells and the solutions for the wet-dry boundary cells are extrapolated from the neighboring wet ones.

Discretization of the sediment transport equation
The partial differential convection-diffusion equations for the sediment concentration are integrated over a vertical layer to get the semi-dicretized equations using the finite volume method. The advection terms and diffusion terms are discretized as the flow governing equations and the firstorder upwind scheme is utilized for the flux calculation.
where the explicit terms , and account for the contributions from the discretization of the horizontal advection, vertical advection and horizontal diffusion terms, respectively. After a finite volume discretization, they have the similar forms as Eqs. (15), (18) and (22).
The time derivative term in the sediment continuity equation was discretized using the backward finite difference scheme. And the flow and bed deformation computations are coupled, meaning that once the flow field and sediment concentration distribution are calculated at time t, the bed changes due to the imbalance of the sediment transport rate are determined for a time interval ∆t. The vertical grid thickness was adjusted according to the water depth variation due to the bed changes. The time step of the sediment computation was selected to guarantee that the bed changes are smaller than 10% of the flow depth.

Linear 3-D standing short wave in deep water
To represent the effects of non-hydrostatic pressure distribution, the linear 3-D standing short wave in deep water is given. The linear 3-D standing short wave oscillates in a 10 m×10 m×10 m closed basin. The horizontal domain is discretized by 40×40 quadrilateral grids. The time step is taken as 0.05 s. The wave amplitude is A=0.1 m and the period is T=3.01 s. Details of analytical solutions of the linear standing wave can be found in Lv et al. (2010).
Figs. 2 and 3 present comparisons of the free-surface elevation between analytical solutions and numerical results with hydrostatic and non-hydrostatic pressure model by applying these two approximations at the top layer. Fig. 4 presents comparisons of the velocity (u, v, w) between analytical solutions and numerical results with non-hydrostatic model by applying the third-order approximation.  LIU Xin et al. China Ocean Eng., 2017, Vol. 31, No. 6, P. 736-746 741 This case demonstrates that the hydrostatic model with thirty vertical layers fails to circulate the standing short wave period and phase. For strong 3-D free surface flows, the hydrostatic pressure approximation is invalid. Simple averaging approximation cause noticeable phase errors with four vertical layers. However, the third-order approximation yields a more accurate simulation of the wave motion, using only four vertical layers.
The average CPU time needed for different number of layers is shown in Table 1. The results confirm the sufficiency of using four vertical layers and 40×40 horizontal grids to simulate strong 3-D free surface flows.

Long wave resonance in a parabolic basin
To validate the ability of the model to capture wet-dry boundary, the long wave resonance in a parabolic basin is given.
The long wave resonance being considered takes place in a basin described by a paraboloid of revolution given by where, h 0 is the still water depth in the center of the basin, r is the distance from the center point, R is the still water radius. Details of analytical solutions of the long wave resonance in a parabolic basin can be found in Ai and Jin (2012). The model is tested using the physical values with h 0 =1 m, R=2500 m, and r 0 =2000 m (Fuhrman and Madsen, 2008). In the numerical experiment, the horizontal grid spacing of ∆x=∆y=80 m and a time step of ∆t=T/240=7.388 s are employed. Only four layers are adopted here.
Figs. 5 and 6 show the time series of the computed and the analytical values of the free surface elevation and velocity, respectively. Fig. 7 shows the comparison of free surface elevation between the numerical and analytical solutions along the cross section y=0, which begins at the second oscillation and covers a half-period. Good agreements are obtained by the present model, demonstrating the accuracy of the wet-dry boundary algorithm. Fig. 8 shows the sketch of a U-shape curved laboratory flume carried out at Changsha University of Science and Technology. The flume consisted of three sections: a 180°b end with a constant centerline radius of curvature 2.00 m, being connected with a 4.00 m long straight channel reaching the upstream and followed by a straight outflow section with the length of 4.00 m at downstream. The cross-section of the flume was rectangular with its width of 0.80 m and its sidewalls (bank) were vertical. The average flow depth in the test case was 0.25 m and the mean flow velocity was 0.25 m/s. The numerical simulations started with a flat sand bed, similar to the experiment, and continued until the equilibrium flow conditions were reached. The simulation replicated the flow condition with 12800 grid cells in which, 80,  (u, v, w) between analytical solutions and numerical results with non-hydrostatic model.  The comparison of the vertical profiles of the streamwise velocity between the measured and simulated results in typical cross sections along the channel bend is demonstrated in Figs. 9 and 10, respectively. The horizontal and vertical axes correspondingly represent the velocity magnitude and water depth. Overall, the predicted velocities show a good agreement with the measurement data at two sections. The simulations reproduce the main features of the velocity distribution over the depth along the channel bend. Fig. 11 shows the secondary flow field being simulated at two typical cross sections along the bend. It can be found the transverse velocity flow towards the outer bank is near the water surface, while the inner bank is near the bed.

S-shape bend flume test
To investigate the applicability of the numerical model in the present study, results simulated with it were compared with the results measured with respect to an S-shape bend flume experiment performed by Onishi et al. (1972). The platform of the flume is shown in Fig. 12. The flume initially had a flat sand bed in the lateral direction. The duration of the flow in the flume was long enough to establish the equilibrium bed topography. The inflow water depth was 0.133 m and the mean inflow velocity was 0.536 m/s. The median diameter, d 50 of bed material was 0.25 mm.
The numerical simulations started similarly to the experiment until the equilibrium bed was reached, which was defined as the topography that does not change in time. In the numerical simulation, it employed 17100 grid cells in which 114, 30, and 5 in the longitudinal, transverse and vertical directions, respectively. The straight reaches up and downstream are expanded to 10.0 m to weaken the influence by boundary conditions. Fig. 13 compares the measured and predicted suspended sediment concentration profiles at two typical cross sections A and B. The overall predictions of the present model are in reasonable agreement with the measurements except for somewhat sharper decay near the bed bottom. Fig. 14 shows the comparisons of the water depth at two Sections A and B around the flume under the equilibrium condition. The bathymetry in bends is characterized by a transverse bed slope, in the first bend, and scour in the outer half of the cross-section (right bank) was more serious than that in the inner half of the cross-section (left bank), while in the last bend, the status was opposite. The simulation catches the characteristics of the bed scouring that are unique in the channel bend. The magnitude of the bed  LIU Xin et al. China Ocean Eng., 2017, Vol. 31, No. 6, P. 736-746 changes and the location of the maximum values matched the measurements from the physical model acceptable.

Conclusions
In this study, a higher-efficient non-hydrostatic finite volume model by solving the Navier-Stokes equations and convection-diffusion equation of the sediment concentration is developed to simulate strong 3-D free surface flows and sediment transport based on the staggered unstructured triangle and quadrilateral grid. An Eulerian scheme for ad- Fig. 9. Comparisons of streamwise velocity at the sections of 60° and 150° between the measured and simulated results.

744
LIU Xin et al. China Ocean Eng., 2017, Vol. 31, No. 6, P. 736-746 vection and diffusion terms is employed in this model and, while it has the attractive property of being conservative. The intermediate velocities are updated by a non-hydrostatic step of which its Possion pressure equation is obtained by employing a divergende-free velocity field. To validate the model, four examples including the lin-ear 3-D standing short wave in deep water, the long wave resonance in a parabolic basin, the U-shape bend flume test and the S-shape bend flume test were performed with a small number of vertical layers. Comparisons between numerical results and analytical or experimental data were carried out. The results for the first test confirm that the thirdorder approximation treatment of the top layer pressure yields a more accurate simulation of strong 3-D free surface flows motion by comparing hydrostatic pressure approximation with simple averaging approximation. The simulation of the long wave resonance in a parabolic basin demonstrates that the non-hydrostatic model capable of accurately capturing the wet-dry boundary. The final simulations for the U-shape bend flume test and S-shape bend flume test    show reasonable results, which manifests that the proposed model can be applied to simulate the velocity distribution over the depth along the channel bend and the subsequent bed deformation calculation. Good agreements between numerical results and experimental data indicate the capability, conservation and efficiency of the model for the simulation of strong three-dimensional free surface flows and sediment transport.
In conclusion, the non-hydrostatic model is accurate and efficient to solve a wide range problems of strong 3-D free surface flows and sediment transport.