*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−1. 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 agreemen<sup>t</sup> between the experiments and simulations is quite good.

**Figure 2.** Comparison of the experimental and simulated long time-averaged vertical liquid velocity profiles.

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.

**Figure 3.** Instantaneous gas fraction from 0 to 0.1 and liquid velocity fields in the test case of ξ3R3 at (**a**) 5 s, (**b**) 10 s, (**c**) 15 s, and (**d**) 20 s.

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.

**Figure 4.** Contours of the gas hold-up over time (half-period, 15 s) for the five sections (Z = −0.02, −0.01, 0, 0.01 0.02).

**Figure 5.** Temporal evolution of the maximum gas fraction locations at Y = 0.32 m and Z = 0; the colors indicate the time.

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 agreemen<sup>t</sup> 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 o ffset characteristic. This is beneficial for understanding the o ffset time, but the time at which the maximum o ffset 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 di fferent 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<sup>−</sup>1) is located in the bubbly flow in the column. High-vorticity magnitude areas (>5 s<sup>−</sup>1) 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.

**Figure 6.** Distribution of various vortices of di fferent magnitudes in the column: (**a**) X-component, (**b**) Z-component.

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 o ffset position, the second is above the reactor, and the worse is above the maximum o ffset position. Therefore, in order to better design the bubble column device, the device height can be reduced e ffectively.

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.

**Figure 7.** Contour of the velocity vector of the central surface in the column.

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 ge<sup>t</sup> the gas–liquid momentum exchange and mass transfer concentrated in this region.

#### *4.4. E*ff*ect of the Gas Volume Flux on the O*ff*set 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.

**Figure 8.** Dimensionless maximum oscillation distance of bubbly flow at the various gas volume fluxes.

**Figure 9.** Oscillation angle of bubbly flow at the various gas volume fluxes.

#### *4.5. E*ff*ect of the Aspect Ratio on the O*ff*set Characteristics of Plume*

Regarding the e ffect of di fferent aspect ratios, simulations were performed at aspect ratios ranging from 2 to 4.5. Figure 10 shows the dimensionless maximum o ffset distance of bubbly flow under these conditions. Under each condition, η is plotted at 6 di fferent 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 di fferent 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.

**Figure 10.** Dimensionless maximum oscillation distance of bubbly flow at the various aspect ratios.

ξ

θ

η

ξ

**Figure 11.** Oscillation angle of bubbly flow at the various aspect ratios.
