Mixing and Fluid Dynamics Effects in Particle Precipitation Processes

Precipitation defined as the rapid formation of moderately soluble crystalline or amorphous solid particles from a liquid solution under high supersaturation conditions is considered. It involves the simultaneous and fast occurrence of primary nucleation and growth of particles together with the secondary processes such as aggregation and breakage. It is shown how the effects of fluid flow and mixing affect the subprocesses forming the overall precipitation process. Examples are presented for reactive precipitation and antisolvent precipitation with supercritical fluids applied as the antisolvents. The elementary subprocesses forming the overall precipitation process (i.e. macro-, meso-, micromixing, chemical reaction, nucleation and growth of particles, aggregation and breakage) are characterized by the related time constants and the application of time constants in the modeling and scale-up of precipitation processes is presented. The effects of turbulence on particle formation are simulated using mechanistic models, CFD and population balances (including the method of moments and the quadrature method of moments). The results of modeling are compared with experimental data. Finally, it is shown in the context of practical applications to what extent the approaches discussed in this paper can be applied to “design” particles.


Introduction
Precipitation refers to the rapid formation of moderately soluble crystalline or amorphous solid particles from a liquid solution, whereas particle formation occurs under high supersaturation conditions. Precipitation involves the simultaneous and fast occurrence of primary nucleation and growth together with the secondary processes such as Ostwald ripening, aggregation and breakage.
Precipitation processes are of importance in technology as they are employed in the production of bulk and fine chemicals including fertilizers, pigments, catalysts, magnetic materials, drugs, etc. The solid product usually has a wide crystal size distribution (CSD) which determines the filtration, washing and settling abilities of suspensions, and thus the quality of the product (crystal size homogeneity, surface area, etc.). Other factors which can strongly influence the product quality are the morphology and purity of the precipitating particles. It is well known that both the CSD and the morphology can strongly de-pend on mixing conditions during the process.
The high supersaturation is usually created either by performing chemical reactions or by mixing the solution with an antisolvent which decreases the solute solubility. Chemical reactions applied in so-called reactive crystallization processes are usually very fast, which means that their course depends on mixing as well. This also means that the role of mixing, including the effects of macro-, meso-and micromixing, can be very important for the creation and dilution of supersaturation. Very important are also effects of the fluid flow as they affect mixing and may affect phenomena of particle agglomeration and coagulation; agglomeration dominates the precipitation processes at high supersaturation, in the case of coagulation only physical forces are involved to hold particles together.
Another aspect of high supersaturation is its effect on the nucleation mechanism. Namely at high supersaturation, there is negligible secondary nucleation resulting from the presence of solute crystalline material but rather primary homogeneous or heterogeneous nucleation dominates the process. Due to the high sensitivity of primary nucleation to supersaturation, the number of primary particles is controlled by nucleation and affected by mixing, whereas crystal growth is less important, because after complete unloading of supersaturation, it is more the number of crystals than the rate of their growth that determines their final size. The growth rate can, however, affect their shape.
In this paper, the effects of mixing and fluid flow observed in reactive precipitation, and antisolvent precipitation with supercritical fluids applied as the antisolvents are considered. The material is organized as follows.
First, the elementary subprocesses forming the overall precipitation process (i.e. macro-, meso-, micromixing, chemical reaction, nucleation, and growth of particles) are characterized by the related time constants and length scales. The application of time constants for the modeling and scale-up of precipitation processes is then presented.
Second, because in most of the precipitation devices the flow and mixing are turbulent, it is shown how to predict interactions of turbulence with particle formation phenomena. In this context, the phenomenological models, mechanistic and based on CFD (including solution of the closure problem for chemical reaction and precipitation) and population balances (including method of moments and quadrature method of moments) are presented and the results of modeling are compared with experimental data.
Third, an attempt to classify aggregation, i.e. agglomeration and coagulation processes, is presented from the point of view of the characteristic length and time scales.
Fourth, suitable kernels for aggregation processes (from orthokinetic to perikinetic, including effects of fluid dynamics and colloidal forces) and breakage processes are presented, and the results of application of these kernels, implemented in some cases in a CFD environment, are presented and compared with experimental data.
This paper finishes with discussion on practical applications where the properties of precipitated particles are of key importance and thus great attention is paid to "particle design". It is also discussed to what extent the approaches presented in this paper can be applied to "design" particles and scale-up processes. This is a review paper and the intention of the author is to review important problems, ideas and challenges related in the author's opinion to the most important aspects of precipitation. The approaches and examples of results presented are based on the author's own work, and supported by the ideas and results of other researchers who investigated this fascinating problem of precipitation.
Despite its practical importance and a long tradition of investigating precipitation, it still requires better understanding and development of theoretical approaches for process interpretation. This is, for example, depicted by a small number of books devoted to precipitation. There are some fundamental books on precipitation including Kinetics of Precipitation by Nielsen A.E. (1964), The Formation and Properties of Precipitates by Walton A.G. (1964) and Precipitation by Söhnel O. and Garside J. (1992); there are also chapters devoted to precipitation in the Crystallization Technology Handbook, edited by Mersmann A. (2001), Handbook of Industrial Crystallization edited by Myerson A.S. (2002), Turbulent Mixing and Chemical Reactions by Bałdyga J. and Bourne J.R. (1999), Supercritical Fluid Technology for Drug Product Development edited by York, P. et al. (2004), Handbook of Industrial Mixing, Science and Practice edited by Paul E.L. et al. (2003) and Formulating Poorly Water Soluble Drugs by Rowe, J.M. and Johnston, K.P. (2012); discussion on computational methods that can possibly be used to model precipitation can be found in Computational Models for Polydisperse Particulate and Multiphase Systems by Marchisio D.L. and Fox R.O. (2013).

