Investigation of Plume Offset Characteristics in Bubble Columns by Euler–Euler Simulation

Abstract

Based on low-cost and easy to enlarge, the bubble column device has been widely concerned in chemical industry. This paper focuses on bubble plumes in laboratory-scale three-dimensional rectangular air-water columns. Static behavior has been investigated in many experiments and simulations, and our present investigations consider the dynamic behavior of bubble plume offset in three dimensions. The investigations are conducted with a set of closure models by the Euler–Euler approach, and subsequently, literature data for rectangular bubble columns are analyzed for comparison purposes. Moreover, the transient evolution characteristics of the bubble plume in the bubble column and the gas phase distribution in sections are introduced, and the offset characteristics and the oscillation period of the plume are analyzed. In addition, the distributions of the vector diagram of velocity and vortex intensity in the domain are given. The effects of different fluxes and column aspect ratios on bubble plumes are studied, and the offset and plume oscillation period (POP) characteristics of bubbles are examined. The investigations reveal quantitative correlations of operating conditions (gas volume flux) and aspect ratios that have not been reported so far, and the simulated and experimental POP results agree well. An interesting phenomenon is that POP does not occur under conditions of a high flux and aspect ratio, and the corresponding prediction values for the conditions with and without POP are given as well. The results reported in this paper may open up a new way for further study of the mass transfer of bubble plumes and development of chemical equipment.

1. Introduction

The bubble column is widely used in the chemical industry [1,2]. It is a typical gas–liquid two-phase flow system [3]. Among them, the water phase exists as a continuous phase, and the air phase as a discrete phase rises from the reactor in the form of bubbles [4]. Bubbles are easy to observe and much attention has been paid to instrument measurement. Bubbles, namely, discrete elements, have always been a hot and challenging topic in both experiments and simulations [5]. Difficulties in bubble columns stem from the fact that the disperse phase is characterized by a complex behavior [6]. Consequently, the design and optimization of bubbly flow equipment present a great challenge to improve the mass transfer efficiency [7]. Although the industrial application of bubble columns has been extensively used, there are still some unsolved problems in their design and amplification, mainly due to the unknown problems in the hydrodynamics of gas–liquid flow [8]. Therefore, we focus our attention on the bubble plume. The bubble plume is the movement of bubbles in a two-phase flow, which is caused by momentum exchange. The bubble plume can also disturb the ambient flow field, and the main objectives are to promote the mixing of the liquids and to enhance the mass transfer. However, due to the complexity of two-phase flow, especially the appearance of bubble motion, the basic properties of gas–liquid hydrodynamics and the characteristics of a bubble plume are still limited [9].

There are many works in the literature on the study of bubble plumes in a bubble column by numerical simulation via bubbly flow observation [10,11,12]. Different numerical methods have been proposed to describe a bubble plume, mainly including the bubble shape, state and plume oscillation [13,14]. The dependence of the period of bubble plume oscillations, global gas hold-up and local phase velocity on the superficial gas velocity and aspect ratio have previously been studied [15]. Diaz, Montes [16] focused on an oscillating plume in terms of the aspect ratio and gas flow rate. It was found that the oscillation period did not depend on the aspect ratio if the latter exceeded a value of 2. They investigated the dependency of the superficial gas velocity and aspect ratio on the time-averaged gas hold-up as well. With the exception of surface tension, Cachaza, Elena Díaz [17] showed that although there are differences in the mechanism of gas–liquid interface properties, the state of bubble flow, including the plume oscillation period, can be unified modeling and analysis. These are the critical problems in the design of bubble column equipment. Bannari, Kerdouss [18] effectively obtain bubble rise and oscillation characteristics, which were in good agreement with experimental results with the aspect ratio of 2.25 and a superficial gas velocity of 0.73 cm/s, and the oscillation characteristics of bubble plume in aspect ratio of 4.5 are analyzed and predicted. Gupta and Roy [19] captured bubble plume oscillation by radioactive particle tracking, and added a series of various interfacial forces (drag, lift, and virtual mass) to the Euler–Euler model and conducted a comparative analysis of different combinations. They emphasized the importance of low superficial gas velocity in the prediction of plume oscillations in transient simulations. Masood and Delgado [20] compared different drag models to bubble plume oscillation in 3D columns. Moreover, the drag model will affect the predicted oscillation period in the aspect ratio of 3. Previous studies have also closely considered interaction and turbulence simulation models to effectively simulate bubble plume, each with a different closure relation model. Liu and Luo [21] found the effect of surface tension on the oscillation behavior of bubble plume is obvious. In addition, the gas velocity is also very important for the prediction of plume oscillation. The velocities and turbulence characteristics of bubble columns with different aspect ratios were studied, and the snapshots of oscillation and turbulent kinetic energy at different flow rates and aspect ratios were provided, and they offered available data for theoretical studies [22]. Fleck and Rzehak [23] compiled an extensive summary of the plume oscillation characteristics and successful prediction of the gas–liquid two-phase characteristics in the bubble column with an aspect ratio of 2.5, and a qualitative study of the plume dynamics. The effects of superficial gas velocity, aspect ratio, reactor and liquid viscosity on low frequency oscillations were studied [24]. When the aspect ratio is greater than 2.0, the plume oscillation frequency remains almost unchanged, and as the aspect ratio decreases from 2.0 to 1.375, the plume oscillation frequency decreases significantly. Different CFD-PBM frameworks have been proposed to describe effectively simulate the plume oscillation period [25,26]. The polydispersed bubble size is used to effectively predict the plume oscillation period with the aspect ratio of 2.5 [27]. Shang, Ng [28] made use of two sets of momentum closures, and it was found that the momentum closure in the numerical simulation has a dramatic effect on the bubble plume oscillation characteristics. Moreover, there have also been reports of experiments to study the oscillation characteristics of bubble plumes [29]. Several available measuring techniques and instruments (i.e., laser Doppler anemometry (LDA), particle image velocimetry (PIV), laser-induced fluorescence (LIF) and computer-aided radioactive particle tracking (CARPT)) to investigate bubble rise and processes of bubble plume oscillations have been applied [11,30,31].

