Monitoring of Batch Industrial Crystallization with Growth, Nucleation, and Agglomeration. Part 1: Modeling with Method of Characteristics.

This paper develops a new simulation model for crystal size distribution dynamics in industrial batch crystallization. The work is motivated by the necessity of accurate prediction models for online monitoring purposes. The proposed numerical scheme is able to handle growth, nucleation, and agglomeration kinetics by means of the population balance equation and the method of characteristics. The former offers a detailed description of the solid phase evolution, while the latter provides an accurate and efficient numerical solution. In particular, the accuracy of the prediction of the agglomeration kinetics, which cannot be ignored in industrial crystallization, has been assessed by comparing it with solutions in the literature. The efficiency of the solution has been tested on a simulation of a seeded flash cooling batch process. Since the proposed numerical scheme can accurately simulate the system behavior more than hundred times faster than the batch duration, it is suitable for online applications such as process monitoring tools based on state estimators.

Introduction

Batch crystallization is an important operation for formulation and separation of high value-added chemicals in crystalline form from liquid solutions in pharmaceutical, food, agriculture, and fine chemical industries. By means of crystallization we are interested in achieving quality targets for the final solid product such as high purity and yield, morphology, and crystal size distribution (CSD). These properties are determined by complex crystallization kinetics of growth, and nucleation, as well as agglomeration and breakage, which in turn depend on the crystallization operating windows (mainly supersaturation, temperature profile, and mixing). Hence, optimal operating conditions strategies should be followed to guarantee a more efficient and environmental-friendly achievement of the specification targets. To attain this task a process engineer can rely on modern computer-aided techniques. To this end, rigorous simulation models of the system at hand are a valid tool for process understanding, process design, design of experiments, evaluation of process alternatives, and online control and monitoring.

From the perspective of the crystallization operation supervision, the accurate online monitoring of industrial crystallizers is still considered a challenge due to the limitations of hardware sensors for the measurements of the key variables such as solute concentration and CSD. These limitations include high investments costs, measurement delays, and measurements errors because of calibration difficulties.[1,2] Hence, there is much to be gained by using online model-based technologies such as state estimators, also known as soft sensors, that are able to predict the performance of the operation within a desired accuracy based on the process model and measurements of secondary variables.

In a recent study on the estimability properties of batch crystallization systems Porru and Özkan[3] found that the performance of the CSD estimation relies on the accuracy of its model since the CSD is not distinguishable neither with measurements of secondary variables such as temperature and volume nor concentration measurements.

In the light of the above-mentioned arguments, the goal of this work is to develop a simulation model of the CSD evolution in industrial crystallizers usable for real time monitoring objectives. The model has to satisfy two requirements:

it must be as detailed as required for the modeling of growth, nucleation, and agglomeration kinetics in order to cover for an appropriate prediction of the undistinguishable dynamics of the CSD

it must be computationally efficient to predict the CSD in real time.

A rigorous process model for crystallization systems consists of material and energy balances for the liquid phase, and material balances for the crystalline phase. The natural framework for the modeling of the crystalline phase in terms of CSD is the population balance equation (PBE).[4−6] The PBE is a partial differential equation (PDE) which describes the dynamics of the particle phase space allowing the complete description of the properties of the particle distribution. More specifically, the PBE is a conservation balance of the number of particles in the crystallizer, which is a function of time and space, where the spatial domain may include an external and internal coordinate. The external coordinate identifies the location of the crystals in the crystallizer. The internal coordinate represents the characteristic size of the particles. Normally, the chosen internal coordinate is the crystal length L rather than the crystal volume v.[7]

It must be pointed out that obtaining the solution of the PBE is challenging because it is analytically intractable.[7] Analytical solutions have been only found for the most idealized systems. For example Gelbard and Seinfeld[8] provide analytical solutions when the seeds distribution belongs to the exponential distribution and the internal coordinate is expressed in terms of v. The analytical solution proposed by Pinar et al.[9,10] is based on the limiting assumption that the CSD does not change its shape during the process.

