VOF-Based Simulation of Turbulent Air-Water Flow over Gravel Beds in Open Channels

1 Introduction

In engineering applications, streams and rivers are considered the primary examples of open channel flows. Understanding how the general characteristics of open channels depend on cross-sectional geometry is essential, as it allows engineers to evaluate flow capacity across different shapes [1]. In mountainous regions, gravel transport is common, and numerous studies have investigated how gravel beds influence flow behavior [2].

Surface or bed roughness directly affects the Reynolds number (RN) by increasing its value. To account for roughness effects in flow channels, scalable wall functions must be applied to various mesh types to assess how roughness influences fluid dynamics and RN [3,4]. Free surface effects can be captured using the Volume of Fluid (VOF) model in computational fluid dynamics (CFD) simulations, with the key element in open channel models being the air-water interface. Evaluating RN is crucial for understanding flow performance in channels with rough beds [4,5]. Turbulence near the free surface is particularly important because vorticity, momentum, heat, and mass transfer are largely governed by surface dynamics. Free surface behavior presents significant challenges in rough bed flows, and knowledge of two-phase interactions remains limited [6].

Roughness is not confined to the bed; it also manifests at the free surface through surface tension between phases, and its magnitude depends on the roughness height [7]. Understanding flow development in rough channels with turbulent near-surface interactions is complex. To capture these effects, the free surface was modeled using the VOF approach coupled with the standard k-ε turbulence model. Roughness impacts velocity magnitude and alters the direction of flow vectors over rough surfaces [1]. Gravel roughness in rectangular channels significantly affects flow resistance and RN, as well as transitions between laminar, transitional, and turbulent regimes [8]. In open channels, the critical Reynolds number marks the transition from laminar to turbulent flow, and surface roughness increases RN while influencing turbulent characteristics [8].

Open channel flows are characterized by a free surface exposed to the atmosphere. Key parameters such as the friction coefficient and shear velocity are essential for quantifying flow resistance; however, accurately estimating these factors remains challenging. Traditional open channel equations often provide imprecise results due to difficulties in determining energy slopes and maintaining steady, uniform flow in human-made channels [9]. Channels can be natural or artificial, with bed roughness manipulated in the latter to control flow. The friction coefficient is a critical parameter for predicting flow velocity and ensuring accurate design. While CFD research often focuses on circular pipe flows, theoretical approaches remain important for open channel experiments due to computational complexity [10].

Recent advances in numerical modeling have demonstrated the utility of CFD in simulating complex open channel and rough bed flows, validated against laboratory data. For example, Yu et al. [11] conducted a comparative 2D CFD study of turbulent flow in open channels with emerged vegetation using multiple RANS turbulence closures, highlighting the influence of model choice on turbulence predictions and mean velocity structures while confirming the applicability of standard k-ε formulations. Xie et al. [12] employed OpenFOAM to simulate rough bed channels, showing strong agreement with experiments and capturing roughness effects on turbulent structures. Large-eddy simulation has been applied to free surface gravel bed flows to better resolve turbulence and bed interaction phenomena beyond RANS capabilities, providing insight into near-bed vortical behavior and free surface dynamics [13]. Additionally, direct numerical simulations of turbulent open channel flows at moderately high Reynolds numbers have documented detailed turbulence statistics and free surface effects, demonstrating the continuing development of high-fidelity computational approaches [14,15].

Determining roughness parameters remains challenging and often requires field measurements or laboratory experiments [3]. The friction factor is widely recognized as a critical variable for predicting flow behavior [16]. In laminar open channel flows, channel shape significantly affects the relationship between the friction factor and RN for non-Newtonian fluids [7]. Bed roughness enhances turbulence, typically reducing velocity near the rough surface while increasing it at the free surface [17]. Open channel studies are used to forecast flow patterns for a variety of scenarios, including stormwater, drainage, sediment transport, and suspended particle tracking in turbulent flows [2,18]. The impact of roughness on flow characteristics varies with channel design and configuration [2].

In this study, experimental investigations were conducted to assess flow resistance in a rectangular open channel with uniformly distributed gravel. Complementary numerical simulations were performed in ANSYS using VOF coupled with the k-ε turbulence model. The objectives were to evaluate the effect of designed roughness on flow, quantify the friction factor-RN relationship, and compare experimental and numerical results for validation. The observed f-RN relationship aligns with previous studies, providing insights into the hydraulic performance of gravel-lined channels. These findings have practical implications for improved flow prediction, optimized channel design, flood risk assessment, sediment transport management, and environmental considerations in hydraulic engineering [1,12,13].

