CFD Contribution to a Design Procedure for Discs and Doughnuts Extraction Columns

0263±8762/98/$10.00+0.00 € Institution of Chemical Engineers Trans IChemE, Vol 76, Part A, November 1998

CFD CONTRIBUTION TO A DESIGN PROCEDURE FOR DISCS AND DOUGHNUTS EXTRACTION COLUMNS M. S. AOUN NABLI, P. GUIRAUD and C. GOURDON Laboratoire de GeÂnie Chimique UMR 5503 CNRS-INPT-UPS, Toulouse, France

T

his paper deals with a Computational Fluid Dynamic contribution to a design procedure for discs and doughnuts extraction columns. The approach can be considered as a numerical alternative to experimental determinations of average hydrodynamic characteristics involved in the phenomenological and transport models used in order to design these contactors. Simulations are carried out by the ESTET calculation code coupled to the k e turbulent model, for steady-state regimes, and for pulsed ¯ ow regimes. In the case of the steady-state regime, correlations between turbulent characteristics of the steady ¯ ow and pressure drop in the column with geometrical parameters and dynamic condition are established. The numerical results are in very good agreement with the experimental ones available. In the case of the pulsed ¯ ow regime, for which 24 different combinations of geometrical and pulsation parameters have been examined, the results presented consist of the averages over a pulsation cycle of hydrodynamic and mixing characteristics correlated to the pulsation amplitude and frequency, the distance between disc and doughnut and the open free area of the internals. The numerical results are also in fairly good agreement with the experimental ones available. Keywords: extraction column design; discs and doughnuts columns; pulsed ¯ ow; CFD

INTRODUCTION

particles (breakage, coalescence, interphase mass transfer, diffusion in solid particles) in a global simulation of the apparatus. Phenomenological laws are coupled with 1D models describing the transport of the phases (plug ¯ ow with axial dispersion for both phases), the evolution of the dispersion being taken into account via population balances (Casamatta3 , Zimmerman et al.4 ). The phenomenological laws commonly used are related to hydrodynamic characteristics such as the average of the turbulent energy and of its dissipation rate, derived from experimental measurements of the pressure drop in the column. The usual way to obtain the different parameters required for the transport models is based on tracer measurements of the Residence Time Distribution (RTD). These experimental determinations present the major disadvantage of being very tedious and expensive. Moreover, this type of experiment, performed on pilot plants, does not solve the problem of scaling up. The numerical study of hydrodynamics and axial mixing appears, since 1980 (Oh 5 ), the most promising way to understand the workings of these columns for different internal design and pulsation conditions. In fact, since the complex experimental investigations of these processes are very long and expensive, a large development of computer-aided modelling was stimulated as a useful tool in the study and design of these contactors (Casamatta3 , Gourdon1 ). The aim of this study is to propose a numerical alternative to these experiments by the means of hydrodynamic simulations within the discs and doughnuts columns. The ¯ ow generated in the column is complex: two-phase, unsteady and turbulent, presenting a double periodicity in

Solvent extraction is one of the most widely used unit operations in the process industry, especially in the treatment of waste in nuclear fuel reprocessing. The operation consists of separating one or several substances (solute) present in a solid or a liquid phase by the addition of another liquid phase in which these substances are transferred preferentially. This mass transfer is often operated in counter-current contactors. One of the most frequently encountered contactors is the discs and doughnuts pulsed column. Discs and doughnuts, alternately arranged along the column, combined with a mechanically imposed pulsation, make the ¯ ow unsteady and turbulent. The major design problem, in the case of an industrial column, lies in the choice of the geometrical parameters of the packing and of the pulsation conditions (Gourdon1 ). The requirements are ® rst, a suf® cient interfacial area production (by means of the turbulence needed to drop breakage), and, secondly, adequate contact time between the phases. Nevertheless, the axial mixing generated by the pulsation intensity has to be kept at the lowest level because it reduces the ef® ciency of the column by decreasing the mass transfer driving force. Neglecting the effect of axial mixing when designing the columns can lead to an overestimation of the mass transfer ef® ciency of about 30% or more (Buratti2 ). Besides, the pressure drop, enhanced by the packing and the pulsation intensity, diminishes the throughput. The standard approach to elaborate a design procedure for the industrial columns is based on the integration of the different phenomena concerning the drops or the solid 951

952

AOUN NABLI et al.

space (discs and doughnuts are alternately arranged) and time (sinusoõÈ dal pulsation). The present work intends to perform a synthesis of the results that have been obtained with single phase predictions. The chronology is as follows: ® rst, the work of the Oh5 relative to the laminar unsteady case, then Angelov et al.6,7 , for the turbulent steady ¯ ow, and then Legarrec8 , again under laminar unsteady ¯ ow conditions, but the ® rst to study axial mixing numerically. Finally, the work of Aoun Nabli9 , for turbulent unsteady ¯ ow with axial mixing. Unfortunately, the major part of the experimental results available have been obtained for the steady state regime, because measurements under unsteady conditions are very dif® cult to perform. As a consequence, and also because the condition of steady state is necessary in order to perform unsteady-state simulations, the ® rst part of this work has been devoted to the simulations of the steady-state regime under different geometrical and dynamic conditions. Correlations, allowing the prediction of hydrodynamic characteristics according to the design (packing geometry) and the dynamic condition, summarize the local results obtained via the Computational Fluid Dynamic approach with the industrial code ESTET, developed at the National Laboratory of Hydraulics (EDF Chatou-France) and commercialized by SIMULOG company. The numerical correlations established are compared with some available experimental ones. In the second part of the work, the same approach has been applied to the oscillating regime. As in the case of the steady-state regime, the predetermination of the local monophasic unsteady pulsed ¯ ow structure allows, on one hand, an understanding of the behaviour of the ¯ ow within the apparatus and, on the other hand provides the appropriate correlations between the turbulent characteristics, pressure drop, axial mixing, and the operating conditions and geometrical parameters of the column.