With the advancement of computer science, the numerical solution of the PBE has become very appealing. However, the accurate numerical solution of the PBE is often challenging due to numerical diffusion,[11] and instability problems.[12] A variety of numerical techniques have been proposed. They can be divided into four main classes: the finite difference method, the method of weighted residuals, Monte Carlo methods, and discretization techniques. A complete overview of these strategies can be found in Mesbah et al.[12] and Costa et al.[13] It is adequate to say that the discretization techniques consist of the discretization of the internal coordinate domain along a grid of length (or volume) classes, and the substitution of the PDE with a system of ordinary differential equations (ODEs). According to Mesbah et al.,[12]

the discretization techniques have the best numerical accuracy,

among these discretization techniques, the finite volume method and the method of characteristics (MOC) are particularly amenable for online applications because computationally efficient algorithms can be built based on them.

Regardless of the numerical technique employed, the majority of the simulation studies presented in the literature only takes into account growth and nucleation kinetics,[12,14−17] while only a few works address the simulation of the agglomeration phenomena. As a result, these studies are often unable to describe industrial crystallization operations in a realistic way. The model of the aggregation of parent crystals[4] to obtain agglomerates has been first developed in volume basis since the total volume of the crystals is conserved in case of pure agglomerating systems.

The simulation of the agglomeration rates presents computation complexities related with the presence of integrals of nonlinear functions. However, when a numerical solution based on a discretization technique is used, it has been found convenient to treat the CSD as a truly discrete distribution and to reformulate the agglomeration mechanism with discrete equations.[18] Many mechanisms have been proposed: Batterham et al.[19] propose an ad-hoc model applicable to a geometric discretization scheme so that the volume classes v obey to v/v = 2. The main drawback in this approach is that this type of grid becomes very coarse at large sizes, which produces a consistent overprediction of CSD at large volumes. Aligned with this approach, the more recent works of Kumar and Ramkrishna[20−22] are remarkable. They first propose a discrete mechanism for the agglomeration, that is applicable to uniform grids,[20] observing that the accuracy of the prediction depends on the density of the grid points, but this demands high computation time. To overcome this problem, they then develop an algorithm that automatically changes the grid: a fine grid is added when needed, and a coarse grid is added at smaller sizes.[21] Finally, they simulate a system undergoing growth, nucleation, and agglomeration by combining the proposed model and the method of characteristics (MOC).[22] Results have a good fit with the available analytical solutions. However, the CPU times required by this technique remain too high[23] to make the algorithm appealing for online monitoring purposes.

Even if the agglomeration mechanism proposed in the literature[21,22] seems to be very accurate, its use is limited to those cases for which the number density can be described on a volume basis. In fact, the number density would be expressed in length units if the kinetics parameters are estimated from CSD data given on a length basis.[24] Indeed, data from focus beam reflectance measurement probes measure a chord length distribution (CLD), which is given in length bases. CLD is then converted in CSD by means of geometric inverse modeling.[14,25] CSD measurements in terms of length also naturally arise from images techniques[25,26] due to the inherent 2D feature of the pictures. In the case of length units used, the agglomeration term has a more complex structure including the integral of a function which is nonlinear in the internal coordinate. This makes the method proposed by Kumar and Ramkrishna[20−22] difficult to extend.

To the best of our knowledge few papers have proposed a numerical solution for agglomeration expressed in length basis. Marchal et al.[27] proposed an approximation of the agglomeration term based on the application of the mean value theorem on frequency. This method does not conserve the particle number and gives poor performance.[20,28] Hounslow et al.[18] proposed a method that consists of a discrete reformulation of the PBE, but it is applicable to grids which obey to the geometric series L/L = 21/3, L being the ith length class. This type of grids becomes very coarse for large sizes. In fact, to cover a length range from 0.1 to 1300 μm, only 42 classes would be modeled. Moreover only 5 classes would cover the range from 500 to 1300 μm: 516, 650, 820, 1032, and 1300 μm, respectively. Hence, the use of the algorithm proposed by Hounslow et al.,[18] would lead to a poor estimation of the number density function at large crystal dimensions. Majumder et al.[23] use an algorithm based on the Lattice Boltzmann method. However, the CPU time required by this method seems not to be compatible with real-time applications, being comparable to the process duration.

