3D Numerical Simulations of Gas–Liquid Two-Phase Flows in

Jul 13, 2017 - The basic flow patterns were predicted successfully by both the CFD and CFD–PBM simulations, but comparatively big differences remain...
0 downloads 5 Views 2MB Size
Article pubs.acs.org/IECR

3D Numerical Simulations of Gas−Liquid Two-Phase Flows in Aluminum Electrolysis Cells with the Coupled Model of Computational Fluid Dynamics−Population Balance Model Shuiqing Zhan,*,† Zhentao Wang,† Jianhong Yang,*,‡ Ruijie Zhao,§ Changfeng Li,† Junfeng Wang,† and Jiemin Zhou∥ †

School of Energy and Power Engineering, ‡Institute of Green Materials and Metallurgy, and §Research Center of Fluid Machinery Engineering and Technology, Jiangsu University, Zhenjiang 212013, China ∥ School of Energy Science and Engineering, Central South University, Changsha 410083, China ABSTRACT: The 3D numerical simulations of gas−liquid twophase flows in three-anode and scaled-up 300 kA cell models using the computational fluid dynamics−population balance model (CFD−PBM) were presented. An Euler−Euler approach incorporated with a dispersed standard k−ε model and bubbleinduced turbulence was used. The CFD−PBM was first validated with referenced experimental data. Then the effect of the electromagnetic force (EMF) was investigated. The basic flow patterns were predicted successfully by both the CFD and CFD−PBM simulations, but comparatively big differences remained in terms of gas bubble distributions. The bubble coalescence behavior played a predominant role under the anode bottom regions. The bubble breakup process became increasingly important once the different-sized bubbles escape from anode edge regions. Moreover, the overall distributions of various flow parameters in the simulations considering EMFs are clearly different from the cases without EMFs in an industrial 300 kA cell. bubble behaviors and two-phase flows in the small-scale electrolytic systems cannot fully represent the actual electrolytic features. The large-scale air−water and low-temperature electrolytic models are both imperfect due to restriction on the experimental conditions, where the generation mechanism of anodic bubbles is completely different from the actual electrochemical reaction. Generally, most of the experimental approaches have revealed that the dynamics and the size distributions of gas bubbles in the cell depend on the spatial positions of the cell. Furthermore, it is clear that some key parameters (the local bubble size, the gas volume fraction, etc.) are very important for the accurate predicting of gas−liquid flows, but most of which remain somewhat limited to qualitative analysis, and some key problems still remain unsolved. To date, with the significant progress in numerical techniques and computing power, more and more studies on the complex gas−liquid flows have been carried out in aluminum electrolysis cells with modern computational fluid dynamics (CFD) methods. Various CFD models for predicting gas−liquid flow

1. INTRODUCTION In the process of high-temperature aluminum electrolysis, the anodic bubbles are mainly generated under the surface of the carbon anode. Then they slide along the anode bottom surface and escape from the top surface of the electrolyte, which can cause several local recirculations.1 The complex electrolyte flow pattern driven by the anode gas in the aluminum electrolysis cells has an obvious effect, not only on the homogenization of distributions of the alumina concentration and the temperature fields, but on the cell design, operational efficiency, and energy saving. Therefore, an in-depth knowledge of detailed gas−liquid flows is very important to the optimal design and scale-up of the electrolysis cell. Due to the extremely high temperatures and hostile environments, it is extremely difficult to conduct the detailed measurements in industrial cells. Therefore, both the experimental and computational approaches have been extensively studied for different types of aluminum electrolysis cells. The experimental approach mainly includes small-scale high-temperature electrolytic models2−5 and large-scale air− water6−10 and low-temperature electrolytic models.11−14 These experimental investigations were mainly focused on obtaining the comprehensive phenomenon of bubble formation, coalescence, breakup, and detachment and of gas-induced liquid flows. However, due to the scale limitations, the spatial © XXXX American Chemical Society

Received: Revised: Accepted: Published: A

April 26, 2017 July 10, 2017 July 13, 2017 July 13, 2017 DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research fields have been proposed on the basis of Euler−Lagrange,15,16 volume of fluid (VOF),17−19 and Euler−Euler approaches.20−24 In the Euler−Lagrange approach, an individual bubble is distinctly traced, which gives a direct physical interpretation. However, it is computationally intensive and therefore may not be appropriate for modeling the large industrial cells with higher gas volume fractions under the anode bottom regions. The advantage of the VOF approach is that it tracks and deals with the bubble interfaces directly. However, the main drawback is that it can consume a lot of computing power due to the need for fine-enough meshing, and hence, it cannot be directly used for industrial-scale cells either. To develop design tools for engineering purposes, the classical Euler−Euler approach is preferred due to its lower computational resources compared to the first two methods. Therefore, the Euler−Euler multiphase flow models have been extensively used and extended for simulating industrial-scale cells. For the twofluid model, a proper solution for the gas−liquid flows is mainly dependent on the correct modeling of interphase forces, turbulence models, and appropriate selection of bubble mean diameter. There is a general agreement that the influence of the interphase drag force is dominant compared with the other nodrag forces (for example, lift, virtual mass, wall lubrication, and turbulent dispersion forces) in the literature for different bubbly flow applications.25 Therefore, most studies numerically simulating gas−liquid flows in different types of aluminum electrolysis cells only considered the widely used drag force model of Ishii−Zuber26 between phases and fixed the bubble diameter as 10 mm20,21,24,27 or 5 mm.28 It is debatable whether these investigations have obtained reasonable gas−liquid flow patterns owing to the lack of experimental data for the validation. The first comprehensive and innovative framework for the emergence of experimental liquid flow patterns made by Cooksey and Yang9 using particle image velocimetry (PIV) is an important step for the validation of such CFD models. On the basis of their experimental equipment and measurement methods, similar experiments and numerical simulations were used to assess the effects of both the drag and turbulent dispersion forces by Feng et al.20,24,29 However, the Ishii− Zuber model and Schiller−Naumann model30 were used in two different investigations, respectively. Unfortunately, they do not explain the circumstances under which the model will show better performance, although the simulation results were found to be in reasonable agreement with experiments. However, it is worth mentioning that some numerical techniques, especially the reliable experimental information on the electrolyte velocity can provide us some new ideas about how to further study this topic. On the basis of the reliable experimental data,20 Zhan et al.22,23 developed gas−liquid two-phase simulations with the Schiller−Naumann model and a constant bubble size of 10 mm for different flow patterns in cells, but a relatively large error still exists. For the first time, a systematic investigation by Zhan et al.31 has been conducted to explore the collective effects of much more drag force and other no-drag force models together on the gas−liquid flows. Simulation results obtained also demonstrate that the drag force is indeed the predominant interphase force and that the lift, virtual mass, and wall lubrication forces can be ignored in the gas−electrolyte flows. Moreover, the Grace drag force model32 and the Simonin turbulent dispersion force model33 should be considered. Another major problem encountered in the simulations is the turbulence model. For the description of a turbulent enclosure