By sorting out a lot of previous research work, we found that more attention has been paid to the oscillation period of a bubble plume, and few studies focus on the offset characteristics of the bubble behavior. Furthermore, we also noticed that there are many comparisons between simulations and experiments, and the emphasis is on constructing numerical simulation methods for bubble plumes. For instance, different models for the description of turbulence have previously been studied [32], which also showed a notable dependence of the κ-ε turbulence model. Yang, Zhang [33] suggested that the resistance correction has a strong ability to simulate the plume oscillation, which should be based on the size distribution of bubbles and the local gas holdup. The application of bubble offset is of great significance to the further study of gas–liquid mixing and mass transfer behavior, while the plume oscillation and offset characteristics correspond to each other. Based on these, the offset characteristics of bubble plume fully developed in bubble column are analyzed in detail.

These issues are illustrated in this work. The main purpose of the present work is to analyze the mixing of gas–liquid and oscillation characteristics of dynamic bubbly flows (periodic bubble plumes). In particular, the dynamic rise process of bubbles and the offset behavior of bubbly flow in a rectangular column are examined. The feasibility of the bubbly flow simulation method is demonstrated with a practical case in rectangular bubble columns. Based on available experimental data and our own simulations, correlations are proposed for plume offset characteristics in dimensionless form. In addition, a mechanism for generating no-plume oscillation period (POP) is given.

In this paper, previous experimental data are first described in Section 2. In Section 3, numerical setup implementation details are provided, all the hydrodynamic equations are derived and computational models are reviewed. Then, the simulation results are presented and compared to experimental data in Section 4. Changes in the plume oscillation position over time and space are provided as well. Finally, conclusions are drawn in Section 5.

2. Experimental Data in the Literature

A series of experiments investigated bubbly flow in rectangular bubble columns, where tap water for the liquid phase and air for the gas phase were used. These experiments are listed in

Table 1. Details of the experimental data in the literature.

Refs.	W-D-H, cm	Aspect Ratio	Gas Distributor	Simplified Structure, mm	Gas Volume Flux, L/h	DB, mm
[10]	20–5–45	2.25	sparger (8 holes)	rectangle (24 × 12)	20–90	1–10
[4]	20–4–45	2.25	single-orifice hole	diameter of 1	48	1–10
[34]	20–5–120	6	sparger (8 holes)	rectangle (18 × 6)	56–296	1–10
[16]	20–4–180	1.25–2.25	sparger (8 holes)	none	69–613	1–10
[13]	20–4–45	2.25	sparger (8 holes)	rectangle (18 × 6)	69–613	1–10
[22]	20–5–120	1.05–4	sparger (8 holes)	none	48–600	1–10
[35]	26.7–1.5–50	1.87	needle	diameter of 0.4	3–12	1.5–2.5

, and experimental data for different column sizes and aspect ratios, reactor structures (including a simplified structure in simulations) and operating conditions were gathered from the literature. Columns [10] were maintained at 25 °C and atmospheric pressure, while most of the experiments were performed on gas spargers with several holes, including pitches. The holes were located at the center position of the sparger. Many of the spargers produced a broad bubble size distribution (1–10 mm) with a mean size of 5 mm. The gas volume flux was measured by volumetric flow meters, varying from 3 to 745.5 liters per hour (LPH). All the experimental data were considered to examine the bubbly flow oscillation and offset characteristics.

By summarizing the above experiments, we have a more specific understanding of the bubble plume in the bubble column. The size of bubble column varies, and the volume flux is also different. However, the range of bubble size in the domain is roughly the same. Of course, it is related to the reactor and flux, which is significant for the selection of bubble size below. In order to study the offset characteristics of plume oscillation more accurately, we control the flux range of this work within the experimental data, which has more practical significance. This is also one of the purposes of gathering experimental data. However, in the experiment, the snapshot of the gas phase is only taken instantaneously, which is not useful in this paper. We extract the period of plume oscillation in the literature and compare it with the results. In summary, we can get the proper flux, bubble size distribution, size of device and distributor, and part data of plume oscillation periods from these experiments.

3. Models and Numerical Details

