Three-Dimensional Numerical Simulation of a Two-Phase Supercritical Open Channel Junction Flow

Abstract

This study investigates the computational fluid dynamics (CFD) modeling of supercritical open channel junction flow using two different turbulence models: k-ω shear stress transport (SST) and k-ω SST scale-adaptive simulation (SAS), in conjunction with Volume of Fluid (VOF) and mixture multiphase models. The efficacy of these models in predicting the intricate free surface fluctuation and free surface elevation in a supercritical junction is evaluated through a comprehensive analysis of time-averaged free surface data obtained from CFD simulations and Light Detection and Ranging (LIDAR) measurements. The dimensionless Reynolds (Re) and Froude (Fr) numbers of the investigated scenario were Fr = 9 and Re = 5.1 × 10⁴ for the main channel, and Fr = 6 and Re = 3.3 × 10⁴ for the side channel. The results of the analysis demonstrated a satisfactory level of agreement with the experimental data. However, certain limitations associated with both CFD and LIDAR were identified. Specifically, the CFD performance was limited by the model’s incapacity to consider small-scale turbulent effects and to model air bubbles smaller than the cell size while the LIDAR measurements were limited by instrument range, inability to provide insight into what is happening below the water surface, and blind spots. Nonetheless, the k-ω SST turbulent model with the VOF multiphase model most closely matched the LIDAR results.

1. Introduction