Driving force for precipitation
Precipitation is a kinetic process and its rate depends on the appropriately defined driving force, called supersaturation in the case considered here. The supersaturation, σ, is defined as the dimensionless difference between the molar chemical potential (the molar Gibbs free energy) μ in supersaturated solution and μ eq in the crystalline state (equal to the potential in saturated solution) that can be expressed using solute activities, a, (Davey R. and Garside J., 2000).
where μ and a are the actual molar chemical potential and the actual activity of the solute in solution, respectively, whereas T is the absolute temperature. The value μ -μ eq = 0 defines equilibrium: for μ -μ eq > 0 spontaneous crystallization may occur, and μ -μ eq < 0 means that an opposite transformation, i.e. dissolution, becomes spontaneous. The sign of Δμ determines the direction of the possible phase transformation, whereas the absolute value of Δμ signifies the distance from the equilibrium. Δμ/RT represents the thermodynamic driving force for crystallization, and both the rate of nucleation and the rate of crystal growth depend on Δμ/RT. The driving force is usually expressed in a simpler form than Eqn. (1): eq a eq eq eq eq ln 1 1 a a a S a a In chemical engineering, concentration is often expressed in the units appropriate for the specific application, e.g., molar concentration, c c c,eq eq a a c S S a c    where a c represents the activity of dissolved component of concentration c, c eq is the equilibrium solubility and S a ≈ S is the supersaturation ratio that is also known as relative supersaturation or saturation ratio. Supersaturation is also sometimes expressed as a concentration difference eq c c c    (4) The solution is supersaturated when S a > 1, S > 1, Δc > 0 or σ > 0. For electrolytes, which dissociate giving v + cations and vanions, v = v + + v -, and where   where D is the diffusion coefficient of the solute, k B is Boltzmann's constant, d is the molecular diameter, N* is the number of molecules or ions forming the critical embryo-nucleus, σ is the surface energy and v is the molecular volume. Nielsen A.E. (1964), Nývlt J. (1971), and Söhnel O. and Garside J. (1992) have shown that over limited ranges of supersaturation, Eqn. (7) may be simplified to where n, the kinetic 'order' of nucleation, varies between 1 and 10. Heterogeneous nucleation results from presence of foreign substances that are always present in a solution and may act as 'catalysts' of nucleation. This catalytic effect results mainly from decrease the energy barrier to nucleation, het G G      thanks to smaller σ for nucleus forming on the surface of the foreign phase. Hence, rate of heterogeneous nucleation (Turnbull D. and Vonnegut B., 1952) can be presented in similar way as homogeneous nucleation, Eqn. (7).
where Ω het < Ω. Heterogeneous nucleation is observed much more frequently than homogeneous nucleation. In true secondary nucleation, nuclei are either formed on the surface of the solid phase as dendrites that are subsequently broken off the crystal or in the surface proximity (adsorption layer, liquid adjacent to the crystal); in contact secondary nucleation, nuclei result from the effects of impact energy of the crystal-crystal and crystalequipment collisions. In the case of precipitation of moderately soluble materials, the secondary nucleation is not an important mechanism of nuclei formation because the particles are too small for the high-impact collisions that would be able to produce secondary nuclei, and the effects of crystal-solution interactions are negligible (Söhnel O. and Garside J., 1992). It is only in the case of rapid agglomeration that the resultant particles can be large enough to induce a sufficiently high impact energy for the generation of nuclei.
Once nuclei are formed in a supersaturated solution, they start to grow into crystals. The rate of crystal growth can be presented as the rate of growth of each crystal face. In the case of precipitation, particles are small and thus a concept of the overall rate of growth, G = dL/dt is used with the crystal volume expressed by v c = k V L 3 , k V being the volume shape factor. The rate of growth is affected by two processes: the mass transfer of growth units to the crystal-solution interface by bulk diffusion and the surface integration of growth units into the crystal lattice. This can be represented by where k D = 2k d /ρ c , k d is the mass transfer coefficient, k R is the surface integration constant that is independent of crystal size and r represents the "order" of the surface integration process. In Eqn. (10), c is the bulk solute concentration, the subscript i describes the concentration at the crystal-solution interface and the subscript eq indicates the concentration at equilibrium. When the resistance of mass transfer to the crystal 1 D k  is much lower than that of the subsequent surface integration step , the surface integration process controls the overall rate of growth, G = k R (c -c eq ) r . In the opposite case, the overall rate of growth is determined by the rate of mass transport from the solution to the crystal G = k D (c -c eq ). When the growth is controlled by combined mechanisms, then the concentration at the crystalsolution interface c i in Eqn. (10) has to be eliminated, in this way relating the rate of growth to the bulk and equilibrium solute concentrations.
A similar approach can be used in the case of more complex systems when there is, for example, a reaction between the constituent ions on the crystal surface before integration into the lattice. As an example, consider precipitation of the barium sulfate BaSO 4 . Vicum L. et al. (2003) applied thermodynamic models for aqueous Ba 2+ , 2 4 SO  , Na + , Clsolutions differing in composition.
The dependence of solubility on composition is then expressed using the thermodynamic solubility product where c i represents the concentration of ions i with γ i being the activity coefficient of the ionic species i, and γ ± represents the mean ionic activity coefficient that depends on the ionic strength.
The mean ionic activity coefficient γ ± can be calculated using semi-empirical forms of the Debye-Hückel law, either that proposed by Bromley L.A., (1973), applicable for the ionic strength up to 6M, or the Pitzer K. (1991) model, valid for higher values of the ionic strength. The equilibrium between ions and undissociated ion pairs is expressed by the dissociation equilibrium, with K a being the dissociation constant The driving force for precipitation, S, is then expressed by: The impact of the composition on the driving force, in-cluding its effect on complex formation, is presented in Fig. 1 with k r = 9.1 × 10 -12 ms -1 (Vicum L. et al., 2003) and k D = 4.0 × 10 -5 m 4 s -1 kmol -1 (Bałdyga J. et al., 1995).

Time constants for precipitation and mixing
It has been well known since years (Bałdyga J. and Bourne J.R., 1999) that the history of supersaturation and the history of supersaturation structure determine the PSD and the crystal morphology in precipitation.
Before starting scale-up or modeling the precipitation processes, it is good policy to assess the importance of mixing in these processes. One should consider the time constants for all mixing and precipitation phenomena that affect the distribution of supersaturation and compare the characteristic time constants for these processes (Bałdyga and Bourne, 1999).

