Mesoregime-Oriented Investigation of Flow ... - ACS Publications

Jul 17, 2019 - School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China. ‡. State Key ..... Depending on the structu...
0 downloads 0 Views 2MB Size
Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Mesoregime-Oriented Investigation of Flow Regime Transition in Bubble Columns Chao Han†,‡ and Jianhua Chen*,‡ †

School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China

Downloaded via GUILFORD COLG on July 25, 2019 at 19:35:07 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



ABSTRACT: This work investigates regime transition in bubble columns through a mesoscale framework. The dominant mechanisms in bubble columns are analyzed according to the concept of mesoregime. Model predictions with the modified stability condition agree well with experimental data. It indicates a scenario of the regime transition. That is, the homogeneous regime is liquiddominated, which is characterized by the dominance of liquid over the movement of small bubbles, while the heterogeneous regime is gas-dominated, which is characterized by the dominance of large bubbles over the liquid movement. Consequently, the transition regime is a mesoregime depending on the principle of compromise in competition. The holdup profiles of the three regimes form an enclosed structure, which is corroborated by experiments in the literature. Moreover, the influence of superficial liquid velocities on the regime transition can be qualitatively captured by the model. Also, a GPU solver is developed to cope with the traversal optimizing procedure.

1. INTRODUCTION Flow regime transition in bubble columns has been widely researched due to its significant impacts on reactor performance.1,2 Generally, three regimes of homogeneous, transition, heterogeneous (or churn-turbulent) flow may appear in bubble columns with different operating conditions.3 If small diameter ( 59.3) n

15