Modeling supercritical open channel flow with computational fluid dynamics (CFD) presents a formidable challenge due to its inherent complexity. The flow in supercritical conditions is characterized by shallow depth, high velocity exceeding the speed of surface gravity waves, and intricate free surface dynamics, making it prone to abrupt changes and complex flow patterns. This flow regime is established when the Froude number (Fr = u/(gL)^0.5; u is local flow velocity, g is the gravitational acceleration, and L is characteristic length or depth of flow, a dimensionless parameter representing the ratio of flow inertia to gravity force, which exceeds unity. Additionally, the high Reynolds (Re = ρuL/μ; ρ is the fluid density, u is the flow speed, L is a characteristic length, μ is the dynamic viscosity) and the Froude numbers associated with supercritical flows result in non-homogeneous and non-stationary behaviors. The dynamics involve the interplay of various factors such as turbulence, air entrainment, and the interaction of water and air at the free surface [1], further complicating the modeling process. The presence of junctions in open channels introduces secondary flows, increased velocities, pulsations, and air bubble entrainment, amplifying the intricacies of the system [2,3]. The flow complexity, combined with the need for high spatiotemporal resolution, poses computational challenges and demands careful consideration in selecting appropriate turbulence models and numerical schemes.

Recent events, such as the hundred-year floods in Slovenia in 2023 [4], have highlighted the importance of the planning and construction of infrastructure to enhance flood resilience and safety measures. The accurate modeling of hydraulic phenomena in supercritical open channel flows is of great importance for the design of hydraulic structures such as high-speed outlets, sewage, and storm drain systems, hydroelectric power plant spillways, drop shafts, etc. [5]. The knowledge of hydraulic phenomena is well advanced in subcritical flows but is still limited in supercritical flows and transitions to multiphase flows [5]. Free-surface flows in hydraulic structures are usually turbulent and characterized by air bubble entrainment due to surface deformation and high local shear stresses [2,6]. Free-surface flows with high Reynolds and Froude numbers (Re > 10⁴ and Fr > 3, respectively) are highly complex, non-homogeneous, and non-stationary with rough water surface and large velocity variations [7,8,9]. Interestingly, such high Reynolds and Froude number flows may form large standing and fluctuating quasi-periodic waves at junctions [10]. Standing waves at the junction can also form concave surfaces that are difficult to capture using the available measurement methods, often preventing the clear determination of the free surface position and height.

The highly dynamic and multiphase nature of junction flows precludes the use of simpler point measurement methods such as manometers [11], wave probes [12,13], point gauges [14], or electromagnetic/ultrasonic devices [15,16] when accurate modeling is desired. The accurate spatial and temporal measurements of complex three-dimensional free surface flows in natural and man-made hydraulic structures are only viable by high-resolution non-contact methods such as LIDAR [17,18], laser triangulation [19] or photogrammetric reconstruction from cameras [20] with overlapping fields of view [21]. While photogrammetric methods are still under development, laser-based methods [10,22,23] have been proven as reasonably accurate in the measurement of such hydraulic phenomena. Considering the progress in experimental methods with a high spatiotemporal resolution, the development of related numerical simulations by CFD is the next logical step in the modeling of junction flows. CFD simulations validated by experimental data may be used to determine several flow-related quantities (velocity, pressure, and air void fraction) with a spatiotemporal resolution unchallenged by experimental methods. This enables a detailed analysis of flow characteristics, offering insights that might be challenging to achieve solely through experimental approaches. Furthermore, an important advantage of CFD simulations over experimental research is the capability for adaptation to different junction geometries without the requirement of constructing a new physical model, given that flow conditions are not extrapolated too far from the experimentally evaluated range.

Feurich and Olsen [24] employed a three-dimensional numerical model (STAR-CCM+) to compute the free surface in several supercritical flow cases, successfully reproducing oblique standing waves and the location of the free water surface. Nevertheless, it is important to emphasize that this specific type of investigation has yet to be conducted for supercritical junction flows. While extensive attention and examination have been dedicated to the three-dimensional numerical studies of subcritical flow structures in channel junctions [25,26,27,28], the exploration of supercritical flow scenarios remains relatively underexplored and less understood.

For the modeling of violent free surface flows, two commonly utilized approaches are the classical finite volume method and the incompressible Smoothed Particle Hydrodynamics (SPH)-based method [29]. Supercritical junction flows have been effectively simulated using a one-dimensional SPH approach, as demonstrated by [30]. This method employs one-dimensional (1-D) Saint-Venant equations, commonly referred to as the shallow water equations [31], which are derived from depth-integrating the Navier–Stokes equation. Conversely, the 3D numerical study of flow structures in a channel junction with high Froude numbers imposes many challenges for the SPH-based methods as they struggle to predict highly aerated flows, rendering them unsuitable for this type of analysis.

Regarding the finite volume approach, different surface treating methods have been developed for simulating three-dimensional phase interactions: the rigid lid, Volume of Fluid (VOF), mixture, and user-specified models [32,33]. However, most studies are restricted to the single-phase approach or flows with low Froude numbers [25,32,34,35]. Notable finite volume-based studies investigated the various aspects of junction flows, such as free surface profiles around the junction [36], the effect of different junction angles on velocity profiles [37], water level and velocity profiles at 90° junction [32], and flow patterns around a vertical circular pier after the junction [38]. These studies, once again, primarily focused on scenarios with low inlet velocity magnitudes and, consequently, low Froude numbers (Fr < 1), resulting in relatively “flat” water surface with minimal air entrainment. In these studies, this allowed for using rigid lid boundary conditions to capture the water surface. Specifically, the rigid lid method assumes a fixed surface and can lead to substantial errors when the actual surface varies. Furthermore, neglecting the behavior of turbulent structures near the water surface, particularly in simulations involving open channels with rough beds, this circumstance can give rise to problematic outcomes [39]. The choice of the correct free surface modeling method plays an important role as it greatly impacts simulation accuracy in open channel junction flows. VOF and mixture methods offer improved accuracy but still have limitations when simulating shallow water flows.

Numerically simulating flows of this nature presents formidable challenges, as it involves the intricate interplay of turbulence and multiphase complexities [40]. Prandtl’s secondary currents add another layer to the intricacy, attributed to turbulence anisotropy [41]. Brown and Crookston [42] discuss the numerical simulation of complex flow patterns in supercritical flows acknowledging challenges in instances where numerical estimations for flow depth can disagree by more than 50%. Contributions from Jia et al. [43] and Song et al. [44] further contribute to the understanding of large-scale vortices’ evolution and recirculation zones, respectively, both critical components in the temporal and spatial exchange of mass and momentum in supercritical open channel flows. Nasif et al. [40] attributed the three-dimensional nature of the turbulent flow and the inherent limitations and high cost of conventional experimental devices, which do not offer quantitative isosurfaces of velocity, pressure, and vorticity as a numerical approach.

Overall, the simulation of supercritical junction flows introduces numerous challenges due to the complexity of the phenomena. These challenges encompass the occurrence of secondary flows, elevated velocities, air entrainment, and increased computational demands due to the necessity of smaller time steps and cell sizes which lead to prolonged calculation time and increased simulation cost. Consequently, accurately simulating such flows requires careful consideration. Acknowledging the intricate complexities, our study investigates the plausibility of such a modeling approach.

Our primary goal has been to determine a feasible combination of multiphase flow and turbulence models to accurately simulate supercritical junction flows with standing waves, which to our knowledge has not been previously accomplished. For this purpose, different RANS k-ω turbulence models were tested with VOF and mixture multiphase models to simulate the highly complex and unstable flow phenomenon at hand.

As the study is pioneering in its examination of supercritical open-channel flows using CFD, our emphasis was on a preliminary assessment rather than an exhaustive analysis of flow dynamics. Our investigation prioritized evaluating the performance of VOF and mixture two-phase models in conjunction with k-ω SST turbulent model in predicting the flow surface and determining which model aligns most closely with the experimental data obtained through Lidar measurements. Additionally, our objective was to assess if the simulation could be further enhanced with the Scale Adaptive Simulation (SAS) model, particularly concerning the k-ω SST mixture model.

2. Materials and Methods

2.1. Simulated Scenario and Reference Experimental Setup

The geometry of the CFD simulation domain was designed to resemble a 90° T-junction junction that had been previously used in hydraulic experiments [10,23] where the topography of the free water surface, namely the standing wave forming at the junction, was measured by LIDAR. In the present study, the existing LIDAR measurements will be utilized for the validation of the CFD simulations under identical boundary conditions.

The reference experimental setup comprised a 90° junction of two sharp-edged horizontal channels with a 0.5 m × 0.5 m square cross section. The main and side channel length upstream of the junction was 1 m, while the total length of the main channel was 6 m. All the channel surfaces were made from glass and with a minimum number of joints to reduce hydraulic roughness and its effect on flow conditions as much as possible. Each channel was supplied with water via a separate pipeline equipped with valves and an electromagnetic flow meter for accurate flow rate setting. Such a setup allowed for the generation of supercritical flows (Froude number up to 20) causing the formation of standing waves with complex geometry and exceeding the channel depth. Desired inflows were provided through a pressure vessel with adjustable height openings, and a free outflow was installed at the end of the main channel. The experimental setup was equipped with rails and an equipment platform for the mounting and precise positioning of the LIDAR device, as shown in Figure 1.

The observed junction flow was turbulent and highly aerated, but the height and shape of the standing wave did not vary significantly over time. This allowed for the time-averaged laser scanning measurements of the free surface topography. The LIDAR device (SICK LMS400, SICK AG, Waldkirch, Germany) was moved along the main channel in 100 mm increments, measuring the average elevation profile of the free water surface at each incremental position along the x-axis. Then, a 3D surface of the standing wave and surrounding flow areas was reconstructed from the averaged transversal free water surface profiles. The LIDAR device operated at a visible red-light wavelength (λ = 650 nm) and was configured with a line scan frequency of 270 Hz and an angular resolution of 0.2°. Each scan line included 350 measurement points within an angular range of 70°. To reconstruct the average transverse profile of the free water surface, a total of 6000 scan lines were employed. A detailed description of the post-processing techniques applied to the LIDAR data can be found in previously published papers [10,23].

The accuracy/uncertainty of the experimental measuring method used depends upon the water flow characteristics, including the entrainment of air bubbles in the water, water drops above the water surface, and water surface fluctuations. Given the complex free water surface topography, state-of-the-art methods face significant limitations in capturing hydraulic phenomena. The estimated measurement error is derived from the measurement uncertainty of the LIDAR scanner, as specified by the manufacturer, and the outcomes obtained through measurement filtering, as detailed in the study by Rak et al. [10]. The measurement uncertainty varies from ±5 mm in the locations with minimal vertical dynamics to ±10 mm in the locations with pronounced vertical surface dynamics [45]. In worst-case scenarios, this deviation may extend to values up to ±25 mm [46].

It is crucial to note that the data acquired in the near-wall region, within a 50 mm distance from the wall, raise concerns about reliability due to light reflection from the wall and should not be dismissed. Nevertheless, we validated the laser scanner measurements through three approaches. Firstly, since the traditional manual measurement method does not allow for the measurement of water surface dynamics, we visually recorded the local water level to gauge the average water surface level and fluctuations in individual points. This was accomplished using a thin slice ruler inserted in the water and a high-speed Casio EX-F1 camera with an image acquisition frequency of 300 fps. The influence of the ruler on the turbulent two-phase flow was deemed negligible [10]. Secondly, we compared the results of laser scanning measurements with the data obtained through laser triangulation [46]. Thirdly, we employed a photogrammetric method [20]. All the comparisons revealed a substantial agreement among the results of the different methods, though some discrepancies arose from variations in the working principles of the measuring methods. Regardless, the results affirm the aforementioned measuring uncertainty.

Our CFD simulation was based on a scenario (Figure 2) where the inlet height for both the main and side channels was set to 15 mm. The volume flow rates were measured at 25.7 ± 1.0 L/s for the main channel and 16.8 ± 0.1 L/s for the side channel. This corresponds to an average flow speed of 3.45 m/s in the main channel and 2.24 m/s in the side channel.

Considering the height of the opening for the characteristic length scale (15 mm), the calculated Froude numbers for the main and side channels were 9 and 6, respectively, while the calculated Reynolds number was 5.1 × 10⁴ for the main channel and 3.3 × 10⁴ for the side channel. The flow is shown from various perspectives in the figure below. For the validation of our CFD model, the measured free surface was then compared to the results obtained from the LIDAR measurements.

2.2. Computational Fluid Dynamics Governing Equations

In CFD simulations, the flow domain is divided into a mesh of computational cells or elements, and the behavior of the fluid is solved numerically within each cell. In the case of Reynolds-averaged Navier–Stokes equations (RANS), the solution variables for the velocity, pressure, and other scalar quantities of the exact Navier–Stokes equations are decomposed into mean and fluctuating components.

Instantaneous variables in the continuity equation and momentum equation are substituted with time-averaged expressions to derive the instantaneous time-averaged momentum equations, as represented by Equations (1) and (2) in tensor notation [47].

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0$$

(1)

$$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_l}{\partial x_l}\right)\right]+\frac{\partial }{\partial x_j}\left(-\rho \overline{u_i^{’}{u}_j^{’}}\right)$$

