Bubble columns have reportedly been used in the chemical, petrochemical, bioprocesses, and pharmaceutical industries. In simple bubble column reactors, the gas phase is dispersed into a liquid or liquid-solid continuous phase. In general, depending on superficial velocities and column diameter, the regime inside the bubble column is either homogeneous, transitional, or heterogeneous (Deckwer, 1992). The former involves simpler physics as compared to the latter and most of the models (interphase, coalescence and breakup) were developed in that regime before being later implemented in the heterogeneous regime. The gas holdup as been reported as the most important design criterion in bubble columns. The latter is related to the bubble size, which ultimately allows determining the interfacial area and ultimately, defines the mass transfer phenomena. In biphasic non-reactive bubbly flows, the bubble size varies due to the gas and liquid velocities, inlet geometry, bubble coalescence, bubble breakup, and bubble growth, hence complicating the hydrodynamic behavior inside the system (Fan, 1989; Yeoh et al., 2014). Furthermore, the gas and liquid flow need closure terms in the interphase momentum transfer equations that depend locally on the velocity profiles, physical properties of phases and on the turbulence parameters that are still under development in the open literature (Ishii and Hibiki, 2011). Hence, a comprehensive understanding of the fluid dynamics is required and the latter would in turn be very useful in many industrial fields. Many researchers have used computational fluid dynamics (CFD) techniques to simulate biphasic bubble columns. The latter are in most cases simulated by the Euler-Euler approach (two-fluid) due to a lesser computational cost when compared to the Euler-Lagrange or volume of fluid (VOF) approaches. In gas-liquid flow, the interface involves both drag and non-drag forces.

The drag force has an influence on the macroscopic structure of the flow. For instance, the radial profiles of velocity and holdup depend on the drag coefficient, Reynolds number, Eotvos numbers, terminal velocity and on the physical properties of the continuous phase (Wang and Yao, 2016). Rzehak et al. (2017) simulated a bubbly flow in different operating conditions and the geometries using the Ishii drag coefficient (Ishii and Zuber, 1979) and the predicted results were reported to be in good agreement with experimental values. This drag correlation is suitable for a wider range of bubbles sizes and covers all flow regimes (homogeneous, transitional or heterogeneous).

In bubble columns, the shape of the radial profiles may change according to the net lateral lift force. According to Tomiyama (2004), small bubbles (d_b < 5.8 mm) have a positive lift coefficient and tend to go toward the reactor wall. However, larger bubbles have a negative lift coefficient and tend to stay at the core of the bubble column. Zhang et al. (2005, 2006) suggested that the inclusion of the Tomiyama lift coefficient could predict a better correlation with experimental values. Nevertheless, Masood and Delgado (2014) and Yamoah et al. (2015) studied the influence of wall lubrication force and found that the Tomiyama correlation (Tomiyama et al., 1995) tends to over-estimate the velocity profiles when compared to the Antal correlation (Antal et al., 1991) that however agrees well with experimental data. Finally, Lucas et al. (2007) developed a 1D-model, studying the effect of wall lubrication and turbulent dispersion forces and suggested that the combination of these non-drag forces provides reasonable results.

Laborde-Boutet et al. (2009) focused on different formulations of the turbulent model and showed that the Renormalization Group (RNG k-ε) predicts higher turbulent dissipation rate as compared to Standard k-ε, which is known to be underestimated (Jakobsen et al., 2005). Most of the coalescence and breakup kernels depend on the local turbulent dissipation rate. Therefore, Laborde-Boutet suggested that the RNG k-ε turbulent model should be used to implement a population balance model (PBM).