ij μ yz 4 Eo·M −0.149jjjj l zzzz j μref z 3 k {

−0.14

Model parameter

H=

Morton number

M=

Eotos number

Eo =

Drag coefficient for single particle

Tomiyama drag model (pure water)

16

μ l4 g(ρl − ρg )

17

ρ l2 σ 3 g (ρl − ρg )db2

l 24 o 0.687 o ) (Re < 1000) o o Re (1 + 0.15Re C D0,p = m o o o o o (Re ≥ 1000) n 0.44 ÄÅ ÅÄ É ÅÅ ÅÅ 16 48 ÑÑÑÑ 8 Eo C D0,b = maxÅÅÅÅminÅÅÅ (1 + 0.15Re 0.687), Ñ, Å ÅÅÇ ÅÇ Re Re ÑÑÑÖ 3 Eo +

18

σ

19

ÉÑ ÑÑ ÑÑ Ñ 4 ÑÑÑÖ

20

Index “i” refers to “S” and “L” for the small bubble phase and the large bubble phase, respectively.

a

and Tomiyama56 (for simplicity, only pure water case in the later model), and the related model equations are listed in Table 2, notice that the correlation of eq 12 is employed to modify the drag coefficient for a swarm of bubbles.34 It is obvious that for larger bubbles (above several millimeters) which are normally seen in bubble columns, the drag coefficients are higher than that of the particles. Thus, Nsurf is calculated according to this difference (refer to eq 7). In addition, Nsurf is assumed to be zero if the bubble drag coefficient is lower. It should be noticed that this ad hoc approach for calculating Nsurf is a temporary expedient and leaves a weakness in the DBS model. A hopeful way is to compute the bubble oscillating behavior and all the energy consumption terms by fully resolved numerical simulation.57,58 It is noteworthy that a critical bubble diameter of dcrit exists in the Grace model, denoting a minimum CD which corresponds to H = 59.3. Similarly, a local minimum CD also exists in the Tomiyama model. In particular, when the Grace model is used, positive Nsurf not only occurs in the right branch of db > dcrit but also in a short-range of the left branch. However, when the Tomiyama model is used in the DBS model, positive Nsurf only appears in part of the right branch, which accounts for its incapability of predicting the jump change, and more details can refer to Chen et al.37 The calculation of Nbreak harnesses the theory of isotropic turbulence applied in bubble columns.59 It is noticed that Nbreak depends on Nturb (eqs 8−11), which is consistent with the consensus of strong coupling between the turbulence and bubble coalescence and breakup.46,60,61 So far, all the three energy terms of Nsurf, Nturb, and Nbreak are calculable depending on the six structural parameters. If these parameters are known, Nsurf can be directly calculated. Consequently, regarding the constant NT, the undetermined Nturb and Nbreak can be iteratively solved according to eq 5.

Though the energy consumption terms are modeled by the aforementioned strategy, the structural parameters are still open variables yet to be closed. Thus, the stability condition (eq 4) is introduced based on the energy resolution. In this sense, the extrema of Nsurf = min and Nturb = min behave as two separate mechanisms representing liquid-dominated and gasdominated, respectively.37 Thus, the stability condition is proposed as the sum of Nsurf and Nturb tending minimum driven by the principle of compromise in competition. Among all the feasible solutions, only the group of structual parameters which lead to the global minimum of the stability criterion (eq 4) is accepted as the real solution. In other words, different from the work of Krishna et al.,62 which takes the small bubble diameter as an input paramter, the DBS model determines both bubble classes by an optimizing procedure. 2.2. GPU Solver for DBS Model. A traversal numerical method was used to solve the DBS model in Yang et al.34 and subsequent works. Three variables of dS, dL, and Ug,S were traversed within a suitable range, leading to a grid with 55 × 55 × 55 nodes. That is, both dS and dL were divided into 55 logarithmic intervals within the acceptable bubble size range of [dmin, dmax], where dmin corresponds to the bubble diameter at the condition of H = 2 according the drag model (refer to eqs 15 and 16), and dmax is the upper bound of the bubble diameter which is arbitrarily specified as 60 mm. While Ug,S was divided into 55 linear intervals within the range of [0, Ug]. Notice that the number of nodes is selected as a rule of thumb due to the limitation of computing load on CPU. Hence a three-dimensional grid of structural parameters was built to initiate the ergodic global search algorithm, with known values of dS, dL, and Ug,S on each node. The conservation equations were solved sequentially to obtain the stability condition value on each node, and then, the minimal Nst was obtained to determine the real state at given operating conditions. The D

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research solving scheme is shown in Figure 2. However, the algorithm was very time consuming due to the three-dimensional

Figure 3. Sketch of GPU solver for DBS model.

Figure 4. Comparison between CPU and GPU calculations (air− water system: ρl = 998 kg/m3; ρg = 1.2 kg/m3; μ = 0.001 Pa s; σ = 0.072 N/m. The experimental data values are from Camarasa et al.3 (Reprinted with permission from Camarasa, E.; Vial, C.; Poncin, S.; Wild, G.; Midoux, N.; Bouillard, J. Influence of coalescence behavior of the liquid and of gas sparging on hydrodynamics and bubble characteristics in a bubble column. Chem. Eng. Process. 1999, 38, 329. Copyright 1999, Elsevier).

Figure 2. Model solving scheme.

optimizing procedure. In order to improve the computing speed, Wang et al.63 developed a simulated annealing (SA) algorithm. Nonetheless, the SA algorithm was found too sensitive to the initial parameters and should be adjusted for different operating conditions. Since the algorithm is easily paralleled, the present work resorts to the GPU acceleration to improve the solver. In the current work, the CUDA 6.0 platform is utilized to perform the algorithm. The solution space is resolved into 64 × 64 × 64 nodes, and the computing threads are conducted as 8 × 8 in the thread blocks and 64 × 64 in the thread grids. By virtue of the intense arithmetic and logic units on the GPU, all the nodes can be computed simultaneously. The coordinates of each node are sent from CPU to GPU, and the conservation equations are solved by the GPU kernel function. Finally, the Nst values are collected on the CPU, and the minimum of Nst is easily obtained, as well as the real solution via the optimized variable set. This procedure is illustrated in Figure 3. The calculation on a NVIDIA Tesla C2050 GPU demands 3−6 s in comparison to 12−15 min on an Intel Xeon E5520 CPU core. Nearly 100 folds of acceleration are achieved without reducing accuracy. The results calculated by CPU and GPU solvers agree well, as seen in Figure 4.

et al.,9 as seen in Figure 5. Inspired by these experimental clues, a possible holdup structure is sketched in Figure 6. The idea of the enclosed structure in Figure 6 is also underpinned by the practice in gas−solid systems.44 That is, the fluiddominated regime leads to an upper voidage profile, whereas the particle-dominated regime leads to a lower voidage profile,

Figure 5. Typical gas holdup profiles in experiments: (a) Zahradnik et al.9 (Reprinted with permission from Zahradnik, J.; Fialova, M.; Ruzicka, M.; Drahos, J.; Kastanek, F.; Thomas, N. H. Duality of the gas−liquid flow regimes in bubble column reactors. Chem. Eng. Sci. 1997, 52, 3811. Copyright 1997, Elsevier). (b) Camarasa et al.3 (Reprinted with permission from Camarasa, E.; Vial, C.; Poncin, S.; Wild, G.; Midoux, N.; Bouillard, J. Influence of coalescence behavior of the liquid and of gas sparging on hydrodynamics and bubble characteristics in a bubble column. Chem. Eng. Process. 1999, 38, 329. Copyright 1999, Elsevier).

3. EXTREMUM TENDENCIES AND DOMINANT MECHANISMS Experimentally, different holdup profiles come out depending on the uniformity of gas distribution. Typical results of an air− water system can be found in Camarasa et al.3 and Zahradnik E

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

to be lower afterward. Nturb = min leads to a monotonously increasing holdup profile, which deviates from the expected lowest one. Moreover, in the case of Nturb = min, the optimized bubble diameter always locates on the upper bound of dmax. Since the drag coefficient declines with reducing bubble diameters in the range of db > dcrit, decreasing dmax can lead to lower holdup. Though further adjustment of dmax can provide a plausible result for Nturb = min, it is not the objective of the current work and hence avoided. To analyze the aforementioned deviation, first the energy consumption rate terms are surveyed, as shown in Figure 8.

Figure 6. Conceived sketch of regime-specific gas holdup profiles.

accordingly the particle-fluid compromising regime shows an intermediate voidage profile which bridges the former two cases. Detailed meaning of the dominant mechanisms in Figure 6 is further illustrated. In order to elicit the dominant mechanisms, investigating various extrema is a good starting point. As stated in the EMMS model for gas−solid flow, the original stability condition of Nst = min was derived from the compromise in competition (CIC) between the two dominant mechanisms of Wst = min and ε = min. The fluid-dominated mechanism of Wst = min leads to a dilute phase, while the particle-dominated mechanism of ε = min leads to a dense phase. The mesoregime of Nst = min shows a voidage profile enclosed by those of Wst = min and ε = min.44 It is expected that the dominant mechanisms of Nsurf = min and Nturb = min play likewise in the gas−liquid flow. That is, Nsurf = min corresponds to a small bubble state, while Nturb = min corresponds to a large bubble state. Furthermore, small bubbles indicate the flow is liquiddominated by means of viscosity and surface tension with relative stable bubble interface, while large bubbles indicate the flow is gas-dominated by means of buoyancy. Therefore, it is expected that Nsurf = min is liquid-dominated with small bubbles, whereas Nturb = min is gas-dominated with large bubbles, and in between is the so-called mesoregime.64 Following this idea, the regime specific behaviors are investigated, and the modified model can be solved efficiently owing to the improved GPU solver. Figure 7 shows the calculated total gas holdup in an air− water system under different stability conditions. Taking the original solution of Nst = min as the benchmark, Nsurf = min leads to a higher holdup before Ug = 0.25 m/s, but it turns out

Figure 8. Dimensionless energy consumption rates.

For all the three cases, the ratio of Nbreak keeps fairly stable. In the case of Nst = min, the ratio of Nsurf and Nturb alternate their preponderance with increasing Ug. On the contrary, the ratio of Nsurf and Nturb always take smaller values under their minimization conditions, respectively. Furthermore, the energy consumption rates in all cases are continuous, which provide little information to analyze the complex shifts in the holdup profiles. Furthermore, detailed structural parameters are shown in Figure 9. As shown in Figure 9a and b, in the case of Nturb = min, dS and dL are always identical to dmax, so that the component holdup of f S and f L cannot be distinguished. In the case of Nsurf = min, dL declines to dcrit before Ug = 0.04 m/s, while dS increases to a plateau until Ug = 0.26 m/s and then jumps to dcrit. That means beyond Ug = 0.26 m/s, dS and dL are also identical for Nsurf = min. For the case of Nst = min, as reported in previous publications of the DBS model, a jump change of dS to dcrit happens, which coincides with the holdup change. The large bubble diameter declines until it arrives at dcrit at very high gas velocities. Figure 9c and d shows f S and f L for the cases of Nst = min and Nsurf = min, respectively. Since it is meaningless to distinguish f S and f L when dS and dL are identical, the results of Nturb = min are omitted, likewise for Nsurf = min beyond Ug = 0.26 m/s and Nst = min beyond Ug = 0.48 m/s. In the case of Nst = min, f S and f L are comparative before the sharp decrease in f S. Then, f S increases with a steep slope, while f L approaches a maximum around Ug = 0.25 m/s and declines until dL = dcrit. The results indicate that during the compromise in competition with increasing Ug, neither bubble phase can prevail before the jump change. However, after the jump change, the small bubble phase taking the form of dS = dcrit is going to prevail, even though smaller f S is found at the beginning. In the case of Nsurf = min, there is little large bubble phase before Ug = 0.1 m/s. Afterward, f S declines and f L

Figure 7. Total gas holdup under different stability conditions. F

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 9. Structural parameters under different stability conditions.

increases, noticing dS < dcrit and dL = dcrit in the gas velocity range of 0.1−0.26 m/s. In short, the results shown in Figure 9 imply that at high gas velocities the bubble diameter tends to be dcrit. It invokes a suspect that db = dcrit may be a substitute for the mechanism of Nturb = min. Obviously, the lowest drag coefficient means the bubbles pass through the liquid bed as soon as possible, thus leading to a minimal gas holdup; in other words, db = dcrit is equivalent to f = min. Thus, if Nturb = min is replaced by f = min as a dominant mechanism and Nsurf = min as another, the stability condition should be reconstructed to express the compromise in competition, which is elucidated in the next section.

Figure 10. Three regimes with different stability conditions: Nsurf = min, f = min, and Nsurf × f = min.

4. MESOREGIME-ORIENTED RESEARCH As a first step, the product of Nsurf × f = min is tested, and with this modification, the calculation of Nturb and Nbreak is dispensable and hence can be omitted. In addition, several remarks should be made here. First, different from Nst = min in the gas−solid flow, where Nst refers to the energy consumption rate for suspending and transporting with respect to unit mass of particles, thus far there is no explicit physical meaning for Nsurf × f in the gas−liquid flow. In fact, the multiplication between Nsurf and f is a feasible combination of the multiobjective optimization components. Second, the bubble phase with dcrit in diameter should be regarded as an equivalence of the “large” bubble phase in terms of high rising velocity rather than the apparent “small” bubble diameter. 4.1. Results and Discussions. Figure 10 shows the calculated gas holdup by different stability conditions. Interestingly, the holdup profiles form an enclosed structure with the upper bound of Nsurf = min, the lower bound of f = min, and in between of Nsurf × f = min. The holdup profile of f = min increases monotonously, while the holdup profile of

Nsurf = min has a weak shift around Ug = 0.1 m/s and merges with the lower bound at Ug = 0.26 m/s. These changes are attributed to the bubble diameter variations, which are clearly seen in Figure 11. For the case of Nsurf × f = min, dS jumps to dcrit at Ug = 0.11 m/s in comparison to Ug = 0.26 m/s for the case of Nsurf = min. Moreover, the case of f = min always leads to a homogeneous bubble phase with dS = dL = dcrit. In general, Nsurf × f = min inherits the feature of Nsurf = min, but it is also subjected to the constraint of f = min which precipitates the transition of the small bubble diameter to the stable value of dcrit. Figure 12 compares the calculated gas holdup profiles and experimental data extracted from Figure 5. Despite the quantitative deviations, the holdup profiles predicted by f = min coincide with the experimental results by nonuniform distribution, while Nsurf × f = min reproduces holdup profiles similar to the experiments of perforated spargers with multiple orifices, featuring a hump with a maximum around Ug = 0.08 G

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

aeration and small disturbance exerted by bubble swarms, so that the macroscale flow is suppressed by the liquid viscosity. Therefore, the homogeneous regime is only observed in the cases of uniform spargers and very low gas velocities, which indicates the dominance of the liquid phase over the small bubbles. Owing to the identification of dominant mechanisms, the mesoregime can be defined via Nsurf × f = min, which roots in the CIC principle. The mesoregime undergoes a shift from the liquid-dominated regime to the gas-dominated regime with increasing gas velocities. The hump in the holdup curve is qualitatively captured by the mesoregime, as expected in the sketch of Figure 6. Furthermore, the effect of liquid velocity on the regime transition is investigated by using this modified model. As is well known, the influence of the liquid velocity on the gas holdup depends on whether it is comparable to the gas velocity.10,20 Low liquid velocities have negligible effects on the gas holdup because the bubble acceleration is weak. However, the gas holdup will be affected remarkably if Ul is high. In general, positive Ul (co-current) reduces the gas holdup, and negative Ul (counter-current) increases the gas holdup. This phenomenon is due to the fact that the liquid motion will accelerate or suppress the bubble rising in the co-current and counter-current modes, respectively. Figure 13 compares the experimental data reported by Besagni and Inzoli20 and the model prediction by this work. In

Figure 11. Calculated bubble diameters by different stability conditions.

Figure 12. Comparison between the calculated and experimental holdup profiles: (a) Zahradnik et al.9 (Reprinted with permission from Zahradnik, J.; Fialova, M.; Ruzicka, M.; Drahos, J.; Kastanek, F.; Thomas, N. H. Duality of the gas−liquid flow regimes in bubble column reactors. Chem. Eng. Sci. 1997, 52, 3811. Copyright 1997, Elsevier). (b) Camarasa et al.3 (Reprinted with permission from Camarasa, E.; Vial, C.; Poncin, S.; Wild, G.; Midoux, N.; Bouillard, J. Influence of coalescence behavior of the liquid and of gas sparging on hydrodynamics and bubble characteristics in a bubble column. Chem. Eng. Process. 1999, 38, 329. Copyright 1999, Elsevier).

Figure 13. Comparison between the experimental data in literature and the model prediction. The experimental data values are from Besagni and Inzoli20 (Reprinted with permission from Besagni, G.; Inzoli, F. Influence of internals on counter-current bubble column hydrodynamics: Holdup, flow regime transition and local flow properties. Chem. Eng. Sci. 2016, 145, 162. Copyright 2016, Elsevier).

m/s. Moreover, the end of the hump agrees well with the second regime transition point. The mechanism of f = min can be understood as gasdominated because the lower gas holdup means a shorter residence time for the gas phase, which is equivalent to the flow system with very large bubbles. In this regime, bubbles are dominated by the buoyant force rather than the interaction exerted by the liquid phase, so that gas passes through the liquid bed with the least resistance. In addition, the merged holdup profiles at high gas velocities confirm the gas-prevailing characteristics of f = min. On the other hand, concerning the nonuniform aeration experimental systems, the concentrated gas stream tends to induce a macroscale circulation in the liquid bulk phase, and consequently, gas rises more easily via the concurrent liquid flow in the central region. Therefore, f = min and nonuniform aeration are both considered to generate a gas-dominated regime. In contrast, the mechanism of Nsurf = min is liquiddominated because it is prone to give small bubbles until Ug = 0.26 m/s. Notice that though dL is larger than dcrit at lower Ug (Figure 11a), f L is always zero before dL arrives at dcrit as seen in Figure 9d, and it has no contribution to the total holdup. In practice, the liquid-dominated regime demands uniform

the experiment of counter-current flow (Figure 13a), with increasing magnitude of liquid velocities (|Ul|), the holdup augments at lower gas velocities, and the transition point from homogeneous to transition regime occurs at lower Ug. The reason was attributed to the compact effects, viz., the downward liquid slows the rise of bubbles and leads to higher holdup. Consequently, the bubbles arrange in a more compact mode which makes an earlier transition. On the other hand, the model predicts similar trends on several aspects, as shown in Figure 13b. First, the holdup profile varies with different |Ul|, showing mutative location and shape of plateau, and the holdup is higher for larger |Ul| at lower gas velocities, e.g., before arriving at the plateau. Second, the calculated plateaus shift to lower Ug with increasing |Ul|. Third, the holdup is positively correlated with |Ul| beyond the plateau, viz., in the gas-dominated regime. It should be noticed that the deviation between the calculation and the experimental data is mainly attributed to the sparger effect. Since the spider type sparger was used in Figure 13a, lower holdup profiles were obtained H

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

According to this observation, we argue that the homogeneous regime is liquid-dominated, which only happens at lower gas velocities combined with a uniform aeration. In contrast, the heterogeneous regime is gas-dominated, which happens in nonuniform aeration or at higher gas velocities in uniform aeration. Consequently, the transition regime is a result of a compromise in competition between the moving tendencies of gas and liquid. As implied by the effects of p, the real bubbly flow depends on various factors, including the aeration uniformity, liquid properties, column geometries, and so forth. Therefore, the transition regime is regarded as a flexible region. Indeed, neither of the dominant mechanisms can exclusively realize their tendencies in the compromising region. Even in the so-called liquid-dominated or gasdominated regimes, the subordinate mechanism may still exert some influence on the flow field. According to this framework, gradual change is more reasonable than abrupt change to characterize the regime transition in bubble columns, albeit the abrupt change predicted by the original DBS model reveals the intrinsic nonlinearity of the gas−liquid flow. As a byproduct, the index of p reconciles the gradual transition mode in the current work and the abrupt mode in the DBS model. The results by varying p prove the robustness of the DBS/EMMS framework. In this sense, it is just an observable appearance of the abrupt or gradual change in the holdup profiles, yet the underlying CIC principle is more essential for understanding regime transition in bubble columns. There is no doubt that the present work should not be overstated for tackling such a complex problem. Though both abrupt and gradual changes can be predicted by adjusting the stability condition, there are still some shortcomings in this modification. First, the liquid turbulence is not adequately modeled, which leads to the distortion of the predicted bubble diameter; that is, the bubble diameters of both phases are no greater than dcrit. Indeed, the holdup profile of f = min approximates the fully developed turbulent bubbly flow, which is based on the fact that high rising velocities appear for both very large bubbles and bubbles in diameter of dcrit. Second, Nsurf is calculated in a phenomenological way which depends on the drag models. A mechanical model is highly appreciated for calculating each part of the energy consumption rates, so that the exact stability condition can be verified. Finally, the oversimplification of neglecting Nbreak weakens the linkage between this heuristic model and the widely used CFD methods. All these shortcomings need to be addressed in the future to consolidate the rigor of physics.

for the nonuniform distribution, whereas the model is zerodimensional without concerning of the sparger and column geometry, so that it cannot reflect the sparger effect directly. Essentially, the model results support the viewpoint that counter-current flow will destabilize the homogeneous flow and lead to earlier transitions as reported in the literature.10,20,41 More generally, with increasing algebraic value of Ul, the plateau on the holdup curve shifts to higher Ug and the gas-dominated regime shows lower gas holdup, as shown Figure 14. So far, it is not very clear about the underlying

Figure 14. Model prediction for the effects of the liquid velocity on the gas holdup and regime transition.

mechanisms of these phenomena; a possible explanation may be that negative Ul compacts the bubble swarms whereas positive Ul loosens the swarms. Therefore, the so-called “gasdominated” regime (viz., gas buoyancy dominates the bubble motion) is advanced in the counter-current mode and delayed in the co-current mode. 4.2. Further Discussions. An imminent question is therefore how to reconcile the abrupt change predicted by the original DBS model and the gradual change by the current work. Therefore, an extended criterion of Nsurf × f p = min is proposed, where the index of p refers to the intensity of competition. Both gradual and abrupt holdup profiles are observed with varying values of the power index, as shown in Figure 15. Higher p values lead to earlier transitions to the gas-

5. CONCLUSIONS Flow regime transitions in bubble columns are believed to relate with the mesoscale structures. The DBS model was capable to predict the second transition point by an abrupt change in the gas holdup profile. However, it deviates from the gradual change widely reported in experiments. To this end, this paper performed a mesoregime-oriented study on this issue and revealed some interesting results. The developed GPU solver is proven to be robust and efficient for the DBS model. Since a direct resolution of the original stability condition cannot eliminate the deviation, the expressions of dominant mechanisms were reconsidered according to the structural parameters. Consequently, f = min ( f denotes the gas holdup) was introduced as the gasdominated mechanism, and Nsurf = min was regarded as the

Figure 15. Variation of holdup profiles versus the power index of p.

dominated regime, which means a strengthened effect of the dominant mechanism of f = min. Moreover, with increasing p, the holdup profile evolves from gradual to abrupt, which indicates a complicated correlation between the structural parameters. I

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

p = power index reflecting the intensity of competition, dimensionless Pb = bubble breakup probability, dimensionless Ug,L = superficial gas velocity for large bubbles, m s−1 Ug,S = superficial gas velocity for small bubbles, m s−1 Ug = superficial gas velocity, m s−1 Ul = superficial liquid velocity, m s−1 UT = terminal rise velocity, m s−1 Wst = energy consumption rate with respect to unit volume of bed, W m−3

liquid-dominated mechanism. The modified stability condition of Nsurf × f = min predicted holdup profiles consistent with the experiments in the literature, featuring a hump with increasing gas velocities. The influence of superficial liquid velocities on the regime transition can be qualitatively predicted. Furthermore, by introducing a power index, both gradual and abrupt changes were observed in the holdup profiles. It indicates that the CIC principle deserves in-depth study when extending mesoscale methods to other systems.





AUTHOR INFORMATION

GREEK LETTERS λ = characteristic size of eddy, m μ = liquid dynamic viscosity, Pa·s ω = collision frequency, s−1 ρ = density, kg m−3 σ = surface tension, N m−1 ε = average voidage in gas−solid flow, dimensionless

Corresponding Author

*E-mail: [email protected]. ORCID

Jianhua Chen: 0000-0002-4068-3956 Notes

The authors declare no competing financial interest.





ACKNOWLEDGMENTS This work is financially supported by the State Key Laboratory of Multiphase Complex Systems (MPCS-2015-A-03) and the National Natural Science Foundation of China (91834303). We are grateful to Prof. Jinghai Li for his plentiful guide of the mesoscale methodology. The helpful discussions with our colleagues, Prof. Ning Yang, Dr. Xiaoping Guan, Prof. Wei Ge, and Dr. Wenlai Huang, are also greatly appreciated.

ABBREVIATIONS CIC = compromise in competition CPU = central processing unit CUDA = compute unified device architecture DBS = dual-bubble-size EMMS = energy minimization multiscale GPU = graphics processing unit SA = simulated annealing





NOMENCLATURE C Db = drag coefficient for a bubble in a swarm, dimensionless CD0,b = drag coefficient for a bubble in a quiescent liquid, dimensionless CDp = drag coefficient for a particle in multiparticle systems, dimensionless CD0,p = drag coefficient for a particle in a quiescent fluid, dimensionless cf = coefficient of surface area augment, cf = f BV2/3 + (1 − f BV)2/3 − 1, dimensionless db = bubble diameter, m dcrit = critical bubble diameter in the drag model, m dL = bubble diameter of the large bubble phase, m dmin = minimum bubble diameter, m dS = bubble diameter of the small bubble phase, m Eo = Eotvos number, dimensionless f = gas holdup, f = f S + f L, dimensionless f BV = breakup ratio of daughter bubble to its mother bubble, dimensionless f L = holdup of the large bubble phase, dimensionless f S = holdup of the small bubble phase, dimensionless g = gravitational acceleration, m s−2 Mo = Morton number, dimensionless Nbreak = energy consumption rate due to bubble breakage with respect to unit mass of liquid, W kg−1 nb = bubble number density, m−3 Nst = sum of Nsurf and Nturb in the gas−liquid flow; also the stability condition in the EMMS model, W kg−1 Nsurf = energy consumption rate due to bubble oscillation with respect to unit mass of liquid, W kg−1 Nturb = energy consumption rate in turbulent liquid phase with respect to unit mass of liquid, W kg−1 NT = total energy consumption rate with respect to unit mass of liquid, W kg−1

SUBSCRIPTS b = bubble g = gas phase l = liquid phase L = large bubble p = particle S = small bubble T = total



REFERENCES

(1) Deckwer, W. D. Bubble Column Reactors; Wiley: Chichester, UK, 1992. (2) Shah, Y. T.; Kelkar, B. G.; Godbole, S. P.; Deckwer, W. D. Design parameters estimations for bubble column reactors. AIChE J. 1982, 28, 353. (3) Camarasa, E.; Vial, C.; Poncin, S.; Wild, G.; Midoux, N.; Bouillard, J. Influence of coalescence behaviour of the liquid and of gas sparging on hydrodynamics and bubble characteristics in a bubble column. Chem. Eng. Process. 1999, 38, 329. (4) Besagni, G.; Di Pasquali, A.; Gallazzini, L.; Gottardi, E.; Colombo, L. P. M.; Inzoli, F. The effect of aspect ratio in countercurrent gas-liquid bubble columns: Experimental results and gas holdup correlations. Int. J. Multiphase Flow 2017, 94, 53. (5) Ruzicka, M. C.; Zahradnik, J.; Drahos, J.; Thomas, N. H. Homogeneous−heterogeneous regime transition in bubble columns. Chem. Eng. Sci. 2001, 56, 4609. (6) Sharaf, S.; Zednikova, M.; Ruzicka, M. C.; Azzopardi, B. J. Global and local hydrodynamics of bubble columns − Effect of gas distributor. Chem. Eng. J. 2016, 288, 489. (7) Besagni, G.; Gallazzini, L.; Inzoli, F. Effect of gas sparger design on bubble column hydrodynamics using pure and binary liquid phases. Chem. Eng. Sci. 2018, 176, 116. (8) Gemello, L.; Plais, C.; Augier, F.; Cloupet, A.; Marchisio, D. L. Hydrodynamics and bubble size in bubble columns: Effects of contaminants and spargers. Chem. Eng. Sci. 2018, 184, 93. J

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (9) Zahradnik, J.; Fialova, M.; Ruzicka, M.; Drahos, J.; Kastanek, F.; Thomas, N. H. Duality of the gas-liquid flow regimes in bubble column reactors. Chem. Eng. Sci. 1997, 52, 3811. (10) Besagni, G.; Inzoli, F. Comprehensive experimental investigation of counter-current bubble column hydrodynamics: Holdup, flow regime transition, bubble size distributions and local flow properties. Chem. Eng. Sci. 2016, 146, 259. (11) Mudde, R. F.; Harteveld, W. K.; van den Akker, H. E. A. Uniform Flow in Bubble Columns. Ind. Eng. Chem. Res. 2009, 48, 148. (12) Ruzicka, M. C.; Drahos, J.; Mena, P. C.; Teixeira, J. A. Effect of viscosity on homogeneous−heterogeneous flow regime transition in bubble columns. Chem. Eng. J. 2003, 96, 15. (13) Besagni, G.; Inzoli, F. The effect of liquid phase properties on bubble column fluid dynamics: Gas holdup, flow regime transition, bubble size distributions and shapes, interfacial areas and foaming phenomena. Chem. Eng. Sci. 2017, 170, 270. (14) Besagni, G.; Inzoli, F.; De Guido, G.; Pellegrini, L. A. The dual effect of viscosity on bubble column hydrodynamics. Chem. Eng. Sci. 2017, 158, 509. (15) Ruzicka, M. C.; Vecer, M. M.; Orvalho, S.; Drahos, J. Effect of surfactant on homogeneous regime stability in bubble column. Chem. Eng. Sci. 2008, 63, 951. (16) Ribeiro, C. P., Jr; Mewes, D. The influence of electrolytes on gas hold-up and regime transition in bubble columns. Chem. Eng. Sci. 2007, 62, 4501. (17) Besagni, G.; Inzoli, F. The effect of electrolyte concentration on counter-current gas−liquid bubble column fluid dynamics: Gas holdup, flow regime transition and bubble size distributions. Chem. Eng. Res. Des. 2017, 118, 170. (18) Besagni, G.; Inzoli, F. Bubble size distributions and shapes in annular gap bubble column. Exp. Therm. Fluid Sci. 2016, 74, 27. (19) Wilkinson, P. M.; Spek, A. P.; van Dierendonck, L. L. Design parameters estimation for scale-up of high-pressure bubble columns. AIChE J. 1992, 38, 544. (20) Besagni, G.; Inzoli, F. Influence of internals on counter-current bubble column hydrodynamics: Holdup, flow regime transition and local flow properties. Chem. Eng. Sci. 2016, 145, 162. (21) Shaikh, A.; Al-Dahhan, M. H. A Review on Flow Regime Transition in Bubble Columns. Int. J. Chem. React. Eng. 2007, 5, 1. (22) Majumder, S. K. Hydrodynamics and Transport Processes of Inverse Bubbly Flow; Elsevier: Amsterdam, 2016. (23) Drahos, J.; Zahradnik, J.; Puncochar, M.; Fialova, M.; Bradka, F. Effect of operating conditions on the characteristics of pressure fluctuations in a bubble column. Chem. Eng. Process. 1991, 29, 107. (24) Vial, C.; Poncin, S.; Wild, G.; Midoux, N. Experimental and theoretical analysis of the hydrodynamics in the riser of an external loop airlift reactor. Chem. Eng. Sci. 2002, 57, 4745. (25) Zuber, N.; Findlay, J. A. Average volumetric concentration in two-phase flow systems. J. Heat Transfer 1965, 87, 453. (26) Shnip, A. I.; Kolhatkar, R. V.; Swamy, D.; Joshi, J. B. Criteria for the transition from the homogeneous to the heterogeneous regime in two-dimensional bubble column reactors. Int. J. Multiphase Flow 1992, 18, 705. (27) Joshi, J. B.; Deshpande, N. S.; Dinkar, M.; Phanikumar, D. V., Hydrodynamic stability of multiphase reactors. In Advances in Chemical Engineering; Academic Press: 2001; Vol. 26, pp 1. (28) Lucas, D.; Prasser, H. M.; Manera, A. Influence of the lift force on the stability of a bubble column. Chem. Eng. Sci. 2005, 60, 3609. (29) Monahan, S. M.; Fox, R. O. Linear stability analysis of a twofluid model for air-water bubble columns. Chem. Eng. Sci. 2007, 62, 3159. (30) Wang, T.; Wang, J.; Jin, Y. Theoretical prediction of flow regime transition in bubble columns by the population balance model. Chem. Eng. Sci. 2005, 60, 6199. (31) Zahradnik, J.; Fialova, M. The effect of bubbling regime on gas and liquid phase mixing in bubble column reactors. Chem. Eng. Sci. 1996, 51, 2491. (32) Orvalho, S.; Hashida, M.; Zednikova, M.; Stanovsky, P.; Ruzicka, M. C.; Sasaki, S.; Tomiyama, A. Flow regimes in slurry

bubble column: Effect of column height and particle concentration. Chem. Eng. J. 2018, 351, 799. (33) Olmos, E.; Gentric, C.; Midoux, N. Numerical description of flow regime transitions in bubble column reactors by a multiple gas phase model. Chem. Eng. Sci. 2003, 58, 2113. (34) Yang, N.; Chen, J.; Zhao, H.; Ge, W.; Li, J. Explorations on the multi-scale flow structure and stability condition in bubble columns. Chem. Eng. Sci. 2007, 62, 6978. (35) Li, J.; Tung, Y.; Kwauk, M. Method of Energy Minimization in Multi-Scale Modeling of Particle-Fluid Two-Phase Flow. In Circulating Fluidized Bed Technology II; Basu, P., Large, J. F., Ed.; Pergamon Press: New York, 1988; p 75. (36) Zhao, H. Multi-scale modeling of gas-liquid (slurry) reactors. Ph.D. Thesis, Institute of Process Engineering, Chinese Academy of Sciences, Beijing, 2006. (37) Chen, J.; Yang, N.; Ge, W.; Li, J. Modeling of Regime Transition in Bubble Columns with Stability Condition. Ind. Eng. Chem. Res. 2009, 48, 290. (38) Yang, N.; Chen, J.; Ge, W.; Li, J. A conceptual model for analyzing the stability condition and regime transition in bubble columns. Chem. Eng. Sci. 2010, 65, 517. (39) Chen, J.; Yang, N.; Ge, W.; Li, J. H. Stability-driven Structure Evolution: Exploring the Intrinsic Similarity Between Gas-Solid and Gas-Liquid Systems. Chin. J. Chem. Eng. 2012, 20, 167. (40) Yang, N. Mesoscale Transport Phenomena and Mechanisms in Gas−Liquid Reaction Systems. In Advances in Chemical Engineering; Marin, G. B., Li, J., Eds.; Academic Press, 2015; Vol. 46, Chapter Five, pp 245. (41) Guan, X.; Yang, N. Modeling of co-current and counter-current bubble columns with an extended EMMS approach. Particuology 2019, 44, 126. (42) Li, J.; Ge, W.; Wang, W.; Yang, N.; Liu, X.; Wang, L.; He, X.; Wang, X.; Wang, J.; Kwauk, M. From Multiscale Modeling to MesoScience; Springer, 2013. (43) Li, J.; Huang, W. From multiscale to mesoscience: Addressing mesoscales in mesoregimes of different levels. Annu. Rev. Chem. Biomol. Eng. 2018, 9, 41. (44) Du, M.; Hu, S.; Chen, J.; Liu, X.; Ge, W. Extremum characteristics of energy consumption in fluidization analyzed by using EMMS. Chem. Eng. J. 2018, 342, 386. (45) Li, J.; Kwauk, M. Particle-Fluid Two-Phase FlowThe EnergyMinimization Multi-Scale Method; Metallurgical Industry Press: Beijing, 1994. (46) Radl, S.; Khinast, J. G. Multiphase flow and mixing in dilute bubble swarms. AIChE J. 2010, 56, 2421. (47) Ziegenhein, T.; Lucas, D. Observations on bubble shapes in bubble columns under different flow conditions. Exp. Therm. Fluid Sci. 2017, 85, 248. (48) Liao, Y.; Ma, T.; Krepper, E.; Lucas, D.; Fröhlich, J. Application of a novel model for bubble-induced turbulence to bubbly flows in containers and vertical pipes. Chem. Eng. Sci. 2019, 202, 55. (49) Shi, W.; Yang, X.; Sommerfeld, M.; Yang, J.; Cai, X.; Li, G.; Zong, Y. Modelling of mass transfer for gas-liquid two-phase flow in bubble column reactor with a bubble breakage model considering bubble-induced turbulence. Chem. Eng. J. 2019, 371, 470. (50) Bhavaraju, S. M.; Russell, T. W. F.; Blanch, H. W. The design of gas sparged devices for viscous liquid systems. AIChE J. 1978, 24, 454. (51) Ge, W.; Chen, F.; Gao, J.; Gao, S.; Huang, J.; Liu, X.; Ren, Y.; Sun, Q.; Wang, L.; Wang, W.; Yang, N.; Zhang, J.; Zhao, H.; Zhou, G.; Li, J. Analytical multi-scale method for multi-phase complex systems in process engineering  Bridging reductionism and holism. Chem. Eng. Sci. 2007, 62, 3346. (52) Ashgriz, N.; Movassat, M. Oscillation of Droplets and Bubbles. In Handbook of Atomization and Sprays: Theory and Applications; Ashgriz, N., Ed.; Springer: Boston, MA, 2011; pp 125. (53) McDougald, N. K.; Leal, L. G. Numerical study of the oscillations of a non-spherical bubble in an inviscid, incompressible liquid. Part I: free oscillations from non-equilibrium initial conditions. Int. J. Multiphase Flow 1999, 25, 887. K

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (54) Zhang, Y. Heat transfer across interfaces of oscillating gas bubbles in liquids under acoustic excitation. Int. Commun. Heat Mass Transfer 2013, 43, 1. (55) Grace, J. R.; Wairegi, T.; Nguyen, T. H. Shapes and velocities of single drops and bubbles moving freely through immiscible liquids. Trans. IChemE 1976, 54, 167. (56) Tomiyama, A. Struggle with computational bubble dynamics. Multiphase Sci. Technol. 1998, 10, 369. (57) Tryggvason, G. Virtual motion of real particles. J. Fluid Mech. 2010, 650, 1. (58) Gumulya, M.; Joshi, J. B.; Utikar, R. P.; Evans, G. M.; Pareek, V. Characteristics of energy production and dissipation around a bubble rising in water. Chem. Eng. Sci. 2019, 193, 38. (59) Luo, H.; Svendsen, H. F. Theoretical Model for Drop and Bubble Breakup in Turbulent Dispersions. AIChE J. 1996, 42, 1225. (60) Han, L.; Gong, S.; Li, Y.; Gao, N.; Fu, J.; Luo, H. a.; Liu, Z. Influence of energy spectrum distribution on drop breakage in turbulent flows. Chem. Eng. Sci. 2014, 117, 55. (61) Cui, Z.; Fan, L. S. Turbulence energy distributions in bubbling gas-liquid and gas-liquid-solid flow systems. Chem. Eng. Sci. 2004, 59, 1755. (62) Krishna, R.; Urseanu, M. I.; van Baten, J. M.; Ellenberger, J. Influence of scale on the hydrodynamics of bubble columns operating in the churn-turbulent regime: experiments vs. Eulerian simulations. Chem. Eng. Sci. 1999, 54, 4903. (63) Wang, Y.; Xiao, Q.; Yang, N.; Li, J. In-Depth Exploration of the Dual-Bubble-Size Model for Bubble Columns. Ind. Eng. Chem. Res. 2012, 51, 2077. (64) Li, J. Mesoscale spatiotemporal structures: opportunities from challenges. Natl. Sci. Rev. 2017, 4, 787.

L

DOI: 10.1021/acs.iecr.9b01715 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX