CFD Simulations of Two Phase Flow in Asymmetric Rotary Agitated

Nov 26, 2018 - In the present work, a computational fluid dynamics (CFD) methodology has been applied to study the hydrodynamics of two phase flow in ...
0 downloads 0 Views 8MB Size
Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. 2018, 57, 17192−17208

CFD Simulations of Two Phase Flow in Asymmetric Rotary Agitated Columns Raosaheb A. Farakte, Nilesh V. Hendre, and Ashwin W. Patwardhan*

Ind. Eng. Chem. Res. 2018.57:17192-17208. Downloaded from pubs.acs.org by TULANE UNIV on 12/19/18. For personal use only.

Department of Chemical Engineering, Institute of Chemical Technology, Matunga, Mumbai-400019, India ABSTRACT: In the present work, a computational fluid dynamics (CFD) methodology has been applied to study the hydrodynamics of two phase flow in an asymmetric rotating disk contactor (ARDC) and asymmetric rotating impeller column (ARIC). The Euler−Euler model for multiphase flow and mixture k−ε model for turbulence have been used. The effect of different drag models on holdup of the dispersed phase has been investigated. The Kumar and Hartland (KH) drag model has been modified to predict the holdup accurately. The results of simulation have been validated with published experimental data. The average error in the prediction of holdup using a modified KH drag model is within ±11% of the experimental values. Hydrodynamics of ARDC and ARIC have been compared using simulations at different agitator speeds. The fraction of static holdup in ARDC and ARIC has been quantified. The developed CFD model has been used to predict the residence time distribution (RTD) of the dispersed phase. single phase simulations.9,10 These researchers used a 2D axisymmetric approach to reduce computational cost. Fei at al.10 measured the velocity profiles in RDC for a single phase system using Laser Doppler Anemometry (LDA). CFD simulations were performed on 2D axisymmetric geometry of RDC and were found to be in good agreement with LDA measurements. On the basis of the CFD results, a modified design for RDC was proposed to improve mass transfer performance. Drumm and Bart11 carried out CFD simulations of RDC using 2D geometry for single and two phase systems. The effect of different turbulent models on the flow profile was studied. The Reynolds stress model (RSM) and realizable k−ε model were found to produce better results for single and two phase systems, respectively. Yadav and Patwardhan12 performed CFD simulations of PSPC using 2D geometry and Euler−Euler methodology. Experimental results of dispersed phase holdup and height of the accumulated layer were used for validation of simulations. Ghaniyari-Benis et al.13 performed 3D simulations of RDC using a population balance model (PBM) to study drop size distribution. CFD simulations were used to study a system of ionic liquid and a mixture of heptanes and toluene in RDC by Onink et al.14 A modified Schiller−Naumann drag model was used to match the higher holdup observed in experiments at high flux of both the phases. Simulations of an industrial RDC of diameter 0.45 m have also been carried out by Drumm et al.15 The results of single phase simulations were in accordance

1. INTRODUCTION Liquid−liquid extraction (LLE) is one of the most important industrial processes. It is widely used in process industries such as refining, pharmaceuticals, mining, nuclear, paints, and dyes. The process of LLE is performed using various liquid−liquid contactors. The selection of the contactor is mainly based on the type of fluids to be handled and their physical and chemical properties. Mixer-settlers, rotating disk contactors (RDCs), asymmetric rotating disk contactors (ARDCs), pulse sieve plate columns (PSPCs), pulse disk and doughnut columns (PDDCs), etc. are the few widely used LLE units. The intimate contact of two liquid phases is achieved by mechanical or pulsed agitation. RDCs are widely used on account of their simple design and high throughput, but their performance is reduced by high axial mixing. ARDC has advantages of high mass transfer rate and reduced axial mixing in the two phases.1 In ARDC, the mixing of both phases is achieved using rotating disks. The disks are mounted on a shaft, located off-center to the column. The entire column is divided into a series of mixing and settling zones (Figure 1). Both phases pass through these zones in a counter-current direction. Several experimental studies have been carried out to understand the performance of ARDC.2−7 These studies have focused on the holdup, drop size distribution, and axial dispersion. Few researchers have developed empirical correlations for these parameters based on a wide range of experimental data.5,8 A computational fluid dynamics (CFD) study of various extraction columns has been reported in the literature. Table 1 represents the recent CFD simulation studies on liquid−liquid contactors such as RDC, PSPC, PDDC, etc. reported in the literature. Earlier studies on LLE columns were limited to © 2018 American Chemical Society

Received: Revised: Accepted: Published: 17192

June 19, 2018 November 24, 2018 November 26, 2018 November 26, 2018 DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 1. Internals and top view of ARDC.

experimental data on ARDC (36 stages) as reported by Kadam et al.5 A good agreement between the simulated results and experimental data of Kadam et al.5 was reported. In the present work, to overcome the issue of large mesh count and to reduce the computational load, simulations were first carried out for a section of ARDC with five stages. By this way, mesh count could be reduced, and case was setup with relatively less computational cost. The developed simulation methodology was then applied for full scale ARDC (15 stages). CFD simulations can provide a comprehensive understanding of the local hydrodynamics, which in turn aids in design and optimization of these extraction columns. In view of this, the objectives of the present work are as follows: (i) to develop a CFD model to study the hydrodynamics of ARDC and ARIC, (ii) to validate the model using experimental holdup data, and (iii) to predict the flow pattern and Peclet number of the dispersed phase (Ped) using CFD.

with the experimental data. However, improvements in simulation methodology were suggested for the two-phase system. Chen et al.16 carried out CFD simulations coupled with population balance modeling and compared the values of continuous phase Peclet number (Pe) with those obtained from the correlation developed by Strand et al.17 Sen et al.18 studied the axial dispersion in single phase in PSPC using CFD. They analyzed the effect of hole diameter and percent free area on the axial dispersion coefficient. Sen et al.19 reported two phase flow simulations using different drag models to predict holdup. It was found that the Kumar and Hartland20 drag model (KH drag model) predicts the holdup better than other drag models. To further improve the predictions, the KH drag model was modified. Yi et al.21 simulated a nonpulsed disk and doughnut column using a modified drag model similar to Onink et al.14 The modified drag model reduced the error in prediction of continuous phase Peclet number. In a recent publication,7 our research group studied the hydrodynamic and mass transfer performance of an asymmetric rotating impeller column (ARIC), in which 45° pitched blade-disk impellers were used for agitation. They developed generalized correlations for drop size and holdup in terms of power consumption per unit volume and compared the hydrodynamic performance of ARDC and ARIC. The flow pattern in ARDC and ARIC are completely different on account of different agitators. Although there has been significant work done on CFD simulations of the conventional rotating disk contactor (RDC), no research work has been published on CFD simulations of ARDC/ARIC. Because of the rotational symmetry, 2-D axisymmetric simulations of RDC are possible and require less computational time. However, in ARDC, the rotating shaft is located off center, because of which, 2-D or 2-D axisymmetric simulations of ARDC are not possible. Also, periodic boundary conditions can only be applied for single phase simulations and not for simulating two phase flow. The only way is to simulate full 3-D geometry of ARDC. Most of the experimental studies reported in the literature have been carried out on pilot scale ARDC. Simulating geometries of such a large scale would be computationally intensive. To the best of our knowledge, Deshmukh22 has carried out CFD simulations of ARDC previously; however, this is an unpublished work. The author has carried out single phase (water) and two phase (water/ toluene) simulations of ARDC. The geometry used for simulations was three stages, and results were compared with

2. EXPERIMENTAL DETAILS The experimental data reported by Hendre et al.7 were used in the present study. Figures 1 and 2 show the internals of the experimental setup and dimensions of the ARDC and ARIC, respectively. A total of 15 number stages were present in each column. The internal assembly (stators and rotors) of the glass column is made of SS316. For hydrodynamic study, aqueous and organic phases were selected as the continuous and dispersed phases, respectively. The PBD impeller (O.D. 5.2 cm) consists of a central disk of diameter 4 cm and four blades (1.2 × 1.2 × 0.2 cm) pitched at an angle of 45°. The experimental setup was fabricated by KVC Pvt. Ltd., Thane. Two liquid−liquid systems of different physical properties were used for the hydrodynamic study. The details of physical properties are shown in Table 2. System 1 is comprised of water as a continuous phase and heavy normal paraffin as a dispersed phase. In system 2, phosphoric acid is the continuous phase, and organophosphorus solvent is the dispersed phase. The organic phase in system 2 is a mixture of TBP, D2EHPA, and HNP. The hydrodynamic parameters such as drop size and holdup at different rotation speeds were used for validation of simulations. The pulse injection technique was employed in determining the residence time distribution (RTD) of the dispersed phase. The RTD was analyzed using a dispersion model to calculate the Peclet number of dispersed phases (Ped). The Peclet number is a measure of relative importance of convective to diffusive transport phenomena. It is inversely proportional to the axial dispersion coefficient. A 17193

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

