Two-phase flow modelling by an error-corrected population balance model

Abstract

High-velocity aerated flow is a common phenomenon in spillways. Its accurate modelling is challenging, mainly due to the lack of realistic physics in the conventional two-phase models. To this end, this study establishes a population balance model (PBM) approach to account for the evolutionary process of air bubbles. The air-water flow in a stepped chute is examined. The model performance is evaluated by statistical metrics: correlation coefficient (CC), root mean squared error (RMSE), and mean absolute error (MAE). Compared with conventional models, the PBM generates improved air-water predictions. However, the flow parameters are still underestimated, particularly in areas with intense air-water interactions. For further development, an error-corrected PBM (EPBM) is proposed by incorporating machine learning (ML) techniques into the PBM. Compared with the PBM, the EPBM leads to a mean augmentation in velocity prediction by 19.8% for the CC, 73.0% for the RMSE, and 77.1% for the MAE. The gains in air concentration estimation are 2.0%, 67.6% and 73.5%, respectively. The EPBM generates the most accurate results, with 99.6% and 89.6% of the velocity and air concentration predictions within a 20% relative error range. The main contributions are establishing a PBM for air-water flows and developing an error-corrected PBM using ML.

KEYWORDS:

• Spillway
• aeration
• two-phase model
• population balance model
• machine learning

1. Introduction

Safe operation of flood discharge structures depends on not only water flows, but also the simultaneous movement of air flows in the system. Due to the significant difference in density, water and air are often separated by a sharp interface. However, if the turbulence level is sufficiently high to overcome both surface tension and gravity effects, intensive air-water mixing (aeration) occurs (Chanson, 1993). Aeration is a highly turbulent process commonly appearing in head-head spillways (Figure 1) and low-level outlets. With the presence of air, the hydraulic performance of the structures is significantly affected by spatial changes in bulk properties, flow depth, turbulence characteristics, etc. (Wood, 1991).

Figure 1. Aeration phenomenon (white water) in a prototype spillway during flood discharge (image by James Yang).

Spillway flows often feature high velocity and strong turbulence, giving rise to surface disturbance and air entrainment. This process transforms the flow into a multiphase flow composed of a mixture of water droplets and waves with dispersed air bubbles (Valero, 2018). Due to its complexity, efforts have been dedicated to understanding the air-water flow behaviours. Laboratory experiments are a dominant approach. The earliest studies date back to Ehrenberger (1926) who physically models air-water flows in a chute. Lane (1939) first defines the onset condition for self-aeration in open-channel flows, i.e. aeration starts at the point the boundary layer reaches the free surface. Later, Bauer (1951) and Halbronn (1952) propose expressions for predicting boundary layer development. Alternative entrainment theories are also developed. Hickox (1945) suggests that the kinetic energy of surface eddies should be strong enough to overcome the surface tension. Soo (1956) argues that the turbulent velocity should be greater than the bubble rise velocity. More laboratory investigations are made by Straub and Anderson (1958), Killen (1968) and Gulliver et al. (1990).

With new measurement and instrumentation technologies (e.g. non-intrusive techniques), the past decades have witnessed progress in air-water flows research. In a spillway model with a fibre-optical probe, Boes and Hager (2003) perform experiments, from which empirical models are introduced to estimate air concentration and flow velocity. By varying the slope of a spillway chute, Ohtsu et al. (2004) investigate the aerated flows, presenting a classification of flow regimes and calculation of flow depth. Felder and Chanson (2009) conduct turbulence measurements in a stepped spillway, in which the validity of the Froude and Reynolds similarities are tested. Felder and Chanson (2016) focus on a stepped spillway placed on an embankment, aiming to provide simple design criteria in terms of dimensionless residual energy. Zhang and Chanson (2017) obtain experimental data from a stepped chute and propose an analytical solution for the estimation of air diffusion. Kramer and Chanson (2018) study the transitional flows in a laboratory spillway, presenting a detailed description of air-water flow properties and an image-based analysis of pool depth fluctuations. Kramer and Chanson (2019) employ a novel filtering technique for optical flow measurements, demonstrating the potential of non-intrusive image-based methods for estimating air-water flow velocity with high spatial and temporal resolution. Kramer et al. (2019) apply the LIDAR technique to spillway flows, concluding that it is applicable in fully aerated flow regions, but not in clear water and rapidly varied flows.

Computational fluid dynamics (CFD) is another approach to study spillway flows. Chen et al. (2002) conduct one of the earliest numerical simulations of stepped spillway overflows, with flow pressure and depth compared with experimental results, showing good agreement. Cheng et al. (2006) evaluate the performance of turbulence models in reproducing spillway hydraulics. Kositgittiwong et al. (2013) compute the velocity profiles of spillway flows numerically, concluding that CFD is a reliable method. Bai et al. (2017) report the two-phase flows in a unique V-shaped stepped spillway, exhibiting the complex flow structures. Liu et al. (2018) simulate the flow in a shaft spillway, analysing its hydraulic features, e.g. flow pattern and air core distribution. Imanian and Mohammadian (2019) investigate the flow properties of a high-head ogee-crested spillway, presenting results of turbulence characteristics, discharge coefficient, and pressure field. Ghaderi et al. (2020) introduce a spillway with trapezoidal labyrinth steps. Numerical results show that this step layout heightens the friction and augments energy losses.

CFD simulations often reproduce mean two-phase flow characteristics, e.g. flow pattern, velocity, pressure, etc. (Cheng et al., 2022; Dou et al., 2021; El Assad et al., 2021; He et al., 2022; Kyriakopoulos et al., 2022). However, estimating a mixture flow remains a challenging task. Valero and Bung (2015) model the self-aeration and air transport processes in spillway flows, witnessing remarkable differences in the void fraction compared to experiments. Wan et al. (2017) present a smoothed particle hydrodynamics (SPH) method to predict the flow with reaeration in a stepped spillway. The simulated results of mixture flow velocity exhibit notable deviation from the experimental ones. A similar conclusion is drawn in the studies by Morovati et al. (2016). Downstream of a spillway aerator, Teng et al. (2016) find that numerical models are incapable of making accurate predictions of air concentration. With a modified two-fluid model, Yang et al. (2019) examine the high-velocity air-water flows in a large chute, concluding that the near-bottom air concentration is overestimated and air detrainment is underestimated. In the work by Güven and Mahmood (2021), one also notices the insufficiency of CFD models in satisfactory estimation of air concentration in chute spillways.