model between the gas and liquid phases in aluminum electrolysis cells, as far as we know, the standard k−ε model was adopted widely in all previous reports.20−24,27−29 Moreover, an additional term of effective turbulent viscosity proposed by Sato et al.34 has been used when considering the influence of the bubble-induced turbulence (BIT) in most investigations.20,21,24,27,28 Although some satisfactory results have been obtained with such an approach, the physical turbulence mechanisms are still not systematically explained and investigated. Therefore, different turbulence models (k−ε, k−ω, and Reynolds stress models) with a secondary phase consideration of a dispersed case have been employed and compared in detail with the experimental data20 in our previous work.31 Significative conclusions driven from the CFD predicted results indicate that the dispersed standard k−ε model is the most appropriate model. After going through previous investigations, it is evident that the CFD model is a useful tool to examine the gas−liquid flows in aluminum electrolysis cells. Moreover, the interphase forces and turbulence models used in the simulations have become a common cognition and can basically meet the requirements, which can provide the guidance and theoretical foundation for the current work. However, the assumption of a constant bubble mean diameter might bring some uncertainty and even unreasonable results for current Euler−Euler two-fluid CFD models. In fact, due to the large highly flattened anodes in modern industrial-scale cells, the bubble coalescence and breakup phenomenon can occur more often between different sized-bubbles, which result in more complex bubble dispersion and heat and mass transfer processes. Further, the presence of magneto-hydrodynamic (MHD) or electromagnetic force (EMF) may cause a series of effects on the electrolyte flows35 and more complex bubble coalescence and breakup behaviors. Therefore, further research on the fundamental theory of bubble size or bubble size distribution (BSD) should be conducted to better understand the performance in the cells. Recently, a well-known population balance model (PBM) was applied widely as an effective tool for calculating the statistical distribution information on a dispersed phase in various different multiphase flows,36−46 and the coupling of computational fluid dynamics−population balance model (CFD−PBM) was proposed to investigate the bulk multiphase flow fields by solving the CFD model as well as the BSD by solving the population balance equation (PBE). The CFD− PBM coupled model has been widely used in chemical engineering36−43 as well as in metallurgical engineering.44−46 At present, however, to our best knowledge, only Einarsrud et al.17,47 developed a multiscale approach to investigate the subgrid bubbles with the PBM method by coupling the VOF method in the two-dimensional (2D) small-scale computational domain, which cannot effectively and successfully acquire the detailed BSD information and macroscopic flow characters (especially when considering the effect of the MHD) in the three-dimensional (3D) computational domain for the modern large-scale industrial cells. In our own recently reported work,48 a CFD−PBM coupled model was used to simulate the BSD only in a one-anode model with a simplified bubble coalescence and breakup model without any direct parameter validation. Moreover, our previous work only considered the Schiller− Naumann drag force correlation (theoretically, the Grace model should be adopted) and ignored the turbulent dispersion force and the effect of the BIT, which is obviously not reasonable and accurate according to our latest research.31 B

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

the total interphase forces only include the drag force and turbulent dispersion force. The drag force reflects the anodic bubbles motion in the surrounding electrolyte and it plays a predominant role in the total interphase forces. The drag force can be formulated as

Actually, these factors mentioned above and the bubble diameter are not independent but interactive with each other, which can further influence the accurate prediction of the BSD and the macroscopic gas−liquid flows together. The primary objective of this study is to present a more reasonable and accurate CFD−PBM coupled model based on analysis of the appropriate coalescence and breakup mechanisms of dispersed anodic bubbles in the cell. For an initial quantitative comparison study, the CFD−PBM simulation considering the BSD and the CFD simulation with a constant bubble size were both performed in a three-anode model and are compared to each other with referenced PIV data.20,49 Then, the further qualitative validation and analysis were also proposed to study the gas−liquid flows and the BSD in the different flow regions of the cell. At last, the CFD−PBM simulation of a scaled-up 300 kA industrial aluminum electrolysis cell was proposed to acquire a better understanding of the impact EMFs on the BSD and gas−liquid flows.

FlD = −FgD =

(6)

where CD is the drag coefficient and db is the Sauter mean bubble diameter. In this work, the drag force experienced by a single bubble was calculated using the well-established Grace model32 based on our previous research,31 and it can be formulated as ⎛ ⎛ g d (ρ − ρ ) ⎞ 4 b l 8 g , ⎟⎟ C D = max⎜⎜min⎜⎜ 2 3 3⎠ ρ U t l ⎝ ⎝ ,

2. MODELING APPROACH 2.1. CFD Modeling of Two-Phase Flow. The Euler− Euler two-fluid approach is used in this work, whereby both the continuous phase (electrolyte) and the disperse phase (anode gas) are considered as interpenetrating media, identified by their local volume fractions. The mass and momentum transport equations for both phases are employed as ∇·(ρk αk uk) = 0

3 CD αgαlρl (ul − u g)|ul − u g| 4 db

Ut =

μg ρg db

Reb =

(1)

⎞ 24 (1 + 0.15Reb 0.687)⎟⎟ Reb ⎠

(7)

Mo−0.149(J − 0.857) (8)

ρ l |u l − u g |d b μl

(9)

where J and H are the empirical functions of the Morton (Mo) and Eötvös (Eo) dimensionless groups

T

∇·(αkρk ukuk) = −αk∇P + ∇·(αkμeff, k (∇uk + (∇uk) ))

⎧ 0.94H 0.757 2 < H < 59.3 J=⎨ ⎩ 3.42H 0.441 H > 59.3

(10)

⎛ μ l ⎞−0.14 4 ⎟ Eo Mo−0.149⎜ ⎝ 0.0009 ⎠ 3

(11)



+ Sk + Fk

(2)



where αk is the volume fraction of the gas or liquid phase, ρk is the phase density, uk is the phase vector velocity, and the subscript k is the phase (l for liquid and g for gas). μeff,k is the effective viscosity of the phase and P is the pressure. Sk is the momentum sources due to external body forces, which mainly include buoyancy and EMFs (the EMFs are ignored for the three-anode model). Fk is the total interphase forces between the two phases. The effective viscosity of the liquid phase μeff,l contains three components: the molecular viscosity (μl), the turbulence viscosity (μt,l), and the extra term due to BIT. μeff,l = μ l + μt,l + μBIT,l

H=

Mo =

Eo =

To account for the effects of gas bubbles on liquid turbulence in the momentum conservation equations mentioned above, an additional term of the turbulent viscosity μBIT,l due to BIT proposed by Sato et al.34 is used and can be calculated as (4)

where Cμb is a constant equal to 0.6. The effective viscosity of the gas phase (μeff,g) can be calculated from the effective viscosity of the liquid phase: μeff,g = μeff,l ρg /ρl

ρl 2 σgl 3

(12)

g(ρl − ρg )db 2 σgl

(13)

where σgl is the anodic gas/electrolyte interfacial tension (taken as 0.134 N/m). For comparison, the Schiller−Naumann drag force model30 was also used in the simulations. The detailed expressions can be found in our previous work.31 Besides the drag force in this work, the bubble experiences a turbulent dispersion force to account for the random influence of turbulent eddies on bubble diffusion. The Somonin averaged turbulent dispersion force derived by Simonin and Viollet33 is used in this study, which can be defined as

(3)

μBIT,l = Cμ bρl αgdb|ul − u g|

gμ l 4 (ρl − ρg )

FTD = −C TDK gl

(5) 31

As pointed out by our previous work, the effects of the lift and wall lubrication forces are almost negligible for the current gas−electrolyte systems. The virtual mass force is not included as the steady state of electrolyte flows and it is been our experience that it reduces model stability markedly. Therefore,