2 Materials and Methods

2.1 Experimental and Numerical Setup

A rectangular flume with a gravel bed was used for the experiments. The flume measured 10 m in length and 0.3 m in width. All experiments were conducted in the Hydraulics Engineering Laboratory of the Department of Civil Engineering, University of Engineering and Technology (UET), Lahore. The channel was glass-sided with a depth of 0.75 m and featured an adjustable tailgate at the downstream end to regulate flow depth.

Discharge was supplied to the channel via a centrifugal pump, passing through a honeycomb flow straightener to ensure uniform and smooth distribution across the channel width. At the downstream end, water entered a sediment tank where sediments were trapped, and the flow then discharged from the main channel. A plan view of the laboratory flume and experimental setup is shown in Fig. 1.

Flow depth measurements were obtained using a hook gauge mounted on top of the flume, allowing precise determination of water surface elevation at multiple locations along the channel. These measurements were used to compute hydraulic parameters such as flow area (A) and hydraulic radius (Rh). Flow velocity (V) was calculated indirectly from the measured discharge (Q) and corresponding flow area, following standard open-channel flow procedures.

Figure 1: Experimental setup: (a) rectangular flume used for open channel experiments, (b) gravel bed segments (1 m each) installed in the flume, (c) 2D numerical model of the channel with the gravel bed surface.

Since the flow was non-uniform along the channel, the bed slope (So) could not be reliably used to estimate Manning’s “n”. The gravel bed (Fig. 1b) was prepared on a horizontal wooden plank measuring 1 m × 0.3 m, following Bazin’s channel bed characteristics [1]. The gravel layer was then installed in the flume to reproduce realistic roughness conditions. A corresponding 2D numerical model (Fig. 1c) was generated in AutoCAD to incorporate gravel bed effects for open channel flow simulations using ANSYS Fluent.

In the VOF approach, a volume fraction α describes the local presence of one phase (e.g., liquid) in a multiphase system. The governing equations (Eqs. (1)–(5)) consist of the conservation of volume fraction, the incompressible continuity equation, and the momentum equation (Navier-Stokes with surface force). ∂α∂t+∇·αu=0(1) ∇·u=0(2) ρα∂u∂t+u⋅∇u=−∇p+∇·μα∇u+∇uT+ραg+fσ(3) fσ=σk∇α(4) k=−∇·∇α∇α(5) where σ is the surface tension coefficient and k is the interface curvature.

In numerical simulations, meshing plays a critical role in accurately predicting fluid behavior. A refined mesh was applied in regions of interest, particularly near the gravel bed and the free surface. Triangular elements with a characteristic length of 0.001 m were employed, resulting in 39,717 cells, 60,225 faces, and 20,509 nodes, following the grid refinement study shown in Fig. 2. The mesh quality was evaluated using skewness and orthogonality criteria to ensure numerical stability and solution reliability. No highly skewed or poorly orthogonal elements were detected. Additional refinement near the gravel bed and free-surface regions enabled accurate resolution of velocity gradients and turbulence effects. Smooth convergence behavior and stable residual reduction confirmed the adequacy of the adopted mesh for the present simulations.

Figure 2: Grid refinement against friction factor (f).

The multiphase simulations involved air and water, with gravitational acceleration of 9.81 m/s2 applied in the vertical direction and a free surface condition specified. Numerical simulations were performed in ANSYS Fluent (version 2020). A steady-state formulation was adopted to investigate mean flow behavior and resistance characteristics under fully developed turbulent flow conditions. The VOF model captured air-water interactions, and turbulence was modeled using the standard k-ε model. Pressure-velocity coupling was handled using the SIMPLE algorithm, and second-order upwind schemes were applied to momentum and turbulence equations for improved accuracy. Convergence was ensured by reducing residuals below 10−5 for continuity and momentum and below 10−6 for turbulence quantities, while monitoring the stabilization of key hydraulic parameters such as velocity and discharge.

The experimental and numerical workflow employed in this study is summarized in the flow chart shown in Fig. 3.

Figure 3: Flow chart illustrating the methodology of the experimental and numerical investigation.

While applying a free surface boundary condition, the implicit method was used, and a no-slip condition was enforced at all channel walls [5,14]. Under the no-slip condition, the wall is assumed stationary, and the fluid velocity at the wall is zero [19]. The hydraulic grade line (HGL) provides information on pressure at various sections, which is essential for the structural design of channels. For steady flow, the flow depth at a given cross-section remains constant (Eq. (6)). Prismatic and man-made channels typically exhibit average steady flow conditions [1,20].