Discs and doughnuts columns are made up of a vertical cylindrical body which contains a packing of discs and doughnuts alternately arranged and maintained by braces of constant length, giving to the ensemble a spatial periodicity. Figure 1 shows a schematic representation of such a column, showing the main geometrical parameters. The two phases (continuous and dispersed) are ¯ owing countercurrently. This counter-current ¯ ow is obtained via the difference of density between a rising light phase and a ¯ owing down heavy phase. Two settlers, with a larger diameter than the active part of the column, are located at each end of the contactor. Under pulsating regimes, the entire ¯ ow system oscillates by means of an external pulsation imposed by a piston acting on the lower settler. Operating Conditions In standard operating conditions, the instantaneous bulk velocity u t (¯ ow rate/section area of column ratio) can be split up into a permanent uo , and a periodic up t component according to (1): uo

up t

up t is the component of the bulk velocity generated by the pulsation. In the case where this pulsation is obtained via a piston, up t can be expressed as a sinusoidal function of the time, with two parameters, the amplitude A and the frequency f (2): up t

pAf cos 2pft

1

2

T

1/f is the pulsation period. In industrial operating conditions, the permanent component is negligible compared to the oscillating one (Buratti2 ). Then, the ¯ ows can be considered as purely oscillating, in such a way: ut

THE DISCS AND DOUGHNUTS COLUMN Column Geometry

ut

Figure 1. Schematic representation of a discs and doughnuts pulsed extraction column.

up t

pAf cos 2pft

3

In such a situation, after each half period of the pulsation, the ¯ ow takes the opposite direction in the column. In the case of a stationary regime, the constant bulk velocity U is taken to be equal to the average over a pulsation cycle of the instantaneous bulk velocity (3) of the oscillating regimes. Therefore: U

T

1 T

up t dt

2Af

4

0

Dimensionless Parameters Dimensionless parameters are chosen in order to put into evidence geometrical and dynamic similitude conditions of the ¯ ow in the column (Legarrec8 ). The typical dimensionless parameters for this kind of contactors are as follows: geometrical parameters h

H/D

Tdi

1

Tdo

: ratio of the cell dimensions; Ddi /D

Ddo /D

2

2

: free area of the disc; : free area of the doughnut.

Trans IChemE, Vol 76, Part A, November 1998

CFD CONTRIBUTION TO A DESIGN PROCEDURE As in the general case (Buratti2 , Oh5 , Legarrec8 ), the doughnut and disc free area are equal in this work: Tdi

Tdo

T

: free area of the packing

dynamic parameters The dynamic parameter is restricted, in the case of the steady-state regime, to the Reynolds number Re U D/n, where n is the kinematic viscosity. For the oscillating regimes, the following dimensionless parameters characterize the pulsation conditions: A f

A/D 2

f D /n

: amplitude parameter; : frequency parameter.

NUMERICAL SIMULATION CONDITIONS Physical Hypothesis and Basic Assumptions The numerical experiments are concerned with steady and unsteady turbulent ¯ ows of an uncompressive Newtonian single phase. The turbulence is represented by a standard k e model with logarithmic wall laws as used in the steady-state case. The choice of a k e model has to be justi® ed. This model is built with two main hypotheses. First, the turbulence is isotropic. This condition is not generally ful® lled in complex ¯ ows, but the fairly good results of a great number of simulations show that only the satisfaction at a small scale of this isotropy condition seems to be necessary. Secondly, the turbulence models developed with a turbulent viscosity concept suppose a local spectral equilibrium. This problem has been investigated for pulsed ¯ ows in rectilinear channels (Ramaprian and Tu 10 , Boudghene et al.11 ). In this case, the authors show that the spectral equilibrium and the turbulence structure of the ¯ ow are affected when the frequency of the oscillations and the characteristic frequency of the turbulence interfere with each other. If the turbulence structure is not modi® ed by the pulsation, the results obtained using the k e model are in good agreement with measurements. Oh5 performed a spectral analysis of velocity measurements obtained by Laser Doppler velocimetry in a discs and doughnuts column, the geometrical parameters being (h 0.156; T 23.5%) and pulsation conditions corresponding to a frequency of oscillations around 1 Hz. At several positions in the column, Oh found a 5/3 slope for the inertial range of the spectrum representing the contribution of each frequency to the turbulence energy (frequency spectrum) showing that this inertial range is not disturbed by the pulsation, certainly because the pulsation frequency (about 1 Hz) is out of the inertial range experimentally observed (5±40 Hz). As a consequence, oscillations of frequencies about 1 Hz do not modify the turbulence structure of the ¯ ow and the k e model of turbulence can be used. The resulting set of equations that has to be solved by the ESTET code to simulate the hydrodynamic ® elds in the column (mean and turbulent ¯ ow) is classical and can be found elsewhere9 . The axial dispersion coef® cient is evaluated from the simulation of the mixing behaviour of a tracer in the column. As in the case of an experimental determination, an Trans IChemE, Vol 76, Part A, November 1998

953