The bubble coalescence and breakup phenomenon requires a Population Balance Equations (PBE) which allows discretization into N-classes of bubble size that can then be coupled with a two-fluid model. A single momentum equation can be solved for all N-classes, an approach called Homogeneous multi size group (MUSIG). In the case of non-homogeneous multi size group (called iMUSIG), multiple momentum equations are solved, making the solutions computationally costly. Krepper et al. (2008) developed and worked on the simulation of a gas-liquid phase using an iMUSIG model and suggested that 2–3 subgroups are sufficient to capture the fluid behavior. It was concluded that although iMUSIG allows simulating the local radial profiles, it is still limited by breakup and coalescence kernels that use an isotropic turbulent approach. Similarly, Xu et al. (2013) simulated a bubble column using both MUSIG, and iMUSIG. The results showed that the former (which includes lift force) and the latter both agreed well with experimental results. Wang et al. (2003, 2005) studied the effect of different coalescence and breakup kernels and results showed that the Luo breakup (Luo and Svendsen, 1996) predicts lower breakup rates, while the Lehr kernel (Lehr et al., 2002) predicts higher breakup rates with regards to empirical values. The key difference in both kernels is in the estimation of the breakup efficiency. Luo's model includes the surface energy constraint, which shows that the break-up could only occur if the kinetic energy of the colliding eddies is higher than the surface energy necessary for bubble breakage. However, Lehr's model only considers the capillary constraint, assuming that the interfacial and inertial forces balance each other. Chen et al. (2005) also studied the effect of different kernels and concluded that the radial profiles were not sensitive as long as the breakup is increased 10 times. Xu et al. (2013) used Luo's model for bubble coalescence and breakup and modified the coalescence parameter to 0.5, generating results that were in good agreement with experimental data.

In literature, the most commonly used kernel for bubble coalescence is the Luo's model while for bubble breakup, the Luo and Lehr model are usually preferred. Luo's coalescence model over-predicts the collision frequency and needs adjustment to reduce the coalescence rate which can be achieved by tuning the coalescence parameter (Wang et al., 2005; Yeoh et al., 2014). Finally, the effect of coalescence parameter in Luo's model was not extensively reported in the open literature and only a handful of studies have been published so far (such as Xu et al., 2013, 2014). Furthermore, the influence of this coalescence parameter on radial profiles of gas holdup and axial liquid velocity using two different bubble breakup models is limited.

In light of this, this work intends to fill the gaps identified in the previously reported approaches with the specific target to fit with industrial applications. Hence, the main objectives of this paper are as follows:

▪ Investigate the influence of the coalescence parameter on radial profiles using a combination of Luo coalescence and Luo breakup models.

The presented results for a biphasic bubbling column were generated on a 2D-axisymmetric geometry and the predicted time-averaged profiles were compared with the literature data of Hills (1974). CFD-simulations were conducted using the commercial software ANSYS-Fluent v.17.2.

All simulations were run on a 2D axis-symmetric geometry. The assumption for the 2D axis-symmetric could be reasonable since experimental data reported by Hills (1974) and Degaleesan (1997) showed that the time-averaged flow field produces a stationary axis-symmetric flow pattern, hence supporting the validity of the 2D model. Simulations were validated using the experimental data published by Hills (1974), which was shown to be robust and is often used by other authors (Krishna et al., 1999; Van Baten, 2000; Ekambara and Joshi, 2005). Hills data has been extensively cited in literature explaining why the model developed in this work was validated using these empirical values. The two-fluids involved in the experiments consisted of air (acting as disperse phase) and water (considered as the continuous phase). The superficial gas velocity was varied between 0.019 and 0.038 m/s, range in which an homogeneous regime could be achieved (Krishna et al., 1999). The diameter and height of the cylindrical column were of 0.138 and 1.38 m, respectively. The static liquid height was 0.9 m and all the experimental observations were performed at a 0.6 m height. The inlet geometry of the experimental setup consisted of a perforated plate with 61 holes which all had a 0.0004 m diameter. Due to the limitation associated with the mesh size and computational cost, the gas was assumed to be introduced uniformly from the bottom of the column. This assumption was supported by Buwa and Ranade (2002) where the influence of the sparger design using a perforated plate (actual experimental inlet with holes) and the sintered plate was investigated and was shown to induce no significant difference with regards to empirical data. They also concluded that a hole diameter of 0.8 mm requires in turn a very fine meshing in the simulations, making it computationally very expensive. Similarly, Chen et al. (2005) also simulated a bubble column using a sintered plate instead of a perforated plate and reported that it is not essential to use the actual experimental inlet configuration.

