FVCOM model simulation of local scouring around bridge pile

Scouring is one of many damages that water can cause. Scouring can occur as a consequence of bridge pile existence. The problem on local scour around single pier will be studied by using FVCOM numerical model. This study objective is to find out how accurate FVCOM model to predict local scour behavior. FVCOM model is based on the finite volume method to solve Navier Stokes, Meyer Peter Muller, and Exner equations. FVCOM computed numerical result then will be verified with computed and measured data in previous numerical (FSUM model) and experimental study. Results from this study show FVCOM model were successfully simulated typical features of local scour around piers such as downflow and wake vortex, but failed to simulate horseshoe vortex. Both computed numerical (FSUM and FVCOM) results are then compared with measured experimental data for its magnitude and time-series of maximum scour depth. FVCOM result shows value 0.99 r-squared correlation and 5.96 percent average error, and FSUM result shows value 0.98 r-squared correlation and 6.82 percent average error. Therefore, it can be deduced that FVCOM successfully predict local scour depth and its time-series and proven that FVCOM is more accurate than FSUM model.


Introduction
Water is the most important commodity for sustainability of human life on earth. Water has numerous benefits for human life, such as irrigation, hydropower, recreation, etc. Besides having countless benefits, water also has potential to damage people's lives. Water damage include flooding, erosion, drought, landslides, land subsidence occurrence, etc. Humans as water users should be able to anticipate the potential destructive force.
One potential water damage is scouring due to river flow on bridge pile. This potential damage should be anticipated since bridge is important for transportation and has implications on the economy and welfare of the community. To be able to anticipate bridge pile local scour, first we need to do an assessment of the magnitude of the scour potential.
FVCOM model is used to simulate local scour around single bridge pile. This model is benchmarked to previous numerical and experimental study conducted by Thanh et al. [1]. Thanh et al. simulated local scour of single bridge pile using FSUM model and benchmarked with experimental results. This study is intended to reveal capability and accuracy of FVCOM model compared to FSUM model to simulate local scour around single bridge pile.

Literature review 2.1 Fluid characteristic
Although the differences in fluid and solid objects can be explained qualitatively based on its molecular structure, a more specific differences is based on how the fluids and solids deform due to external forces. The fluid is defined as a substance that continually deforms (flows) when given shear force. Although solids also can deform (value of the deformation is very small), but a solid object deformation is not continuous (no flow) [2]. Some materials, such as slurries, tar, putty, toothpaste, and so on, are not easily classified into solid objects or fluids. This is because the material behaves like a solid when given a small shear force, but if the stress exceeds some critical value, the substance will flow. Studies of this kind is called rheology and are not included in the discussion of the fluid in this study. Munson et al. [2] explained that laws that apply in the fluid is Newton's laws of motion, conservation of mass, and the first and second laws of thermodynamics. Thus, there are similarities between the approaches used in the analysis of fluid motion with the approach used in solid objects.

Scouring and sediment transport
Annandale [3] defines scour as an excessive erosion event at a particular location. Erosion is transported sediment material event due to a specific flow phenomena, such as deflection of flow around the bridge pile, narrowing flow at bridge abutments, supercritical flow from dam spillway, and so forth. Basically, the information necessary to analyze the phenomenon of scours include: sediment ability to resist, sediment transport capacity of flowing water, and critical value associated with sediment transport capacity when sediment begins to move, known as incipient motion criteria.
Incipient motion occurs when value of sediment transport capacity due to water flow just exceeds value of sediment ability to resist erosion. This event is the beginning of erosion, and will continue until a new equilibrium is reached. A new equilibrium is a state where the sediment transport capacity due to the flow of water becomes equal to or lower than sediment ability to resist erosion. In such that condition, maximum scouring has been reached.
Quantification of value of sediment transport is not an easy task. In general, the approach is to use indicators or parameters of flow with ties or influence on sediment transport. Indicators or parameters used such as shear stress, average velocity, and stream power. The methods used to quantify the amount or value of the sediment transport has a very high variability and yield inconsistent trends. This causes practical problems for researchers and practitioners to analyze scour phenomenon.

