CFD analysis of laboratory scale phase equilibrium cell operation

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Jama, Mohamed Ali; Nikiforow, Kaj; Qureshi, Muhammad Saad; Alopaeus, Ville


I. INTRODUCTION
In the chemical engineering related industry, many processes are multiphase systems, where components are transferring across interfaces between the phases. Thus, mass transfer has been subjected to extensive studies, experiments, and modeling since the dawn of the whole chemical engineering discipline more than a hundred years ago. Many of these studies during the past two decades have used computational fluid dynamics (CFD) coupled with mass transfer. In the field of drinking water disinfection, Cockx et al. (1999), Talvy et al. (2011), and Gong et al. (2007) have modeled the ozonation process to predict ozone dissolving. Laakkonen et al. (2006a;2006b) have investigated mass transfer in agitated fermenter vessels considering local population balance distribution. In separation process analysis for liquidliquid extraction in micro and macroscale (Buffo et al., 2016;Tsaoulidis and Panagiota, 2015;Xu et al., 2013;Chasanis et al., 2010;and Jildeh et al., 2012) or distillation process (Jiang et al., 2013;Sun et al., 2007;Rahimi et al., 2012;Wang et al., 2004;Rahimi et al., 2006;and Noriler et al., 2010), CFD is used to predict flow field, mass transfer, concentration profiles, and finally efficiencies. Importantly, all these studies a) Author to whom correspondence should be addressed: mohamed.jama@ aalto.fi. Tel.: +358 50 434 7372. Fax: +358 94 702 2694. b) Email: kaj.nikiforow@aalto.fi c) Email: muhammad.qureshi@aalto.fi d) Email: ville.alopaeus@aalto.fi have a common assumption of thermodynamic equilibrium at the interface between the two fluid phases. Accurate prediction of this equilibrium is necessary to specify driving forces for mass transfer correctly.
In the classical equilibrium stage design of separation processes, phase equilibrium is assumed at each ideal stage (Seader andHenley, 1998 andKing, 1980). In order to predict the number of ideal stages correctly for the desired operation, the requirements for accurate phase equilibrium data are extremely high, typically much higher than for transport or other physical properties. Šoóš et al. (2003) have compared results from designing rectification distillation column using VLE data from measurements, Aspen HYSYS database with those predicted from the UNIFAC (functional-group activity coefficients) method. The predicted distillation column parameters, i.e., reflux ration, flows, and composition of the phases on the stages were different when using the three different VLE data sets. This indicates that the design results are highly sensitive to changes in the equilibrium data (Xin and Whiting, 2000;Whiting et al., 1999;Dong et al., 2005;Mathias, 2014;and Wakeham et al., 2001). This stresses the importance of accurate equilibrium data.
As CFD coupled with mass transfer is becoming the state of the art modeling tool for bench, pilot, and large scale of two phase systems in chemical engineering, little or zero attention has been paid to laboratory scale apparatuses where chemical equilibrium data are measured. The quality of the data depends on the experimental procedure and on the apparatus used.  Qureshi et al. (2015) have conducted measurements for infinite dilution activity coefficients (IDAC) for various bio-oil related components and fitted thermodynamic model parameters against the experimental results. IDAC measurements are carried out in the region where the solute is in trace level in the solution. In the IDAC cells, gas is bubbled through the liquid phase, which is also mixed to improve the contact between gas and liquid and to level out concentration differences in the liquid phase. In the work of Qureshi et al. (2015), equilibrium was assumed based on an earlier work by Richon et al. (1980). Approach to equilibrium was further assessed by sensitivity analysis, for example, by observing the variation in the measured activity coefficient with respect to gas flow rates. This experimental validation approach is, however, sensitive to small experimental errors. Furthermore, the approach to equilibrium depends strongly on the solutes being measured, and thus separate validation runs should be done for each new chemical system to be measured. If the equilibrium is not reached, the experiments have to be stopped and modification to the apparatus done, which delays the experimental campaign. The first step to improve this design and operation purely based on empiricism is to use the mass transfer model based on idealized flow pattern as in the work by Richon et al. (1980). However, with these idealized mass transfer models, it is not possible to assess if some of the bubbles are carried with the liquid flow to the surface in a shorter time than required for proper equilibration; this by-passing would not be compensated by bubbles remaining in the cell longer than required as no further mass transfer takes place for them. Only a proper fluid flow modeling would fully answer this.
The most important physical phenomena affecting the operation of the IDAC and several other equilibrium cell designs are relatively close to the operation of a separation stage in distillation or absorption. This suggests that the same state-of-the art modeling tools should be used to analyze both. In this work CFD, Euler-Lagrangian approach is used for the first time to model a practical laboratory scale equilibrium apparatus (dilutor cell), such as used by Richon et al. (1980). The RTD (residence time distribution) of the bubbles inside the dilutor is studied using CFD and compared with the residence time from ideal mixing, plug flow, and the time required for the bubbles to reach equilibrium in the dilutor cell. With this analysis, the fraction of bubbles reaching equilibrium can be estimated to give an insight into whether the equipment works properly, how the equipment should be designed in case of poor equilibration, and whether simpler models are sufficient for this kind of analysis.