17194

PSPC

nonpulsed disk doughnut column

Sen et al.19

Yi et al.21

PSPC

Sen et al.18

RDC

Onink et al.14

RDC

RDC

GhaniyariBenis et al.13

Chen et al.16

sieve plate column and PSPC

Yadav and Patwardhan12

RDC

RDC

Drumm and Bart11

Drumm et al.15

RDC, modified RDC with perforated disk at stator opening

type of column

Fei et al.10

reference

Dt = 0.05 m Hc = 0.1 m DO = 0.002 m Z = 0.5 m Dt = 0.219 m Dd = 0.11 m Ds = 0.152 m Hc = 0.0865 m Dt = 0.06 m Dd = 0.04 m Ds = 0.06 m Hc = 0.032 m Dt = 0.450 m Dd = 0.270 m Ds = 0.315 m Hc = 0.090 m Dt = 0.150 m Dd = 0.090 m Ds = 0.105 m Hc = 0.030 m Dt = 0.075 m Dh = 0.003 m (triangular pitch) Hc = 0.050 m Dt = 0.075 m Dh = 0.003 m (triangular pitch) Hc = 0.050 m Dt = 0.0725 m Ddisk = 0.0624 m Ddoughnut = 0.0369 Hc = 0.019 m

Dd = 0.90 m Ds = 0.105 m Hc = 0.030 m

Dt = 0.29 m Dd = 0.174 m Ds = 0.19 m Hc = 0.078 m Dt = 0.150 m

dimensions

2D, two phase, steady state, RTD realizable k−ε modified drag model

2D, two phase Euler−Euler k−ε turbulent model modified KH drag model

2D, single phase, unsteady state RTD Euler−Euler k−ε model

2D, two phase, steady state, unsteady state RTD Euler−Euler realizable k−ε mixture Schiller−Naumann drag population balance model

2D, single and two phase, steady state, different turbulent models, Euler−Euler, Schiller−Naumann drag model

2D, two phase, steady state Euler−Euler, realizable k−ε, modified schiller−Naumann drag model

3D, two phase, steady state Euler−Euler, k−ω turbulence model, Musig model, Schiller−Naumann drag

Euler−Euler, k−ε turbulent model, Schiller−Naumann drag

steady state Euler−Euler different turbulent models Schiller−Naumann drag 2D, two phase, unsteady state

2D, single and two phase

2D, single phase steady state k−ε turbulent model

simulation methods

Table 1. Literature Survey of CFD Simulations of Liquid−Liquid Contactors

sulfuric acid-alamine 336 and isodecanol

TBP−nitric acid

water

mixture (of water and EG) and nheptane

water

ionic liquid: mixture of n-heptane and toluene

n-butanol-succinic acid−water

toluene−water

mixture (of water and EG)- n-heptane

water

system

drop size from experiments

drop size from experiments

single phase simulations

population balance model

drop size from experiments

drop size from experiments

population balance model

drop size from correlation

drop size from experiments

single phase simulations

input of drop size major findings

modified drag model improves prediction of Pe

good agreement between experimental and simulated results of holdup

modified drag model predicts the holdup in good agreement with experimental data

axial dispersion coefficient can be predicted for single phase flow in PSPC

developed a pseudo-3D model for RDC to combine the benefits of both 2D and 3D models

single phase results are in agreement with experiments; holdup profile needs improvement, whereas velocity profile predicted accurately by simulations

predicted hydrodynamics characteristics of RDC for extraction of toluene

good agreement between experimental holdup with simulation data; average error is below 11%

holdup and height of accumulated layer below dispersed phase is predicted accurately

good agreement between simulated and experimental velocity profiles measured by PIV single phase: RSM is best turbulent model two phase: k−ε -RNG is best turbulent model

good agreement between simulated and experimental velocity profiles measured by LDA; proposed new RDC geometry

Industrial & Engineering Chemistry Research Article

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

modified drag model improves holdup prediction for both ARDC/ ARIC

large value of Pe indicates convection dominated flow, and a small value of Pe indicates diffusive flow.

3. SIMULATION METHODOLOGY For simulation of two-phase flow, the Euler−Euler two fluid model was adopted. As per the Euler−Euler method, the two phases were considered as interpenetrating continua. The concept of phasic volume fraction (αi) was utilized to specify the volume fraction of the respective phase. An interphase exchange coefficient (Kij) was used to model the exchange of momentum between phases. Continuity and momentum eqs (eqs 1 and 2, respectively) were solved for both continuous and dispersed phases. ∂ (αiρi ) + ∇·(αiρi Ui ) = 0 ∂t

(1)

drop size from experiments

∂ (αiρi Ui ) + ∇·(αiρi UU i i) ∂t

n

= −αi∇p − αi∇·τi + αiρi g +

∑ R ij i=1

HNP-water, solventphosphoric acid

+ (Fi + Flift, i + Fvm, i)

(2)

3D two-phase, steady state, RTD standard k−ε modified drag model

The pressure (p) was shared by both the phases. The lift forces (Flift,i) and virtual mass forces (Fvm,i) can be neglected for liquid−liquid system.11,12 The term τ̿i is the stress tensor, and it is composed of two components, i.e., viscous stress τ̿i,laminar and turbulent stress τ̿i,turbulent. These terms are expressed as follows: τi̿ ,laminar = μi (∇Ui + ∇ UiT )

(3)

τi̿ ,turbulent = μt , m (∇Ui + ∇ UiT )

(4)

The closure of the above equations was achieved by the definition of interphase exchange force (∑R̅ ij) as

∑ R ij = ∑ K ij(Ui − Uj)

(5)

where Kij was defined as K ij =

3αiαjρj C D|Ui − Uj| 4d p

(6)

Dt = 0.103 m Dd = 0.04, 0.052 m Hc = 0.05 m prediction of static holdup RTD simulation of dispersed phase

where CD is the drag coefficient and dp is the drop diameter of the dispersed phase. The universal drag model and KH drag model were used in the present study. The universal drag model has been previously used for calculating the drag coefficient in a gas−liquid flow.23 As per the universal drag model, for a drop Reynolds number (Re) between 1 and 1000, CD was defined as follows:

ARDC/ARIC

CD =

24 (1 + 0.1Re 0.75) Re

(7)