In uniform flow, the flow depth, velocity, acceleration, slope, and channel cross-section remain constant along the distance between two sections of the channel (Eq. (7)) [15].

In steady-uniform flow, characteristics such as discharge, velocity, and flow area remain constant over both time and distance. In contrast, steady-non-uniform flow describes a regime where these properties remain consistent over time but vary along the channel length [19]. Unsteady-uniform flow is characterized by fluid properties that remain consistent along a specific distance but vary over time. Conversely, unsteady-non-uniform flow refers to a regime in which fluid properties vary with both time and distance [1].

Although the flow regime depends on the roughness configuration of the channel, hydraulic dynamics are also governed by fluid properties. Flow is categorized as laminar or turbulent based on the relative influence of viscous and inertial forces [20]. According to RN, flow is considered laminar when RN < 500. For pipe flow, the laminar-to-turbulent transition occurs approximately between RN = 2000 and 50,000, while in open channels, the transition range extends from RN = 500 to 12,500 [21]. The relationship between viscous and inertial forces can be expressed in terms of RN (Eq. (8)). RN=VLv(8) where V is the average flow velocity (m/s) and L is the characteristic length (m). For open channels, L is taken as the hydraulic radius (Rh) of the conduit, while for pipe flow, it is equal to the pipe diameter. The kinematic viscosity of water is denoted by ν (m2/s). The fluid is assumed to be ideal, incompressible, and non-viscous. The coefficient C is a function of flow forces and the relative roughness of the channel (Eq. (9)) [22,23]. V=CRS(9) Rh=AreaWettedPerimeter(10) where C is Chezy’s roughness coefficient, S is the channel bed slope, and Rh is the hydraulic radius (Eq. (10)). This formula is particularly suitable for lined channels. In open channel flows, Manning’s n is determined using Eq. (10). Factors such as channel inclination, vegetation, gravel, sediments, and flow velocity contribute to increased disturbances in the flow. Determining the value of n is challenging, as each type of open channel exhibits a different value [8,17].

According to Manning’s formulation, the horizontal component of gravitational force is balanced by the shear resistance of the bed [22,24]. The Manning equation provides an empirical relationship between flow area, channel slope, and flow velocity, with Manning’s coefficient representing the friction exerted by channel bed roughness (Eq. (11)). As channel roughness increases, Manning’s n correspondingly increases [1,15,18]. Q=1nARh23 S12(11) where n is Manning’s roughness coefficient, and S ≈ So is the channel bed slope. This formula is preferred for unlined channels. The different coefficients, including n, C, and f, were calculated using the following relationships (Eqs. (12) and (13)) [12,13].

2.2 Validation of Numerical Results

The Colebrook equation provides a computational method for determining f in fully developed pipe flows or ducts, and it can also be applied to approximate open channel flow. This equation is primarily used for turbulent flows, including high-intensity fluid regimes approaching their upper limits. The mathematical formulation for f is given in Eq. (14) [14,25].

The second method used for numerical validation is the Haaland approximation [16,26]. This method is also applicable to turbulent flow conditions, although turbulent flows are restricted to a specific RN range. The RN for which this approximation is valid extends from a minimum of 4000 to a maximum of 5 × 108. The corresponding mathematical expression is given in Eq. (15). In addition to empirical and semi-empirical equations, data-driven approaches such as artificial neural networks have also been applied to predict friction factor behavior in open channel flow with satisfactory accuracy [27,28].

2.3 Statistical Indicators

Several statistical tests were performed to evaluate the performance of the experimental and numerical results, as summarized in Table 1. Here, Ni represents the numerical values, while O denotes the observed data. The symbol O ¯ indicates the mean value of the observed data.

3 Results and Discussions

The observed inverse relationship between f and RN (Fig. 4) reflects the progressive dominance of inertial forces over viscous effects as flow intensity increases. At lower RN, flow resistance is strongly influenced by bed roughness and near-bed turbulence, resulting in higher f. As RN increases, the relative effect of roughness diminishes, and the flow becomes fully turbulent, reducing frictional resistance. This behavior is consistent with classical hydraulic theory and well-established empirical formulations, including Bazin’s curves and the Colebrook and Moody relationships [1,14]. Comparable inverse trends between f and RN have also been reported for vegetated and rough open channel flows under turbulent conditions [26,29].