The discrepancy is presumably attributable to the insufficient physical representations of the two-phase flow physics (Yang et al., 2019). In the mixture flow models, air bubbles are assumed to be spherical in shape and remain constant in size, regardless of the collision, aggregation and breakage processes (Lei et al., 2021), which is not the case in reality (Sarhan et al., 2016, 2018a, 2018b). To rectify these issues, the population balance model (PBM) provides an alternative. It is an Eulerian-Eulerian approach capable of simulating the evolutionary process of mixture particles. The model has primarily been used for gas-liquid flows (e.g. in fermenters) in chemical and other industrial applications (Castellano et al., 2018; Raesi & Maddahian, 2022). Its use in high-speed aerated flows is limited. A hydraulic jump study by Xiang and Tu (2016) shows some promising results. However, errors of void fractions still exist in the free-surface regions, necessitating improvements.

Motivated by the arguments mentioned above, the objective of the present study is twofold: (a) to assess the performance of the PBM in predicting high-speed two-phase flows, and (b) to explore the potential of soft computing approaches in improving the accuracy of the PBM. Due to high efficiency and adaptability, soft computing techniques (e.g. machine learning) are receiving a growing interest in spillway hydraulics, with successful examples in estimating energy loss and flow discharge capacity (Roushangar et al., 2014; Roushangar et al., 2019). The main contribution of the work is (a) to implement a novel CFD algorithm for a better physical description of air-water flows and (b) to develop a machine learning assisted CFD approach for improved prediction of two-phase flow properties. The study aims to establish an accurate model for modelling two-phase air-water flows in chute spillways, which also provides a solution for flow modelling in similar flow environments.

2. Modelling methodology

The study consists mainly of two parts: evaluations of the two-phase models and the PBM, and development of a hybrid approach for accurate two-phase flow prediction. The emphasis is placed on the PBM and its improvement with error correction.

2.1. Turbulence model

A review of previous works shows that the k-ϵ turbulence model is a reliable method for simulating spillway flows (Chen et al., 2002; Kositgittiwong et al., 2013; Teng et al., 2016; Valero & Bung, 2015; Yang et al., 2019). It has several variants, including the standard, the renormalisation group (RNG), and the realisable k-ϵ model. The RNG model is chosen for this study based on its satisfactory accuracy (Bombardelli et al., 2011; Morovati et al., 2016; Valero & García-Bartual, 2016). The turbulence kinetic energy, k, and its rate of dissipation, ϵ, are obtained from the following equations (1) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho k\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho k{\overrightarrow{v}}_i\right)& =\frac{\mathrm{\partial }}{\mathrm{\partial }xj}\left({\alpha }_k{\mu }_{\text{eff}}\frac{\mathrm{\partial }k}{\mathrm{\partial }xi}\right)\\ & +G_k+G_b-\rho ϵ-Y_M+S_k\end{array}$(1) (2) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho ϵ\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho ϵ{\overrightarrow{v}}_i\right)& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left({\alpha }_{ϵ}{\mu }_{\text{eff}}\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }{x}_j}\right)\\ & +C_{1ϵ}\frac{ϵ}{k}\left(G_k+C_{3ϵ}{G}_b\right)\\ & -C_{2ϵ}\rho \frac{{ϵ}^2}{k}-R_{ϵ}+S_{ϵ}\end{array}$(2) where G_b (G_k) = turbulence kinetic energy due to buoyancy (mean velocity gradients), Y_M = contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, S_k, S_ϵ = user-defined source terms, μ_eff = effective viscosity, σ_k, σ_ϵ = turbulent Prandtl numbers, and C_1ϵ, C_2ϵ, C_3ϵ = constants.

The drag force dominates the movement of air bubbles in the dispersed flow. It acts on bubbles as a resisting force caused by the non-uniform pressure distribution of the surrounding water. The drag coefficient is determined based on the theory by Schlichting and Kestin (2017). One is referred to Sarhan et al. (2016, 2018a, 2018b) for detailed model descriptions.

2.2. Two-phase flow models

In combination with the PBM, three conventional two-phase models are used in the study: the Volume of Fluid (VOF), Two-Fluid and Mixture models.

2.2.1. VOF model

The VOF model is a simple and easy-to-implement method widely used in industrial applications such as two-phase flow in pipelines, droplet dynamics, and multiphase flow in chemical reactors. The VOF model is particularly useful for simulating problems involving the movement and mixing of immiscible fluids with a free surface. In this model, the fluids share a single set of momentum equations, and each phase is represented by its VOF in a computational cell. Water is usually treated as a primary phase and air as a secondary one. In a cell, the water and air have the same flow velocity ($\overrightarrow{v}$), flow pressure (p) and other turbulence properties (Bai et al., 2017; Güven & Mahmood, 2021). The momentum equations are expressed by (3) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho \overrightarrow{v}\right)+\mathrm{\nabla }\cdot \left(\rho \overrightarrow{v}\overrightarrow{v}\right)\\ & =-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left[\mu \left(\mathrm{\nabla }\overrightarrow{v}+\mathrm{\nabla }{\overrightarrow{v}}^T\right)\right]+\rho \overrightarrow{g}+{\overrightarrow{F}}_1\end{array}$(3) where t = time, g = gravity, ρ = density, μ = dynamic viscosity, and ${\overrightarrow{F}}_1$ = interfacial force. In each cell, the sum of the volume fraction α_q for the liquids is equal to one: ${\sum }_{q=1}^n{\alpha }_q=1$ (n is the number of phases). α_q is calculated by (4) $\begin{array}{r}\frac{\mathrm{\partial }\left({\alpha }_q\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\overrightarrow{v}}_q{\alpha }_q\right)=0\end{array}$(4)