amount of tracer is injected in a compartment where the turbulent pulsed ¯ ow blends it. The concentration ® eld is determined by solving a convection-diffusion equation, the turbulent diffusivity of the tracer being set equal to the turbulent viscosity. Calculation Domain and Simulation Strategy The ¯ ow is axisymmetrical and two-dimensional (Oh5 , Angelov et al.6 ). The spatial periodicity of the packing permits the simulation of the ¯ ow using only a small number of compartments8 . A compartment represents a spatial period of the packing, constituted by two consecutive cells. The ¯ ow calculation domain used in this work contains three compartments (6 cells), each one being limited by two successive discs. The axisymmetry of the ¯ ow allows a 2D calculation to be performed between the axis and the wall of the column (Figure 2), in cylindrical coordinates (r, z). The calculation domain is discretized using a structured and irregular mesh. The mesh, in the case of the lower cell dimension-ratio (h 0.156; T 23.5%), consists of: · 33 lines in the radial direction, 6 of which are in the free area of the disc, and 14 in the free area of the doughnut; · in the axial direction, 14 lines in the space separating two successive baf¯ es, and 2 relative to the baf¯ e thickness. The calculation strategy for every case of simulation is as follows: · the steady-state regime, corresponding to the timeaveraged amplitude of the oscillating ¯ ow rate, is ® rstly simulated for every case; · using these steady-states as initial conditions, the entry ¯ ow-rate begins to be pulsed until it reaches a locally periodic ¯ ow; at the end of this step, the hydrodynamics characteristics can be established; · then, the tracer is introduced into the central compartment of the computational domain. In the course of the following pulsations, it is mixed into the surrounding ¯ uid under the action of the mean ¯ ow convection and of the turbulent diffusion.

Figure 2. Calculation domain.

954

AOUN NABLI et al.

The choice of the time-step D t is performed with respect to the following classical conditions of stability and consistency of the numerical solution: · the convection number of Friedrichs and Levy is less than unity: Nc V D t/D 1 < 1; 1 2 · the diffusion number: Nd YD t/D 1 < 2; In all simulated cases, D t lies between 1 ´ 10 3 and 2 ´ 10 3 seconds. A period of pulsation of 1 second is described by 500 to 1000 time-steps. Two to four periods of pulsation are necessary to reach the convergence of the oscillating regime. Simulated Operating Conditions The geometrical parameters (h , T ) of the different simulated column designs vary in the ranges (0.156±0.406) and 17±40%) respectively. For all the cases, the column diameter D and the thickness of the internals are kept constant (D 0.288 mm and e 0.003 m). For the steady-state simulations, the dynamic parameter Re varies between 15,000 and 52,000 (see Table 1). The ® rst value of the table corresponds to the pilot column studied by Oh5 and Laulan10 . For the pulsed regime, f and A lie between 55,406 and 124,582, and 0.057 and 0.200, respectively. These variation ranges correspond, for a diameter D of 0.288 m, to frequencies between 0.5 and 1.5 Hz, and amplitudes between 1.5 and 6 cm. The conditions used for each simulation case are reported in Table 2. SinusoõÈ dal versus time pro® les (constant in the case of the steady-state regimes) are imposed for the axial velocities at each edge of the calculation domain. Hydrodynamic results are extracted from the central compartment where the ¯ ow is established, as has been proved by increasing the number of stages being simulated (Aoun Nabli9 ). STEADY-STATE RESULTS Local Hydrodynamic Structures Numerical results concerning the local hydrodynamic structures (mean and turbulent ¯ ow) have been presented in a previous work9 . They show that the steady-state ® eld does not depend on the Reynolds number Re, but only on the design parameters. Legarrec 8 reached similar conclusions when the Reynolds number is larger than a critical value Rec (Rec 8640 for h 0.234 and T 30%). The mean ¯ ow presents vortices down stream of doughnuts and discs (one vortex per cell). The recirculation volumes increase with the increase of the cell dimensionratio h , and slightly decrease with the increase of the open free area T . The same behaviour has been numerically observed by Angelov et al.6 . The zones, where the turbulent energy is most intensive, are located in the annular aperture between the discs and the column wall; in the aperature of the doughnuts and just beneath the discs (in the case of an upward ¯ ow). These locations have been detected by other authors5,6 . Even though the discs and doughnuts have the same open free area T , the turbulent energy level is more important downstream of a doughnut than downstream of a disc. Qualitatively, the ¯ ow is accelerated downstream of

a disc, and decelerated downstream of a doughnut, leading to a decrease and an increase of the turbulent level, respectively. These velocity variations are simply due to the evolution of the radial area free for the ¯ ow with the radial position. This area decreases when the ¯ ow goes from the wall to the axis of the column under a doughnut, and increases when the same ¯ ow returns from the axis to the wall under a disc. Space-Averaged Dynamic Flow Characteristics The averaged characteristics of the ¯ ow have been calculated from the local value ® elds over the compartment, for every simulated dynamic and geometrical condition (h , T , Re). Table 1 presents the turbulent characteristics (k, e), as well as the pressure gradient and the pressure drop coef® cient calculated between two consecutive discs. Each dimensionless space-averaged characteristic X has been correlated under the following form: j h b T l Red

X

5

The values of j, b, l and d are given in Table 2. Turbulent characteristics-correlations The correlations show that K kL and K e L are respectively proportional to the square and the cube of the Reynolds number. The proportional factor only depends on the geometrical parameters of the internals. The validity domain of the dimensionless parameters is as follows: 15, 000 # 0.15 # 17 #

h # T #

Re #

52, 000

0.41 40%