K gl =

Dt,gl ⎛ ∇αg ∇α l ⎞ ⎜⎜ ⎟ − σgl ⎝ αg αl ⎠⎟

3 CD αgαlρl |ul − u g| 4 db

(14)

(15)

where CTD is the turbulent dispersion coefficient with a value of 0.231 and Kgl is the drag momentum exchange coefficient. C

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Geometry and grid schemes of the three-anode model: (a) experimental geometry, (b) boundary conditions, and (c) grid schemes.

2.2.1. Bubble Coalescence Model. The closures for the PBM in the current work is the closure for the bubble coalescence and breakup mechanisms proposed by Luo and Svendsen,51,52 considering the bubble coalescence only caused by various turbulence eddies and the bubble breakup only due to the isotropic turbulence. The coalescence rate [c(di,dj)] is usually calculated as the product of the collision frequency [ωC(di,dj)] and the coalescence efficiency [PC(di,dj)]

The density ratio of electrolyte and anodic gas in the actual cell is larger than 5000 and the main flow regions contain one continuous phase (electrolyte) and a dilute secondary phase (anodic gas) in the form of anodic bubbles. It is exactly based on this consideration that the dispersed turbulence model is an appropriate model. Our prior numerical research experience31 indicates that the predicted results from the dispersed standard k−ε model could give better agreement with the experimental results than other dispersed turbulence models (standard k−ε, RNG k−ε, standard k−ω, SST k−ω, and the RSM models). Hence, the dispersed standard k−ε model will be considered in the following CFD−PBM simulation. The details of the relative equations and theories can be referred to ANSYS Fluent 14.550 and will not be covered in this work. 2.2. Population Balance Model (PBM). As mentioned above, only a constant bubble mean diameter can be considered in conventional Euler−Euler two-fluid model. Apparently, this is actually an obviously oversimplified calculation method, at least for the regions under the anode bottom. Due to the large highly flattened anodes in industrial cells, the bubbles gradually slide and gather under the anode bottom surfaces, resulting in a variety of possible bubble coalescence and breakup phenomenon. Then such phenomenon may contribute to bubble evolution with different numbers, sizes, and spatial distributions, which in turn affects two-phase flow patterns. In order to solve the problem, the formulation of the PBM was put forward in our current investigation. The PBM is a well-known approach for calculating the bubble size distribution and accounting for the coalescence and breakup effects in gas−liquid bubbly flows. For coalescence and breakup only, the steady-state population balance equation (PBE) can be expressed as ∇·(u gni) = BB − DB + BC − DC

c(di ,dj) = ωC(di ,dj) PC(di ,dj)

(17)

where di and dj are the diameter of bubble groups i and j, respectively. The collision frequency resulting from turbulent fluctuation with some modified forms proposed by Wang et al.53 is considered and expressed as ωC(di ,dj) =

αg,max π 2 (di + dj)2 (di 2/3 + dj 2/3)1/2 4 αg,max − αg ε1/3

(18)

The coalescence efficiency of the bubble groups i and j is calculated as ⎧ ⎫ [0.75(1 + xij 2)(1 + xij 3)]1/2 ⎪ ⎪ 1/2 ⎬ We PC(di ,dj) = exp⎨−c1 ij 1/2 3 ⎪ ⎪ (ρg /ρl + 0.5) (1 + xij) ⎭ ⎩ (19)

where C1 = 1.0 and xij = di/dj. Weij is the Weber number, which can be calculated as Weij = ρLdiui̅ j2/σl, where ui̅ j = (ui̅ 2 + ui̅ 2)1/2 and ui̅ = (1.69εidi)1/2. 2.2.2. Bubble Breakup Model. The bubble breakup model proposed by Luo and Svendsen52 is considered in this research. This model has been widely used in different previous studies of bubbly flows and is based on the theories of isotropic turbulence and probability. Assuming that binary breakup occurs between bubbles, the breakup rate of bubbles of size dj into bubbles of size di can be obtained as follows:

(16)

where ni is the number density function of group i bubbles. BB, DB, BC, and DC are the bubble birth and death rates caused by different bubble coalescence and breakup behaviors, respectively. The detailed theories can also be referred to ANSYS Fluent 14.5.50 D

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 2. Geometry and grid schemes of the 300 kA cell model: (a) geometry and boundary conditions and (b) grid schemes.

gB(di ,dj) = 0.9238(1 − αg)εl1/3di−2/3 exp( −bξ

∫ξ

min

1

(1 + ξ)2

that the PBM can effectively obtain the d32 assumption in the CFD model and the CFD model can offer the necessary gas− liquid two-phase flows inversely.

11/3

ξ

−11/3

) dξ

(20)

3. NUMERICAL DETAILS Two computational geometries and meshes, such as a threeanode model (see Figure 1) and a 300 kA aluminum electrolysis cell model (see Figure 2), were built in this work for different purposes. The computational domain of the threeanode models have the same dimensions as the experimental apparatus described by Feng et al.20 and Yang and Cooksey,49 respectively, which have been described carefully in our recent work.31 The CFD and CFD−PBM simulations will be validated with the experimental electrolyte velocities at the four locations of A, B, C, and D20 and the experimental electrolyte turbulent kinetic energies at location B49 (see Figure 1a). Unlike the interanode channel width of 40 mm in the experimental apparatus,20 the interanode channel width is 20 mm in the experimental apparatus.49 The 300 kA cell model is used as an application and the scaled-up case to investigate the relative influence of the gas−liquid flows induced by the EMFs on the BSD. The dimensions of the 300 kA cell model can also be referred from our previous work.23 A local grid refinement approach was presented for the anode−cathode distance (ACD) and the interanode channels (see Figures 1c and 2b). A number of hexahedral cells around 311 070 for meshing the three-anode model (for details, please refer to our previous study31) and 1 926 250 for meshing the 300 kA cell model was found to be sufficient to characterize the two-phase flow parameters after further grid refinements, which will be adopted in the present work. All the simulations are implemented using the commercial CFD software Fluent 14.5 with all the models described above. The double-precision model was selected for accurate simulation of the PBM. A mass-flow-inlet boundary condition with an average gas volume fraction of 0.5 was applied at the anode bottom surface for each anode on the basis of a

where gB(di,dj) is the breakup rate for a parent bubble with volume Vi into a daughter bubble with volume Vj. ξ = λ/d is the dimensionless size of eddies in the inertial subrange of isotropic turbulence. ξmin = λmin/d, f = dj3/di3, β = 2.047, and b = 12[f 2/3 + (1 − f)2/3 − 1]σlρl−1εl−2/3di−5/3β−1. 2.3. CFD−PBM Coupled Model. The class method (CM) dealing with different bubble classes (bubble groups) is adopted to solve the PBEs and can be implemented in the current CFD approach. Then, the steady-state PBEs for the ith bubble class can be expressed as ∇·(αgρg u gfi ) = BBi − DBi + BCi − DCi

(21)

By defining the volume fraction of group i bubbles as f i = αi/ αg, the sum of volume fractions for all different bubble groups equals the local gas volume fraction αg based on the assumptions that all size groups share the same density and velocity. So the number density function ni related to αg can be obtained from the following equation:

αgfi = niVi

(22)