2.2.2. Mixture model

The Mixture model is a simplified two-phase model. The method simulates the phases (fluid or particulate) by solving the momentum, continuity and energy equations for the mixture, the volume fraction equations for the secondary phases, and algebraic expressions for the relative velocities. It predicts two-phase flows where the phases move at different velocities but assume local equilibrium over short spatial length scales. Alternatively, it simulates homogeneous two-phase flows with strong coupling and the phases moving at the same velocity. The model can also be used to calculate non-Newtonian viscosity. Typical applications include sedimentation, cyclone separators, particle-laden flows with low loading, and bubbly flows where the gas volume fraction remains low. The Mixture model differs from the VOF model in three respects: it allows the phases to interpenetrate, to move at different velocities, using the concept of slip velocities and considering inter-phase interaction of mass, momentum and energy transfer. (Viitanen et al., 2020). The continuity equation is given by (5) $\begin{array}{r}\frac{\mathrm{\partial }\left({\rho }_m\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_m{\overrightarrow{v}}_m\right)=0\end{array}$(5) where ρ_m = mixture density and ${\overrightarrow{v}}_m$ = mass-averaged velocity, calculated by (6) $\begin{array}{r}{\overrightarrow{v}}_m=\sum _{q=1}^n\frac{{\alpha }_q{\rho }_q{\overrightarrow{v}}_q}{{\rho }_m}\end{array}$(6) The momentum equation for the mixture is obtained by summing the individual ones for all phases. (7) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\rho }_m{\overrightarrow{v}}_m\right)+\mathrm{\nabla }\cdot \left(\rho {\overrightarrow{v}}_m{\overrightarrow{v}}_m\right)\\ & =-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left[{\mu }_m\left(\mathrm{\nabla }{\overrightarrow{v}}_m+\mathrm{\nabla }{\overrightarrow{v}}_m^T\right)\right]+{\rho }_m\overrightarrow{g}+{\overrightarrow{F}}_2\\ & +\mathrm{\nabla }\cdot \left(\sum _{q=1}^n\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_{dr,q}{\overrightarrow{v}}_{dr,q}\right)\right)\end{array}$(7) where μ_m = dynamic viscosity of the mixture, ${\overrightarrow{F}}_2$ = body force, and ${\overrightarrow{v}}_{dr,q}$ = drift velocity for secondary phase q.

2.2.3. Eulerian model

