**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]. Di fficulties 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 grea<sup>t</sup> challenge to improve the mass transfer e fficiency [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–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 agreemen<sup>t</sup> 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 o ffset is of grea<sup>t</sup> significance to the further study of gas–liquid mixing and mass transfer behavior, while the plume oscillation and o ffset characteristics correspond to each other. Based on these, the o ffset 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 o ffset 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 o ffset 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, and experimental data for di fferent 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 o ffset characteristics.


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

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 di fferent. 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 o ffset 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 ge<sup>t</sup> 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**

Di fferent 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. The surface tension and material properties of tap water and air at atmospheric pressure and temperature are summarized in Table 3. Numerical details suitable for this work are introduced later.


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


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. Six values are chosen for the aspect ratio (ξi) and gas volume flux (Ri), 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.


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

**Figure 1.** Schematic of the column geometry: (**left**) 3D representation and (**right**) 2D view.

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:

$$
\eta = \frac{2L}{W} \tag{1}
$$

$$\Theta = \operatorname{artan} \frac{L}{h} \tag{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.
