Evaluation of Three-Dimensional Environmental Hydraulic Modeling in Scour Hole

Abstract

The main goal of this study was comparing the performance of an open-source code OpenFOAM and a commercial software Ansys Fluent in simulating the turbulent flow through a scour hole developed in a sand bed channel, which helps to give a hint in choosing the appropriate calculating tool. Both models were set with the same mesh and as similar as possible numerical settings, with RANS turbulence modeling, applying the k-ωSST model, in transient simulations. The results of flow pattern, velocity, and turbulence properties were collected and compared with laboratory experimental data. The analyzed results showed that, although both of the two models cannot perfectly reproduce the values from a laboratory experiment, they can quite well capture the flow in scour hole near the wall, with a bit higher performance coming from the OpenFOAM model application.

1. Introduction

Due to the rapid development of computing technology, currently using computers have the super computational capacity that leads to a boost in modeling many different issues or simulating complicated systems. This method is widely accepted as a less expensive option compared to the laboratory experiment or real-field operation when related to large-scale studies. Instead of building a big and high-cost apparatus for testing new concepts, numerical modeling can reduce the work down to some effective computer codes. The reliable software can help to investigate or provide more insight on the behavior of large systems that are limited in measurement options [1].

This method fits well while studying large environmental systems such as rivers, estuaries, or oceans. A scour hole developed on a sandy bed channel near hydraulic structures, such as a bed protection of weir, has a considerable effect not only on the safety of the structure, but also on the habitat. In general, physical studies of the scour hole in a hydraulic laboratory take more time and effort than numerical studies, while the former is much more accurate if the mathematical equations for the physical mechanism are not well developed. The flow characteristics are complicated in the scour hole with recirculation, sediment transport, and effect of boundary layer.

Both physical and numerical modeling studies were conducted by many previous researchers in scour holes. In the most of numerical modeling studies, the results and accuracy depend on selection of numerical model and method due to the complexity of physical presses. There are many options for numerical modeling. A free computational fluid dynamic (CFD) code named OpenFOAM and the commercial software Ansys Fluent were chosen to compare in this present study. OpenFOAM was developed by OpenCFD Ltd. and OpenFOAM Foundation [2], and the company provides it as an open-source code; therefore, the users can modify or add more codes for the part they want. OpenFOAM also has a wide range of applications of fluid dynamics, with no limitation for parallel computing. These characteristics make it become one of the most widely used CFD packages recently in both industrial and academic sectors [3]. Meanwhile, ANSYS Fluent is a commercial software model developed by ANSYS Inc [4]. Besides the strong capacities for CFD modeling, this option also comes with a powerful graphical user interface and various supporting services, such as the integration with several grid-generating and post-processing software, and a typical support team [5]. However, the users have to pay for the license and cannot interfere with its inside codes and functions. This study focused on a comparison of the accuracy of these widely used tools and especially their performances in hydraulic modeling of flow in the scour hole near the bottom wall.

Previously, several related research studies have been done to compare the performance of these two CFD codes in simulations of specific targets. Lysenko et al. [3] conducted an experiment using OpenFOAM and Ansys Fluent for simulating the turbulence separated flows, and their analyzed data show essentially equal results. A possibility to archive the agreeing results from both models was proved by Ambrosino and Funel [6] when they examined the exterior flow field around simplified passenger sedan geometry. However, research result from Balogh, et al. [7] has shown some differences. In the study by Balogh, et al. [7], the authors used Reynolds Average Navier Stokes (RANS) with k-ε turbulence model for both OpenFOAM and Ansys Fluent codes for the simulations. They concluded that the flow velocity was more accurately predicted by OpenFOAM, but the turbulent kinetic energy was more accurately predicted by Ansys Fluent. For the flow near the wall, such as flow in scour hole, there also are many studies conducted by OpenFOAM [8,9,10,11,12] or Ansys Fluent [13,14,15,16,17]. Still, there are not many works focusing on the two models’ comparison of the flow near the wall according to the authors’ knowledge. Therefore, in the present study, simulation of the turbulence flow through a scour hole in a sand bed channel was conducted by both OpenFOAM-v1712 (released in 2017) and Ansys Fluent 19.1 (in 2018) for comparing their accuracy and performance in simulation. Both models are utilized for 2D modeling of laboratory-scale experiments taken from our own research.