In the case of Re > 1000, the following expression is used for the calculation of drag coefficient (CD). 18/7 y2 ij σ yz zz 2 ijj d p yzzjij 1 + 17.67(1 − αi) jj zz z zzjjj = λ ; C D = jj jj z z z jj g Δρ zzz 3 jk λRT z{k 18.67(1 − αi)3 ij { { k

present work

type of column reference

Table 1. continued

dimensions

simulation methods

system

input of drop size

major findings

Industrial & Engineering Chemistry Research

0.5

(8)

where λ is the Rayleigh−Taylor instability wavelength. 17195

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 2. Internals and top view of ARIC. n

Table 2. Physical Properties of Liquid−Liquid Systems Used for the Study

μm =

∑ αiμi i=1

phase system 1 system 2

water HNP phosphoric acid organophosphorus solvent

density (kg/ m3)

density difference (kg/m3)

998.2 770 1230

248.2

860

390

viscosity (Pa s) 0.001 0.0024 0.002

interfacial tension (N/m)

n

ρm =

∑ αiρi

÷◊÷ n ∑i = 1 αiρi Ui ÷÷÷÷◊ Um = n ∑i = 1 αiρi

0.035

i=1

0.02

0.0036

For evaluation of Re, effective viscosity based on volume fraction of dispersed phase was utilized as per eqs 9 and 10. ÷◊÷ ÷◊ ÷ ρj d p|Ui − Uj| Re = μeff (9)

k2 ε

(17)

There are many aspects such as turbulence modeling and droplet interactions that can affect the accuracy of predictions of CFD models as discussed above. In the two fluid model, monodispersed drop size is assumed throughout the entire computational domain. In reality, there exists a drop size distribution in LLE columns because of several mechanisms like drop breakage and coalescence. Hence, assumption of monodispersed drop size could lead to inaccuracy in simulations. The assumption of monodispersed drop size in CFD simulations provides satisfactory results at less computational power. A realistic approach can be made by coupling CFD with population balance models, which accounts for breakage and coalescence of droplets, but it requires huge computational power in case of 3-D simulations. The simulation results are also affected by the choice of turbulence models. As listed in Table 1, there are many models of turbulence available, which have been used for simulations of LLE columns (Table 1). However, it is not feasible to address all issues together that could otherwise lead to increased computational load. A reasonable solution can be obtained by combining all these uncertainties in one or two terms of governing equations, and the error between predicted and experimental results can be reduced by modifying these terms. There are few reports in which this approach has been followed for the two phase CFD simulations, wherein the drag model has been modified to approximate the predictions of simulations with experimental results.19,24,25

(10)

24 yz i zz(1 + Aαi B) C D = jjj0.53 + (11) Re { k The value of empirical constants A and B are 4.56 and 0.73, respectively. The k−ε mixture model was used for modeling the turbulence. The following equations of k and ε were solved for the mixture instead of all phases. ÷÷÷÷◊ ∂ (ρm k) + ∇·(ρm k Um) ∂t ijji μt , m zy yz zz∇k zzz + G − ρ ε = ∇·jjjjjjjμm + k ,m zz z m j σ k (12) { { kk

The general form of the KH drag model is shown in eq 11:

÷÷÷÷◊ ∂ (ρ ε) + ∇·(ρm ε Um) ∂t m ijij μt , m yz yz zz∇εzzz + ε (C G − C ρ ε) = ∇·jjjjjjjμm + 2ε m zz z k 1ε k , m j σ ε { kk {

(16)

The expression for generation of turbulent kinetic energy due to mean velocity (Gk,m) is shown in eq 18, ÷÷÷÷◊ ÷÷÷÷◊ ÷÷÷÷◊ Gk , m = μt, m (∇ Um + (∇ Um)T ):∇ Um (18)

μj (1 − αi)2.5

(15)

The turbulent viscosity of the mixture is given by eq 17: μt, m = ρm Cμ

μeff =

(14)

(13)

The viscosity, density, and velocity of the mixture is calculated as per the following equations: 17196

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

flow.27−29 This phenomenon has been attributed to the interaction between droplets and turbulent flow field. The higher value of A in the drag model for the pulsed column than for the spray columns could be related to the droplet interactions and extent of turbulence prevailing in the two columns due to a difference in geometry and mode of agitation. In addition, a literature study shows that the prediction of drag coefficient for a particular system has always defied researchers. In the past, several empirical drag models were developed for prediction of CD for different liquid−liquid systems. Typical plots of CD vs Re are available in the literature for various liquid−liquid systems.30 It has been observed that in the region of high Reynolds number, the predicted values of CD at a similar Reynolds number are quite different. Thus, there is always some uncertainty in predicting CD for a specific system. All the simulations in the present work were performed under steady state conditions. The CFD code Fluent 17.0 was used for carrying out simulations. For two phase simulations, use of velocity inlet boundary conditions for the inlets and the outlets of the LLE column has been reported.11 The velocity inlet with negative velocity is used at the outlets of LLE column to ensure the respective phase is moving out of the solution domain. The simulations with two velocity inlets and two pressure outlets were performed. However, this set of boundary conditions resulted in unrealistic flow rates at the outlets. The flow rate through the aqueous outlet was observed to be much higher than the inlet flow rate, and the flow rate at the outlet of the organic phase was observed to be lower than at the inlet. Thus, this set of boundary conditions yielded incorrect results. At a steady state, inlet and outlet flow rates for respective phases should be equal. Thus, boundary conditions for both outlets were chosen as mass flows with equal flow rates as the inlet with direction normal pointing outside of the solution domain. Rotation of the disk was modeled using frame motion of the cylindrical domain around all disks. Moving wall BC was given for the wall of disks and shaft. Gauge pressure at the organic outlet was fixed to zero as the organic phase exits the column at atmospheric pressure. To model static pressure in the column, specified operating density was fixed to 0 kg/m3. This enabled the calculation of static pressure in the column and U-tube, which helped in maintaining the interface. In liquid−liquid extraction columns, a steady state is assumed when the interface position in the column does not change with time. This condition is reached when the static pressure at the bottom of the column is balanced by the pressure in the U-loop. In the beginning, when the two phase simulations were being carried out, it was observed that the interface position always changed, and a converged solution was not obtained because of a pressure imbalance despite having provided a U-loop. This issue was solved when the specified operating density option was activated and its value was set to 0 kg/m3. This resulted in accurate calculation of static pressure, and the steady interface position in the column was maintained. Numerical solutions were obtained using a phase-coupled SIMPLE solution scheme along with a second order discretization methodology for momentum, k and ε. The volume fraction equation was discretized using the QUICK scheme. The holdup of the dispersed phase in the column was monitored, and a converged solution was obtained after attaining a steady value of hold-up. The dispersed phase was assumed as monodispersed, and an experimentally measured Sauter mean diameter was provided

The drag force is an important factor, which decides the rise or settling velocity of dispersed particles/droplets in a continuous medium through which it is moving. The Schiller−Naumann drag model is commonly used to simulate two phase flow in LLE columns. This drag model has been derived for a solid sphere immersed in an inifinite pool of liquid. Thus, this model is applicable for the case of very lean dispersions. In addition, modeling of the drag coefficient in a liquid−liquid dispersion is more complex due to lower interfacial tension and presence of impurities at the interface. The large droplets of droplets wobble during the flow, thus changing the form of drag acting on them. All these phenomena are not covered in the Schiller−Naumann drag model. Hence this drag model has not been used for the simulations of the present case. The drag models used in the present work consider the effect of the volume fraction of the dispersed phase on the drag coefficient (CD). In the universal drag model, the effect of the volume fraction has been incorporated in terms of effective viscosity of the continuous phase (eq 10). The region with a higher volume fraction will have higher viscosities and lower Re. Thus, as per eq 7, CD will be higher in the region of high volume fraction, which could lead to better prediction of the holdup in the column. The expression of the KH drag model (eq 11) shows that CD is proportional to the holdup of the dispersed phase. The KH drag model was originally developed for estimation of holdup and slip velocity for liquid−liquid dispersion systems in spray columns. The first term in this equation represents the CD value for a single particle system. The second factor (1 + A·αiB) in eq 11 accounts for the effect of particle interaction and wall effect on drag coefficient in multiparticle systems. The term “particle interaction” refers to the hindrance effect caused by presence of neighboring particles on the settling velocity of individual particles. The term “wall effect” refers to retardation force exerted by the container walls. This reduces the terminal velocity of particles. The term 1 + A·αiB was used by Famularo and Happel26 as a correction factor in the calculation of settling rates by Stoke’s model to account for the effect of resistance on settling in dilute suspensions. The constant “A” is termed as a concentration coefficient, and its value was estimated to be 1.3 for random suspension in the vicinity of a settling particle. The value of B equals 1/3 in the expression 1 + A·αiB is obtained by mathematical rearrangement. Barnea and Mizrahi25 estimated the value of A to be 1.0 for a solid/liquid multiparticle system. Kumar-Hartland20 evaluated the values of A and B for measurement of slip velocity in spray columns for liquid−liquid systems. The best fit values of constants A and B in the drag coefficient equation were found to be 4.56 and 0.73. For liquid−liquid dispersions, the extent of hindrance to the settling of dispersed phase droplets also depends on the physical properties of the system. Sen et al.19 modified the values of A and B for more accurate prediction of holdup in pulsed sieve plate columns. With an increase in A and decrease in the value of B, an increase in holdup was reported. The holdup was found to be less sensitive to B as compared to A. In spray columns, the magnitude of droplet size and holdup are primarily decided by sparger design and velocities of continuous and dispersed phases. In pulsed columns, the magnitude of turbulence is higher as compared to spray columns because of extra turbulence induced by pulsing action. In the literature, it has been reported that there is an increase in drag coefficient with an increase in turbulence for two phase 17197

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

used for the tracer equation. A total of 2 mL of tracer was introduced into the organic inlet of the computational domain. Courant number for the initial flow time (up to 5 s) was maintained below 0.8. The time step was increased periodically without deviating residual criteria. Concentration of the tracer was monitored at the organic outlet. Simulations were carried out for a time period close to thrice the experimental mean residence time of the tracer. The Peclet number values are determined from the two moments of the RTD curve, mean residence time, and standard deviation. These two moments are evaluated at different operating conditions from the outlet concentration profile of the tracer, and the Peclet number was calculated as follows:31

as input for simulations. The Sauter mean drop size at respective rotation speeds in ARDC and ARIC has been reported in our previous publication.7 These drop sizes were provided as input for the simulations of respective rotation speeds, and the values are given in Table 3. Using this case Table 3. Drop Size for System 1 and 2 at Various Rotation Speeds (Vc = Vd = 1 mm/s) column

system

RPM

drop size, mm

ARDC

2

200 400 600 800 400 800 100 150 200 250 150 250

3.03 2.58 1.36 0.79 3.0 0.99 2.34 1.6 1.0 0.92 2.2 0.6

1 ARIC

2

1

σ2 2Pe + 8 = 2 2 tm Pe + 4Pe + 4

(20)

4. COMPUTATIONAL MODEL AND MESH The flow profile in ARDC/ARIC is not axisymmetric, due to the asymmetric location of the shaft and rotors. It is not possible to simulate a 2D or pseudo 3D model of ARDC/ ARIC. The grid independency study and case setup were performed on five-stage ARDC instead of a complete 15-stage ARDC in order to reduce computational time. A clearance of 0.5 mm was maintained between the stator and inner surface of the glass column in simulations, in accordance with the experimental setup. The modeled geometry is shown in Figure 3. A cylindrical domain was created around all the disks to simulate the rotation of disks as per the moving reference frame (MRF) model. In the modeled geometry, the size of clearance between the wall and stator is around 0.5 mm, which makes it difficult to mesh the column geometry with hexahedral cells. The rotor domain is a 3-D cylinder and can be easily meshed with hexahedral elements. In order to reduce mesh count and computational cost, the rotor domain was meshed with hexahedral cells while the remaining part was meshed with tetrahedral elements. The diameter and length of the rotor domain was 5.4 and 2.5 cm, respectively.

setup and five-stage ARDC geometry, simulations were carried out to match the experimental holdup profile. The simulation methodology developed for 5-stage geometry was further applied to 15-stage geometry. The developed flow profile in ARDC (15 compartments) was used for the simulation of RTD of the dispersed phase. The RTD simulations were carried out with the help of the species transport model. Following the simulations of two phase flow, a tracer transport eq (eq 19) was solved for tracer material. ÷◊ ÷ ∂ρYi ⃗ + ∇·(ρUY j i ) = −∇· Ji (19) ∂t A tracer material having the same properties as the dispersed phase was created. The tracer equation was solved using a second order upwind scheme. A residual criterion of 10−5 was

Figure 3. Geometrical model and mesh of five-stage ARDC. 17198

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research A grid independency test was performed on the model geometry consisting of five stages. Three grid sizes were used with cell counts of 1 444 119, 1 698 086, and 1 923 439, respectively. These grids were called low, medium, and high cell count grids. The axial velocity profile of the continuous phase was monitored along a line in the third stage of the column, at a distance of 1.2 cm below the disk. Simulations were performed on these grids for a rotation speed of 500 rpm. System 2 was selected for simulations, and superficial velocities of continuous (Vc) and dispersed phases (Vd) were kept at 1 mm/s. The superficial velocities Vc and Vd have been defined based on column cross section calculated using column inner diameter. Drag force was modeled using the universal drag model. Figure 4 shows the aqueous phase velocity profile in the mixing stage. Higher velocities can be observed near the

Figure 5. Holdup profile for system 2. (a) Universal drag model. (b) Original Kumar Hartland drag model: (Vc = Vd = 1 mm/s).

The effect of different turbulence models like standard k−ε, reliazable k−ε, RNG k−ε, and the k−ω model on the dispersed phase holdup was investigated at 400 rpm in ARDC. It has been found that there is minimal change in the prediction of holdup using these models. Also, the model of Burns et al.32 was introduced to account for the effect of turbulent dispersion force. The turbulent dispersion force term appears in the interphase momentum transport equation. The influence of turbulent dispersion on the simulated results was found to be quantitively insignificant. It also has been previously shown that for simulations of liquid−liquid dispersions, the turbulent dispersion term has a negligible contribution to the overall interphase momentum transport.11,33 Therefore, for further simulations, turbulent dispersion force was not considered. In order to reduce the error between experimental and simulated results, it was thought to modify the two abovementioned drag models that would increase the value of the drag coefficient and increase the resistance to rising droplets which would lead to higher holdup in the column. The drag coefficient in liquid−liquid dispersion is affected by the presence of impurities, surface active contaminants, and rigidity of the interface and walls of container. Traces of impurities or contaminants and drop deformations diminish internal circulations, which significantly increase the drag on dispersed phase particles. The effect of all these parameters on the drag coefficient is not considered in the equation of drag models. In order to select a suitable drag model for modification, the dependency of the drag coefficient on the volume fraction was analyzed for these two drag models. The volume averaged slip velocities were calculated at 400 and 800 rpm from simulations performed using the universal drag

Figure 4. Axial velocity of continuous phase (system 2) for three different grids (Vc = Vd = 1 mm/s).

rotating shaft. The velocities close to the shaft are different for a mesh with low grid count than that for high and medium grid counts. The velocity profiles of medium and high cell count grids are identical. Thus, a grid with a medium cell count was used for further simulations. On the basis of the grid independency test, a 1.5 mm cell size was used in the cylindrical domain around the disk, whereas a cell size of 3 mm was used in the remaining part of the model. The holdup values obtained from simulations were 7.91%, 7.87%, and 7.78% for low, medium, and high cell count grids, respectively. The holdup values for all three grids are almost comparable. The generated mesh is shown in Figure 3. The 15-stage geometry of ARDC and ARIC were meshed with the same cell sizes as those of the fie-stage geometry.

5. RESULTS AND DISCUSSION 5.1. Simulations of Five-Stage ARDC. Figure 5 depicts the comparison of simulation and experimental variation of holdup with agitator speed in the case of ARDC for system 2. Figure 5a and b shows the simulated holdup values using the universal drag model and KH drag model, respectively. Below 500 rpm, the holdup values obtained from simulations are in good agreement with the experimental results in the case of both the drag models. Above 500 rpm, simulation shows very low holdup as compared to experimental results. 17199

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

change in holdup was less sensitive to changes in the value of B; it was decided to change the value of A only. Simulations were carried out for the values optimized by Sen et al.,19 and the predicted results were observed to be far from the experimental one. The reason can be inaccurate prediction of drag force in ARDC, which leads to error in simulated results. As has been discussed in the previous section, the drag coefficient is significantly affected by the extent of turbulence prevailing in the two phase system. In the case of mechanically agitated extraction columns like ARDC, the degree of turbulence can be considered significantly higher than that in spray and a pulsed column on account of turbulence generated by rotating impellers in the former. Also, the flow pattern in ARDC is quite different from that observed in the pulsed column. In the pulsed column, the two phase flow is virtually counter-current, whereas in ARDC, the two phases follow a helical path in the column because of the presence of horizontal stators and baffle arrangement. Thus, the flow behavior of pulsed columns and ARDC are thoroughly different, which significantly affects the hydrodynamics. Therefore, the optimized values of A and B in eq 11 of the KH drag model of Sen et al.19 needs to be modified for adequate prediction of the drag coefficient in ARDC. The constant A in the KH drag model (eq 11) was increased from its original value of 4.56. Simulations with the modified KH drag model were carried out for system 2 at 800 rpm in ARDC. Figure 7 shows the variation in the values of holdup

model. These values of slip velocities along with drop size at respective rotation speed were used to calculate Re, which in turn was used to calculate CD from the universal and KH drag models (eqs 7 and 11). Figure 6a and b depict the variation of CD with respect to the holdup of the dispersed phase at 400 and 800 rpm,

Figure 6. Relation between holdup and CD for different drag models (a) 400 rpm, (b) 800 rpm (Vc = Vd = 1 mm/s).

respectively. An increase in the value of CD is observed with an increase in holdup values for both drag models. However, in the case of the KH drag model, the increase in CD is much higher in comparison with the universal drag model. For simplicity of calculations, a constant slip velocity was assumed for all holdup values, at a particular agitation speed. In reality, an increase in holdup lowers the slip velocity,20 which leads to a further increase in CD. These results show that CD predicted by the KH drag model is more sensitive to the holdup of the dispersed phase. 5.2. Modification in KH Drag Model. The KH drag model was developed by Kumar and Hartland20 for the gravity settling of liquid−liquid dispersions in spray columns by analyzing experimental data for several liquid−liquid systems. The general form of a KH drag model has been given in eq 11. The empirical constants in this equation, A and B, are fitted parameters, and their values are 4.56 and 0.73, respectively. This drag model was modified by Sen et al.19 to bring the predictions of simulations closer to the experimental holdup values. As can be seen from eq 8, to increase the holdup, the value of B needs to be decreased. However, the modification in B is limited to positive values. Selecting a negative value of B would show the effect of lowering CD with an increase in holdup. Furthermore, Sen et al.19 have also shown that a

Figure 7. Effect of KH drag model constant (A) on holdup; N = 800 rpm (drop diameter = 0.79 mm, Vc = Vd = 1 mm/s).

with an increase in the value of A. A steady rise in holdup is observed with an increase in value of A up to 19. In this region of moderate slope, the maximum holdup is less than 20%, much lower than the experimental value. With a further increase in A, a rapid rise in the holdup is observed. At A equal to 19.4, the holdup is 28%. The rapid increase in holdup at a value of A greater than 19 can be attributed to the synergistic effect of an increase in holdup and decrease in slip velocity. Simulations for A > 19.4 did not attain a steady state. The holdup was observed to increase with progress in simulation. The CD value is a function of holdup, slip velocity, and parameter A. At higher values of A (>19.4), the calculated value of CD based on holdup and slip velocity is such that it further escalates holdup and reduces slip velocity. It was observed that even after a large number of iterations, the CD value did not attain a steady value, causing a continuous increase in dispersed phase holdup. This causes increasing holdup for the value of A above 19.4. As the value of simulated 17200

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research holdup is close to experimental holdup at A equal to 19.4, further simulations were carried out using the modified KH drag model with a value of A equal to 19.4. 5.3. Simulations of ARDC and ARIC Using Modified KH Drag Model. The prediction of dispersed phase holdup at various RPMs using the modified KH drag model for system 2 is shown in Figure 8a. Initially, simulations were performed in

Table 4. Drop Size for System 2 at Various Rotation Speeds in ARDC column

system

Vc

Vd

RPM

drop size, mm

ARDC

2

3 3 3 1 1 1

1 1 1 3 3 3

200 400 600 200 400 600

2.84 2.4 1.0 2.85 2.37 0.99

Figure 9. Holdup vs RPM for different Vc and Vd in ARDC; system 2.

holdup at different Vc and Vd values in ARDC at 400 rpm. It can be seen that effect of Vd on holdup is much pronounced than that of Vc, as observed in experiments. The average deviation in prediction of holdup for the effect of Vc and Vd is 20% and 33%, respectively. The simulation methodology developed for ARDC was checked for ARIC. The effect of drag model modification is shown in Figure 10 for ARIC. Figure 10a depicts the variation in holdup with RPM for system 1, for the original and modified drag models. Both modified and original drag models predicted similar hold-up values. Figure 10b shows a comparison of holdup for system 2 between experimental and simulation results at four different impeller speeds, i.e., 100, 150, 200, and 250 rpm. At lower speeds, holdups predicted using both original and modified drag models are close to experimental values. At 250 rpm, the original drag model underpredicts the holdup, whereas the modified drag model shows a holdup very close to the experimental value. This shows that the simulation methodology developed for ARDC is also applicable to ARIC geometry. The average error in the prediction of dispersed phase holdup in ARDC and ARIC using CFD simulations is ±11% for both systems. 5.4. Comparison of Hydrodynamics of ARDC and ARIC. In this section, dispersed phase volume fraction and axial velocities of both phases in ARDC and ARIC have been compared near and above the critical rotor speed (NCR). Hendre et al.7 have reported the values of NCR for system 2 in ARDC and ARIC as 427 and 127 rpm, respectively. Therefore, the above-mentioned parameters have been compared at 150 rpm of ARDC and 400 rpm of ARIC, respectively. Furthermore, these parameters have been compared above NCR at 800 rpm of ARDC and 250 rpm of ARIC, respectively. In addition, the power per unit volume (P/V) for these cases can be found in our previous publication.7 The value of P/V at 400 rpm of ARDC using system 2 is 34 W/m3, as compared to

Figure 8. Holdup prediction using modified KH drag model. (a) System 2. (b) System 1 (Vc = Vd = 1 mm/s).

the five-stage geometry to check the performance of the modified drag model at all RPMs. It can be seen that the simulated values of holdup are in good agreement with the experimental data. Similar holdup profiles were obtained for 15-stage geometry using the modified drag model. The maximum error in the prediction of dispersed phase holdup for five-stage and 15-stage geometry are 20.7% and 21.2% at 600 rpm. The 15-stage model of ARDC was further tested for system 1. The results of simulated and experimental values of holdup are shown in Figure 8b. The maximum error in the prediction of dispersed phase holdup is 21.7%. Thus, the modified drag model is suitable for predicting the hydrodynamics of ARDC at all agitation speeds up to 800 rpm. The effect of modified A (from 4.56 to 19.4) on holdup at lower rotational speed can be assessed by comparing Figure 4b and Figure 7. A maximum 12% increase in holdup is observed at 500 rpm. Thus, there is no significant change in holdup at lower rotational speed when the value of A is increased. Thus, the constant A is well within the region of smaller slope of Figure 7 if drawn for respective rotational speeds. Simulations were carried out for ARDC to investigate effect of Vc and Vd on the dispersed phase in ARDC. The corresponding values of the Sauter mean diameter given as input to simulations have been added in Table 4. Figure 9 depicts the effect of RPM on 17201

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 10. Comparison of simulated and experimental holdup profiles. (a) System 1 and (b) system 2 in ARIC (Vc = Vd = 1 mm/s).

Figure 11. Contours of the volume fraction of the dispersed phase for system 2 (a) 400 rpm ARDC, (b) 800 rpm ARDC, (c) 150 rpm ARIC, (d) 250 rpm ARIC (Vc = Vd = 1 mm/s).

48 W/m3 at 150 rpm in ARIC. Similarly, P/V at 800 rpm of ARDC is 165 W/m3, as compared to 226 W/m3 at 250 rpm in ARIC. The P/V values near NCR in both columns are comparable, whereas above NCR, the difference in P/V is modest. However, a comparison can be done to understand the variation in two phase flow above and below NCR for ARDC and ARIC. 5.4.1. Contours of Volume Fraction of Dispersed Phase. Contours of the volume fraction of the dispersed phase (system 2) on a central vertical plane of the column are shown in Figure 11. The vertical plane is chosen such that it passes through both mixing and settling zones. The mixing and settling zones are marked in the figure for better visualization. At low rotation speeds, it is observed that a major fraction of the dispersed phase is present below the stators in mixing and settling zones in both ARDC and ARIC (Figure 11a and c). Volume fraction in this region is close to 1, which indicates accumulation of dispersed phase droplets below the stators. The holdups in the remaining regions of the column are extremely low, indicating flow of the dispersed phase from one compartment to another, without getting completely dispersed

in the mixing region. The accumulated holdup below the stator in the mixing zone is more in ARIC than that in ARDC because of the downflow nature of the PBD impeller. This is reflected in the value of holdup as high as 10% in ARIC as compared to 6% in ARDC, at 150 and 400 rpm, respectively. At lower agitation speeds, volume fractions in the settling region are similar in both columns. A major part of the holdup in this region is located below the small stator. A fraction of the dispersed phase can be observed above the clearance between stators and the column wall (Figure 11a and c) due to the flow of the dispersed phase through the clearances. With an increase in agitation speed, the volume fraction throughout the column is increased, in both the columns (Figure 11b and d). The fraction of the dispersed phase accumulated below the stators is observed to decrease with an increase in agitation speed. This indicates that the dispersed phase is evenly distributed in the mixing region at higher rotation speeds in ARDC and ARIC. 5.4.2. Contours of Axial Velocities of Continuous Phase. To understand the flow behavior of continuous phase in ARDC and ARIC, contours of axial velocities on the central 17202

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 12. Contours of axial velocity (m/s) of continuous phase for system 2 (a) 400 rpm ARDC, (b) 800 rpm ARDC, (c) 150 rpm ARIC, (d) 250 rpm ARIC (Vc = Vd = 1 mm/s).

Figure 13. Contours of axial velocity (m/s) of the dispersed phase overlapped with velocity vectors for system 2 (a) 400 rpm ARDC, (b) 800 rpm ARDC, (c) 150 rpm ARIC, (d) 250 rpm ARIC (Vc = Vd = 1 mm/s).

vertical plane are analyzed (Figure 12). Figure 12a and c show the axial velocities in ARDC and ARIC at 400 and 150 rpm, respectively. The PBD impeller produces downward velocities near the blade, while the radial flow is generated by a disk. The continuous phase velocity near the column wall in the settling zone shows axially upward velocities due to the flow of the dispersed phase through the clearance. At 400 rpm, the presence of upward velocity below the disk and the downward velocity near the wall and vertical baffle indicates the presence of recirculation zones. Similarly, the recirculation zones are also observed at 150 rpm in ARIC. The radial flow generated by the disks splits into upward and downward flow near the vertical baffle and column wall in the mixing zone. The magnitude of upward velocities in the region near the vertical baffle and column wall in ARIC is higher than that observed in ARDC. Figure 12b and d shows the axial velocities at 800 and 250 rpm in ARDC and ARIC, respectively. Strong recirculation

zones are observed in the mixing region of both the columns. The magnitude of axial velocities in ARIC is higher than that in ARDC. At 250 rpm, axially downward velocities generated by the PBD impeller cover the entire length of the mixing zone. At 150 rpm, axial velocity in the region below the big stator is 5 times less than the maximum axial velocity. This reduction in axial velocities is due to the presence of an accumulated layer of dispersed phase below the big stator. Thus, it can be said that the PBD impeller provides higher agitation in the mixing zone at 250 rpm than at 150 rpm, leading to increased breakage of droplets. 5.4.3. Contours of Axial Velocities of Dispersed Phase. Contours of axial velocity overlapped with velocity vectors of the dispersed phase for system 2 are depicted in Figure 13. Figure 13a shows the axial velocity of the dispersed phase at 400 rpm in ARDC. The axial velocity profile of the dispersed phase is completely different than that of the continuous phase at 400 rpm. The dispersed phase has axially upward velocities 17203

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 14. RTD of the dispersed phase (system 2): (a) 400 rpm (ARDC), (b) 600 rpm (ARDC), (c) 150 rpm (ARIC), (d) 250 rpm ARIC (Vc = Vd = 1 mm/s).

the displacement method, all inlet and outlet valves are closed after attaining a steady state in the column. Rotation is continued for some time to avoid accumulation of the dispersed phase below stators and disks. The fall in interface level from a steady state and after shutdown provides an estimate of holdup in the column. In other methods, for evaluating static holdup, all the inlet and outlet valves along with rotation are shut off. This causes only dynamic holdup to move out of the active section of the column. The fall in interface corresponds to the dynamic holdup. However, in the case of ARDC, alternate mixing and settling zones as well as the presence of vertical baffle damps the droplet velocity after rotation is stopped. This results in accumulation of a large part of the dynamic holdup below the stators. Thus, the experimental procedure to measure static holdup can lead to erroneous results. The high volume fraction of the dispersed phase below the stators at 400 and 150 rpm has very low axial velocities, indicating that these are the regions of static holdup. Such regions of high volume fraction and low velocities are less at 800 and 250 rpm in ARDC and ARIC, respectively. An attempt has been made in the present work to quantify the static holdup. The region having a volume fraction above 70% has been considered as the static holdup region. The values of static holdup at 400 rpm (ARDC) and 150 rpm (ARIC) are 57% and 65%, respectively, of the total holdup. Thus, a large fraction of total holdup exists as static holdup below NCR. At 800 rpm, the static holdup in ARDC is reduced to 30% of the total holdup, whereas static holdup in ARIC is only 20% of the total holdup at 250 rpm. Although the contribution of regions of static holdup to mass transfer is less, it is helpful for the coalescence of droplets. The coalesced droplets further break down in the following mixing stage. This process of coalescence and breakage aids in surface renewal, which in

close to 0.15 m/s near the tip of the disk and above clearances. At 400 rpm, a higher drop size is observed, resulting in buoyancy driven flow of the dispersed phase. Similarly, the flow of the dispersed phase is driven by buoyancy of the drops at 150 rpm in ARIC (Figure 13c). Despite the strong downflow generated by the impeller, the dispersed phase shows axially upward flow in most of the regions. Only a small region above the impeller shows downward velocities. The region below the big stators show axial velocities close to zero at both 150 and 400 rpm in ARIC and ARDC, respectively, because of the accumulated dispersed phase in these regions. At 800 rpm in ARDC (Figure 13b), drop size is small, leading to low buoyancy force. Therefore, downward velocities of the dispersed phase are observed at 800 rpm. In addition, the flow of the dispersed phase through clearance is negligible in both mixing and settling zones. At 250 rpm in ARIC (Figure 13d), the flow pattern of the dispersed phase shows a strong downward flow near the impeller. The flow through clearance is also much less. The contours of axial velocity of the dispersed and continuous phases at 250 rpm in ARIC and 800 rpm in ARDC portray similar flow patterns in the mixing region. The results of axial velocities of the dispersed phase are important to understanding the flow behavior near and above the critical rotor speed (NCR). Below NCR, the flow of the dispersed phase is buoyancy driven, whereas with an increase in agitator speed beyond NCR, the dispersed phase starts following a similar velocity pattern to that of the continuous phase. This is beneficial as it causes better mixing of the dispersed phase with the continuous phase and leads to an increase in performance of extraction columns. 5.5. Static Holdup. The experimental procedure for determining static holdup is like the displacement method used for measurement of overall holdup. This method has been previously used for quantification of static holdup in RDC.34 In 17204

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

Figure 15. Tracer flow vs time at different locations in the column. (a) Various locations where tracer flow was monitored, (b) 400 rpm, (c) 800 rpm (inset figures: RTD at inlet of first stage, tracer flow (kg/s × 105) vs time (s)).

length of the column with respect to time (Figure 15a). At 400 rpm (Figure 15a), the RTD shows an initial peak at the inlet of first stage within 1 s after injection. This peak is followed by the nearly constant flow rate of the tracer from 1.5 to 2 s. The later plateau is observed because of the dispersion of the fraction of the tracer in the region between the sparger and first-stage inlet. At 400 rpm, the tracer injected in the column enters the first stage within a period of 4 s. However, at 800 rpm (Figure 15b), ∼25 s are required for the total amount of tracer to enter the first stage. As compared to 400 rpm, the tracer peak at 800 rpm is much lower. At 800 rpm, the region between the sparger and the first stage behaves as a highly backmixed region as the RTD is similar to that of a CSTR. At high agitation speeds, drops of smaller size are present in the system, having low buoyancy force and following the flow pattern of the continuous phase leading to backmixing. The tracer flow was monitored at a plane in stage 10. A split peak is observed at stage 10, due to the bypass of the dispersed phase through the clearance present between stators and the wall of the column. A similar split peak is not observed at the outlet of the last stage of the column. The maximum flow rate at the tracer outlet is observed to be lower than the flow rate at the outlet of the last stage. The existence of an accumulated organic phase between the last stage and tracer outlet acts like a CSTR. Thus, the entire column can be modeled as tanks in series followed by a CSTR having a volume of the accumulated organic phase. The increase in mean residence time at 800 rpm as compared to 400 rpm can be attributed to the decrease in rise velocities of drops on account of small size and higher holdup. In addition, RTD at 800 rpm at stage 10 and the last stage

turn enhances the mass transfer coefficient, leading to improved extraction performance.30 5.6. RTD of Dispersed Phase. Simulations were carried out to study the RTD of the dispersed phase of system 2. The drag model used for simulations was the modified KH drag model. The results of simulations were compared with the experimental exit age distribution (Figure 14) in terms of RTD functions (E(t) vs t). The RTD simulations were also carried out using the original drag model. The effect of modification in the drag model is more pronounced at higher RPMs. RTD results are compared for 250 rpm in ARIC. Figure 14d shows the comparison of E(t) vs t for 250 rpm in ARIC using the original and modified drag models with experimental results. The modified drag model predicts the E(t) curve satisfactorily as compared to the original drag model. The original drag model predicts a lower holdup in the column, which results in lower residence time. The other RTD simulations have been carried out using the modified KH drag model. It can be seen that simulations predict the RTD curve satisfactorily with accurate peak time. The coefficient of determination (R2) for the simulated RTD profile was compared with the experimental data. The value of R2 at 400 and 600 rpm in ARDC is 0.87 and 0.90, respectively. In ARIC, the values of R2 at 150 and 250 rpm are 0.92 and 0.85, respectively. These indicate that the use of the modified drag model in CFD simulations predicts satisfactory results of dispersed phase RTD. 5.6.1. Tracer Flow Profile at Different Locations in the Column. In order to understand the flow pattern of the dispersed phase at different locations in the column, the tracer flow profile was monitored at different locations along the 17205

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research

3-D simulations of pilot scale geometries is very timeconsuming and requires large computational resources. In addition, CFD simulations have helped in analyzing the flow pattern of the dispersed phase near and above critical agitator speed (NCR) in ARDC and ARIC. Below NCR, the flow of the dispersed phase is governed by the buoyancy of droplets, whereas above NCR, the flow generated by a rotating disk/ impeller regulates the flow pattern of the dispersed phase. The presence of recirculation zones, bypassing of the dispersed phase through clearances at different operating conditions, has been observed through simulations. In addition, CFD simulations have been used to quantify the static and dynamic holdup in the column, which is difficult to measure through experiments. Simulations predicted the effect of rotational speed and impeller geometry on variations of static holdup. The static holdup has been found to be reducing with an increase in rotational speed. Below NCR, the contribution of static holdup to the total holdup is ∼60%. Above NCR, the contribution of static holdup reduces to 20 to 30% of total holdup. Reduction in static holdup at higher agitation speed leads to an increase in extraction efficiency. The methodology developed here can be further used to study the effect of various geometrical parameters like compartment height, eccentricity, and impeller type on the recirculation, bypassing, and static holdup in the column. This will be useful in the optimization of geometrical parameters during design and scale up of such extraction columns. Although the developed CFD model provides satisfactory agreement between predicted and experimental holdup values, significant improvement in the model for prediction of Ped is required. Simulations predicted a continuous increase in Ped with an increase in agitator speed in both the columns, as against the experimental data. This certainly limits the general applicability of the developed CFD model. The discrepancy between simulated and experimental results can be due to the assumption of monodispersed drop size in simulations. Coupling of CFD with population balance models with suitable breakage and coalescence kernels for prediction of Ped and the development of mechanistic formulations for drag can be a way forward in this topic to improve the CFD predictions. Despite such limitations, the developed CFD approach will encourage further research in this area. The RTD study has revealed reduction in the bypass of the dispersed phase through the clearance between the stator and wall at higher rotational speed. Such observations can be utilized to optimize the clearance and rotation speed while designing the LLE columns. The flow patterns of the continuous and dispersed phase can be examined for different agitator designs that will lead to a better understanding of local hydrodynamics. This offers the capability of predicting the hydrodynamic behavior of such extraction columns to investigate the influence of operating and design parameters without actually constructing the experimental set up. The effect of geometrical parameters like vertical baffle position and stage height on the local hydrodynamics can be investigated. This understanding will assist in reducing the uncertainties in design and scale up of such extraction columns.

outlet do not show multiple peaks, indicating the reduction in bypass of the dispersed phase through clearances at higher RPM. The increase in spread of the RTD curve at 800 rpm as compared to 400 rpm signifies an increase in axial mixing at higher RPM. The tracer flow profiles in ARIC (at both 150 and 250 rpm) do not show a split peak at stage 10. Also, the stage 10 profiles at 250 rpm in ARIC and 800 rpm in ARDC are observed to be similar. Thus, the flow of the dispersed phase through clearance in ARIC is negligible or well mixed with the main streamflow. 5.6.2. Peclet Number of Dispersed Phase Flow (Ped). Table 5 compares the values of Ped obtained through experiments Table 5. Dispersed Phase Peclet Number (Ped) of ssytem 2 Ped type

N (RPM)

experiment

simulation

ARIC ARIC ARDC ARDC

150 250 400 600

13.06 13.99 12.11 12.75

8.68 20.84 5.58 8.07

and simulations in both columns for system 2. The values of Ped obtained through experiments ranges from 12 to 14, indicating that they do not show much variation with agitator speed. However, simulations show an increase in Ped with an increase in agitation speed. The increase in Ped from simulations indicates that dispersed phase flow approaches plug flow behavior. This behavior can be attributed to the assumption of monodispersed drop size in simulations. In reality, there exists a broad range of drop size distribution in extraction columns due to the coalescence and breakage of drops. These processes are intensified by the presence of mixing and settling zones. The drops of different size exhibit a wide range of rise velocities. The presence of larger drops results in forward mixing that cannot be predicted by the assumption of monodispersed drop size. To overcome this limitation, a population balance method can be incorporated in CFD, wherein the drop size distribution is predicted based on coalescence and breakage phenomena.

6. CONCLUSIONS The present paper has provided a methodology to simulate two phase flow in pilot scale ARDC and ARIC, which can reduce the efforts of experimental studies. The effect of two different drag models which account for the effect of dispersed phase holdup on the drag coefficient has been evaluated for the prediction of holdup. It has been observed that the original KH drag model is not able to predict the holdup at a high rotational speed. Increasing the value of the constant parameter in the drag model reduced the error between the predicted and experimental results. The unreliable prediction of holdup with the original KH drag model could be attributed to the complex geometry of the column, wall effects, and droplet interactions in the dispersed phase. To overcome this difficulty, the KH drag model was modified to minimize the error in predictions. This technique lumps all the uncertainties in the modified empirical constants. A value of parameter A equal to 19.4 was observed to give the best fit with an average error of 11% for both ARDC and ARIC. However, a major limitation of this method is that the modified drag model is valid only for the respective test setup and simulation method. Also, this process of optimizing parameters of drag models in



AUTHOR INFORMATION

Corresponding Author

*Tel.: 91-22-33612018. Fax: 91-22-33611020. E-mail: aw. [email protected]. 17206

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research ORCID

PBD = Pitched Blade-Disk PBM = Population Balance Method PDDC = Pulsed Disk Doughnut Column PSPC = Pulsed Sieve Plate Column RDC = Rotating Disk Contactor RPM = Rotations Per Minute RSM = Reynolds Stress Model RTD = Residence Time Distribution

Ashwin W. Patwardhan: 0000-0001-9893-1249 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS R.A.F. wishes to acknowledge University Grants Commission for providing a Ph.D. fellowship. N.V.H. acknowledges DAEICT for providing a Ph.D. fellowship.





NOMENCLATURE A, B = KH drag model constants CD = drag coefficient DO = Orifice diameter (m) Dt = diameter of column (m) Dd = diameter of disk in RDC/ARDC (m) Ds = stator opening diameter (m) Dh = hole diameter in PSPC (m) Ddisk = diameter of disk in PDDC (m) Ddoughnut = diameter of doughnut in PDDC (m) dP = dispersed phase drop diameter (m) F̅ i = external body force (N) Flift,i = lift force (N) Fvm,i = virtual mass force (N) g = acceleration due to gravity (m/s2) Gk,m = generation of turbulent kinetic energy due to mean velocity (m2/s2) Hc = compartment height (m) k = turbulent kinetic energy (m2/s2) Kij = interphase exchange coefficient L = characteristic length (m) N = agitation speed (RPM) NCR = critical rotation speed (RPM) p = pressure (Pa) Pe = Peclet number (−) Ped = Peclet number of dispersed phase (−) R2 = coefficient of determination (−) R̅ ij = interphase exchange force (N) Re = Reynolds number (−) Ui = velocity of ith phase (m/s) Z = column height (m)

REFERENCES

(1) Torab-Mostaedi, M.; Asadollahzadeh, M. Mass transfer performance in an asymmetric rotating disc contactor. Chem. Eng. Res. Des. 2015, 94, 90−97. (2) Haman, J.; Mišek, T. Analysis of extractor performance. 1. Longitudinal mixing in a single-phase system. Chem. Eng. J. 1988, 39, 1−9. (3) Haman, J.; Mišek, T. Analysis of extractor performance. II. Longitudinal mixing of the continuous liquid in a two-phase system. Chem. Eng. J. 1988, 39 (1), 11−16. (4) Mišek, T.; Haman, J.; Rod, V. Analysis of Extractor Performance. III. Longitudinal Mixing in the Dispersed Phase. Chem. Eng. J. 1988, 39, 69−77. (5) Kadam, B. D.; Joshi, J. B.; Patil, R. N. Hydrodynamic and mass transfer characteristics of asymmetric rotating disc extractors. Chem. Eng. Res. Des. 2009, 87 (5), 756−769. (6) Torab-Mostaedi, M.; Asadollahzadeh, M. Mass transfer performance in an asymmetric rotating disc contactor. Chem. Eng. Res. Des. 2015, 94, 90−97. (7) Hendre, N. V.; Venkatasubramani, V.; Farakte, R. A.; Patwardhan, A. W. Hydrodynamics and Mass Transfer Characteristics of Asymmetric Rotary Agitated Columns. Ind. Eng. Chem. Res. 2018, 57, 1630. (8) Kumar, A.; Hartland, S. A Unified Correlation for the Prediction of Dispersed-Phase Hold-Up in Liquid-Liquid Extraction Columns. Ind. Eng. Chem. Res. 1995, 34 (11), 3925−3940. (9) Angelov, G.; Journe, E.; Line, A.; Gourdon, C. Simulation of the flow patterns in a disc and doughnut column. Chem. Eng. J. 1990, 45 (2), 87−97. (10) Fei, W. Y.; Wang, Y. D.; Wan, Y. K. Chem. Eng. J. 2000, 78, 131−139. (11) Drumm, C.; Bart, H. Hydrodynamics in a RDC Extractor : Single and Two-Phase PIV Measurements and CFD Simulations. Chem. Eng. Technol. 2006, 29, 1297−1302. (12) Yadav, R. L.; Patwardhan, A. W. CFD modeling of sieve and pulsed-sieve plate extraction columns. Chem. Eng. Res. Des. 2009, 87 (1), 25−35. (13) Ghaniyari-Benis, S.; Hedayat, N.; Ziyari, A.; Kazemzadeh, M.; Shafiee, M. Three-Dimensional Simulation of Hydrodynamics in a Rotating Disc Contactor using Computational Fluid Dynamics. Chem. Eng. Technol. 2009, 32 (1), 93−102. (14) Onink, F.; Drumm, C.; Meindersma, G. W.; Bart, H.; de Haan, A. B. Hydrodynamic behavior analysis of a rotating disc contactor for aromatics extraction with 4-methyl-butyl-pyridinium·BF4 by CFD. Chem. Eng. J. 2010, 160, 511−521. (15) Drumm, C.; Hlawitschka, M. W.; Bart, H.-J. CFD Simulations and Particle Image Velocimetry Measurements in an Industrial Scale Rotating Disc Contactor. AIChE J. 2011, 57 (1), 10−26. (16) Chen, H.; Sun, Z.; Song, X.; Yu, J. Key Parameter Prediction and Validation for a Pilot-Scale Rotating-Disk Contactor by CFD− PBM Simulation. Ind. Eng. Chem. Res. 2014, 53, 20013−20023. (17) Strand, C. P.; Olney, R. B.; Ackerman, G. H. Fundamental Aspects of Rotating Disk. AIChE J. 1962, 8 (2), 252−261. (18) Sen, N.; Singh, K. K.; Patwardhan, A. W.; Mukhopadhyay, S.; Shenoy, K. T. CFD Simulations of Pulsed Sieve Plate Column: Axial Dispersion in Single-phase Flow. Sep. Sci. Technol. 2015, 50 (16), 2485−2495. (19) Sen, N.; Singh, K. K.; Patwardhan, A. W.; Mukhopadhyay, S.; Shenoy, K. T. CFD simulation of two-phase flow in pulsed sieve-plate



GREEK SYMBOLS αi = volume fraction of secondary phase i ρi = density of ith phase (kg/m3) τ̿i = stress−strain tensor μ = viscosity (N s/m2) ε = turbulent kinetic energy dissipation rate (m2/s3) μeff = effective viscosity (N s/m2) μj = viscosity of primary phase (N s/m2)



ABBREVIATIONS ARDC = Asymmetric Rotating Disk Contactor ARIC = Asymmetric Rotating Impeller Contactor BC = Boundary Conditions CFD = Computational Fluid Dynamics CSTR = Continuous Stirred Tank Reactor KH = Kumar and Hartland LDA = Laser Doppler Anemometry LLE = Liquid Liquid Extraction LPM = Liters Per Minute MRF = Moving Reference Frame 17207

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208

Article

Industrial & Engineering Chemistry Research column - Identification of a suitable drag model to predict dispersed phase hold up. Sep. Sci. Technol. 2016, 51 (17), 2790−2803. (20) Kumar, A.; Hartland, S. Gravity settling in liquid/liquid dispersions. Can. J. Chem. Eng. 1985, 63 (3), 368−376. (21) Yi, H.; Wang, Y.; Smith, K. H.; Fei, W. Y.; Stevens, G. W. CFD Simulation of Liquid-Liquid Two-Phase Hydrodynamics and Axial Dispersion Analysis for a Non-Pulsed Disc and Doughnut Solvent Extraction Column. Solvent Extr. Ion Exch. 2016, 34 (6), 535−548. (22) Deshmukh, S. Computational and Experimental Fluid Dynamics: Design of Liquid-Liquid Multiphase Reactors/Contactors; Ph.D. Thesis, Institute of Chemical Technology, University of Mumbai, 2008. (23) Nygren, A. Simulation of Bubbly Flow in a Flat Bubble Column, Master’s Thesis, Lund University, 2014. (24) Rusche, H. Computational Fluid Dynamics of Dispersed TwoPhase Flows at High Phase Fractions; Imperial College London (University of London), 2003. (25) Barnea, E.; Mizrahi, J. A generalized approach to the fluid dynamics of particulate systems. Chem. Eng. J. 1973, 5 (2), 171−189. (26) Famularo, J.; Happel, J. Sedimentation of dilute suspensions in creeping motion. AIChE J. 1965, 11 (6), 981−988. (27) Brucato, A.; Grisafi, F.; Montante, G. Chem. Eng. Sci. 1998, 53 (18), 3295−3314. (28) Khopkar, A. R.; Kasat, G. R.; Pandit, A. B.; Ranade, V. V. Ind. Eng. Chem. Res. 2006, 45, 4416−4428. (29) Khopkar, A. R.; Ranade, V. V. CFD Simulation of Gas − Liquid Stirred Vessel: VC, S33, and L33 Flow Regimes. AIChE J. 2006, 52 (5), 1654. (30) Laddha, G.; Degaleesan, T. E. Transport Phenomena in Liquid Extraction; McGraw-Hill Companies, 1978. (31) Fogler, H. S. Elements of Chemical Reaction Engineering, second ed.; Pearson Education, 2006. (32) Burns, A. D.; Frank, T.; Hamill, I.; Shi, J. Favre Averaged Drag Model for Turbulent Dispersion in Eulerian Multi-Phase Flows. 5th International Conference on Multiphase Flow, ICMF ‘04, Yokohama, Japan, May 30−June 4, 2004; International Conference on Multiphase Flow, 2004; no 392, pp 1−17. (33) Wang, F.; Mao, Z. Ind. Eng. Chem. Res. 2005, 44, 5776−5787. (34) Molavi, H.; Bahmanyar, H. Local and Average Static Holdups in a Countercurrent Flow Rotary Disc Contactor. Chem. Eng. Commun. 2011, 198, 977.

17208

DOI: 10.1021/acs.iecr.8b02720 Ind. Eng. Chem. Res. 2018, 57, 17192−17208