Modelling and Simulation of Particle Size Distribution of Precipitates in Continuous Tubular Crystallizers

This paper presents a population balance model for describing the temporal evolution of the particle size distribution of a precipitate produced in a laboratory scale tubular crystallizer, including nucleation, growth and agglomeration of crystals. The quadrature method of moments is used to calculate the moments of the crystal size distribution. The gamma probability density function with the sixth order Laguerre polynomials are used to reproduce the particle size distribution from moments. A sensitivity analyses is carried out by simulation that can be helpful when designing a crystallizer. Influence of the design and process parameters on the particle size distribution of product is analysed.


Introduction
Precipitation is a widely used method to synthesize and purify materials.The basic concept of a precipitation process is to obtain compounds, ordinary by fast chemical reactions or salting-out phenomena, which are characterized by low solubility compared with the reactants.Under these conditions, a chemical compound appears in the fluid phase in much higher concentration that would be thermodynamically permitted.Due to this thermodynamic instability of the fluid phase a crystallization process occurs leading to a solid disperse precipitate and at the same time to desupersaturation and stabilization of the liquid phase.
In modern chemistry, beside purity such crystalline properties like the mean particle size and size distribution of particles also obtain significant attention [1][2][3][4].Particle size dis-tribution (PSD) can be tailored by operations like grinding and granulation.These secondary operations are needed because a lot of processing properties, for instance the surface specific properties like adsorption capacity, solubility, rate of dissolution or filtration and sedimentation, depend on PSD.Therefore, it seems to be rather useful to obtain product with the desired particle size distribution directly from reaction precipitation.
Micromixers of different types, such as impinging jet mixers [5,6], vortex-type mixers [7,8] and Y-mixers [9,10] are often used in controlled reaction and anti-solvent precipitation which ensure fast and good mixing of reactant solutions and may produce precipitates with desired particle size distributions [2,11].Micromixers, not depending on geometries, usually are combined with tubular reactors ensuring the subsequent plug flow transport of the crystalline suspension.Although, these devices are often used and studied there is still a lack of detailed mathematical models by means of which it would be possible to predict, control and optimize the precipitation processes.
The aim of the paper is to present a model based approach to predict the PSD of a precipitate produced in a laboratory scale tubular reactor following a Y-micromixer which ensures fast and complete mixing of reactant streams.The chemical reaction is considered instantaneous and the process is assumed to 59(2), pp.138-144, 2015 DOI:10.3311/PPch.2181Creative Commons Attribution b be isothermal [12].Both primary and secondary nucleation, growth and agglomeration of crystals are included.To visualize the effects of constitutive and process parameters of the crystallizing substance and reactor setup on the PSD of precipitate a detailed sensitivity analysis is presented and analysed.

Mathematical model
The precipitation process to be modelled is considered according to the scheme in Fig. 1.The Y-micromixer is followed by a tube which transports the homogenous stream produced in micromixer into a collective vessel.In the micromixer, an instantaneous chemical reaction takes place with primary nucleation A+B→C↓ while the subsequent processes, i.e. further primary nucleation, crystal growth, secondary nucleation and agglomeration takes place in the tube in which plug flow of suspension is assumed.The output of the tube is a solid-fluid suspension having a specific PSD.
As a consequence, the mathematical model consists of the following equations.
Mass balance equation for the Y-micromixer having volume V m : subject to the initial condition c m (0) = 0, where c A denotes the input concentration of reactant A, S = c m /c s is the supersaturation ratio with the solubility c s , and L n denotes the linear size of nuclei, q denotes the flow rate and M is the molecular mass of the crystallizing species.In Eq.(1) B p stands for the primary nucleation rate written in the form while the second term on the right hand side denotes the rate of consuming of the precipitating substance C by nucleation.
The mass balance equation for species C in the tubular crystallizer, assuming plug flow of suspension in that, is written as the partial differential equation subject to the c(x, 0) = c 0 (x) initial and c(0,t) = c m (t) boundary conditions.Here, x denotes the axial coordinate of the tube, v is the linear velocity of plug flow and k V µ 3 denotes the total volume of crystals in a unit volume of suspension, expressed using the third order moment of the linear crystal size In Eq.(4), n(x,L,t) denotes the population density function of crystals population by means of which n(x,L,t)dL provides the number of crystals from the interval of size (L,L+dL) at time t and at axial coordinate x in a unit volume of suspension.
The population balance equation governing the temporal evolution of the population density function is written in the form subject to the initial condition n (x, L, 0) = 0. Here, the boundary conditions are written as where B s is the secondary nucleation rate, written as depending on the total surface A c (x) of the precipitate in a unit volume of suspension and denotes the growth rate of crystals.In Eq. ( 9) the total surface A c (x) is expressed by means of the second order moment of the linear crystal size with the surface shape factor k A of crystals.
, , , The second and third terms on the left hand side of Eq.( 5) describe, respectively, the rates of change of the population density function due to flow of crystals along the tube and crystal growth while the two terms on the right hand side represent the rates of decreasing and increasing of the number of crystals because of crystals agglomeration with kernel function β.