(2)

Here ρ is the density, x_i is the coordinate in the i direction, u_i is the velocity magnitude in the i direction, μ is the dynamic viscosity, and δ_ij is the Kronecker delta.

To solve Equation (3), Reynolds stress $-\rho \overline{u_i^{’}{u}_j^{’}}$ needs to be modeled. This is performed with the use of the Boussinesq hypothesis which relates Reynolds stresses to mean velocity gradients:

$$-\rho \overline{u_i^{’}{u}_j^{’}}={\mu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\left(\rho k+{\mu }_t\frac{\partial u_k}{\partial x_k}\right){\delta }_{ij}$$

(3)

Depending on the definition of turbulent viscosity μ_t, several turbulence closure models exist. These models incorporate additional transport equations for turbulence variables. For k-ω based turbulence models, turbulent viscosity is defined as a ratio of turbulent kinetic energy (k) times density (ρ) to specific dissipation rate ω:

$${\mu }_t={\alpha }^{*}\frac{\rho k}{\omega }$$

(4)

In the k-ω SST model [48], the transport behavior of Reynolds stresses is achieved by introducing a limiter to the formulation of eddy viscosity. Additionally, the k-ω SST-SAS model dynamically adapts to resolved structures by integrating the von Kármán length scale into the turbulence scale equation, resulting in behavior similar to the more computationally demanding Large Eddy Simulations (LES) as detailed in [49].

2.3. Multiphase Flow Models (VOF and Mixture)

The complex dynamics of multiphase flows demand the use of unique techniques to accurately predict the location of a free surface. Among these techniques, the VOF and mixture models are widely employed due to their versatility and reliability [50].

The VOF model tracks the shape and position of the interface of two or more immiscible fluids. On the other hand, the mixture model assumes fluid mixing at the interface, calculating the density and velocity of the resulting mixture. These models find applications in simulating various phenomena, including free surface flows, bubbly flows, and turbulent mixing [47].

In numerical simulations, the volume fraction (α) plays a pivotal role, representing the presence and proportion of each phase within a computational cell (0 to 1 range). For the VOF method, the interphase boundary is established by solving the continuity equation for the volume fraction:

$$\frac{1}{{\rho }_q}\left[\frac{\partial }{\partial t}\left({\alpha }_q{\rho }_q\right)+\nabla ·\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)\right]=\frac{1}{{\rho }_q}\left[{\displaystyle \sum }_{p=1}^n\left({\dot{m}}_{pq}-{\dot{m}}_{qp}\right)\right]$$

