Ind. Eng. Chem. Res. 2007, 46, 3041-3047
3041
Transformation of the Chord-Length Distributions to Size Distributions for Nonspherical Particles with Orientation Bias† Nandkishor K. Nere,‡ Doraiswami Ramkrishna,*,‡ Bruce E. Parker,§ Willis V. Bell III,| and Pankaj Mohan§ School of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47907; Global Manufacturing SerVicessProcess Engineering Center, Eli Lilly and Company, Indianapolis, Indiana 46285; and Process Engineering Center, Eli Lilly and Company, Tippecanoe Laboratories, Lafayette, Indiana 47909
The applicability of the focused-beam reflectance measurement (FBRM) device for particle-size measurements depends on an appropriate transformation of the measured chord lengths to particle sizes. This problem has received considerable attention in the recent past. However, most of this effort has been directed toward either spherical particles or nonspherical particles of uniformly distributed orientations. In actual practice, the particle orientation must be biased by the flow of the fluid in which the particles are suspended. Such orientation bias becomes particularly serious for highly nonspherical particles such as those of needle shape that occur frequently in the pharmaceutical industry. This paper addresses the significant problem of transforming chord lengths of such particles with orientation bias to their sizes. Our approach is predicated on the solution of an inverse problem for extracting the orientation bias by measuring the chord-length distributions of particles whose size distributions are known. As the flow in most process equipment will be a function of position, the orientation bias will depend on the site of measurement. Thus, such a calibration approach will be specific to particle shape and measurement site because the flow conditions will vary with location. The proposed methodology is the first of its kind, which we expect to be a promising approach to fully exploit the simplicity of the FBRM device. 1. Introduction It is a pleasure to contribute to a festschrift in felicitation of Professor Man Mohan Sharma, a chemical engineer extraordinaire, who flourished as a researcher and educator, remarkably percipient on diverse issues of the chemical process industry. As his basic contributions have steadfastly focused on relevance to chemical process technology in one way or another, it was our view that this article must address a practical issue regardless of its mathematical base. Consequently, we consider here the application of a device that has had promise for online measurement of particle-size distributions in various dispersed phase processes. Determination of particle-size distribution (PSD) is commonly exercised in various chemical processes and unit operations, namely, polymerization, precipitation, crystallization, filtration, milling, and conveying, in which the quality of the production depends on the measurement and control of the PSD. Furthermore, PSD also plays a vital role in governing the reactivity in general and bioactivity of pharmaceutical molecules in particular. The measurement of chord lengths as an index of the particle size has become increasingly popular in recent times. Interpretation of such a chord-length distribution (CLD) is of paramount importance in order to extract the actual PSD, which can be significantly different from the measured CLD. Incorrect transformation of the measured data in terms of size distributions * To whom correspondence may be addressed. Phone: 1-765-4944066. E-mail:
[email protected]. Fax: 1-765-494-0805. † A part of this work was first presented at the 2005 AIChE annual meeting at Cincinnati, Ohio. ‡ Purdue University. § Global Manufacturing ServicessProcess Engineering Center, Eli Lilly and Company. | Process Engineering Center, Eli Lilly and Company, Tippecanoe Laboratories.
may adversely affect derived kinetic models and process monitoring, leading to flawed process design and control. It is the correct transformation of the CLD to PSD that is addressed in the present paper. In formulating the inverse problem of the determination of PSD given the CLD, we resolve the dependence of particle size on chord length into two sets of variables. The first refers to particle orientation as described by Eulerian angles between a fixed Cartesian frame of reference relative to the measurement site and another Cartesian frame attached to the particle; the second refers to geometric variables necessary to uniquely identify the chord length from particle size for a given orientation. The orientation bias is described by a distribution function for the Eulerian angles and is to be determined by solving the inverse problem to be defined presently. The geometric coordinates are assumed to be uniformly distributed, an assumption consistent with the use of a device that essentially measures the size distribution averaged over the measurement volume without recognition of gradients within the volume. The orientation bias so extracted is a more reasonable input to predicting particle size from chord-length measurements as long as the number densities are not so large as to affect the orientation bias by modifying the fluid flow. 2. FBRM Device and Past Work The focused-beam reflectance measurement (FBRM) technique has attracted considerable attention in the particulate process industry because of its simplicity and capability for online measurement. The schematic of the FBRM probe and the mechanism of the measurement of the chord are reproduced from ref 1 in Figure 1 for the ready reference. It consists of a rotating laser beam, which when passed through the suspension of particles receives reflectance for the duration of obstruction by the particle. The time is then converted using the linear laser
10.1021/ie0609463 CCC: $37.00 © 2007 American Chemical Society Published on Web 11/30/2006
3042
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
3. Mathematical Development
Figure 1. Schematic of FBRM probe (Worlitschek1).
velocity to get the length, which essentially is detected as a chord of a particle. Obviously, the chord is not a direct measure of the actual particle size, and must be transformed appropriately to particle size. The CLD measured by FBRM is of course a function of a measurement location. This is clearly due to the localized nature of the hydrodynamics prevailing at the measurement location. The location will decide CLD depending on the orientation of the particles that are presented to the probe. The transformation of FBRM chord-length measurements to particle size has of course received considerable attention in the past. However, most of these efforts have been directed toward either spherical particles or nonspherical particles of uniformly distributed orientations.1-3 The literature approach can be concisely described by the following generalized model for the CLD measured by FBRM.
fL(l) )
∫0∞ fL|X(l|x) fX(x) dx
(1)
where fL is the chord-length distribution, fL|X is the conditional probability distribution of chord length given a particle size assuming uniform orientation distribution, and fX is the size distribution. Alternatively, the measured CLD is assumed to be comprised of the CLDs of each of the particle sizes. The CLD for a given size is obtained by assuming uniform orientation distribution of a particle. Obviously, the transformation of the measured CLD to PSD will require a model for CLD for a given size, which is obtained by assuming that the location at which the chord is cut is uniformly distributed. While the CLD can be analytically derived for spheres and ellipsoid, one has to adopt a numerical approach to obtain the CLD for highly nonisotropic particles such as octahedral crystals with high aspect ratio and prismatic-shaped needles that are frequently encountered in pharmaceutical industries. Clearly, in the proposed literature approaches, the particle orientation effect on the CLD is not taken into account. However, in actual practice, the particle orientation must be biased by the flow of a fluid in which the particles are suspended. Such orientation bias becomes particularly serious for nonspherical particles with high aspect ratios. The CLD would be specific to the orientation, and this CLD is to be integrated over only the possible orientations determined by the local flow. It is the determination of distribution of orientations of particles caused by flow that is addressed in the present paper.
The major motivation of the present work is the need for converting chord-length measurements into particle sizes where the particles are significantly different from spheres with a high degree of anisotropy. The FBRM method of determining chord lengths depends on rotating a laser beam around an axis parallel to the laser beam, thus sweeping a cylindrical surface with length equal to that traversed by the beam in the apparatus and radius equal to the distance of the axis of rotation from the laser beam. The particles that happen to be on the cylindrical surface swept by the laser will be intercepted by the rotating laser beam. It is assumed here that the laser beam will intercept the particles sufficiently fast so as to keep the orientation of the particle constant during this period. As the laser beam intercepts the particle, it subtends a chord whose length is equal to the time of interception multiplied by the linear velocity with which the laser beam moves. The use of this FBRM technique is predicated on the assumption that there will be enough information available to calculate the particle size from the chord length. Obviously, the calculation of the particle size from the chord length depends on the particle shape and what one construes to be the size from that shape. Once this is decided upon, it is assumed that the orientation of the particle should suffice to calculate the geometrical relationship between the particle size and the subtended chord. The orientations of the particles, however, depend upon the local flow characteristics of the fluid in which they are suspended. It must be clear at this stage that the prediction of particle-size distribution from the chord-length distribution requires knowledge of the orientation distribution. The work proposed here is toward identifying the orientation distribution of particles subjected to an FBRM device at a particular location with flow conditions supposedly remaining invariant with time, at least in the averaged sense. In order to describe the orientation of the particle, we will assume the following: (a) The particle has an axis of symmetry whose orientation could be used to describe the orientation of the particle itself. (b) For the sake of generality, the specification of orientation of the particle on the cylinder surface somehow allows the exact calculation of the particle size from the chord length and vice versa. (c) The orientation of the axis of symmetry (and, hence, that of the particle) is described by three Eulerian angles (Θ, Φ, Ψ), which obviously can vary from 0 to 2π. (d) The subtended chord length can be calculated using geometrical arguments from specification of the following: particle size X, its orientation (Θ, Φ, Ψ), and a few other geometrical details, abstractly represented by vector Y (which can be anywhere in a domain ΩY of finite volume measure VY). An example of the vector Y may be given by the location of the particle on the cylindrical surface swept by the rotating laser beam, the point of intersection of the axis of symmetry with the cylinder in terms of its distance from the “centroid”, and so on. Let X denote the particle size and L denote the chord length that the laser beam subtends by intercepting the particle of orientation (Θ, Φ, Ψ) and location on the surface given by vector Y. The assumption that specification of the particle size and the particle orientation allows the calculation of the chord length is equivalent to postulating the existence of a function G(X,Θ,Φ,Ψ,Y) such that the random variables L and X are related as follows L ) G(X,Θ,Φ,Ψ,Y). We let the probability density for the distribution of X, i.e., the size distribution, be given by fX(x) and that for the
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3043
Figure 2. Algorithm to estimate the chord-length distribution function, fL|ΘΦΨX for a specified orientation and size.
distribution of L be given by fL(l), i.e., the chord-length distribution. Further, we let the conditional orientation distribution for a particle of size x be given by fΘΦΨLY|X(θ,φ,ψ,l,y|x). Then, from the total probability theorem, it follows that
fL(l) )
∫0∞ dx∫02π dθ∫02π dφ∫02π dψ × ∫Ω fΘΦΨLY|X(θ,φ,ψ,l,y|x)fX(x) dy
(2)
Y
writing fΘΦΨLY|X(θ,φ,ψ,l,y|x) as
fΘΦΨLY|X(θ,φ,ψ,l,y|x) ) fL|XΘΦΨY(l|x,θ,φ,ψ,y)fΘΦΨY|LX(θ,φ,ψ,y|l,x) Clearly, the orientation distribution is dependent only on the size of the particle and we may further assume that
fΘΦΨY|LX(θ,φ,ψ,y|l,x) ) fΘΦΨ|X(θ,φ,ψ,|x)fY|X(y|x)
(3)
so that we obtain from eqs 2 and 3 the following:
fL(l) )
∫0∞ dx∫02π dθ∫02π dφ∫02π dψ × ∫Ω fΘΦΨLY|X(θ,φ,ψ,l,y|x)fX(x) dy
(4)
Y
Figure 3. Computational algorithm schematic of two-step procedure to transform the chord-length distribution to size distribution for nonspherical particles with orientation bias: (a) Step I, determination of orientation bias function; (b) Step II, determination of PSD from CLD knowing the orientation bias.
From geometric considerations, we clearly have
fL|XΘΦΨY(l|x,θ,φ,ψ,y) ) δ(l - G(x,θ,φ,ψ,y))
(5)
If chord-length distributions fL(l) are available for a variety of situations in which size distributions fX(x) are also simultaneously obtained by, say, image analysis, then the inVerse problem would be to determine the conditional orientation distribution fΘΦΨ|X(θ,φ,ψ|x), provided we are able to estimate the conditional density fY|X(y|x). We assume that this distribution is uniform in the domain ΩY. More detailed considerations in regard to this assumption are not within the scope of this theory, but clearly its imposition at this level is much more reasonable than at the level of the orientation of particles. The assumption enables the calculation of fL|ΘΦΨ(l|θ,φ,ψ,x), which is demonstrated by an example below. 3.1. Determination of fL|ΘΦΨX(l|θ,O,ψ,x). The following subsections progressively discuss the geometrical considerations in the determination of the conditional probability distribution function in question. 3.1.1. Laser Surveillance Zone. We consider an absolute Cartesian coordinate frame of reference OsX1X2X3, with the X3 pointing along the direction of the laser beam and the origin
O being at the axis of rotation and also where the laser beam enters the surveillance zone. In other words, the laser enters the surveillance zone at X3 ) 0. We also entertain cylindrical coordinates (r, R, z) obtained through transformation as follows.
()
r ) xX12 + X22, R ) tan-1
X2 , z ) X3 X1
(6)
The laser surveillance zone is regarded as the region 0 < X3 < lz and in the “vicinity” of the cylindrical surface r ) RL (traced by a laser with a radius of rotation RL) that depends on the optical properties of the suspension and orientation and size of the particle, respectively. 3.1.2. Model for the Particle Shape and Size. The particles are assumed to be geometrically similar with ξ as the single characteristic size dimension. We assume that the particle shape can be described by an equation for the particle, given by
f(x,a) ) 0
(7)
where x are coordinates of a point on the particle surface and a is a parametric vector governing the particle size. The interior
3044
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
of the particle is denoted by the set I ≡ {x:f < 0}, while the exterior is given by E ≡ {x:f > 0}. We let the surface be denoted by the set S ≡ {x:f ) 0}. The origin of the coordinates associated with the particle is (symmetrically) inside the particle and is defined by
f(0,a) ) -1
(8)
The geometric similarity implies that the foregoing function is invariant when the particle coordinates and the parameter vector are scaled by the same scalar. In other words, f(λx,λa) ) f(x,a). An example is the superquadrics4 defined by
f(x1,x2,x3;a1,a2,a3) ≡
() () () x1 a1
m
+
x2 a2
m
+
x3 a3
m
-1)0
(9)
The particle size ξ is defined to be some functional of the function f and is represented by ξ ) Ξ(f). An example is the diameter of a sphere that is equivalent to the particle in that it has the same surface-to-volume ratio, which can be written as a (nonlinear) functional of f. It is also of interest to describe the particle domain in terms of absolute coordinates. We let the absolute coordinates of the particle centroid be represented by (X1,o, X2,o, X3,o) or by (ro, Ro, zo) related to the Cartesian coordinates by the above transformation (eq 6). The orientation of the particle coordinate axes (x1, x2, x3) relative to absolute coordinates (X1, X2, X3) is given in terms of Eulerian angles θ, φ, and ψ. Clearly the transformation is given by
[] [ ]
x1 X1 - X1,o x2 ) L X2 - X2,o , xi ) Lij(Xj - Xj,o) x3 X3 - X3,o
[
(10)
L≡ c(φ)c(ψ) s(θ)s(φ)c(ψ) + c(θ)s(ψ) -c(θ)s(φ)c(ψ) + s(θ)s(ψ) -c(φ)s(ψ) -s(θ)s(φ)s(ψ) + c(θ)c(ψ) c(θ)s(φ)s(ψ) + s(θ)c(ψ) s(φ) -s(θ)c(φ) c(θ)c(φ)
]
where s and c are sine and cosine of their arguments, respectively, and the absolute coordinates of the origin of the particle coordinate axes are (X1,o, X2,o, X3,o). In eq 10, we have used the summation convention for repeated indices. 3.1.3. Equations for Particle Surface in Terms of Absolute Coordinates. Since the relationships between absolute and relative (particle) coordinates are known above, we can write the equation of the particle shape in absolute coordinates. Alternatively, we may write
f(x) ≡ F(X - Xo) ) 0 where we have dropped the parametric vector a from the notation for simplicity. If we further transform to cylindrical coordinates of points on the particle surface, we let
Xo f (ro, Ro, zo) and X f (r, R, z) Further, we let g(r,R,z,ro,Ro,zo) ≡ F(X - Xo). It is now of interest to see how a particle with a fixed orientation in space can intersect a cylindrical surface of specified radius from the X3- or (z-) axis. Let the cylindrical surface be oriented with its central axis along the z- axis. The cylinder surface is defined by r ) RL, 0 < z < lz (the laser surveillance zone). By moving the centroid of the particle along the line (Ro, zo), the particle is
Figure 4. (a) Two test size distributions: Case I, s; Case II, - - -; (b) and (c) imposed ODF contours at θ ) 0c and θ ) 2.4c, respectively.
being displaced radially with its orientation fixed. This is accomplished by varying the value of ro by keeping the values of Ro and zo fixed. First, we consider the case of ro > RL, and we examine the equation F(X - Xo) ) 0, restated as g(r,R,z;ro,Ro,zo) ) 0. As we radially move the particle without changing its orientation, the particle surface eventually touches the cylinder surface. This obviously occurs at some point with r ) RL. Let us denote the corresponding values of R and z by R1 and z1, respectively. Thus (RL, R1, z1) represents the coordinates of the “upper” point of tangency between the particle and the cylindrical laser sweep
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3045
Figure 5. Simulated chord-length distributions: (a) Case I and (b) Case II.
Figure 6. Inverted size distribution along with the input PSD: (a) Case I; (b) Case II; input PSD (s); after inversion (O).
surface and, consequently, ro ) R h . On the other hand, by moving the centroid radially inward along (Ro, zo), i.e., ro < RL, we can find ro ) R at the “lower” point of tangency between the particle and the cylindrical surface. We can calculate R and R h in a single effort by asking that r ) RL, and we can solve for two different solutions of ro. To derive the required equations for this step, recognize that the particle surface and the surface r ) RL must be tangential at the required points. At the points of tangency, the normal vectors to the two surfaces must be a multiple of each other. Thus, ∇g ) β∇(r - RL), which yields two equations, ∂g/∂z ) 0 and ∂g/∂R ) 0, that are clearly free of unknown β. The third equation is obviously given by g ) 0. The solution to these three equations will yield directly the two values R and R h for the variable ro along with the values of (z1, R1) and (z2, R2) that belong to the two points where the particle will intersect the cylindrical surface swept by the laser beam. Once the values of R and R h are known, the range
specified orientation of the particle since R and R h depend on the particle orientation. We call this range of the coordinates as Y-coordinates. Clearly this is the spatial region that a particle can take in order to yield a detection of chord by the laser. In the absence of the means of exact determination of the location of particle that gives rise to a detection of chord, we will assume that Y-coordinates are uniformly distributed in the domain given by ΩY in eq 11. 3.1.4. Calculation of Chord Length. Knowing the geometric coordinates (ro, Ro, zo) of the particle centroid, which clearly is an element of the set of Y-coordinates that can be identified for a specified orientation angles in accordance with the algorithm described in an earlier subsection, it is straightforward to identify the end points of the chord (RL, R1, z1) and (RL, R2, z2). These in turn can be easily used to calculate the chord length, which is given by RL|∆R| where ∆R is the difference between the angular solutions obtained for the two intersections. Thus, given the orientation angles, the geometric coordinates of the particle centroid can be generated to yield the chordlength distribution, fL|ΘΦΨX(l|θ,φ,ψ,x). The algorithm for the determination of the conditional chord-length distribution has been shown in Figure 2.
ΩY ) {R < ro < R h ; 0 < Ro < 2π; 0 < zo < ls}
(11)
of the geometric coordinates, ro, Ro, and zo, is now completely known. It should be noted that this domain depends on the
3046
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
Figure 7. Contours of arbitrary ODF chosen to observe its effect on the inversion of CLD to PSD: (a) θ ) 0c and (b) θ ) 2.4c.
Thus, the overall model for the FBRM chord-length distribution can be succinctly stated as follows,
∫Ω
∫0∞ fL|OX (l|o,x) fO|X(o|x)fX(x) dx
(12)
Figure 8. Inverted size distribution along with the input PSD: (a) Case I; (b) Case II; input PSD (s); after inversion (×).
where fL|OX is the chord-length distribution for a given size with a specific orientation and fO|X is the probability density function for the orientation of the particle with a specified size. Furthermore, o and Ωo are the orientation vector and the domain of the orientation vector, respectively. The proposed model obviously suggests a two-step strategy to transform the chord-length distributions to size distributions at a given location of measurement, which is depicted in Figure 3. The first step would be to fix the FBRM location and measure the CLD for the suspension of precharacterized PSD. It must be ensured that the particles do not undergo any physical/ chemical change during the measurement. The characterization of the PSD in terms of its size distribution and shape of the particle may be done a priori by, say, image analysis. Perform the inversion to get the ODF. Now this ODF can be used in subsequent determination of unknown PSD given the CLD at the same location of the measurement. It is clear that the first step is to be performed only once as long as the location of measurement remains the same.
can be performed following usual methods reported in literature.1,2 All the calculations were carried out using Matlab.5 The two test cases of widely different size distributions for ellipsoidal particles with the shape parameters in eq 9 as a1/ a2/a3 ) 1:1:6 and m ) 2 are depicted in Figure 4a. The arbitrary normal ODF chosen for the purpose of verification is also shown in the same figure. In the absence of the actual experiment, the chord-length distribution was simulated for both the test cases using the assumed ODF (parts b and c of Figure 4). The simulated CLDs have been shown in Figure 5. Thus, the FBRM data was imitated for a known particle size distribution (PSD) and the ODF was obtained by the inversion as per the procedure outlined earlier. The extracted ODF was then used along with the simulated CLD (pretended as if it was measured by FBRM device) to obtain the PSD. The PSDs resulting out of the inversion have been shown in Figure 6, which indicate good agreement with the PSD. Parts a and b of Figure 8 show the PSDs obtained after the inversion of CLD by assuming some arbitrary ODF depicted in parts a and b of Figure 7. It is clear that, without the knowledge of correct ODF, one is likely to get erroneous PSD.
4. Verification of the Methodology
5. Summary The two-step methodology for the transformation of the chord length to size distribution for the nonspherical particles sus-
fL|OX (l) )
do o
The orientation bias function can be expressed in terms of three-dimensional Fourier series expansion, and the inversion
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3047
pended in the flow biasing the orientation of the particles has been developed and successfully demonstrated for the arbitrary size distributions with known orientation distributions. The importance of the first step in which the ODF is to be determined for a given location of measurement is illustrated through the comparison of the inverted PSDs from the CLDs with different ODFs. Acknowledgment Sponsorship by Eli Lilly and company for the project on the modeling of crystallization kinetics is gratefully acknowledged. Notations fL ) chord-length distribution fL|X ) conditional probability distribution of chord length given a particle size assuming uniform orientation distribution fL|ΘΦΨX ) conditional probability density function for chord length given the orientation for a given size fΘΦΨ|X ) probability distribution of orientation angles for a particle of given size fX ) size distribution fL|OX ) chord-length distribution for a given size with a specific orientation fO|X ) probability density function of the orientation of a particle for a specified size fY ) Y-coordinate distribution function g ) particle surface RL) radius of the cylinder swiped by laser (i.e., laser rotation) (ro, Ro, zo) ) cylindrical coordinates of the particle center l ) chord length lz ) length of laser surveillance zone L ) rotation matrix
o ) orientation vector consisting of the Eulerian angles R h ) radial coordinate of the particle center lying outside the swiped cylinder still touching the cylinder R ) radial coordinate of the particle center lying inside the swiped cylinder still touching the cylinder Greek Letters R ) angular coordinate in cylindrical coordinates β ) constant in the equality of the gradients in equation: ∇g ) β∇(r - RL) (θ, φ, ψ) ) Eulerian angles ξ ) characteristic particle size defined by ξ ) Ξ(f) Ωo ) space of orientation angles ΩY ) Y-coordinate space Literature Cited (1) Worlitschek, J. Monitoring, Modeling and Optimization of Batch Cooling Crystallization. Ph.D. Dissertation, Swiss Federal Institute of Technology Zurich, Switzerland, 2003. (2) Li, M.; Wilkinson, D. Determination of Non-spherical Particle Size Distribution from Chord Length Measurements. Part 1: Theoretical Analysis. Chem. Eng. Sci. 2005, 60, 3251. (3) Tadayyon, A.; Rohani, B. Determination of Particle Size Distribution by ParTec_100: Modeling and Experimental Results. Part. Part. Syst. Charact. 1996, 15, 127. (4) Barr, A. H. Superquadrics and Angles-Preserving Transformations. IEEE Comput. Graphics Appl. 1981, 1, 11. (5) MATLAB, release 13; The Mathworks: Natick, MA, copyright 19842002.
ReceiVed for reView July 20, 2006 ReVised manuscript receiVed October 5, 2006 Accepted October 9, 2006 IE0609463