Solution by means of the quadrature method of moments
The complex system of equations ( 1), ( 3) and ( 5) is solved by using the quadrature method of moments rewriting those, taking into account the time scales of the process components, into quasi-stationary form.
According to the assumption of instantaneous chemical reaction with high rate primary nucleation due to very high supersaturation in the Y-micromixer, we can write dc m /dt ≈ 0 so that from Eq. ( 1) we obtain simply the algebraic equation In this way, a well stirred suspension characterized with high supersaturation is fed into the plug flow tubular crystallizer therefore we can assume also quasi-stationary state in the tube thus Eqs. ( 3) and ( 5) are written, respectively, in the form subject to the initial condition c(τ = 0) = c m , as well as in the form subject to the initial and boundary conditions Note that the second term on the right hand side of Eq. ( 14) denotes the disappearing of particles due to agglomeration and the third term takes into consideration that in each agglomeration step a new particle, agglomerate is formed.
Introducing the quadrature form of moments of the linear size of crystals and substituting (19) into Eq.(18) we obtain which together with Eq. ( 13) forms a closed set of ordinary differential equations.Note that the quadrature weights (w) and abscissas (L) have no physical meaning.Those are used only to calculate the moments of the particle size.
For describing the effects of turbulence in the tubular crystallizer the turbulent shear agglomeration kernel [18] is used where ε denotes the dissipation rate of turbulent kinetic energy and ν denotes the viscosity of suspension.The dissipation rate is expressed as ,1, 2... ) , , 0,1, 2... (20) where p denotes the pressure, A is the cross-section area of the tube.The viscosity of the sus-pension is given as where φ denotes the volume fraction of the solid phase.

Simulation results and discussion
The set of ordinary differential equations formed by Eq. (13) and by Eqs.(20) for the first six leading moments was solved in Matlab environment.A detailed sensitivity analysis was carried out studying the effects of the kinetic parameters (k p0 ,b p ,k s0 ,b s ,k g ,g) and the parameters of the reactor (X, d, q) on the crystal size distribution of the resulting precipitate.The basic collection of parameters used in the sensitivity analyses is listed in Table 1 and alteration of parameters the effects were actually studied were modified with respect to these values.
Table 2 contains the limiting values of parameters between which those were altered.In each case five simulation runs were performed changing the values of parameters uniformly between the extreme values.In each case the basic value is in the middle of the interval so in each case one of five simulations was performed by using the basic parameter set.This method made the sensitivity analysis results more comparable.
Note that the kinetic parameters used in the current sensitivity analysis typically hold for the reaction precipitation.
The PSD can be approximated knowing the leading moments of distribution by means of different techniques.In this study, the gamma density function corrected by the Laguerre polynomials was used [19].Accordingly, the population density function is given as where while the n th order associated Laguerre polynomial is expressed as Figure 2 illustrates the effects of changing the coefficient of crystal growth rate on the PSD computed on the basis of Eqs.(24)-( 27) from the moments of particle size.Increasing the growth rate coefficient shifts the support interval of the crystal size distribution toward the larger crystals but the dispersion is likely also increased.Computation of the mean value and standard deviation confirm these observations as it is indicated in Table 3.This alteration caused an approximately linear increase in the mean particle size but this apparent linear tendency could not be generalized because the system is highly nonlinear.
This behaviour is illustrated also in Fig. 3 showing the effects of increasing the supersaturation dependency of the secondary nucleation rate.The first two steps did not cause considerable change in the PSD but the third step already did as it is seen clearly.Starting from the fourth step the effects of the same step size became rather significant.Explanation of this behaviour can be found in the nonlinear characteristics of the process.While the lower values of b s were applied the primary nucleation and crystal growth were the determining phenomena but above a certain level a relatively small change of b s started to play more important role in determining the process.
Table 1 The basic values of model parameters used in the sensitivity analysis  Figure 4 presents evolution of the number of crystals, represented by the zero order moment μ 0 of the crystal size along the length of the tubular crystallizer depending on the growth rate coefficient.It is seen that the relatively high crystal growth rate compared with the nucleation rate may prevent further nucleation of crystals.It is also seen that primary nucleation appeared to be the dominant crystal production process in the tubular crystallizer although in some cases observable secondary nucleation arose at the final stage of the process.Only a weak effect of agglomeration might be observed indicating that the agglomeration rate proved to be too small regarding the conditions provided by the remaining kinetic and process conditions.It is also seen in Fig. 4 that after a certain length the number of crystals remains approximately constant in the reactor so no further nucleation occurs.The reason is that the supersaturation, i.e. the driving force of crystallization is consumed by the growth of particles.We see the same phenomenon when the mean particle diameter is investigated as it is illustrated in Fig. 5.
Figure 5 presents evolution of the Sauter mean diameter (m 3 / m 2 ) along the length of the crystallizer.Comparing the diagrams of number of crystals in Fig. 4 and of the Sauter mean diameter in Fig. 5 indicates that both the production of crystals and the Sauter mean diameter achieve their stationary states almost at the same crystallizer length.This means that under the kinetic and process conditions applied in simulation the supersaturation ratio at these moments becomes so small that even further crystal growth becomes negligible.
Influence of the model parameters studied in simulation on two characteristic parameters of the crystal size distribution, i.e. on the mean size and standard deviation is presented in Table 3 in a qualitative manner.Here, the meaning of the symbols are: "+": slight increase; "++": strong increase; "-": slight decrease, "--": strong decrease.
Table 3 illustrates well that the parameters of crystallizer significantly smaller influence on the PSD than the kinetic and constitutive parameters of the precipitating species.The effects of the design of crystallizer on the PSD manifest mainly in increasing agglomeration of crystals which, as it seen in Table 3, is a function of both the contact time and flow regime.

Conclusions
A model based sensitivity analysis was presented for a continuous precipitation reactor by computer simulation.The investigated characteristics were, in principle, the particle size distribution which was described by using a population balance model of the crystallizer including chemical reaction, nucleation, growth and agglomeration of crystals.The chemical reaction was assumed to be instantaneous and the model was built

Table 3
Qualitative assessment of the effects of model parameters on the PSD

Fig. 2 3 Fig. 5 Fig. 4
Fig. 2 Effects of variation of the growth rate coefficient on the PSD Fig. 3 Effects of variation of the exponent b s of secondary nucleation rate on the PSD

Table 2
The limiting values of model parameters used in the sensitivity analysis