II. APPARATUS AND VLE MEASUREMENT
The dilutor cell, which is subject to modeling in this work, has been designed and used for IDAC measurements by Qureshi et al. (2015), from where the details of the apparatus and its operating principles can be found. The dilutor cell has a diameter of 24 mm and working volume of either 18 cm 3 or 27 cm 3 , corresponding to liquid levels of either 40 mm or 60 mm, respectively. The mixing in the dilutor is induced by magnetic stirrer bar. There are no baffles in the  (Qureshi et al., 2015) dilutor. Nitrogen, which is used as the carrier gas, is introduced in the cell by evenly spaced capillaries, with 0.2 mm diameter. Inert gas, in this case nitrogen, is used for stripping the dilute solute from a solution and solute concentration is detected from the inert gas stream exiting the apparatus by using gas chromatography (GC). Six bio-oil based compounds (n-butanol, valeraldehyde, furan, crotonaldehyde, toluene, and acrolein) are used as solutes in the work of Qureshi et al. (2015) and are tested from a numerical point of view in this study. The measured infinite dilute activity coefficients for the above-mentioned compounds are shown in Table I.

III. MODELING
The steady state solver with the Euler-Lagrangian (E-L) approach was used to model the gas-liquid fluid flow in the dilutor, explained in Secs. III A and III B. The Euler-Euler (E-E) approach with post-processing for the bubble tracking was also investigated and compared with the E-L approach.
In the Euler-Euler (E-E) approach: first, steady-state twophase calculations were conducted for fluid flow. Then UDF (user-defined function) code of the equation for particle age (Liu and Tilton, 2010;Baléo and Cloirec, 2000;and Simcik et al., 2012) is used to calculate the age of the bubbles (a) in the dilutor cell. The mean age function and its multiphase scalar transport representation are (Simcik et al., 2012) where D T is the flow dispersivity and Sc is the turbulent Schmidt number, with default value of 0.7. A passive scalar transport equation, in dispersed phase, is solved in a transient simulation. The value of the scalar (the age of the bubbles) at gas inlet stream is set to zero. At the walls and outlet, the flux of the scalar is set to zero. Then the area weighted average value is traced at the liquid level (the top of the computational domain) to get the mean residence time. Very similar mean residence times were obtained by the E-E and E-L methods, see Table II. This indicates that the methods are equally accurate in this case and that the results are supporting each other.