(5)

Here ${\dot{m}}_{qp}$ is the mass transfer from phase q to p and ${\dot{m}}_{pq}$ is the mass transfer from phase p to phase q. The volume fraction equation is not solved directly; instead, the volume fraction of the primary phase is calculated based on the following constraint:

$${\displaystyle \sum }_{q=1}^n{\alpha }_q=1$$

(6)

For the mixture model, we are solving the momentum, continuity, and energy equation for the mixture and algebraic expressions for the relative velocities. The volume fraction equation in this case is obtained from the continuity equation for secondary phase p:

$$\frac{\partial }{\partial t}\left({\alpha }_p{\rho }_p\right)+\nabla ·\left({\alpha }_p{\rho }_p{\overrightarrow{v}}_m\right)=-\nabla ·\left({\alpha }_p{\rho }_p{\overrightarrow{v}}_{dr,p}\right){\displaystyle \sum }_{q=1}^n\left({\dot{m}}_{qp}-{\dot{m}}_{pq}\right)$$

(7)

Solving the volume fraction equation increases the computational cost and time, making both methods very costly. Although this approach is very effective in capturing complex free surface phenomena, it can produce significant errors in the vicinity of the interface [51]. For instance, the VOF approach typically yields a location with an error approximately three times that of the cell size [32].

