Monte Carlo simulation of drop dispersion behaviour

Chemical Engineering and Processing 39 (2000) 201 – 206 www.elsevier.com/locate/cep

Monte Carlo simulation of drop dispersion behaviour M. Balmelli, L. Steiner * Swiss Federal Institute of Technology (ETH), Department of Industrial and Engineering Chemistry, 8092 Zurich, Switzerland Received 5 March 1999; received in revised form 10 May 1999; accepted 10 May 1999

Abstract To check the equations for breakage and coalescence rates in turbulent dispersions a simulation program based on the Monte Carlo method has been written. The equations of Coulaloglou and Tavlarides (C.A. Coulaloglou, L.L. Tavlarides, Drop size distribution and coalescence frequencies of liquid-liquid dispersions in flow vessels, A.I.Ch.E. J. 22 (1976) 289) were compared to experimental data and the free parameters were recalculated by optimisation. Special values were necessary for each liquid system and setting. © 2000 Elsevier Science S.A. All rights reserved. Keywords: Turbulent dispersions; Breakage; Coalescence; Drop population balance; Monte Carlo simulation

1. Introduction The turbulent dispersions are produced by intensive agitation of two immiscible liquids. One of the liquids is dispersed into drops, which are prevented from settling by continuous input of the stirring energy. The drops may break down into smaller ones or, especially at higher hold-ups of the dispersed phase, coalesce with other drops. The drops size distribution is determined by the balance of the breakage and the coalescence processes, which are in equilibrium in steady state. For a quantitative description drop population balances are used as introduced by Coulaloglou and Tavlarides [1]. However, the corresponding integro-differential equations are complicated and require a numerical solution, which is tedious even with modern desktop computers. The time consumption is prohibitive especially if multistage processes are simulated or values of free parameters evaluated by optimisation. A possible alternative to the numerical solution of the population balance equations is the Monte Carlo simulation, based, as the name reveals, on random numbers. It assumes that a simulation, repeated many times with random values of the parameters, produces a good approximation of the exact solution. The main advantage is the simplicity of the calculation algorithm. * Corresponding author. Present address: Tannrietlistrasse 8, 8166 Niederweningen, Switzerland.

The deviation from the exact solution lies usually between 5 and 10% (Sobol [2]), and may be further decreased with additional computing effort. Spielman and Levenspiel [3] adopted this technique for simulation of chemical reactors, Ramkrishna [4] indicated the connection between the Monte Carlo method and the drop population balance. A sample of the dispersion is considered, observing all events of influence. In a batch stirred cell filled with two immiscible liquids, both breakage and coalescence of drops have to be traced. The probability of the considered breakage and coalescence processes has to be known. There are two different approaches for a Monte Carlo simulation of drop dispersions (Hsia and Tavlarides [5,6]). The first one, ‘time-oriented’, splits the total simulated time into equal time intervals Dt, an event may take place in each interval. If there is no event, the next interval is examined. The second approach, ‘event-oriented’, estimates the time interval between events, shifting the watch ahead of necessary amount of time until the next event occurs. The second approach saves a significant amount of simulation time but requires a guess of the time elapsing between the different events, which changes during the simulation. Rod and Mı´sˇek [7] adopted the time-oriented approach. Bapat et al. [8] investigated the second possibility, however, without discretising the drops into classes. Their procedure allowed only a tracing of 300–500 drops so a simulation of a significant increase of the

0255-2701/00/$ - see front matter © 2000 Elsevier Science S.A. All rights reserved. PII: S 0 2 5 5 - 2 7 0 1 ( 9 9 ) 0 0 0 7 4 - 4

202

M. Balmelli, L. Steiner / Chemical Engineering and Processing 39 (2000) 201–206

Table 1 Breakage and coalescence frequencies Breakage frequency for class j

Overall breakage frequency

Coalescence frequency between classes i and j

Overall coalescence frequency for class i

bj =

Nj g(dj ) V

(2)

fb =

1 M % Nj g(dj ) Vj = 1

(3)