Similarly, the positive relationship between C and V (Fig. 5) indicates an increase in flow conveyance efficiency with higher velocity, as increased momentum allows the flow to overcome bed-induced resistance more effectively. The agreement between experimental observations and numerical predictions demonstrates that the adopted modeling approach can capture the dominant hydraulic mechanisms governing gravel-bed open channel flow, despite the inherent simplifications of depth-averaged and two-dimensional simulations.

As discharge (Q) increases, the flow area (A) also increases, which interacts with more roughness, ultimately raising Manning’s n (Table 2). The relationship between V and n for a gravel bed was found to be exponential and inverse (Table 2 and Table 3), with higher velocities corresponding to lower n values. This occurs because increased velocity reduces flow depth, causing the flow to approach supercritical conditions, as evident from the specific energy diagrams [10].

During experiments with the gravel bed, several interesting behaviors were observed. The relationship between f and RN was linear and inverse (Fig. 4), showing strong correlation with Bazin’s curve [1]. Notably, the C value corresponding to the 18th experimental observation (Table 2) was unusually high compared to the rest of the dataset. This outlier is attributed to measurement uncertainty in flow depth at different channel locations during that run. The resulting very small energy slope produced an unrealistically high C value, which does not represent the actual hydraulic behavior of the gravel-bed channel. This data point is retained for completeness and transparency, but it is not considered representative of the overall flow trend.

Figure 4: Relationship between the f and RN from experimental measurements in a gravel-bed open channel.

Figure 5: Relationship between C and V from experimental measurements in a gravel-bed open channel.

Fig. 4 illustrates a negative correlation between f and RN. As RN increases, indicating a transition from laminar to turbulent flow, f decreases. This behavior is typical in fluid dynamics, as turbulent flow, associated with higher RN, generally produces lower frictional losses. The fitted trend line highlights this inverse relationship, emphasizing how increased turbulence reduces f.

The C value for the 18th gravel-bed measurement was unusually high, recorded as 100. This anomaly is attributed to measurement errors in flow depth at multiple channel locations. Specifically, the energy line slope for this reading was extremely low, suggesting negligible head loss, which is physically unrealistic for a gravel-bed channel.

Similarly, Fig. 5 shows a positive correlation between C and V. As V increases, C also increases, as indicated by the upward trend of the fitted line. This relationship appears to be linear, with the majority of data points following a consistent increasing trend, suggesting that higher flow velocities correspond to higher Chezy coefficients.

The numerical results were evaluated using ANSYS Fluent. In the simulations, the free surface was defined by specifying the fluid depth in the modeled open-flow channel. Air was assigned as phase-1, and water as phase-2, allowing computation of f, C, and RN. The effects of gravel-bed roughness were analyzed within the phase-2 region, while phase-1 served as the boundary layer for the free surface.

For comparison with experimental data, V was plotted against C, as shown in Fig. 6. The numerical results indicate greater stability in C values relative to the experimental measurements. Fig. 6 demonstrates a positive linear relationship between C and V, with both variables increasing consistently. As V increases, C also tends to rise, as indicated by the upward-sloping trend line. However, the data points exhibit some variability, particularly at higher V values, suggesting that while a general increasing trend exists, deviations occur. Overall, the plot indicates that C is directly influenced by V, though with a certain degree of fluctuation.

Figure 6: Relationship between C and V obtained from numerical simulations for a gravel-bed open channel.

Fig. 7 illustrates a negative relationship between f and RN. As RN increases, f decreases, as shown by the downward-sloping trend line. This behavior aligns with fluid dynamics principles, where higher RN, typically associated with more turbulent flow, corresponds to lower frictional losses. The plot confirms a clear inverse correlation, indicating that as flow becomes more turbulent, the friction factor diminishes.

Figure 7: Relationship between the f and RN obtained from numerical simulations for a gravel-bed open channel.

Table 4 presents statistical indicators for three variables: f, RN, and C. The correlation coefficient (r) indicates that RN exhibits the strongest positive correlation, followed by f and C. The percentage error (PE) shows that RN has the lowest error, whereas C and f display higher errors. In terms of mean absolute error (MAE), RN has the highest value, reflecting larger prediction deviations compared to f and C. Similarly, the root mean square error (RMSE) is greatest for RN, suggesting that RN experiences the largest prediction discrepancies, followed by C and f. Overall, RN demonstrates strong correlations but also exhibits substantial prediction errors.

Computational studies inherently involve some degree of error, as numerical models attempt to replicate real-world physical systems but cannot achieve exact correspondence with experimental conditions. The differences between the experimental and numerical results are illustrated in Fig. 8a,b.