The purpose of this work is focusing on the performing of simulations rather than analyzing the inner workings of those two models. This study tries to find out which software shows better simulation for the flow near the bottom wall. The results from this study will contribute as a reference for choosing the appropriate modeling tool and scheme for further research, similar to numerous other previous works [3,6,7]. The computational mesh was created by OpenFOAM and then converted to the Ansys Fluent mesh input form. Both simulation models and boundary conditions were set as similar as possible between the two models. The results of water flow behaviors such as streamlines, velocity profiles, and turbulence kinetic energy, were analyzed and compared to our own laboratory experiment data.

2. Methodology

2.1. Experimental Study

The results from laboratory experiments were used for comparing the two numerical models. In this study, a laboratory experiment with the setting as presented in Figure 1 was conducted to investigate the hydraulic properties of a scour hole developed in a sand bed channel. This physical study was performed by using a 12.5 m long and 0.6 m wide channel made of two parallel vertical glasses. The input water runs into the channel controlled by a hydraulic pump, and a volumetric flow rate (Q) of 0.35 m³/s was maintained. A tail gate was set at the end of the channel to control the water depth of 0.144 m. The detail conditions of the experiment are as shown in

Table 1. Physical experiment setup.

Case	Q (m3/s)	h0 (m)	u (m/s)	d50 (m)	Re
exp.	0.035	0.144	0.405	0.0012	58,300

, where h₀ denotes the water height (m) and d₅₀ presents the mean particle size by weight (mm). The inlet flow velocity u (m/s) can be calculated as the division of the volumetric flow rate Q (m³/s) and the area of the inlet boundary, A (m²):

$$u=\frac{Q}{A}$$

(1)

The reference length D (m) of the water height at the inlet (h₀) was used to calculate the Reynolds number as:

$$Re=\frac{uD}{\nu }$$

(2)

and $\nu =1\times {10}^{-6}$ is the kinematic viscosity of water (m²/s).

For reproducing the riverbed protection part, the experimental channel bed was set up with two different parts as shown in Figure 1. A 4.5 m long acrylic bed was used to produce a smooth condition at the upstream, while the downstream was of 8 m length filled with sand. In this work, the mean particle size by weight (d₅₀) was to 1.2 mm. This laboratory experiment took up to 12 days to perform until it reached the equilibrium state, defined as there being no further increase in the scour hole. At the final stage, an Acoustic Doppler Velocimeter (ADV) system was used to measure the water flow properties and bed profile.

2.2. Numerical Study

2.2.1. Governing Equation

Both models of Ansys Fluent and OpenFOAM use the RANS formulation in their codes, and are solved in the ensemble-averaged form. The governing equations of continuity and its RANS formula are as follows:

$$\frac{\partial u_i}{\partial x_i}=0$$

(3)

