Chemical Engineering Science On the prediction ... - Semantic Scholar
Chemical Engineering Science 64 (2009) 686 -- 696
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s
On the prediction of crystal shape distributions in a steady-state continuous crystallizer Christian Borchert a , Nandkishor Nere b , Doraiswami Ramkrishna b,∗ , Andreas Voigt c , Kai Sundmacher a,c a b c
Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, D-39106 Magdeburg, Germany Purdue University, School of Chemical Engineering, 480 Stadium Mall Dr., West Lafayette, IN 47907, USA Otto von Guericke University, Process Systems Engineering, PSF 410, D-39106 Magdeburg, Germany
A R T I C L E
I N F O
Article history: Received 26 November 2007 Received in revised form 27 March 2008 Accepted 15 May 2008 Available online 22 May 2008 Keywords: Crystallization Modeling Morphology Multidimensional population balance
A B S T R A C T
Zhang and Doherty [2004. Simultaneous prediction of crystal shape and size for solution crystallization. A.I.Ch.E. Journal 50, 2101–2112] have provided a one-dimensional analysis of crystallization based on the assumption that the relative face-specific growth rates of a (2-D) crystal are independent of supersaturation and hence invariant with time. Subsequent work by these authors [Zhang, Y., Sizemore, J.P., Doherty, M.F., 2006. Shape evolution of 3-dimensional faceted crystals. A.I.Ch.E. Journal 52, 1906–1915) consider shape evolution of single three-dimensional crystals with morphological changes. In this work, we present a multidimensional population balance approach accounting for dependence of the relative face-specific growth rates on supersaturation, a situation more commonly encountered. For example, Joshi and Paul [1974. Effect of supersaturation and fluid shear on habit and homogeneity of potassium dihydrogen phosphate crystals. Journal of Crystal Growth 22, 321–327] and Mullin and Whiting [1980. Succinic acid crystal-growth rates in aqueous solution. Industrial & Engineering Chemistry Fundamentals 19, 117–121] report face-specific growth rates with different dependence on the supersaturation. Thus it has been observed that there exists significantly different crystal shapes in a crystallizer [Yang, G., Kubota, N., Sha, Z., Louhi-Kultanen, M. Wang, J., 2006. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Crystal Growth and Design 6, 2799–2803]. Consequently, the population of crystals at any instant will have widely varying crystal shapes and sizes depending upon the initial crystal shape and size distribution. Computations are presented for the shape distributions of the crystal population emerging from a steady-state continuous crystallizer for two cases: (1) feed without crystals including nucleation for the formation of new crystals, and (2) feed with seed crystals of known shape, with suppressed nucleation. In the range of mean residence times investigated, the calculated crystal volume distributions for the first case show geometrically dissimilar shapes without morphological variations. However, in the second case, because the feed crystals of the chosen shape were susceptible to morphological changes, the volume distributions display this feature with shape and size distributions for each of a number of different morphologies. By varying operating conditions such as the flow rate, the inlet supersaturation, and the shapes of feed crystals, the proposed model can clearly be used to manipulate the crystal shape and size distributions and their morphologies. © 2008 Elsevier Ltd. All rights reserved.
1. Introduction We embark in this paper on an analysis of crystal populations in the habitat of a crystallizer, for the latter is the source of crystalline products in industrial practice, where the shape of faceted crystals is a major quality factor. In the chemical industries in general and in the speciality and fine chemical industries, in particular, the crystal shape and its distribution are of considerable interest because it can remarkably influence the product downstream processing and
∗
Corresponding author. Tel.: +1 765 494 4066; fax: +1 765 494 0805. E-mail address: [email protected] (D. Ramkrishna).
0009-2509/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.05.009
usability. For industrial applications it is important to predict, optimize and control shape distributions. Since the growth of single crystals depend on supersaturation in the environment, which must change due to the cumulative growths of all crystals around, the analysis must account for the reciprocal effects of crystal growth and supersaturation dynamics. Indeed population balances represent a natural implement for such an analysis. Since the potential complexity of a full multidimensional approach calls for some sense of parsimony in regard to the extent of detail, it behooves us to begin with an examination of scenarios for simplification. In this connection, the work of Zhang and Doherty (2004) is a landmark that represents in some way the most desirable simplification, viz., that crystal
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
growth may be described in a unidimensional space of one kind or another. We argue that, in describing crystal growth the authors have stretched the most out of a one-dimensional framework and that no improvement is possible without transgression into higher dimensions. More specifically, the authors assume that the relative face-specific growth rates are independent of the supersaturation and hence of time. This assumption traces the evolution of crystal growth to be along a rectilinear path in crystal shape space as determined by the vector, h ≡ (h1 , h2 , . . . , hn )T , representing distances of the various faces from the center of the crystal. This is a line crystal shape space, described by some parameter along it, which could either be the volume, mass or age of the crystal or a component of belonging to a secure face. Strict unidimensionality of a crystal population must confine all crystals to this line, which must necessarily pass through the origin of shape space. Clearly, all crystals lying on this line must be geometrically similar, characterized by a scalar property of one kind or another. Thus no new crystal formed by nucleation outside of this line will be geometrically similar to the set of crystals on the line. As Zhang et al. (2006) observe, faces of a crystal can appear or disappear during the course of growth signifying a morphological change. This change in morphology is indeed an interruption of its geometric similarity with its past although crystal growth will persist along the same direction with respect to hi 's representing unchanged faces. There are thus two kinds of outcasts from the geometrically similar crystals along the line of unidimensionality defined by constant relative face growths. The first are those that are born by nucleation outside of the line, and the second are those that undergo morphological changes. Strict unidimensionality cannot of course accommodate either kind. With the assumption of constant relative face growths, unidimensionality is still possible in an approximate sense. Fresh crystals formed by nucleation outside of the one-dimensional line can never be on it as it must proceed parallel to it from outside it. However, the shape of a growing crystal, as long as it does not change morphology, will become progressively insensitive to its initial perturbation from the line. Thus outcasts of the first kind, in course of time, may be regarded as geometrically similar in an approximate sense. Outcasts of the second kind also can similarly be inducted into the approximately similar set if morphological transitions terminate. Zhang and Doherty (2004) show that the geometric area factor of crystals becomes eventually constant as the population evolves through various morphological changes into an asymptotic shape. At this stage, a one-dimensional description of the crystal population becomes possible as the scalar property can be inverted to construct the full shape of each crystal. The discussion above must indicate that Zhang and Doherty (2004) have extracted the most out of a one-dimensional model. To seek improvement, one must admit multidimensionality (a fact recognized by Zhang and Doherty, 2004), making full use of crystal symmetry to minimize the increase in dimension. Before we set out to address our approach, however, we note one additional feature of the work of Zhang and Doherty (2004). While the assumption of constant relative face growth confines crystal growth to a line independently of supersaturation, the speed with which crystals move is regulated by supersaturation so that coupling between the two processes is accommodated. Our multidimensional framework, however, couples supersaturation dynamics intimately with growth of the crystal population. Crystal growth occurs along curved lines in shape space. Due to the increase in complexity, we confine our effort to relatively simple situations. For our treatment, we assume that the crystals are surrounded by a continuous phase in which physico-chemical processes occur leading to supersaturation. This work will ignore the path to supersaturation by assuming its direct availability.
687
Due to crystal growth, dissolved material is transferred from the continuous phase to the dispersed phase, more specifically: to the crystal faces. We assume that the crystal is built of only one species. The driving force for growth is the supersaturation. In the subsequent text we use the relative supersaturation defined by
=
− sat , sat
(1)
where is the mass concentration of the solute while sat its saturation concentration. The paper is organized as follows. In Section 2 we introduce a single crystal shape model along the lines of the work of Zhang et al. 2006). However, this work departs in one significant way from that of Zhang and Doherty (2004) in that the relative growth of different crystal faces, which was assumed constant in earlier work, is allowed to depend on supersaturation. The coupling of the single crystal with the continuous phase is discussed within the framework of population balance which accounts for the joint effects of the growth of all crystals in the population. Furthermore, we establish the mass balance in terms of supersaturation. Section 3 shows how one obtains a conditional equation for the steady-state supersaturation of a mixed suspension mixed product removal (MSMPR) crystallizer. In Section 4, the results of model simulations using the conditional equation for the steady-state supersaturation are discussed. The supersaturation description has been combined with (i) the classical age distribution of ideally mixed systems and (ii) with the single crystal shape evolution model from Section 2 to obtain shape distributions. In the first example, we consider an arbitrary cuboid crystal, which is a rather simple system while the second example shows how the developed tool can be applied to a more complex case (paracetamol) using the growth rate measurements by Ristic et al. (2001). Finally, we conclude by summarizing the outcomes of the present work in brief. 2. Model formulation Crystallizers contain a large number of crystals, called a crystal population. The major task in crystallization is to grow crystals. Therefore we focus on growth throughout the paper and neglect breakage and aggregation. The framework that allows the investigation of the dynamics of a crystal population is the population balance equation. However, in order to describe the behavior of a crystal population properly one should fully understand the dynamics of a single crystal. Thus we introduce at first the shape model for a single crystal followed by the discussion of some of the important properties. Finally, the population balance model for the suspension is introduced. 2.1. Shape model We consider crystals with a convex shape composed of n flat faces. A point O is chosen, for example the crystal center, to be the origin of a Cartesian coordinate system in which the crystal form is described. The natural starting point to decompose a convex crystal is to consider each face separately. The orientation of the ith face is given by its unit normal ni , which is obtained from transforming Miller indices to normal vectors in Cartesian coordinates. The distance of the ith face to the origin O is denoted by hi . The ith face of a crystal is a finite sector on the plane ni · r = hi ,
i = 1 . . . n,
(2)
where r is a point in 3-D physical space. Clearly, hi is a property of the crystal that serves to describe its geometry. Thus hi is interpreted as the ith entry of the geometrical state vector h = [h1 . . . hn ]T . The geometry of a crystal is
688
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
altered when h changes. Therefore we will consider the change of h due to face-specific crystal growth. Note that the orientation ni of a face does not change due to growth but only hi . Microscopically, a crystal grows by the attachment of growth units to its surface. We regard the growth units to be orders of magnitude smaller than the crystal so that the discrete statistical processes of crystal growth can be approximated by a continuous model. The growth of a crystal is usually anisotropic, i.e. the growth rates of the faces follow different laws, see e.g. Sangwal (1998). The dynamics of h(t) is described by the vector differential equation: dh ˙ = H(h, ), dt
h(0) = h0 ,
(3)
˙ = [H ˙ n ]T , and H ˙ is the perpendicular growth rate of the ith ˙ , ... , H H 1 i ˙ is also called the face velocity and is face, see Zhang et al. (2006). H i in general a function of the crystal state and supersaturation . The crystal shape S is defined by the planes containing crystal faces as
S = {r : ni · r < hi , i = 1 . . . n},
(4)
which is a convex polyhedron. The shape evolution is given by t ˙ S(t) = r : Nr < h(t), h(t) = h0 + )), ()) d , H(h( (5) 0
where ⎡ T⎤ n1 ⎢ . ⎥ ⎥ N=⎢ ⎣ .. ⎦
For many crystals we find, at least for some idealized situation a symmetrical crystal shape. Obviously it is then sufficient to consider only one of the symmetrical faces. Hence, we can reduce the dimension of the geometrical state vector h remarkably when we exploit crystal symmetry. There may be other properties of interest, for example the crystal age, color, density or composition. In this paper we will have only one additional property: crystal age . Hence, the full property state vector is x = [hT ]T .
(6) 2.2. Dynamics of real and virtual faces
nTn is the n by three matrix containing the unit normals of the faces. Shapes S1 and S2 are defined to be similar when hS = shS , 2 1 where s is a scalar. Usually, the shapes at different instants are not ˙ the influsimilar due to anisotropic effects in growth. For t?hi,0 / H i ence of the initial geometry h0 on h is negligible compared to the influence of the growth. Further, if > 0 is at steady state, denoted ˙ does not by ss , as, for example, in a continuous crystallizer and if H i ˙ depend on the geometry (jHi / jhj = 0), we find by replacing time t with age , a variable more appropriate for the dynamics of individual ˙ . That means s = crystals in a steady-state crystallizer that h() = H 2 / 1 is a scaling factor to relate crystal shapes at different instants: h(1 ) = sh(2 ),
Fig. 1. In case that the face with the thick line has a remarkably higher growth rate than the neighboring faces, its surface area becomes lower until it grows out (top) or its increases (bottom).
(7)
which means that the shapes are similar. If successive shapes are similar under certain growth conditions (for instance supersaturation) then we speak of the steady-state geometry. It is opportune at this stage to discuss the precise difference between this work and that of Zhang and Doherty (2004). In the discussion below, we will refer to the space to which the vector h belongs as the crystal shape space or simply shape space. Zhang and Doherty (2004) assume the relative growths of different faces to remain constant in time regardless of the supersaturation. Thus their view traces the evolution of shape along a rectilinear trajectory in shape space, independently of supersaturation dynamics in the continuous phase. Of course insofar as supersaturation does influence the growth of each face it dictates the speed with which shape evolution occurs without changing its direction. In contrast, in our work supersaturation not only determines the speed with which shape evolution occurs but also its direction. Thus in the steady state envisaged above although the present work will feature crystals that are geometrically similar, the shape and size will vary with residence time due to variation in the steadystate supersaturation. On the other hand, Zhang and Doherty (2004) will only perceive differences in the shape for a single crystal. Their work ignores the feedback of crystal growth on supersaturation.
Since crystal faces grow at different rates, it can happen that one face is much faster than its neighbors so that it finally grows out, see Fig. 1. Therefore one generally has only n < n faces present on the crystal surface, whereas there are n = n − n faces which are virtual (likely faces). In the following we assign to the distance and unit normals of virtual faces a double-prime to identify them as such, for example hi and ni . In case that we explicitly refer to a real face we provide the distance and unit normal with a single prime, hi and ni . We stipulate that a virtual face is always tangential to the crystal, see Fig. 2. Since no two faces have the same orientation, the polyhedron S has an extreme point Vijk in the direction n of the virtual face . The point Vijk is a vertex defined by the intersection point of at least three real faces i, j, k. The coordinates of the intersection are obtained by vijk = N−1 h , ijk ijk
where N−1 = [ni nj nk ]T , ijk
hijk = (hi hj hk )T .
(8)
Note that there exists an extreme line (edge) Vijk Vijl on the crystal when n ·(vijk −vijl )=0. Then it is sufficient for our purposes to iden tify only one point on Vijk Vijl , e.g. Vijk to be the extreme point of S with respect to n . The distance h is dictated by vijk because the vir tual face is tangential to S. Then h is the projection of vijk onto the direction n , that is h = n · vijk . By the virtue of Eq. (8) the distance of the virtual face is a linear function of the distances of real faces: h = q · hijk ,
q = (n )T N−1 . ijk
(9)
Constructing an n -dimensional real geometrical state vector h which contains only the distances of the real faces and a virtual geometrical state vector h (n -dimensional), we can write h = Q h ,
(10)
where the n by n matrix Q contains on its rows the entries of q on the corresponding position. The matrix Q changes its entries and
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
689
Fig. 2. Different morphologies of the same crystal exhibiting different faces. Even though the faces which cut the edges of the octahedron (left) are not real (green faces) we have to keep track of them due to the fact that they may appear (right). We call them virtual faces (blue) and they are always tangential to the crystal shape.
dimension depending on the geometrical configuration of the crystal. It is especially necessary to track the number of real and virtual faces n and n . That is, for qualitatively different geometrical configurations as shown in Fig. 2 we find different matrices Q . An increase in dimension in the real part of the shape space (h space) induced by the conversion of a virtual face into a real one must be viewed as a step change in the morphology as it will signal a discontinuous change in the number of vertices, edges or planes. This new morphology, however, is deemed to remain invariant until the next step change occurs. ˙ =dh /dt, have to be determined from Growth rates of real faces, H i i crystal growth models and verified by experiments. We assume here ˙ is a function of supersaturation without loss of generality that H i only. There exist several crystal growth theories in the literature, see for instance Sangwal (1998) and Prywer (2005). The displacement rate of a virtual face is—like its position—a function of growth rates of real faces. Differentiation of Eq. (10) with respect to time yields the velocity vector of virtual faces: ˙ = Q H ˙ (), H
˙ = dh where H dt
and
˙ () = dh . H dt
(11)
An algorithm for the calculation of shape evolution when faces can change their state from real to virtual and vice versa can be found in Zhang et al. (2006). For the present, we abstractly observe that step changes can occur in the evolution when transitions occur between virtual faces/edges and real faces/edges. Such transitions are characterized by equalities involving the growth vectors and the motion of vertices, which together define some hyper surface in the joint space of supersaturation and crystal shape. These hyper surfaces define the domain of operation for restricting the variation in morphology. 2.3. Mass balance for a single crystal
For practical calculations we decompose the polyhedron into pyramids Si whose base plane is the ith face and apex is the origin of the coordinate system O. Pyramids are again decomposed into tetrahedra Sij whose vertices are the origin and the vertices of the ith face, see Fig. 3. We number the coordinates of the vertices of the ith face from vi1 , . . . , vi . Then we can find i − 2 tetrahedra per face. i The vertex coordinates of the jth tetrahedron of the ith face is stored in the matrix j = 1 . . . (i − 2).
(13)
n i −2
V(S) =
V(Sij ).
(15)
i=1 j=1
In the following we represent the volume of a crystal, V = V(S) by V(h). Assuming that the density of a crystal, s , is constant, the rate of change of the crystal mass due to growth, Gm is given by n
jV d ˙ ˙ = ( ∇ V) · H, H V(h) = s s h dt jhi i
(16,17)
i=1
where ∇h =(j/ jh1 , . . . , j/ jhn ). This mass transfer term will later be used in the mass balance for the continuous phase. We can calculate the volume change as follows. The change of the volume requires the evaluation of the derivative of the tetrahedra V(Sij ), see Eq. (14): 1 d d V(Sij ) = | det(Tij )|. dt 6 dt
(18)
The derivative of the norm of a function f is d|f | = sign(f ) df . For the derivative of a determinant of a matrix A we find (see Petersen and Pederson, 2006): d det(A) = det(A) tr(A−1 d A), where tr is the operator delivering the trace of the matrix. The derivative of Tij is denoted by T˙ and contains the vertex velocity vectors v˙ which are ij
ij
obtained from differentiation of Eq. (8):
The volume of the tetrahedron Sij is readily found to be V(Sij ) = 16 | det(Tij )|.
Summation over all tetrahedra of a real face and summation over all n real faces yields the volume of the crystal
Gm = s
The volume of a crystal or rather the polyhedron S is given by V(S) = dV. (12) S
Tij = [vi1 vi(j+1) vi(j+2) ],
Fig. 3. Decomposition of the volume under the ith face into tetrahedra.
(14)
T˙ ij = [v˙ i1 v˙ ij+1 v˙ ij+2 ]T ,
j = 1, . . . , i − 2.
(19)
690
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
Then we get 1 d V(Sij ) = sign(V(Sij )) det(Tij ) tr((Tij )−1 T˙ ij ) dt 6 1 = sign(V(Sij )) tr((Tij )adj T˙ ij ), 6
(20a) (20b)
where adj refers to the adjugate or adjoint (transpose of the cofactor).
3. Steady-state analysis
2.4. Multidimensional population balance equation So far we have discussed the single crystal model. In the following we account for a whole population of crystals. The crystal state vector x = [hT , ]T is a point in the continuous property state space R(n+1) . Here we consider the domain of state space, to be of our interest. Particles entering originate at one point xnuc ∈ j with the rate of nucleation B0 . It is assumed that all nuclei have the same state xnuc . Following Ramkrishna (2000), the population balance equation for a population of crystals undergoing convective growth in a wellmixed continuous crystallizer is
j˜ ˙ f˜ (x, t) = − 1 f˜ (x, t), f (x, t) + ∇x · X 1 jt 1 1
(21a)
I.C.: f˜1 (x, 0) = f˜1,0 ,
(21b)
˙ · nx )f˜ (x, t) = B (x − xnuc ), B.C.: (−X 1 0
x ∈ j,
R.C.: f˜1 (x → ∞, t) = 0,
for the supersaturation using the preceding definitions in connection with Eqs. (16), (23) and (24) in Eq. (22) becomes
1 d d = (in + 1)Fin − ( + 1) Fout − ( + 1) ( Vr ) dt Vr dt s ˙ · (∇ V) dV . − (f˜ H) (26) h h sat 1
We assumed in Eqs. (21) that all crystals nucleated in the crystallizer have the same state xnuc ∈ j. At steady-state operation, all crystals see the same environmental conditions, namely the (steadystate) supersaturation ss . Therefore all crystals undergo the same shape evolution according to Eq. (5). That is, all crystals evolve along the same trajectory h(, ss ) through shape space. The time which a crystal spends within the crystallizer is the age . Suppose we know the steady-state supersaturation ss , then we can calculate the trajectory h(, ss ) for an individual crystal using Eq. (3). That means we know the shape as a function of crystal age, S(, ss ), and thus the volume V(, ss ) for all crystal ages. We consider the population balance equation only for the crystal age since shape is then recoverable purely from . Integration of Eq. (21) over h yields:
(21c)
j j 1 f (, t) + f (, t) = − f1 (, t), jt 1 j 1
(27a)
(21d)
I.C.: f1 (, 0) = f1,0 ,
(27b)
B.C.: f1 (0, t) = B0 (),
(27c)
R.C.: f1 ( → ∞, t) = 0,
(27d)
where f˜1 is the crystal number density: [f˜1 ] = #/ i [xi ]/m3 . The crystallizer volume is Vr and is the mean residence time. Vector nx is an outer unit normal on the boundary j. It is obvious that for the case of continuous seeding the equation will require an additional source term on the right-hand side.
where f1 (, t) = f˜1 (h, , t) dVh . The steady-state age distribution h is given by
j 1 f () = − f1,ss (), j 1,ss
(28a)
The continuous phase is in general characterized by a state vector y but we restrict our analysis to mass m of the crystal solute. Let ˙ in be the mass flow rate of the saturated solution into and m ˙ out m ˙ disp be the mass be the mass flow rate out of the crystallizer, and m consumption of the dispersed phase. The mass balance follows:
B.C. : f1,ss (0) = B0 (),
(28b)
R.C. : f1,ss ( → ∞) = 0,
(28c)
dm ˙ in − m ˙ out − m ˙ disp . =m dt
f1,ss () = B0 () exp(−/ ),
2.5. Mass balance for crystal population
(22)
The mass transfer of crystal solute from the continuous phase to the crystalline phase is given by ∞ ˙ disp = Vr Gm (h)f˜1 (h, , t) dVh d, (23) m 0 where dVh is an infinitesimal volume measure in h-space and Vr is the total volume of the slurry. Vr is the sum of the volume of the continuous phase and the volume of the dispersed phase: Vr = Vc + Vdisp . The mass transfer rate Gm is given by Eq. (16). Furthermore we write m = Vr ,
˙ in = in Fin , m
˙ out = Fout , m
(24)
where is the mass concentration of dissolved material, Fin and Fout is the volumetric flow rate into and out of the crystallizer and is the volume fraction of the liquid phase defined as
=
Vdisp Vc =1− . Vr Vr
(25)
The relative supersaturation is given by Eq. (1). We assume that the saturation concentration sat is constant. The differential equation
whose solution is readily found to be: (29)
where 0 for 0, () = 1 otherwise,
(30)
is the Heaviside step function. According to Eq. (26), the steady-state mass balance in terms of supersaturation for constant volume, dVr /dt = 0, and negligible crystal volume, = 1, is ∞ j 1 f1,ss () V(h) d. (31) 0 = [in − ss ] − s
sat
j
0
Integration by parts yields: ∞ j f1,ss () V(h()) d 0
j
= [f1,ss ()V(h(, ss ))]∞ 0 −
∞ 0
V(h(, ss ))
j f d. j 1,ss
(32)
Using Eqs. (28a) and (28c) and the assumption that V( = 0) → 0, we rewrite the preceding equation as ∞ ∞ j 1 f1,ss () V(h(, ss )) d = − V(h(, ss )) f1,ss () d. (33) 0
j
0
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
691
Then we obtain for Eq. (31): 0 = [in − ss ] −
s sat
∞ 0
V(h(, ss ))f1,ss () d
(34)
and with Eq. (29) 0 = [in − ss ] −
s B sat 0
∞ 0
V(h(, ss )) exp(−/ ) d.
(35)
This equation in conjunction with a single crystal shape model allows us to determine the steady-state supersaturation ss . Note that the trajectory h(, ss ) depends on the steady-state supersaturation ss , which is to be solved iteratively. In what follows we present two case studies illustrating the application of the developed framework for a simple case of cuboid crystals followed by the complex case of paracetamol crystals. 4. Model simulations 4.1. Steady-state distributions for cuboid-shaped seed crystals
Fig. 5. steady-state concentration ss versus mean residence time .
For a first case study we consider a well-mixed crystallizer which has a mean residence time of a few seconds. We shall study the influence of the mean residence time on averaged geometrical properties. It is also shown that for this specific system we find a unique solution for the steady-state supersaturation that is multiple steadystates do not occur. The crystals under consideration posses a cuboid shape with face distances, measured from the crystal center, h1 , h2 and h3 , see Fig. 4. The crystals are symmetric about three planes, therefore the three distances are sufficient to describe the shape with six faces. Hence, the geometrical state vector is ⎛
⎞ h1 h = ⎝ h2 ⎠ . h3
(36)
The crystal volume is given by V(h) = 8h1 h2 h3 .
(37)
Assume that the size of a nucleus is zero (hi,0 = 0, i = 1 . . . 3), then the volume as a function of the crystal age is ˙ H ˙ ˙ 3 V(, ) = 8H 1 2 H3 .
(38)
Let the growth rates depend on the supersaturation only and be of the form of a power law ˙ = k gi . H i g,i
(39)
in = 2 s = 100 kg m−3 sat = 1 kg m−3
1 kb = 48 × 1016 # m−3 s−1 b=0 kg,1 = kg,2 = kg,3 = 10−6 m s−1 g1 = 1 g2 = 12 g3 = 3
The system parameters do not refer to a specific physical system.
According to Eq. (35) the steady-state concentration can be determined by solving the following equation ∞ b+g¯ ss = in − kb k¯ g ss 3 exp(−/ ) d, (41) 0
where = s / sat , k¯ g = 8 kg,i and g¯ = gi . The integral can be solved analytically: ∞ 3 exp(−/ ) d = 4 3!. (42) 0
Thus, ¯
b+g ss = in − 6 kb k¯ g 4 ss .
The nucleation follows the law: B 0 = k b b .
Table 1 Values for model parameters
(40)
(43)
Essentially, one has to solve the equation:
ss + ss − in = 0, = 6 kb k¯ g 4 > 0, ¯
= b + g, 0 < ss < in .
(44)
The derivative of the function g(ss ) = ss + ss − in is dg
−1 = ss + 1. dss
Fig. 4. Crystal comprising of six faces. The distances of the opposite faces to the origin of the coordinate system are assumed to be equal.
(45)
It follows that g(ss ) is strictly increasing for ss > 0 as > 0 and > 1. Therefore it has only one root (for ss > 0). Thus, Eq. (43) has only one solution in the relevant domain. We can compute the zeros of g(ss ) using an appropriate numerical method, e.g. the Newton–Raphson method. The steady-state supersaturation as a function of the mean
692
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
residence time for in = 2, = 1, = 4.5 is shown in Fig. 5. The parameters for growth and nucleation are given in Table 1. The residence time distribution f1,ss 1 = exp(−/ ) E() = ∞ f d 0 1,ss
In the next example we illustrate the application of the population balance framework and steady state analysis to simulate the shape distributions of the population of more complex paracetamol crystals.
(46)
˙ ss ( )) can be used in conjunction with the trajectory, h(, ) = H( to obtain average geometrical properties as a function of the mean residence time: ∞ ¯ ) = h( h(, )E() d. (47) 0
4.2. Steady-state distributions for paracetamol crystals We consider an MSMPR crystallizer in which monoclinic paracetamol is crystallized from aqueous solution. We predict the steadystate crystal shape distributions for different mean residence times and inlet concentrations. Monoclinic paracetamol has space group symmetry P21 /a with unit cell parameters a = 126.51 nm, b = 88.87 nm, c = 72.36 nm, = 114.848 nm, see Boerrigter et al. (2002). For the following case study ¯ {0 1 1} and {0 0 1} are selected. the most likely forms {1 1 0}, {2 0 1}, Further we assume full symmetry of the crystal shape. This leads to a four-dimensional geometrical state vector: ⎛ ⎞ h{1 1 0} ⎜h ⎟ ⎜ {2 0 1} ¯ ⎟ ⎟. (48) h=⎜ ⎜ ⎟ ⎝ h{0 1 1} ⎠
Fig. 6 shows the variation in the distances of three different faces at the steady state as a function of mean residence time. We use this information to reconstruct mean crystal geometries for different mean residence times. At = 0.36 we find that crystals have a needle-like habit, see Fig. 7(a). For = 1 and 2.5 one obtains cubes and a plate-like habit, respectively, as depicted in Fig. 7(b) and (c). Thus one can use the residence time as a control parameter to obtain one form preferably over the other. This is especially important for crystalline products which have considerably different dissolution behavior on different faces. Thus one can control the product dissolution behavior by varying the mean residence time of the crystallizer. Also it is necessary to control crystal shapes with respect to downstream processability. For instance the pressure drop in filter cakes containing mainly plate-like crystals is larger than in filter cakes containing more compact forms.
Ristic et al. (2001) have investigated the face-specific growth behavior of rather large crystals. The seed crystals with a volume of about 0.125 cm3 were grown to a crystal of about 2.5 cm3 . We have fitted their experimental data (see Fig. 8) over the range, 0 < < 0.2
Fig. 6. Mean face distances h¯ 1 (solid), h¯ 2 (dash-dot) and h¯ 3 (dotted) as a function of mean residence time .
¯ (dashed), {0 1 1} (dash-dot), and Fig. 8. Growth rates for forms {1 1 0} (solid), {2 0 1} {0 0 1} (dotted).
h{0 0 1}
Fig. 7. Average crystal geometries according to Fig. 6. = 0.36 s (needle-like habit), = 1 s (cube), = 2.5 s (plate-like habit) (from left to right).
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
693
x 10−5 ¯ (011)
( 011)
1
z
( 001) ¯ (110)
0
( 110) −1 ¯ (201)
−2 x 10−5
0 −1 x 10−5
0 y
1
x
2
Fig. 10. Paracetamol seed. Geometrical state is h0 = (1, 1, 1, 1) × 10−5 m. T
Fig. 9. Steady-state morphologies for = 0.2, = 0.15, = 0.1, = 0.05 (from left to right).
to obtain the following growth expressions: ⎧ 1.04 × 10−12 for < 0.1076, ⎪ ⎪ ⎨ −12 2 2.13 × 10 ˙ H {1 1 0} = ⎪ −2.14 × 10−11 ⎪ ⎩ −4.37 × 10−12 otherwise,
(49a)
−14 3 + 1.85 × 10−12 2 ˙ H ¯ = − 5.88 × 10 {2 0 1}
− 1.13 × 10−12 ,
(49b) Fig. 11. Steady-state supersaturation ss as a function of mean residence time for various inlet concentrations in = 0.2, 0.15, 0.1, 0.05 (top to bottom).
−14 3 + 7.46 × 10−13 2 ˙ H {0 1 1} = − 2.24 × 10
+ 5.94 × 10−12 , ⎧ −1.27 × 10−14 2 ⎪ ⎪ ⎪ +2.43 × 10−12 ⎪ ⎨ −14 3 ˙ H {0 0 1} = ⎪ −5.70 × 10 ⎪ ⎪ +1.87 × 10−12 2 ⎪ ⎩ −5.98 × 10−12
(49c) for < 0.0533, (49d) otherwise.
It is assumed here that these growth expressions are true even for small crystals, that is, the growth rates are independent of the crystal ˙ / jh = 0. Fig. 9 shows size (or in general the crystal geometry): jH i j the steady-state morphologies for crystal age → ∞, that is in the sense that the initial geometry can be ignored, see Eq. (7). It is further reported in Ristic et al. (2001) that for < 0.2 nucleation could be avoided in their growth cell. We assume that this is also true for the crystallizer under consideration. Let the seed crystals be added to the crystallizer at a rate of B0 = 1010 #/m3 /s. The seeds have the geometrical state ⎛ ⎞ 1 ⎜1⎟ ⎜ h0 = ⎝ ⎟ × 10−5 m. 1⎠ 1
(50)
The corresponding crystal shape is shown in Fig. 10. Using the conditional equation for the steady-state supersaturation (35) in connection with the trajectory of a single crystal h(, ss ), we can determine the steady-state supersaturation as a function of mean residence time as shown in Fig. 11. Here we observe the variation in the steady-state supersaturation for various inlet concentrations of in = 0.05, 0.1, 0.15, 0.2, which of course decreases for higher values for . At ss = 0.0533 and 0.1076 we find discontinuities, which reflect the discontinuities in the growth kinetics, Eq. (49), which clearly influence the mass consumption of a single crystal Gm , Eq. (16), and thus the mass consumption of the whole population, see Eq. (23). We analyze the crystallizer with mean residence time, =3×103 s. The steady-state supersaturation ss (in , ) for various inlet concentrations can be inferred from Fig. 11 as: ss (0.2, 3 × 103 s) = 0.1158, ss (0.15, 3 × 103 s) = 0.1107, ss (0.1, 3 × 103 s) = 0.0864, ss (0.05, 3 × 103 s) = 0.0457. We see that the increase of the inlet supersaturation from in = 0.05 to 0.1 leads to an increase in the steady-state supersaturation by ss = 0.0407. In comparison, the increase from in = 0.15 to ss = 0.2 yields an increase of the steady-state super-
694
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
Fig. 12. Volume distribution for in = 0.2, ss = 0.1158 and corresponding crystal shapes.
Fig. 13. Volume distribution for in = 0.05, ss = 0.0458 and corresponding crystal shapes.
saturation of only ss = 0.0051. The observation can be readily attributed to the morphology transitions and hence the consumption of supersaturation due to anisotropic growth. We now consider crystal distributions in a crystallizer with the residence time, = 3 × 103 s. From the trajectory h(, ), one can calculate the volume of a single crystal using Eq. (15). This allows the calculation of the volume distribution as follows: V(, )f1,ss () q3 (, ) = ∞ . 0 V(, )f1,ss () d
(51)
Fig. 12 shows the volume distribution q3 () for = 3 × 103 s with in = 0.2 plotted against crystal age. Since crystal age is not a measurable quantity the volume distribution against the sphere equivalent diameter is additionally shown in Fig. 14. The crystal shapes for = 5 × 103 , 104 , 1.5 × 104 s are also shown in Fig. 12. We observe that the shapes are similar as per Eq. (7). Fig. 13 depicts volume distributions for the same mean residence time but for a different inlet concentration (in = 0.05). It also shows crystal shapes at various ages, namely = 2.5 × 103 , 5 × 103 , 104 s. Clearly, we observe that the crystal shape changes remarkably with its age. This can be attributed to the fairly smaller growth rates at ss = 0.0457 which makes the initial geometry, Eq. (50), the deciding factor in the final shape. This can also be seen in Fig. 15. Here the mass growth rate Gm is plotted as a function of the crystal's surface area. If the steadystate shape for a given supersaturation is reached, that is shapes at different instants are similar in the sense of Eq. (7), then the mass growth rate is proportional to the surface area A, see dashed lines in Fig. 15. However, if the shape is considerably influenced by the initial geometry, which means that it is not at steady state, then the mass growth rate can depart significantly from linear behavior. For the growth of paracetamol at ss = 0.1158 (as discussed in Fig. 12) the mass growth is higher than for the steady-state geometry. This is because faster growing faces have a larger surface area compared to the steady-state geometry. As the crystal grows, it approaches the steady-state geometry and a linear growth law. Paracetamol growth at ss = 0.0457 shows particularly interesting behavior because the
Fig. 14. Volume distribution for in = 0.05, ss = 0.0458 versus sphere equivalent diameter.
mass growth rate Gm decreases with larger surface area in a first growth period, see Fig. 15. This comes from the fact that the form ¯ is relatively fast compared to surrounding forms (compare to { 2 0 1} ¯ faces, Figs. 8 and 10). Hence, a lot of material is attached to the {2 0 1} but the complementary effect of the high growth rate is that these faces grow out—therefore their area becomes continuously smaller. Because of that one observes that the growth rate increases with rising surface area as soon as the forms which will grow out have only marginal influence on the overall mass growth. Certainly, inlet concentration and mean residence time are not the only control variables we may play around in order to engineer the crystal morphology. The seed rate is one such control variable. In the earlier illustrations it was set to B0 = 1010 #/m3 /s. It is our
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
τ = 3⋅104
10−16 Mass Growth Rate Gm [kg s−1]
695
σss = 0.1158 τ=0
10−17
σss = 0.0457 τ=0
τ = 3⋅104 Faces 201 disappear
10−18 10−9
10−8 Surface Area A
[m2]
Fig. 15. Mass growth rate versus surface area for the cases shown in Figs. 12 (top, solid, crystal age = 0 . . . 3 × 104 ) and 13 (bottom, solid for crystal age = 0 . . . 3 × 104 , dotted for > 3 × 104 ). Mass growth rates for crystals possessing steady-state geometry in the sense of Eq. (7) are proportional to the surface area (dashed lines). The mass consumption rate of crystals with shapes being not in steady state exhibit a nonlinear dependency of the surface area.
Fig. 17. Volume distribution for in = 0.2, B0 = 109 , yielding ss = 0.1398 and corresponding crystal shapes for selected crystal ages .
Fig. 16. steady-state supersaturation ss as a function of mean residence time for different seed rates B0 = 1011 , 1010 , 109 , 108 #/m3 /s (left to right).
objective to illustrate its effect. Fig. 16 shows the steady-state supersaturation as a function of mean residence time for various seed rates in a range B0 = 109 . . . 1011 , where in = 0.2. One can recognize that the steady-state supersaturation obviously depletes faster as the seed rate is increased. This is obvious as the increase in crystal population increases the consumption of the supersaturation due to cumulative growth. The volume distribution for = 3 × 103 s, yielding a steady-state concentration of ss = 0.1398 is shown in Fig. 17. Also the shapes as a function of crystal age have been shown in the figure that indicate the significant effect of the crystal age on the shape evolution. 5. Conclusion The multidimensional population balance framework of this paper has been able to extract shape distributions in a steady state
continuous crystallizer for two different cases. One involves seed crystals in the feed with no nucleation in the crystallizer and the other with no crystals in the feed but nucleation providing crystals in the system. The effects of (i) the mean residence time, (ii) seeding rate of crystals, and (iii) inlet supersaturation have been investigated with consequences that are readily understood. All of these affect the shape and morphology distributions significantly. Although energy balance was not included in this analysis, its coupling with the mass and population balance to describe the complete dynamics of the crystallizer does not present a particularly challenging computational problem for the cases at hand. Other processes such as crystal breakage and aggregation that occur in real crystallizers perhaps deserve inclusion for a more realistic analysis. There also are other complicating sources such as the effect of impurities, surfactants, hydrodynamic forces, and so on that could affect crystal growth that have been omitted in this effort. In light of all of the foregoing factors, our current effort must be regarded as a modest beginning in spite of its focus on a more realistic crystal growth. Of course, in conclusion, it is fair to state that the quality of predictions from such models is contingent upon the quality of information on crystal facespecific growth rates whether they are obtained from experiments (e.g. Mullin and Whiting, 1980) or theory (e.g. Liu et al., 1995; Winn and Doherty, 1998; Boerrigter et al., 2002).
Notation B0 f1 (x, t) F Gm h hi ˙ H
nucleation rate, # m−3 s−1 number density, # [xi ]−1 m−3 s−1 flow rate, m3 s−1 mass transfer rate, kg s−1 geometrical state vector distance of ith face to the origin O, m velocity in h-space
696
˙ H i ˙ m n ni N Q r
S Tij vijk V Vc Vijk Vr x
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
velocity of ith face mass flow rate, kg s−1 number of crystal faces unit normal of ith crystal face matrix containing unit normals morphology matrix coordinates in 3-d physical space crystal shape (point set in 3-d physical space) tetrahedra matrix coordinates of Vijk crystal volume, m3 volume of continuous phase, m3 vertex at which faces i, j, k meet volume of entire system, m3 crystal state vector
Greek letters
(x) (x) i sat in ss ∇h
Dirac delta [x]−1 void fraction Heaviside step mean residence time, s number of vertices of ith face concentration of solute, kg m−3 saturation concentration of solute, kg m−3 supersaturation inlet supersaturation steady-state supersaturation crystal age, s domain nabla operator in shape space
References Boerrigter, S.X.M., Cuppen, H.M., Ristic, R.I., Sherwood, J.N., Bennema, P., Meekes, H., 2002. Explanation for the supersaturation-dependent morphology of monoclinic paracetamol. Crystal Growth and Design 3, 357–361. Joshi, M.S., Paul, B.K., 1974. Effect of supersaturation and fluid shear on habit and homogeneity of potassium dihydrogen phosphate crystals. Journal of Crystal Growth 22, 321–327. Liu, X.Y., Boek, E.S., Briels, W.J., Bennema, P., 1995. Prediction of crystal-growth morphology based on structural-analysis of the solid-fluid interface. Nature 374, 342–345. Mullin, J.W., Whiting, M.H.L., 1980. Succinic acid crystal-growth rates in aqueous solution. Industrial & Engineering Chemistry Fundamentals 19, 117–121. Petersen, K.B., Pederson, M.S., 2006. The matrix cookbook. Technical University of Denmark http://www2.imm.dtu.dk/pubdb/p.php?3274 . Prywer, J., 2005. Kinetic and geometric determination of the growth morphology of bulk crystals: recent developments. Progress in Crystal Growth and Characterization of Materials 50, 1–38. Ramkrishna, D., 2000. Population Balances. Academic Press, San Diego. Ristic, R.I., Finnie, S., Sheen, D.B., Sherwood, J.N., 2001. Macro- and micromorphology of monoclinic paracetamol grown from pure aqueous solution. Journal of Physical Chemistry 105, 9057–9066. Sangwal, K., 1998. Growth kinetics and surface morphology of crystals grown from solutions: recent observations and their interpretation. Progress in Crystal Growth and Characterization of Materials 36, 163–248. Winn, D., Doherty, M.F., 1998. A new technique for predicting the shape of solutiongrown organic crystals. A.I.Ch.E. Journal 44, 2501–2514. Yang, G., Kubota, N., Sha, Z., Louhi-Kultanen, M., Wang, J., 2006. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Crystal Growth and Design 6, 2799–2803. Zhang, Y., Doherty, M.F., 2004. Simultaneous prediction of crystal shape and size for solution crystallization. A.I.Ch.E. Journal 50, 2101–2112. Zhang, Y., Sizemore, J.P., Doherty, M.F., 2006. Shape evolution of 3-dimensional faceted crystals. A.I.Ch.E. Journal 52, 1906–1915.
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s
On the prediction of crystal shape distributions in a steady-state continuous crystallizer Christian Borchert a , Nandkishor Nere b , Doraiswami Ramkrishna b,∗ , Andreas Voigt c , Kai Sundmacher a,c a b c
Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, D-39106 Magdeburg, Germany Purdue University, School of Chemical Engineering, 480 Stadium Mall Dr., West Lafayette, IN 47907, USA Otto von Guericke University, Process Systems Engineering, PSF 410, D-39106 Magdeburg, Germany
A R T I C L E
I N F O
Article history: Received 26 November 2007 Received in revised form 27 March 2008 Accepted 15 May 2008 Available online 22 May 2008 Keywords: Crystallization Modeling Morphology Multidimensional population balance
A B S T R A C T
Zhang and Doherty [2004. Simultaneous prediction of crystal shape and size for solution crystallization. A.I.Ch.E. Journal 50, 2101–2112] have provided a one-dimensional analysis of crystallization based on the assumption that the relative face-specific growth rates of a (2-D) crystal are independent of supersaturation and hence invariant with time. Subsequent work by these authors [Zhang, Y., Sizemore, J.P., Doherty, M.F., 2006. Shape evolution of 3-dimensional faceted crystals. A.I.Ch.E. Journal 52, 1906–1915) consider shape evolution of single three-dimensional crystals with morphological changes. In this work, we present a multidimensional population balance approach accounting for dependence of the relative face-specific growth rates on supersaturation, a situation more commonly encountered. For example, Joshi and Paul [1974. Effect of supersaturation and fluid shear on habit and homogeneity of potassium dihydrogen phosphate crystals. Journal of Crystal Growth 22, 321–327] and Mullin and Whiting [1980. Succinic acid crystal-growth rates in aqueous solution. Industrial & Engineering Chemistry Fundamentals 19, 117–121] report face-specific growth rates with different dependence on the supersaturation. Thus it has been observed that there exists significantly different crystal shapes in a crystallizer [Yang, G., Kubota, N., Sha, Z., Louhi-Kultanen, M. Wang, J., 2006. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Crystal Growth and Design 6, 2799–2803]. Consequently, the population of crystals at any instant will have widely varying crystal shapes and sizes depending upon the initial crystal shape and size distribution. Computations are presented for the shape distributions of the crystal population emerging from a steady-state continuous crystallizer for two cases: (1) feed without crystals including nucleation for the formation of new crystals, and (2) feed with seed crystals of known shape, with suppressed nucleation. In the range of mean residence times investigated, the calculated crystal volume distributions for the first case show geometrically dissimilar shapes without morphological variations. However, in the second case, because the feed crystals of the chosen shape were susceptible to morphological changes, the volume distributions display this feature with shape and size distributions for each of a number of different morphologies. By varying operating conditions such as the flow rate, the inlet supersaturation, and the shapes of feed crystals, the proposed model can clearly be used to manipulate the crystal shape and size distributions and their morphologies. © 2008 Elsevier Ltd. All rights reserved.
1. Introduction We embark in this paper on an analysis of crystal populations in the habitat of a crystallizer, for the latter is the source of crystalline products in industrial practice, where the shape of faceted crystals is a major quality factor. In the chemical industries in general and in the speciality and fine chemical industries, in particular, the crystal shape and its distribution are of considerable interest because it can remarkably influence the product downstream processing and
∗
Corresponding author. Tel.: +1 765 494 4066; fax: +1 765 494 0805. E-mail address: [email protected] (D. Ramkrishna).
0009-2509/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.05.009
usability. For industrial applications it is important to predict, optimize and control shape distributions. Since the growth of single crystals depend on supersaturation in the environment, which must change due to the cumulative growths of all crystals around, the analysis must account for the reciprocal effects of crystal growth and supersaturation dynamics. Indeed population balances represent a natural implement for such an analysis. Since the potential complexity of a full multidimensional approach calls for some sense of parsimony in regard to the extent of detail, it behooves us to begin with an examination of scenarios for simplification. In this connection, the work of Zhang and Doherty (2004) is a landmark that represents in some way the most desirable simplification, viz., that crystal
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
growth may be described in a unidimensional space of one kind or another. We argue that, in describing crystal growth the authors have stretched the most out of a one-dimensional framework and that no improvement is possible without transgression into higher dimensions. More specifically, the authors assume that the relative face-specific growth rates are independent of the supersaturation and hence of time. This assumption traces the evolution of crystal growth to be along a rectilinear path in crystal shape space as determined by the vector, h ≡ (h1 , h2 , . . . , hn )T , representing distances of the various faces from the center of the crystal. This is a line crystal shape space, described by some parameter along it, which could either be the volume, mass or age of the crystal or a component of belonging to a secure face. Strict unidimensionality of a crystal population must confine all crystals to this line, which must necessarily pass through the origin of shape space. Clearly, all crystals lying on this line must be geometrically similar, characterized by a scalar property of one kind or another. Thus no new crystal formed by nucleation outside of this line will be geometrically similar to the set of crystals on the line. As Zhang et al. (2006) observe, faces of a crystal can appear or disappear during the course of growth signifying a morphological change. This change in morphology is indeed an interruption of its geometric similarity with its past although crystal growth will persist along the same direction with respect to hi 's representing unchanged faces. There are thus two kinds of outcasts from the geometrically similar crystals along the line of unidimensionality defined by constant relative face growths. The first are those that are born by nucleation outside of the line, and the second are those that undergo morphological changes. Strict unidimensionality cannot of course accommodate either kind. With the assumption of constant relative face growths, unidimensionality is still possible in an approximate sense. Fresh crystals formed by nucleation outside of the one-dimensional line can never be on it as it must proceed parallel to it from outside it. However, the shape of a growing crystal, as long as it does not change morphology, will become progressively insensitive to its initial perturbation from the line. Thus outcasts of the first kind, in course of time, may be regarded as geometrically similar in an approximate sense. Outcasts of the second kind also can similarly be inducted into the approximately similar set if morphological transitions terminate. Zhang and Doherty (2004) show that the geometric area factor of crystals becomes eventually constant as the population evolves through various morphological changes into an asymptotic shape. At this stage, a one-dimensional description of the crystal population becomes possible as the scalar property can be inverted to construct the full shape of each crystal. The discussion above must indicate that Zhang and Doherty (2004) have extracted the most out of a one-dimensional model. To seek improvement, one must admit multidimensionality (a fact recognized by Zhang and Doherty, 2004), making full use of crystal symmetry to minimize the increase in dimension. Before we set out to address our approach, however, we note one additional feature of the work of Zhang and Doherty (2004). While the assumption of constant relative face growth confines crystal growth to a line independently of supersaturation, the speed with which crystals move is regulated by supersaturation so that coupling between the two processes is accommodated. Our multidimensional framework, however, couples supersaturation dynamics intimately with growth of the crystal population. Crystal growth occurs along curved lines in shape space. Due to the increase in complexity, we confine our effort to relatively simple situations. For our treatment, we assume that the crystals are surrounded by a continuous phase in which physico-chemical processes occur leading to supersaturation. This work will ignore the path to supersaturation by assuming its direct availability.
687
Due to crystal growth, dissolved material is transferred from the continuous phase to the dispersed phase, more specifically: to the crystal faces. We assume that the crystal is built of only one species. The driving force for growth is the supersaturation. In the subsequent text we use the relative supersaturation defined by
=
− sat , sat
(1)
where is the mass concentration of the solute while sat its saturation concentration. The paper is organized as follows. In Section 2 we introduce a single crystal shape model along the lines of the work of Zhang et al. 2006). However, this work departs in one significant way from that of Zhang and Doherty (2004) in that the relative growth of different crystal faces, which was assumed constant in earlier work, is allowed to depend on supersaturation. The coupling of the single crystal with the continuous phase is discussed within the framework of population balance which accounts for the joint effects of the growth of all crystals in the population. Furthermore, we establish the mass balance in terms of supersaturation. Section 3 shows how one obtains a conditional equation for the steady-state supersaturation of a mixed suspension mixed product removal (MSMPR) crystallizer. In Section 4, the results of model simulations using the conditional equation for the steady-state supersaturation are discussed. The supersaturation description has been combined with (i) the classical age distribution of ideally mixed systems and (ii) with the single crystal shape evolution model from Section 2 to obtain shape distributions. In the first example, we consider an arbitrary cuboid crystal, which is a rather simple system while the second example shows how the developed tool can be applied to a more complex case (paracetamol) using the growth rate measurements by Ristic et al. (2001). Finally, we conclude by summarizing the outcomes of the present work in brief. 2. Model formulation Crystallizers contain a large number of crystals, called a crystal population. The major task in crystallization is to grow crystals. Therefore we focus on growth throughout the paper and neglect breakage and aggregation. The framework that allows the investigation of the dynamics of a crystal population is the population balance equation. However, in order to describe the behavior of a crystal population properly one should fully understand the dynamics of a single crystal. Thus we introduce at first the shape model for a single crystal followed by the discussion of some of the important properties. Finally, the population balance model for the suspension is introduced. 2.1. Shape model We consider crystals with a convex shape composed of n flat faces. A point O is chosen, for example the crystal center, to be the origin of a Cartesian coordinate system in which the crystal form is described. The natural starting point to decompose a convex crystal is to consider each face separately. The orientation of the ith face is given by its unit normal ni , which is obtained from transforming Miller indices to normal vectors in Cartesian coordinates. The distance of the ith face to the origin O is denoted by hi . The ith face of a crystal is a finite sector on the plane ni · r = hi ,
i = 1 . . . n,
(2)
where r is a point in 3-D physical space. Clearly, hi is a property of the crystal that serves to describe its geometry. Thus hi is interpreted as the ith entry of the geometrical state vector h = [h1 . . . hn ]T . The geometry of a crystal is
688
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
altered when h changes. Therefore we will consider the change of h due to face-specific crystal growth. Note that the orientation ni of a face does not change due to growth but only hi . Microscopically, a crystal grows by the attachment of growth units to its surface. We regard the growth units to be orders of magnitude smaller than the crystal so that the discrete statistical processes of crystal growth can be approximated by a continuous model. The growth of a crystal is usually anisotropic, i.e. the growth rates of the faces follow different laws, see e.g. Sangwal (1998). The dynamics of h(t) is described by the vector differential equation: dh ˙ = H(h, ), dt
h(0) = h0 ,
(3)
˙ = [H ˙ n ]T , and H ˙ is the perpendicular growth rate of the ith ˙ , ... , H H 1 i ˙ is also called the face velocity and is face, see Zhang et al. (2006). H i in general a function of the crystal state and supersaturation . The crystal shape S is defined by the planes containing crystal faces as
S = {r : ni · r < hi , i = 1 . . . n},
(4)
which is a convex polyhedron. The shape evolution is given by t ˙ S(t) = r : Nr < h(t), h(t) = h0 + )), ()) d , H(h( (5) 0
where ⎡ T⎤ n1 ⎢ . ⎥ ⎥ N=⎢ ⎣ .. ⎦
For many crystals we find, at least for some idealized situation a symmetrical crystal shape. Obviously it is then sufficient to consider only one of the symmetrical faces. Hence, we can reduce the dimension of the geometrical state vector h remarkably when we exploit crystal symmetry. There may be other properties of interest, for example the crystal age, color, density or composition. In this paper we will have only one additional property: crystal age . Hence, the full property state vector is x = [hT ]T .
(6) 2.2. Dynamics of real and virtual faces
nTn is the n by three matrix containing the unit normals of the faces. Shapes S1 and S2 are defined to be similar when hS = shS , 2 1 where s is a scalar. Usually, the shapes at different instants are not ˙ the influsimilar due to anisotropic effects in growth. For t?hi,0 / H i ence of the initial geometry h0 on h is negligible compared to the influence of the growth. Further, if > 0 is at steady state, denoted ˙ does not by ss , as, for example, in a continuous crystallizer and if H i ˙ depend on the geometry (jHi / jhj = 0), we find by replacing time t with age , a variable more appropriate for the dynamics of individual ˙ . That means s = crystals in a steady-state crystallizer that h() = H 2 / 1 is a scaling factor to relate crystal shapes at different instants: h(1 ) = sh(2 ),
Fig. 1. In case that the face with the thick line has a remarkably higher growth rate than the neighboring faces, its surface area becomes lower until it grows out (top) or its increases (bottom).
(7)
which means that the shapes are similar. If successive shapes are similar under certain growth conditions (for instance supersaturation) then we speak of the steady-state geometry. It is opportune at this stage to discuss the precise difference between this work and that of Zhang and Doherty (2004). In the discussion below, we will refer to the space to which the vector h belongs as the crystal shape space or simply shape space. Zhang and Doherty (2004) assume the relative growths of different faces to remain constant in time regardless of the supersaturation. Thus their view traces the evolution of shape along a rectilinear trajectory in shape space, independently of supersaturation dynamics in the continuous phase. Of course insofar as supersaturation does influence the growth of each face it dictates the speed with which shape evolution occurs without changing its direction. In contrast, in our work supersaturation not only determines the speed with which shape evolution occurs but also its direction. Thus in the steady state envisaged above although the present work will feature crystals that are geometrically similar, the shape and size will vary with residence time due to variation in the steadystate supersaturation. On the other hand, Zhang and Doherty (2004) will only perceive differences in the shape for a single crystal. Their work ignores the feedback of crystal growth on supersaturation.
Since crystal faces grow at different rates, it can happen that one face is much faster than its neighbors so that it finally grows out, see Fig. 1. Therefore one generally has only n < n faces present on the crystal surface, whereas there are n = n − n faces which are virtual (likely faces). In the following we assign to the distance and unit normals of virtual faces a double-prime to identify them as such, for example hi and ni . In case that we explicitly refer to a real face we provide the distance and unit normal with a single prime, hi and ni . We stipulate that a virtual face is always tangential to the crystal, see Fig. 2. Since no two faces have the same orientation, the polyhedron S has an extreme point Vijk in the direction n of the virtual face . The point Vijk is a vertex defined by the intersection point of at least three real faces i, j, k. The coordinates of the intersection are obtained by vijk = N−1 h , ijk ijk
where N−1 = [ni nj nk ]T , ijk
hijk = (hi hj hk )T .
(8)
Note that there exists an extreme line (edge) Vijk Vijl on the crystal when n ·(vijk −vijl )=0. Then it is sufficient for our purposes to iden tify only one point on Vijk Vijl , e.g. Vijk to be the extreme point of S with respect to n . The distance h is dictated by vijk because the vir tual face is tangential to S. Then h is the projection of vijk onto the direction n , that is h = n · vijk . By the virtue of Eq. (8) the distance of the virtual face is a linear function of the distances of real faces: h = q · hijk ,
q = (n )T N−1 . ijk
(9)
Constructing an n -dimensional real geometrical state vector h which contains only the distances of the real faces and a virtual geometrical state vector h (n -dimensional), we can write h = Q h ,
(10)
where the n by n matrix Q contains on its rows the entries of q on the corresponding position. The matrix Q changes its entries and
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
689
Fig. 2. Different morphologies of the same crystal exhibiting different faces. Even though the faces which cut the edges of the octahedron (left) are not real (green faces) we have to keep track of them due to the fact that they may appear (right). We call them virtual faces (blue) and they are always tangential to the crystal shape.
dimension depending on the geometrical configuration of the crystal. It is especially necessary to track the number of real and virtual faces n and n . That is, for qualitatively different geometrical configurations as shown in Fig. 2 we find different matrices Q . An increase in dimension in the real part of the shape space (h space) induced by the conversion of a virtual face into a real one must be viewed as a step change in the morphology as it will signal a discontinuous change in the number of vertices, edges or planes. This new morphology, however, is deemed to remain invariant until the next step change occurs. ˙ =dh /dt, have to be determined from Growth rates of real faces, H i i crystal growth models and verified by experiments. We assume here ˙ is a function of supersaturation without loss of generality that H i only. There exist several crystal growth theories in the literature, see for instance Sangwal (1998) and Prywer (2005). The displacement rate of a virtual face is—like its position—a function of growth rates of real faces. Differentiation of Eq. (10) with respect to time yields the velocity vector of virtual faces: ˙ = Q H ˙ (), H
˙ = dh where H dt
and
˙ () = dh . H dt
(11)
An algorithm for the calculation of shape evolution when faces can change their state from real to virtual and vice versa can be found in Zhang et al. (2006). For the present, we abstractly observe that step changes can occur in the evolution when transitions occur between virtual faces/edges and real faces/edges. Such transitions are characterized by equalities involving the growth vectors and the motion of vertices, which together define some hyper surface in the joint space of supersaturation and crystal shape. These hyper surfaces define the domain of operation for restricting the variation in morphology. 2.3. Mass balance for a single crystal
For practical calculations we decompose the polyhedron into pyramids Si whose base plane is the ith face and apex is the origin of the coordinate system O. Pyramids are again decomposed into tetrahedra Sij whose vertices are the origin and the vertices of the ith face, see Fig. 3. We number the coordinates of the vertices of the ith face from vi1 , . . . , vi . Then we can find i − 2 tetrahedra per face. i The vertex coordinates of the jth tetrahedron of the ith face is stored in the matrix j = 1 . . . (i − 2).
(13)
n i −2
V(S) =
V(Sij ).
(15)
i=1 j=1
In the following we represent the volume of a crystal, V = V(S) by V(h). Assuming that the density of a crystal, s , is constant, the rate of change of the crystal mass due to growth, Gm is given by n
jV d ˙ ˙ = ( ∇ V) · H, H V(h) = s s h dt jhi i
(16,17)
i=1
where ∇h =(j/ jh1 , . . . , j/ jhn ). This mass transfer term will later be used in the mass balance for the continuous phase. We can calculate the volume change as follows. The change of the volume requires the evaluation of the derivative of the tetrahedra V(Sij ), see Eq. (14): 1 d d V(Sij ) = | det(Tij )|. dt 6 dt
(18)
The derivative of the norm of a function f is d|f | = sign(f ) df . For the derivative of a determinant of a matrix A we find (see Petersen and Pederson, 2006): d det(A) = det(A) tr(A−1 d A), where tr is the operator delivering the trace of the matrix. The derivative of Tij is denoted by T˙ and contains the vertex velocity vectors v˙ which are ij
ij
obtained from differentiation of Eq. (8):
The volume of the tetrahedron Sij is readily found to be V(Sij ) = 16 | det(Tij )|.
Summation over all tetrahedra of a real face and summation over all n real faces yields the volume of the crystal
Gm = s
The volume of a crystal or rather the polyhedron S is given by V(S) = dV. (12) S
Tij = [vi1 vi(j+1) vi(j+2) ],
Fig. 3. Decomposition of the volume under the ith face into tetrahedra.
(14)
T˙ ij = [v˙ i1 v˙ ij+1 v˙ ij+2 ]T ,
j = 1, . . . , i − 2.
(19)
690
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
Then we get 1 d V(Sij ) = sign(V(Sij )) det(Tij ) tr((Tij )−1 T˙ ij ) dt 6 1 = sign(V(Sij )) tr((Tij )adj T˙ ij ), 6
(20a) (20b)
where adj refers to the adjugate or adjoint (transpose of the cofactor).
3. Steady-state analysis
2.4. Multidimensional population balance equation So far we have discussed the single crystal model. In the following we account for a whole population of crystals. The crystal state vector x = [hT , ]T is a point in the continuous property state space R(n+1) . Here we consider the domain of state space, to be of our interest. Particles entering originate at one point xnuc ∈ j with the rate of nucleation B0 . It is assumed that all nuclei have the same state xnuc . Following Ramkrishna (2000), the population balance equation for a population of crystals undergoing convective growth in a wellmixed continuous crystallizer is
j˜ ˙ f˜ (x, t) = − 1 f˜ (x, t), f (x, t) + ∇x · X 1 jt 1 1
(21a)
I.C.: f˜1 (x, 0) = f˜1,0 ,
(21b)
˙ · nx )f˜ (x, t) = B (x − xnuc ), B.C.: (−X 1 0
x ∈ j,
R.C.: f˜1 (x → ∞, t) = 0,
for the supersaturation using the preceding definitions in connection with Eqs. (16), (23) and (24) in Eq. (22) becomes
1 d d = (in + 1)Fin − ( + 1) Fout − ( + 1) ( Vr ) dt Vr dt s ˙ · (∇ V) dV . − (f˜ H) (26) h h sat 1
We assumed in Eqs. (21) that all crystals nucleated in the crystallizer have the same state xnuc ∈ j. At steady-state operation, all crystals see the same environmental conditions, namely the (steadystate) supersaturation ss . Therefore all crystals undergo the same shape evolution according to Eq. (5). That is, all crystals evolve along the same trajectory h(, ss ) through shape space. The time which a crystal spends within the crystallizer is the age . Suppose we know the steady-state supersaturation ss , then we can calculate the trajectory h(, ss ) for an individual crystal using Eq. (3). That means we know the shape as a function of crystal age, S(, ss ), and thus the volume V(, ss ) for all crystal ages. We consider the population balance equation only for the crystal age since shape is then recoverable purely from . Integration of Eq. (21) over h yields:
(21c)
j j 1 f (, t) + f (, t) = − f1 (, t), jt 1 j 1
(27a)
(21d)
I.C.: f1 (, 0) = f1,0 ,
(27b)
B.C.: f1 (0, t) = B0 (),
(27c)
R.C.: f1 ( → ∞, t) = 0,
(27d)
where f˜1 is the crystal number density: [f˜1 ] = #/ i [xi ]/m3 . The crystallizer volume is Vr and is the mean residence time. Vector nx is an outer unit normal on the boundary j. It is obvious that for the case of continuous seeding the equation will require an additional source term on the right-hand side.
where f1 (, t) = f˜1 (h, , t) dVh . The steady-state age distribution h is given by
j 1 f () = − f1,ss (), j 1,ss
(28a)
The continuous phase is in general characterized by a state vector y but we restrict our analysis to mass m of the crystal solute. Let ˙ in be the mass flow rate of the saturated solution into and m ˙ out m ˙ disp be the mass be the mass flow rate out of the crystallizer, and m consumption of the dispersed phase. The mass balance follows:
B.C. : f1,ss (0) = B0 (),
(28b)
R.C. : f1,ss ( → ∞) = 0,
(28c)
dm ˙ in − m ˙ out − m ˙ disp . =m dt
f1,ss () = B0 () exp(−/ ),
2.5. Mass balance for crystal population
(22)
The mass transfer of crystal solute from the continuous phase to the crystalline phase is given by ∞ ˙ disp = Vr Gm (h)f˜1 (h, , t) dVh d, (23) m 0 where dVh is an infinitesimal volume measure in h-space and Vr is the total volume of the slurry. Vr is the sum of the volume of the continuous phase and the volume of the dispersed phase: Vr = Vc + Vdisp . The mass transfer rate Gm is given by Eq. (16). Furthermore we write m = Vr ,
˙ in = in Fin , m
˙ out = Fout , m
(24)
where is the mass concentration of dissolved material, Fin and Fout is the volumetric flow rate into and out of the crystallizer and is the volume fraction of the liquid phase defined as
=
Vdisp Vc =1− . Vr Vr
(25)
The relative supersaturation is given by Eq. (1). We assume that the saturation concentration sat is constant. The differential equation
whose solution is readily found to be: (29)
where 0 for 0, () = 1 otherwise,
(30)
is the Heaviside step function. According to Eq. (26), the steady-state mass balance in terms of supersaturation for constant volume, dVr /dt = 0, and negligible crystal volume, = 1, is ∞ j 1 f1,ss () V(h) d. (31) 0 = [in − ss ] − s
sat
j
0
Integration by parts yields: ∞ j f1,ss () V(h()) d 0
j
= [f1,ss ()V(h(, ss ))]∞ 0 −
∞ 0
V(h(, ss ))
j f d. j 1,ss
(32)
Using Eqs. (28a) and (28c) and the assumption that V( = 0) → 0, we rewrite the preceding equation as ∞ ∞ j 1 f1,ss () V(h(, ss )) d = − V(h(, ss )) f1,ss () d. (33) 0
j
0
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
691
Then we obtain for Eq. (31): 0 = [in − ss ] −
s sat
∞ 0
V(h(, ss ))f1,ss () d
(34)
and with Eq. (29) 0 = [in − ss ] −
s B sat 0
∞ 0
V(h(, ss )) exp(−/ ) d.
(35)
This equation in conjunction with a single crystal shape model allows us to determine the steady-state supersaturation ss . Note that the trajectory h(, ss ) depends on the steady-state supersaturation ss , which is to be solved iteratively. In what follows we present two case studies illustrating the application of the developed framework for a simple case of cuboid crystals followed by the complex case of paracetamol crystals. 4. Model simulations 4.1. Steady-state distributions for cuboid-shaped seed crystals
Fig. 5. steady-state concentration ss versus mean residence time .
For a first case study we consider a well-mixed crystallizer which has a mean residence time of a few seconds. We shall study the influence of the mean residence time on averaged geometrical properties. It is also shown that for this specific system we find a unique solution for the steady-state supersaturation that is multiple steadystates do not occur. The crystals under consideration posses a cuboid shape with face distances, measured from the crystal center, h1 , h2 and h3 , see Fig. 4. The crystals are symmetric about three planes, therefore the three distances are sufficient to describe the shape with six faces. Hence, the geometrical state vector is ⎛
⎞ h1 h = ⎝ h2 ⎠ . h3
(36)
The crystal volume is given by V(h) = 8h1 h2 h3 .
(37)
Assume that the size of a nucleus is zero (hi,0 = 0, i = 1 . . . 3), then the volume as a function of the crystal age is ˙ H ˙ ˙ 3 V(, ) = 8H 1 2 H3 .
(38)
Let the growth rates depend on the supersaturation only and be of the form of a power law ˙ = k gi . H i g,i
(39)
in = 2 s = 100 kg m−3 sat = 1 kg m−3
1 kb = 48 × 1016 # m−3 s−1 b=0 kg,1 = kg,2 = kg,3 = 10−6 m s−1 g1 = 1 g2 = 12 g3 = 3
The system parameters do not refer to a specific physical system.
According to Eq. (35) the steady-state concentration can be determined by solving the following equation ∞ b+g¯ ss = in − kb k¯ g ss 3 exp(−/ ) d, (41) 0
where = s / sat , k¯ g = 8 kg,i and g¯ = gi . The integral can be solved analytically: ∞ 3 exp(−/ ) d = 4 3!. (42) 0
Thus, ¯
b+g ss = in − 6 kb k¯ g 4 ss .
The nucleation follows the law: B 0 = k b b .
Table 1 Values for model parameters
(40)
(43)
Essentially, one has to solve the equation:
ss + ss − in = 0, = 6 kb k¯ g 4 > 0, ¯
= b + g, 0 < ss < in .
(44)
The derivative of the function g(ss ) = ss + ss − in is dg
−1 = ss + 1. dss
Fig. 4. Crystal comprising of six faces. The distances of the opposite faces to the origin of the coordinate system are assumed to be equal.
(45)
It follows that g(ss ) is strictly increasing for ss > 0 as > 0 and > 1. Therefore it has only one root (for ss > 0). Thus, Eq. (43) has only one solution in the relevant domain. We can compute the zeros of g(ss ) using an appropriate numerical method, e.g. the Newton–Raphson method. The steady-state supersaturation as a function of the mean
692
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
residence time for in = 2, = 1, = 4.5 is shown in Fig. 5. The parameters for growth and nucleation are given in Table 1. The residence time distribution f1,ss 1 = exp(−/ ) E() = ∞ f d 0 1,ss
In the next example we illustrate the application of the population balance framework and steady state analysis to simulate the shape distributions of the population of more complex paracetamol crystals.
(46)
˙ ss ( )) can be used in conjunction with the trajectory, h(, ) = H( to obtain average geometrical properties as a function of the mean residence time: ∞ ¯ ) = h( h(, )E() d. (47) 0
4.2. Steady-state distributions for paracetamol crystals We consider an MSMPR crystallizer in which monoclinic paracetamol is crystallized from aqueous solution. We predict the steadystate crystal shape distributions for different mean residence times and inlet concentrations. Monoclinic paracetamol has space group symmetry P21 /a with unit cell parameters a = 126.51 nm, b = 88.87 nm, c = 72.36 nm, = 114.848 nm, see Boerrigter et al. (2002). For the following case study ¯ {0 1 1} and {0 0 1} are selected. the most likely forms {1 1 0}, {2 0 1}, Further we assume full symmetry of the crystal shape. This leads to a four-dimensional geometrical state vector: ⎛ ⎞ h{1 1 0} ⎜h ⎟ ⎜ {2 0 1} ¯ ⎟ ⎟. (48) h=⎜ ⎜ ⎟ ⎝ h{0 1 1} ⎠
Fig. 6 shows the variation in the distances of three different faces at the steady state as a function of mean residence time. We use this information to reconstruct mean crystal geometries for different mean residence times. At = 0.36 we find that crystals have a needle-like habit, see Fig. 7(a). For = 1 and 2.5 one obtains cubes and a plate-like habit, respectively, as depicted in Fig. 7(b) and (c). Thus one can use the residence time as a control parameter to obtain one form preferably over the other. This is especially important for crystalline products which have considerably different dissolution behavior on different faces. Thus one can control the product dissolution behavior by varying the mean residence time of the crystallizer. Also it is necessary to control crystal shapes with respect to downstream processability. For instance the pressure drop in filter cakes containing mainly plate-like crystals is larger than in filter cakes containing more compact forms.
Ristic et al. (2001) have investigated the face-specific growth behavior of rather large crystals. The seed crystals with a volume of about 0.125 cm3 were grown to a crystal of about 2.5 cm3 . We have fitted their experimental data (see Fig. 8) over the range, 0 < < 0.2
Fig. 6. Mean face distances h¯ 1 (solid), h¯ 2 (dash-dot) and h¯ 3 (dotted) as a function of mean residence time .
¯ (dashed), {0 1 1} (dash-dot), and Fig. 8. Growth rates for forms {1 1 0} (solid), {2 0 1} {0 0 1} (dotted).
h{0 0 1}
Fig. 7. Average crystal geometries according to Fig. 6. = 0.36 s (needle-like habit), = 1 s (cube), = 2.5 s (plate-like habit) (from left to right).
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
693
x 10−5 ¯ (011)
( 011)
1
z
( 001) ¯ (110)
0
( 110) −1 ¯ (201)
−2 x 10−5
0 −1 x 10−5
0 y
1
x
2
Fig. 10. Paracetamol seed. Geometrical state is h0 = (1, 1, 1, 1) × 10−5 m. T
Fig. 9. Steady-state morphologies for = 0.2, = 0.15, = 0.1, = 0.05 (from left to right).
to obtain the following growth expressions: ⎧ 1.04 × 10−12 for < 0.1076, ⎪ ⎪ ⎨ −12 2 2.13 × 10 ˙ H {1 1 0} = ⎪ −2.14 × 10−11 ⎪ ⎩ −4.37 × 10−12 otherwise,
(49a)
−14 3 + 1.85 × 10−12 2 ˙ H ¯ = − 5.88 × 10 {2 0 1}
− 1.13 × 10−12 ,
(49b) Fig. 11. Steady-state supersaturation ss as a function of mean residence time for various inlet concentrations in = 0.2, 0.15, 0.1, 0.05 (top to bottom).
−14 3 + 7.46 × 10−13 2 ˙ H {0 1 1} = − 2.24 × 10
+ 5.94 × 10−12 , ⎧ −1.27 × 10−14 2 ⎪ ⎪ ⎪ +2.43 × 10−12 ⎪ ⎨ −14 3 ˙ H {0 0 1} = ⎪ −5.70 × 10 ⎪ ⎪ +1.87 × 10−12 2 ⎪ ⎩ −5.98 × 10−12
(49c) for < 0.0533, (49d) otherwise.
It is assumed here that these growth expressions are true even for small crystals, that is, the growth rates are independent of the crystal ˙ / jh = 0. Fig. 9 shows size (or in general the crystal geometry): jH i j the steady-state morphologies for crystal age → ∞, that is in the sense that the initial geometry can be ignored, see Eq. (7). It is further reported in Ristic et al. (2001) that for < 0.2 nucleation could be avoided in their growth cell. We assume that this is also true for the crystallizer under consideration. Let the seed crystals be added to the crystallizer at a rate of B0 = 1010 #/m3 /s. The seeds have the geometrical state ⎛ ⎞ 1 ⎜1⎟ ⎜ h0 = ⎝ ⎟ × 10−5 m. 1⎠ 1
(50)
The corresponding crystal shape is shown in Fig. 10. Using the conditional equation for the steady-state supersaturation (35) in connection with the trajectory of a single crystal h(, ss ), we can determine the steady-state supersaturation as a function of mean residence time as shown in Fig. 11. Here we observe the variation in the steady-state supersaturation for various inlet concentrations of in = 0.05, 0.1, 0.15, 0.2, which of course decreases for higher values for . At ss = 0.0533 and 0.1076 we find discontinuities, which reflect the discontinuities in the growth kinetics, Eq. (49), which clearly influence the mass consumption of a single crystal Gm , Eq. (16), and thus the mass consumption of the whole population, see Eq. (23). We analyze the crystallizer with mean residence time, =3×103 s. The steady-state supersaturation ss (in , ) for various inlet concentrations can be inferred from Fig. 11 as: ss (0.2, 3 × 103 s) = 0.1158, ss (0.15, 3 × 103 s) = 0.1107, ss (0.1, 3 × 103 s) = 0.0864, ss (0.05, 3 × 103 s) = 0.0457. We see that the increase of the inlet supersaturation from in = 0.05 to 0.1 leads to an increase in the steady-state supersaturation by ss = 0.0407. In comparison, the increase from in = 0.15 to ss = 0.2 yields an increase of the steady-state super-
694
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
Fig. 12. Volume distribution for in = 0.2, ss = 0.1158 and corresponding crystal shapes.
Fig. 13. Volume distribution for in = 0.05, ss = 0.0458 and corresponding crystal shapes.
saturation of only ss = 0.0051. The observation can be readily attributed to the morphology transitions and hence the consumption of supersaturation due to anisotropic growth. We now consider crystal distributions in a crystallizer with the residence time, = 3 × 103 s. From the trajectory h(, ), one can calculate the volume of a single crystal using Eq. (15). This allows the calculation of the volume distribution as follows: V(, )f1,ss () q3 (, ) = ∞ . 0 V(, )f1,ss () d
(51)
Fig. 12 shows the volume distribution q3 () for = 3 × 103 s with in = 0.2 plotted against crystal age. Since crystal age is not a measurable quantity the volume distribution against the sphere equivalent diameter is additionally shown in Fig. 14. The crystal shapes for = 5 × 103 , 104 , 1.5 × 104 s are also shown in Fig. 12. We observe that the shapes are similar as per Eq. (7). Fig. 13 depicts volume distributions for the same mean residence time but for a different inlet concentration (in = 0.05). It also shows crystal shapes at various ages, namely = 2.5 × 103 , 5 × 103 , 104 s. Clearly, we observe that the crystal shape changes remarkably with its age. This can be attributed to the fairly smaller growth rates at ss = 0.0457 which makes the initial geometry, Eq. (50), the deciding factor in the final shape. This can also be seen in Fig. 15. Here the mass growth rate Gm is plotted as a function of the crystal's surface area. If the steadystate shape for a given supersaturation is reached, that is shapes at different instants are similar in the sense of Eq. (7), then the mass growth rate is proportional to the surface area A, see dashed lines in Fig. 15. However, if the shape is considerably influenced by the initial geometry, which means that it is not at steady state, then the mass growth rate can depart significantly from linear behavior. For the growth of paracetamol at ss = 0.1158 (as discussed in Fig. 12) the mass growth is higher than for the steady-state geometry. This is because faster growing faces have a larger surface area compared to the steady-state geometry. As the crystal grows, it approaches the steady-state geometry and a linear growth law. Paracetamol growth at ss = 0.0457 shows particularly interesting behavior because the
Fig. 14. Volume distribution for in = 0.05, ss = 0.0458 versus sphere equivalent diameter.
mass growth rate Gm decreases with larger surface area in a first growth period, see Fig. 15. This comes from the fact that the form ¯ is relatively fast compared to surrounding forms (compare to { 2 0 1} ¯ faces, Figs. 8 and 10). Hence, a lot of material is attached to the {2 0 1} but the complementary effect of the high growth rate is that these faces grow out—therefore their area becomes continuously smaller. Because of that one observes that the growth rate increases with rising surface area as soon as the forms which will grow out have only marginal influence on the overall mass growth. Certainly, inlet concentration and mean residence time are not the only control variables we may play around in order to engineer the crystal morphology. The seed rate is one such control variable. In the earlier illustrations it was set to B0 = 1010 #/m3 /s. It is our
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
τ = 3⋅104
10−16 Mass Growth Rate Gm [kg s−1]
695
σss = 0.1158 τ=0
10−17
σss = 0.0457 τ=0
τ = 3⋅104 Faces 201 disappear
10−18 10−9
10−8 Surface Area A
[m2]
Fig. 15. Mass growth rate versus surface area for the cases shown in Figs. 12 (top, solid, crystal age = 0 . . . 3 × 104 ) and 13 (bottom, solid for crystal age = 0 . . . 3 × 104 , dotted for > 3 × 104 ). Mass growth rates for crystals possessing steady-state geometry in the sense of Eq. (7) are proportional to the surface area (dashed lines). The mass consumption rate of crystals with shapes being not in steady state exhibit a nonlinear dependency of the surface area.
Fig. 17. Volume distribution for in = 0.2, B0 = 109 , yielding ss = 0.1398 and corresponding crystal shapes for selected crystal ages .
Fig. 16. steady-state supersaturation ss as a function of mean residence time for different seed rates B0 = 1011 , 1010 , 109 , 108 #/m3 /s (left to right).
objective to illustrate its effect. Fig. 16 shows the steady-state supersaturation as a function of mean residence time for various seed rates in a range B0 = 109 . . . 1011 , where in = 0.2. One can recognize that the steady-state supersaturation obviously depletes faster as the seed rate is increased. This is obvious as the increase in crystal population increases the consumption of the supersaturation due to cumulative growth. The volume distribution for = 3 × 103 s, yielding a steady-state concentration of ss = 0.1398 is shown in Fig. 17. Also the shapes as a function of crystal age have been shown in the figure that indicate the significant effect of the crystal age on the shape evolution. 5. Conclusion The multidimensional population balance framework of this paper has been able to extract shape distributions in a steady state
continuous crystallizer for two different cases. One involves seed crystals in the feed with no nucleation in the crystallizer and the other with no crystals in the feed but nucleation providing crystals in the system. The effects of (i) the mean residence time, (ii) seeding rate of crystals, and (iii) inlet supersaturation have been investigated with consequences that are readily understood. All of these affect the shape and morphology distributions significantly. Although energy balance was not included in this analysis, its coupling with the mass and population balance to describe the complete dynamics of the crystallizer does not present a particularly challenging computational problem for the cases at hand. Other processes such as crystal breakage and aggregation that occur in real crystallizers perhaps deserve inclusion for a more realistic analysis. There also are other complicating sources such as the effect of impurities, surfactants, hydrodynamic forces, and so on that could affect crystal growth that have been omitted in this effort. In light of all of the foregoing factors, our current effort must be regarded as a modest beginning in spite of its focus on a more realistic crystal growth. Of course, in conclusion, it is fair to state that the quality of predictions from such models is contingent upon the quality of information on crystal facespecific growth rates whether they are obtained from experiments (e.g. Mullin and Whiting, 1980) or theory (e.g. Liu et al., 1995; Winn and Doherty, 1998; Boerrigter et al., 2002).
Notation B0 f1 (x, t) F Gm h hi ˙ H
nucleation rate, # m−3 s−1 number density, # [xi ]−1 m−3 s−1 flow rate, m3 s−1 mass transfer rate, kg s−1 geometrical state vector distance of ith face to the origin O, m velocity in h-space
696
˙ H i ˙ m n ni N Q r
S Tij vijk V Vc Vijk Vr x
C. Borchert et al. / Chemical Engineering Science 64 (2009) 686 -- 696
velocity of ith face mass flow rate, kg s−1 number of crystal faces unit normal of ith crystal face matrix containing unit normals morphology matrix coordinates in 3-d physical space crystal shape (point set in 3-d physical space) tetrahedra matrix coordinates of Vijk crystal volume, m3 volume of continuous phase, m3 vertex at which faces i, j, k meet volume of entire system, m3 crystal state vector
Greek letters
(x) (x) i sat in ss ∇h
Dirac delta [x]−1 void fraction Heaviside step mean residence time, s number of vertices of ith face concentration of solute, kg m−3 saturation concentration of solute, kg m−3 supersaturation inlet supersaturation steady-state supersaturation crystal age, s domain nabla operator in shape space
References Boerrigter, S.X.M., Cuppen, H.M., Ristic, R.I., Sherwood, J.N., Bennema, P., Meekes, H., 2002. Explanation for the supersaturation-dependent morphology of monoclinic paracetamol. Crystal Growth and Design 3, 357–361. Joshi, M.S., Paul, B.K., 1974. Effect of supersaturation and fluid shear on habit and homogeneity of potassium dihydrogen phosphate crystals. Journal of Crystal Growth 22, 321–327. Liu, X.Y., Boek, E.S., Briels, W.J., Bennema, P., 1995. Prediction of crystal-growth morphology based on structural-analysis of the solid-fluid interface. Nature 374, 342–345. Mullin, J.W., Whiting, M.H.L., 1980. Succinic acid crystal-growth rates in aqueous solution. Industrial & Engineering Chemistry Fundamentals 19, 117–121. Petersen, K.B., Pederson, M.S., 2006. The matrix cookbook. Technical University of Denmark http://www2.imm.dtu.dk/pubdb/p.php?3274 . Prywer, J., 2005. Kinetic and geometric determination of the growth morphology of bulk crystals: recent developments. Progress in Crystal Growth and Characterization of Materials 50, 1–38. Ramkrishna, D., 2000. Population Balances. Academic Press, San Diego. Ristic, R.I., Finnie, S., Sheen, D.B., Sherwood, J.N., 2001. Macro- and micromorphology of monoclinic paracetamol grown from pure aqueous solution. Journal of Physical Chemistry 105, 9057–9066. Sangwal, K., 1998. Growth kinetics and surface morphology of crystals grown from solutions: recent observations and their interpretation. Progress in Crystal Growth and Characterization of Materials 36, 163–248. Winn, D., Doherty, M.F., 1998. A new technique for predicting the shape of solutiongrown organic crystals. A.I.Ch.E. Journal 44, 2501–2514. Yang, G., Kubota, N., Sha, Z., Louhi-Kultanen, M., Wang, J., 2006. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Crystal Growth and Design 6, 2799–2803. Zhang, Y., Doherty, M.F., 2004. Simultaneous prediction of crystal shape and size for solution crystallization. A.I.Ch.E. Journal 50, 2101–2112. Zhang, Y., Sizemore, J.P., Doherty, M.F., 2006. Shape evolution of 3-dimensional faceted crystals. A.I.Ch.E. Journal 52, 1906–1915.