The correlations, obtained with a mean relative error of about 2% for K kL and 3.5% for K e L , show that the turbulent energy level increases with the decrease of the cell dimension-ratio and/or of the open free area. Combining the correlations of K kL and K e L leads to expressions concerning the mean turbulent scales: 3/2 l tm K k L /K e L being the mean macro-scale (Taylor integral scale) which represents the dimension of the large vortices structures; gm n3 /K e L 1/4 being a mean Kolmogorov micro-scale. The mean macro-scale is practically independent of the cell dimension-ratio h and increases with the transparency T . The turbulent macro-scale size, in the permanent case, is ® xed by the open free area T and not by the cell dimension-ratio h . The turbulent micro structure-scale increases with the increase of geometrical parameters and decreases with the increase in the ¯ ow rate in the column. At a ® xed geometry, the ratio of the mean macro-scale over the mean micro-scale is proportional to Re3/4 . It can be noted that locally, the ratio l t /g is equal to Re3/4 t , where Ret is a local turbulent Reynolds number de® ned as Ret k1/2 l t /n.

Pressure drop law in steady-state regime The pressure gradient calculated between two successive discs represents the pressure drop due to discs and doughnuts singularity, and friction in the column. The Trans IChemE, Vol 76, Part A, November 1998

CFD CONTRIBUTION TO A DESIGN PROCEDURE

955

Table 1. Mean dynamic characteristics of the steady-state ¯ ow.

h* 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.233 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.406 0.406

T * (%)

Re

102 K kL (m2 s 2 )

23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 17.0 23.5 30.0 40.0 23.5 23.5 23.5 23.5 23.5 23.5 23.5

15000 20000 25000 30000 35000 40000 45000 30000 40000 40000 40000 40000 20000 25000 30000 35000 45000 30000 52000

1.4 2.4 3.8 5.5 7.5 9.8 12.3 4.3 7.1 6.1 5.5 4.2 1.5 2.3 3.4 4.7 7.7 2.8 8.6

K eL (m2 s 3 )

D P 2r H e (m s 2 )

Cpdc (m 1 )

0.052 0.122 0.241 0.411 0.651 0.966 1.357 0.276 0.711 0.496 0.419 0.226 0.062 0.120 0.209 0.335 0.702 0.158 0.822

1.5 2.7 4.3 6.2 8.2 10.7 13.6 4.1 9.2 5.5 4.3 2.3 1.4 2.1 3.1 4.3 7.0 2.3 7.1

553 560 571 571 555 555 557 378 477 285 223 119 290 279 286 291 287 212 218

simulation results show that the pressure gradient is proportional to the square of the Reynolds number: Table 2. Coef® cients of the steady-state correlations.

D2 K k L n2 D4 K k L n3 l tm D gm D Cpdc

j

b

l

d

0.62

0.68

0.57

2

3.2

0.96

1.16

3

0.15

0.06

0.3

0

3/4

0.24

0.29

3/4

2.9

0.95

1.54

0

D P rL

Cpdc

n 2 2 Re D

6

This result has also been shown experimentally by many authors5,12,13,14 . The proportional factor Cpdc , called the pressure drop coef® cient, depends only on the geometrical parameters of the internals. From all the numerical results, the correlation for Cpdc with geometrical parameters has been obtained with a mean relative error of about 2.6%. This correlation in the steady-state regime shows that the pressure drop increases with the decrease of h and/or T , in other words with the con® ning of the ¯ ow by the internals. Figure 3 shows that the numerical results obtained in this study are in very good agreement with the experimental ones of Milot13 and Leroy 14 in the steady-state regime in the same type of column. Bracou15 has studied the pressure drop in a column equipped with truncated discs. The lower value of the pressure drop coef® cient obtained experimentally by this author in comparison with the numerical data shows that the truncated discs dissipate less energy than discs and doughnuts packing. OSCILLATING REGIMES RESULTS Local Hydrodynamic Structures

Figure 3. Variation of the dimensionless pressure drop coef® cient Cpdc D according to h (T 23.5%).

Trans IChemE, Vol 76, Part A, November 1998

Numerical results relative to the time evolution of the local hydrodynamic structures (mean and turbulent ¯ ow) have been detailed in a previous work9 . Two types of hydrodynamic structures have been observed. A regime characterized by unstable large vortices in the case of low values of h , and another, with permanent stable vortices in each stage of the column, when h is larger. By comparison, the other parameters show a slight in¯ uence on the ¯ ow structure. These two different hydrodynamic regimes have been also observed by Oh5 and Legarrec8 .

956

AOUN NABLI et al.

Space-averaged Flow Characteristics over the Pulsation Cycle As has been done for the steady-state case, the mean characteristics of the ¯ ow, averaged over a compartment, have been derived from the instantaneous local values, for every simulated case. This concerns: the turbulent kinetic energy and its dissipation rate which are the classical parameters of the drop breakage modelling, the vorticity ® ®V®(calculated over a cell limited by a doughnut V 12 rot and by a disc) which represents the importance of the ¯ ow shearing, and the instantaneous pressure gradient calculated between two successive discs. Figure 4 shows the time evolution of K k L and K e L with A 0.173; f 82944; h 0.156 and T 23.5%. The initial time (t 0) corresponds, for all the cases, to the maximum of the ascending ¯ ow rate. The time evolution of the space-averaged values (over a compartment) of K kL and K e L are periodic, with a period equal to the half of the pulsation cycle. In fact, the ¯ ow structure in a compartment (two cells) during the ® rst half of the pulsation period is symmetrical to the one during the second half of the pulsation period. The time evolutions of the space-averaged values over only one cell of k and e (not represented here) are periodic with the same period as that imposed by the piston. The maxima of K kL and K e L are observed during the deceleration phases of the ascending or descending ¯ ow, the minima being during the acceleration phases. This is in agreement with experimental results about oscillating ¯ ows in channels (Hino et al.16 , Sergreev17 ) and with laser Doppler velocimetry measurements performed at some points in the column (Zaharieva et al.18 ), and con® rms that turbulence is enhanced during the deceleration phases. The time evolution of the one cell-averaged vorticity (the cell lower limit being a doughnut) is periodic with the same period as the one imposed by the pulsation (Figure 5). In the case of the immediately adjacent cell, a similar time evolution is obtained, with a shift equal to half a pulsation. The vorticity presents two maxima. Each of them occurs at the beginning of the deceleration phases (deceleration of the ascending ¯ ow or of the descending ¯ ow), after the ¯ ow rate has reached its maximum. The maxima can be easily related to the size of the vortices generated within the cell downstream of the discs and of the doughnuts. The larger the vortice, the more important is the vorticity: as a