The coupling between the CFD and PBM is achieved through the calculation of the bubble Sauter mean diameter d32, which is obtained as follows: N

d32 =

N

∑ nidi 3/∑ nidi 2 i=1

i=1

i = 1, 2, ..., N (23)

The d32 is then used in the calculation of the interphase forces between the gas and liquid phases and the BIT, thus completing the coupling relationship between the macroscopic CFD and mesoscopic PBM approaches. Both the CFD and the PBM can benefit from each other, with the advantage being E

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research hypothesized estimated average covering that fluctuated around 50% under the normal cell operation.6,8,54 The gas mass flow rates (i.e., 2.388 × 10−3 kg/s for three-anode model and 3.083 × 10−2 kg/s for 300 kA cell model) were obtained by a known Faraday’s law. A degassing boundary condition was considered at the top surface of the electrolyte, which means an outlet for the dispersed gas phase. A no-slip boundary condition was used for the electrolyte and a free-slip boundary condition was used for the gas bubbles at the wall, respectively. The physical properties of the gas and electrolyte can be found in our previous reports.23,31 On the basis of various electrolytic tests under investigated operating conditions,54,55 the average equivalent sizes of the initial bubbles that moved away from nucleation sites and the largest bubbles that escaped from the detachment positions were about 1 and 35 mm, respectively. As required by the solution strategy of the CM, the effects of the number of bubble-size groups (6, 12, and 18 groups) was studied with a user-customized method with sizes ranging from 1 to 35 mm. The results indicted that there was no significant difference between the 12 and 18 groups for the predicted gas volume fractions in the ACD regions. To reduce the computational resources and times, the subdivision of the different bubbles sizes into 12 groups with diameters ranging from 1 to 35 mm was used in the following CFD−PBM simulations, as shown in Table 1. A number of experiments indicted that the initial

For a quantitative validation, both the CFD and CFD−PBM simulations of gas−liquid two-phase flows in the three-anode model have been validated for the experimental conditions of Feng et al.,20 as shown in Figure 3. The predicted electrolyte velocity at three different sections (z = 0.08, 0.12, and 0.16 m) in four different vertical locations were compared with the experimental data of Feng et al.20 As expected from Figure 3, the simulated results of both the CFD and CFD−PBM simulations using the Schiller−Naumann model are not in particularly good agreement with the experimental data. However, the Grace model using nonspherical bubbles from the two simulations above can both further improve the predicted results. This is mainly due to the fact that the sizes and shapes of local bubbles can vary greatly and have an important effect on the two-phase flows, which heavily depend on the drag force coefficient models. The Schiller−Naumann model is only suitable for spherical rigid bubbles, while the Grace model can be largely applicable to the bubble regimes consisting of spherical and hemispherical caps and ellipse, which will be applied to the rest of this work. In addition, the results obtained from the CFD−PBM modeling slightly agree better with the PIV measurements compared to those obtained from CFD modeling, which indicates that the proposed CFD− PBM simulation could be suitable for predicting gas−liquid two-phase flows in aluminum electrolysis cells. An unexpected observation is that CFD−PBM modeling is more or less insensitive to the electrolyte velocity profiles. This might be because the gas volume fractions involved are very small in the side channels where the BSD is very narrow, and hence, the bubble coalescence and breakup are not significant in those regions. Figure 4 shows the quantitative comparison between the CFD-predicted and experimental electrolyte turbulent kinetic energy distributions at two different sections (z = 0.06 and 0.18 m) at vertical location B49 for the CFD and CFD−PBM simulations with the Grace model. Clearly, high turbulent kinetic energy mainly remains in the regions near the top surface of the electrolyte and anode ends related to the complex turbulence interactions between the gas and electrolyte phases. That is because of the strong momentum transfer from the rising bubbles to the electrolyte. However, it can be seen that the comparison is relatively not so good at some points, especially for the profiles in the center channel near the cell ends. This might be attributed to the complicated turbulence interactions and bubble-induced turbulence as well as their connections with bubble coalescence and breakup mechanisms, which needs to be addressed in future studies. It should be noted that, in reality, the bubbles have different average diameters at various regions of the cells (i.e., the side, interanode, and central channels and anode bottom regions). For instance, the cell more likely produces larger bubbles under the anode bottom regions due to the high coalescence rates. In contrast, the cell contains smaller bubbles in the side channels. In other words, if we can obtain more experimental data, such as the electrolyte velocity profiles under the anode bottom regions and even the bubble size information, the impact of the PBM may be significant. Further comparisons and discussions will be conducted in the following several sections. 4.2. Comparisons between the CFD and CFD−PBM Simulations for the Three-Anode Model. To understand more deeply the effects of incorporating the PBM, Figure 5 shows detailed comparisons of gas−liquid flows at a horizontal plane (z = 0.03 m) in the ACD regions for two different CFD

Table 1. Sizes of Discrete Bubble Groups group

bubble diameter/mm

group

bubble diameter/mm

1 2 3 4 5 6

1 2 3 4 5 7

7 8 9 10 11 12

10 13 17 22 29 35

bubbles generated under the anode bottom surface were very small and generally uniform in size, so the inlet bubble mean diameter was set to 1 mm (that is, the volume fraction of group 1 with a bubble size of 1 mm is 1.0). For comparison, the CFD simulation was also carried out with a constant bubble mean size of 7 mm. The mesh-to-mesh solution interpolation for EMFs in the 300 kA cell model and the coalescence and breakup kernels were incorportated through user-defined functions (UDFs).

4. SIMULATION RESULTS AND DISCUSSION 4.1. Quantitative Comparison with the Experimental Values. As stated in our previous work,31 the agreements with the experimental data were improved significantly with the drag force models of Grace, Tomiyama, and Ishii−Zuber compared with that of Schiller−Naumann (SN). Moreover, the predicted electrolyte velocity from the Grace model was shown to be relatively better than that from both the Ishii−Zuber and Tomiyama models. However, we should admit that the drag force models have a great relationship with the BSD. Therefore, In the current numerical study, we first might find it necessary to examine the effect of the two different drag force models of Schiller−Naumann and Grace on the accuracy of electrolyte flows predicted by the CFD simulations with a constant bubble mean size of 7 mm and by the CFD−PBM simulations considering bubble coalescence and breakup. F

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 3. Comparison of vertical electrolyte velocity profiles at different locations for different CFD models: (a) location A, (b) location B, (c) location C, and (d) location D.

the velocity magnitude of the flows predicted by CFD and CFD−PBM simulations, especially for the velocity magnitudes of the gas phase. This is perhaps to be expected, since the horizontal anode acts in the direction perpendicular to the flow and hence promotes the strong lateral movements of the bubbles. The above effects further confirmed that the great difference can be seen for gas volume fraction distributions with different simulations, although the gas layer morphology is quite similar. This has been also reported from the CFD simulations on the distributions of the gas volume fraction.27 The gas volume fraction is higher for CFD−PBM simulation than that for CFD simulation in the central regions. The amount of gas tends to move toward the edges of the anode and then escape upward to the top surface of the electrolyte. To further demonstrate and compare the gas bubble layer more directly, the gas volume fractions at the centerlines of four