Figure 8: Comparison between experimental and numerical results; (a) f versus RN; (b) C versus V.

The observed discrepancies arise from modeling errors, which occur when converting a real physical system into a numerical model. These errors represent deviations between experimental and numerical results. In practice, it is not possible to replicate the exact geometry of the gravel bed within the modeling software; instead, the numerical model approximates the bed based on available simulation data.

The results generally follow the expected trends: lower f values correspond to higher RN, consistent with Bazin’s curve, which shows that as RN increases, friction decreases. Similarly, the relationship between C and V remains direct, while the f versus RN relationship is inverse. However, due to the inherent complexity of the flow and the approximations in the numerical model, the results cannot perfectly match the exact physical behavior.

Validation of Experimental and Numerical Results

Effective presentation of results is crucial for conveying the findings clearly and ensuring a proper understanding of the investigated problem. In the numerical simulations, contour plots are used to visualize the flow patterns and solution distribution. These plots serve as a conceptual framework for the model, guiding the interpretation of results during the post-processing stage and facilitating accurate analysis of flow behavior.

The velocity magnitude contours indicate that the fluid follows the prescribed flow direction and boundary conditions, as illustrated in Fig. 9. Accurate interpretation of the VOF model requires careful examination of these contours to ensure consistency with the defined computational domain. In the VOF contours, the red regions represent the water phase, while the blue regions correspond to the air phase, as shown in Fig. 10.

Figure 9: Velocity magnitude contour obtained from numerical simulations.

Figure 10: VOF contour for the water phase obtained from numerical simulations.

In the open channel flume, the low-density fluid is air, designated as phase-1 in the simulations to represent the upper region of the flume. In the VOF contours, red indicates the air phase, while blue represents the water phase. The VOF distribution for air is shown in Fig. 11.

Figure 11: VOF air contour for the air phase obtained from numerical simulations.

Investigating flow in an open channel is a challenging task, particularly for non-prismatic channels. To assess the accuracy of the numerical results, validation was performed using two approximation techniques: the Colebrook free-surface equation and the Colebrook-Moody approximation (Fig. 12). Both methods confirmed the f-RN relationship, demonstrating the expected decrease in f as RN increases.

Figure 12: Relationship between f and RN for the gravel-bed open channel using: (a) Colebrook free-surface equation, (b) Colebrook approximation by Moody’s method.

4 Conclusions

This study combined laboratory experiments with numerical modeling to better understand flow resistance in a gravel-bed open channel and to examine how classical hydraulic relationships perform when applied using modern CFD tools.

• 1.The main objective was to investigate the relationship between the friction factor and Reynolds number in a gravel-bed open channel and to assess how bed roughness influences flow resistance under turbulent conditions.
• 2.Experimental tests were carried out in a rectangular flume with a uniformly distributed gravel bed, where flow depths were measured using a hook gauge and hydraulic parameters were obtained from measured discharge and channel geometry.
• 3.A two-dimensional CFD model was developed using the Volume of Fluid (VOF) approach coupled with the standard k-ε turbulence model to represent free-surface flow and turbulence effects.
• 4.The results consistently showed that the friction factor decreases as the Reynolds number increases, indicating that at higher flow intensities the influence of bed roughness becomes less dominant.
• 5.A clear increasing trend was also observed between flow velocity and Chezy’s coefficient, suggesting improved flow conveyance with increasing velocity.
• 6.The numerical predictions showed good agreement with the experimental data, with a strong correlation for the Reynolds number (r ≈ 0.99) and reasonable correlations for the friction factor (r ≈ 0.75) and Chezy’s coefficient (r ≈ 0.71).
• 7.The average percentage differences between numerical and experimental results were about 17% for Reynolds number, 25% for the friction factor, and 12% for Chezy’s coefficient, which are acceptable considering experimental uncertainty and modeling simplifications.
• 8.Validation using the Colebrook equation and Moody’s approximation further supported the observed inverse relationship between friction factor and Reynolds number.
• 9.The overall findings indicate that the combined experimental-numerical approach used in this study is suitable for analyzing gravel-bed open channel flows and can support practical hydraulic assessments and preliminary channel design.
• 10.Future work could extend this research by applying three-dimensional simulations, testing different gravel sizes and arrangements, and incorporating more advanced turbulence models to better resolve local flow structures.

In summary, the study confirms that classical resistance relationships such as those proposed by Bazin, Colebrook, and Moody remain applicable when supported by CFD modeling, providing a reliable framework for studying gravel-bed open channel flows.