A. Hydrodynamics modeling
A single phase hydrodynamics is solved for water, representing the dilute aqueous mixture. The Lagrangian based discrete phase model (DPM) is used to model the RTD of the nitrogen bubbles in the dilutor. In the Lagrangian framework, each particle is tracked individually. The residence time for each bubble is detected on the liquid surface. In this approach, the system is two-way coupled. The continuous phase flow field is affected by the discrete phase (nitrogen bubbles) and vice versa. Due to the very dilute solvent used in the present experimental setup, the physical properties of the continuous phase are taken to be those of water. Nitrogen bubbles with constant size of 1 mm or 2 mm are used in the discrete model. Due to the small scale and a low number of bubbles in the system, application of traditional bubble breakage and coalescence models is not taken into account (Bakker and Van den Akker, 1994;Lane et al., 2005;Liu et al., 2013;Moilanen et al., 2005;and Moilanen et al., 2008).
Single phase Eulerian mass and momentum conservation equations were used, where the interface momentum exchange, the source term (S i ) is determined by the bubble force balance equation of the Lagrangian approach. The modeling process is carried out as follows: (1) the velocity of the continuous phase is determined by solving the Euler conservation equations, (2), the velocity of the dispersed phase is solved by the force balance equation, and (3) the source term for the momentum equation (5) is updated. Thus, the numerical iteration between these steps takes into account the interaction between phases.

B. Discrete phase model
The discrete phase model of the Fluent 16 software (Fluent, Inc., 2006) was used to track bubbles in the Lagrangian reference frame. The trajectory of each bubble is determined by the force balance equation, where the subscript i denotes the bubble number, F D is the drag function, u p the bubble velocity, u l the continuous phase velocity, g is the gravitational acceleration constant, and ρ p and ρ l are the bubble and continuous phase densities, respectively. The two first right-hand side terms are drag and gravity forces, where the last one F i represents additional case specific forces per unit bubble mass.
The drag function is defined as where µ l is the molecular viscosity of the continuous phase, D p is the bubble diameter, Re is the Reynolds number [see Eq. (9)], and C D is the drag coefficient. Bubbles are considered as monosized and incompressible in the drag coefficient correlation (Morsi and Alexander, 1972), shape deformation is neglected by using shape factor equal to one, where a 1 , a 2 , and a 3 are empirical constants.
In the two-way coupled model, bubbles are likewise affecting the continuous phase flow in the dilutor. To model accurately, this mutual interaction between the continuous and discrete phases, more forces are introduced to the force balance equations. The virtual mass force is considered as significant in the Lagrangian bubble tracking calculations. This assumption is supported by Ziegenhein et al. (2015).
The added formula for virtual mass force is Bubble trajectories are assumed to be significantly affected by the pressure gradient of the continuous fluid. The used pressure gradient force formula is written as Two-way turbulence coupling is also employed to model turbulent production and damping between phases (Fluent, Inc., 2006). Bubbles are tracked using the stochastic tracking (random walk model) approach of the Fluent 16 software (Fluent, Inc., 2006). Thus, the turbulent dispersion affecting on bubbles is taken into account. Bubbles with the same injection position may take different paths due to the turbulent fluctuations. The bubble trajectories are affected by the instantaneous and fluctuation velocities of the continuous phase, as follows: whereŪ is the mean velocity and u is the fluctuation velocity of the continuous phase (Fluent, Inc., 2006). The force balance equation is integrated to calculate the instantaneous bubble velocity. Thereafter, the trajectory of each bubble is determined by the kinematic