models. Both the streamlines for the electrolyte and gas phases have been plotted to identify the overall flow patterns. From Figure 5, it is clearly shown that the basic distributions of the flow patterns (i.e., the scales and positions of vortexes) are very similar to each other and are predicted successfully by the two different simulations. The electrolyte flows into the ACD regions from the wider tap-end channel and then flows out into the narrower duct-end channel, which may have to do largely with the flow field structures of the gas phase. Besides, the flow patterns of both the electrolyte and gas under the bottom regions of anode B1 and anode B2 are quite different from those under the anode B3 bottom regions. These predicted results indicate that the channel width has a significant influence on the flow patterns. This feature is well-consistent with the previous results in the literature.24,29,49 However, this research also indicates that there are significant differences for G

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

simulation and CFD−PBM simulation is 6.21 and 7.16 mm, respectively, which is comparatively higher than some measurements (the maximum bubble layer thickness is about 4−6 mm) in the literature.3−6,55 This may be due to the difference for the scale and shape characteristic of the anode and the bubble dynamic behaviors between different electrolytic systems, which also needs to be addressed in future studies. It is obvious from the preceding discussions that there existed different impacts of the PBM on the two-phase flows in the side channels and ACD regions. The further explanations of the gas volume fractions at different vertical planes have been shown in Figure 7. These planes are plane A (a vertical plane through the middle of anode B2), plane B (a vertical plane through the middle of interanode channel between anode B2 and anode B3), and plane C (a vertical plane through the middle of anodes along the width direction), respectively. Also, it is clear from Figure 7 that there are significant differences for the magnitude of the gas volume fractions with different simulations, though they show similar gas volume fraction distributions. The gas volume fraction is lower in the side and central channels (see Figure 7a) than that in the interanode channels (see Figure 7b,c). The main reason maybe that the bubbles prefer to travel along the shortest path (anode width direction) because of the strong buoyancy effect of the rising bubbles, and hence, gas bubbles escape from the interanode channels more quickly than from the side and central channels (anode length direction). Moreover, a careful observation reveals that the gas volume fraction is very low in the side, interanode, and central channels, especially in the side and central channels, compared to that under the anode bottom regions. Therefore, the bubble behaviors like collision, coalescence, and breakup can be significantly promoted under the anode bottom regions, where the BSD should be very wide and hence the impact of the PBM will be higher. As previously discussed, the liquid turbulence is the primary mechanism responsible for bubble coalescence and breakup. In a given gas−liquid system, bubble coalescence and breakup processes occurring within the cell determine the BSD and are

Figure 4. Comparison of electrolyte turbulent kinetic energy profiles at location B for different CFD models with the Grace drag force model.

different z-sections (z = 0.035, 0.03, 0.025, and 0.02 m) along both the length and width directions of anode B2 for different CFD model predictions are shown in Figure 6. For each zsection, the predicted gas volume fraction profile has nearly a symmetric parabola shape for each direction and is significantly higher in the CFD−PBM simulations than in the CFD simulations. The gas volume fraction becomes lower when the fluid regions are farther away from the anode bottom surface. When the coalescence and breakup processes of bubbles are considered, different-sized bubbles are just formed under the anode bottom regions. These bubbles are more likely to carry lateral movements due to forces (including the magnitude and the direction) acting in the lateral direction for the large horizontal anode. More bubbles would not move and escape effectively from the edges of the anodes from CFD− PBM simulations. The predicted equivalent bubble layer thickness (the ratio of the accumulated bubble volume to the total areas of three anodes in ACD regions) for the CFD

Figure 5. Comparison of gas−liquid flows over a horizontal plane (z = 0.03 m) for different CFD models: (a) CFD simulation and (b) CFD−PBM simulation. H

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. Comparison of gas volume fractions at the centerlines of different z-sections along both the (a) length and (b) width directions of anode B2 for different CFD models.

Figure 7. Comparison of gas volume fractions over different vertical planes for different CFD models: (a) plane A, (b) plane B, and (c) plane C.

(see Figure 8a), which makes the smaller bubbles more likely to form larger bubbles, mainly due to the higher coalescence rate and lower breakup rate. That is why the electrolyte velocity profiles are relatively insensitive to the PBM implementation in the side and central channels (as pointed out above). It can be also observed from Figure 8 that the BSD is relatively uniform in the interanode channels, except in the regions near the top surface of the electrolyte. This can be attributed to the appropriate gas volume fractions (Figure 7b,c) and electrolyte turbulent energy dissipation rates (Figure 8a), which cause an approximate equilibrium of bubble coalescence and breakup processes.

primarily influenced by a combination of the local gas volume fraction and the liquid kinetic energy dissipation rate ε. The coalescence rates of smaller bubbles would increase for the higher local gas volume fraction, while the breakup rates of larger bubbles would increase for the higher local dissipation rate. Figure 8 illustrates the predicted electrolyte turbulent energy dissipation rate and bubble Sauter mean diameter d32 at different vertical planes for the CFD−PBM simulations. One can clearly observe that the local mean bubble size under the anode bottom regions is much larger than that in the side and central channels. Under the anode bottom regions, the gas volume fraction is higher (see Figure 7a), whereas the electrolyte turbulent energy dissipation rate is much lower I

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 8. (a) Electrolyte turbulent energy dissipation rate and (b) bubble Sauter mean diameter at different vertical planes for CFD−PBM simulations.

Figure 9. (a) Electrolyte turbulent energy dissipation rate and (b) bubble Sauter mean diameter at a horizontal plane (z = 0.03 m) in the threeanode model.

and the dissipation rate is also high, only the smaller bubbles can be found and the BSD is relatively uniform (mainly around 8−11 mm in diameter). Such phenomena seems expected but justified and has its theoretical basis. As a matter of fact, due to the horizontal large anode, the coalescence process of small bubbles often occurs during the lateral movements of bubbles from the central region to the boundary region, mainly along the width direction, as illustrated earlier in Figure 5b for the gas streamlines. Therefore, the larger bubbles are accumulated in regions near the interanode and duct-end channels, where the chance for bubble coalescence can occur continuously between the earlier large bubbles (also due to the earlier coalescence behavior), while the breakup behavior is not significant because of the low dissipation rate. Conversely, smaller bubbles (also due to the initial stages of the coalescence behavior from the smaller bubbles) mainly accumulated in central regions with low coalescence rate, although the gas volume fraction is very high. This is largely due to the fact that the initial mean size of

It is well-known that the BSD under the anode bottom regions have a great impact on the gas−liquid two-phase flows and cell design and scale-up. Unfortunately, there is no quantitative experimental BSD information available or proper bubble mean diameter for validation at present that we know of. Therefore, there exists a need to explore the detailed BSD and its evolution under the anode bottom regions, especially for the large anode configuration in modern industrial cells. The predicted electrolyte turbulent energy dissipation rate and bubble Sauter mean diameter at a horizontal plane (z = 0.03 m) are shown in Figure 9. From Figure 9, as discussed earlier, the BSD is very wide under the whole anode bottom region. Moreover, it is interesting to notice that the larger bubbles (mainly larger than 15 mm in diameter) are located at the boundary regions near the interanode and duct-end channels, where the gas volume fraction is very low (see Figure 5b) and the turbulent energy dissipation rate is also relatively low. For the central regions, although the gas volume fraction is higher J

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