Time scales for precipitation
A comparison of the characteristic time constants for precipitation and mixing should yield the rate of the controlling mechanism. In the following section, the time constants that can be used to characterize precipitation  Vicum L. et al., 2003) are presented. The first one is the time constant for the n th -order chemical kinetics of homogeneous reaction, where k n is the rate constant of the chemical reaction of order n and c A0 is the initial or feed concentration of reactant A. This time constant characterizes the rate of consumption of reactant A and thus the rate of creation of supersaturation provided that the chemical reaction controls this process. For a more general definition and more exact identification of τ R , the method given by Bałdyga J. and Bourne J.R. (1999) is recommended; for example, for reaction of the type For the primary nucleation, the time constant can be approximated by the induction period, i.e. the time that elapsed between mixing on the molecular scale two reacting solutions and the appearance of the first crystals of precipitate; the induction period is considered to be inversely proportional to the rate of nucleation (Dirksen J.A. and Ring T.A., 1991).
where K is a constant. The time constant τ N represents a characteristic time that is necessary to form a characteristic, significant number of nuclei per unit volume. The characteristic time for crystal growth can be expressed in relation to the rate of product concentration decrease caused by the crystal growth where c c is the concentration of precipitating substance, and A g is the specific surface area per unit volume of suspension.
Note the different role of processes characterized by the time constants presented above; the time constants τ R and τ Gcr characterize processes that in a similar way to mixing directly affect the distribution of supersaturation in the system, but only τ Gcr depends on supersaturation itself. On the other hand, the nucleation process that is characterized by τ N significantly depends on the supersaturation distribution but does not directly affect this distribution.

Time scales for mixing
The time constants for precipitation should be correlated with the process time and time constants for mixing.
The process time represents the batch time in the case of a batch precipitator and can be identified with the mean residence time τ in the case of the continuous flow systems.
Chemical reactions, nucleation and crystal growth are essentially the molecular-molecular-level process, which means that mixing on the molecular scale can directly influence the course of precipitation. However, the largerscale mixing mechanisms often indirectly affect the precipitation process by changing the environment for local mixing. Hence the complete sequence of mixing processes, namely macromixing, mesomixing and micromixing, needs to be considered.
Macromixing can be identified as fluid blending on the scale of the precipitator. The characteristic length scale for this process is thus the size of the precipitator. Macromixing in the tank of volume V t can be identified with the macroscopic flow pattern; e.g. the blending time τ b or cir- can be used as simple characteristics of macromixing in a stirred tank. In Eqn. (19), N represents the stirrer frequency and D S is the stirrer diameter.
There are two mechanisms of mixing on a meso-scale.
The first mechanism of mixing on the meso-scale refers to the coarse-scale turbulent exchange between the fresh feed and its surroundings. The time constant τ D for meso-scale mixing by turbulent diffusion can be expressed by where Q is the feeding rate, D t represents turbulent diffusivity, and ū is the local, average fluid velocity. The related length scale can be estimated either as (Q/ū) 1/2 for high feeding rates or can be represented by the feed pipe diameter d f for low feeding rates, when (Q/ū) 1/2 < d f . The second mechanism of mesomixing refers to disintegration of large eddies of a size larger than the Kolmogorov microscale, λ K = v 3/4 /ε 1/4 , but smaller than the integral scale of concentration fluctuations, Λ C . Inertial-convective disintegration proceeds without being influenced by molecular-scale mixing, however, it itself affects the micromixing by changing the local environment. The time scale τ S for meso-scale mixing resulting from the inertial-convective disintegration of the large eddies reads where A is the proportionality constant equal roughly to 1.2. Micromixing processes are driven by the mechanism of the viscous-convective deformation of fluid elements, followed by molecular diffusion. A characteristic time constant for viscous convective mixing, i.e. for eddies smaller than the Kolmogorov microscale scale, λ K is known as the engulfment time constant (Bałdyga J. and Bourne J.R., 1989): that for a high Schmidt number can be replaced by the time constant for the viscous-convective and viscousdiffusive mixing Based on these characteristic times, Bałdyga J. and Bourne J.R. (1999) proposed a classification of precipitation by comparing the characteristic times for mixing and reactive precipitation. Using τ M to denote the characteristic mixing time (this can be any mixing time listed above), one can classify the precipitation subprocesses from the point of view of competition between precipitation and mixing. Therefore the reaction can be instanta- . The effects of mixing should be observed in the first two cases.
It should be noted that the chemical reaction can be classified as slow compared to micromixing, but fast or instantaneous when compared to mesomixing or macromixing. For an instantaneous, mixing-controlled reaction (i.e. for τ R much smaller than any mixing time scale) but with τ N > τ VCVD , the influence of molecular diffusion on nucleation can be neglected. However, when τ N ≈ τ E , the viscous-convective micromixing (engulfment) affects precipitation. Usually, τ N << τ Gcr -thus mixing affects precipitation mainly through the nucleation step.
Note also that for a sufficiently slow feeding τ >> τ C , one can assume a uniform concentration of the bulk and just model the reaction zone using meso-and/or micromixing modeling. If precipitation is slow compared to meso-and micromixing, the assumption of ideal mixing can be used. Considerations similar to the ones presented above can be helpful in choosing the model of mixing that is most suitable for the modeling and scaling up of precipitation processes.

