Simulation of the Dynamic Behavior of Continuous Crystallizers - ACS

Chapter 12 Simulation of the Dynamic Behavior of Continuous Crystallizers 1

2

1

Herman J. M. Kramer , Sjoerd de Wolf , and Johan Jager 1

2

Downloaded by UNIV OF AUCKLAND on May 8, 2015 | http://pubs.acs.org Publication Date: September 21, 1990 | doi: 10.1021/bk-1990-0438.ch012

Laboratory for Process Equipment and Laboratory for Measurement and Control, Delft University of Technology, Mekelweg 2, Delft, Netherlands

The simulation of a continuous, evaporative, crystallizer is described. Four methods to solve the nonlinear p a r t i a l d i f f e r e n t i a l equation which describes the population dynamics, are compared with respect to their applicability, accuracy, efficiency and robustness. The method of lines transforms the partial differential equation into a set of ordinary differential equations. The Lax-Wendroff technique uses a finite difference approximation, to estimate both the derivative with respect to time and size. The remaining two are based on the method of characteristics. It can be concluded that the method of characteristics with a fixed time grid, the LaxWendroff technique and the transformation method, give satisfactory results in most of the applications. However, each of the methods has its own particular draw-back. The relevance of the major problems encountered are dicussed and it is concluded that the best method to be used depends very much on the application. The observed t r a n s i e n t s o f the c r y s t a l s i z e d i s t r i b u t i o n (CSD) o f i n d u s t r i a l c r y s t a l l i z e r s are e i t h e r caused by p r o c e s s d i s t u r b a n c e s o r by i n s t a b i l i t i e s i n the c r y s t a l l i z a t i o n process i t s e l f (1). Due to the introduction o f an on-line CSD measurement technique (2), the c o n t r o l o f CSD*s i n c r y s t a l l i z a t i o n p r o c e s s e s comes i n t o sight. Another requirement to reach t h i s g o a l i s a dynamic model f o r the CSD i n i n d u s t r i a l c r y s t a l l i z e r s . The dynamic model f o r a continuous c r y s t a l l i z a t i o n process consists o f a nonlinear p a r t i a l d i f f e r e n c e e q u a t i o n coupled to one or two ordinary d i f f e r e n t i a l equations (3.tiL) and i s completed by a set o f algebraic r e l a t i o n s f o r the growth and nucleation k i n e t i c s . The k i n e t i c r e l a t i o n s are empirical and contain a number o f p a r a m e t e r s w h i c h have t o be e s t i m a t e d from t h e experimental data. S i m u l a t i o n o f the experimental data i n combination with a n o n l i n e a r parameter e s t i m a t i o n i s a p o w e r f u l l technique to determine the k i n e t i c parameters from the experimental

0097-6156/90A)438-0159$06.00A) © 1990 American Chemical Society

In Crystallization as a Separations Process; Myerson, A., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1990.

CRYSTALLIZATION AS A SEPARATIONS PROCESS


160