three different directions (x-direction, y-direction, and zdirection). From Figure 10, it can be seen that the x-direction EMFs (Fx; see Figure 10a) are easily taken with the antisymmetrical distributions in the two sides (duct-end and tap-end channels) of the cell, with a higher value in the end regions and a lower value in the central regions. The distributions of y-direction EMFs (Fy; see Figure 10b) are unsymmetrical with a higher value in the upstream regions (pointing to the downstream side) and a lower value in the downstream regions (pointing to the upstream side). Further, it can be seen that the values of z-direction EMFs (Fz; see Figure 10c) are close to zero, except near the two very narrow side channels. The characteristics of the EMFs distributions suggested above could significantly affect the two-phase flows and hence further influence the bubble coalescence and breakup behaviors. The predicted overall distributions of the gas−liquid twophase flows and bubble Sauter mean diameter at a horizontal plane (z = 0.03 m) in the 300 kA aluminum electrolysis cell with and without EMFs are presented in Figure 11. It can be clearly seen that the velocity distributions of both the electrolyte and gas phases show significantly different flow patterns when considering the EMFs or not. Specifically, since the bubbles generated under each anode bottom surface move in various separate regions, the electrolyte flow patterns present a series of small eddies surrounding each anode when ignoring the EMFs (see Figure 11a). But when the EMFs are taken into account, the overall electrolyte flow patterns mainly present two large approximately symmetrical eddies with a few small eddies around several anodes near the central regions of the cell (see Figure 11b). The flow patterns of the anode gas show the same features as the electrolyte flow patterns for both the cases with and without the EMFs. And it noted that the magnitude for the gas velocity is larger near the interanode channel regions. In general, as discussed in the previous discussion, the EMFs do bring a serious influence on the horizontal two-phase flows, but they have almost no effect on the vertical flow patterns, mainly due to the distribution characteristics of EMFs mentioned before. Therefore, in the ensuing discussion, we will mainly explore the distributions of the gas volume fraction and electrolyte turbulent energy dissipation rate and their influence on the BSD at horizontal planes in the ACDs. The predicted distributions of the gas volume fraction and electrolyte turbulent energy dissipation rate are also shown in Figure 11. The gas volume fraction and the electrolyte turbulent energy dissipation rate in the case with the EMFs (Figure 11b) could be clearly different from the case without the EMFs (Figure 11a). The basic distribution characteristics of these two parameters under each anode bottom region are quite similar to each other and are essentially the same for the results in the three-anode model discussed earlier when not considering the EMFs. Once the EMFs are considered, the predicted gas volume fraction and electrolyte turbulent energy dissipation rate could change obviously, especially for the results in the regions near the two large symmetrical eddies. In these regions, the gas volume fraction is relatively low and the turbulent energy dissipation rate is relatively high, which will significantly influence the coalescence and breakup phenomenon of different-sized bubbles and hence may have an important effect on the BSD, as also shown in Figure 11. Obviously, one can observe from the contours in Figure 11 that the BSD under each anode bottom region is similar (for instance, the locations of the large and small bubbles) in the

all bubbles is assumed to be 1 mm (the smallest bubble size group, as shown in Table 1) at the anode bottom surface in the CFD−PBM simulation. As a result, the proportions of the volume fraction of several smaller-size groups (i.e., 1, 2, and 3 mm and so on) are relatively high, resulting in relatively small bubble diameters in central regions under the high gas volume fraction condition. Therefore, the CFD−PBM-simulated BSD in this paper is acceptable and the complex bubble coalescence and breakup behaviors are considered to be satisfactory. In summary, we can infer that the coalescence behavior is more predominant than the breakup behavior because of the higher gas volume fractions and lower turbulence during the bubble movements under the anode bottom regions and that the breakup of bubbles has become increasingly important during the transporting of bubbles from the anode edges to the top surface of the electrolyte. However, due to the restriction of existing experimental instruments or conditions as emphasized in the earlier sections, the validation presented in this paper did not extend to comparisons with bubble size measurements, which will be further investigated in future work. 4.3. Extending the CFD−PBM Simulation for a 300 kA Aluminum Electrolysis Cell Model. The good agreements of the current CFD−PBM simulations with experimental electrolyte velocity distributions for the three-anode model give us confidence that the further investigation of the hydromechanics and BSD performance in the larger-scale industrial cells will also be fruitful. As an extension of the three-anode model, the 3D steady-state CFD−PBM simulations of gas−liquid flows in a 300 kA aluminum electrolysis cell were performed. It is generally known that both the EMFs and bubble-driven forces have a big effect on the two-phase flows in real industrial cells.21,23,28,44,56,57 The EMFs are included through the Lorentz force in the cell, which are calculated from the magnetic and electric fields obtained from the metal-pad model using ANSYS. Figure 10 shows the predicted distributions of the EMFs in

Figure 10. Predicted distributions of the EMFs in three different directions in a 300 kA aluminum electrolysis cell: (a) x-direction, (b) y-direction, and (c) z-direction. K

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 11. Gas−liquid two-phase flows, electrolyte turbulent energy dissipation rate, and bubble Sauter mean diameter at a horizontal plane (z = 0.03 m) in a 300 kA aluminum electrolysis cell model: (a) without EMFs and (b) with EMFs.

bubble coalescence and breakup mechanisms. The predicted results from the CFD−PBM simulation considering the BSD and the CFD simulation with a constant bubble mean size have been compared and analyzed in details. Further, the CFD− PBM simulation of a scaled-up 300 kA industrial aluminum electrolysis cell was proposed to acquire a in-depth understanding of the impact of the EMFs on the gas−liquid flows and the BSD. To sum up, several conclusions can be drawn: (1) The simulations with the Grace model were proven to be better than those with the Schiller−Naumann model. The CFD− PBM simulation can show better predictions than the CFD simulation by simply assuming a constant bubble mean size of 7 mm, from which we can infer that the CFD−PBM simulation can give a relatively excellent prediction of the gas−liquid twophase flows. (2) The basic distributions of the flow patterns under the anode bottom regions are predicted successfully by both the CFD and CFD−PBM simulations. However, the gas volume fraction distributions predicted by the two different simulations are markedly different and the gas volume fraction is higher for CFD−PBM simulations in the central regions. (3)

case without the EMFs. More specifically, the larger bubbles can be observed in the edge regions of the anode near the interanode channels and small bubbles appear in the central regions and the edge regions of the anode near the two side channels, which are also the same as the results mentioned in the three-anode model above. The possible reasons for these results will not be repeated here. In this situation, when considering the EMFs, the gas volume fraction is lower, while the electrolyte turbulent energy dissipation rate is higher in the regions near the two large symmetrical eddies, which may easily result in smaller bubbles because of the higher breakup rate and lower coalescence rate. By far, although there have been no available data information on bubble mean size or BSD for the present systems, the CFD−PBM-predicted BSD appears to be in accordance with expectances.

5. CONCLUSIONS In the present work, a 3D CFD−PBM simulation for predicting the gas−liquid two-phase flows in two aluminum electrolysis cells of different scale was established by considering proper L

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