Different from the VOF concept, the Eulerian approach, often called the Two-Fluid model, is a two-phase two-fluid model. It treats the two immiscible yet penetrating phases separately in a cell – they do not share the same velocity. However, the pressure is the same for the phases. The model framework is based on ensemble-averaged mass and momentum transport equations for each of the two phases. Air bubbles are treated as spherical particles; an average bubble size needs to be specified by the user. The model includes lift force, turbulence dispersion force, virtual mass force, drag force and wall-lubrication force. The continuity equation is expressed by (8) $\begin{array}{r}\frac{\mathrm{\partial }\left({\alpha }_q{\rho }_q\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)=0\end{array}$(8) The momentum equation is given as (9) $\begin{array}{rl}& \frac{\mathrm{\partial }\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q{\overrightarrow{v}}_q\right)\\ & =-{\alpha }_q{\mathrm{\nabla }}_p+\mathrm{\nabla }\cdot {\overline{\overline{\tau }}}_q+{\alpha }_q{\rho }_q\overrightarrow{g}\\ & +\left({\overrightarrow{F}}_q+{\overrightarrow{F}}_{\text{lift},q}+{\overrightarrow{F}}_{\text{wl},q}+{\overrightarrow{F}}_{\text{vm},q}+{\overrightarrow{F}}_{\text{td},q}\right)\end{array}$(9) where ρ_q = density of phase q, ${\overrightarrow{v}}_q$ = velocity of phase q, and ${\overline{\overline{\tau }}}_q$ = stress–strain tensor of q^th phase, ${\overrightarrow{F}}_q$ = external body force, ${\overrightarrow{F}}_{\text{lift},q}$ =  lift force, ${\overrightarrow{F}}_{\text{wl},q}$ = wall lubrication force, ${\overrightarrow{F}}_{\text{vw},q}$ = virtual mass force, ${\overrightarrow{F}}_{\text{td},q}$ = turbulent dispersion force. A detailed model description is found in Teng et al. (2016) and Yang et al. (2019).

2.3. Population balance model

The PBM is a statistical formulation to describe the evolution of a population of particles in a multiphase flow (Wang et al., 2005). A spillway flow always involves an air phase. The size distribution of air bubbles and droplets evolves in conjunction with transport in the two-phase system. The evolutionary processes can be a combination of different phenomena like dissolution, growth, dispersion, aggregation, breakage, etc. In addition to momentum, mass, and energy balances, a balance equation is required to describe the particle population changes. This balance is generally referred to as the population balance, taking into account bubble breakup and coalescence. The PBM is based on the mechanisms of breakup of dispersed bubbles caused by shear forces and particle coalescence resulting from flow turbulence. The air-water interaction is controlled by the critical Weber and Capillary numbers. Besides general flow features, such a model leads to the identification of bubble-size distribution in different flow regimes and adds more realism to the modelling. The PBM approach is an advanced tool for modelling two-phase flow issues. The available literature shows that the approach is mainly used in nuclear, thermal, chemical and water treatment fields. The use in the hydraulic area is somewhat limited. The difference from the other applications is the much higher flow velocity in spillway flows.

A balance equation is employed to describe the changes in the particle population, as follows (10) $\begin{array}{r}\frac{\mathrm{\partial }\left(n_i\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(n_i{u}_i\right)=B_{B,i}-D_{B,i}+B_{C,i}-D_{C,i}\end{array}$(10) where n_i = bubble number density, u_i = bubble velocity, B_B,I (B_C,i) = source term of the bubble i generated by the breakage of all bubbles larger (smaller) than the diameter d_i, and D_B,I (D_C,I) = vanishing term caused by bubble disappearance (aggregation). Bannari et al. (2008) present a method for computing these parameters. The bubble aggregation rate Ω_ag is defined by (11) $\begin{array}{r}{\mathrm{\Omega }}_{ag}\left(V_i,V_j\right)=w_{ag}\left(V_i,V_j\right)P_{ag}\left(V_i,V_j\right)\end{array}$(11) where ω_ag = collision frequency, P_ag = aggregation efficiency, and V_i, V_j = bubble volume. Using the energy theory, Luo and Svendsen (1996) provide a procedure for calculating the breakage probability and collision frequency. A detailed description is given by Lei et al. (2021) and Raesi and Maddahian (2022).

2.4. Error-corrected CFD procedure

CFD is useful in the prediction of two-phase spillway flows. However, relevant studies report significant discrepancies in the predicted mixture properties (Güven & Mahmood, 2021; Morovati et al., 2016; Valero & Bung, 2015). Researchers have made efforts to improve simulation results, with machine learning (ML) as a promising approach (Hanna et al., 2020; Vinuesa & Brunton, 2022; Zhao et al., 2020), which is due to its robust learning and prediction performance.

Motivated by the potential of ML, this study is intended to develop an improved hybrid approach by coupling CFD with the ML technique to achieve enhanced results. Figure 2 illustrates the established CFD procedure that is based on error corrections. The workflow of the framework consists of three phases: CFD modelling, residual mapping, and simulation reconstruction.

Figure 2. Illustration of the error-corrected CFD approach.

Phase 1: CFD modelling. For the simulation of the bubbly flows, the Eulerian and Mixture models are first examined and coupled with the PBM, leading to the PBM-Mixture and PBM-Eulerian models. The VOF is used separately as it cannot be combined with the PBM. The numerical results (five models) are then compared with experimental data, and the residuals are determined. The best-performed CFD models are thus selected for further improvement.

Phase 2: residual mapping. The CFD simulation errors are estimated by Bayesian optimised ML models, namely support vector machine (SVM) and Bagging. They represent two distinct categories of algorithms: conventional ML and ensemble learning. The choice is based on their satisfactory performance in similar hydraulic studies. One is referred to Li et al. (2021) and Ahmadianfar et al. (2021) for detailed model descriptions and their performance.

In ML models, hyperparameter optimisation is the key to obtain accurate results. The conventional trial and error method is time-consuming and cannot guarantee global optimum. To this end, the Bayesian optimisation (BO) is integrated into the ML models to search for the optimal architecture by tuning hyperparameters. Based on Bayes’ conditional probability rule, the BO evaluates the results from the previous iterations and chooses values for the next iteration. Its primary principles are: (a) to establish a probabilistic surrogate model based on the Gaussian process, and (b) to apply an acquisition function to determine the next observation location, where the observation property is expected to be the optimum. The choice of the BO rests on its fast computation and successful applications in similar problems (Wu et al., 2019).

Phase 3: simulation reconstruction. The CFD results and ML-predicted residuals are combined to construct error-corrected models. Their performance is evaluated statistically, and the most accurate ones are generated for practical application.

2.5. Evaluation metrics

Comparing the experimental and numerical results involves a large amount of data. As a result, Statistical metrics are introduced to quantify model skills. The indices often refer to correlation coefficient (CC), root mean squared error (RMSE), and mean absolute error (MAE). The CC is expressed by (12) $\begin{array}{r}\text{CC}=\frac{{\sum }_{i=1}^N\left(O_i-\overline{O}\right)\left(S_i-\overline{S}\right)}{\sqrt{{\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2\sum _{i=1}^N{\left(S_i-\overline{S}\right)}^2}}\end{array}$(12) where O_i = i^th observation, S_i = i^th simulation,$\overline{O}$ = mean of observations, $\overline{S}$ = mean of simulations, and N = number of data. The CC ranges from −1 to 1. CC = 1 and −1 indicate a perfect positive and negative correlation, with CC = 0 implying zero correlation. The RMSE measures the discrepancy between modelled and observed values, defined by (13) $\begin{array}{r}\text{RMSE}=\sqrt{\frac{{\sum }_{i=1}^N{\left(O_i-S_i\right)}^2}{N}}\end{array}$(13) The MAE computes the mean of all individual errors, expressed as (14) $\begin{array}{r}\text{MAE}=\frac{1}{N}\sum _{i=1}^N\left(|O_i-S_i|\right)\end{array}$(14) Both the RMSE and the MAE vary from 0 to +∞. The closer to zero, the higher accuracy.

3. Experimental set-up

To validate the developed numerical models, the study adopts the experimental results from a physical model of a stepped spillway (Felder & Chanson, 2009). As shown in Figure 3, the test rig comprises a large feeding tank, an upstream broad-crested weir, a stepped spillway chute, and a tailwater channel. The model is constructed with Perspex sidewalls and smooth marine plywood staircase invert. The chute width is constant, 1.0 m. The weir is 0.6 m long, including the upstream rounded corner (radius 0.057 m). The chute steps are identical, totalling 20 in number, with step height H = 0.05 m and length L = 0.125 m. A point gauge marks the clear-water depth in the chute, and a calibrated broad-crested weir records the water discharge Q. A double-tip conductivity probe (scan rate: 20 kHz per probe sensor, and sampling duration: 45 s) measures the two-phase flow properties of air-bubble sizes and air concentration.

Figure 3. Sketch of the stepped spillway model with an overflow weir followed by 20 steps.

A pump circulates the water in the test rig, with Q = 0.0519 m³/s. The critical flow depth over the crest is d_c = 0.066 m, with the corresponding Reynolds number R = ρ_wVD_h/μ = 2.1 × 10⁵ (ρ_w = water density, V = flow velocity, D_h = hydraulic diameter). Under this condition, the flow is in skimming regime and free-surface aeration starts from step edge 7. The data of V and air concentration C, totalling 378 groups, are collected in the aerated region. A coordinate system is defined: the x-axis runs through the pseudo-bottom formed by step edges, and the y-axis is perpendicular to the pseudo-bottom. The origin is at the aeration onset point (step edge 7 in this study).

4. Results and analysis

4.1. CFD simulations

Together with the conventional two-phase models, the contribution from the PBM is evaluated in predicting air-water flow properties. The most representative parameters are mixture velocity and air concentration, and they are chosen for model appraisals and comparisons.

4.1.1. Mesh and boundary conditions

Figure 4 shows the numerical model of the stepped chute. The computational domain, containing the reservoir and the tailwater channel, is meshed with structured grids. The total number of cells is 0.65 million, with a min. and max. grid size of 5 and 10 mm (in the chute), respectively. For boundary conditions, the vertical upstream face is divided into two portions delimited by the water surface. The upper portion is set as pressure inlet of air, and the lower is one velocity inlet of water. The top of the domain is defined as pressure inlet of air and the exit in the tailwater pressure outlet. The spillway flow is deemed symmetrical; only half of the chute is simulated, with the central plane treated as symmetry. The remaining boundaries are assigned to walls.

Figure 4. Numerical model with boundary conditions and local mesh.

In the Mixture and Eulerian model, it is required to define a mean diameter d_a as the characteristic value for air bubbles. Only a limited number of studies have reported the d_a. In the study of spillway aerator flows, Yang et al. (2019) examine the effects of the d_a within the range 0.5–4.0 mm, suggesting that d_a = 0.5 mm leads to better estimations for C. Consequently, this value is used in the simulations. The PBM based model allows one to define the distribution of bubble diameters. A physical test indicates that 85% of the air bubbles fall within d_a = 0.5–4 mm (Yang et al., 2019). Based on this, the PBM parameters are defined as: min. d_a = 0.1, particle volume coefficient = π/6 (default value), number of particle size bins = 10, and ratio exponent = 2, which leads to d_a = 0.1–0.64 mm.

4.1.2. Grid independence

The grid quality needs to be checked to guarantee grid independent solutions. Three grid resolutions are examined, namely, 0.48 million (coarse), 0.65 million (medium) and 1.02 million (fine). All the grids are structured. The VOF method is used for the evaluation of mesh sensitivity. Along the chute, the flow pressure is compared; Figure 5 shows the results at step no. 10. The results from the three grids differ insignificantly, with a max. difference below 10%. Considering the computational efficiency, the medium-sized mesh is adopted.

Figure 5. Effects of grid density on simulation results: pressures acting on step no. 10. (a) along the horizontal face and (b) along the vertical face. p_m = max. pressure on the step, and l and h = horizontal and vertical distance from the step corner vertex.

4.1.3. CFD-predicted velocity

The air-water interfacial velocity is an emblematic parameter in mixture flows, with its distribution exhibiting typical flow characteristics. At each step edge, a power-law model approximates the mixture velocity distribution: (15) $\begin{array}{r}\begin{array}{ll}\frac{V}{V_{90}}={\left(\frac{y}{y_{90}}\right)}^{1/\text{β}},& \text{0}\le \frac{y}{y_{90}}\le 1\\ \frac{V}{V_{90}}=1,& \text{1}\le \frac{y}{y_{90}}\le 2.5\end{array}\}\end{array}$(15) where y₉₀ = y at C = 0.9, V₉₀ = characteristic air-water velocity at y = y₉₀, and β = constant. Felder and Chanson (2009) indicate that β = 10 leads to the best fit for the present experimental datasets.

Figure 6 compares the experimental, CFD and empirical velocity distributions at typical step edges. Seemingly, all the models capture the overall trend of the velocity. However, some models exhibit evident discrepancies. For example, at step 7, the VOF shows notable overestimations, demonstrating its inability to reproduce intense air-water interactions in the developing region. At step 14, the Mixture model underestimates the velocity near the free-surface area. The empirical correlation reasonably approximates the experimental results. However, its use might be restricted in practical applications. For instance, y₉₀ should be known to predict the velocity, which depends in return on the availability of C. The C measurement is a time-consuming task. It would be advisable to directly measure V instead of estimating it from the C data.

Figure 6. Comparisons of velocity among the experimental, CFD and empirical approaches at typical steps.

For a coherent evaluation of the models, Table 1 summarises their performance statistics for velocity estimation. Compared with the PBM based approach, the conventional models are less accurate in estimating interfacial velocities, with their CC = 0.44–0.50. The PBM module generates improved predictions, with the PBM-Eulerian leading to the highest accuracy (CC = 0.56). Compared with the conventional models, the PBMs improve, on average, by 18.9% in the CC, 3.8% in the RMSE, and 8.1% in the MAE. The enhancement is likely attributed to better physical descriptions of the system, e.g. bubble evolution.

Table 1. Statistical performance of the CFD models for interfacial velocity estimation.

Model	CC	RMSE	MAE	Model	CC	RMSE	MAE
VOF	0.444	0.145	0.095	PBM-Mixture	0.510	0.160	0.089
Mixture	0.403	0.144	0.093	PBM-Eulerian	0.560	0.140	0.088
Eulerian	0.503	0.179	0.101	–	–	–	–

4.1.4. CFD-predicted air concentration

Immediately downstream of the inception point, substantial free-surface aeration occurs and remains sustained. In the skimming flow, C at a given step edge is approximated using an analytical model based on the advective diffusion of air bubbles (Chanson & Toombes, 2002): (16) $\begin{array}{r}C=1-{\tanh }^2\left(K-\frac{y/y_{90}}{2D}+\frac{{\left(y/y_{90}-1/3\right)}^3}{3D}\right)\end{array}$(16) where K = integration constant and D = function of C_m. They are estimated by (17) $\begin{array}{r}K=\frac{1}{2D}-\frac{8}{81D}+0.3275\end{array}$(17) (18) $\begin{array}{r}{C}_m=0.7622\times \left(1.0434-e^{-3.614D}\right)\end{array}$(18) where C_m = depth-averaged C, obtained by (19) $\begin{array}{r}{C}_m=\frac{1}{y_{90}}\times {\int }_0^{y_{90}}\left(1-C\right)\times dy\end{array}$(19) At typical step edges, Figure 7 presents the comparisons of C among the experiments, CFD and empirical model. The VOF and Mixture models fail, on the whole, to generate accurate results, with appreciable underestimations at some locations. The predictions by the Eulerian and PBMs agree well with the experimental results. The empirical correlation provides seemingly accurate estimations. However, despite its close fit onto the physical data, it cannot be used for practical purposes. For the determination of C, a sufficient number of measurements are required in Equation (19). In addition, predictions at large y/y₉₀ lead to significant errors (excluded in the figures).

Figure 7. Comparison of C among the experimental, CFD and empirical results.

Table 2 presents the error statistics of the models in the C estimations. Among the conventional models, the Eulerian method exhibits the highest accuracy, and the VOF shows the lowest. The PBM approach generates the most accurate predictions, with CC > 0.97. The PBM module significantly improves the effectiveness of the conventional models. The average improvement of the PBM is 6.7% in the CC, 56.5% in the RMSE, and 53.8% in the MAE.

Table 2. Statistical performance of the CFD models for air concentration estimation.

Model	CC	RMSE	MAE	Model	CC	RMSE	MAE
VOF	0.857	0.283	0.165	PBM-Mixture	0.977	0.082	0.050
Mixture	0.942	0.149	0.095	PBM-Eulerian	0.979	0.076	0.048
Eulerian	0.952	0.113	0.058	–	–	–	–

4.2. Error-corrected PBM

The PBM module exhibits superior efficacy over the conventional models in estimating V and C. However, notable insufficiency remains. To this end, the BO optimised ML models are incorporated into the PBM approach for error correction.

4.2.1. Model set-up

The goal of the practice is to establish an error-corrected PBM, denoted as EPBM. For this purpose, the SVM and Bagging are used to estimate the residuals. With a 7:3 ratio, the residuals are divided into a calibration (265) and a validation set (113). The predictors are the original CFD results (V or C) and their corresponding coordinates (x and y). The ML-predicted errors are combined with the PBM results to generate new predictions (V or C). To guarantee high performance, the BO approach optimises the ML architecture. In the SVM, the critical parameters for optimisation are box constraint, kernel scale, and epsilon. In the Bagging, they are the learning cycle, decision split, and leaf node. For their definitions, one is referred to Li and Yang (2022). As optimisation is a stochastic process, the outputs from iterations might differ. Thus, the models are run five times, and the one with the highest accuracy is saved. For each time, the BO operates 30 iterations.

4.2.2. EPBM-predicted velocity

Figure 8 presents the optimisation process of the model parameters. After 20 iterations, the errors reach a relatively low level for all the models. The optimal values of the parameters are shown in Table 3. The optimised models are then used for constructing the EPBM.

Figure 8. Optimisation of hyperparameters for different models in V estimation by the EPBM.

Table 3. Optimal values for hyperparameters in V estimation by the EPBM.

	SVM			Bagging
Hyperparameters	PBM-Mixture	PBM-Eulerian	Hyperparameters	PBM-Mixture	PBM-Eulerian
Box constraint	218.96	318.18	Learning cycle	143	86
Kernel scale	5.0671	4.9567	Decision split	156	264
Epsilon	4.8654 × 10^−4	5.6267 × 10^−3	Leaf node	2	1

Figure 9 compares the V results between the experiments and the EPBM approach at typical steps. All the EPBM models generate accurate results; the differences from the experiments are insignificant. Table 4 shows the statistical performance of the EPBM for velocity estimation. The satisfactory goodness-of-fit indexes and small errors confirm the enhanced prediction accuracy.

Figure 9. Comparison between the experimental and the EPBM-generated velocity profiles.

Table 4. Statistical performance of the EPBM for velocity estimation.

Model	CC	RMSE	MAE
PBM-Mixture-SVM	0.632	0.040	0.021
PBM-Eulerian-SVM	0.589	0.042	0.021
PBM-Mixture-Bagging	0.691	0.038	0.019
PBM-Eulerian-Bagging	0.642	0.041	0.020

4.2.3. EPBM-generated air concentration

In the ML models, the hyperparameters are first optimised. Figure 10 exhibits the optimisation process, showing that, after a number of iterations, all the models achieve an optimal architecture, guaranteeing good prediction performance. Table 5 presents the optimised results.

Figure 10. Optimisation of hyperparameters for different models in estimation of C by the EPBM.

Table 5. Optimal values for hyperparameters in the estimation of C by the EPBM.

	SVM			Bagging
Hyperparameters	PBM-Mixture	PBM-Eulerian	Hyperparameters	PBM-Mixture	PBM-Eulerian
Box constraint	410.96	187.25	Learning cycle	487	377
Kernel scale	1.6145	1.0713	Decision split	171	54
Epsilon	5.81 × 10^−4	3.96 × 10^−4	Leaf node	1	2

Figure 11 compares the C results between the experiments and the EPBM approach at typical steps. At step no. 19, the PBM-Eulerian-Bagging method shows some insufficiency, which is presumably attributed to the high turbulence level. Apart from this, all the EPBM models can accurately capture the C distributions. Table 6 presents the performance statistics of the models. The EPBM approach generates sufficiently accurate predictions, with CC = 0.993–0.998, RMSE = 0.013–0.041 and MAE = 0.005–0.021.

Figure 11. Comparisons of C between the experiments and the EPBM approach.

Table 6. Statistical performance of the EPBM approach for C estimation.

Model	CC	RMSE	MAE
PBM-Mixture-SVM	0.998	0.020	0.010
PBM-Eulerian-SVM	0.999	0.013	0.005
PBM-Mixture-Bagging	0.997	0.028	0.016
PBM-Eulerian-Bagging	0.993	0.041	0.021

4.3. Model comparison

Comparisons are made to evaluate the efficiency of the conventional and modified PBM models. Figure 12 presents the V contours from the experiments, the PBM and the EPBM. Due to the roughness effects, small V values near the pseudo-bottom are expected. However, immediately downstream of the inception point, the experiments exhibit some suspiciously high velocities, which are presumably attributable to measurement errors (Felder & Chanson, 2008). In the transition from water to air, the mixture velocity is reasonably consistent in the y direction, with large values appearing near the surface. In the PBM approach, the near-bottom and near-surface V results are mostly underestimated. The PBM yields reasonable estimations in the transitional regions. In comparison with the PBM, the EPBM gives rise to better estimations. The velocity, both close to the bottom and in the transitional area, agrees well with the laboratory results. Near the surface, the EPBM slightly underestimates the V.

Figure 12. Velocity contours from (a) experiment, (b) PBM-Mixture, (c) PBM-Eulerian, (d) PBM-Mixture-SVM, (e) PBM-Eulerian-SVM, (f) PBM-Mixture-Bagging, and (g) PBM-Eulerian-Bagging.

Figure 13 presents C contours from the experiments, the PBM and the EPBM, where the x-axis runs through the pseudo-bottom (Figure 3). In the experiments, intense air-water mixing occurs some 20 mm above the pseudo-bottom. The PBM results indicate that the air-water mixing zone develops further down towards the pseudo-bottom (smaller y). The EPBM reproduces the aeration extent more satisfactorily, including the turbulent air-water interactions. In the highly aerated regions, both the PBM and the EPBM capture well the general pattern of C distributions.

Figure 13. Air concentration contours from (a) experiment, (b) PBM-Mixture, (c) PBM-Eulerian, (d) PBM-Mixture-SVM, (e) PBM-Eulerian-SVM, (f) PBM-Mixture-Bagging, and (g) PBM-Eulerian-Bagging.

Compared to the PBM models, the improvements of the corresponding EPBM models are statistically presented in Table 7. For velocity estimation, the EPBM considerably improves the accuracy of the PBM, with the CC enhanced by up to 35.5%, the RMSE lowered by up to 76.3%, and the MAE curtailed by up to 78.7%. The average drop in error is 73.0% for the RMSE and 77.1% for the MAE. For aeration prediction, the EPBM approach also demonstrates remarkable enhancement, with an augment in the CC by 1.4–2.1%, and a decline in the RMSE and the MAE by 46.1–82.9% and 56.3–89.6%, respectively. The mean reduction in error is 67.6 % for the RMSE and 73.5% for the MAE.

Table 7. Improvement of the EPBM as compared with the PBM in modelling of V and C.

	V/V90			C
Model	ΔCC (%)	ΔRMSE (%)	ΔMAE (%)	ΔCC (%)	ΔRMSE (%)	ΔMAE (%)
PBM-Mixture-SVM	23.9	−75.0	−76.4	2.1	−75.6	−80.0
PBM-Eulerian-SVM	5.2	−70.0	−76.1	2.0	−82.9	−89.6
PBM-Mixture-Bagging	35.5	−76.3	−78.7	2.0	−65.9	−68.0
PBM-Eulerian-Bagging	14.6	−70.7	−77.3	1.4	−46.1	−56.3
Average	19.8	−73.0	−77.1	1.9	−67.6	−73.5

The Taylor diagram is a mathematical tool for quantifying the degree of correspondence between the modelled and the observed behaviour. Figure 14 presents the Taylor diagrams for V and C. For each model, three statistical parameters are plotted. The CC refers to the azimuthal angle (blue contours). In the simulated field, the centred RMSE is proportional to the distance from the point on the horizontal axis identified as ‘observed’ (pink contours). The standard deviation of the simulated pattern is proportional to the radial distance from the origin (green contours). In the diagram, a model with higher accuracy is expected to generate a prediction closer to the experiment. For the V results, the PBM leads to significant variations and high error levels. For the C results, the PBM indicates a similar issue. In comparison, for both V and C, the EPBM results collapse nearest the experimental data, showing their improvements. It can be stated that the EPBM introduces considerable enhancement of predictive skills.

Figure 14. Taylor diagrams for (a) velocity and (b) air concentration.

In the form of scatter plots, Figure 15 compares the experimental and predicted V and C. For the PBM models, their estimated V results show significant dispersion, and underestimation exists in both models. Similar issues appear in the C estimations. The PBM generates satisfactory predictions only at the low and high C values, i.e. at C = 0–0.3 and C = 0.9–1. For 0.3 < C < 0.9, overestimations are evident at most locations. Although the PBM approach is imprinted with physical rationales, it fails to make accurate predictions for the chute flows examined, especially in the flow regions with intense air-water mixing, which necessitates improvements. The EPBM provides enhanced predictions for the flow parameters. The results agree well with the experimental data, without obvious under- or overestimations. The biases in the PBM become significantly lower in flow regions with both weak and intense air-water interactions.

Figure 15. Scatter plots of the dimensionless interfacial velocity (a and b) and air concentration (c and d).

Figure 16 presents the cumulative frequency (CF) of the relative error (RE = |S − O|/O). For V, the EPBM exhibits comparative accuracy and outperforms the PBM. For instance, 76.5% of the PBM-predicted results show a RE below 10%. Within the same error range in the EPBM, the proportion of data increases to 98.0%. For RE < 20%, the CF is 85.4% for the PBM and 99.6% for the EPBM. Their max. RE is 75.9% and 32.1%, respectively. For C, the EPBM demonstrates superior performance over the PBM. 63.0% of the predictions by the PBM-Eulerian and 64.8% by the PBM-Mixture fall within 0 < RE < 10%. The PBM-Eulerian-SVM results in the highest accuracy, with 95% of its estimations falling within 0 < RE < 20%. For RE < 10% and RE < 20%, the average improvement in CF of the EPBM is 20% and 15%, respectively.

Figure 16. Cumulative frequency (CF) of the relative error (RE) for (a) velocity and (b) air concentration

5. Discussion

High-speed flows in hydraulic structures lead to intense air-water interactions. For safe design and operation, CFD is often used to reproduce the aeration phenomenon. One challenging task in its modelling is to account for the evolutionary process of air bubbles reasonably. The conventional two-phase models (e.g. VOF, Mixture and Eulerian) lack such rationales, thus resulting in limited accuracy. The PBM has shown potential in bubbly flow predictions (Castellano et al., 2018; Wang et al., 2005). The major difference from other industrial applications is the much higher velocity and turbulence level in spillway flows. For improved predictions considering more flow physics, this study incorporates the PBM approach into the conventional two-phase models. The resulting PBM procedure generates predictions of air-water mixing with higher accuracy, which is consistent with previous studies (Lei et al., 2021; Xiang & Tu, 2016). This is plausibly due to the realistic representation of the size distribution and evolution of air bubbles (i.e. breakup, coalescence, growth, dispersion, etc.).

However, despite the enhancement, the PBM still exhibits some insufficiency in producing accurate results in flow regions with strong air-water exchanges, which is also reported by Xiang and Tu (2016). To amend this, the PBM is integrated into the ML technique to create an error-corrected procedure. The resulting EPMB significantly improves the efficacy of the PBM in terms of mixture flow properties. Using ML techniques to achieve error correction for CFD simulations provides a solution for accurately modelling complex systems, as demonstrated by Hanna et al. (2020).

This study develops an improved hybrid model for accurate modelling of air-water flows. The main findings are twofold: (a) the PBM gains in physical rationales and generates reasonable predictions of high-velocity two-phase flows, and (b) ML techniques are a powerful tool to enhance CFD simulations. This study provides a novel solution for establishing error-free digital twins. Due to the complexity of physical systems and sometimes the insufficiency of numerical models, simulation errors are inevitable. Using the ML technique to map the residuals helps create high-performance mathematical models that generate accurate predictions. The robustness of the developed models is validated in terms of flow velocity air concentration. The goal is to achieve reliable calculations for the entire flow field. In future studies, it is desirable to extend the research and construct an EPBM framework that accurately captures all the flow features.

6. Conclusions

A spillway flow is characterised by large discharge, high velocity, strong turbulence and intense air-water mixing, which differs significantly from flows in many other industrial disciplines. For accurate modelling, this study first establishes a PBM approach and compares it with the conventional two-phase models. For further improvement, the PBM is then incorporated into the ML technique to establish the EMPB, in which the BO algorithm is employed for optimisation. The conclusions from the study are summarised as follows.

• The conventional two-phase models (VOF, Mixture and Eulerian) can capture the overall trend of flow velocity and air concentration in the stepped spillway flows. However, compared with the experimental results, the simulations produce appreciable errors if air-water interactions are intense.

• Incorporating the PBM into the conventional two-phase models improves the modelling accuracy. For the flow velocity, the average improvement is 18.9% for the CC, 3.8% for the RMSE, and 8.1% for MAE. For the air concentration, the corresponding indexes are 6.7%, 56.5% and 53.8%, respectively. Despite the enhancement, the flow velocity and air concentration are often underestimated.

• With the error-corrected procedure using machine learning techniques, the EMPB significantly boosts the performance of the PBM. For the flow velocity, the average augmentation is 19.8% in the CC, 73.0% in the RMSE, and 77.1% in the MAE. For the air concentration, the improvement is by 1.9%, 67.6% and 73.5%, respectively.

• The EPBM approach generates the most accurate results for the stepped chute flow. For the flow velocity, 99.6% of the predictions fall within 0 ≤ RE ≤ 20%, and 98.0% within 0 ≤ RE ≤ 10%. For the PBM approach, the corresponding results are 85.4% and 76.5%. For the air concentration, 89.6% of the EPBM estimations are within 0 ≤ RE ≤ 20% and 84.0% within 0 ≤ RE ≤ 10%. The corresponding numbers for the PBM are 81.7% and 74.7%, respectively.

The established EPBM approach is a novel and reliable method for the computation of high-velocity two-phase flows in chute spillways. It also provides a solution procedure for error-corrected CFD modelling. Further research should be devoted to obtaining error-free simulations for the entire flow features. Keeping in mind the complexity in the mathematical formulation of air-water mixing subjected to high velocity, the PBM should be evaluated with other spillway configurations. Due to the highly dynamic nature of the spillway flow, eddy simulations can also be assessed instead of the conventional turbulence models.

Acknowledgement

As part of the research project Quality and trust of numerical modelling of water-air flows for safe spillway discharge (VKU14151), this study is funded by the Swedish Hydropower Centre (SVC). SVC is established by the Swedish Energy Agency, Energiforsk and Svenska Kraftnät, together with the Royal Institute of Technology (KTH), Luleå University of Technology (LTU), Uppsala University (UU) and Chalmers University of Technology (CTH). Participating companies and industry associations include AFRY, Andritz Hydro, Boliden, Fortum Generation, Holmen Energi, Jämtkraft, Karlstads Energi, LKAB, Mälarenergi, Norconsult, Rainpower, Skellefteå Kraft, Sollefteåforsens, Statkraft Sverige, Sweco Energuide, Sweco Infrastructure, Tekniska verken i Linköping, Uniper, Vattenfall R&D, Vattenfall Vattenkraft, Voith Hydro, WSP Sverige and Zinkgruvan. Holger Ecke of Vattenfall R&D, Anders Ansell of KTH, and Carolina Holmberg, Emma Hagner, Lennart Kjellman and Stina Åstrand of SVC are acknowledged for coordination of diverse issues.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Swedish Hydropower Centre [grant number VKU14151].