Numerical methods
Numerical modeling of water flow phenomena is an analysis of water flow system and process based on the computer simulations. Versteeg and Malalasekara [4] describes the use of numerical methods several advantages compared to the experimental approach in the case of a fluid system, including: • Lower cost and modeling time.
• Capable to simulate the system where controlled experiments are difficult or impossible to implement (e.g. very large systems). • Ability to model the system with a high level of danger or outside the limits that can be retained by the system (for example, is the study of security and disaster scenario). • Unlimited levels od results details.
A numerical model is an algorithm calculations performed by the computer with the input or the input specified by the user. Numerical modeling is generally divided into three processes, namely pre-pocessor, solver and post-processor.
Pre-processor is input from water flow problems into numerical model consists of: definition of model domain, division of domain into small parts with simple geometry that do not overlap (also called a grid, mesh, cell, or element), the selection and discretization governing equation that describes the phenomenon of water flow, determining or defining property of water flow, and selection of flow boundary conditions that can represent actual water flow in reality.
There are three main techniques or methods to solve governing equations that describes phenomenon of water flow, finite difference method; finite volume method; and finite element method. Post-processor is the extraction process output, visualization of flow parameters, and interpretation of modeling results. Many computer applications can facilitate extraction process and visualization of modeling output. Visualization can be done by displaying animated and dynamic result on the computer screen.

Finite volume method
Mazumder [5] states that in finite volume method, governing equations in the form of partial derivative equations solved at each volume control. The first step is to divide computational domain into cells, as can be seen in Figure 1. Cells has simple shape and can be sized differently from one another. Boundary line between two cells known as cell surface. Points on cells called nodal. Flow parameters information are stored in each cell, called center of the cell. In Figure 1, cell center is a small box that is contained in the triangular cell, whereas nodal point is thick at both ends of the line that limits one cell to another. Unlike finite difference method, in finite volume method, the equation is not immediately resolved, but integrated into each cell in computational domain then approximation and completion are made. Because there are no cells outside the domain boundary, the boundary conditions can not be directly resolved.
The most important property in finite volume method is called conservation property. The amount of variable transfer rate per unit area of the cell surface called a flux. In case of diffusion, the amount of flux is proportional to the gradient of variable magnitude being simulated. Fundamentally, finite volume method is a flux balance equation. Thus, conservation characteristic is always automatically maintained in finite volume method.

FVCOM model
FVCOM is a 3-D model of primitive equations for free surface flow simulation with finite volume approach on irregular grid developed by Chen and Beardsley [6]. FVCOM model has been equipped with sediment calculation module refers to ROMS model consisting of suspended sediment and bed sediment tramsport. ROMS is a model of sediment transport with a structured grid which is then modified to accommodate irregular grid (unstructured) on FVCOM model.
In FVCOM finite volume method, triangular mesh is applied with an irregular grid. Each element consists of three nodes, one centroid, and three sides. Size of mesh formed in domain is defined by user. To improve accuracy of computational results, vectorial variables u and v (speed x and y) placed at elements centroid (center point), while scalar variables are placed at nodal points.

Local scour around bridge pile
Local scour is bed scouring that occur locally around the buildings or obstructions in a stream of water. Due to the buildings or other obstructions, the flow was interrupted, and the parameters of the flow changes.
Adikesuma [7] states that basic mechanisms that cause scours is a vortex that occurs at the bed, known as Horseshoe Vortex. Ameson et al. [8] gives a schematic representation of the dimensions of the phenomenon of Horseshoe Vortex at bridge piers as can be seen in Figure  2 and Figure 3.
Whirlpool also occur in downstream of building or barrier known as Wake Vortex as can be seen in the illustration in Figure 3. The intensity of Wake Vortex quickly diminished with increasing distance downstream of buildings or obstructions. This then leads to the disposition of the base material on the downstream structure.

Hydrodynamic Equations
By ignoring Coriolis parameter and non-hydrostatic pressure, the governing equations used in FVCOM model is continuity and momentum equation as follows: (4) is continuity equation where x, y, and z are the Cartesian coordinate system, while u, v, and w are the components of velocity in directions of x, y, and z, ρ is density of water, P is hydrostatic pressure, g is gravitational acceleration, Km is vertical component of Eddy viscosity, and Fu, Fv represents horizontal force. Illustration of an orthogonal coordinate by Chen et al. [9] can be seen in Figure 4, where total water column is calculated using equation (5), H is bed elevation relative to z = 0, and ζ is elevation of free surface water relative to z = 0.

Equations (1) -(3) is momentum equations, and equation
= + (5) where p = pa + pH + q and pH meet the following conditions: Values of u, v, and w parameters on surface and bottom boundary condition at the absence of evaporation, transpiration and bottom fluxes can be seen as in equation (7) and (8). = + + at z = ζ (x, y, t) = − − at z = -h (x, y) To solve the problems of irregular bathymetri, use the coordinates of sigma (σ) as defined as:  where C is a constant parameter, and Ω u is the area of the momentum control element. Am value varies depending on model resolution and the horizontal velocity gradient. Am value decreases as the grid size or horizontal velocity gradient decreases.