The boundary condition involves a uniform inlet bubble size which was calculated from Kumar's correlation and was obtained for diameters of 3.6 mm and 4.5 mm at superficial gas velocities of 0.019 and 0.038 m/s, respectively (Kumar et al., 1976). The outlet and wall include atmospheric pressure and non-slip boundary conditions, respectively. Gas was the only mixture introduced from the inlet (α_p = 1), while the 0.9 m static column height involved α_q = 1 and α_p = 0. Above this level (free board), the gas and water phase fractions were α_q = 0 and α_p = 1, respectively. The bubble volume of each class was calculated from the following formula (vi+1/vi=2r), where r is the ratio factor which equals to 1, 2 … n. For all the simulations (except for the mesh sensitivity analysis), a third order upwind scheme was used to discretize the continuity equation while the rest of the transport equations were solved by a second order scheme (see Table 1). The mesh sensitivity analysis was performed using a first order scheme due to convergence issues that were faced when solving the transport equations with higher order schemes for finer mesh. Convergence problems were encountered when an adaptive time step approach was used. In such cases, solutions tended to diverge due to the variation of the time step, especially at the initial flow time. The fixed time step was well consistent in term of convergence. Hence, 1E-04 s time steps were used and guaranteed that the courant number for air and water velocities was <1. Once a statistically steady state was reached, a time-averaged sampling was calculated for 30 s.

In this work, 2D-axisymmetric simulations of a bubbling column were performed and the simulated time-averaged radial profiles were compared with the empirical data obtained from Hills (1974). The investigated superficial gas velocities were 0.019 and 0.038 m/s, covering the homogeneous bubbly regime. The developed model consisted of a two-fluid model coupled with a PBM. The former included the gas-liquid interface that considered the drag, lift, wall lubrication and turbulent dispersion forces. The latter involved the Luo bubble coalescence model as well as two different bubble breakup models: Luo and Lehr. From this work, the following conclusions could hence be formulated:

The sensitivity analysis of the bubble classes was performed using 14, 20, and 22 bins and it was shown that the solution was independent of the bubble classes. Hence, the lowest number of bubble classes were selected to reduce computational cost.

Verification of the numerical schemes was performed using first and second orders and results showed that numerical schemes had no significant influence on the predicted radial profiles of gas holdup and axial liquid velocity. However, at the center of the column, a slight discrepancy was observed, which might be related to numerical diffusion.

Mesh sensitivity was conducted on 1.5 mm (finer), 3 mm (medium), and 6 mm (coarse) mesh sizes. The predicted axial liquid profile of coarse (6 mm) mesh size differed from medium and fine mesh size, and hence was ignored. The fine meshing showed a non-realistic behavior near the wall without the inclusion of wall lubrication force, which might be related to the y⁺ value of 12.5 and once introduced, there was no significant difference between fine- and medium-sized mesh. In addition, the predicted axial liquid velocity slightly improved at the core of the bubble column.

The combination of the Luo coalescence and Luo breakup kernels (Luo-Luo) was shown to under-predict the gas holdup both at a 0.019 and 0.038 m/s superficial gas velocity. The gas holdup was increased to a maximum when the coalescence parameter was reduced. However, at the lowest Ug and the c₀ (=0.1) values, the predicted velocity profile was far away from the experimental values. It is thus recommended to tune the coalescence parameter when using the Luo-Luo kernels and both the holdup and axial liquid profiles should be considered for validation purposes.

Simulations using a combination of the Luo coalescence and Lehr breakup kernels (Luo-Lehr) predicted a closer holdup both for the 0.019 and 0.038 m/s superficial velocities when compared with experiments. Scaling of coalescence parameter, in combination with the Lehr model leads to no significant improvement in the gas holdup. Furthermore, a decrease of the coalescence parameter significantly influences the axial liquid profile that under-predicts the profile compared to experiments. Results have shown that it is better to use Luo-Lehr kernels without any modification of the coalescence parameter.

AS: Wrote and drafted the article. MB, TM: Helped in developing the model and provided the technical guidelines. JL: Reviewed and approved the article, managed the research, reviewed the results and provided the technical guidelines.