Scale-up principles for precipitation
Scaling up industrial apparatus represents a fundamental step in the realization of industrial plants. The scale-up procedures are, however, not obvious in the case of complex processes, including precipitation. To scale up the process, one can use an empirical approach, i.e. carry out the process at a number of scales and gain enough information to make an empirical prediction of the system performance on a larger scale, i.e. extrapolate the process using principles of similarity. Dimensionless similarity criteria can be derived using the available governing equations for mass, momentum, species, energy, and population balances; when they are not available, the Buckingham's Pi Theorem can be used instead (Zlokarnik M., 2002).
In the geometrically similar systems, a complete similarity occurs if all necessary dimensionless criteria derived either from the differential equations or obtained using the Pi Theorem are equal. For example, to have a similar velocity distribution in a single-phase system, it is enough to apply the same Reynolds number in both systems. In complex precipitation processes, such complete similarity is impossible and in fact not necessary. We usually want to reproduce the product quality, for example, the particle size and morphology at the larger scale and obtain identical not similar products from the systems differing in scale. Thus, one should apply a partial similarity and in this way obtain a reduced number of the most important similarity criteria for the process.
A more fundamental approach to scale-up is based on modeling with the goal of predicting the effect of scale-up and of understanding the process in sufficient detail.
The strategy recommended here is to start by performing an analysis of time constants to identify which sub-processes are less important on the large scale, and gain information on which requirements related to such processes can be neglected. Afterwards, a physical model employing CFD should be applied to predict the scale-up effects more correctly and finally to optimize the process.
As an example of scaling up with a combined application of time scale analysis and CFD, let us consider an antisolvents precipitation method known as the solution-enhanced dispersion by supercritical fluids process (SEDS TM ) that has been proposed by York P. and Hanna M. (1996). In the SEDS TM process, the supercritical fluid (SCF) acts as an antisolvent and the substrate solution is co-introduced into a nozzle and mixed intensively with the antisolvent in a nozzle chamber, from where they both flow into the precipitation vessel at high velocity as a turbulent jet, as schematically depicted in Fig. 2. The process is carried out above the mixture critical pressure, so the antisolvent (CO 2 ) is completely miscible with the solvent (ethanol). As mentioned earlier, the distribution of supersaturation and the supersaturation history of fluid elements determine the PSD and the crystal morphology in precipitation, which means that macro-, meso-, and micromixing processes can affect the particle size distribution.
To scale a process up or down, one needs invariability of all characteristic times: τ = idem, and τ i = idem from Eqn. (15) to (23), which is of course impossible.
Details on the scale-up of the SEDS process can be found in Bałdyga J. et al. (2010). Interpreting the precipitation vessel shown in Fig. 2 as a pipe-in-pipe system, with the outer tube fluid velocity equal to zero, one can characterize a similarity of the flow pattern using the Curtet number (Becker H.A., et al., 1963;Craya A. and Curtet R., 1955), which characterizes turbulent selfentrainment.
where d 0 is the nozzle orifice diameter and D v is the diameter of the precipitation vessel. This means that under process conditions, there is a recirculation eddy induced as shown in Fig. 3, and the flow similarity requires d 0 /D v = const. The recirculation eddy causes the exhausted, residual fluids to dilute the fresh supersaturated solution, which significantly reduces the nucleation rate.
The supersaturation is then unloaded mainly due to the particle growth.
Using the jet nozzle diameter d 0 for the length scale, and the nozzle velocity u 0 for the velocity scale, we express the rate of energy dissipation by represents the dimensionless distributions of the rate of energy dissipation.
This means that both the mean residence time τ and the characteristic convection time τ C can be expressed by and the fluid elements should have the same dimensionless positions 0 x d  after the same dimensionless time period t/τ C .
The time constants for meso-scale mixing then read The time constants for the viscous-convective and viscous-diffusive micromixing processes scale as ε -1/2 see Eqns. (22) and (23), which means that when the fluid properties are the same in both scales considered. Hence, the ratio of any mesomixing time constant and micromixing time constant scales as . This shows that at large Reynolds number values, τ S >> τ E , with mesomixing controlling the process and negligible micromixing effects. Based on this, Bałdyga J., et al. (2010) proposed the following: (1) for developed turbulence, characterized by the high Re 0 values, micromixing effects are negligible, the scale-up criterion is then (d 0 /u 0 ) 1 = (d 0 /u 0 ) 2 and thus for weakly turbulent flows micromixing dominates, and the scale-up cr iter ion is then τ E1 = τ E2 or and thus (u 0 ) 2 /(u 0 ) 1 = X 1/7 and (d 0 ) 2 / (d 0 ) 1 = X 3/7 . As a safe scale-up method, the criterion based on high Reynolds number flows has been recommended.