In this work, we propose a novel numerical solution with the method of characteristic (MOC) in order to obtain an accurate and computationally efficient solution for PBE accounting for the growth, the nucleation, and the agglomeration phenomena together expressed on a length basis. The proposed method includes an accurate computation of the integral constituting the agglomeration term. The method is applicable to any CSD shape, and does not require any specific grid for the length domain, overcoming the limitations of the algorithm proposed by Hounslow et al.[18] With the proposed method we are able to simulate the system behavior more than hundred times faster than the batch run duration (with a computer with Intel(R) Core (TM)2 Quad CPU 2.83 GHz and MATLAB 2016a), which, differently from other methods found in the literature,[23] makes it appealing for online use for monitoring purposes.

This paper is organized as follows. Section presents the formulation of the CSD dynamic model by means of the PBE including growth, nucleation, and agglomeration kinetics. The MOC and the discretization scheme for the CSD are described respectively in sections and 4. In section the numerical solution of the PBE with the MOC and the proposed discretization scheme is addressed. This section also explains step by step how to handle the complexity due to the agglomeration term. The proposed algorithm is tested on the simulation of a seeded batch flash cooling crystallization in section . Conclusions are given in section .

The Population Balance Equation and Crystallization Kinetics

The most detailed model for describing the evolution of the CSD is the population balance equation (PBE) (eq ). The PBE describes the conservation of crystal identities in the crystallizer. The crystals can be distributed in time and space. The spatial domain may include an external (the position of the entity in the crystallizer) and an internal (the characteristic size of the identity) coordinates. In batch industrial crystallization applications, the crystallizer is often assumed to be well mixed, and the CSD is not a function of the external coordinate, but a function of the internal coordinate x and the time instant t:In (eq ) stands for the number density function which defines the number of the crystals per unit of crystallization volume and with size ranging from x to x + dx. is the growth rate, is the nucleation rate, B(x, t) and D(x, t) are the birth and death rates due to agglomeration and breakage phenomena. V is the crystallization volume, x0 is the dimension of the crystal nuclei. Often, the chosen internal coordinate is the crystal length L. Accordingly, the number density function has units [#/μm m3].

The crystallization is commonly operated in seeded, batch, or semibatch mode where the introduction of crystal seeds prevents undesirable primary nucleation phenomena that negatively influence the CSD. The secondary nucleation, generally due to attrition generated by the crystal–crystal and crystal–crystallizer collisions, produces small crystals of length L0. The agglomeration mechanism which involves the cohesion of parent crystals is often observed in industrial setups, while the breakage is less common at sufficiently low agitation speeds. Finally, the changes of the number density due to volume variation are usually negligible compared to the other mechanisms. Under the above-mentioned conditions the PBE (eq ) reduces towith left boundary condition (eq ) and initial condition (eq )where n(L,t) denotes the number density function in [#/μm m3], G denotes the growth rate in [μm/s] and B0 is the nucleation rate in [#/s m3]; D and B denote the death and birth rates due to agglomeration only. L0 is the dimension of the crystal nuclei. nseeds denotes the initial CSD of the seeds in terms of number density function.

The PBE is coupled with the mass and energy balances in the liquid phase of the crystallization system by means of the kinetics G, B, D, B0 which are functions of time with the supersaturation. Supersaturation, on the other hand, depends on temperature, and concentration of the solute. The crystallization kinetics considered in this work are listed in (eqs ). The size independent power law kinetics for crystal growth G (eq ) is widely used in crystallization modeling[12,16] because of its simplicity. The secondary nucleation B0 (4b) is modeled through the Evans kinetics[29] with only crystal-impeller collisions considered. The birth B and death D functions due to agglomeration phenomena are modeled according to Hounslow et al.,[18] as shown in (eq ) and (eq ). Note that the modeling of the agglomeration phenomena is a source of nonlinearities for the system. The agglomeration kernel β is considered size-independent and calculated according to the empirical expression (eq ).In (eqs ) k, k, Lmin, g, g, and a are kinetic parameters to be estimated. C and Csat(T) are the solute concentration and the solute concentration at saturation, respectively, N, N, and ε are the impeller parameters: flow and power numbers and specific power input, respectively.

The models for crystals agglomeration in length basis (eq ) and (eq ) were originally derived by Hounslow et al.[18] by converting the mechanism expressed in volume basis[4] (eq ) and (eq ) to a length-based form. In (eqs ) the functionality with time is omitted.where the prime is employed to emphasize that the variables and parameters are expressed in volume basis rather than length basis. The agglomeration kernel β′(v – ϵ, ϵ) is a measure of the number of collisions between particles of volume v and ϵ that successfully merge to produce a particle of volume v + ϵ. One can note that the solution of the agglomeration terms is more challenging when expressed in length basis compared with the volume basis because the former includes a nonlinear functionality of the internal coordinate.

For the sake of clarity, it must be pointed out that the kinetic laws (eq ), (eq ), and (eq ) reported in this section do not influence the derivation of the proposed methodology. Hence, they can be replaced with different ones on the condition that length-independent growth rate laws are chosen.

Method of Characteristics

As it is pointed out in the Introduction, the crystallization model (eq ) in length basis and with kinetics (eqs ) cannot be solved analytically without making certain assumptions. Thus, we have chosen the numerical solution with discretization methods. Among the various techniques,[12,13] the one based on the method of characteristic (MOC) is very appealing because of its simplicity and accuracy.

The MOC is a mathematical procedure that identifies characteristic curves. On these curves, the partial differential equation (PDE) becomes an ordinary differential equation (ODE). In the particular case of the population balance equation (PBE) at hand, this method eliminates the convection term through the transformation of (eq ) in[22]The system (eq ) is defined by the boundary condition (eq ) and the initial condition (eq ) inherited from the PDE (eq ):plus the extra initial condition for the internal coordinates:In (eqs ), the integral of the first equationcorresponds to the characteristic curve along on which the ODE for n is solved. In other words, the set of differential equations (eqs ) represents the evolution of the number density as seen by an observer moving with the growth velocity.

The use of the MOC is particularly convenient in the case of size-independent crystal growth, since in this case the system (eq ) simplifies as (eq ).which is accompanied by the boundary condition:and the initial conditions (eq ) and (eq )

Discretization of the Distributed Domain of the Internal Coordinate L and the Distributed Variable n

To solve the PBE with the MOC (eq ) in a simulation environment, the internal coordinate L has to be discretized along a certain number of grid points at every time interval, as shown in Figure . In the following we refer each point of the grid with the expression “length class” L. The crystal classes are collected in the vector variable Λ = [L0, L1, L2, ..., Lmax]T with dimension dim(Λ) = γ. The number of classes γ has to be tuned as a trade-off between accuracy of the numerical solution and computational time, and it may increase at every iteration as it will be explained in subsection . According to the discretized grid, the variable number density function n(L) distributed in the continuous domain of characteristic length L is treated as a vector N = [n(L0), n(L1), n(L2), ..., n(Lmax) ]T (or in short N = [n0, n1, n2, ..., nmax]T) associated with the vector of the classes Λ, as shown in Figure . Thus, dim(N) = dim(Λ) = γ.

Figure 1

Discrete crystal length and size distribution.

On the basis of this notation, the number of crystals M0, per unit volume in the size range [L, L] is computable with good approximation using the trapezoidal rule, provided that L – L is sufficiently small, that is to say:

Numerical Solution of the PBE by Using the MOC

To understand the solution of the PBE by using the MOC, let us decompose the general problem into three subproblems. Let us analyze the behavior of a system accounting for (i) only crystal growth, (ii) crystal growth and nucleation, and (iii) only agglomeration. The subsequent step will be the combination of the findings to simulate a crystallization process accounting for the three kinetics, and the improvement of the computational efficiency of the proposed algorithm.

Crystal Growth Only

In this section, the use of the MOC is addressed for the case of only crystal growth. In this case the system (eq ) reduces tofrom which it is easy to conclude that the number density maintains the shape of the initial seeds nseeds(L), i = 1, ..., γ while the length axis travels at a velocity equal to the growth rate, or in other words:every point of the initial profile is displaced in time t by a quantity Gt in the positive direction of L. The number of classes γ remains constant during the simulation duration. Practically speaking, an iterative numerical scheme to solve (eq ) is given bywhere Nsamp is the number of iterations needed to simulate the batch run, and t – tin is the sampling interval. tin (or tfin) is the beginning (or end) of the time interval. In this scheme one can use chemical and physical properties, and kinetics calculated at the beginning of the time interval tin (e.g., G = Gin) to obtain an approximation of the integral of the ODEs (eq a), provided that the time interval is small enough. If numerical instability or inaccurate solution of the ODEs are experienced, we recommend to shorten the sampling interval or to choose a different discretization scheme taking into account a balance between computational complexity (i.e., number of functions evaluations) and online implementation capability.

Crystal Growth and Nucleation

In the presence of crystal growth and nucleation the discretized model of the system within the MOC becomesPractically speaking, in this case the L axis is shifted in the positive direction due to the growth phenomena. In the mean time, the smallest modeled crystal class L1 becomes larger than the characteristic dimension of the nuclei L0. Thus, the distribution of seeded crystals does not undergo changes in its shape, and additionally, a number of crystals equal to B0/G is continuously added. From a numerical point of view, one can analyze the phenomena during a single sampling time: at every sampling time the L axis is shifted due to growth, and a new crystal class (or in other words a new element in the vector Λ) has to be added to accommodate the nuclei. Accordingly, the dimension of the grid γ increases with time, and the grid may become nonhomogeneously spaced. The behavior of the L axis can be visualized in Figure .

Figure 2

Schematic evolution of the length classes due to growth and nucleation.

In parallel, the number density of the crystals existing at time tin is conserved. In addition, a number of nuclei equal to B0/G is generated by the nucleation phenomena. A corresponding numerical scheme can beIt must be pointed out that the necessity of incorporating to the mesh grid a new length class L0 at every time interval makes the grid nonuniform when the growth rate is not constant with time. Indeed, the spacing between the first class L1 and its adjacent class L2 at any time tfin is L2, – L1, = L2, – L0 = G(tfin – tin).

The major drawback of dealing with the nucleation by means of MOC is that for long simulation durations the dimension of the vectors N and Λ becomes intractable, and the simulation time becomes comparable with (or longer than) the operation time, making the model difficult to use for any real time application. However, this problem can be easily overcome for batch processes by setting an appropriate sampling interval tfin – tin, which guarantees a good trade-off between speed and accuracy. In case this would not be enough, a regularization of the mesh grid by elimination of length classes in the region of the fines could be performed. According to the discretization scheme used in this paper, the elimination of the cells does not affect the shape of the distribution, while the number of the crystal is conserved by means of (eq ).

Agglomeration Only

The literature is limited when it concerns the simulation of the agglomeration term in length basis due to the inherit complexity of the solution of the PBE in this case. Thus, most of the models proposed are not suitable to describe industrial cases. In this section a numerical solution of the net agglomeration rate B – D compatible with the MOC is presented. The net agglomeration rate is handled by numerically solving the integral terms in (eq ) and (eq ) according to the discretization:The key step in the numerical solution of the integral is the evaluation of the number density n at the lengths (L3 – λ3)1/3, j = 1, ..., i, by using a linear extrapolation between values at computed nodes of the grid.

To test the accuracy of the proposed solution by comparing it with the available analytical solution in literature,[18] a pure agglomerating system is considered first. Under pure agglomeration, the system (eq ) becomesThe adopted numerical procedure is described in the following, highlighting the challenges and the strategies used to overcome them.

Let us recall the discretization scheme adopted as shown in Figure . The mesh grid is described by the vector Λ = [L0, L1, L2, ..., Lmax]T of crystal lengths. The number density function is described by the vector N = [n0, n1, n2, ..., nmax]T. The CSD is univocally defined by the pair [Λ, N].

∀L ∈ Λ:

1. The evaluation of (eq ) involves the computation of the elements of the vector Λ*and γ* is the dimension of the subset of crystal classesThis vector Λ* provides a set of lengths for which the number density has to be evaluated. The drawback is that, according to the proposed discretization, the number density n(Λ*) is usually not available because the crystal classes described by Λ* may not belong to the set Λ. In the example below, this concept is shown by using a vector of crystal lengths Λ = [Λ1, ···, Λγ]T, Λ ∈ Rγ and γ = 11 with elements defined as Λ1 = 0; Λγ = 100; Λ = Λ+10. Λ (or Λ) denotes the ith (or (i+1)th) element of the vector Λ. The corresponding vector Λ* computed at L = 100 is Λ100* = [100; 99.97; 99.73; 99.09; 97.82; 95.65; 92.21;86.13; 78.73; 64.71;0]. None of the elements of Λ100* but Λ100,1* = 100 and Λ100,11* = 0 belong to Λ (see Figure for further clarifications).

Figure 3

Comparison between the elements of the vectors Λ– and Λ* with a numerical example.

2.The number density n(Λ*) is calculated by linear extrapolation involving two adjacent length classes L, L such that L < (L3 – L3)1/3 < L, according to the simple formula:

3. Let us now recall the integral term in (eq )This integral is evaluated by computing the area under the curve defined by the matrix of points [Λ–, F] withby means of a trapezoidal rule:F is the vector of the integrand function evaluated at discrete points.

4. The birth rate by agglomeration is then given by

5. The evaluation of the death rate D(L) is much easier because it involves the computation of the first moment M0 of the CSD approximable withwhere M0, is the number of crystals per unit volume in the size range [L, L] computed with a trapezoidal rule according to (eq ). Thus,It must be pointed out that in the description of the numerical method adopted for the simulation of the agglomeration the variable time t is omitted, but the net agglomeration rate B(L) – D(L) is time dependent and must be recomputed at every sampling time of the simulation.

After the evaluation of the agglomeration rate B(L) – D(L), the system (eq ) must be integrated over time. For a small enough time interval (tfin – tin), the following approximation can be used:Efficiency and performance of the proposed solution scheme has been compared with the analytical solution of the agglomerating system proposed by Gelbard and Seinfield,[8] and recalled below (eqs –29)where τ = N0β t, and which is valid for the idealized case of an exponential CSDSince (eq ) is given in a volume basis, n′(v) must be converted in a length basis for comparison purposes with the proposed numerical solution by means of

Figure shows the time behavior of the CSD obtained with the analytical (solid line) and the numerical (dot line) solution. One can notice that the proposed numerical solution is in good agreement with the analytical one. Small deviations are observable in the region of the first half of the size domain, while no deviations are present in the region of the larger lengths.

Figure 4

Comparison between the analytical solution (solid curves) in length basis (eqs –29) and the numerical solution (dot line). The initial CSD (eq ) has parameters N0 = 0.285 and v0 = 3.50. The agglomeration kernel is taken constant at the value β = 1 during the five iterations.

Table shows the conservation of the third moment M3 of the CSD which is important to verify the accuracy of the algorithm, since M3 does not change under pure agglomeration kinetics. In particular, Table reports the third moment of the analytical solution and the numerical solution, and its loss due to the numerical computation. One can notice that the proposed numerical scheme is accurate since M3 is very well preserved. It has been noticed that the numerical errors are more abundant especially when the crystals size is distributed around 1, while they become negligible for crystal sizes preferentially distributed at L > 1. Therefore, the authors discourage the use of scaled length domains while adopting the proposed algorithm.

Table 1

Conservation of the III Moment of the CSD

	M3	M3
time	analytical solution	numerical solution	loss at every iteration
0	1.9099	1.9099
1	1.9099	1.8779	1.66%
2	1.9099	1.8533	1.31%
3	1.9099	1.8333	1.08%
4	1.9099	1.8166	0.91%
5	1.9099	1.8023	0.78%

To demonstrate this, another simulation is carried out, with seeds lognormally distributed around the length 74 μm and with a spread of 1.21. The behavior of the CSD due to pure agglomeration in this case is shown in Figure . Results are obtained applying (i) the proposed algorithm accompanied by a mesh grid for the crystal length uniformly distributed covering 300 μm with 300 crystal classes; (ii) agglomeration kernel-taken constant at the value of 1.508 × 10–12; (iii) 3600 s of crystallization; (iv) sampling time of 60 s. Under these conditions the simulation takes 25.1 s. The computation of the third moment reveals that in this latter case the volume of the crystals is almost perfectly conserved, with a total mismatch between M3 at 0 s and its value at 3600 s equal to 0.006%. By virtue of its accuracy and the limited computational effort, the proposed numerical scheme seems adequate for real time simulation for monitoring purposes. Moreover, it must be pointed out that the proposed method is very flexible since it is independent of the grid, which can be even inhomogeneous, and thus perfectly fits with the numerical technique based on MOC adopted for growth and nucleation simulation in this paper.

Figure 5

One hour agglomeration of a crystal population with initial log-normal distribution around 74 μm with a spread of 1.21.

Crystal Growth, Nucleation, and Agglomeration

When the combination of growth, nucleation, and agglomeration has to be simulated with the MOC approach, one has to consider the system (eq ) together with the algebraic equation (eq ) accounting for nucleation at the L0 length, and the initial conditions (eq ) and (eq ). This system is recalled below:Since all the phenomena previously discussed are involved here, a suitable combination of the schemes proposed in previous subsections must be implemented. For example, we suggest the following scheme at every sampling time,Here (eq ) and (eq ) upgrade the mesh grid due to nucleation and growth. In (eq ) the number of particles per unit volume n* is due to the nucleation and growth phenomena according toOne can note that this numerical scheme comes from the combination of (eqs ) and (eq ).

In the following sections, the proposed numerical scheme is applied to simulate a crystallization process achieved by means of the flash cooling technology.

Results and Discussion: Seeded Flash Cooling Crystallization

The method of characteristics (MOC) has been applied to simulate the behavior of a seeded batch flash cooling crystallization of a chemical species in water. In the flash cooling crystallization the driving force (the supersaturation) is achieved by the combination of (i) evaporation of the solvent and (ii) cooling of the solution, obtained by means of vapor extraction. The crystals dynamics are assumed to be determined by crystal growth, secondary nucleation due to attrition, and agglomeration, and the corresponding kinetics manifest a dependency with other states of the system, namely temperature and solute concentration.

The model of the process consists of material and energy balances for the liquid and solid phases. The dynamics of the temperature T (eq ), the concentration C (eq ), and the volume V (eq ) are described by ODEs, while the particulate feature of the solid product is modeled with the PBE.[5] Under the assumptions of perfect mixing, size-independent crystal growth rate, absence of crystals and solute in the vapor flow, dilute solution, and negligible effect of the volume on the dynamics of the CSD, the crystallizer model (eq ) is provided as follows.whereThe meaning of the symbols is presented in the nomenclature section. The values of the kinetic and impeller parameters, and crystal properties are reported in Table . The set of crystallization kinetics (eqs 4) has been employed to calculate G, B0, B, and D. Kinetics parameters are chosen to give a fast and growth driven crystallization with moderate nucleation and agglomeration. The flow rate of vapor Fw is an input of the system chosen to guarantee a low and almost constant supersaturation level over the batch run.

Table 2

Kinetic and Impeller Parameters

parameter	value	unit
kg	4.64 × 10–5	m/s
kci	exp(12.54)	#/m3 s
Lmin	100	μm
gg, gn	1
a	3.016 × 10–12
NQ	1.6
NP	2
ε	2.1	m2/s3
ρc	1600	kg/m3

The applied vapor profile generates a variation (i) of 22.5 K between the initial and final temperature, (ii) of 141.9 kg/m3 between the initial and final concentration and (iii) of 0.3 m3 between the initial and final volume, during 6720 s of batch run.

We assume that the seeds are lognormally distributed with location parameter μ = 74 μm and spread of σ = 1.6. The initial mesh grid is uniformly distributed and consists of 300 lengths with length interval of 2 μm. The nuclei have characteristic length L0 = 0.1 μm. The behavior of the crystallization system is simulated in MATLAB R2016a. The ODEs are solved with ODE45 within a time interval of 60 s. At the end of this time interval, the crystallization kinetics are provided and used to upgrade the CSD based on the numerical schemes (eqs 31–32) previously discussed. The value of the calculated second moment M2 is fed back to the ODEs for the next iteration. It must be pointed out that simulations conducted with a shorter time interval than 60 s return the very same dynamic behavior, but take a much longer simulation time. In the following, simulation results of the CSD dynamics due to (i) crystal growth, (ii) crystal growth and nucleation, and (iii) crystal growth, nucleation, and agglomeration are presented and discussed.

CSD Time Behavior under Crystal Growth Only

For the case of crystallization driven by the growth kinetics alone, the behavior of the number density n during the batch run is shown in Figure a. In this graph the CSD is plotted every 10 min. One can appreciate that the rate of shift of the CSD curve is decreasing during the batch run according to the growth rate which is depicted in Figure b. The latter depends in turn on the temperature and the concentration in the crystallizer.

Figure 6

(a) CSD in terms of number density n for a pure growth crystallization system at every 10 min of operation. (b) Evolution of the growth rate during the batch run under pure crystal growth.

For the sake of completeness, one can be interested in observing the behavior of the CSD in the more intuitive volume fraction function vf [m3/μm m3] which is presented in Figure . Compared with the number density n, the volume fraction vf retains the information on the total volume of the crystal phase according to vf = n kL3. For this reason the area under vf is increasing with time.

Figure 7

CSD in terms of volume fraction vf for a pure growth crystallization system at every 10 min of operation.

CSD Dynamics under Growth and Nucleation

The presence of the nucleation term implies that crystal nuclei of length L0 are produced. Thus, the number density evolves with a tale on its left side. The growth of the seeds is negligibly perturbed by the nucleation phenomena, as shown in Figure a. It must be pointed out that the nucleation rate (depicted in Figure b) is taken small, and consequently the growth rate is not significantly different than the one of Figure b, obtained for the case of growth only. If the nucleation rate were taken higher, the growth rate would have been lower because nucleation and agglomeration rates compete for the same driving force (the supersaturation). Due to the small nucleation rate, which generates a relative small left-tale for the number density n, the plot of the volume fraction vf is similar to the one obtained for the case of growth only, and for the sake of conciseness it is not repeated here. However, if the nucleation rate is larger, the volume fraction vf would develop a second mode in the region of the fines. The simulation is obtained with 60 s of sampling interval and it is as accurate as the one obtained using a sampling time of one second, but faster.

Figure 8

(a) CSD in terms of number density n for a crystallization system undergoing crystal growth and moderate nucleation, at every 10 min of operation. (b) Evolution of the nucleation rate due to crystal-impeller attrition during the batch run.

Growth, Nucleation, and Agglomeration

When all the crystallization phenomena are taken into account, the shape of the CSD in terms of number density (see Figure a) becomes very interesting. The average dimension of the crystals becomes considerably larger due to agglomeration, and the initial log-normal shape is not retained.

Figure 9

(a) CSD in terms of number density n for a crystallization system undergoing crystal growth, moderate nucleation, and attrition, at every 10 min of operation. (b) Net agglomeration rate B – D as a function of the crystal length, at the beginning of the batch (blue), after 5 (yellow), and after 10 (red) minutes of operation.

Figure b shows the net agglomeraton rate at the initial time t0 (blue line) and after 5 (yellow line) and 10 (red line) minutes. One can notice that (i) at the initial time crystals of length smaller than 60 μm are consumed to generate larger crystals, (ii) at time 5 and 10 min, the maximum consumption of crystals due to agglomeration has moved toward larger and larger crystal lengths, which is possible because large length classes are continuously generated to the detriment of crystals belonging to the finest classes which are consumed, and (iii) the agglomeration also involves fine crystals produced by nucleation.

Figure shows the effect of the phenomena on the volume fraction vf. One can notice that the volume fraction evolves three modes that can be associated with (from the right to the left) (i) agglomeration and growth of the initial seeds, (ii) growth of seeds, and (iii) nucleation and agglomeration of the nuclei.

Figure 10

CSD evolution in terms of volume fraction vf for a crystallization system undergoing crystal growth, moderate nucleation, and attrition, at every 10 min of operation.

Finally it must be pointed out that the simulation of 6720 s (112 min) of batch crystallization accounting for growth, nucleation, and agglomeration takes only 62 s, which makes the proposed numerical scheme very attractive for online use for real time monitoring purposes.

Conclusions

In this paper the problem of simulating seeded batch crystallization systems accounting for growth, nucleation, and agglomeration for online model-based estimation technologies has been addressed. The crystal size distribution (CSD) evolution has been modeled within the framework of the population balance equation (PBE). The method of characteristics (MOC) has been chosen for the discretization of the PBE. This method, coupled with a numerical integration with the upgrade of the CSD every 1 min, offers accurate solution and fast computational time. The proposed simulation model is thus an appropriate estimation model to be used within monitoring tools based on state estimators driven by plant measurements.