$$\rho \left(\frac{\partial U_i}{\partial t}+\frac{\partial }{\partial x_j}\left(U_i{U}_j\right)\right)=-\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}\left(2\mu S_{ji}-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$$

(4)

In these equations, u, U, and $u^{\prime }$ are the instantaneous, mean, and fluctuating velocities of the flow (m/s), respectively; t denotes the time in second (s); $\rho$ is the density (kg/m³), P is the pressure (Pa), $\mu$ is the molecular viscosity (kg/m·s), and $S_{ji}$ is the mean strain-rate tensor.

2.2.2. Turbulent Model

Based on the previous research [8,9], the $\mathit{k}-\mathit{\omega }\mathit{S}\mathit{S}\mathit{T}$ turbulence scheme, which was proved to be the best to investigate the turbulence properties of the flow in the scour hole, was employed in the present study. Both OpenFOAM and Ansys Fluent models use the $\mathit{k}-\mathit{\omega }\mathit{S}\mathit{S}\mathit{T}$ formula that was developed from Menter [19]. The idea of this approach is blending or switching the $\mathit{k}-\mathit{\epsilon }$ model with the $\mathit{k}-\mathit{\omega }$ model. This task can be done by the addition of a blending function, which can switch from one in the near-wall region to activate the $\mathit{k}-\mathit{\omega }$ model to zero while in the far field to change to the $\mathit{k}-\mathit{\epsilon }$ model. The standard $\mathit{k}-\mathit{\omega }$ model is a low Reynolds number model and can be used where the boundary layer is thick and the viscous sublayer can be resolved. On the other hand, the standard $\mathit{k}-\mathit{\omega }$ is best used for near-wall treatment. Meanwhile, the $\mathit{k}-\mathit{\epsilon }$ model is ideal for predicting the flow away from the wall. Therefore, the combination is less sensitive to the free stream conditions compared to the standard $\mathit{k}-\mathit{\omega }$ model and also avoids build-up of excessive turbulent kinetic energy near stagnation points as compared to the standard $\mathit{k}-\mathit{\epsilon }$ model.

In general, the transport equations for turbulent kinetic energy, k, and specific dissipation rate $\omega$ can be written in the simple forms as:

$$\frac{\partial }{\partial t}\left(\rho \mathit{k}\right)+\frac{\partial }{\partial x_i}\left(\rho u_i\mathit{k}\right)=\frac{\partial }{\partial x_j}\left[\left(\frac{{\mu }_t}{{\sigma }_{\mathit{k}}}+\mu \right)\frac{\partial \mathit{k}}{\partial x_j}\right]+G_{\mathit{k}}-Y_{\mathit{k}}+S_{\mathit{k}}$$

(5)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_i}\left(\rho u_i\omega \right)=\frac{\partial }{\partial x_j}\left[\left(\frac{{\mu }_t}{{\sigma }_{\omega }}+\mu \right)\frac{\partial \omega }{\partial x_j}\right]+G_{\omega }-Y_{{}_{\omega }}+S_{\omega }+D_{\omega }$$

(6)

where $G_{\mathit{k}}$ and $G_{\omega }$ are the generations due to the mean velocity gradient, and $Y_{\mathit{k}}$ and $Y_{\omega }$ are the dissipations due to turbulence of k and $\omega$, respectively; $S_{\mathit{k}}$ and $S_{\omega }$ denote the moduli of the mean rate-of-strain tensors, while ${\sigma }_{\mathit{k}}$ and ${\sigma }_{\omega }$ present the diffusion rates of turbulent kinetic energy and specific dissipation rate, respectively; $D_{\omega }$ is the cross-diffusion term; ${\mu }_t$ is the dynamic turbulent viscosity.

2.2.3. Numerical Setup

To achieve the best comparison between the two models, the numerical study settings in OpenFOAM and Fluent were selected as similar (as shown in

Table 2. Setup of numerical models.

Model	Turbulent Model	Pressure Scheme	Velocity/Turbulence Scheme	Pressure-Velocity Coupling
OpenFOAM	$\mathit{k}-\omega SST$	second-order (GAMG)	second-order (smoothSolver)	PISO
FLUENT		second-order (PRESTO)	second-order (-)

). A computational mesh that resembles the laboratory experiment was created by OpenFOAM and converted to mesh input file of Ansys Fluent (Figure 2).

The domain of numerical study was set same as the laboratory experiment setup, which was 12.5 m in length and 0.6 m in width. To optimize the calculation results near the wall, the mesh was gradually set to refine toward the bottom. Additionally, the grid size in the scoured area, where the flow is supposed to be complicated, was also created smaller to enhance the results. Moreover, a study about the grid size and the time step convergence were conducted to select the appropriate parameters, as is always recommended [8,9]. This step helps to reduce the errors from the discretization. To get an optimal combination of the two parameters, the simulation results from three different grid sizes (accounting for the first cell from the wall) of $\Delta z$ = 1, 1.5, and 3 mm ($\Delta z/h_0$ = 4.7 × 10⁻³, 7.0 × 10⁻³, 1.4 × 10⁻²) and three different time steps increase from $\Delta t$ = 1.0 × 10⁻⁵ to 1.0 × 10⁻⁴ and 1.0 × 10⁻³ s, which are compared in Figure 3. Here, the horizontal velocities u_x were collected at 4.5 m downstream. The result indicates that the time step of Δt = 1.0 × 10⁻⁴ s and grid size of minimum Δz = 1.5 mm (close to the wall boundary) and maximum of 10 mm (near the atmosphere boundary) are sufficient to produce the results within 0.1% of the higher resolution cases.