Figure 5. Time evolution of the mean over one cell of the vorticity and of the pressure gradient over a compartment; (h 0.156; T 23.5%; A 0.174; f 82944).

consequence, the maximum of the vorticity observed at the start of the ascending ¯ ow deceleration phase (t 0.1T, recirculation developed at the edge of the doughnut) is larger than the one observed at the start of the descending ¯ ow deceleration phase (t 0.6, recirculation developed at the edge of the disc of the same cell). Two minima at the beginning of the ¯ ow acceleration phases can also be observed. They also can be related to the size of the vortices present in the cell at these instants. The shear-stress generated by the doughnuts is more important than the one generated by the discs. This explains the fact that the turbulent intensity is more important down stream of the doughnuts than down stream of the discs, as has been shown in the local study9 . The time evolution of the pressure gradient is periodic with the same period than the one imposed by the pulsation (Figure 5). Near the maximum ¯ ow rate instants, disturbances of pressure have been found. In fact, similar disturbances have been experimentally observed by Leroy 14 and Bracou15 in similar columns, but with pneumatic pulsation which is not rigorously sinusoidal. Space and Time Averaged Flow CharacteristicsÐ Correlations with Dynamic and Geometrical Parameters The search for correlations concerning the evolution of the averaged characteristics according to the geometrical and dynamic parameters represents a very important step for the design and scale-up of extraction columns. They are established, using the simulation results obtained, in the variation range of parameters presented in Table 3. Table 3 contains the average over one pulsation cycle of the dynamic characteristics for all the simulated oscillating cases. All obtained correlations have the same following form: Y

ah b T c A d f

e

7

where Y is the dimensionless space and time averaged ¯ ow characteristic; a, b, c, d and e are the coef® cients of the correlations reported in Table 4. These coef® cients are identi® ed by means of a ® tting method. Figure 4. Time evolution of the space-averaged values of k and e; (h 0.156; T 23.5%; A 0.174; f 82944).

Turbulent characteristics The average of the turbulent dissipation rate K e L , Trans IChemE, Vol 76, Part A, November 1998

CFD CONTRIBUTION TO A DESIGN PROCEDURE

957

Table 3. Mean dynamic characteristics and axial dispersion coef® cients of the oscillating regimes.

(m s )

K kL (m2 s 3 )

K VL (s 1 )

D P T 2r H e (m s 2 )

(m2 s 1 )

1.2 2.0 3.0 4.2 5.5 6.9 8.5 1.4 2.5 6.2 8.6 11.4 3.5 5.2 4.1 3.8 2.7 1.6 2.1 2.5 3.4 5.1 2.2 5.7

0.046 0.088 0.166 0.270 0.403 0.572 0.783 0.056 0.128 0.476 0.775 1.163 0.169 0.360 0.219 0.187 0.104 0.051 0.071 0.097 0.162 0.299 0.077 0.348

2.6 3.9 5.0 5.9 6.9 7.7 8.4 3.5 4.6 7.1 8.4 8.9 5.6 6.9 5.5 5.4 4.1 3.3 3.8 4.1 5.0 6.3 4.3 6.3

1.4 2.9 4.0 5.2 6.7 8.3 10.1 2.0 3.3 7.6 10.4 13.7 3.3 5.4 3.3 2.6 1.7 1.4 1.8 2.1 2.8 4.3 1.7 4.0

0.60 0.83 1.00 1.16 1.32 1.48 1.60 0.70 0.94 1.36 1.56 1.72 2.27 4.41 4.63 4.61 4.41 2.69 3.11 3.63 4.11 5.00 ± 7.65

102 K kL cases

h

T (%)

A

O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O15 O16 O17 O18 O19 O20 O21 O22 O23 O24

0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.233 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.313 0.406 0.406

23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 17.0 23.5 30.0 40.0 23.5 23.5 23.5 23.5 23.5 23.5 23.5

0.057 0.077 0.096 0.115 0.134 0.154 0.173 0.115 0.115 0.115 0.115 0.115 0.115 0.154 0.154 0.154 0.154 0.077 0.096 0.115 0.134 0.173 0.115 0.200

2

f 82944 82944 82944 82944 82944 82944 82944 55406 69258 96961 110813 124582 82944 82944 82944 82944 82944 82944 82944 82944 82944 82944 82944 82944