1 Ni Nj 2V 2 h(di,dj )l(di,dj )

(4)

M 1 Ni % Nj 2 2V j=1 h(di,dj )l(di,dj )

(5)

cij =

ci =

M 1 M % Ni % Nj 2V 2i = 1 j = 1 h(di,dj )l(di,dj )

Overall coalescence frequency

fc =

(6)

breakage and coalescence functions g(dj ), h(di, dj ) and l(di, dj ) of Coulaloglou and Tavlarides [1] being related to one drop, it is necessary to form values valid for the respective drop classes and balance zones (Table 1). The simulated dispersion volume is given by V=

p M % N d3 6Xj = 1 j j

The factor 12 in the coalescence terms avoids counting the events twice. A special routine has been written to eliminate self-coalescences. The probabilities of breakage and coalescence are Pb =

fb fc and Pc = 1− Pb = fb + fc fb + fc

2. Description of the simulation procedure

n

% −ln(Ri )

ln(R) i=0 Dt = − while lim n n“ l

=1

(1)

where R is a random number and l is the average number of events for time unit. The simulation is based on a linearly dicretised drop diameter grid. Only binary events are allowed (breakage into two drops and coalescence between two drops). Breakage may produce daughter drops of different size, coalescence is possible between drops with different diameters. The classical

bj and Pb, ji = b(dj,di ) fb

(9)

In case of coalescence, probabilities that class j is affected and interacts with class i are Pc, j =

A guess for the time interval between events randomly distributed in time, i.e. when every event can take place anytime, independent from the previous event, is given by Gordon [9] as

(8)

In case of breakage, probabilities that a drop of class j is affected and that the first daughter drop belongs to class i are Pb, j =

stirrer speed was not possible, since the drop sample would grow too large. In this work drop discretisation is combined with event oriented time handling.

(7)

ci c and Pc,ji = ij fc fc

(10)

The generated random number R is uniformly distributed between 0 and 1. It is adapted to a function with density f(x) by

&

x

f(x) dx

,&

−



f(x) dx = R

(11)

−

For discrete case, di, j \ 0 and the total number of classes M we have d

% f(di )Ddi i=0

,

M

% f(di )Ddi

i=0

Ddi = Dd

=

d

,

M

% f(di ) % f(di )=R i=0

i=0