All of OpenFOAM simulations in this present study were operated with the pisoFoam module, a transient solver for incompressible flow with many choices for turbulence simulation. The Gaussian integration with different interpolation schemes was applied for the spatial discretization of differential operators. The second-order linear interpolation was used for the gradient terms. The linear corrected interpolation method was applied for Laplacian terms. The limited linear interpolation was employed for the divergence terms. The GAMG solver with Gauss Seidel smoother was used for pressure, and the smooth Solver 177 with Gauss Seidel smoother was applied for velocity and turbulence. The pressure–velocity coupling was calculated by the Pressure Implicit with Splitting of Operators (PISO) method, which is an efficient approach to solve the Navier–Stokes equations in unsteady problems.

In Fluent simulations, the PRESTO scheme was used for pressure, while the second-order scheme was applied for velocity and turbulence. Similar to OpenFOAM, the PISO method was employed for the pressure–velocity coupling calculation. The under-relaxation factor of 0.5 was used for pressure, momentum, and turbulent quantities to avoid the divergence of numerical simulations and archive the full convergence.

3. Results and Discussion

3.1. Streamline

Flow patterns of case Q20 h120 d12 (Q = 0.020 m³/s, h₀ = 0.12 m, d₅₀ = 1.2 mm) in the stabilized/equilibrium scoured hole at 245 h later were revealed in the previous study by Park [18]. From the results, stream-wise flow velocity was faster in the vertically upper part (where z > initial bed elevation) than in the lower part (where z < initial bed elevation). In the lower part of upstream scour slope, a reverse and circulating movement of sediment particles that were relatively much slower than the depth-averaged flow velocity was captured (Figure 4).

The flow circulation was also presented in the first plot of Figure 5. Both OpenFOAM and Ansys Fluent showed their ability of capturing this property of the flow by producing a circulation flow right after the transition point of x/h₀ = 0, where the channel bed was set abruptly to change in the roughness. This flow separation phenomenon plays key role in moving or transporting the sediment particles of the channel bed and leads to the development of the scour hole.

As shown in the Figure 5, there is a backflow occurs near the bed in a region called deceleration zone. The size of flow circulation was measured to reach x/h₀ = 6.9 in the laboratory experiment case. However, the numerical results from both models show a little over-prediction, with the separation point located at x/h₀ = 9.7 in the OpenFOAM case and x/h₀ = 10.4 for the Ansys Fluent simulation, which are 1.4 and 1.5 times longer than the experimental case, respectively. This result suggests that both OpenFOAM and Ansys Fluent seem to estimate higher velocities of the flow in this circulation zone than the experimental data. Though both models fail in calculating the exact circulation size, they are in good agreement with each other, with the difference of only 6.7%.

3.2. Velocity

Distribution of the time-averaged stream-wise flow is plotted in Figure 6. Flow velocity was faster in the vertically upper part (z > initial bed elevation) than in the lower part (z < initial bed elevation) in the scoured hole. In the lower part of upstream scour slope, a reverse and circulating flow was relatively much slower than the primary flow in the channel (Figure 6). The distribution of stream-wise flow velocity caused a larger value of velocity gradient in the turbulent shear layer.

Additionally, while the above streamline results are needed to show the flow patterns, the velocity profiles are important for analyzing the flow properties. In this study, although the streamline visually presents the circulation zone, the figures cannot give the exact measurement data of that zone ending point. Therefore, in this case, the velocity profile is normally used to distinguish flow separation by determining the final point of backward velocity vectors. The comparison of numerical simulation results from OpenFOAM and Ansys Fluent with the laboratory experiment data is presented in Figure 7.

In Figure 7, the velocity profiles of the flow at x/h₀ = 0 (means right at the transition point), x/h₀ = 2.1, and x/h₀ = 4.2 are presented. First of all, as can be seen here, the simulation results from OpenFOAM and Ansys Fluent showed a very good agreement in predicting the flow velocity in the scour hole. The difference of their averaged results is only 0.0138 calculated by RMSE value. The velocity value predicted by Ansys Fluent is a bit smaller than that from OpenFOAM. Moreover, both OpenFOAM and Ansys Fluent clearly show that the minus velocity values represent the backward velocity vectors in the region near the wall, which cause the circulation flow. However, results from both of these models show some difference compared to the laboratory experiment data, with the overall RMSE values of 0.04 and 0.07 for OpenFOAM and Ansys Fluent, respectively. For the free stream flow above the fixed bed level (z/h₀ = 0), the models and experiment data fit well. However, the numerical simulation under-predicted the flow velocity inside the scour hole. Additionally, while the results right at the transition point x/h₀ = 0 (where the channel bed roughness changed abruptly) are quite similar with a logarithmic curve of velocity for all measurements, the difference increases when the flow goes downstream into the scour hole.