correlated to the column parameters with a mean relative error of about 3.1%, increases with the decrease of the cell dimension ratio h and of the open free area T . This correlation also shows that the dissipation rate increases more quickly with the pulsation frequency than with the pulsation amplitude. The turbulent kinetic energy K kL varies with the column parameters in the same manner as its dissipation rate. The mean relative error being about 2.1%. The correlations for K k L and K e L show clearly that A and f have a different contribution to the turbulence level in the contactor. Nevertheless, in the works concerning the pulsed columns, the majority of authors had studied the pulsation parameters without distinction between the part due to A and that due to f (Oh 5 , Angelov et al.7 ).

a K kL 2af K kL 2af

d

e

2.4 ´ 10

3

0.68

0.69

0.28

0.45

1.4 ´ 10

3

3

1.29

1.34

0.37

0.62

8.4 ´ 10

2

0.27

0.31

0.05

0.05

3.05

0.32

0.33

0.66

0.90

2.1

0.40

0.54

0.91

0.13

K VL f

Cpdc

c

2

l tm D gm D

D P rL

b

T

3.3 ´ 10

2

1.22

1.30

0.39

0.37

2.7 ´ 10

2

1.22

1.30

0.39

0.37

Trans IChemE, Vol 76, Part A, November 1998

103 Dis

These authors proposed then a correlation according to the pulsation intensity (Af ) with the following form: K kL

Ck Af

2

8

where Ck is a coef® cient which only depends on internals design. Oh5 proposes, from an experimental study of the turbulent ® eld in the column (h 0.156; T 23.5%), a coef® cient Ck 24.9 obtained for a pulsation intensity in the interval 5 # Af # 21 ´ 10 3 m s 1 . Using the correlation established in this work (Table 4) with the same geometrical conditions, and with dynamic parameters giving a pulsation intensity Af 16.5 ´ 10 3 m s 1 (in the range studied by Oh5 ), the following values of Ck are obtained: Ck Ck

Table 4. Coef® cients of the oscillating regimes correlations.

2

20.2 with A 33.6 with A

0.115 and f 0.057 and f

41472; 82944.

These results re¯ ect well the experimental one of Oh5 (Ck 24.9) but also put into evidence the respective in¯ uence of A and f that are globally accounted for via the pulsation intensity Af in the previous work. Angelov7 proposes, in a work based on numerical simulations of the oscillating regime with a quasi-steady ¯ ow hypothesis, the following correlation concerning Ck : Ck

1.85

1 h

0.67

1 T

0.65

1

9

obtained for the following variation range of geometrical parameters: 0.1 # h # 0.4 and 20 # T # 40% and for a pulsation intensity varying in the interval 8.7 ´ 10 3 # Af # 26 ´ 10 3 m s 1 . Applied to the experimental conditions of Oh5 , the correlation proposed by Angelov7 gives

958

AOUN NABLI et al.

Ck 9.3. This underestimated value in comparison to Oh’ s value and to the ones presented in the present work can be attributed to the quasi-steady ¯ ow hypothesis he used. Also, it can be noticed that h and T have practically the same exponent in the expression of K kL as the present study (Table 4) and in the relation (9) developed by Angelov7 . This remark can be explained by the fact that the geometrical parameters (h , T ) appear with almost the same exponents in the correlation obtained for steadystate regimes (Table 2) and in the one developed for oscillating regimes (Table 4), giving less importance to the quasi-steady ¯ ow hypothesis of Angelov. As has been done for the steady-state regimes, a macro scale of turbulence l tm K k L 3/2 /K e L (Taylor scale) has been calculated from the kinetic turbulent energy and its dissipation rate. The exponents are reported on Table 4. The Taylor scale really only depends on the geometrical parameters in the investigated variation range of dynamic parameters. This scale increases with the increase of h and/or T . It can be noticed that the integral scale computed for oscillating cases is more sensitive to the cell dimensionratio h that the one obtained in the case of permanent regimes. In the steady state ¯ ow, the size of the vortices generated at the edge of the discs and at the edge of the doughnuts seems to be related to the radial dimension of these packing elements. On the contrary, in the case of oscillating ¯ ows, the largest vortices are generated when the ¯ ow rate is vanishing in the column. At this instant, vortices occupy the entire volume of the cells; as a consequence, their sizes are related to the distance h separating the packing elements. A Kolmogorov micro-scale gm n3 /K e L 3/4 is also given by a correlation in Table 4. For ® xed geometrical conditions, this micro-scale decreases with the increase of dynamic parameters, whereas the macro-scale remains constant. By increasing the pulsation intensity, the inertial part of the turbulent spectrum reaches smaller structures. In the case of the micro-scale, the effect of geometrical and dynamic parameters are practically similar in permanent and oscillating regimes. Mean ¯ ow-vorticity From the simulation results, the dimensionless correlation concerning the averaged vorticity K V L over one cell and over a pulsation cycle has been developed with an accuracy better than 4%. This relation (Table 4) shows that the shearing induced by the mean ¯ ow increases with the increase of the dynamic parameters and the decrease of geometrical parameters of the column, as is the case for the turbulent characteristics of the ¯ ow. Pressure drop The average over a pulsation cycle of the absolute value T of the pressure gradient D P /rL T 1/T 0 D P /rL dt calculated between two successive discs is correlated to the column parameters (Table 4) with an accuracy better than 3%. This relationship shows that the pressure drop in the column increases with the increase of the dynamic parameters and the decrease of the geometrical parameters, as was the case in the permanent regime. The instantaneous pressure gradient can be written as follows by application of the fundamental principle of

dynamics between two consecutive discs: D P rL

c t

g

D Hf , s L

10