(12) An example in Fig. 1 demonstrates the choice of the event (breakage or coalescence, (Fig. 1a) and the following choice of the interacting class (Fig. 1b). In case of coalescence, the procedure is similar, using Pc,j and Pc,ji respectively.

Fig. 1. Selecting a process with random numbers. (a) Choosing a process (in this case breakage), (b) selecting the interacting class (in this case class 3).

M. Balmelli, L. Steiner / Chemical Engineering and Processing 39 (2000) 201–206 Table 2 Simulations with the Coulaloglou and Tavlarides [1,11] model

3. Experiment

Average deviations [%] between simulated an experimental characteristic diameters in steady state

Toluene, original constants Toluene, modified constants o-Xylene, modified constants

d32

d30

d10

49.0 11.5 15.7

43.0 9.3 25.7

35.2 18.4 46.2

The first daughter drop in a breakage is selected from I−1

I

% b(dj,di )

i=1 M

% b(dj,di )

i=1

% b(dj,di ) BR3 5

i=1 M

203

(13)

% b(dj,di )

i=1

The corresponding computer program, schematicallyshown in the appendix, was written, where different model functions for breakage, coalescence and daughter drop distribution can be inserted to examine the existing models. The program stops when the simulated distribution reached a steady state characterised by the breakage and coalescence frequencies remaining within 2% for a longer time interval. To prevent instabilities and overflows, the number of simulated drops was constricted to be between 50 and 250 times the number of drop classes, outside of this range the simulated volume was adjusted. Even with an older computer of 75 MHz a simulation of a transient process, such as the approach to a new steady state after a change of the stirrer speed, was computed in few minutes. The simulation program was checked with the model of Rod and Mı´sˇek [7], for which the analytical solution is available. A sample comparison of analytical and simulated drop size distributions is shown in Fig. 2a, a simulation of a transient state is in Fig. 2b. Generally, the difference between simulated and analytical characteristic diameters was under 0.5%.

A stirred cell with a diameter of 150 mm and a height of 70 mm, shaped like a stage of an agitated extraction column of the Ku¨hni type, was constructed. It was operated completely filled with the liquids, the six-blade shrouded-turbine agitator had a diameter of 70 mm and a height of 10 mm. The cell was equipped with measuring devices for determination of drop size distribution and its variation with hold-up and stirrer speed. The dispersion was photographed in a 1–3 mm wide space between two glass cylinders inserted into the cell, the drop sizes were evaluated manually on a digitising tablet. More than 750 drops were counted for each steady state, assuring a good accuracy. Additionally, the specific interfacial area was monitored by a laser-light-attenuation method to indicate that the system was in steady state conditions during the photographing. The cell is schematically shown in Fig. 3, full details have been given elsewhere (Balmelli [10]). The experiments were carried out with two liquid systems, i.e. either toluene or o-xylene dispersed in water. The hold-up of the dispersed phase varied between 5 and 20%, the stirrer speed from 200 to 450 r.p.m. Repeated experiments under the same conditions have proved that the reproducibility was between 2 and 4%. To eliminate the effect of the system ageing (changes in mean diameter up to 15% within a week) all experiments were performed with fresh liquids of analytical purity.

Fig. 3. Layout of the stirred cell with drop size distribution measurement devices. 1, cell with liquid – liquid dispersion; 2, stainless steel plates; 3, glass cylinders; 4,5, laser device; 6, stirrer; 7, computer for evaluation of laser output; 8, mirror; 9, glass prism; 10, camera; 11, flash.

4. Results

Fig. 2. Comparison of Monte Carlo and analytical solutions. Examples of a simulated drop size distribution (a) and of an approach to a new steady state after a change of the speed of agitation (b).

The model of Coulaloglou and Tavlarides [1,11] that is still frequently used as a reference was checked by the Monte Carlo technique and the results compared to experimental data. The following equations describe the breakage and coalescence rates in liquid–liquid dispersions:

M. Balmelli, L. Steiner / Chemical Engineering and Processing 39 (2000) 201–206

204

Breakage frequency: g(6) = CI6 − 2/9D 2/3 s



r CIIs(1 + X) exp − 2 1+X rD6 5/9D 4/3 s r

n

modified constants, Eq. (18) (curve C).

2

(14)

Coalescence frequency: h(6,6%)l(6,6%) = CIII(6 2/3 +6%2/3)(6 2/9 +6%2/9)1/2 × D 2/3 s − CIV

r exp (1+X)





hCrCD 2s r3 6 1/36%1/3 2 3 1/3 s (1+X) 6 +6%1/3

n 4

(15)

Size distribution of daughter drops after a breakage: b(6%,6)=



n

2.4 (26 − 6%)2 exp −4.5 6% (6%)2

(16)

The model contains four free parameters (CI to CIV) which should be determined from experimental data. The original authors have published the following values: CI = 0.40 CIII = 2.8(10 − 6)

CII = 0.08 CIV = 1.83(109)

(17)

Influences of physical properties of the phases, speed of agitation and the hold-up being defined explicitly, the free parameters should allow for particular properties of the experimental set-up, such as deviations from the standard stirred tank dimensions. They certainly should remain constant for a given liquid system and experimental arrangement. However, in our evaluations, a reasonable simulation of the measured steady state drop size distributions was possible only with the four parameters evaluated for each particular experiment individually. Their values deviated from the published ones and varied also with the hold-up and the speed of agitation. By a (perforce) averaging for the toluene system we found:

Fig. 4. Comparison of simulated drop size distributions with experimental data. System toluene – water, 450 r.p.m. A, original constants from Eq. (17), B, constants optimised for this experiment only, C, averaged constants for all experiments with the given system (Eq. (18)).

Typically, though the mean diameters may be within acceptable accuracy, the drop size distribution is properly described with the special constants only. Still worse was the situation when the constants optimised for the toluene system were applied for the experiments performed with o-Xylene. The drop size distributions were shifted to smaller diameters and especially the mean Sauter diameter, which is strongly affected by the largest drops, was not correctly reproduced. A typicalcomparison is shown in Fig. 5. The simulation of the approach to a new steady state after a change of the speed of agitation proved to be the most difficult one. Even with the special constants optimised to provide the best possible fit of the final steady state drop size distribution the simulation of the transients was not satisfactory. An example in Fig. 6 indicates that the

CI = 0.036 (−91%) CII = 0.12 ( +50%) CIII = 4.2(10 – 7) (− 85%) CIV = 1.16(108) ( − 94%) (18) Using the original and the modified constants in simulations the agreement with experimental data shown in Table 2 was obtained. The characteristic diameters dxy are generally defined as [S(nid xi )/S(nid yi )]1/(x – y), d32 is the Sauter (volume/surface) diameter, d30 is the diameter of an equivalent sphere, d10 is the arithmetic mean of the measured diameters. Fig. 4 shows comparison of the drop size distributions calculated using the original constants, Eq. (17) (curve A), the constants optimised for the special experiment in question (curve B) and the

Fig. 5. Comparison of a simulated drop size distribution (line) with experimental data (bars). O-Xylene in water, constants optimised for toluene in water.

M. Balmelli, L. Steiner / Chemical Engineering and Processing 39 (2000) 201–206

Fig. 6. Simulation of the approach to a new steady state after an abrupt decrease in the speed of the agitation. Toluene in water, model parameters optimised to best fit the final drop size distribution.

diameters. Fig. 4 shows comparison of the drop size distributions calculated using the original constants, Eq. (17) (curve A), the constants optimised for the special experiment in question (curve B) and the modified constants, Eq. (18) (curve C). Typically, though the mean diameters may be within acceptable accuracy, the drop size distribution is properly described with the special constants only. Still worse was the situation when the constants optimised for the toluene system were applied for the experiments performed with o-Xylene. The drop size distributions were shifted to smaller diameters and especially the mean Sauter diameter, which is strongly affected by the largest drops, was not correctly reproduced. A typical comparison is shown in Fig. 5. The simulation of the approach to a new steady state after a change of the speed of agitation proved to be the most difficult one. Even with the special constants optimised to provide the best possible fit of the final steady state drop size distribution the simulation of the transients was not satisfactory. An example in Fig. 6 indicates that the model predicted much faster transition than measured experimentally. Such behaviour was typical for all measured data.

5. Discussion The developed Monte Carlo simulation of drop distributions combines the event oriented time handling with discretisation of drops into pre-set classes of uniform broadness. With this technique it is possible to investigate the breakage and coalescence phenomena without solving the population balance equations and thus saving a considerable amount of the computer time. The comparison to a model with available analytical solution has shown that a result accurate to 2–5% may be obtained in several minutes, even with a rela-

205

tively slow desktop computer. The method is therefore suitable not only for the simulation but also for simpler optimisations and evaluation of free parameters from experimentally measured drops size distributions. The check of the well-known equations of Coulaloglou and Tavlarides indicated a very limited applicability. The four free parameters depend not only on the arrangement of the experimental equipment but also on hold-up, speed of agitation, physical properties of the liquid system and its purity. Replacing toluene with o-xylene has changed the results so that the constants evaluated for the former system could not be used for the latter one. Even within a single system and in the same apparatus no constant values could be averaged. For a reasonable fit the ‘constants’ would have to be replaced by empirical functions of (at least) hold-up and speed of agitation which, however, would destroy the theoretical background of the original equations for the breakage and the coalescence rates. We believe that a general prediction of the breakage and coalescence rates from semi-theoretical equations still remains very unreliable and that special empirical equations should be evaluated for each liquid system and experimental arrangement. Rather than modifying the constants in the above mentioned equations it is more advantageous to use simple power functions to characterise the variation of the breakage and the coalescence rates with drop diameter, intensity of agitation and hold-up as described by Laso et al. [12]. Such equations enable a realistic simulation for a given system and experimental conditions and may be used for extrapolation and a limited scaling-up. Notation bj ci,cij C D d10 d30 d32 Ds fb, fc g(6), g(d) h(6,6%), h(d, d%) M n N P

breakage frequency for class j (m−3 s−1) coalescence frequency for class j and between classes i and j (1/m−6 s−1) Empirical constant drop diameter (m) mean drop diameter (m) diameter of equivalent sphere (m) Sauter mean drop diameter (m) Stirrer diameter (m) overall breakage frequency (m−3 s−1) overall coalescence frequency (m−6 s−1) breakage frequency functions (m−6s−1) drop collision frequency function (s−1) number of drop classes numerical drop concentration (m−3) number of drops Probability

206

M. Balmelli, L. Steiner / Chemical Engineering and Processing 39 (2000) 201–206

Appendix A. h l l(6,6%), l(d,d%) n(6) (n(6)) r s r R t Dt 6 V X Subscripts i,j b,c C,D

Greek symbols b(6,6%) (b(d,d%)) h l l(6,6%), l(d,d%) n(6) (n(6)) r s

Viscosity (kg/ms) Average number of events (s−1) Collision efficiency between two drops Number of daughter drops in breakage density (kg/m3) Interfacial tension (N/m) stirrer speed (s−1) Random number time (s) time interval (s) drop volume (m3) volume of the balance zone (m3) hold-up of dispersed phase drop classes Breakage and coalescence respectively Continuous and dispersed phases respectively

drop size distribution of first daughter drops in breakage Viscosity (kg/ms) Average number of events (s−1) Collision efficiency between two drops Number of daughter drops in breakage density (kg/m3) Interfacial tension (N/m)

References [1] C.A. Coulaloglou, L.L. Tavlarides, Drop size distribution and coalescence frequencies of liquid-liquid dispersions in flow vessels, A.I.Ch.E. J. 22 (1976) 289. [2] I.M. Sobol, Die Monte Carlo Methode, Deutscher Verlag der Wissenschaften, Berlin, 1991. [3] L.A. Spielman, O. Levenspiel, A Monte Carlo treatment for reacting and coalescing dispersed phase systems, Chem. Eng. Sci. 20 (1965) 247. [4] D. Ramkrishna, Analysis of population balance IV, the precise connection between Monte Carlo simulation and population balances, Chem. Eng. Sci. 36 (1981) 1203. [5] A.M. Hsia, L.L. Tavlarides, A simulation model for homogeneous dispersions in stirred tanks, Chem. Eng. J. 20 (1980) 225.

[6] M.A. Hsia, L.L. Tavlarides, Simulation analysis of drop breakage, coalescence and micromixing in liquid – liquid stirred tanks, Chem. Eng. J. 26 (1983) 189. [7] V. Rod, T. Mı´sˇek, Stochastic modelling of dispersion formation in agitated liquid – liquid systems, Trans. Inst. Chem. Eng. 60 (1982) 48. [8] P.M. Bapat, L.L. Tavlarides, G.W. Smith, Monte Carlo simulation of mass transfer in liquid – liquid dispersions, Chem. Eng. Sci. 38 (1983) 2003. [9] G. Gordon, System Simulation, Prentice-Hall, Englewood Cliffs (NJ), 1969. [10] M. Balmelli, Simulation turbulenter dispersionen in ein—und mehrstufigen Extraktionsanlagen. Ph.D. Thesis No. 12597, ETH Zu¨rich (1998). [11] C.A. Coulaloglou, L.L. Tavlarides, Description of interaction processes in agitated liquid – liquid dispersions, Chem. Eng. Sci. 32 (1977) 1289. [12] M. Laso, L. Steiner, S. Hartland, Dynamic simulation of liquid– liquid agitated dispersions, Chem. Eng. Sci. 42 (1987) 2429.