1. Introduction
Agglomeration and layering growth in fluidized beds are efficient particle formulation and size enlargement processes for solid materials, due to good particle mixing and high mass and energy transfer rates. The possibility to influence product particle properties such as size and porosity directly by a variety of different process variables [
1] and create a high-quality product is important from an industrial point of view. Therefore, these processes are applied widely for the production of agricultural, chemical, pharmaceutical and food products such as fertilizers, milk powder and fine chemicals. In the layering growth process a solution is sprayed continuously on the surface of the fluidized particles. After the solvent has dried a new solid layer remains whereby the particle size increases slowly during the course of the process and an onion-like structure is formed as presented schematically in
Figure 1. In contrast to this the size of the particles increases rapidly due to agglomeration of particles. When two wet particles collide in the turbulent fluidization air stream, there is a chance of aggregation due to viscous and capillary forces acting between the wet surfaces [
2]. After the drying of the solution a solid connection between the two particles is formed as pictured in
Figure 2. Layering growth is usually applied for the solidification of wet materials while agglomeration is desired if the size of particles has to be increased quickly, however depending on the process configuration and setting of process parameters both subprocesses can occur simultaneously [
2,
3].
For large scale production, these processes are operated continuously and often equipped with a sieve-mill cycle which serves for the separation, grinding and recycling of oversized particles. This configuration is more economic due to the reuse of the oversized particles. However, the recycle may also introduce instability in the form of self-sustained oscillations. These have been reported for layering growth in simulation studies of well-mixed process models [
4,
5] and models with zone-formation [
6] and confirmed experimentally [
5,
7,
8]. In all of the latter contributions porous sodium benzoate particles were produced over the course of 25 and 36 h respectively resulting in undamped oscillations of the mean particle size between the milled particle size and a maximal value. From a practical point of view such oscillations have to be avoided for a stable process operation by choosing proper process parameters or by implementing stabilizing process control methods [
9,
10,
11]. The occurrence of oscillations is mainly related to the size of the milled particles (or the milling grade) and also seems to be a result of the control of the total bed mass [
12] which, however, is indispensable from the operational point of view. In fluidized bed configurations where agglomeration is the dominant size enlargement mechanism stable oscillations have not been observed to the best knowledge of the authors. The question regarding processes where both size enlargement mechanisms occur has been investigated by the authors in the conference contribution Otto et al. [
13] by means of numerical bifurcation analysis. It has been shown in different simulation scenarios, that generally, shifting size enlargement from pure layering growth to agglomeration, tends to dampen oscillations and stabilize the whole process. The aim of the present contribution is to give a more detailed description of the underlying process model, to extend the previous results to a larger set of bifurcation parameters and give explanations of the observed phenomena.
The second and third section are concerned with the description of the population balance process model and the parameter bifurcation algorithm, respectively. The contribution concludes with the presentation and discussion of bifurcation plots for various parameters in the fourth chapter and a short summary in section five.
2. Process Model
In this contribution a plant configuration with sieve-mill-recycle is investigated where particles are added to and withdrawn continuously from the fluidization chamber. The withdrawn particles are separated into three fractions by two sieves. Particles meeting the desired product size specification are withdrawn while oversized particles are milled and re-fed to the process chamber together with the undersized fraction. A schematic representation of the plant with an additional external particle feed is presented schematically in
Figure 3.
A well-established framework for the mathematical description of agglomeration and layering growth processes is population balance modeling [
14]. Here the evolution of a particle population, distributed with respect to some internal or external properties
x, such as size, shape or spatial coordinates, is described by a nonlinear partial differential equation, the so-called population balance equation (PBE). The current state of the population is given by the number density function (NDF)
, describing the number of particles at time
t in the infinitesimal interval
of the internal coordinates. For well-mixed systems such as fluidized beds it is reasonable to assume that there are no spatial gradients in the bed, therefore spatial coordinates can be neglected. An important internal granule property is the characteristic size, since many other important properties such as the angle of repose or the porosity are correlated to it. In fluidized bed processes it is convenient to assume that all particles are spherical, thus their respective size can be characterized by either the particle diameter
L or the according volume
v, given by
While layering growth is most commonly described in terms of the particle diameter [
15,
16,
17], the agglomeration term is predominantly given with respect to the volume coordinate [
14,
18]. In the present contribution we choose to transform the layering growth term and model the complete PBE in volume coordinates, i.e.,
and
. The solution of the agglomeration term requires the application of efficient numerical techniques guaranteeing conservation of total mass and correct computation of particle number. For those approaches the description of the particle density with respect to particle volume is convenient, whereas the layering growth term can be discretized using a first-order upwind scheme regardless of the choice of the internal coordinate.
The PBE describing the evolution of
is obtained by balancing all the particles streams into and out of the fluidization chamber according to
Figure 3 and adding sink and source terms reflecting layering growth and agglomeration, yielding
In the subsequent paragraphs the right-hand side terms of Equation (
2) are described in detail in accordance with Radichkov et al. [
4], Heinrich et al. [
15] and Neugebauer et al. [
5] as well as Otto et al. [
13].
The diameter-based layering growth is described by the advection type equation
where
with
describes the diameter-independent growth rate and is given by
which can be found e.g., in Heinrich et al. [
15]. Here,
describes the injection rate of binder mass,
the binder density and
the total particle surface defined by
It is assumed that the binder solution is distributed uniformly on all particles and that the solvent evaporates completely. In order to transform Equation (
4) to a volume-based formulation
the volumetric growth rate
G has to be computed from
for which the relationships
as well as Equation (
1) and (
4) are utilized, resulting in
Notably, the growth rate is inversely proportional to the total particle surface, resulting in slow particle growth for beds with high surface.
The particle size enlargement by agglomeration is modeled using the following integral terms
where the first term on the right-hand side describes the ‘birth’ of new agglomerates due to coalescence and the second term describes the ‘death’ of the two parent particles [
14,
18]. The kinetic rate of this second order process is given by the agglomeration kernel
describing the rate of successful agglomeration events depending on the volume of the involved particles. Due to the complex interaction of different mechanical and thermodynamical micro-processes leading to the aggregation of two agglomerates, the function
depends on various process, plant and material parameters. For example, the viscosity of the binder affects the inter-particle forces during particle collision and the temperature of the fluidization medium affects the drying behavior of the binding agent on the particle surface and thereby the formation of solid bridges between the primary particles [
1,
2,
19]. Hence, modeling of
is a challenging task with different approaches of varying complexity to be found in the literature [
2,
20,
21,
22,
23]. In order to investigate the influence of different agglomeration kinetics on the stability and at the same time keep the model simple, the rather general kernel function
is applied in this contribution. With the free parameters
a,
b and
, the latter also called agglomeration efficiency, different shapes of the kernel function and hence different types of agglomeration kinetics can be approximated [
24]. In the nominal case
the kernel function becomes
, corresponding to a volume-independent agglomeration rate.
The feed of external nuclei, also called primary particles is model by
under the assumption the total number of feed particles
is distributed normally around the nominal volume
.
It is assumed that particles are withdrawn from the fluidization chamber without internal separation, i.e.,
where
K denotes the withdrawal rate. The computation of
K will be explained at the end of this section. The withdrawn particles are sieved externally using two sieves with opening diameters
in order to separate particle within the desired product size range
from over- and undersized ones. The non-ideal separation functions characterizing the sieves are modeled by cumulative Gaussian functions
where
determines the selectivity of separation. By introducing Equation (
1), the volume-dependent separation functions
and
are obtained. The fines particle fraction, given by
is fed back directly into the fluidization chamber without time delay for further size enlargement. The product fraction
is removed from the process, while the oversized fraction
is milled into fine particles to serve as new nuclei. For simplicity it is assumed that the milled particles are distributed according to the Gaussian mill function
with milling grade
and variation
. From the number conservation the volume-based mill distribution is obtained as
The number density distribution of milled particles is obtained by equating the mass of oversized and milled particles [
5,
15]:
From a practical point of view it is crucial to keep the total particle bed mass within certain boundaries in order to maintain the fluidization state. Hence, the system is described as a differential-algebraic system with the additional conservation equation
where
is the desired bed mass and
is the particle density. The integral term on the right-hand side describes the first moment of the particle size distribution and corresponds to the total particle volume in the fluidization chamber. Differentiating Equation (
20) and inserting the PBE (
2) yields the ideal withdrawal rate
compensating for the increase in bed mass by binder and external nuclei addition.
The model presented above has been validated without the agglomeration term and some small adjustments in e.g., Neugebauer et al. [
5], the agglomeration process model has been validated experimentally but without recycle in Otto et al. [
25].
3. Bifurcation Analysis
This section is concerned with the bifurcation analysis extending the results from Radichkov et al. [
4], Dreyschultze et al. [
6] and Neugebauer et al. [
5] for the process with agglomeration. For the bifurcation analysis the DAE system defined by Equations (
2) and (
20) has to be solved numerically, typically achieved by application of a moment method [
26] or the discretization of the volume coordinate [
27,
28]. In this contribution we choose the finite volume scheme presented in Singh et al. [
28] which achieves high accuracy of the NDF compared to the moment methods and is straightforward to implement as well as numerically efficient compared to more complex discretization methods. The resulting
N-dimensional system of ordinary differential equations is implemented in Matlab and solved using the build-in function
ode15s for stiff DAE systems.
Since the aim of this contribution is to investigate the stability behavior of the process with layering growth and agglomeration, relevant bifurcation parameters have to be identified. Previous contributions have established the mean size of milled particles
to be the most important variable regarding the occurrence of limit cycles in the process, therefore we adopt it as the first bifurcation parameter. In order to investigate the influence of particle agglomeration it is recognized that the rate at which aggregation events between two granules occur is directly determined by the agglomeration efficiency
. The influence of other plant and process parameters is evaluated by investigating three samples in each case. The two-dimensional bifurcation analysis is reduced to a series of one-dimensional problems by discretizing
logarithmically in the interval
. For every value of
a parameter continuation [
29] w.r.t.
in the interval
is conducted. The parameter continuation starts in the stable region with
, where the initial steady state NDF is obtained by solving the process model on the time interval
. Afterwards
is decreased by a small step
and a prediction for the corresponding steady state is given using a tangent predictor [
29]. The prediction is corrected by using it as initial guess for minimizing the sum of squares of the time derivatives by a Levenberg–Marquardt algorithm implemented in the Matlab function
lsqnonlin. The procedure, where the NDF function is constrained to positive values, is repeated iteratively until
is reached. In order to determine the stability of the steady state, the process model is linearized and eigenvalues are computed. In the case of at least one real part greater than zero, the corresponding steady state NDF is unstable and the existence of a stable limit cycle is tested by time integration. The bifurcation parameters are presented in
Table 1.
4. Results
As a basis for subsequent investigations, the stability analysis is carried out for the nominal case, where the plant and layering growth related parameters presented in
Table 2 correspond to those given in Radichkov et al. [
4], Dreyschultze et al. [
6] and Bück et al. [
30].
The stability plot for the nominal parameter case is presented in
Figure 4. The mill grade interval is chosen as
. Larger particles are already part of the product size interval and smaller particles are undesired since in a practical fluidized bed they might be removed form the chamber by the exhaust air. In the area colored white, one stable steady state solution exists while in the black area the PDE solution converges to a stable limit cycle around an unstable equilibrium. For
, 1 × 10
−14, the amount of agglomeration events is negligibly small and the bifurcation results correspond to the results presented for the layering growth only case in Radichkov et al. [
4] and Dreyschultze et al. [
6], where a zone of instability is limited by two supercritical Hopf-bifurcation points located at
and
. For increasing values of
both bifurcation points are shifted towards smaller mill grades until the left one reaches its minimum mill grade value at
. Further increasing
, the two bifurcation points are moving towards each other and collapse at
. For even larger values of
no self-sustained oscillations are observed.
In the following an explanation for the vanishing of the self-sustained oscillations with increasing agglomeration rates is developed. To this end, the oscillation of the NDF in the layering growth only case is described in three steps in accordance with Radichkov et al. [
4]:
Beginning from a narrow particle size distribution, presented in
Figure 5 (left), with a high fine particle mass fraction, the average particle diameter and the particle growth rate are small due to the high particle surface. Since the oversized mass fraction is small, a limited amount of particles is milled and only few new nuclei are formed. Therefore, the particle size distribution stays narrow.
Due to the particle growth the size distribution shifts towards larger particles and the average particle diameter rises. The fines mass fraction decreases while the product mass fraction increases (
Figure 5, middle). Due to the narrow shape of the particle size distribution a large amount of particles reaches product size in a short time interval. Since the total bed mass is conserved only some of them are withdrawn and the particles grow further into the oversized range.
The oversized particles are withdrawn and milled and a large amount of new nuclei is formed in a short time interval (
Figure 5, right). Since the total particle surface increases quickly the growth rate decreases accordingly and the narrow particle size distribution of the first phase is formed.
For the oscillation cycle described above it is crucial that the number density function stays narrow during the size enlargement, leading to two distinct “phases” where a uni-modal NDF either mainly grows or is mainly milled and fed back to form a new peak around the milling diameter. When is larger than the upper Hopf-bifurcation point, the oversized diameters and are too close for the emergence of distinct phases. For values of below the lower bifurcation point, the amount of new nuclei generated by the milling and the accordingly large total surface slow down the growth such that a large amount of particles can be withdrawn as product before reaching the oversized range. Thus, the formation of a second mode of nuclei is prevented and a stable equilibrium emerges.
In contrast to narrow particle size distributions, described above, wide particle size distributions can simultaneously grow in the small diameter region and be milled and fed back in the large diameter region leading to a flow equilibrium and therefore a stable distribution. It is well-known that the different mechanisms of layering growth and agglomeration promote narrow and wide size distributions, respectively. Due to the convection term describing the layering growth (
6) the particle size distribution largely maintains its shape and only shifts to the right. The agglomeration term (
9) however allows, depending on the specific agglomeration kernel, for the combination of two parent particles forming an agglomerate twice the volume quickly leading to a stretch of the initial number density function and therefore, in general, to a stabilization of the process dynamics. The larger the agglomeration efficiency is, the wider the NDF becomes. For
and
however, agglomeration leads to the emergence of oscillations. In the granulation only case, the high total surface of milled particles slows down the growth and enables a stable equilibrium, here particle aggregation quickly reduces the particle surface and allows quick growth even if the milled particles are small.
In [
11,
31] the dynamics of a rotary drum granulation setup with sieve-mill recycle and size-enlargement due to layering and agglomeration have been investigated. The process model is quite similar to the present one, however an additional spatial coordinate along the rotational axis is introduced. Although no detailed bifurcation analysis with respect to the agglomeration efficiency has been conducted, self-sustained oscillations have been observed for process configurations where the milled particle size (crusher gap) is below some critical value. The proposed explanations are in accordance with the results presented above.
In the following, the influence of additional process and plant parameters on the occurrence and stability of equilibria is investigated. By varying the Kapur kernel parameters
a and
b the shape of the agglomeration kernel is changed. The agglomeration kernels depicted in
Figure 6 show that an increase of the parameter
b between
and
results in preferential agglomeration of small particles by increasing their agglomeration rate. The resulting bifurcation plots are presented in
Figure 7, showing that the region of instability shrinks slightly. Since the variation of the parameter
b increases the total agglomeration rate over all volumes, the results support the general finding that increased agglomeration dampens oscillations.
The layering growth rate
G is directly proportional to the binder injection rate
. Therefore, increasing the latter results in faster particle growth and more particles in the oversized fraction, preventing the formation of the stable equilibrium for mill diameters smaller than
due to the mechanism described above described above. Decreasing
has the inverse effect as presented in
Figure 8. The lower growth rate allows for the withdrawal of most particles in the product region, in turn minimizing the oversized fraction such that no new mode of milled particles is created.
The feed of external nuclei creates new nuclei for agglomeration and layering. The main reason for the occurrence of self-sustained oscillations is that only few new nuclei are added in the “growing phase” and the NDF maintains its narrow shape. Therefore, reducing the number of new nuclei by decreasing the feed rate
promotes the occurrence of limits cycles and the expansion of the unstable region as presented in
Figure 9.
The quality of the product particles is mainly determined by their average diameter and its variation which can be improved by varying the sieving diameter if required. Hence, a change of the sieve diameter
around its nominal value
is investigated and presented in
Figure 10. It can be seen that, in general, larger values shift the unstable region towards larger mill diameters. Due to the reduced product size range, oscillations are possible for larger mill grades and the upper Hopf-bifurcation points are shifted to the right. On the other hand, the stable region left to the lower bifurcation points is also more pronounced. One reason for this might be the increased time span milled particles have before they reach product size allowing for the NDF to broaden by addition of external nuclei.