with c : global acceleration of the ¯ uid mass c dU t /dt 2p 2 Af 2 sin 2pf t , D Hf , s : the instantaneous pressure drop due to the friction on the column wall and to the singularities of the internals. Most of the authors (Oh 5 , Leroy 4 , Laulan12 , Milot13 ) estimate the pressure drop in the oscillating case in the same way as for the permanent case, by means of a pressure drop coef® cient Cpdc . The results obtained for Cpdc in the permanent regime are assumed to be valid in the oscillating cases at each instant of the period (quasi-steady ¯ ow hypothesis): Cpdc L ut ut g

D Hf , s

11

These authors suppose then that Cpdc is independent of time and of the pulsation parameters, and only depends on the geometrical parameters as is the case in the steadystate ¯ ow regime. Under these assumptions and by integrating the absolute value of the pressure gradient (10) over a pulsation cycle, the following relation is obtained: D P rL

Cpdc T

p2 2 Af 8

2

12

Using the correlation for the pressure gradient in Table 4 and the previous relation (12), a correlation can be derived for the pressure drop coef® cient Cpdc in the oscillating case (the coef® cients are reported in Table 4). It is also clear that Cpdc cannot be considered as a constant and depends on the geometrical and, also, on dynamic parameters. This result invalidates the hypothesis used in previous works about discs and doughnuts columns5,12 14 . The authors supposed that the Cpdc value obtained in the permanent regime is valid in the case of the oscillating regime, and, as a consequence, is not a function of A and f . Moreover, the value obtained for Cpdc for oscillating regimes is larger than the one obtained for the permanent regime. In the case of the lower cell dimension-ratio (h 0.156; T 23.5%), Cpdc of the oscillating regime varies in the interval [777; 1200 m 1 ] for the studied range of dynamic parameters. Meanwhile, for the permanent case, the application of the numerical correlation (Table 2) or of the experimental correlation of Leroy 14 gives a coef® cient Cpdc ~ 500 m 1 . In fact, the oscillating turbulent regime dissipates more energy than the permanent turbulent regime. It is important to note that the only experimental result performed for the oscillating regime (Laulan12 ) gives a value of Cpdc 820 m 1 , in the range obtained by the simulation for the same geometry. Axial mixing of the continuous phase The axial mixing in the column has been simulated by solving a conservation equation for a tracer injected in the central compartment and dispersed under the action of turbulent pulsed ¯ ow. From the time evolution of the tracer cloud (concentration ® elds), an instantaneous axial dispersion coef® cient has been computed by taking the time derivative of the variance of the axial mass distribution of Trans IChemE, Vol 76, Part A, November 1998

CFD CONTRIBUTION TO A DESIGN PROCEDURE

Figure 6. Variation of Dis with pulsation frequency (A h 0.156).

0.115;

the tracer. This instantaneous axial dispersion coef® cient oscillates with a period equal to the half pulsation cycle, around a constant value. This value is the axial dispersion coef® cient Dis involved in the plug ¯ ow with the axial dispersion model used to describe the transport phenomenon in the population balance based simulations (Zimmermann et al.4 ). Details about the procedure used to estimate Dis are reported elsewhere (Aoun Nabli19 ). This numerical procedure has been applied to 24 combinations of the geometrical and pulsation parameters. The obtained values of Dis are reported in Table 1, and are correlated to the geometrical and dynamic parameters, the coef® cients of the correlation being reported in Table 4. This correlation, obtained with a mean relative error less than 4.5%, shows that the axial dispersion coef® cient is independent of the open free area, increases with the cell dimension-ratio and the pulsation parameters, and is more dependent on the frequency than on the amplitude of the pulsation. Buratti2 experimentally observed the same tendencies using an electrolytic drawing method, especially the independence of Dis with regard to the open free area. Figure 6 shows a comparison between the numerical values of Dis obtained in this work, with Oh5 ’ s and Buratti2 ’ s experimental ones as a function of the frequency parameter. Even though numerical experimentation seems sometimes to underestimate the axial mixing intensity, the results are in acceptable agreement with the experimental values. The accuracy is suf® cient for the design and scale-up of extraction columns. CONCLUSIONS In this study, a numerical contribution to a design procedure for discs and doughnuts extraction columns was presented. This CFD based approach can be considered as a numerical alternative to long and tedious experimental measurements. Knowledge about the in¯ uence of the geometrical and dynamic parameters on the hydrodynamics of the discs and doughnuts pulsed columns has been gained by this means. At the present time, it would not be possible to extract so much information from experiments. The numerical results obtained in this work have been successfully compared with the available experimental ones, as far as design is concerned. Improvement of the simulation to ® t better with the experimental results is dif® cult to propose without a more precise experimental knowledge of Trans IChemE, Vol 76, Part A, November 1998

959

the hydrodynamics within these columns. Particle Image Velocimetry seems to be an interesting technique to use in order to perform a time and space resolute investigation of the local ¯ ow velocities which will allow a precise comparison with simulation results. The averaged values of the turbulence characteristics and of the axial dispersion coef® cient have been extracted from local ® elds and are then available to be implemented into the phenomenological laws describing drop break-up or coalescence and in 1D transport models, via the correlations presented in Table 4. This integration allows a pre-design of industrial columns. In a second stage of the design, for which precise simulations are needed, the numerical procedure presented here can be used again with the actual dynamic and geometrical parameters of the apparatus in order to re® ne the values of turbulence and of axial mixing generated by the equipment. Obviously, this procedure can also give very useful knowledge about existing columns. The next phase of our research concerning the ¯ ow behaviour in extraction columns consists of the study of the comportment of the dispersed phase by Lagrangian simulations of droplets trajectories. Always with a design objective, this will allow the prediction of the mean residence time and the axial dispersion coef® cient of the dispersed phase, taking into account properties of the dispersed phase such as the size distribution and the density difference with the continuous phase. NOMENCLATURE A A a, b, c, d, e Ck Cpdc D Ddi Ddo Dis e f f g H h H/D k L 2H e Nc V D t/D l Nd YD t/D l 2 Re U D/n Ret k1/2 lt /n r T t Tdi 1 Ddi /D Tdo Ddo /D 2 T Tdi Tdo U 2 Af