2.4. Numerical Grid, Boundary Conditions, and Wall Functions

For high Reynolds number flows, smaller cell sizes are necessary to effectively capture the intricate behavior of turbulent flows as complex flow structures associated with high Reynolds numbers occur at smaller length scales. By reducing the cell size, more grid points are available to capture the fluid flow characteristics, such as vortices, boundary layer behavior, and turbulence effects more accurately. However, in transient simulations, reducing the cell size imposes restrictions on the chosen time step of the simulation. The time step must be sufficiently small to satisfy the Courant–Friedrichs–Lewy (CFL) condition, which ensures numerical stability [52]. When dealing with high Reynolds number flows, a trade-off between cell size and time step arises. Achieving an optimal balance between cell size, time step, and computational resources is crucial for obtaining accurate and efficient simulation results.

To reduce the calculation time and cost of the simulation, the geometry of the channel was adjusted accordingly as can be seen in Figure 3. The inflow length for both the main and side channels was reduced to a length of 0.5 m, while the channel downstream of the junction was reduced to 1.5 m. Considering that the standing wave height in the investigated scenario was relatively low, the height of the entire domain was adjusted to 300 mm. The boundary conditions were specified as follows: a pressure outlet boundary condition (red) was imposed at the channel exit and the top of the channel, for channel walls a no-slip wall boundary condition (grey) was assumed, and a constant velocity inlet boundary condition (blue) was prescribed for both the main and side channels.

A numerical grid was created in Ansys ICEM CFD. The medium mesh shown in Figure 3 was a structured grid consisting of 3D blocks with cells of only hexahedral shapes. To accurately capture both inlet velocity profiles, the mesh was refined at the bottom of the channel. For the employment of wall functions, the size of the first wall cell in the entire domain was set to 3 mm. The employment of wall functions means the viscous sublayer is not resolved. This compromise is justified by the relatively insignificant role that the specifics of the near-bed flow dynamics play in these types of flows [53]. The y+ values were monitored throughout the numerical simulation and checked at post-processing. Near the upper boundary of the domain, the numerical mesh was coarser as air circulation in this region had relatively lower significance.

2.5. Solution Methods Settings

The performance of the k-ω SST and k-ω SST-SAS turbulent models was assessed. For both employed multiphase models, the implicit formulation was used. The pressure-implicit with splitting of operators (PISO) algorithm was used for the coupling of pressure and velocity. Regarding spatial discretization, the following schemes and algorithms were used: Least Squares cell-based method for gradients, Presto! algorithm for pressure, Second-Order Upwind scheme for momentum equations, Compressive scheme for volume fractions, Second-Order Upwind scheme for turbulent kinetic energy, and specific dissipation rate. For the transient formulation, the Second-Order Implicit scheme was chosen.

Table 1. Values of simulation parameters.

Parameter	Unit	Value
Main channel water velocity	m/s	3.45
Main channel opening	mm	15
Side channel water velocity	m/s	2.24
Side channel opening	mm	15
Main channel Reynolds number	-	51,516
Side channel Reynolds number	-	33,532
Main channel Froude number	-	9
Side channel Froude number	-	6
Operating pressure	Pa	101,325
Temperature	°C	15