Scouring and sediment transport equation
FVCOM model have been equipped with module to calculate suspended sediment and bed sediment transport, and there are areas near the bed where there is no clear distinction between these two sources. On implementation, suspended sediment and suspended sediment are calculated separately and added together to produce total sediment transport.

Suspended sediment transport equation
Suspended sediment model use concentration approach with the following evolution equation: where Ci is sediment concentration i, Ah is horizontal eddy viscosity, Kh is vertical eddy viscosity, and wi is the settling velocity. On the surface there is no-flux boundary condition, while at the bottom, sediment flux is the difference between the deposition and erosion. Both of these boundary conditions are mathematically presented in equation (16) where Qi is the erosive fluxes, Pb is the porosity of sediments, and τci are critical shear stress of sediment i. Sediment eroded when the local shear stress that occurs is greater than critical shear stress. Balanced advection, vertical diffusion, sediments transported and deposited sediment will generate a sediment concentration profile. Settling term (4th term in the evolution equation sediment concentration) need to be considered with caution, because there is a sharp profile of concentration gradient near the base. FVCOM model used flux-limited scheme. This scheme use min-mod limiter and eliminated second-order accuracy in extreme gradients. Settling flux on the bed is stored as a deposition flux in the equation (17).

Bed sediment transport and morphology equations
To calculate thickness of eroded sediment, Exner mass balance equation is used, such as the following: Where η is a bottom elevation, t is time, ε0 is porosity, and ▽ qs is bed sediment flux gradient in space. Sediment flux is computed by using Meyer-Peter and Muller schemes. In this scheme, the sediment flux is calculated based on the difference of actual stress with a critical stress, such as the following: * = ( * − * ) 3 where qs is volumetric sediment discharge per unit width, u* is shear velocity, s is specific weight of sediment, g is gravitational acceleration, D is water depth.

Finite volume discretization
Similar to the finite element method, numerical computation domain is divided into several horizontal triangle cells. Cells triangle consists of three nodes, one centroid, and three sides ( Figure 5).
Integration of equation (13) to the area of triangular cells provide: where vn is velocity component normal to the cell edge and s is three sides of a triangle. Equation (23) are integrated numerically by the 4 th order Runge-Kutta scheme, with following procedures: where k = 1, 2, 3, 4 and (α1, α2, α3, α4) = (1/4, 1/3, 1/2, 1). n represents n-th time-step. Ω is area bounded by the midpoint of the side of a triangle around the nodal point.
With the same way integration, momentum equations and sediment transport can be discretized using finite volume scheme. Space discretization is based on flux balance on each cell, and time discretization using 4 th order Runge-Kutta method.

Results and discussion
Domain and physical parameters used in this FVCOM model simulation refer to experimental study conducted by Thanh et al. [1]. Domain is 6 meters long and 1.1 meters wide. 10 cm x 10 cm rectangle pile is placed right in the middle. Water flows at constant discharge 52 liters per second with 22 cm water depth, and initial bed sediment thickness is 15 cm.
Simulation results show that bed roughness value affected rate and maximum scour that occurred around bridge pile. The greater value of bed roughness, the greater and faster the scouring that occurs. Figure 6-7 show current velocity vectors and magnitudes on vertical cross section upstream of bridge pile using FVCOM model. Figure 8 shows vertical vortex that occurs in the upstream pillar based on a study conducted by Thanh et al. [1]. Figure 6-7 show current velocity vectors and magnitudes on horizontal cross section around bridge pile using FVCOM model.       [1] in Figure  8, this is caused by the velocity vector of the vertical direction is a result of the movement of the axis of sigma and not the actual vertical direction movement. Even though, FVCOM model simulation is able to capture horizontal vortex phenomena (wake vortex) on downstream of the pile very well, as can be seen in Figure  9-10.
Scouring that occurs in the FVCOM model simulation along with comparisons with experimental and FSUM model simulation conducted by Thanh et al. [1] presented in a time-series graph which can be seen in Figure 11. Pearson correlation value and the error rate of FVCOM model results compared with the results of experimental studies by Thanh et al. [1]. Correlation value and the error rate is then compared with a correlation value and the error rate resulted by FSUM numerical modeling by Thanh et al. [1]. These values can be seen in Table 1.