(5) Zhao, Z. B.; Wang, Z. W.; Gao, B. L.; Feng, Y. Q.; Shi, Z. N.; Hu, X. W. Anodic bubble behavior and voltage drop in a laboratory transparent Aluminum electrolytic cell. Metall. Mater. Trans. B 2016, 47, 1962−1975. (6) Fortin, S.; Gerhardt, M.; Gesing, J. A. Physical modeling of bubble behaviour and gas release from aluminum reduction cell anodes. Light Met. 1984, 721−741. (7) Chen, J. J. J. Some physical model studies of gas-induced flows in aluminum cells. JOM 1994, 46, 15−20. (8) Vekony, K.; Kiss, L. I. Morphology of two-phase layers with large bubbles. Metall. Mater. Trans. B 2010, 41, 1006−1017. (9) Cooksey, M. A.; Yang, W. PIV measurements on physical models of aluminum reduction cells. Light Met. 2006, 359−365. (10) Wang, Y. F.; Zhang, L. F.; Zuo, X. J. Fluid flow and bubble behavior in the aluminum electrolysis cell. Light Met. 2009, 581−586. (11) Qian, K.; Chen, Z. D.; Chen, J. J. J. Bubble coverage and bubble resistance using cells with horizontal electrode. J. Appl. Electrochem. 1998, 28, 1141−1145. (12) Xue, Y. Q.; Zhou, N. J.; Bao, S. Z. Normal temperature analogue experiment of anode bubbles behavior in aluminum electrolysis cells. Chin. J. Nonferrous Met. 2006, 16, 1823−1828. (13) Alam, M.; Yang, W.; Mohanarangam, K.; Brooks, G.; Morsi, Y. S. Investigation of anodic gas film behavior in Hall-Heroult cell using low temperature electrolyte. Metall. Mater. Trans. B 2013, 44, 1155− 1165. (14) Das, S.; Weerasiri, L. D.; Jegatheesan, V. Bubble flow in a static magnetic field. Light Met. 2015, 789−793. (15) Li, X. P.; Li, J.; Lai, Y. Q.; Zhao, H. Q.; Liu, Y. X. Mathematical simulation of the gas induced bath flow in a drained aluminum reduction cell. Trans. Nonferrous Met. Soc. China 2004, 14, 1221−1226. (16) Doheim, M. A.; Elkersh, A. M.; Ali, M. M. Computational modeling of flow in aluminum reduction cells due to gas bubbles and electromagnetic forces. Metall. Mater. Trans. B 2007, 38, 113−119. (17) Einarsrud, K. E.; Johansen, S. T. Modelling of bubble behaviour in aluminium reduction cells. Prog. Comput. Fluid Dyn. 2012, 12, 119− 130. (18) Zhang, K. Y.; Feng, Y. Q.; Schwarz, M. P.; Wang, Z. W.; Cooksey, M. Computational fluid dynamics (CFD) modeling of bubble dynamics in the aluminum smelting process. Ind. Eng. Chem. Res. 2013, 52, 11378−11390. (19) Zhao, Z. B.; Feng, Y. Q.; Schwarz, M. P.; Witt, P. J.; Wang, Z. W.; Cooksey, M. Numerical modeling of flow dynamics in the aluminum smelting process: comparison between air−water and CO2cryolite systems. Metall. Mater. Trans. B 2017, 48, 1200−1216. (20) Feng, Y Q.; Yang, W.; Cooksey, M. A. Development of bubble driven flow CFD model applied for aluminum smelting cells. The Seventh International Conference on CFD in the Minerals and Process Industries; CSIRO, 2009; pp 1−8. (21) Li, J.; Xu, Y. J.; Zhang, H. L.; Lai, Y. Q. An inhomogeneous three-phase model for the flow in aluminium reduction cells. Int. J. Multiphase Flow 2011, 37, 46−54. (22) Zhan, S. Q.; Zhou, J. M.; Li, M.; Zhou, Y. W.; Yang, J. H. Numerical simulation of gas−liquid two-phase flow in aluminum reduction cells with perforated anodes. CIESC J. 2013, 64, 3612−3619. (23) Zhan, S. Q.; Li, M.; Zhou, J. M.; Yang, J. H.; Zhou, Y. W. CFD simulation of effect of anode configuration on gas−liquid flow and alumina transport process in an aluminum reduction cell. J. Cent. South Univ. 2015, 22, 2482−2492. (24) Feng, Y. Q.; Cooksey, M.; Schwarz, M. P. CFD modeling of electrolyte flow in aluminium reduction cells. Light Met. 2007, 339− 344. (25) Oey, R. S.; Mudde, R. F.; Van den Akker, H. Sensitivity study on interfacial closure laws in two-fluid bubbly flow simulations. AIChE J. 2003, 49, 1621−1636. (26) Ishii, M.; Zuber, N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 1979, 25, 843−855. (27) Yang, S.; Zhang, H. L.; Xu, X. J.; Zhang, H. H.; Zou, Z.; Li, J.; Lai, Y. Q. Effects of slot cutting at prebaked anodes on bubble

The bubble coalescence behavior is more predominant than the bubble breakup behavior under the anode bottom regions due to the higher gas volume fractions and lower turbulence. The bubble breakup process has become increasingly important, while the bubbles with different sizes transport from the anode edges to the top surface of the electrolyte. (4) The EMFs have a serious influence on the horizontal two-phase flows but have almost no effect on the vertical flows, mainly due to the distribution characteristics of EMFs. Moreover, the overdistributions of the gas volume fraction, the turbulent energy dissipation rate, and the bubble Sauter mean diameter in the case with the EMFs are clearly different from the case without the EMFs, mainly in the regions near the two large symmetrical eddies induced by the horizontal EMFs. In general, the above results demonstrate the reasonableness of the CFD−PBM simulations in this work and their ability to predict the complex gas−liquid two-phase flows in the scale-up industrial cells. Future work will be focused on improving the models for more appropriate coalescence and breakup mechanisms of the bubbles and the turbulence modification for the case of a high gas volume fraction with different bubble groups. Moreover, the transient irregular anode bottom geometry related with the electrolyte−metal interface deformation in a real cell should also be further discussed. Therefore, even with the current results, we expect that the implementation of PBM into the conventional CFD method is feasible and essential.



AUTHOR INFORMATION

Corresponding Authors

*S.Z. e-mail: [email protected]. *J.Y. e-mail: [email protected]. ORCID

Shuiqing Zhan: 0000-0002-9759-6511 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors are grateful for the financial support of the National Natural Science Foundation of China (11502097, 51509110), the Natural Science Foundation of Jiangsu Province (BK20170551, BK20150511, BK20171301), the Natural Science Foundation of Higher Education Institutions of Jiangsu Province (17KJB450001), the Foundation of Senior Talent of Jiangsu University (2015JDG158), the China Postdoctoral Science Foundation (2016M591781), the “Qing Lan” Project Foundation of Jiangsu Province and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). Our special thanks are due to anonymous reviewers for insightful suggestions on this work.



REFERENCES