data under instationary conditions (5.). Experimental d a t a from our c r y s t a l l i z e r , indicate that for a correct description of the CSD of the c r y s t a l l i z e r , a simple model with s i z e independent c r y s t a l growth i s not s a t i s f a c t o r y and t h a t o t h e r phenomena, l i k e growth dipersion and s i z e dependent growth or c r y s t a l a t t r i t i o n have to be i n c o r p o r a t e d i n t o the model. Implementation of these phenomena i n the s i m u l a t i o n however i n c r e a s e s t h e r e q u i r e d CPU t i m e f o r s i m u l a t i o n , and more parameters are to be e s t i m a t e d . N o n l i n e a r multiparameter o p t i m i z a t i o n techniques combined w i t h s i m u l a t i o n r e q u i r e s a l a r g e amount o f computing power. I t i s t h e r e f o r e important to choose the o p t i m a l s i m u l a t i o n a l g o r i t h m f o r t h e parameter e s t i m a t i o n program. In t h i s paper a comparison between four a l g o r i t h m s , which have been implemented and t e s t e d i n our simulation programs i s presented. I n d u s t r i a l c r y s t a l l i z e r s are o f t e n equiped w i t h an a n n u l a r zone f o r f i n e s removal. T h i s fines removal system gives r i s e to a mass accumulation i n the large annular zone (4_), which a f f e c t s the p r o c e s s dynamics a t l e a s t under unsteady conditions. I t i s shown, t h a t implementation o f t h i s f i n e s d e s t r u c t i o n system i n t o t h e simulation program has some implications for the algorithm used. The i n f l u e n c e o f the mass accumulation i n the a n n u l a r zone on t h e process dynamics i s also discussed. The Model In t h i s s e c t i o n the model for a continuous evaporative c r y s t a l l i z e r i s discussed. The c r y s t a l l i z e r i s o f the d r a f t tube b a f f l e d (DTB) type and i s equiped with a f i n e s removal system c o n s i s t i n g of a large annular zone on the outside o f the c r y s t a l l i z e r (see F i g u r e 1) . In o r d e r to v a r y the dissolved fines flow without changing the c u t - s i z e of the fines removal system, the flow i n the a n n u l a r zone i s kept constant and the flow i n the d i s s o l v i n g system i s varied by changing the recycle flow r a t e . The model assumptions are: - The c r y s t a l l i z e r i s well mixed. - The c r y s t a l l i z e r has a constant temperature and volume. - Only secondary nucleation takes place at s i z e L=0. - The feed flow i s c r y s t a l free. - A l l fines are dissolved i n the dissolved fines flow. - The f i n e s flow i n the a n n u l a r zone i s a p l u g f l o w and t h e c l a s s i f i c a t i o n o f the fines flow i s immediate at the entrance of the annular zone. - The removal of the fines and the product i s s i z e dependent and can be described by the product o f a s i z e dependent c l a s s i f c a t i o n f u n c t i o n h ( L ) and t h e p o p u l a t i o n d e n s i t y n ( L , t ) i n t h e crystallizer. The model f o r t h i s c r y s t a l l i z e r c o n f i g u r a t i o n has been shown to consist of the well known population balance (4_), coupled w i t h an ordinary d i f f e r e n t i a l equation, the concentration balance, and a set of algebraic e q u a t i o n s f o r the vapour flow r a t e , the growth and n u c l e a t i o n k i n e t i c s (4_). The p o p u l a t i o n balance i s a f i r s t - o r d e r hyperbolic p a r t i a l d i f f e r e n t i a l equation:


12. KRAMER ETAL.

Dynamic Behavior of Continuous Crystallizers


t •—Vapour-zone Dissolving vessel

t Phf«» , Of

Annular-zone

. Phc

Figure 1. Schematic picture of the continuous, evaporative crystallizer with fines removal system.


161


162

L

3n(L,t) 3(G(L,t)n(L,t)) at aL Q ( t ) . h (L) +

n(L,t-x)

v

+ K

V * > y > * W "

L

> '

K

t)

9

-

=0

(1)


with the boundary condition n(0,t) = B(t)/G(0,t) with G the c r y s t a l growth r a t e , B the nucleation r a t e , n the population density, Q the volume flow r a t e , h the c l a s s i f i c a t i o n function, V the c r y s t a l l i z e r volume.