3.3. Turbulent Kinetic Energy

Values of turbulent kinetic energy were calculated and presented in Figure 8. Mostly, at 2 cm from the scoured hole, the values were estimated as excessively increased due to the larger values of velocity gradients and flow circulation. According to previous studies [20,21], high turbulence intensities are possible in decelerating turbulent flows, due to the formation of layers with great velocity gradients in abrupt expansions of water depth. Hoffmans and Booij [21] described turbulent mixing (or shear) layer in the scour hole with experimental data and numerical approaches. As they proposed, re-circulation (or reverse movement) of flow in the upstream scour slope occurs near the bed, and a mixing layer develops between the transient flow and the recirculating flow [20,21].

The non-dimensional turbulent kinetic energy profiles produced by the two numerical models at several cross-sections are plotted as in Figure 9 in comparison with the laboratory experimental data. For a turbulent flow, such as the flow of water running from a smooth bed through a sand bed channel in this study, the turbulent kinetic energy k is normally used to investigate the turbulence characteristic of that flow, since it is defined as the mean kinetic energy per unit mass associated with eddies in turbulent flow.

As can be seen here, both OpenFOAM and Ansys Fluent well captured the happenings of the high turbulent region in the scour hole. The turbulent kinetic energy value increases dramatically from the layer close to the bed to reach the maximum in the center of the circulation flow, and then decreases in the free stream upper region, to finally reach a constant value, which is similar to the value of k at the transition point of x/h_{0 Figute}= 0. This is due to the change in horizontal, vertical, and lateral velocity fluctuations in these regions. While velocities near the bed are relatively low due to the vertical flow separation, the circulation causes back-flow and induces the formation of small eddy scales in the scour hole that increase the velocity fluctuation, and therefore increases the turbulent kinetic energy of this zone. At the beginning of the scour hole and in the free stream region above the fixed bed level, the results of k from both models are almost identical for both model results. However, they show some difference inside the circulation zone. In this region, the predicted turbulent kinetic energy from Ansys Fluent model is smaller than that from OpenFOAM. Overall, the difference between these two results is RMSE = 2.5 × 10⁻⁴, but it reaches RMSE = 3.5 × 10⁻⁴ in the scour hole. Moreover, though both the OpenFOAM and Ansys Fluent models can capture the overall behavior of turbulent kinetic energy profile of the flow, they are still not so accurate when compared with the experimental data. In this study, the $\mathit{k}-\omega SST$, which was developed to simulate both the low and high Reynolds number or close and far from the wall, was used in all simulations. However, overall, the value of differences compared to laboratory experiment case reached RMSE = 0.0058 and 0.0062 for OpenFOAM and Ansys Fluent, respectively.

4. Conclusions

This study focused on comparing the performance of two numerical models—the open-source package OpenFOAM and the commercial software Ansys Fluent—when simulating the near-wall flow, or the turbulent flow in a scour hole to be specific. Both codes were set up similarly in conditions and schemes. The results of flow pattern streamline, velocity profiles, and turbulent kinetic energy were collected and compared to previous laboratory experimental data. Overall, both models showed a good ability of capturing the flow behavior, such as circulation flow that happens in a scour hole, as well as the development of velocity and turbulence. Although they cannot perfectly archive the results as the laboratory experimental data, the performances by the two models are reasonable, and especially quite similar for the first two categories, while only showing a small difference in the last one.

In conclusion, the calculation of OpenFOAM is better than Ansys Fluent for the flow in a scour hole near the wall as shown in this study. According to the energy cascade principle, the higher turbulent kinetic energy produced by OpenFOAM can result in creating more vortices and therefore can lead to the higher mixing in the scour hole region than by Ansys Fluent. This result can be a useful note when make a choice between the models.