Different TFMs are repeatedly cited in many previous works. Therefore, a concise description is provided here. References to the original works are presented as well. The conservation equations of two-phase flow are listed in

Table 2. Summary of the governing equations for two-phase flow.

Models	Equations	Refs./Remarks
Conservation equations	$\frac{\partial }{\partial t}({\rho }_i{\alpha }_i)+\nabla \cdot ({\rho }_i{\alpha }_i{\mathbf{u}}_i)=0$ $\frac{\partial }{\partial t}({\rho }_i{\alpha }_i{\mathbf{u}}_i)+\nabla \cdot ({\rho }_i{\alpha }_i{\mathbf{u}}_i{\mathbf{u}}_i)=\nabla \cdot ({\alpha }_i{\mathbf{T}}_i)-{\alpha }_i\nabla {\mathrm{p}}_i+{\alpha }_i{\rho }_i\mathbf{g}+{\mathbf{F}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ ${\mathbf{T}}_i={\mu }_i^e\left(\nabla {\mathbf{u}}_i+{\left(\nabla {\mathbf{u}}_i\right)}^T\right)$	$i=g,l$
Turbulence Equations for water phase	$\frac{\partial }{\partial \mathrm{t}}\left({\mathrm{k}}_l{\rho }_l\right)+\nabla \cdot {(\mathrm{k}}_l{\rho }_l{\alpha }_l{\mathbf{u}}_l)=\frac{\partial }{\partial x_i}\left[{\alpha }_l\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+{\alpha }_l\left(G_k+G_b-\rho \epsilon -Y_M\right)$ $\begin{array}{l}\frac{\partial }{\partial \mathrm{t}}\left({\epsilon }_l{\rho }_l\right)+\nabla \cdot ({\epsilon }_l{\rho }_l{\alpha }_l{\mathbf{u}}_l)=\frac{\partial }{\partial x_i}\left[{\alpha }_l\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial {\epsilon }_l}{\partial x_i}\right]+{\alpha }_l{G}_{1\epsilon }\frac{{\epsilon }_l}{k_l}\left(G_k+C_{3\epsilon }{G}_b\right)\\ -{\alpha }_l{C}_{2\epsilon }\rho \frac{{\epsilon }_l^2}{k_l}\end{array}$${\mathsf{\mu }}_{\mathrm{l}}^{\mathrm{e}}={\mathsf{\mu }}_{\mathrm{t},\mathrm{l}}+{\mathsf{\mu }}_{\mathrm{l}\mathrm{a}\mathrm{m},\mathrm{l}}=\mathsf{\rho }{\mathrm{C}}_{\mathrm{u}}\frac{{\mathrm{k}}_{\mathrm{l}}{}^2}{{\mathsf{\epsilon }}_{\mathrm{l}}}+{\mathsf{\mu }}_{\mathrm{l}\mathrm{a}\mathrm{m},\mathrm{l}}$	$\begin{array}{l}{\mathrm{C}}_{1\mathsf{\epsilon }}=1.44,{\mathrm{C}}_{2\mathsf{\epsilon }}=1.92,\\ {\mathrm{C}}_{3\mathsf{\epsilon }}=0.09,{\mathsf{\sigma }}_{\mathrm{k}}=1.0,\\ {\mathsf{\sigma }}_{\mathsf{\epsilon }}=1.3\end{array}$
Turbulence viscosity of air phase	${\mathsf{\mu }}_{\mathrm{g}}^{\mathrm{e}}={\mathsf{\mu }}_{\mathrm{t},\mathrm{g}}+{\mathsf{\mu }}_{\mathrm{l}\mathrm{a}\mathrm{m},\mathrm{g}}={\mathsf{\mu }}_{\mathrm{t},\mathrm{l}}\times {\mathsf{\rho }}_{\mathrm{g}}/{\mathsf{\rho }}_{\mathrm{l}}+{\mathsf{\mu }}_{\mathrm{l}\mathrm{a}\mathrm{m},\mathrm{g}}$	-
Interphase forces	${\mathbf{F}}_l^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}=-{\mathbf{F}}_g^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}},{\mathbf{F}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}={\mathbf{F}}_{Di}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}+{\mathbf{F}}_{Li}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}+{\mathbf{F}}_{Wi}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}+{\mathbf{F}}_{Ti}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}+{\mathbf{F}}_{VMi}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$	$i=g,l$
Drag forces	$\begin{array}{l}{\mathbf{F}}_{\mathrm{D}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}=\frac{3}{4}{\mathsf{\alpha }}_{\mathrm{g}}{\mathsf{\rho }}_{\mathrm{l}}{\mathrm{C}}_{\mathrm{D}}\frac{1}{{\mathrm{d}}_{\mathrm{g}}}\left|{\mathbf{u}}_{\mathrm{g}-}{\mathbf{u}}_{\mathrm{l}}\right|\left({\mathbf{u}}_{\mathrm{g}-}{\mathbf{u}}_{\mathrm{l}}\right)\\ {\mathrm{C}}_{\mathrm{D}}=\mathrm{M}\mathrm{A}\mathrm{X}\left[\frac{24}{\mathrm{R}\mathrm{e}}\left(1+0.15\mathrm{R}{\mathrm{e}}^{0.687}\right),\frac{8}{3}\left(\frac{\mathrm{E}\mathrm{o}}{\mathrm{E}\mathrm{o}+4}\right)\right]\end{array}$	[36]
Transverse lift forces	$\begin{array}{l}{\mathbf{F}}_{\mathrm{L}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}=-C_l{\rho }_l\left({\mathbf{u}}_g-{\mathbf{u}}_l\right)rot{\mathbf{u}}_l\\ {\mathrm{C}}_{\mathrm{l}}=\{\begin{array}{c}\mathrm{m}\mathrm{i}\mathrm{n}\left[0.288\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(0.121\mathrm{R}\mathrm{e}\right),\mathrm{f}\left(\mathrm{E}{\mathrm{o}}_{\mathrm{d}}\right)\right]\\ \mathrm{f}\left(\mathrm{E}{\mathrm{o}}_{\mathrm{d}}\right)\\ -0.27\end{array}\hspace{1em}\hspace{1em}\begin{array}{c}\mathrm{E}{\mathrm{o}}_{\mathrm{d}}<4\\ 4<\mathrm{E}{\mathrm{o}}_{\mathrm{d}}<10.7\\ \mathrm{E}{\mathrm{o}}_{\mathrm{d}}>10.7\end{array}\\ \mathrm{f}\left(\mathrm{E}{\mathrm{o}}_{\mathrm{d}}\right)=0.00105\mathrm{E}{\mathrm{o}}_{\mathrm{d}}{}^3-0.0159\mathrm{E}{\mathrm{o}}_{\mathrm{d}}{}^2-0.0204\mathrm{E}{\mathrm{o}}_{\mathrm{d}}+0.474\\ \mathrm{E}{\mathrm{o}}_{\mathrm{d}}=\frac{\mathrm{g}\left({\mathsf{\rho }}_{\mathrm{l}}-{\mathsf{\rho }}_{\mathrm{g}}\right){\mathrm{d}}_{\mathrm{H}}^2}{\mathsf{\sigma }}\end{array}$	[37]
Wall lubrication forces	$\begin{array}{r}{\mathbf{F}}_{\mathrm{W}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}=-\left({\mathrm{C}}_{{\mathrm{w}}_1}+{\mathrm{C}}_{{\mathrm{w}}_2}\frac{{\mathrm{D}}_{\mathrm{s}}}{{\mathrm{y}}_{\mathrm{w}}}\right)\frac{{\left[\left({\mathsf{\alpha }}_{\mathrm{g}}{\mathsf{\rho }}_l{\mathbf{u}}_{\mathrm{g}}-{\mathbf{u}}_{\mathrm{l}}\right)-\left(\left({\mathbf{u}}_{\mathrm{g}}-{\mathbf{u}}_{\mathrm{l}}\right)\cdot {\mathbf{n}}_{\mathsf{\omega }}\right){\mathbf{n}}_{\mathsf{\omega }}\right]}^2}{{\mathrm{D}}_{\mathrm{s}}},\\ {\mathrm{C}}_{{\mathrm{w}}_1}=0.0064.,{\mathrm{C}}_{{\mathrm{w}}_2}=0.016\end{array}$	[6]
Turbulent dispersion forces	${\mathbf{F}}_{\mathrm{T}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}={\mathrm{C}}_{\mathrm{T}\mathrm{D}}{\mathrm{C}}_{\mathrm{D}}\frac{{\mathsf{\mu }}_{\mathrm{t}\mathrm{g}}}{{\mathsf{\sigma }}_{\mathrm{t}\mathrm{g}}}\left(\frac{\mathsf{Ñ}{\mathsf{\alpha }}_{\mathrm{l}}}{{\mathsf{\alpha }}_{\mathrm{l}}}-\frac{\mathsf{Ñ}{\mathsf{\alpha }}_{\mathrm{g}}}{{\mathsf{\alpha }}_{\mathrm{g}}}\right),{\mathrm{C}}_{\mathrm{T}\mathrm{D}}=1,{\sigma }_{\mathrm{t}\mathrm{g}}=0.9$	[38]
Virtual mass forces	${\mathbf{F}}_{\mathrm{V}\mathrm{M}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}={\mathrm{C}}_{\mathrm{V}\mathrm{M}}{\mathsf{\rho }}_{\mathrm{l}}{\mathsf{\alpha }}_{\mathrm{g}}\left(\frac{\mathrm{D}{\mathbf{u}}_{\mathrm{g}}}{\mathrm{D}\mathrm{t}}-\frac{\mathrm{D}{\mathbf{u}}_{\mathrm{l}}}{\mathrm{D}\mathrm{t}}\right)$	CVM=0.5
Inlet	$\mathbf{u}=x\mathbf{i}+y\mathbf{j}+z\mathbf{k}=\mathrm{C}\mathbf{j}$	C is Constant value
Outlet	$\frac{\partial u}{\partial y}=0,\frac{\partial p}{\partial y}=0$	-
Wall	$\mathbf{u}=0$	-
Convergence criterion	$\left|\frac{\mathsf{\phi }\left(\mathrm{n}+1\right)-\mathsf{\phi }\left(\mathrm{n}\right)}{\mathsf{\phi }\left(\mathrm{n}+1\right)}\right|\le 0.001$	n, n + 1 are the steps of iterations

. The surface tension and material properties of tap water and air at atmospheric pressure and temperature are summarized in

Table 3. Details of the physical properties and environmental conditions.

Materials	Density, kg/m3	Viscosity, Pa·s	Temperature, °C	Surface Tension, N/m
air	1.225	1.79 × 10−5	25	0.0725
tap water	998.2	1.01 × 10−3

. Numerical details suitable for this work are introduced later.

With regard to the conservation equations, each phase is assumed to be incompressible, and the volume fractions of the two phases always add up to 1 under all circumstances. Interphase forces are exchanged for momentum transfer between the two phases. For the liquid phase, the turbulent dynamic energy and energy rate are calculated from the standard κ-ε turbulence (SKE) model. To obtain a reliable solution of the turbulence equations, the dispersed method is adopted, and standard wall functions are applied to describe the wall–turbulence interactions.

The simulations refer to a bubble column [10] 0.2 m wide, 1.8 m high and 0.04 m deep, as previously tested [13]. The height of the apparatus can be arbitrarily adjusted. The column is filled with water up to a certain height (such as 0.45 m) from the bottom, and the water volume is regarded as the computational domain. Air is injected through a sparger in the center position of the bottom of the domain containing eight holes (diameter: 1 mm) and a 6 mm pitch. Moreover, the sparger is simplified as a rectangular region of 18 × 6 mm size in the simulations. An inlet boundary condition is prescribed on the simplified area, while the holes are treated as a whole. An outlet boundary condition is applied at the top of the domain, and the pressure at the outlet is uncertain. On the remainder of the domain boundary area, a free-slip condition is employed for the air phase and a no-slip condition for the water phase. The water phase is fully turbulent (at a turbulence intensity equal to 5%), and the hydrodynamic diameter is set as the width of the column. For the dispersed phase (i.e., the air phase), the gas velocity is converted by the gas volume flux through the inlet region, and no backflow occurs in the domain. The bubble size is set to 5.0 mm for the sparger on the basis of previous experiments.

The convergence criterion is set to 0.001 for the simulation cases to ensure a high calculation accuracy. Spatial discretization schemes are critical for the solution process. Therefore, the second order upwind scheme is adopted for κ and ε. The time step is set to 0.01 s, which guarantees that the Courant–Friedrichs–Lewy (CFL) number always remains below 0.5. All the simulations are performed with commercial software ANSYS FLUENT 16.0 and conducted on a 2.10-GHz platform (56-core CPUs) with 192 GB of RAM.

The simulated bubble column is shown in Figure 1. The position of the maximum lateral displacement is indicated by the yellow bubble. The partial working conditions chosen for the simulations are listed in

Table 4. Details of the partial working conditions adopted in the simulation cases.

Test Case (Partial)	Aspect Ratio	Gas Volume Flux, l/h	DB, mm	dB, mm	Ug, mm/s
ξ1R1	2	136	1–10	5	3.77
ξ2R2	2.5	226	1–10	5	6.26
ξ3R3	3	317	1–10	5	8.78
ξ4R4	3.5	407	1–10	5	11.28
ξ5R5	4	497	1–10	5	13.78
ξ6R6	4.5	588	1–10	5	16.29

. Six values are chosen for the aspect ratio (ξ_i) and gas volume flux (R_i), which yields 36 groups of test cases by cross-combination. The dynamic flow properties during the plume oscillation period are closely considered in all simulations.

The position of the knee point of the first oscillation during bubbly flow is used to illustrate the offset characteristics of the plume. The maximum offset distance and angle of the bubble plume in dimensionless form are calculated by using the following equations:

$$\mathsf{\eta }=\frac{2L}{W}$$

(1)

$$\mathsf{\theta }=\mathrm{arttan}\frac{L}{h}$$

(2)

where L is the maximum offset distance of the bubble, h is the vertical distance of the maximum offset, and η and θ are the maximum offset distance and angle, respectively. As a result, the calculated η value is not more than 1.

4. Simulation Results and Discussion

All simulations are presented in this section, and they have been run in the transient mode across the full 3D domain of the column. The results are calculated over a simulated physical time of 300 s, which is long enough to ensure full domain development, and the resulting domain no longer exhibits notable fluctuations. Since the different studies described in Section 2 provide different sizes and parameters, the gas inlet type must adopt the volume flux to ensure that the simulations are meaningful while enabling a comparison to experiments which is as accurate as possible.

4.1. Model Validation

As mentioned above, a grid independence study is required to obtain the optimal balance between the computational time and numerical accuracy. Non-uniform hexahedral grids are considered in this work. The results of the grid independence study are not provided, and further details can be found in [13]. Calculations were conducted on fine grids to ensure an adequate resolution. The simulation case against which the approach was validated [10] compared velocity profiles to experimental profiles. The POP data are also assessed via a comparison to the simulation results.

Figure 2 shows the results obtained by using the test case compared to the experimental data with R = 48 L h⁻¹. It is worth noting here that the depth of the experiment is 0.05 m, and the depth of this model is approximately the same. Regardless of whether the model is an experimental or simulation model, it is not a real 3D model (since the depth of the column is very small). Therefore, a comparison is feasible. The agreement between the experiments and simulations is quite good.

The predicted values of the middle position are slightly different, but the overall trend is consistent, which is due to the intense momentum exchange between the gas and liquid phases in the intermediate region. Furthermore, the POP obtained with this gas volume flux is 14.6 s, and the prediction model captures this feature, while the calculated POP almost agrees with the predicted values in previous studies.

4.2. Transient Evolution of Bubble Plume

In handling the dynamics of bubbly flow, in contrast to the time-averaged flow problem, the rise process of bubbles is required for two-phase flow. Few studies have addressed the variations in bubble plume oscillation with simulations and measurements. Furthermore, most are two-dimensional dynamic studies. To better understand the continuous behavior at different times, the performance of the bubble plume has been determined. Predicted values are obtained through simulations with the use of three-dimensional rectangular columns.

The instantaneous gas fraction ranging from 0 to 0.1 is shown in the test case of ξ3R3, and the liquid velocity fields are also shown in Figure 3. The POP in this case is 4.9 s, and the POP start time is 12.5 s. A clear change is observed between the results of the four different physical moments, which are obtained at 5, 10, 15 and 20 s. Due to an insufficient calculation time, the bubbly flow rises vertically, as shown in Figure 3a. An asymmetric gas fraction is observed in the bubble plume. Transient oscillation occurs and is implicitly contained in the gas fraction contour, as shown in Figure 3b. Thereafter, POP occurs, and a periodic state is attained. Two vortices emerge on the two sides of the bubbly flow. The momentum inequality between the two vortices results in the development of an asymmetric flow. Compared to the results of adjacent periods, more bubble plume oscillation details are similar, as shown in Figure 3c,d. Note that the gas fraction changes over time and space. Through the observation of the gas fraction, it is found that the bubbly flow spirals to the outlet of the top domain. Moreover, the gas fraction distribution is not fully developed before the onset of POP. The bubbly flow shape varies since the bubbles flowing in and out of the domain are not in dynamic equilibrium. Some bubbles become separated from the main bubbly flow zone. A steady bubble plume oscillation phenomenon can be observed in the bubble column after 12.5 s, and the gas fraction distribution continuously changes during this period.

In a word, from the instantaneous gas fraction of bubbles, we can clearly see the oscillation and offset characteristics of bubbles, especially in the period of oscillation.

In the first 10 s, the offset of bubble plume is very weak, and the gas–liquid mixing is uniform, but as the plume oscillation period begins, the offset of bubble increases, resulting in the flow field disorder and the gas–liquid mixture is inhomogeneous, but the scope and intensity of momentum exchange increase constantly.

In the test case of ξ3R3, contours of the gas fraction are generated at five sections, as shown in Figure 4, which are located at different Z-coordinates. On both sides of the sections (Z = −0.02, 0.02), the simulated contours are not notable. This effect becomes more pronounced at higher positions, where notable deviations are observed from the central plane. At the other sections (Z = −0.01, 0, 0.01), the contours in this case agree very well, except for a small difference at the sparger. With increasing distance from the central plane, the gas fraction notably decreases. However, the profile simultaneously changes from a full bubbly flow to a partial bubbly flow. The symmetrical profiles always remain the same. The higher the bubbly flow position is, the smaller the gas volume fraction is. This illustrates that the central surface (Z = 0) suitably represents the column to analyze bubble plume oscillation.

In the test case of ξ3R3, the POP was acquired by averaging the time interval of the X-component of the water velocity at a certain point in the center plane at a height of y = 0.25 m. The trajectory of the maximum offset position is shown in Figure 5 at heights of Y = 0.32 m and Z = 0 since the maximum offset can be determined in terms of the height from the two-dimensional perspective of the central surface during POP. In addition, we also discover that the Z-component of the locations of the maximum gas fraction always occurs at the center position (Z = 0) at different times, except for the moments at 2 and 11 s (the Z-components are −0.0016 and −0.0012, respectively). This may occur due to bubble plume instability, but most of the moments are observed in the center. A spline curve is applied to connect all the moments sequentially, which only reveals the trend.

Before the onset of POP, an unstable fluctuation of the maximum gas fraction is observed since the fluid fields are still developing. Thereafter, the locations exhibit regular changes, and the agreement is good for each POP, especially in regard to the peaks. Moreover, it is suggested that this happens due to the balance between bubble coalescence and breakup. This explanation matches the appearance of the bubble plume.

We can find that the period of bubble oscillation is synchronous with the offset characteristic. This is beneficial for understanding the offset time, but the time at which the maximum offset position is located is not necessarily the start or end of the oscillation period. This provides a basis for determining the location and time of mass transfer and gas–liquid mixing.

4.3. Analysis of the Vorticity Distribution and Velocity

Similarly, we choose the test case of ξ3R3 as an example for analysis. The vorticity indicates the velocity and orientation of local rotation, which is adopted to represent the vortex characteristics. To determine the influence of bubbly flow in the column, it is of significance to investigate the vortex intensity and vorticity distribution. Figure 6 shows the distribution of various vortices of different magnitudes in the column at 15 s under the condition of ξ = 3. It is clear that the number and intensity of the vortices in the Y-component are very low, and we thus ignore them here. As shown in Figure 6a, consistent with the distribution of the liquid velocity field, a low-vorticity magnitude area (≤5 s⁻¹) is located in the bubbly flow in the column. High-vorticity magnitude areas (>5 s⁻¹) are located near the corner, top surface and maximum oscillation positions. The results also indicate that the vorticity distribution of the X-component is highly dispersed along the column, which corresponds to long-distance two-phase flow in the X-direction. Figure 6b shows the vorticity magnitude distribution of the Z-component. Larger-scale vortex structures are captured in the column. The results indicate that the position in the Z-direction and momentum exchange are limited, but consistent with the gas fraction distribution, more intense vortices are clearly encountered. High-vorticity magnitude areas occur near the walls and at the inlet.

From the distribution of vortex intensity in directions, the greater the vortex intensity, the more obvious the mass transfer. The best mixing is around the maximum offset position, the second is above the reactor, and the worse is above the maximum offset position. Therefore, in order to better design the bubble column device, the device height can be reduced effectively.

Figure 7 shows the velocity vector of a particular three-dimensional simulation. The dynamic movement of the bubble plume is similar to the experiments. Vortices in the opposite direction are also shown surrounding the bubbly flow. The result of the section clearly captures bubble oscillation, and the velocity magnitudes are mainly concentrated in the middle of the domain. This is explained by the turbulent viscosity of the continuous phase, which decreases the bubble movement, and the size of one vortex is slightly larger than that of the other.

Vectors of velocity and vortex distribution are interrelated as well. In Figure 7, we find that the gas–liquid mixing at the vortex on both sides of the bubble flow is very good, and the liquid phase is surrounded by the gas phase in the middle area, we get the gas–liquid momentum exchange and mass transfer concentrated in this region.

4.4. Effect of the Gas Volume Flux on the Offset Characteristics of the Plume

The offset characteristics of bubbly flow are the main feature of the bubble column, and are captured by acquiring model predictions as a function of η and θ. Moreover, the central surface (Z = 0) is applied to analyze the offset characteristics since it effectively represents the whole domain and the column depth is very small. The characteristics are approximately the same in each POP, and the 20th period is consequently chosen for unified analysis. It is worth noting that there is no oscillation period in the case of a high gas flux and aspect ratio, but the first oscillation of bubbly flow reaches the left or right wall of the column within a short time, and the height of the knee point of the first oscillation tends to remain unchanged. Hence, it has no impact on our data. An interesting phenomenon occurs at a dimensionless maximum oscillation distance of 1, which is the condition in this situation, since the column width is relatively small, and the gas flux is sufficiently high that the bubbly flow cannot be fully expanded, resulting in the bubbly flow not exhibiting a POP.

To investigate the effect of the gas volume flux on the bubbly flow offset behavior in columns of different aspect ratios (ξ), test simulations were performed at different aspect ratios and six gas fluxes: 136 LPH (lowest flow rate), 226 LPH (relatively low flow rate), 317 LPH (moderately flow rate), 407 LPH (relatively moderate flow rate), 497 LPH (relatively high flow rate), and 588 LPH (highest flow rate).

Figure 8 shows the dimensionless maximum offset distance at the various gas fluxes, and Figure 9 shows the angles under the corresponding conditions. At the low fluxes, the bubble plume exiting the sparger follows a specific path, and bubbly flow is developed. Depending on the initial condition, this bubbly flow can repeatedly change course within a certain range. However, at the high fluxes, bubble oscillation intensifies, and the domain becomes disordered. It is also observed that the area occupied by the bubble plume rapidly increases until bubbles flow along the column wall. Therefore, a typical no-POP mode of bubbly flow is observed at the high fluxes and aspect ratios. Figure 9 indicates that the liquid surrounding the bubble plume remains relatively stable at the low fluxes, and an increasing trend of the angle is captured. However, at the high fluxes, large angles are observed, and bubbles flow upward close to the containing wall, even along the column.

4.5. Effect of the Aspect Ratio on the Offset Characteristics of Plume

Regarding the effect of different aspect ratios, simulations were performed at aspect ratios ranging from 2 to 4.5. Figure 10 shows the dimensionless maximum offset distance of bubbly flow under these conditions. Under each condition, η is plotted at 6 different levels. The results indicate that at a low aspect ratio, η gradually increases, which essentially confirms that bubbly flow occurs along a single path and that its trajectory does not touch the two walls. One can determine the occurrence of POP by the plume distance. At the upper level, the distance continues to increase and even remains unchanged under some conditions. This happens because bubbly flow is limited in the domain by the column width and cannot fully develop, which implies that the momentum exchange between gas and liquid occurs forcefully. Figure 11 shows the plume angle of bubbly flow at the different aspect ratios. Here, the results indicate that at the low levels, the angle decreases. Two main vortices are formed around the bubbly flow in which liquid flows in the opposite direction. It is also found that the bubbly flow progressively moves toward the column walls with increasing aspect ratio. This is observed because while the two vortices are formed, their sizes are not uniform. At the high levels, it can be inferred that some bubble plumes continuously move toward the walls and exhibit POP. Moreover, other bubble plumes move along the walls with no POP, and the liquid flows violently.

4.6. Correlations

As already mentioned, correlations for the offset characteristics are derived based on the data in this work and under turbulent flow conditions. In Figure 8, Figure 9, Figure 10 and Figure 11, the predicted lines are presented as fitting curves. At low aspect ratios (ξ = 2, 2.5, 3), η is linearly fitted with the flux, as shown in Figure 8. The minimum adjusted R-square value is 0.93, and the value for the predicted red line is 0.95. At high aspect ratios (ξ = 4, 4.5), η and the flux are suitably fitted with a lognormal curve, and the minimum adjusted R-square value is 0.96, while the correlation coefficient of the purple line is 0.99. At an aspect ratio of 3, the change in data is irregular in the transition stage, and there is no suitable fitting relationship. However, the fitting correlation between the flux and angle is uniform, both of which exhibit exponential growth. The minimum correlation coefficient is 0.85, and the red line in Figure 9 has a value of 0.99. As shown in Figure 10, the lognormal profile is adapted to the correlation of the aspect ratio and η. The correlation coefficient of the red line is as high as 0.99, and the lowest value is not smaller than 0.95. Finally, at the highest flux (R = 588 LPH), the aspect ratio and angle reveal an exponential decay relationship (adjusted R-square = 0.94). Under the other conditions, the angle change is very slight, and a suitable fitting curve cannot be established.

4.7. Plume Oscillation Period

Comparing the previous data in terms of the POP at the different values of the aspect ratio and gas volume flux, the simulation results obtained agree well with the predicted data of [16]. It was previously noted [22] that no POP occurred at a high flux and high aspect ratio. This conclusion is confirmed by the present simulations, with more detailed results provided. Figure 12 reveals that in some cases, plume oscillation is not observed. In other words, this means that at a high aspect ratio of 4 or 4.5, the bubbly flow does not oscillate at these high fluxes, and according to a particular flow path along the walls, the bubble plume is almost completely governed by the narrow width (and cannot fully develop). POP only happens under certain conditions.

5. Conclusions

To sum up, we use numerical simulation method to study the offset characteristics of a bubble plume in a bubble column, in which momentum exchange induced by gas–liquid interaction leads to bubble offset. The detailed flow field shows that the characteristics of bubble offset need time to fully develop. It is similar to the period of bubble oscillation, and also presents periodic range changes. In addition, from the height surface (Y-coordinates unchanged) of the maximum offset position, the offset characteristic exhibits periodic fluctuations in each section. The vortices are mainly concentrated around the bubble flow and on the wall in the domain. Furthermore, we apply the maximum offset position to construct the quantitative analysis of the offset characteristics of the bubble plume, and successfully draw the correlations between the dimensionless offset distance and angle in the oscillation characteristics of different flux and aspect ratio. The fitting degree is good as well. With the increase of volume flux, due to the limited area space, the dimensionless offset distance has an upper limit, so that the value will not change after 1, and the offset angle gradually increases; with the increase of aspect ratio, the dimensionless offset distance shows the same trend, but the angle presents a decreasing trend when the volume flow is large, and other situations fluctuate greatly without obvious regularity. Considering the effect of a limited closed region on the offset of bubble plume, the conditions without oscillation period are found as well.

A series of closure models are applied to the case of bubbly flow, especially to examine bubble plume oscillation based on previous works, with a focus on the dynamic bubble plume process. The simulation results agree quite well with the experimental data, including the velocity profile and POP. Some of the observed deviations could occur because a fixed bubble diameter was used in the simulations, and the change in bubble size was neglected. The models can be further optimized, because the resulting effects are not yet truly reliable. The dynamic details of the bubble plume, including the velocity, volume fraction, vortex intensity and dynamic bubbly flow variations, are observed in detail. It is found that the column is not a real three-dimensional model because its depth is very small. Therefore, the intermediate plane (Z = 0) is chosen to describe the dynamic behavior of bubble oscillation. Furthermore, the influence of the aspect ratio and gas volume flux on the oscillation characteristics is examined based on our simulation results, and qualitative fitting correlations are provided. It is noteworthy that when the dimensionless offset distance (η) is 1, no POP is observed. In other words, POP does not occur in bubble columns at a high gas volume flux and aspect ratio. A more detailed investigation of the specific conditions resulting in the absence of POP has a great potential. Additional research on bubble plume oscillation is required for the development of accurate bubble interaction and turbulence models. The gas–liquid mixing and mass transfer can be enhanced around the maximum offset position. Gas–liquid mixing is mainly concentrated in the middle region. Furthermore, the maximum offset position is related to the aspect ratio and volume flux. These should have a wide application range in the chemical industry.