(2) [m/s] [#/(m s)] [#/m*] [m /s] 3

3

3

[m ]

The i n d i c e s r e f e r to the f i n e s ( f ) , the product (p) and the fines recycle (r) flows. Note that part of the removed f i n e s flow i s f e d back i n the c r y s t a l l i z e r with a delay T , which depends on the volume of the annular zone and on the fines flow r a t e . The n u c l e a t i o n r a t e i s detemined by the w e l l known power law (3.»^U • c l a s s II system the growth rate i s calculated from the concentration b a l a n c e F

o

r

a

(i). 0

/ G ( L , t ) n ( L , t ) L * d L = [Q (t)C (t)+Q (t)C (t)+Q (t)C (t) i

i

d

d

r

r

-

C{Q (t) Q -Q (t)(l-c (t))}]/(p -C)3k V p

with

+

f

r

r

c

v

(3)

3

C p c k

the concentration, [kg/m ] the density, [kg/m ] the f r a c t i o n c l e a r l i q u i d of the s l u r r y , a shape f a c t o r . 3

v

The i n d i c e s i , d , c r e f e r to the feedflow, the dissolved fines flow and c r y s t a l substance respectively. The energy balance, then gives a r e l a t i o n f o r the vapour flow r a t e . The solution of the population balance i s strongly a f f e c t e d by growth r a t e f u n c t i o n used i n the m o d e l . I t has been shown t h a t s i z e dependent growth r a t e i s important to describe the d i f f u s i o n l i m i t e d c r y s t a l growth. Moreover i t has been shown that the s i z e reduction by a t t r i t i o n can also be described as a n e g a t i v e s i z e dependent growth r a t e (6.). F o r the s i m u l a t i o n s described here two types of s i z e dependent growth rates are used: G ( L , t ) = G(0,t)( k

for

1 + aL)

(4)

simple s i z e dependent c r y s t a l growth and G(L,t) = G ( L , t ) + G ( L , m , t ) k

a

3

to describe €he c r y s t a l a t t r i t i o n .


(5)

12.


KRAMER ETAL.


Simulation Algorithms. Simulation of the dynamic behaviour of i n d u s t r i a l c r y s t a l l i z e r s i s mainly concerned w i t h the s o l u t i o n o f the p a r t i a l d i f f e r e n t i a l equation, which i s solved e i t h e r d i r e c t l y o r by a t r a n s f o r m a t i o n method. Four a l g o r i t h m s have been i n c o r p o r a t e d i n our simulation program. Here a comparison w i l l be given of the d i f f e r e n t algorithms i n terms of: applicability e f f i c i e n c y with respect to c a l c u l a t i o n time accuracy robustness F i r s t the d i f f e r e n t methods w i l l be discussed b r i e f l y . Method o f l i n e s . T h i s m e t h o d , w h i c h was i n t r o d u c e d i n the c r y s t a l l i z a t i o n research by Tsuruoka and Randolph (J) , transforms the population balance i n t o a set of ordinary d i f f e r e n t i a l equations by d i s c r e t i z a t i o n of the s i z e axis i n a fixed number of g r i d p o i n t s . The d i f f e r e n t i a l 3 { G ( L , t ) n ( L , t ) } / 3 L i s then approximated by a f i n i t e d i f f e r e n c e scheme. As has been shown (J) a f o u r t h - o r d e r accurate scheme using an equidistant g r i d (i= l . . . . z ) r e s u l t s i n : 3F dL

F

i-2

-

8 F

i-l

(6)

+8F. , - F . )/12AL + O(AL) l+l 1+2'' ' 0

x

where F^= G(L t)n(L t). At the boundaries t h i s c e n t r a l difference E q u a t i o n 6 cannot be used because the g r i d points L _ ^ , L Q , L » ^ and , do not e x i s t . A f i v e point Lagrange i n t e r p o l a t i o n i s used to e i t i i mate these g r i d points, leading to the following equations: i-2. 3F -6F F . ) / 12AL = (-3F -10F l 8 F (7) dL i=z-l, 3F (8) (-F 6F _ -18F 1 0 F + 3 F ) / 12AL 3L ±9

it

z +

T

U

±

± - 3

i=z,

3F dL

( + 3 F

+

±

i-4 "

l 6 F

+

i + 1

2

± + 2

± - 1

+ 3 6 F

i-3

+

i-2

+

+ 3

±

±+1

-48F _ ±

+25F ) / 12AL ±

1

(9)

w i t h F = G ( 0 , t ) n ( 0 , t) . As t h e p o p u l a t i o n d e n s i t y a t L = 0 i s determined by the boundary condition, the scheme leads to a set of z-1 ordinary d i f f e r e n t i a l equations, which are solved u s i n g a Runge K u t t a a l g o r i t h m . Due to the fixed simulation g r i d and the numerical approach, t h i s scheme i s general applicable i n the sence, that there are no l i m i t a t i o n s i n the c r y s t a l l i z e r configuration or the c r y s t a l growth r a t e f u n c t i o n . However, t h e method i s s e n s i t i v e t o d i s c o n t i n u i t i e s i n Gn as i s shown i n the next s e c t i o n . To examine these o s c i l l a t i o n s i n more d e t a i l , also lower order versions of t h i s scheme have been implemented. The second-order version i s given by: 3F dL L . 3F 3L L /

=

+

/ ILIL^ " 2TL< i l , t - i , t >

(

1

5

)

(

1

6

)

+

F

F

i+. .t+.5 " i - . 5 , f . 5 > 5

This implementation i s s e c o n d - o r d e r a c c u r a t e w i t h r e s p e c t to the time and the s i z e s t e p . The scheme i s g e n e r a l a p p l i c a b l e and as shown i n the next s e c t i o n , t h i s scheme i s a l s o s e n s i t i v e f o r d i s c o n t i n u i t i e s i n Gn as caused by the R-Z model for fines removal. The o s c i l l a t i o n s are however l e s s severe than f o r the method o f l i n e s . A l s o f o r t h i s method a f i r s t - o r d e r scheme was implemented. Here the s o - c a l l e d Lax scheme was chosen (8): n

j.

t +

i 4 L*

(20)

T h i s removal f u n c t i o n g i v e s r i s e t o a d i s c o n t i n u i t y i n t h e p o p u l a t i o n d e n s i t y a t the c u t s i z e o f the f i n e s . The n u c l e a t i o n parameters are given i n equation 19. In Figure 3 the responses a r e shown o f the p o p u l a t i o n d e n s i t y at 120 pm and of the growth rate after a step i n the heat input to the c r y s t a l l i z e r from 120 to 170 kW for three simulation algorithms. The cut-size of the fines L was 100 jim, a s i z e dependent growth r a t e was used as d e s c r i b e d by Equation 4 with a= -250 and the number of g r i d points was 400. When the s i m u l a t i o n was performed with the method o f l i n e s , s e v e r e o s c i l l a t i o n s are present i n the response of the population density at 120 pm, which dampen out rather slowly. Also the response o f the Lax-Wendroff method shows these o s c i l l a t i o n s to a lesser extend. f



12.

KRAMER ETAL.

Dynamic Behavior of Continuous Crystallizers167

Figure 2. Response of the growth rate G on a negative step in the heat input of the crystallizer simulated with the Lax-Wendroff (line) and the MCFS (dots) algorithm.




168

Figure 3. Response of the population density at 120 /im (top) and of the crystal growth rate (bottom) on a step in the heat input of the crystallizer simulated with method of lines (line), Lax-Wendroff (dots), and the MCFS (dashes) algorithm.



12.

KRAMER ETAL.


This unstable behaviour i s not seen i n the r e s p o n s e s o f the i n t e g r a t e d p r o c e s s v a r i a b l e s l i k e the moments and the growth rate (see Figure 3B). In the responses of the two implementations o f the methods o f c h a r a c t e r i s t i c s , no u n s t a b l e behaviour was found. The response of the MCFS, however again shows the dynamic e r r o r , which r e s u l t s from the v a r i a b l e timestep. The delay i n the responses could almost completely be removed by a d o u b l i n g o f the number o f g r i d p o i n t s (not shown). To investigate whether the o s c i l l a t i o n s i n the response of the population density a f t e r a s t e p i n the heat i n p u t are caused by the o r d e r o f accuracy of the algorithm, simulations were performed using lower order implementations o f the method o f l i n e s and of the Lax-Wendroff algorithm. The conditions were s i m i l a r to those i n Figure 3» The r e s u l t s are given i n Figure 4 and 5» while i n T a b l e I I I the stationary values of m and the growth rate before and 35000 seconds after the step i n the neat i n p u t are g i v e n . The r e s u l t s i n d i c a t e t h a t t h e r e i s a c l e a r r e l a t i o n s h i p between the i n t e n s i t y of the o s c i l l a t i o n s and the order of the algorithm. In the r e s p o n s e o f t h e f i r s t - o r d e r method o f l i n e s s i m u l a t i o n s no o s c i l l a t i o n s are p r e s e n t . However the response ( F i g . 4) and the s t a t i o n a r y v a l u e s (Table I I I ) d e v i a t e somewhat from the higher o r d e r s o l u t i o n , which cannot be compensated c o m p l e t e l y by an increase i n the number of g r i d points used, as i s shown i n Figure 6. The f i r s t - o r d e r Lax scheme g i v e s too h i g h d e v i a t i o n s and i s not reliable. The s t a t i o n a r y values as given i n Table I I I give a much higher spread than for the MSMPR c r y s t a l l i z e r , which i s due to a combination o f the d i s c o n t i n u o u s f i n e s removal f u n c t i o n and the difference schemes used i n the d i f f e r e n t methods. The r e l e v a n c e o f the o s c i l l a t i o n s seen i n the method of l i n e s and the Lax-Wendroff techniques are not easy to d e s c r i b e . In most cases one i s n o t i n t e r e s t e d i n the population density at a s i n g l e c r y s t a l s i z e , but i n integrated proces v a r i a b l e s l i k e the number o f c r y s t a l s i n a c e r t a i n s i z e range, the moments or the c r y s t a l growth r a t e , i n which these o s c i l l a t i o n s a r e n o t p r e s e n t . On t h e o t h e r h a n d , t h e o s c i l l a t i o n s w i l l affect the accuracy of these responses. The lower order schemes of the method of l i n e s form a reasonble a l t e r n a t i v e to decrease the e f f e c t o f the unstable behaviour. The two methods of c h a r a c t e r i s t i c s do not show t h i s u n s t a b l e b e h a v i o u r . However, the MCFS technique suffers from large dynamic errors and i s l i m i t t e d i n i t s a p p l i c a t i o n . The MCFT seems to form a good a l t e r n a t i v e . fi

As d i s c u s s e d i n the previous section and summarized i n Table I I , a drawback of the MCFT method i s that the mass a c c u m u l a t i o n i n the f i n e s removal system cannot be simulated, therefore we examined whether t h i s mass accumulation has a noticable e f f e c t on the process dynamics. In the s i m u l a t i o n the fines removal i s simulated with a cut s i z e s of 150 pm. The f i n e s flow r a t e Q and the r e c y c l e flow r a t e Q were 1.25 and .75 l i t e r per second. The r e s u l t s are shown i n Figure 7. I t i s c l e a r that the mass a c c u m u l a t i o n has indeed an e f f e c t on t h e p r o c e s s dynamics. Even on the mean s i z e o f the c r y s t a l s a c l e a r s h i f t i n the response i s seen. I t appears t h a t the e f f e c t i s s t r o n g l y dependent on the value of the recycle flow rate (not shown). The conclusion from these r e s u l t s i s t h a t the e f f e c t s o f mass accumulation i n the fines system are present, and can only be neglected at low cutsizes and low fines recycle r a t e s . f

p


169


170


Figure 4. Response of the population density at 120 /*m on a step in the heat input of the crystallizer simulated with method of lines, fourth-order (line), second-order (dots), and first-order (dashes).

Figure 5. Response of the population density at 120 /im on a step in the heat input of the crystallizer simulated with Lax-Wendroff, second-order (line), and the first-order Lax scheme.



KRAMER ETAL.

Dynamic Behavior of Continuous Crystallizers171

Figure 6. Effect of the number of grid points on the response of the population density at 120 /im on a step in heat input; fourth-order (line), first-order 400 grid points (dots), and first-order 800 grid points (dashes).




172

Figure 7. Effect of the mass accumulation on the response of the L ^ (A) and on that of the growth rate G and the nucleation rate B (B) on a step in the heat input of the crystallizer. The cutsize of the fines was 150 pm.


12. KRAMER ETAL.


Table I I I . S t a t i o n a r y v a l u e s o f G a n d m f o r a DTB c r y s t a l l i z e r operated with a fines removal system and a s i z e dependent growth regime. The values are given before and 35000 second a f t e r a s t e p i n the heat input of the c r y s t a l l i z e r from 120 to 170 kW. The s i m u l a t i o n has been p e r f o r m e d w i t h 400 gridpoints 2


Order

Method of Lines tt tt tt tt Lax-Wendroff tt tt MCFS MCFT

4 2 1 2 1 1 1

before step mO *1.E10 G *l.E-7

.5620 .5624 .5688 .5628 .5332 .5587 .5582

.8416 .8381 .8531 .8370 .8740 .8368 .8361

35000 s a f t e r step G *l.E-7

mO *1.E10

.5363 .5364 .5226 .5369 .5055 .5328 .5294

1.491 1.485 1.5H 1.484

1.552 1.482 1.478

Conclusions. F o u r s i m u l a t i o n a l g o r i t h m s have been d i s c u s s e d i n terms o f a p p l i c a b i l i t y , accuracy, c a l c u l a t i o n time and robustness. Two of the methods d e s c r i b e d are based on numerical methods, while the other two are based on the method of c h a r a c t e r i s t i c s . To simulate a simple MSMPR c r y s t a l l i z e r , a l l four algorithms give reasonable r e s u l t s . The methods were a l s o a p p l i e d t o more c o m p l e x c r y s t a l l i z e r c o n f i g u r a t i o n s and on d i f f e r e n t growth rate k i n e t i c s . Results show that the best method to be used depends on the a p p l i c a t i o n , w h i l e a l l the methods s t u d i e d are unsatisfactory i n one or more aspects s t u d i e d . The most important drawbacks o f t h e method o f l i n e s a l g o r i t h m are i t s long c a l c u l a t i o n time and the unstable behaviour i n the response o f p o p u l a t i o n d e n s i t y i n the p r e s e n c e o f d i s c o n t i n u i t i e s i n G n . Lower o r d e r implementations can be an a l t e r n a t i v e . In t h i s case more g r i d p o i n t s w i l l be n e c e s s a r y to maintain the required accuracy. The Lax-Wendroff method i s l e s s a c c u r a t e , but i s f a s t e r than the method o f l i n e s . The scheme i s also sensitive f o r d i s c o n t i n u i t i e s i n Gn, although the o s c i l l a t i o n s are less severe. The f i r s t - o r d e r Lax scheme i s not r e l i a b l e . The MCFS method i s simple to use and f a s t , The technique i s however not general applicable and shows severe dynamic e r r o r s . The MCFT algorithm i s the most favorable of the s t u d i e d s i m u l a t i o n methods. The only draw-back of the method i s that i t cannot be used to s i m u l a t e the e f f e c t s o f the mass accumulation i n t h e f i n e s removal system. In t h i s paper i t has been shown however, that the effects of the mass accumulation i n the fines removal system on the p r o c e s s dynamics can not be neglected, unless low a cutsize f o r the fines removal and low fines recycle rates are used.


173


174 Acknowledgments

The authors would l i k e to express t h e i r gratitude to the Netherlands Technology Foundation (STW), AKZO, DOW, DSM, Dupont de Nemours, Rhone Poulenc , Suiker Unie and Unilever for t h e i r f i n a n c i a l support of the reseach program.


References 1. Randolph, A.D., Beckman, J.R., Kraljevich, Z.I., (1977) AIChE J.,23, 500-510 2. Jager, J., De Wolf, S., De Jong, E . J . and Klapwijk, W. in Industrial Crystallization 87 , Nyvilt, J., Zacek, S. Eds, Elsevier Amsterdam, 1989 p 415 418. 3. Randolph, A.D. and Larson, M.A. Theory of particulate processes; analysis and techniques of continuous crystallization; Academic Press New York, 1971. enz 4. De Wolf, S., Jager, J., Kramer, H.J.M., Eek, R. and Bosgra, O.H. Proc. IFAC Symp. on Dynamics and Control of Chemical Reactors, Destination Columns and Batch Processes, Maastricht, The Netherlands, 1989. 5. Frank, R., David, R, Villermaux, J . and Klein, J.P. Chem. Engng. Sci., 43, 69-77 6. Polisch, J . and Mersman, A. Chem. Engng. Sci. Technol.,1988, 11, 40-49 7. Tsuruoka, C. and Randolph, A.D., AIChE Symposium series, no. 253,1987, vol.83, pp. 104-109 8. Vemuri, V. and Karplus, W.J. (1981) Digital computer treatment of partial differential equations. Prentice Hall, Engelood C l i f f s , New Jersey. 9. De Leer, B. G. M., Ph.D. Thesis Delft University of Technology, Delft 1981. RECEIVED May 12, 1990


Simulation of the Dynamic Behavior of Continuous Crystallizers - ACS

Recommend Documents