includes an overview of the input and operating parameters alongside the properties for both water and air. The standard initialization method was chosen, where the initial values were computed from all the zones. The volume fraction of all the cells was set to zero initially, representing the domain filled only with air. The values for k and ω were determined as 0.03173 m²/s² and 1.77∙105 s⁻¹, respectively. For a more in-depth understanding of the determination of initial values, please refer to Ansys Fluent Theory Guide [47].

2.6. Convergence Criterion

The convergence criterion for all the simulations was established by setting the residuals below 10⁻⁴, which is widely accepted by the CFD community for this magnitude of the transient multiphase flow. Due to the turbulent nature of the investigated fluid flow, mass fluctuations were present in the domain. This could be easily solved by extending the channel further and allowing the flow to settle. However, such an extension would lead to an impractically excessive length and would render the simulation economically unviable. Therefore, to tackle this issue while still ensuring accurate and reliable results, we monitored mass fluctuations in the domain, calculated as the difference between the water mass inflow from both inlets and the water mass outflow. In addition, the following criterion was specified: the integral of mass imbalance on the investigated time-averaged interval of 10 s must be less than half of the maximum amplitude of the mass fluctuations. This criterion allows for a rigorous assessment of the convergence and ensures that the model accurately captures the dynamic behavior of the system while not introducing any mass accumulation in the domain.

2.7. HPC Cluster

The simulation was run at the Faculty of Mechanical Engineering, University of Ljubljana, on an HPC cluster with two AMD Epyc 7402 processors and 120 Gb Ram per node. For our simulation, three to five nodes were used for an individual case. To advance the simulation, a time step of 0.00025 s (medium mesh) was used with a maximum of 20 iterations per time step, with real-time data storage for direct simulation control. The simulation was run for a total of 20 s of flow time.

3. Results

3.1. Flow Development

At the onset, the entire fluid domain was filled only with air. The development of flow over the first 10 s is depicted in Figure 4. Here, the free surface with a volume fraction of 0.5 is colored-coded by the instantaneous maximum velocity magnitude. Initially, the flow front in both channels moves at the specified inlet velocity until the two flow structures collide at the junction at approximately 0.21 s. At approximately 1.0 s, the water structure gradually morphs into a standing wave (no hydraulic jump was observed in the upstream channel sections). At 2.5 s, the bulk liquid reaches the outlet of the domain, also discussed in Section 4.1. and shown on the graph of mass fluctuations.

Over the next few seconds, the flow structures begin to stabilize. The standing wave height gradually diminishes with the increasing distance from the junction. The last image in Figure 4 shows the free surface of a fully developed quasi-steady state junction flow at a flow time of 10 s. Beyond this point, another 10 s of flow time was simulated to investigate the quasi-steady-state junction by time-averaging the variables of interest. This timeframe has proven adequate in similar studies [54] as it effectively captures the multiple characteristic oscillation periods of the variables of interest. The time-averaged results were then compared against the LIDAR measurements.

3.2. Grid Independence Test

The grid independence test was conducted using four meshes with different levels of refinement. The number of computational cells for each mesh is showcased in

Table 2. Grid independence test for the k-ω SST VOF model.

Model					Point 1 (1000, 0, 30)		Point 2 (2000, 100, 20)
k-ω SST VOF	No. of Cells (-)	Max Cell Length (mm)	Time Step (ms)	Calculation Time (Core Hours)	Velocity Magnitude (m/s)	RMSE (m/s)	Velocity Magnitude (m/s)	RMSE (m/s)
Coarse mesh	1.42 × 105	35.3	0.505	1478	2.00	0.12	1.87	0.05
Medium mesh	2.130 × 106	7.35	0.250	14,224	2.53	0.08	1.85	0.16
Fine mesh	3.923 × 106	5.7	0.220	29,690	2.52	0.18	1.94	0.24
Ultra fine	8.004 × 106	4.7	0.183	89,368	2.55	0.20	1.92	0.20

. Since wall functions were employed, when refining the mesh, the distance between the first cell and the wall remained unchanged for all three grids. The y+ value in the initial cell at the wall was maintained at 49, respectively, while the maximum reported y+ value during the simulation was 263. The selected parameter of interest for the grid independence test was the mean velocity magnitude and its corresponding root mean square error at two designated points, point 1 (1500, 0, 30) and point 2 (2500, 100, 20). The point coordinates refer to the actual geometry of the physical model and not to the geometry of the computational mesh. To allow comparison, the Courant number needed to be the same for all three meshes which corresponded to different time step values, showcased in Table 1. Throughout most of the simulation, the Courant number consistently remained in the range of 0.3 to 0.5. Notably, in specific intervals, it approached values higher than 0.9 but remained below the critical threshold of 1.