(1) Liu, Y. X.; Li, J. Modern aluminum electrolysis; Metallurgical Industry Press, Beijing, China, 2008. (2) Aaberg, R. G.; Ranum, V.; Williamson, K. The gas under anodes in aluminum smelting cells. Part II: gas volume and bubble layer characteristics. Light Met. 1997, 341−346. (3) Cassayre, L.; Plascencia, G.; Marin, T.; Utigard, S. T. Gas evolution on graphite and oxygen-evolving anodes during aluminium electrolysis. Light Met. 2006, 379−383. (4) Yang, S. H.; Xie, B. R.; Yang, F. L.; Li, M. Z. The behavior of anode bubbles in aluminum electrolytic cell. Nonferrous Met. Sci. Eng. 2013, 4, 20−24. M

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research elimination in aluminum reduction cell. J. Cent. South Univ. 2012, 43, 4617−4625. (28) Wang, Q.; Li, B.; He, Z.; Feng, N. Simulation of magnetohydrodynamic multiphase flow phenomena and interface fluctuation in aluminum electrolytic cell with innovative cathode. Metall. Mater. Trans. B 2014, 45, 272−294. (29) Feng, Y. Q.; Schwarz, M. P.; Yang, W.; Cooksey, M. Two-phase CFD model of the bubble-driven flow in the molten electrolyte layer of a Hall−Héroult aluminum cell. Metall. Mater. Trans. B 2015, 46, 1959−1981. (30) Schiller, L. A.; Nauman, Z. A drag coefficient correlation. Ver. Dtsch. Ing. 1935, 77, 135−147. (31) Zhan, S. Q.; Yang, J. H.; Wang, Z. T.; Zhao, R. J.; Zheng, J.; Wang, J. F. CFD simulation of effect of interphase forces and turbulence models on gas−liquid two-phase flows in non-industrial aluminum electrolysis cells. JOM, 201710.1007/s11837-017-2327-5. (32) Grace, J. R.; Wairegi, T.; Nguyen, T. H. Shapes and velocities of single drops and bubbles moving freely through immiscible liquids. Institution Chem. Eng. 1976, 54, 167−173. (33) Simonin, O.; Viollet, P. L.; Hewitt, F. G. Modelling of turbulent two-phase jets loaded with discrete particles. Phase-Interface Phenom. Multiphase Flow. 1990, 26, 259−269. (34) Sato, Y.; Sekoguchi, K. Liquid velocity distribution in two-phase bubble flow. Int. J. Multiphase Flow 1975, 2, 79−95. (35) Baiteche, M.; Taghavi, S. M.; Ziegler, D.; Fafard, M. LES turbulence modeling approach for molten aluminium and electrolyte flow in aluminum electrolysis cell. Light Met. 2017, 679−686. (36) Cheung, S. C. P.; Yeoh, G. H.; Tu, J. Y. On the numerical study of isothermal vertical bubbly flow using two population balance approaches. Chem. Eng. Sci. 2007, 62, 4659−4674. (37) Cheung, S. C. P.; Yeoh, G. H.; Tu, J. Y. Population balance modeling of bubbly flows considering the hydrodynamics and thermomechanical processes. AIChE J. 2008, 54, 1689−1710. (38) Liang, X. F.; Pan, H.; Su, Y. H.; Luo, Z. H. CFD−PBM approach with modified drag model for the gas−liquid flow in a bubble column. Chem. Eng. Res. Des. 2016, 112, 88−102. (39) Gao, Z. M.; Li, D. Y.; Buffo, A.; Podgórska, W.; Marchisio, D. L. Simulation of droplet breakage in turbulent liquid−liquid dispersions with CFD−PBM: comparison of breakage kernels. Chem. Eng. Sci. 2016, 142, 277−288. (40) Chen, H.; Sun, Z.; Song, X. F.; Yu, J. G. Key Parameter prediction and validation for a pilot-scale rotating-disk contactor by CFD−PBM simulation. Ind. Eng. Chem. Res. 2014, 53, 20013−20023. (41) Chen, H.; Sun, Z.; Song, X. F.; Yu, J. G. A pseudo-3D model with 3D accuracy and 2D cost for the CFD−PBM simulation of a pilot-scale rotating disc contactor. Chem. Eng. Sci. 2016, 139, 27−40. (42) Xie, L.; Luo, Z. H. Multiscale computational fluid dynamics− population balance model coupled system of atom transfer radical suspension polymerization in stirred tank reactors. Ind. Eng. Chem. Res. 2017, 56, 4690−4702. (43) Yan, W. C.; Li, J.; Luo, Z. H. A CFD−PBM coupled model with polymerization kinetics for multizone circulating polymerization reactors. Powder Technol. 2012, 231, 77−87. (44) Zhan, S. Q.; Li, M.; Zhou, J. M.; Yang, J. H.; Zhou, Y. W. CFD simulation of dissolution process of alumina in an aluminum reduction cell with two-particle phase population balance model. Appl. Therm. Eng. 2014, 73, 805−818. (45) Chen, G. J.; He, S. P.; Li, Y. G.; Wang, Q. Modeling dynamics of agglomeration, transport, and removal of Al2O3 clusters in the Rheinsahl−Heraeus reactor based on the coupled computational fluid dynamics-population balance method model. Ind. Eng. Chem. Res. 2016, 55, 7030−7042. (46) Liu, Z. Q.; Qi, F. S.; Li, B. K.; Cheung, S. C. P. Modeling of bubble behaviors and size distribution in a slab continuous casting mold. Int. J. Multiphase Flow 2016, 79, 190−201. (47) Einarsrud, K. E.; Eick, I.; Bai, W.; Feng, Y. Q.; Hua, J. S.; Witt, P. J. Towards a coupled multi-Scale, multi-physics simulation framework for aluminium electrolysis. Appl. Math. Model. 2017, 44, 3−24.

(48) Zhan, S. Q.; Li, M.; Zhou, J. M.; Yang, J. H.; Zhou, Y. W.; Zhou, C. Q. A CFD−PBM coupled model predicting anodic bubble size distribution in aluminum reduction cells. Light Met. 2014, 777−782. (49) Yang, W.; Cooksey, M. Effect of slot height and width on liquid flow in physical models of aluminum reduction cells. Light Met. 2007, 451−456. (50) ANSYS Fluent User Manual, Ansys Inc., Canonsburg, PA, 2011. (51) Luo, H. Coalescence, breakup and liquid circulation in bubble column reactors (Doctoral thesis), Norwegian Institute of Technology, Trondheim, Norway, 1993. (52) Luo, H.; Svendsen, H. F. Theoretical model for drop and bubble breakage in turbulent dispersions. AIChE J. 1996, 42, 1225−1233. (53) Wang, T. F.; Wang, J. F.; Jin, Y. Population balance model for gas-liquid flows: influence of bubble coalescence and breakup models. Ind. Eng. Chem. Res. 2005, 44, 7540−7549. (54) Poncsak, S.; Kiss, L. I.; Toulouse, D.; Perron, A.; Perron, S. Size distribution of the bubbles in the Hall−Heroult cells. Light Met. 2006, 457−462. (55) Cassayre, L.; Utigard, T. A.; Bouvet, S. Visualizing gas evolution on graphite and oxygen-evolving anodes. JOM 2002, 54, 41−45. (56) von Kaenel, R.; Antille, J.; Romerio, M. V.; Besson, O. Impact of magnetohydrodynamic and bubbles driving forces on the alumina concentration in the bath of an Hall-Héroult cell. Light Met. 2013, 585−590. (57) Das, S.; Brooks, G.; Morsi, Y. Theoretical investigation of the inclined sidewall design on magnetohydrodynamic (MHD) forces in an aluminum electrolytic cell. Metall. Mater. Trans. B 2011, 42, 243− 253.

N

DOI: 10.1021/acs.iecr.7b01765 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX