An investigation of the effect of initial bubble diameter on the bubble trajectory in the flotation column cell using CFD simulation

The effect of initial bubble diameter on the bubble motion pattern in a flotation column has been studied by the twophase computational fluid dynamics (CFD) method. The two-phase simulations have been done using the volume of a fluid (VOF) model in ANSYS® Fluent® software. The computational field was a square cross-section column with a width of 0.1 m and a height of 1 m into which air was interred as a single bubble from the lower part of the column by an internal sparger. An experimental test has been also performed and the simulated results have been validated using the values obtained for the bubble rise velocity. A comparison of the simulation and the experimental results has confirmed that CFD can predict the bubble rise velocity profile and its value in the flotation column less than 5% relative to the experimental values. Then the simulations have been repeated with a 20% decrease and increase in the initial bubble diameter to investigate the effect of bubble diameter on the bubble flow pattern. The investigations have shown that as the bubble diameter increases, the velocity decreases and the bubble rises in a more zigzag direction as a result of two counter-rotating trailing vortices behind the bubble increasing.


Introduction
With the decline of minerals in the present age, the maximum usage of natural resources and mines has become very important and this has given mining and its related industries a special place in the economy. Therefore, increasing the efficiency of processing operations and optimizing related systems is of particular importance. Understanding the mechanism and manner of performing physical processes is effective in improving the performance of a system. Fluid science and in particular "Computational Fluid Dynamics (CFD)" is the tool that can be used in this regard. Process efficiency in the column flotation is strongly influenced by hydrodynamic components and CFD as a numerical method can help analyze the effect of these components on the flow system. So far, examples of CFD simulations have been performed on processing equipment with different approaches and objectives (Wang et al., 2018). Table 1 lists the titles and summarizes the studies performed on the flotation column collection area using the CFD simulation method.
In these research studies, most attempts have been made to investigate the overall flotation operations. However, due to the special role of bubbles in floating the particles and transporting them to the froth phase and for a better understanding the effect of hydrodynamic components, it is necessary to conduct more studies on them.
The bubble trajectory (the rising path) is also one of the most prominent variables, which affects the hydrodynamics of flotation (Levich, 1962

Not mentioned
To study the dynamic characteristics of oil-water separation in loop flotation.
while some prominent aspects of bubble dynamics in the flotation process, such as the effect of initial bubble diameter on bubble rising trajectory and rise velocity still need to be researched in future studies. Even though the effect of bubble diameter on the rise velocity has been investigated before, the interaction of initial bubble size, rise velocity and bubble trajectory has not yet been investigated. The present study aims to investigate the trajectory and rise velocity of different bubble diameters in a single-bubble system by the use of a video processing technique and CFD simulation to establish an insight into the bubble dynamics in the flotation process. Therefore, to achieve a reliable result in simulation, firstly a complete understanding of the governing physics phenomenon is a need, and secondly, a selection of suitable physical models (multiphase model, turbulence model, etc.) is required to use with simulation tools (e.g. Fluent software). Therefore, this study has To measure the cross-section and axial section of gas-liquid two-phase flow field in a lab-scale cyclonic flotation column, particle image velocimetry (PIV) is combined with an endoscopic measurement and phase discrimination technique. been defined and carried out concerning the undeniable importance of improving the efficiency of flotation as the most efficient separation unit and increasing the sciences related to flotation in both experimental and simulation. In the first step, to simplify the problem, the motion of a single bubble was investigated experimentally and numerically. First, the experiment, its details and then the method used in the simulation and its steps are explained. Fi-nally, experimental and simulation results are presented, and by comparing the two, a conclusion is made.

Experiments
A laboratory flotation column was constructed and an experimental test was carried out to obtain necessary experimental data to validate the simulation predictions. In the following subsections, the flotation column characteristics and setup are explained.

Column flotation cell characteristics and setup
Column dimensions: the flotation column with a height of 1.00 m and a square cross-section with a crosswise of 0.100 m and a height to diameter ratio of 7.14 was used, which is shown schematically in Figure 1. The square cross-section makes it easy to record experimental observations. The column material was selected from clear Plexiglas with a thickness of 0.005 m; because it is strong and light and its transparency allows for the observation of the process. There are four holes in the column; two holes at the bottom at a height of 4 cm and two holes at the top of the column at a height of 0.940 m.
This system is equipped with an internal sparger, including a pneumatic hose with an inner diameter of 0.002 and an outer diameter of 0.004 m, which is located in the lower hole of the column, and an automatic air pump (ACD-500) with a power of 5 W, a pressure of 0.012 MPa and the capacity to inject 3 litres of air per minute was employed for aeration.
In the present work, a bubble entered the system individually; in fact, the time interval between the entries of two consecutive bubbles into the column was adjusted in a way that the bubbles' velocities were not affected by their frequency. Since the inlet airflow to the column is very small in these conditions, it is not possible to use conventional airflow measuring instruments. Therefore, to adjust the inlet airflow to the column, a needle valve was used after the pump.
The experiment was performed at atmospheric pressure and an air temperature of 30±0.5 °C. The purpose of the experiment was to investigate the hydrodynamic conditions in the flotation column in the presence of two phases, gas and liquid. Therefore, first, the column was filled with water to a height of 0.900 m (considering the possible flow fluctuations) and then the air valve was opened and air entered the column. After ensuring hydrodynamic stability, imaging and recording of the results were performed.

Image analysis
After ensuring the stability of the flow inside the column, imaging of the column and recording the results were performed using a Canon EOS 500d. The video recording speed was 30 frames per second and the continuous shooting speed was 3.4 frames per second. Figure 2 summarizes the imaging process and the purpose of doing it, which is explained in the following: • Measuring the bubble diameter: for this purpose, the image of the bubble was recorded in the inlet, and then in Digimizer 1 defined a unit of measurement from the ruler installed on the column wall, and the diameter of the bubble was measured from the image. To calculate the bubble diameter, a volumetric equivalent diameter was calculated according to Equation 1 (Clift and Grace, 1978). (1) Where: d eq -the volumetric equivalent diameter (m); h -the smallest diameters of the bubble (m); b -the largest diameters of the bubble (m).
• Measuring the bubble displacement in two continuous images and calculating the bubble rise velocity: to measure the bubble rise velocity, the camera's continuous shooting mode settings were used. By measuring the displacement of a bubble in two consecutive images and considering the time interval between two images (1/3.4 = 0.29 s), the velocity of the bubble was calculated. For this purpose, photos were taken in four sections from 0 to 0.20 m, 0.20 to 0.45 m, 0.45 to 0.70 m, and 0.70 to 0.90 m. • Counting the number of bubbles in the column and calculating the gas holdup: an overall view of the column was taken and the number of bubbles inside the column was counted. By multiplying the number of bubbles that are in the column at the same 1 Digimizer is a flexible software package that is very useful for analyzing images and allows you to manually accurately measure the components of images. time in the volume of one bubble, the volume of air in the column was obtained, and by dividing it by the volume of water in the column, the gas holdup was calculated. • Counting the number of interred bubbles to the column per second: to measure the inlet air flow rate to the column, in stable conditions, the location of the sparger was filmed and then the number of interred bubbles to the column was counted per unit time and the inlet air flow rate was calculated by considering the volume of bubbles.

Modelling and Simulation
The geometry and meshing of the computational domain were created using Gambit ® 2.4.6 and then imported into ANSYS ® Fluent ® 15 for flow calculations by a computer system with a Xeon processor consisting of 42 cores used in parallel mode. The details of numerical models, column geometry, mesh generation, choice of the optimal mesh convergence criteria, and simulation conditions are discussed in the following subsections.

Governing Equations
In the present work, the Volume Of Fluid (VOF) model was used as a multiphase model to simulate a bubble motion in the water, which allows the construction of the interface to become part of the solution based on the same grid system. The movement of the gas-liquid interface was tracked based on the distribution of α g , the volume fraction of gas in a computational cell, where α g =0 in the liquid phase and α g =1 in the gas phase. Therefore, the gas-liquid interface existed in the cell where α g lies between 0 and 1. The interfacial mass transfer was neglected here. The mass and momentum equation for the two-phase flow are represented as Equations 2 and 3: (2) Where: ρ -density (kg/m 3 ); -the velocity vector (m/s); p -the scalar pressure value (Pa); τ -the viscous stress tensor (Pa); -the gravitational acceleration (m/s 2 ); -the surface tension source term (N/m 2 ); t -time (s). In a two-phase system, the phases are represented by the subscripts g and l, the density and viscosity in each cell are given by Equations 4 and 5 (Brennen, 2005) as follows: Where: ρ -density (kg/m 3 ); α g -the volume fraction of gas-phase (%); ρ g -the gas density (kg/m 3 ); α l -the volume fraction of liquid-phase (%); ρ l -the liquid density (kg/m 3 ); µ -viscosity (Pa.s); µ g -the gas viscosity (Pa.s); µ l -the liquid viscosity (Pa.s).
The tracking of the interface between the gas and liquid is accomplished by the solution of a continuity equation for the volume fraction of gas, which is in Equation 6: Where: α g -the volume fraction of gas-phase (%); -the velocity vector (m/s); t -time (s). The liquid volume fraction is calculated based on Equation 7: Where: α g -the volume fraction of gas-phase (%); α l -the volume fraction of liquid-phase (%).

Geometry
The solution geometry is based on a laboratory flotation column used for experiments to validate the CFD simulations. According to the experimental observations in the present work, the movement of a single bubble and its rise path is almost located in the middle of the column; therefore, the flow is assumed to be linear and the simulations are performed in two dimensions. The schematic of the geometry and the boundary conditions are shown in Figure 3.

Meshing
After generating geometry in Gambit software, a meshing operation is performed. The desired mesh is quad. The criterion of meshing was located 20 to 30 cells in the bubble diameter direction. Since the range of transverse oscillation of the bubble is (2* (d bubble + 0.5d bubble )), so for meshing, the width of the column is divided into three parts: the middle part of the column, where the  bubble is interred and begins to oscillate in an upward motion; and two parts on the sides. In the bubble rise path -the middle part of the column-the grid is smaller (Cano-Lozano et al., 2014). Figure 4 shows the mesh in the middle of the column, around the bubble, and the sides.

Simulation and boundary conditions
In this section, choosing and setting simulation and boundary conditions, such as transient flow and time step, two-phase model, convergence criteria, and the discretization scheme are mentioned.
• Transient flow: the purpose is to study the movement of a bubble during its ascent, so all simulations are done time-dependently. • Time step: in solving the problem transiently, it is necessary to divide the processing time into small parts. A good way to determine the amount of time step is to observe the number of iterations done in each time step to achieve convergence. The appropriate number of repetitions in each time step is 10 to 30 (Ansys Fluent Theory Guide, 2015). In this research, the variable time step has been used. Thus, the calculations start with the time step 10 -8 and based on the value determined for the Courant number (0.5) can be increased to 10 -4 . It is necessary to mention two points: firstly, the mentioned values are used only based on the author's experience, and secondly, it has the maximum accuracy that can be used according to the available hardware for the authors. • Laminar flow: since the motion of a single bubble is studied, the flow is laminar and the turbulence model is not applied (Van der Pijl, 2005). • Two-phase flow: all simulations are performed by considering two phases: water as the primary phase and air as the second phase. The properties of fluid used in the simulation at a temperature of 30±0.5 ºC (the experimental temperature) are shown in Table  2. Also, the amount of surface tension is considered to be 0.0712 (N/m). It should be noted that since the ratio of minimum diameter to maximum diameter is less than one percent, so the interred bubble is assumed to be spherical (Clift and Grace, 1978).
tion less than 10 -6 , were considered as the convergence criteria. • The discretization scheme: the discretization of the momentum equations was performed using the first-order upwind scheme, whereas the QUICK scheme was used for the volume fraction equation. The PISO algorithm was utilized to couple the pressure and the velocity fields.

Mesh verify
It is important to use enough computational cells in the solution domain to solve the equations. The effect of the mesh size on the accuracy of the solution is investigated by studying the bubble rise velocity. In this way, the number of meshes along the bubble diameter and consequently other parts of the computational domain increases to the extent that the values obtained from the simulation have the least difference with the experimental results. Table 3 shows three different mesh numbers used for this purpose. The results of these simulations are compared to the bubble rise velocity profile in Figure 5. • Gravitational acceleration: in all simulations, a gravitational acceleration of 9.81 m/s 2 is considered in the y-direction. • Convergence criteria: the normalized residual values, for continuity less than 10 -3 and volume frac- As can be seen in Figure 5, meshes 2 and 3 predict the velocity pattern more similar to the experimental results than mesh number 1. The numerical values also show that by decreasing the mesh number and using mesh number 3 in the simulations, the velocity values slightly change. Therefore, with an approximation of less than 2%, mesh number 2 has been selected and used in the simulations.

Results and discussion
An experimental test and the CFD simulations were performed in the mentioned conditions. After validating the numerical results, simulations were repeated with a 20% increase and decrease in bubble diameter. The result of the experiment, validation details, and discussion about the result obtained are presented in the following subsections.

Experimental results
During the experiment, the hydrodynamic variables of the column were measured and calculated according to the algorithm presented in Figure 2 and the results are shown in Table 4.

Dimensionless numbers
Reynolds number: quantity without a significant unit used in fluid mechanics to predict flow pattern. This number is the ratio of inertial force to viscous force. Fluids with lower Reynolds numbers tend to have flow with a laminar, layered pattern, while fluids with high Reynolds numbers have more turbulent flows. The Reynolds number value is calculated by Equation 8: Where: Re -Reynolds number. ρ 1 -the liquid density (kg/m 3 ); u b -the bubble velocity (m/s); d b -the bubble diameter (m); µ l -the liquid viscosity (kg/m.s). Morton number: used in fluid dynamics with the Eo number to describe the shape of a bubble moving in a fluid which is calculated by Equation 9: Where: Mo -Morton number; G -the gravitational acceleration (m/s 2 ); µ 1 -the liquid viscosity (kg/m.s); ρ 1 -the liquid density (kg/m 3 ); σ -surface tension (Pa). Eotvos number: The ratio between the internal force to the surface stress acting on the surface between the liquid and the gas. According to the definition, its value is obtained by Equation 10 (Binder, 1973): Where: Eo -Eotvos number; G -the gravitational acceleration (m/s 2 ); ρ 1 -the liquid density (kg/m 3 ); ρ 1 -the gas density (kg/m 3 ); d b -bubble diameter (m); σ -surface tension (Pa). The values of dimensionless numbers for the bubble diameter of this work are provided in Table 5.

Validation
The next step is comparing the simulation results with the experimental data to validate the simulation results. In the experiment, the velocity of the bubble was calculated by measuring its displacement in two continuous shots, and due to the time interval between two continuous shots, a mean velocity was achieved. Also in the simulation, a centre of mass was defined for the bubble and the value of mean velocity was calculated according to its displacement in defined height per unit time. Figure 6 shows the bubble mean velocity profile for a bubble with an initial diameter of 3.23 mm based on numerical simulations and experimental measurements. As can be seen, the bubble mean velocity profile has two parts, in the first part the bubble velocity increased immediately, and in the second part it decreased with a decreasing acceleration. The axial component of the bubble velocity was zero at the beginning of the movement,

Effect of initial bubble diameter on bubble path rise
The domain of bubble rise path depends on bubble characteristics such as bubble size, deformation, changes in rise velocity, and the presence of surfactants To investigate the effect of initial bubble diameter on the bubble shape, the simulations were repeated for the other two diameters with a 20% increase and decrease.
It was observed that the bubble shape and path changed from the bottom to the top of the column. For a bubble with a diameter of 3.23 mm, at the entering of the bubble to the column and the beginning of the movement, as long as the bubble maintained its spherical shape, it moved in a straight upward direction, but at the same time it accelerated and increased velocity, and the bubble shape changed from sphere to ellipse. As a result Also, the velocity values recorded in the experimental test and the predicted values in the simulation, and the difference between the two are given in Table 6. It shows the numerical values obtained from the simulations predicted the velocity values with a difference of 5% compared to the experimental values. Furthermore, the location of the bubble in the Clift diagram (the pattern of deformation bubble rising in the liquid column) was investigated as a qualitative validation, which is shown in Figure 7. In the bubble size range of the present work, with increasing bubble diameter, Reynolds increased and the surface tension force decreased, other forces (viscosity, inertia, and compressive dynamics) increased and the bubble shape changed from spherical to elliptical and then a higher increase in bubble diameter made a spherical cap bubble shape (Clift and Grace, 1978). of increasing the level of the forehead and flattening the shape of the bubble, its increase in velocity also decreased, and as a result, the shape of the bubble became closer to a spherical shape again. It is the alternating changes in the shape and velocity of the bubble that caused the bubble to ascend in a zigzag and spiral manner (Zhang, Y.; Sam, A.; Finch, J.A., 2003). Actually, as Figure 8 shows, when a bubble rises in water, it experiences a lateral force and torque along its path and two counter-rotating trailing vortices appear behind the bubble, which is believed to be the main causes of path transitions (deviation) from straight to zigzag. The rising trajectories for bubbles of 1.92, 2.56, and 3.23 mm diameter are illustrated in Figure 9. The time interval between two bubbles in each shape is 0.2 seconds. A comparison of bubble trajectory in different diameters indicates while the diameter of the bubble increased between 1.92 mm to 3.2 mm, the motion changed from straight to zigzag. In fact, with an increase in bubble size, its deviation from the straight vertical path as well as its zigzag motions increased. tation cell using a two-phase VOF model was able to predict the bubble rise velocity values well with an error of less than 5% compared to the experimental values. Finally, the bubble rise path in three different diameters was presented and investigated. The results show that the morphology and pattern of bubble motion were affected by the initial bubble diameter, and by increasing the bubble size in the range of 1.92 mm to 3.23 mm, its deviation from the straight vertical path increased as well.
While in flotation, as a method of particle separation using a bubble, on the one hand, the zigzag motion of the bubble increases the probability of bubble-particle collision and on the other hand, because of these oscillations, the probability of the particle slipping off the bubble and separating is more probable, so it is very important to choose the optimal bubble range produced by a sparger in flotation.
These findings can help researchers better understand the bubble properties and their dynamics, which are fundamental to the performance of the flotation process, and how they are affected by the initial bubble diameter. In addition, since the relationship between the bubble dynamics and the performance of true flotation is one of the least understood aspects of the flotation, it is worthwhile to do further theoretical investigations to develop a more comprehensive model, which is the current work of this research team.

Conclusion
In this research, the single bubble motion in the flotation column was investigated in two experimental and numerical sections. First of all, an experiment was performed and the bubble hydrodynamic components, including its rise velocity, were measured using image processing. After that, the motion of a single bubble in the flotation column was simulated in two phases using a VOF model. To mesh verify, simulations were performed for three mesh numbers; the results were compared to the bubble rise velocity profile and finally, with an approximation of less than 2%, mesh number 2 by 752500 elements was selected and used in the simulations. Then the bubble rise velocity profile was plotted and compared with the experimental values. The results show that the bubble motion simulation in a column flo-