Fig. 3
The structure of the flow in the SEDS particle formation vessel, represented by the stream lines (Bałdyga J. et al., 2010(Bałdyga J. et al., , copyright 2010 Elsevier Ltd., reprinted with permission). To validate this criterion, CFD simulations and experimental investigations for acetaminophen precipitation were carried out for three equipment scales: a laboratory scale with d 0 = 2.0·10 -4 m, a pilot scale with d 0 = 4.0·10 -4 m, and a small manufacturing plant (SMP) scale with d 0 = 9.0·10 -4 m. Fig. 4 shows the importance of the geometric similarity; the particle size distributions obtained for the same process conditions, but different vessel diameters and thus different C t values, are shown in Fig. 4.  Fig. 5 shows the results of application of the scale-up criterion (d 0 /u 0 ) 1 = (d 0 /u 0 ) 2 . The micromixing effects decrease with increasing system scale, following the relation Re 02 /Re 01 = (d 02 /d 01 ) 2 , and are negligible in the case of both the pilot plant and the SMP results presented in Fig. 5.  Fig. 6 shows an example of the CFD prediction of the supersaturation distribution in the small manufacturing plant; a similar distribution was observed in the case of the pilot plant, whereas in the case of the laboratory scale, the supersaturation was about 7 % higher due to a too small Reynolds number.
The importance of the time constants, either for the di-rect modeling-based scaling up or for the formulation of scale-up rules, has been emphasized in several publications. For example, Zauner R. and Jones A.G. (2000) have shown that the calcium oxalate precipitation can be scaled up by using the time constants for micro-and mesomixing in combination with a simple segregated feed model. An interesting discussion on using time constants is presented by Gavi E. et al. (2007), in the context of the scale-up of Confined Impinging Jet Reactors (CIJR) based on the results of Marchisio D.L. et al. (2006).

Concept of population balance
The balancing of populations is related to classic problems of statistical mechanics that were defined in the 19 th century, whereas the fundamentals of the methods of treating populations of aggregating particles were introduced by Smoluchowski M. (1917 where x i are the external coordinates representing positions in the physical space and x j are the internal coordinates representing properties of particles, n is the population number density, u pi refers to the component of the rate of change of particle position in the physical space, m is the number of particle properties, and w pj refers to the component of the rate of change of the particle property, and B and D are the birth and death functions at the point in the phase space. Very often, the property of the particle is identified with the particle size, L, and then Eqn. (29) takes the form Eqn. (30) can be used under assumption that the shape of the particles remains unchanged during particle growth, and the concept of constant shape factors can be used, with the particle volume v p = k v L 3 , the particle surface area a p = k a L 2 , and the characteristic size L chosen so that the ratio of the shape factors k a /k v = 6. However, this  is not always the case. When, for example, the particle shape is determined by the growth of two faces, Eqn. (29) should be used as Eqn. (29) should be solved together with the mass, momentum, energy, and species balances to determine the particle velocity in the physical space as well as the rates of change of the particle properties. This creates numerical difficulties because the population balance equation should be solved in (m + 3)-dimensional space, whereas all other equations are defined in the 3-dimensional physical space.
To reduce the dimensionality, the volume-averaged population balance can be used by applying either a mechanistic model of micromixing, (Bałdyga J. et al., 1995) or another macroscopic domain including the entire equipment, e.g. the crystallizer (Randolph A.D. and Larson M.A., 1971). This reduces the number of coordinates to the number of particle properties, m.
In the literature, for m > 1, one can find mostly publications dealing with such volume-averaged populations. Ramkrishna D. and Mahoney A.W. (2002) considered the application of the multi-dimensional population balances to represent the dynamics of particle morphology. They discussed the possibility of finding precipitation kinetics from the dynamics of the particle size distribution by the inversion of the multidimensional population balances (i.e. by solving the inverse problem). Efficient numerical techniques to treat the problem numerically were considered by Qamar S. and Seidel-Morgenstein A. (2009). Briesen H. (2006) tried to simulate the crystal size and shape by means of the population balance model based on the crystal volume and the shape factor as internal coordinates. Ochsenbein, D.R. et al. (2015ab) considered the growth and agglomeration of particles including effects on their shape.
To reduce the problem to the 3-dimensional physical space, one can use the moment transformation of the population balance. For example, for particles described by one internal coordinate L, Eqn. (30) can be used, so for m = 1, one gets the set of equations represented by Eqn.
for i = 1, 2, 3 and j = 0, 1, 2, ... and   , , n x t L  represents the crystal size distribution at the time t and position x  , and B fj and D fj are the birth and death terms for the moments of order j. The averaged birth and death terms for the moments of the order j are not closed as the form of the CSD is not known in advance. This creates the first "closure problem" that can be solved using an idea of the presumed distribution. A good example here is the quadrature method of moments (McGraw R., 1997).

Application of CFD
Consider the process of precipitation of a solid product from two liquid ionic solutions A and B A n+ + B n-→ C↓ involving an instantaneous chemical reaction that produces a supersaturated solution of the product C, and its subsequent crystallization-precipitation, in the systems with turbulent flow.
A quantitative description of this process requires a complete knowledge of the spatial distribution of the instantaneous local concentration of ions A n+ and B n-, product C, and other ions present in the system to properly express the driving force. Identification of these concentrations requires the solving of a system of partial differential equations describing the momentum, mass, species, and population balances in the system.
The Reynolds averaging procedure, typical for treating turbulent flows, produces more unknowns than the number of the balance equations, and the set of the partial differential equations becomes unclosed. We thus have the second "closure problem". This "closure problem" is usually solved by forming a "closure hypothesis". The "closure hypothesis" introduces additional information based on the experimental data, investigator's experience, or intuition. They are included in CFD models such as the k-ε model, Reynolds Stress Model or Large Eddy Simulation model. To simulate the course of precipitation, CFD should be combined with the averaged population balance equation. In the simplest case it can be the standard k-ε turbulence model (Launder B.E. and Spalding D.E., 1972) which employs the hypothesis of gradient diffusion for momentum and species transport.
Eqn. (32) employs the local, instantaneous values of the velocity and concentration. Application of Eqn. (32) to the case of turbulent flow by utilizing the Reynolds decomposition and the Reynolds averaging (denoted with overbar) was obtained by Bałdyga J. and Orciuch W. (2001).
After introducing a concept of gradient diffusion for particles and employing a concept of mixture fraction f that can be interpreted as a volume fraction of the nonreacting fluid originating in the A stream when A and B solutions are mixed, one gets Eqn. (34) in the closed form where D pT is the local value of the particle turbulent diffusivity coefficient. The statistics of fluid elements is expressed using the probability density functions, PDF, Φ( f ), for example the presumed beta distribution (Bałdyga J. and Bourne J.R., 1999). The distributions of the average concentrations f and the inert tracer concentration variance, are necessary to define the presumed PDF and can be calculated using CFD and one of the mixing models: the multiple-time-scale model of the turbulent mixer (Bałdyga J., 1989) or the spectral relaxation model (Fox R.O., 1995). Eqn. (35) should be solved together with the partial differential equations describing the species balances.
The results of application to the SEDS process of the model presented above were already shown in Fig. 6. More results are presented in Fig. 7. The method of SEDS modeling is described in detail by (Bałdyga J. et al., 2004(Bałdyga J. et al., , 2008a and (Henczka M. et al., 2005). Fig. 7 explains how one can control the particle size by varying the flow rate (expressed here by Re); this effect resulting from the competition between the creation of supersaturation by mixing fresh fluids and its decrease by dilution with residual fluids was considered in Section 2.3 for scaling up precipitation processes. Fig. 7 also illustrates the competition between meso-and micromixing. The comparison of CFD models, the first model neglecting and the second including micromixing effects shows that micromixing dominates the process at small Reynolds number turbulent flows, whereas mesomixing dominates at high Re. It is interesting that predictions of the mechanistic E-model (Shekunov B.Y. et al., 2001) are very good and illustrate very well the effects of creation and dilution of supersaturation due to mixing.
Other methods applying the presumed PDF are based on the application of groups of the Dirac delta functions representing different environments (Piton D. et al., 2000, Marchisio D.L. et al., 2001, or combining DQMOM with the IEM model of micromixing (Gavi E. et al., 2007). A full PDF method was applied by Di Veroli G. and Rigopoulos S. (2010). The comparison of these three approaches, based on the presumed beta-PDF, multienvironment, and the full PDF, was carried out by Fissore D. et al. (2002). The following examples present the effects of mixing on the precipitation of BaSO 4 in the single-feed and double-feed stirred tank reactors operating in semibatch mode. Fig. 8 presents a visualization of the mixing zone for the single-feed simulations, which is regarded as the zone of the intensity of segregation, larger than 0.01. The intensity of segregation represents the normalized variance of concentration of the inert tracer. Clearly, the micromixing is not instantaneous and not localized at the feed point. Fig. 9 compares the measured and predicted effects of the initial concentration on the particle size for the singlefeed mode. The effect of concentration is predicted well, especially at the higher concentration values.
In the case of the double-feed precipitation, the effects of mixing are more significant. Segregation of the fresh feeds A and B prevents the creation of supersaturation, and thus micromixing decreases segregation, increases supersaturation and the nucleation rate close to the feed pipe.
At a further distance from the feed point, micromixing with the residual solution from the bulk decreases supersaturation and limits nucleation. The local intensity of the segregation of fresh fluids A and B is presented in Fig. 10, showing the zone where mixing affects the precipitation process. Fig. 11 presents a comparison of two models: model I neglecting fluctuations of species concentrations, thus neglecting the effect of the inertial-convective and viscousconvective concentration fluctuations, and based on the local average concentrations, with complete model II including the effects of micromixing. There is a satisfactory agreement of the experimental data with the model predictions only in the case of model II.
Neglecting the micromixing (model I) results in a very large amount of very small particles, due to the fact that very high supersaturation values are created in the region of the impeller which-due to the strong nonlinearity of Eqn. (7)-results in a very large amount of small particles. The time constant for nucleation (about 10 -5 s), is much smaller than the time constant for micromixing that-close to the impeller-is of the order of milliseconds. This clearly shows that at high supersaturation, the micromixing effects need to be included to match the simulation results with the experimental data. Fig. 9 Comparison of the model predictions with the experimental results for single-feed mode with the feed introduced above the impeller (Bałdyga J. et al., 2006, reprinted with permission). Model II accounts for the micromixing effect. Fig. 8 Visualization of the mixing zone (Bałdyga J. et al. 2006) for the single-feed mode for N = 5s -1 and t f = 1230 s. (Bałdyga et al. 2006, reprinted with permission).

Fig. 10
Mixing in the feed zone of a double-feed precipitator operating at N = 3s -1 and q f = 0.12 × 10 -6 m 3 s -1 3D representation of the feed zone for the intensity of segregation I S ≥ 0.01 (Bałdyga et al. 2006, reprinted with permission).

Fig. 11
The effect of the stirrer speed on the mean crystal size for q f = 0.12 × 10 -6 m 3 s -1 , , feed concentrations c A0 = c B0 = 0.234 kmol m -3 and α V = 50. Solid line-predicted by Model II accounting for the micromixing effects; dashed line-predicted by Model I neglecting the micromixing effects (Bałdyga et al. 2006, reprinted with permission).

Time and length scales for aggregation
According to Söhnel O. and Garside J. (1992), aggregation denotes all processes of clustering of separate particles to form larger particles, whereas in the agglomeration process, strong secondary particles are formed that are held together by crystalline bridges or physical forces. In the case of coagulation, only very small particles are involved and only physical forces hold them together. Let us now discuss what we mean by small particles, and let us try to classify them based on their size. In the case of solidliquid particulate systems, particles of a size smaller than 1 μm are called colloidal, whereas the larger ones-suspended (Elimelech M. et al., 1995).
In the case of a dispersion subjected to turbulent flow, the classification should be based on the comparison of the particle size expressed by the particle radius a with the characteristic length scales for turbulence.
We apply two characteristic scales of turbulence: the Kolmogorov microscale λ K = v 3/4 /ε 1/4 , the smallest scale of turbulent flow, and L e = (u′) 3 /ε, the scale of the large, energycontaining eddies; here u′ is the root-mean-square velocity fluctuation.
Particles of the size a, smaller than the Kolmogorov length microscale, λ K , are classified as small. The very small particles are groups of the small particles for which  , so their movements are dominated by the Brownian motion, which leads to perikinetic collisions and perikinetic aggregation. The Péclet number is defined as the ratio of the diffusion time

B
a D  to the convection time 1/E. E represents the rate of deformation here and B D  is the coefficient of Brownian diffusion for a particle separated from other particles.
In the group of large particles (a > λ K ), there are inertialconvective effects; aggregation of very large particles, a > L e , is hardly possible in practice. The time-scale analysis is also important to make the description more complete as it considers more mechanisms which affect aggregation. The time scales characterizing particles are: the diffusion time τ B , and the relaxation time τ p , respectively: where 12 D  is the mutual Brownian diffusion coefficient for large particle separation, C D is a friction factor (also known as a drag coefficient), ρ p is the particle density, the Reynolds number, Re p , is defined using the relative velocity fluid-particle and the equal volume sphere diameter. They can be compared with the time scales for turbu-lence: the Kolmogorov time microscale, τ K , and the Lagrangian integral time scale, T L that can be estimated using Eqn. (40) (Hinze J.O., 1975).
where Λ f is the longitudinal integral scale of turbulence, u′ is the r.m.s. (root-mean-square) velocity fluctuation, and β ≈ 0.4. For τ p < τ K , there is no correlation between particle motions if τ B < τ K and the motions are correlated for τ B > τ K . For τ K < τ p < T L and τ B > T L there is a partial correlation of particle motions; for τ p > T L and τ B > T L there is no correlation between motions of the particles.
One can see that it is not enough that particles are small a < λ K , to understand the mechanism of aggregation because either perikinetic or orthokinetic aggregation can be observed depending on the relation between time scales or related values of the Péclet number. Similarly in the case of large particles, a > λ K , with partial or no correlation of particle motions, the appropriate kernels including inertial effects should be applied.

Aggregation
The aggregation rates can be considered either in the Lagrangian or the Eulerian frame of reference. The Lagrangian approach is based on the trajectory analysis of the particles and the Eulerian approach is based on the concept of the pair distribution function that mimics in a continuous way the particle concentration distribution. It enables identification of the probability of finding another particle (and its orientation if not spherical) in the vicinity of the considered particle.
In the Lagrangian approach, instantaneous values of the relative particle velocity are calculated from Newton's second law, including or not including Brownian motions. If they are included then the model takes the form of a stochastic Langevin-type differential equation with solutions representing stochastic trajectories. In practice, the concept of limiting trajectories is most often applied, which means that the random movements are neglected and the limiting trajectories define the collision cross-section. To include random movements, the methods of dynamic simulations can be used; a detailed description of these methods can be found in the book by Elimelech M. et al. (1995).
In the Eulerian approach, the Fokker-Planck equations are applied. Consider small particles with τ p << τ K . The particles are conveyed by the fluid, diffused by the Brownian motion, and are subjected to colloidal DLVO-type forces. They can be in the perikinetic τ B < τ K or orthokinetic τ B > τ K regime, or in the regime in between the two τ B ≈ τ K .
The collision rate constant, β c , can be calculated from the pair probability function, c pair , distribution: where f u  is the relative particle velocity resulting from the fluid flow, int u  is the relative particle velocity resulting from the colloidal interactions, S C is the entire collision surface and n  is the normal to this surface (Melis S. et al., 1999). Eqn. (41) can be simplified (Bałdyga J. and Orciuch W., 2001) by taking an ensemble average of Eqn. (41) and using a concept of the gradient turbulent diffusion, expressed in terms of a generalized form of the stability ratio W.
where V is the interaction potential, G(r,λ) is the hydrodynamic function describing the squeezing of the liquid out of the space between particles, r is the center-to-center distance, λ is the particle size ratio, and R D represents the overall (turbulent and Brownian) diffusivity ratio. Fig. 12 shows the stability ratio values calculated using the DLM model to predict the surface potential.
The predicted values of the surface potential are shown in Fig. 13. They have been calculated using either the Double Layer Model, DLM, Löbbus M. et al. (1998), or the Gel Layer Model, GLM, (Löbbus M. et al., 1998;Adler J.J. et al., 2001). One can see strong repulsion effects for pH > 6, with strong effects on the stability ratio, and almost no repulsion below pH = 6. Fig. 14 shows that at Pe → 0 (because ε → 0), we observe perikinetic aggregation with no effect of the rate of energy dissipation, whereas for large ε, the orthokinetic aggregation controlled by convection (Pe >>> 1) is observed. The effects of colloidal forces are observed at the perikinetic limit only and are negligible when the orthokinetic mechanism dominates. When the relaxation time is large enough, τ p ≥ τ K , the inertial effects cannot be neglected and then the rate of collision depends on fluctuations of the relative velocity of particles w, which depends on the fluid velocity fluctu-ations and relaxation time (Kruis F.E. and Kusters K.A., 1996) Because not every collision leads to agglomeration, the concept of the probability or efficiency of agglomeration, P S , can be used .
where c t represents the average time necessary to build a strong enough bridge between the particles and i t is the average interaction time.
The time period c t can be interpreted as the time nec- Fig. 12 The effect of the process conditions on the stability ratio, W, for 10 nm particles (Bałdyga J. et al., 2012(Bałdyga J. et al., , copyright 2012 Elsevier Ltd., reprinted with permission). where f(λ) represents the shape function as proposed by David R. et al. (1991). The results of the simulations show that at constant supersaturation, both an increase of the energy input and an increase of the particle size lower the probability of agglomeration during collisions. The agglomeration kernel β(L 1 ,L 2 ) can then be expressed as the product of the collision kernel, denoted here by β c , and the agglomeration probability, P S : The results of Bałdyga J. and Krasiński A. (2005) illustrate this phenomenon for the reactive, agglomerative precipitation of benzoic acid in a semibatch, stirred tank reactor. Supersaturation was created in this case by the homogeneous chemical reaction between sodium benzoate and hydrochloric acid. Figs. 15 and 16 illustrate main results.
The model correctly predicts the trends of the effect of the concentration and the stirrer revolutions observed in the experiments.
The microscopic pictures of the primary particles and  the aggregates confirm that the model correctly predicts both the size of the primary particles forming agglomerates and the size of non-agglomerated particles. Notice the effects of supersaturation on both: size of primary particles due to particle growth and size of agglomerates due to building strong bridges between primary particles forming agglomerates.

Breakage
Disintegration processes play a significant role in the chemicals industry. There are many types of devices (e.g. rotor-stators, high-pressure nozzles, stirred tanks) that can be applied to generate stresses. There are of course different mechanisms involved in stress generation and thus different mechanisms of deagglomeration. For example in the case of high-pressure systems, disintegration results from hydrodynamic and cavitation stresses, whereas in the case of the rotor-stator system, only hydrodynamic stresses are active. To select the most suitable equipment for the desired parameters of the final product, one can use relatively simple mechanistic modeling in combination with CFD. The modeling presented in this paper includes the effects of agglomerate structure and suspension rheology on the disintegration processes.
The strength of the agglomerate depends on its structure, which in turn determines the mechanism of its breakage. A typical agglomerate structure of Aerosil fumed silica particles is represented by the cluster-fractal model as shown on Fig. 17. Large clusters of size L i consist of smaller primary aggregates of a size L a , which are composed of primary particles of size L 0 . They are characterized by the fractal dimension D f . The primary particles forming primary aggregates are connected by strong chemical bonds. Such aggregates are strong, stable and resistant against mechanical forces. They are characterized by an average fractal dimension D f0 . The clusters of size L i , connected by adhesion forces, are relatively unstable, and they can be disintegrated by mechanical forces. Fractal geometry is used to estimate the tensile strength of agglomerates.
There is experimental evidence (Baldyga J. et al., 2008bc) that the break-up of Aerosil 200 V agglomerates occurs through erosion in a rotor-stator mixer, and through shattering and rupture in both a high-pressure nozzle disintegrator and the ultrasonic device. Break-up processes control the agglomerate size distribution; details of break-up depend on the agglomerate structure strength and the stresses generated in the fluid. On the other hand, the rheological properties of the suspension depend on the concentration of particles, their size distribution and morphology.
To describe all these phenomena and simulate the process, a model based on the population balance was developed. In the model, the population balance equations are supplemented with adequate breakage kinetics and solved  using the quadrature method of moments (McGraw, R., 1997) in the CFD environment. The commercial CFD code Fluent (the k-ε model of turbulence) is supplemented with UDFs for suspension rheology. Suspension rheology is described with the Buyevich Yu. A. and Kabpsov S.K. (1999) model. Models considered in this work were compared with experimental data in publications (Baldyga J. et al., 2008bc). In the simulations, two mechanisms of stress generation were considered: hydrodynamic stresses generated by the flow with frequency Γ h (L), and stresses generated by cavitation events with cavitation itself generated either by the flow or by the ultrasonic waves, at frequency Γ c (L) for agglomerates of size L. In the simulations, a superposition of breakage mechanisms is assumed.
Comparison of these time constants with residence time values in the region of stresses higher than the tensile strength of agglomerates characterizes the extent of disintegration. Typical results for the rotor-stator, high-pressure nozzle and ultrasonic device are presented below.
The geometry of the rotor-stator device is shown in Fig. 18. Fig. 19 shows the distribution of hydrodynamic stresses and particle size in the rotor-stator mixer.
The stresses are not strong enough to disperse agglomerates in one passage through the system, and thus to disperse agglomerates, the effective residence time should be increased and this is possible by multiplying the number of passages through the rotor-stator as shown in Fig. 20.
Much higher stresses (hydrodynamic and induced by cavitation) are generated in the high-pressure nozzle device (Fig. 21), and the system is so effective that one passage through the system results in effective dispersion of agglomerates (Fig. 22). Comparing the model with the experimental results, one can exclude the rupture mechanism; because the stresses are much higher than those in the rotor-stator mixer, one can expect that the breakage mechanism is shattering. More precisely this is represented by the fragmentation number defined in Rwei S.P. et al. (1990Rwei S.P. et al. ( , 1991; for the stress τ, it can be expressed in a more general form (Bałdyga et al., 2008b)  where τ is the stress, Ha is the effective Hamaker constant and H represents the surface-to-surface distance. Fractal geometry is used to estimate the tensile strength of agglomerates applied in Eqn. (49). The number of bonds connecting aggregates in an agglomerate can be related to the agglomerate porosity by means of the follow-ing equation (Rumpf H., 1962;Tang S. et al., 2001): where F represents the bonding force between aggregates

Fig. 19
Hydrodynamic stresses τ h and mean particle size L 30 distributions in the rotor-stator system for Q = 0.6 dm 3 /s, N = 9000 rpm (150 s -1 ) and 5 wt.% of Aerosil, for conditions applied by Bałdyga et al. (2008c). in the agglomerate that is calculated using the classic DLVO theory. In the present case, F is expressed by the van der Waals' attractive forces F A that are calculated assuming that all primary particles have equal diameters,  To complete the description, one should define the stresses τ for each mechanism of generation.
In the case of turbulent flow for aggregates smaller than the Kolmogorov microscale, λ K = (v 3 /ε) 1/4 , Eqn. (52) works for breakage kernel Γ with deformation rate defined by For turbulent flows, and for aggregates larger than the Kolmogorov microscale, L > λ K , one has When cavitation is induced by the flow then the breakage kernel is recalculated from the condensation rate R C using the model by (Singhal A.K. et al., 2002).
The rate of collapse of the bubbles serving as the breakage kernel is defined as a ratio of the mass condensation rate to the mass of a bubble: with C bc being an empirical constant and ρ v the vapor density. The bubble radius value, R 0 , is expressed using the Blake equation for the critical bubble diameter (Brennen, C.E., 2005). The stresses generated by bubble implosion (generated by shock waves and high-velocity microjets) can be estimated using the semi-empirical equation (Crum, L., 1988) where c is the velocity of compressional waves in the liquid, V j is the microjet velocity and α is a constant that varies from 0.41 to 3.0. The breakage kernel for ultrasonic-wave-induced cavitation and for ultrasonic frequency f, reads: for the ultrasonic wave propagation, we use the wave equation (Saez V. et al., 2005) Fig. 23 shows the experimental system, whereas comparisons of model predictions with experimental data is given in Fig. 24. The model predicts well the small effect of Aerosil concentration on the average agglomerate size. Fig. 25 shows the distribution of cavitation stresses generated in the ultrasonic cell. They are many orders of magnitude larger than the stresses generated in the rotorstator (Fig. 19a), and the deagglomeration process has either a rupture or shattering mechanism. Experimental results and model predictions for the deagglomeration of silica nanoparticle agglomerates (Aerosil 200 V) confirmed that the fragmentation number can be used to predict the mechanism of deagglomeration. Deagglomeration is possible for Fa > 1; for Fa < 100 erosion is observed, rupture is for Fa > 100 and shattering dominates for Fa > 1000.  (Bałdyga et al., 2008b(Bałdyga et al., , copyright 2008 The Institution of Chemical Engineers, reprinted with permission.)

Fig. 24
Comparison of model predictions with experimental data for one passage through the ultrasonic cell. Effect of Aerosil concentration for energy density E v = 300 GJm -3 . (Bałdyga et al., 2008b(Bałdyga et al., , copyright 2008 The Institution of Chemical Engineers, reprinted with permission.)

Final discussion
There are many practical applications where the properties of precipitated particles are of key importance. The most important reason for this is that there is an increasing emphasis on high added value specialty chemicals such as pharmaceuticals, dyes, pigments, paints, and printing inks regarding the crystal size distribution, particle morphology and product purity. As presented by Karpinski (2006), many precipitated pharmaceutical organic products can form a range of polymorphs or solvatomorphs. There is also the effect of the precipitate properties on the downstream processes such as filtration or drying. This explains why attention is paid to "particle design".
The results presented in this paper show that particle size distribution (Figs. 4, 5, 7, 9 and 11), particle morphology and many other properties depend on the spatial distribution of the supersaturation (Fig. 6) and the timing of its evolution. Important is also the evolution of fluid elements in the supersaturated zone which depends on the flow pattern (Fig. 3) and the rate of mixing (Fig. 10). Mixing controls both creation and dilution of supersaturation. However, not only the level but also the chemical structure of supersaturation is important, which can be shown using precipitation diagrams (Söhnel O., 1991); a simple example of such a diagram is shown in Fig. 1. One can plot on the precipitation diagram such important information on the type of precipitate that results from the solution composition represented by points on the diagram. The composition of the solution affects the chemical reaction and resulting properties of the solid phase. One can thus show those regions where a solid phase of uniform composition, morphology, color or crystal habit exists.
Another effect resulting from the solution composition is related to electrostatic interactions (Fig. 13) between colloidal particles that may control the stability of colloidal suspensions (Fig. 12) and particle aggregation (Figs.  14, 15 and 16).
Another important distribution is the distribution of hydrodynamic properties such as the velocity gradients, deformation and rotation tensors, and resulting hydrodynamic stresses. The relation between fluid deformation and rotation affects the efficiency of micromixing, thus also affecting the creation rate of supersaturation. Very important is the distribution of stresses generated in the system, including hydrodynamic stresses (Fig. 19) and the stresses resulting from cavitation (Fig. 25). Direct interaction between the distribution of hydrodynamic properties, the distribution of supersaturation and the distribution of the particle populations results in the size distribution and morphology of agglomerates as shown in Fig. 15 and Fig. 16. There is of course the direct effect of stress distribution on deagglomeration which is presented in Fig. 20, Fig. 22 and Fig. 24.
Precipitation is a very complex process, which means that when designing and scaling up precipitation processes, complete similarity is impossible. In fact when scaling up, complete similarity is not necessary because as explained earlier, we want to obtain identical products from the systems differing in scale. As a first attempt to design the process, the time and length scale analysis is recommended in this paper, which is illustrated with examples. The comparison of time and length scales shows which of them are important for the process. This information can be used for scaling up and for formulating adequate models. The precipitation models have a multiscale character and are very complex. It is thus very important to choose an adequate model of mixing and to decide at which scale modeling is not necessary.