The results of the grid independence test are presented in Table 2. Upon examination, it is evident that at point 1, the variance in velocity between the medium and fine grid is less than 0.4%. However, at point 2, further downstream where computational errors tend to accumulate, the error increases to 4.7%.

4. Discussion

4.1. Result Credibility Assessment

The continuity equation residuals reached the values of 10⁻⁴, whereas all the other residuals were lower, typically on the order of 10⁻⁶ or less. The credibility of the results was verified through the inspection of mass fluctuations which are present in the channel due to the periodicity of the flow. In Figure 5, mass fluctuations which are calculated as the cumulative mass flow from both inlets minus the outlet are presented. In Figure 5a, the fluctuations for all three investigated models, calculated on a medium-size mesh, were less than +/−4 kg/s, equivalent to less than +/−9.3% of the nominal water flow rate. Further, the integral of mass fluctuations on the 10 s interval was less than half of the maximum amplitude, suggesting that there was no significant mass accumulation occurring in the channel that would undermine the credibility of the simulation results. In Figure 5b, the mass fluctuations for different grid sizes are presented for the k-ω SST VOF model. The most pronounced difference is in the value of the fluctuation period between the coarse and medium mesh. Importantly, additional refinement beyond the medium mesh does not significantly impact the mass fluctuations in the channel.

4.2. CFD and LIDAR Comparison

Below, we present the free surface profiles of the junction using both the (CFD) k-ω SST-VOF model and LIDAR measurements. These profiles are showcased at three distinct cross sections along the x-axis, specifically at distances of 1500 mm, 1700 mm, and 2200 mm from the main channel inlet (with distances referencing the physical measuring line).

The CFD-derived free surface profiles, featuring a volume fraction of 0.5, are illustrated with solid black lines, while the color-coded profiles represent the root mean square deviation (RMSD) of volume fraction, α_RMSD, providing insight into the mean fluctuation of volume fraction throughout the time averaging process. This presentation effectively showcases the dynamic and unsteady nature of the free surface. The LIDAR-measured free surface profiles are represented with dotted red lines.

It is important to note that the data obtained from CFD were averaged on a time interval of 10 s, while the LIDAR results were averaged on a shorter 2 s interval. Despite this difference, the CFD results proved to be in good agreement with the LIDAR measurements. The selected model allowed for the identification of several entrainment zones and air pockets which were visually observed but could not be measured by LIDAR.

As demonstrated by Figure 6, there is good agreement between the CFD-calculated and LIDAR-measured free surface profiles near the channel center where the maximum depth of about 200 mm was attained, and in its rightmost 100 mm. The CFD results indicate the presence of air pockets or intense air entrainment zones from z = 70 mm to 240 mm and from z = 280 mm to 360 mm. However, due to the inherent limitations of the LIDAR method, such entrained air structures could not be experimentally measured, thus preventing us from verifying their size and position using this measurement technique. The α_RMSD values for the volume fraction are higher on the uppermost left side, where the flow experiences greater turbulence.

At x = 1700 mm, Figure 7, the comparison of the CFD and LIDAR free surface profiles once again shows a better agreement in the central and right sections of the channel (y > 150 mm). The highest measured water level is observed on the leftmost side which holds true for both the LIDAR and CFD data. However, the LIDAR measurements were about 20 mm higher which can be attributed to the refraction of light near the wall. As we move away from the wall towards the center of the channel, the disparities between the two profiles become more pronounced, primarily due to the scattering of water droplets in the vicinity. At a distance of 225 mm from the left wall, the data from the LIDAR and CFD calculations closely align at a height of 140 mm, respectively. The agreement on the far-right side is notably good. In comparison to the profile at x = 1500 mm, the α_RMSD values are higher, indicating greater fluctuations in this region.

At x = 2200 mm, as shown in Figure 8, where the flow is less turbulent, the CFD profile coincides with the LIDAR profile in a relatively wide central segment and remains quite close to it in other areas. The water surface height varies from 50 to 95 mm across the entire cross section. Notably, the α_RMSD values remain high.