IV. TEST CASES AND NUMERICAL DETAILS
The modeled operating conditions are set equal to those in the laboratory (Qureshi et al., 2015). Modeling is done for impeller speeds of 300-1000 rpm and inert gas flow rates of 2-10 cm 3 /min. Based on visual observations, the size of the bubbles is chosen as 1 mm or 2 mm. 1 mm represents the most typical bubble size observed in the system, whereas 2 mm represents the largest one with the highest mass transfer limitations. In the modeling, the size of the bubbles is assumed constant. Although this assumption may not be precisely valid in practise, a constant bubble size assumption still gives reasonable estimate for residence times and characteristic time scale for mass transfer. The residence time of each bubble size is modeled separately as described earlier. The effect of the liquid level on the residence time distribution is also investigated by using a liquid level of either 40 mm or 60 mm. The process properties and the physical properties of the solvent and the solutes are shown in Tables III, IV and V. The commercial software of ANSYS Fluent R16.0 is used in the hydrodynamic modeling. The Fluent solver package is based on the finite volume method to solve transport equations   (Fuller et al., 1966) (Wilke and Chang, 1955) numerically (Patakar, 1980 andVegteeg andMalalasekera, 1995). The GAMBIT software is employed to conduct the 3D model of the dilutor and mesh generation. The generated mesh is unstructured and contains 130 000 computational cells. The stirrer is modeled as a zero thickness wall of the cross section, no slip velocity boundary condition with standard wall functions was assumed (Fig. 1). Also the capillaries, which distribute the gas bubbles in the dilutor, are modeled as cylindrical elements with 1 mm height, see Fig. 2. The effect of the capillary tubes on the fluid flow is considered negligible. The gas inlet of the sparger is located on the surface pointing downwards, corresponding to the tip of the capillary. The interior geometry of the dilutor is divided into two parts. One part with a cylindrical shape is located at the center of the dilutor and contains the impeller with its surroundings. The other part consists of the rest of the dilutor, i.e., the regions close to the dilutor walls.
The multiple reference frame model was implemented for rotational movement in the cell. The Reynold-Averaged Navier-Stokes with the k ε turbulence model is employed in modeling turbulent flow in the cell. The stirrer Reynolds number is 5000 for 1000 rpm, which is in the transition flow regime. In order to further assess this problem, fluid flow was also predicted by assuming laminar flow. These two approaches resulted in almost identical residence time distributions for the bubbles, indicating that the original model formulation can be used with trust. For the continuous phase, the top boundary (liquid surface) was set as the pressure outlet and the rest were set as the non-slip wall boundary.
The particle type in the discrete phase model of the Fluent R 16.0 was set as an inert particle with the physical properties of nitrogen. The stream injection is defined as a surface. All dilutor walls are defined as reflect boundaries and the liquid surface as the escape boundary for the particles.

V. MASS TRANSFER
In the dilutor, the solute is being stripped from the solvent into the carrier gas stream. Mass transfer limitations may exist either in gas phase, in liquid phase, or in both phases. To study the rate of mass transfer separately in both phases, a quantity called the approach solute concentration is used. This quantity is defined as the ratio of solute concentration to equilibrium solute concentration, For proper operation of the dilutor, the approach solute concentration should be relatively close to unity. In this work, an approach solute concentration equal to τ(t) = 0.999 is used as a reference and the time needed to reach this limit is called the mass transfer limiting residence time (t LIM ). That is, the mass transfer limiting residence time is the minimum gas bubble residence time required (dictated by the mass transfer rate) to reach 99.9% of equilibrium solute concentration.
In the case of mass transfer limitation being in the gas phase, the approach solute concentration is calculated employing an approach by Alopaeus (2000), Time dependent Fourier number is where t is the residence time of the bubble and the time averaged Sherwood number is The solute-carrier gas binary diffusion coefficients, shown in Table V, are calculated by the correlation of Fuller et al. (1966) where the subscript i denotes the solute, j denotes the carrier gas, and ϑ is the group contribution value for atoms, groups, and structural features. Both the molecular weight and group contribution value are averaged by considering the humidity in gas phase (Qureshi et al., 2015).
In the case of mass transfer limitation being in the liquid phase, the approach solute concentration is calculated employing an approach proposed by Richon et al. (1980), where p vap is the solute vapor pressure and γ ∞ is the infinite dilution activity coefficient of the solute in solvent. A correlation by Weber (1975) is used for calculating the liquid phase mass transfer coefficient (k liq ), The model by Wilke and Chang (1955) is used to calculate diffusion coefficients in liquid phase. The residence time weighted approach equilibrium concentration of gas leaving the dilutor is calculated as follows: where τ(t) is the approach solute concentration, E(t) is the residence time distribution, and t is the residence time.