uo up t V z

pulsation amplitude, m dimensionless pulsation amplitude coef® cients of the oscillating-regime correlations coef® cient, equation (9) pressure drop coef® cient, m 1 diameter of the column, m diameter of a disc, m diameter of a doughnut, m axial dispersion coef® cient, m2 s 1 thickness of the baf¯ es, m pulsation frequency, Hz dimensionless pulsation frequency gravity, m s 2 disc-doughnut spacing, m cell dimension-ratio turbulent kinetic energy, m2 s 2 disc-disc spacing, m convection number of Friedrichs and Levy diffusion number Reynolds number turbulent Reynolds number radial coordinate, m pulsation period, s time, s free area of discs free area of doughnuts packing free area average over a pulsation cycle of the instantaneous bulk velocity, m s 1 bulk velocity: ¯ ow rate/section area of column ratio, ms 1 permanent component of the bulk velocity, m s 1 periodic component of the bulk velocity, m s 1 velocity scale, m s 1 axial coordinate, m

Greek letters j, b, l, d Y

coef® cients of the permanent-regime correlations diffusion coef® cient, m2 s 1

ut

2

960

AOUN NABLI et al.

D Hf , s D P D t D l e g gm lt l tm n V r c

instantaneous pressure drop, m pressure gradient, Pa time step, s space step, m dissipation rate of the turbulent kinetic energy, m2 s Kolmogorov micro-scale, m averaged Kolmogorov micro-scale, m turbulent macro-scale (Taylor integral scale), m averaged turbulent macro-scale, m kinematic viscosity, m2 s 1 vorticity, s 1 density, kg m 3 global acceleration of the ¯ uid mass, m s 2

Subscripts di do o p t m

disc doughnut permanent component periodic permanent component turbulent characteristics mean turbulent characteristics

Operator K L

3

averaged value over a compartment in the permanent regime, and over a compartment and a pulsation period in the oscillating regime.

REFERENCES 1. 2. 3. 4.

Gourdon, C., 1989, Ph D Thesis (INP Toulouse, France). Buratti, M. F., 1988, Ph D Thesis (INP Lorraine, France). Casamatta, G., 1981, Ph D Thesis (INP Toulouse, France). Zimmerman, A., Gourdon, C., Joulia, X. and Gorak, A., 1992, Chem Eng J, 57: 229±236. 5. Oh, W. Z., 1983, Ph D Thesis (INP Toulouse, France). 6. Angelov, G., Journe , E., Line , A. and Gourdon, C., 1990, Chem Eng J, 45: 87±97.

7. Angelov, G., Gourdon, C. and Line , A., 1993, Solvent extraction in the process industries, ISEC’.93, York, England, 2: 1183±1190. 8. Legarrec, S., 1993, Ph D Thesis (CNAM, Paris). 9. Aoun Nabli, M. S., 1995, Ph D Thesis (INP Toulouse, France). 10. Ramaprian, B. R. and Tus, S. W., 1983, J Fluid Mech, 137: 59±81. 11. BoudgheÁne-Stambouli, M. H., Ha-minh, H., 1993, Colloque de meÂcanique des ¯ uides numeÂrique de Toulouse, O13.1±O13.4. 12. Laulan, A., 1980, Ph D Thesis (INP Toulouse, France). 13. Milot, J. F., Duhamet, J., Gourdon, C. and Casamatta, G., 1990, Chem Eng J, 45: 111±122. 14. Leroy, P., 1991, Ph D Thesis (INP Lorraine). 15. Bracou, H., Borda, G., Masbernat, O. and Gourdon, C., 1996, In¯ uence of plate wettability upon dispersion transport in a liquid-liquid extraction column, ISEC’ 96, Melbourne, Australia, 2: 1149±1154. 16. Hino, M., Sawamoto, M. and Takasu, S., 1976, J Fluid Mech, Part 2, 75: 193±207. 17. Sergeev, S. I., 1966, Fluid Dynamics, 1: 121±122. 18. Zaharieva, A., Masbernat, O., Aoun Nabli, M., Guiraud, P. and Gourdon, C., 1996, Experimental and numerical 2-D turbulent ® eld in a pulsed extraction column, ISEC’ 96, Melbourne, Australia, 2: 1233±1238. 19. Aoun Nabli, M. S., Guiraud, P. and Gourdon, C., 1997, Numerical experimentations: a tool to calculate the axial dispersion coef® cient in discs and doughnuts pulsed solvent extraction columns, Chem Eng Sci, 52: 2353±2368.

ADDRESS Correspondence concerning this paper should be addressed to Dr P. Guiraud, Laboratoire de Ge nie Chimique, UMR 5503, CNRSINPT-UPS, 18 Chemin de la Loge, F 31078, Toulouse Cedex 4, France. E-mail: [email protected] The manuscript was communicated via our International Editor for Continental Europe, Professor M. Roques. It was received in July 1996 and accepted for publication after revision 28 July 1998.

Trans IChemE, Vol 76, Part A, November 1998