In Figure 9, the lateral flow profile along the channel centerline with k-ω SST VOF on medium mesh is presented. The flow structures exhibit pronounced three-dimensional characteristics. Note that since the LIDAR free surface profiles were measured transversally, the longitudinal measurement data resolution is somewhat sparse, with 10 cm between the neighboring measurement points (i.e., equal to the distance between the transversal LIDAR profiles). Nevertheless, the CFD longitudinal profile is captured relatively well with the LIDAR data.

Additionally, the calculation for k-ω SST mixture and k-ω SST-SAS mixture on a medium-sized mesh was run. A comparison of the free surface profiles at x of 1500, 1700, and 2200 mm is showcased in Figure 10. Near the junction, all three models predict the position of a free surface relatively consistently. However, the further we move downstream, particularly at x = 2200 mm, both mixture models produced an excessively perturbed free surface, which we argue is due to the wide distribution of the particulate phase. Furthermore, the employment of the computationally more intensive SAS model gave no additional benefit. Since SAS was about 25% more computationally expensive than the SST, there is no advantage of the SAS over the SST model in the present flow setup. The position of the obtained free surface profile was significantly affected by the choice of the employed multiphase model. Depending on whether the VOF or mixture two-phase model was used, there were notable variations in the position of the time-averaged free surface.

Upon closer examination of the RMSD values for the selected three models, it becomes evident that the k-ω SST VOF model exhibits higher RMSD values. This indicates a more pronounced periodicity, particularly in the near-wall region where the turbulent mixing of water and air occurs due to wave impact with the wall, in comparison to the corresponding mixture models. This observed periodicity aligns with the experimental observations.

Examining the cross sections at x = 1500 and 1700 mm reveals significant agreement among all three CFD models and the LIDAR data, particularly when factoring in the inherent uncertainty associated with the experimental method and the complexity of the phenomenon at hand. However, at x = 2200 mm, the agreement diminishes for both mixture models employed. Both models project a more turbulent free surface, whereas the VOF model portrays a more uniform free surface, aligning more closely with the LIDAR data. By further inspection of the RMSD values for the selected three models, it becomes evident that the k-ω SST VOF model exhibits higher RMSD values, indicating a more pronounced periodicity, specifically in the near-wall region, where wave impact happens, compared to the corresponding mixture models. This periodicity is supported by the experimental observations. Assessing the improved consistency with the experimental results downstream from the standing wave, it can be inferred that the k-ω SST VOF model stands out as a superior choice in the modeling approach. The accuracy of the CFD results can reduce the need for physical model testing which promises to significantly streamline and expedite the design process of hydraulic structures featuring open channel junctions.

5. Conclusions

This study presents a CFD model of a supercritical open channel junction flow on a T-shaped junction. The selected flow case is characterized by high Reynolds (3.4 × 10⁴ and 5.2 × 10⁴, respectively) and Froude numbers (6 and 9, respectively), the formation of a standing wave at the junction, and the absence of upstream hydraulic jumps. The CFD model is able to estimate the water depth in channel junctions where 3D standing waves of significant heights can form. Quality CFD models can potentially eliminate the need for constructing physical models, thus streamlining and accelerating the design of hydraulic structures with open channel junctions. To our knowledge, no prior CFD studies have successfully undertaken such an approach for supercritical junction flows.

The comparison of the CFD results with the reference LIDAR measurements reveals a good agreement between the predicted and measured free surface elevations, especially when the shortcomings of the LIDAR method in certain measurement zones and situations are considered. The CFD models with appropriate turbulence and multiphase flow models were capable of accurately predicting the free surface fluctuation and free surface elevation in the supercritical junction. The results suggest that a good convergence can be achieved in a given domain geometry employing a simple structured hex mesh with near-wall refinement, although at a high computational expense due to the large number of computational cells and small transient timestep. Comparing the multiphase flow models, both VOF and mixture were able to detect deep air pockets and other entrainment zones that were visually observed but could not be measured by LIDAR due to being below the uppermost free surface. The VOF model proved itself more accurate than the mixture model, particularly in the regions downstream of the junction where the flow fluctuations partly subsided. Comparing the turbulence models, the k-ω SST model was demonstrated to perform well under such simulation conditions and no advantage was gained by using the more computationally intense SAS model. Overall, this study provides valuable insights into the performance and limitations of the employed turbulence and multiphase models for simulating supercritical junction flows.