VI. RESIDENCE TIME DISTRIBUTION
The time available for the mass transfer depends on the residence time of the carrier gas bubbles in the dispersion. A central factor affecting this time is the rate of mixing. Two limiting cases can be identified: (1) in case of perfect mixing, the residence time distribution is that of ideal mixing, and (2) in case of no mixing, bubble residence time is constant and can be calculated with the plug flow assumption.
The ideal mixing residence time distribution E ideal (t) is (Fogler, 2006) In this work, the ideally mixed system mean residence time (t M ) is taken to be the same as that calculated with CFD, where E model is the residence time distribution from CFD modeling. This assumption was made to make comparisons easier, although simple literature correlations for gas hold-up could have been used as well.
In the plug flow calculations, the residence time is constant due to the constant size dependent slip velocity and equal bubble path length for all the bubbles, which is equal to the distance between dilutor liquid surface and capillary tip. The plug flow residence time is defined as follows: where the bubble rise velocity (v S ) is calculated by using the same drag coefficient correlation as in the CFD solution, namely, Morsi and Alexander (1972).

A. Residence time distribution
The effect of the bubble size on the approach equilibrium concentration was studied considering two different bubble diameters: 1 mm and 2 mm. Figure 3 shows the RTD calculated with CFD and also the plug flow and ideal mixing RTD for two  Because the maximum and average observed bubble sizes in the laboratory were between 1 and 2 mm, respectively, the true RTD is close to that of 1 mm. The mean residence time was studied by keeping the agitation speed constant and changing the inert gas flow rate between 2 and 10 cm 3 /min. For 1 mm bubbles, Fig. 4 shows that the higher inert gas flow rate, the more bubble driven hydrodynamics is in the dilutor. In other words, the flow field in the dilutor cell is affected by the rising bubbles so that the mean residence time is shortened (mean RT of the model and the plug flow RT are getting closer). For the 2 mm bubbles, the CFD model mean RT predictions and plug flow RT are equal. This indicates that the gas hold-up is the same. The same drag model used in the CFD model and in plug flow calculation explains the overlapping between the mean RT and plug flow RT ( Table V).
The mass transfer limiting RT time [i.e., the time required to reach 99.9% of equilibrium concentration, calculated with Eqs. (15) and (19)] for both bubble sizes and for all compounds are shown in Tables VI and VII. The operating conditions are the same as those in the case shown in Fig. 3. The results show that mass transfer in the gas phase is relatively fast. Mass transfer in the liquid phase is also fast except for furan and toluene, which are the least soluble compounds. Figure 5 shows both the RTD calculated with CFD and the mass transfer limiting RT (calculated assuming mass transfer resistance totally either in the gas phase or in the liquid phase). For both bubble sizes, the RTD with 4 cm and 6 cm liquid levels is calculated using CFD. The inert gas flow rate was 10 cm 3 /min and agitation speed 300 rpm. Considering mass transfer resistance in gas phase only, the equilibrium is reached and is independent of the studied bubble sizes and liquid level [see Fig. 5(a)]. The influence of the water level on approach concentration of toluene is illustrated in Fig. 6. The bubble slip velocity in plug flow is used to make Eqs. (15) and (19) liquid level dependent. The results show that the equilibrium would be reached with liquid level less than 4 cm (see Fig. 6).
Considering mass transfer resistance in liquid phase only, equilibrium is also reached, except for toluene and furan. The influence of 50% increase in the liquid level, corresponding 50% increase in mean RT, on the equilibrium is negligible as shown in Fig. 5(b). Compared with other compounds, toluene and furan require longer time to reach equilibrium in the liquid phase due to high infinite dilution activity coefficients, see Tables VI and I. Figure 5(b) shows that mass transfer limited RT is fairly longer with respect to time in model RTD, indicating that the leaving gas is not in equilibrium with the liquid. The results presented in Fig. 6 also indicate low concentration at 4 cm, showing that with the larger bubble size, over 20 cm liquid level is needed to reach 99.9% equilibrium concentration. Calculations are based on plug flow conditions [combining Eqs. (19) and (24)].
The approach equilibrium concentration is a function of approach solute concentration and RTD [see Eq. (21)]. Therefore, in order to study quantitatively if the leaving gas is in equilibrium, both the bubble RTD from CFD modeling and approach solute concentration (Fig. 7) should be considered. The approach solute concentration is calculated by using   Figure 7 shows the time needed to reach equilibrium for each bubble size, which is 0.68 s and 1.03 s for large and small bubbles, respectively. The approach equilibrium concentration of all compounds is shown in Table VIII, supporting the low equilibrium concentration shown in Fig. 5(b). Figure 8 shows the effect of agitation speed on the gas bubble mean residence time. The inert gas flow rate is kept constant at 10 cm 3 /min, and the agitation speed is varied in the range 300-900 rpm. The results show that increasing agitation speed increases bubble mean RT. This trend is sharper for the 1 mm than for the 2 mm bubbles. For smaller bubbles, the influence starts to grow after 700 rpm. Figure 9 shows that the shape of the RTD curve of the bubbles starts to change at the agitation speed of 800 rpm. This indicates that small bubbles start to get trapped behind the stirrer, see Fig. 10. This might lead to the formation of larger bubbles due to the coalescence process.
Gas volume fraction increases with increasing mean RT. One of the main assumptions of the Eulerian DPM in the Fluent software is low gas volume fraction (<12%) in the flow. Gas average volume fraction in the 1000 rpm, constant bubble size of 1 mm and 10 cm 3 /min is 3.7%. Thus, the influence of high   bubble volume fraction, i.e., the bubble-bubble interaction, is not taken into account in the model. Mixing speed strongly affects the mean RT (Figs. 8-10). To reach equilibrium concentration, toluene and furan require higher impeller speeds, especially for the larger bubble size. Therefore, the effect of mixing on the approach equilibrium concentration can be studied by changing the agitation speed (see Table IX). Table IX shows that to reach an equilibrium concentration for toluene and furan, the agitation speed should be well above 800 rpm. All other compounds require much less agitation speed, which is around 300-400 rpm, see Table VIII. Agitation improves the approach equilibrium concentration, although ideal mixing is worse than 300 rpm case (Table VIII). Clearly, increased bubble hold-up (trapping of the bubble) overcomes the danger of by-passing.
At agitation speeds above 800 rpm, the flow region in the dilutor is not bubble-driven anymore, as seen in Figs. 9 and 10. The coalescence process takes place in the stirrer vicinity due to high gas volume fraction. Because of this limitation, the maximum approach equilibrium concentration for toluene and furan in the liquid level of 4 cm, Q = 10 cm 3 /min, and bubble sizes of 1 and 2 mm are 0.948 and 0.986, respectively. The effect of the liquid level is studied by changing the liquid level from 4 cm to 6 cm, where the agitation speed is kept at 800 rpm and the inert gas flow rate at 10 cm 3 /min.
The influence of the changing liquid level from 4 cm to 6 cm on the approach equilibrium concentration is similar to that of increasing agitation speed from 800 rpm to 900 rpm for the 4 cm liquid level (Table X). As there is a clear danger of entrapped bubble coalescing above 800 rpm, it is advisable to increase the liquid level instead of impeller speed in this case.

VIII. CONCLUSIONS
Two-phase CFD was used to model mass transfer and fluid flow in a dilutor apparatus designed to measure infinite dilution activity coefficients. Based on the detailed flow pattern and residence time distribution of bubbles, it was possible to evaluate the approach to equilibrium in a more rigorous manner than with traditional mass transfer models assuming ideal flow patterns. Results showed that the approach to equilibrium is in some cases component dependent: toluene and furan require relatively higher agitation speed or significantly higher liquid levels, compared with the other studied bio-oil compounds, for proper equilibration. Results indicate also an upper limit for the mixing rate, originating from the coalescence process. This limit is around 800 rpm for the studied diluter cell. After this, the mixing element starts to collect gas around it, potentially leading to a coalescence of bubbles with reduced mass transfer area and longer times required for equilibration. CFD as a design tool predicts the influence of the compound of interest and changing the dimension on reaching equilibrium in the dilutor cell. CFD was shown to predict potential by-passing in different operation conditions that traditional mass transfer models are not capable of.

ACKNOWLEDGMENTS
The authors would like to acknowledge the financial support from the Academy of Finland. (ExtraFin-PORLIS project).