ARTICLE pubs.acs.org/IECR
In-Depth Exploration of the Dual-Bubble-Size Model for Bubble Columns Yuhua Wang,†,‡ Qi Xiao,†,§ Ning Yang,*,† and Jinghai Li† †
State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, P.O. Box 353, Beijing, 100190, People’s Republic of China ‡ Graduate University of Chinese Academy of Sciences, Beijing, 100049, People’s Republic of China § State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, People’s Republic of China ABSTRACT: The Dual-Bubble-Size (DBS) model is an extension of the Energy-Minimization Multi-Scale (EMMS) model, which was originally proposed for gassolid fluidization to gasliquid systems. This model is featured by a stability condition in addition to conservative equations for two bubble classes. The stability condition is mathematically formulated as a minimization tendency of microscale energy dissipation and physically reflects the compromise between different dominant mechanisms. This work attempts to use the Simulated Annealing (SA) method to solve the nonlinear optimization problem of the model. It is found that the SA method could greatly reduce the computational cost while capturing the jump change of gas holdup more accurately for the DBS model. The jump change reflects the transition from homogeneous and transitional regimes to heterogeneous regime, and essentially arises from the inflection of the trajectory of global minimum point within the space of the three structure parameters. The DBS model then is extended to Triple-Bubble-Size (TBS) and Multiple-Bubble-Size (MBS) models by introducing three or more bubble classes. We find that the TBS and MBS model prediction is reduced to that of DBS model and the two characteristic bubble classes are distinct in the calculation, even though the gas is resolved into three or more bubble classes. This implies that gasliquid flow in bubble columns are essentially dominated by two bubble classes rather than multiple bubble classes, and the DBS model may be an intrinsic model for reflecting the compromising mechanisms in the system and thereby describing the system evolution and gasliquid interaction. Its integration with computational fluid dynamics (CFD) simulation may offer a more reasonable way to model the complex gasliquid flow in bubble columns.
1. INTRODUCTION Although bubble columns have found many applications in the process industry, the hydrodynamics in bubble columns is complicated to understand, because of the formation and dynamic variation of multiscale structures. For example, the distribution of gas holdup and bubble size may be pertinent to those phenomena occurring at mesoscales such as bubble wake, bubble coalescence and breakup, as well as bubble-induced turbulence. Empirical correlations are effective for some cases but may not be suitable for process scaleup and optimization in a more general sense. On the other hand, computational flluid dynamics (CFD) simulation is still far from mature, because of its difficulty in submodels for drag, lift, and turbulence models, and other issues, such as modeling bubble breakup and coalescence, grid resolution, and numerical algorithm,1 especially for higher gas holdup and higher gas flow rate existing in industrial processes. In our previous work, we attempted to propose a DualBubble-Size (DBS) model, featuring a stability condition used to close the mass and momentum equations sets.24 The stability condition was mathematically formulated with an extremum tendency, which physically reflects the compromise between two dominant mechanisms. The model calculation for total gas holdup showed a jump change when increasing superficial gas velocity, which is coincident with some experimental findings on regime variation from homogeneous and transitional regimes to heterogeneous regime.5,6 This jump change of gas holdup was found to be induced by the shift of the global minimum point of r 2011 American Chemical Society
so-called “microscale energy dissipation” within the space of structure parameters.7,8 This may offer a physical interpretation and reflect some intrinsic nature of regime transition. Although the DBS model offers a theoretical framework for analyzing the structure evolution of global gasliquid systems, there are some unsettled issues or open questions that require further investigation. For example, the computation cost is higher when using the global searching (GS) algorithm to seek the global minimum point of microscale energy dissipation formulated as a multivariation function. This article tries to employ a Simulated Annealing (SA) approach to fasten the computation and search for more-accurate solutions. Second, two-bubble classes were introduced in the original DBS model, but the gas phase in practical cases may appear in the form of a wide bubble size distribution. By using the SA approach, we find that the distinction of two bubble classes still prevails even though the gas phase is resolved into more classes in modeling, that is, extending the DBS model to TBS (Triple-Bubble-Size) or MBS (MultipleBubble-Size) models.
Special Issue: Nigam Issue Received: April 1, 2011 Accepted: August 24, 2011 Revised: August 9, 2011 Published: September 12, 2011 2077
dx.doi.org/10.1021/ie200668f | Ind. Eng. Chem. Res. 2012, 51, 2077–2083
Industrial & Engineering Chemistry Research
ARTICLE
2. MODEL DESCRIPTION The DBS model is an extension of the Energy-Minimization Multi-Scale (EMMS) method9,10 that was originally proposed for gassolid fluidization to gasliquid systems. In the EMMS model, the heterogeneous structure was resolved into a particle-rich dense phase and a gas-rich dilute phase, so that eight structure parameters were introduced but only six hydrodynamic equations could be established. A stability condition was then proposed to close the model (that is, the minimization of energy consumption for suspending and transporting particles per unit mass). The stability condition was thought to physically reflect the compromise of two dominant mechanisms, that is, the tendencies for gas to pass through the fluidized bed with least resistance and for particles to maintain the least gravitational potential, and, hence, mathematically formulated as a mutually constrained extremum tendency. Therefore, it supplies another constraint to the complex system, in addition to mass and momentum conservative equations. This strategy is adopted in the DBS model for gasliquid flow in bubble columns. The gas phase can be resolved into two-bubble classes as bimodal bubble-size distribution has been found to exist in experiments,11 and six structure parameters are introduced, that is, bubble diameters (dS, dL), volume fractions (fS, fL), and superficial gas velocities (UgS, UgL). Three conservative equations can be formulated as shown below: n
∑i Ugi ¼ Ug " fi FL g ¼
ði ¼ S, LÞ
ð1Þ
# fi π 2 1 Ugi Ul 2 C d ði ¼ S, LÞ D, i i 4 2 fi 1 fb ðπ=6Þdi 3
ð2Þ In bubble columns, the total rate of energy input fed into the system (NT) cannot be fully dissipated through liquid viscosity (Nturb). One portion is dissipated due to the slip of liquid along bubble surfaces and shape oscillation (Nsurf), and the remaining portion (Nbreak) may be stored temporarily as surface energy in the process of bubble breakage, released in bubble coalescence, and finally dissipated at microscale through the first two paths. While the first two portions are regarded as microscale energy dissipations, the last one could be taken as a type of mesoscale energy dissipation in a relative sense. Hence, the stability condition can be established as the minimization of microscale energy dissipation:24
∑
i ¼ S, L
Nsurf , i þ Nturb ¼ min
ð3Þ
It is equivalent to the maximization of mesoscale energy dissipation. The stability condition essentially reflects the compromise between two dominant mechanisms, namely, the tendency of gas-dominating and, hence, the least energy dissipation via liquid viscosity (N turbfmin), and the tendency of liquid-dominating and, hence, the weakest bubble surface response to liquid (Nsurffmin). The formulation for the three portions of energy dissipation is given below, and
more details about this model can be found in our previous publication.2,7,8 ! CD, p Nsurf ¼ 1 ð4Þ NT CD, b Z Nbreak ¼
db λmin
Z 0
0:5
ωðdb , λÞ Pb ðdb , λ, fBV Þcf πdb 2 σ dfBV dλ ð1 fb ÞFl þ fb Fg
ð5Þ Nturb ¼ NT Nsurf Nbreak
ð6Þ
The effective drag coefficient of a bubble can be calculated from the correlation of terminal velocity of Grace et al.12 With a given superficial gas velocity Ug, the model equations can be solved by separating three open variables (UgS, dS, dL) and searching the global minimum point of eq 3 within the space of the three structure parameters.
3. NUMERICAL METHODS In our previous work, the model was solved by using a global searching (GS) method. The model equations were solved for each point in the searching domain of the combination of the three structure parameters, and the global minimum point is discriminated among all the points. The model calculation could capture the jump change of total gas holdup and the cause inducing the jump change is found to be the shift of the global minimum point between the two local minima when increasing superficial gas velocity.2 However, the computational cost is higher for the GS method, and because the landscape of solution domain is fairly complicated and may involve multiple local minima, the accuracy of the final solution is dependent on the grid resolution. In this work, we try to use the Simulated Annealing (SA) method to overcome this problem. The SA method is a stochastic algorithm developed for the global optimization of combinatorial problems, and it is extensively used in chemical engineering such as optimization of heatexchanger networks13 and reactive distillation processes.14 Its working principle is borrowed from the physical annealing thermal process in solid melting through heating and then slow cooling and crystallization into a minimum free-energy state, allowing transitions from low to high energy levels when controlling the cooling rate. The Metropolis algorithm is often used to simulate this process. For a high energy state Ei at a temperature Ti, a new energy state Ej can be stochastically generated and then accepted as a new state if Ej e Ei. But Ej can also be accepted even if Ej > Ei, with a probability given by Ei Ej PðΔEi Þ ¼ exp KB T i where KB is the Boltzmann constant and Ti is the current system temperature. This can be achieved in calculation by generating a random number Ri for a given Ti and comparing Ri with the probability P(ΔEi) (i.e., the inner loop). Ej is accepted if the probability is greater than Ri. The next energy state Ej+1 will then be stochastically generated within the neighborhood of the last accepted state Ej. Apparently, the acceptance or rejection of the new state is also dependent on Ti. The probability approaches 1 at higher temperature and decreases upon cooling. This implies that the searching range is larger for higher temperature, and the 2078
dx.doi.org/10.1021/ie200668f |Ind. Eng. Chem. Res. 2012, 51, 2077–2083
Industrial & Engineering Chemistry Research
ARTICLE
Figure 1. Comparison of gas holdup calculated from the DBS model with Global Searching (GS) or Simulated Annealing (SA) methods (airwater system).
Figure 2. Trajectory of global minimum point in the space of structure parameters.
Table 1. Schemes for Searching the Minimum Point in the DBS Model node number
range of jump
(UgS dS dL)
change (m/s)
case
method
1
GS
100 1002
0.1090.110
2
SA
1.3 105
0.1010.102
3
GS
100 5002
0.1030.104
4
GS
100 20002
0.1010.102
5 6
GS GS
100 50002 100 200002
0.1010.102 0.1010.102
accepting criterion is more stringent when cooling, narrowing the searching range. The inner step can be performed for several numbers of iterations for a given temperature, and then Ti can be adjusted according to some cooling schedule (i.e., the outer step). A fast cooling schedule may drive the calculation toward a poor local optimum (i.e., a metastable state), whereas a slow schedule may take much computational time while approaching the global optimum. In this work, we use the following function to control the cooling rate:15 1=N
TðiÞ ¼ T0 αi
ð7Þ
α represents the speed of the cooling process and is usually defined in the interval [0.7, 1.0]. We adopt values of 0.95 for α and 2 for 1/N. Similarly, the initial temperature T0 is also important for the optimization. Higher-quality solutions may be obtained with higher initial temperature since the number of moves accepted is larger and the parameter space could be searched more thoroughly, but the computational time is longer. A criterion of balance between the quality of solution and computational cost can be established as suggested by Rayward-Smith et al.,16 starting the calculation with a higher T0 and then cooling it rapidly until 50% of bad solutions are accepted. T0 can be determined from the following formula: ΔE 0:5 ¼ exp ð8Þ T0 Cardoso et al.17 gave a detailed review for the SA method, and we only briefly introduce the basic idea of this approach in this article. Also, the control parameters of the algorithm for solving the DBS model in this work may not be optimal, and, here, our
Figure 3. Relative frequency of the bubble classes with increasing Ug.
single-minded purpose is to accelerate the calculation and, meanwhile, obtain a more-accurate solution. Another motivation of using the SA method is that, when the DBS model is extended to the TBS or MBS model, the GS method cannot be used, because the computational cost is extremely high with the rapid increase of free variables. Therefore, the SA method could be used for TBS or MBS models.
4. SOLVING THE DBS MODEL WITH THE SA METHOD The model equations must be solved for each node within the solution space (UgS, dS, dL) when using the global searching (GS) method, and the solution accuracy and the computational cost is dependent on the number of nodes. In our previous calculation, we used 100 5002 nodes to capture a jump change of the total gas holdup; this requires ∼3 h of computer time for one gas velocity with a personal computer that is equipped with an Intel 3G CPU, and ∼23 days to capture the jump change on the gas holdup curve. Figure 1 indicates the jump change occurring at superficial gas velocities between 0.103 m/s and 0.104 m/s, which is qualitatively coincident with the experimental findings of Zahradnik et al.,5 Camarasa et al.,6 and Ruthiya et al.18 They found that the second regime transition point appears beyond 0.12 m/s of Ug, signifying the transition from homogeneous and transitional regimes to heterogeneous regime. We then use the SA method to recalculate the gas holdup to shorten the computational time. However, the jump change point moves to the left, falling within the range of 0.1010.102 m/s. Compared to the GS method, the calculation takes only
