Aggregate Structure Evolution for Size-Dependent Aggregation by Means of Monte Carlo Simulations †

Aggregation during crystallization and precipitation processes often leads to complex-shaped particle aggregates. As an alternative to low-dimensional deterministic population balance models, where assumptions on the particle shape must be made, stochastic or so-called Monte Carlo methods can be employed. In previous work a hierarchical characterization of aggregates has been proposed (Briesen, AIChE J., 52, 2436-2446, 2006), which allows the use of different levels of detail for modelling the different rate processes as primary particle growth or particle aggregation. With that hierarchical characterization, the detailed geometry of aggregates becomes accessible for rate process modelling and product characterization. Here, this framework is extended to investigate size-dependent collision rates and aggregation ef ficiencies. The results show that the aggregate structures can be modelled by the interplay of shear rate and the growth rate at the particle necks in a mechanistic way. Future work will address the comparison with experimental data and alternative model formu-


Introduction
Aggregation during cr ystallization and precipitation processes often leads to complex-shaped particle aggregates.Even though the primary particle shape and aggregate structure affect downstream processing properties such as filterability, flowability or mechanical stability, the detailed shape is rarely considered.The actual shape is often so complex that a detailed modelling which accounts for the explicit shape is infeasible due to the associated computational effort.To avoid this effort, one-dimensional population balance formulations may be used (e.g. 1) ).There, the particles are typically characterized by one single characteristic size.However, the aggre-gates and their behaviour can hardly be described by only one characteristic length.If the dimensionality is extended, size and shape descriptors can be used as inner coordinates.But even then, a rigorous characterization of complex-shaped aggregates is limited.Additionally, the consistent formulation of rate expressions for multivariate formulations of the population balance equation becomes very tedious.As an alternative to deterministic population balance models, stochastic or so-called Monte Carlo methods have gained increasing popularity for the simulation of par ticulate processes [2][3][4][5][6] . Tere, the evolution of a set of representative, discrete particles is tracked.The interaction of the particles is realized by means of probabilistically selected events.Monte Carlo methods feature three principal advantages: Simplicity of implementation, capability of dealing with high-dimensional problems and ease of representing complex behaviour.Monte Carlo methods are highly promising as the dimensional restrictions of deterministic population balance modelling are a major obstacle to the mechanistic modelling of particulate processes.Despite these principal advantages, Monte Carlo methods are not yet commonly employed because of their usually high computational requirements.Nevertheless, the use of Monte Carlo methods seems to be a logical consequence especially in the light of unprecedented and ever increasing computational power.Stochastic methods have frequently been used to investigate the behaviour of fractal aggregates.Based on the pioneering work by Witten and Sander 7) as well as Meakin 8) , many studies along these lines have appeared in the literature in which the mechanisms are categorized in diffusion-limited, reaction-limited or ballistic cluster aggregation.All these mechanisms lead to different fractal dimensions that provide information on the openness of the aggregates.However, this fractal structure only becomes evident if a very large number of primary particles per aggregate is considered.Here, the number of primary particles per aggregate is too low to characterize the aggregates by a fractal dimension only.Therefore, the detailed structure of the aggregates must be reflected in the characterization to permit appropriate modelling.The implementation of Monte Carlo schemes for solving population dynamics problems is straightforward and far less restricted in the dimensionality of particle characterization.However, their accuracy is of course limited by the number of particles considered in the discrete set.As the particle number has to be large to achieve a satisfactor y accuracy, computation of the single steps in the Monte Carlo algorithm has to be as efficient as possible.Although the geometric operations involved to per form an aggregation event for complex-shaped particles are possible in principle, the consideration of a large set of these particles makes this fully detailed approach prohibitive.Instead, ways of reducing the complexity of the systems characterization need to be found.In previous work 9) , the author has presented a hierarchical particle characterization which uses different levels of detail for representing the different particle phenomena.For aggregation, a substitution system of seven point masses is constructed which preserves the main geometr y of the particles.The main features of this reduced characterization will be discussed briefly in the next section.With this type of characterization, it was possible to qualitatively model complex aggregate structures with only a moderate demand on computer power.As the focus of that contribution was on the hierarchical characterization, many assumptions have been introduced which lead to an unrealistically simple problem description.The present contribution extends the previous approach towards a more realistic modelling of the rate processes accounting for size and structure-dependent behaviour.To rigorously consider size dependency, the aggregation rate for each combination of particles in the Monte Carlo set would need to be computed.To avoid this large effort, the particle set is grouped into several classes for which identical collision behaviour is assumed.In the Monte Carlo scheme, initially two classes are probabilistically selected.A collision is executed for two particles randomly selected from the respective classes.Whether the collision is successful in forming a new aggregate is determined explicitly as a function of the system parameters (physical properties of the solid and the fluid, shear rate, linear growth rate).To determine the success of a certain collision, a dimensionless group suggested by Hounslow et al. 10) is used, which has proven to be very successful in correlating their data.

Hierarchical characterization
In this section, the main features and the construction of substitute systems for aggregation are introduced briefly.A hierarchical particle characterization allows selection of a complexity level which is just appropriate for a certain rate phenomenon.A reduced characterization can be used to perform aggregation, whereas the fully detailed characterization is still available.Thus, the detailed characterization can always be used to derive alternative reduced representations which are more suitable for rate events other than aggregation or also for a more detailed aggregation model.The discussion of the reduced characterization is separated into the characterization of the primary and the aggregate particles.For each of these systems, a fully detailed and a reduced characterization are introduced.Each of the primary particles is characterized by three dimensions.This allows arbitrary cuboids to be represented.Additionally, each primary particle within an aggregate structure is associated with a 4×4 matrix that specifies the relative position of each of the primary particles within the aggregates.The geometrical complexity of this characterization, however, is too high to permit evaluation of aggregation events in a reasonable computation time.To reduce this complexity, substitution systems comprising a set of point masses are introduced.In the substitution system, each primary or aggregate particle i is represented by seven point masses.For the primary particles, the construction of the set of point masses is trivial.If symmetry is exploited and conditions to guarantee the preservation of the moments of inertia with respect to the coordinate axes are imposed, the point mass system can be determined explicitly.The construction of a reduced representation by means of seven point masses for an arbitrar y aggregate structure is more demanding.First, a principal component analysis (PCA) is performed for the set of point masses comprising all point masses of the primar y particles.This operation yields the main orientation of the aggregate.The actual position and mass value for the point masses is constructed by preserving the moments of inertia with respect to the principal component directions.Additionally, the position of the point masses must reflect the extensions of the original aggregate.Fig. 1 shows the full and the reduced representation by means of point masses of a primary and an aggregate particle.In this work, only the position of the mass points is relevant.This position is represented by a matrix in which the rows correspond to the coordinates of the mass points in a coordinate system which has its origin in the centre of mass and whose axes are given by the principal components directions: Note that this reduced system, instead of the full characterization, is used to perform the actual aggregation event.The detailed information is nevertheless still available and can be used to assess the simulation results.Fur ther details of the par ticle characterization and its construction can be found in previous work 9) .During the simulation, the aggregates are sorted into size classes to study size-dependent behaviour.This is similar to a sectional method when discretizing a deterministic population balance equation 11,12) .The size of the aggregate is determined only on the basis of the substitute system.Consider a pseudo-ellipsoid spanned by the axes Li,1 to Li,6.As all the axes have a different length, the volume of the pseudoellipsoid comprises 8 different sections.The size of an aggregate is considered to be the diameter of a sphere with the equivalent volume to the pseudoellipsoid: Note that contrary to standard approaches, it is not the actual aggregate volume, comprised by the total volume of the primar y particles, that is employed.The structural properties of the aggregate can thus be considered in the evaluation of the collision rate (i.e.open structures with the same actual volume have a larger collision diameter).The number of size classes is fixed at 20 for all the simulations presented below.Subsequent intervals are enlarged by a factor of 1.07, respectively.The interval with the largest size values is chosen such that the aggregation efficiency is zero for any value in that inter val.The effective aggregation rate for two aggregates which belong to the size classes i and j defined by the interval (di, di+1] and (dj, dj+1], respectively, is given below (see equation ( 9)).

Monte Carlo scheme
The following steps constitute the basic Monte Carlo scheme: 1. Initialize system.
2. Determine the time between two aggregation events.3. Select two particles to be aggregated.Particles may be primary or aggregate particles.4. Perform aggregation event with the selected particles.5. Check termination criterion.If 'false', go back to step 2. Step1: The system is initialized by specifying an initial number of primar y particles.The initial particles are chosen to be cubes with a specified side length.Additionally, the initial solids volume fraction in the suspension is chosen.
Step 2: The inter-event time strongly depends on the cur- rent state of the system; namely, the size distribution of the current set of particles.To determine the interevent time, all possible effective aggregation rates between particles of two size classes need to be evaluated.The effective aggregation rate constant βij between particles from size class i and j is a product of the collision frequency Cij and the collision efficiency Ψij: For the collision rate, the classical expression for aggregation in a laminar shear flow as proposed by Smoluchowski 13) is used: where γ ˙. is the shear rate and di is the representative size of size class i.The aggregation efficiency is evaluated following the contribution by Hounslow et al. 10) .They concluded from micromechanical considerations that a collision is only effective if the aggregate has enough time to grow a stable neck.The most influential parameters are therefore the growth rate at the neck and the shear rate trying to disrupt the aggregate.They introduced a dimensionless group M which was effective in correlating their simulation and experimental results: where σ* is the material-dependent apparent yield stress and G is the linear growth rate.The kinematic viscosity is given by μ, and d is the diameter of a sphere from a mono-disperse system.Plotting the efficiency determined from their previous simulation studies showed that the effects of supersaturation, shear rate and particle diameter could be perfectly correlated with the single dimensionless group M. Hounslow et al. 19) only give a graphical representation of that correlation.To obtain a functional representation Ψ＝ feff(M), a sigmoid-type function, which is based on the graphical data and shows very good agreement, is used here.The relation used is plotted in Fig. 2.
One major assumption made in the derivation of this dimensionless group is that the system is considered to be mono-disperse.For the diameter d, Hounslow et al. 10) used a mean particle size of the whole distribution.Here, the explicit size dependence shall be considered.As the analytical derivation can only be performed for mono-disperse spheres, a pragmatic approach is chosen to at least capture the qualitative dependence for a size-dependent dimensionless group: Accordingly, the aggregation efficiency is also evaluated in a size-dependent way for each pair of possible sizes.The aggregation kernel for a particular combination of sizes di and dj is then given by As mentioned above, the particles are grouped into size classes.Hence, representative values for these classes need to be determined because equation ( 6) only provides the aggregation rate for two distinct sizes.Accordingly, an averaging of the rates over each size class interval is performed: A graphical representation of the aggregation rate kernel is given in Fig. 3 for a particular set of parameters.If the particles are very small their collision efficiency will be 1.Accordingly, the effective aggregation rate increases with increasing particle size because it increases the cross-sectional area and consequently the probability of a collision.For a further increase of the particle size, the effective rate decreases because of the strongly decreasing collision efficiency for large particles.The actual rate of effective collisions between two particles of the size classes i and j can then be evaluated by where Ni and Nj represent the volumetric number densities of the respective size classes.The aggregation kernel βi, j depends only on the particles' size classes.Thus, the term defined by equation ( 7) can be pre-processed before the actual Monte Carlo run.
From the rates of effective collisions, an average time between two aggregation events can be determined by ∆t = 1 where V is the current volume domain of the Monte Carlo scheme.
Step 3: The selection of a combination of size classes from which particles should be aggregated is governed by the effective collision rates for the respective size class pairs.Obviously, preference should be given to a combination with a high collision rate instead of a combination with a very low rate.For easier notation, consider the collision rates ri,j to be arranged in a vector rˆ with the indices given by k = (i−1)・nj + j : The index ksel for the selected pair of size classes (isel, jsel) is determined from the where RND is a random number on the interval [0,1].
Step 4: From each of the selected particle size classes, one particle is randomly sampled.The following discussion of step 4 will always refer to these two particles denoted by indices 1 and 2. The aggregation of the particles is realized with a procedure similar to that proposed earlier 9) .There, two mass points of the substitution system are randomly selected to coincide.Also, they are oriented such that they share a common axis.As selection of the mass points is completely arbitrary, the structural properties of the two aggregates are not properly reflected.Although a fully rigorous model would exceed the scope of this contribution, a strategy is developed which at least qualitatively reflects the expected structural evolution behaviour.To achieve this, each combination of mass points from the selected particles is associated with a probability of being part of the axis along which the new aggregate is formed.As the aforementioned size classes only reflect a spherical diameter, this averaged size is not suitable for representing the structural properties.Responsible for the collision rate and the collision efficiency is the cross-sectional area which is exposed to the shear flow.After collision, the aggregate tumbles in the shear flow.The orientation of maximum normal force must generally be expected to be different from the orientation at collision.Therefore, different representative values are used for the collision diameter di,coll and the efficiency diameter di,eff.For two selected mass points e and f out of the possible 6 outer mass points represented in equation ( 1), the effective collision rate can be evaluated: (11) The collision diameter is chosen as the circle diameter which has the same cross-sectional area as a pseudo-ellipse.The pseudo-ellipse is spanned by the four axes which do not fall on the same coordinate axis as the selected substitution mass.As the general formulation hinders readability, only the collision diameter for e = 1 is given: For other values of e, the collision diameter is computed analogously.The rationale behind this is that if the particles are orientated in a way that the crosssectional area of the pseudo-ellipse is exposed to the flow, the particles are likely to collide with the remaining axis.
For the efficiency diameter, the doubled value of the axis of the attaching mass point is chosen: The idea behind this choice is to prefer the formation of compact aggregates as the particle size reaches a value where the collision efficiency drops from 1.
The formation of aggregates at the mass point with a small axis should yield more compact aggregates.
Note that these assumptions only allow a qualitative analysis of the structural behaviour.For a quantitative representation, more rigorous examinations of the local micro-mechanical behaviour are needed.Also, it must be fundamentally investigated whether the micro-hydrodynamics of complex-shaped aggregates can be represented by equivalent sphere diameters in the first place.The presented choice nevertheless already allows a first mechanistic estimation of the qualitative behaviour.

Case Study
A closed aggregating system with constant shear and growth rate is investigated as a case study.Different shear rates and growth rates are studied with the presented algorithmic framework to assess particle size and structure.Note that the growth rate will only be used to reflect the growth of the neck of attaching particles.The particles themselves are not considered to grow.The variation of shear rate affects both the collision frequency and the collision efficiency, whereas the growth rate only affects the aggregation efficiency.The values used in the present study are shown in Table 1.The case numbers are used to identify certain combinations of shear and growth rate when reporting the results.For the simulations a nominal, initial particle number of 10000 is chosen.The side length of the initial cubes is 15μm.The initial solids volume fraction is set to 3.3E-5.For the apparent yield stress, the value suggested by Hounslow et al. 10) for calcium oxalate monohydrate is chosen: σ*= 150 GPa.The viscosity μ is set to 1E-3 Pa s.To obtain comparable results, a final average number of primar y particles per aggregate is specified as the termination criterion for the simulation.In other words, the simulation is stopped when the average number of primary particles is 10.This means that the total simulation time may vary strongly.On the other hand, the structural properties of the resulting particles ensemble can be compared on a common basis.

Results
Fig. 4 shows the total number concentration as a function of time for the 9 investigated cases.The number concentration decreases faster for increasing shear (case sequence: 1→4→7, 2→5→8, 3→6→9) and growth rate (case sequence: 1→2→3, 4→5→6, 7 →8→9).Except for the cases 4 and 7, all simulations reach the termination criterion where on average an aggregate comprises 10 primary particles.Cases 4 and 7 reach a state where the effective aggregation rates are zero for all particle combinations before fulfilment of the termination criterion.In other words, the shear rate is already too large to allow the formation of larger aggregates.As the time step between two successful aggregation events increases dramatically in the final phase of these cases, the time axis is cut off at 15000 s, whereas the final time for cases 4 and 7 is in the order of 10 12 s.As a qualifier for the morphology of the particles, the compactness of the particle is introduced.The compactness ω here is defined as the ratio of the actual aggregate volume Vagg and the volume of the cuboid spanned by the substitute system: For the primary particle, the compactness is 1 and decreases quickly to smaller values as the primary particle number in the aggregate increases.Fig. 5 shows a histogram to assess the final state of the particle ensemble.The aggregates are grouped according to their number of primar y particles and their  compactness.The low growth rate in case 1 hindered the appearance of larger aggregates.To reach the final average number of 10 primary particles per aggregate, all initial particles disappeared and formed aggregates with a narrow variation in size and compactness.In case 7, the shear rate was even too large to reach the final state.Similar to case 1, a narrow distribution of particle sizes and compactness is obtained.For the large growth rates, the aggregation efficiency remains large even for large aggregates.As the collision rate for the larger aggregates increases, many initial primary particles remain un-aggregated.
As the growth rate of the neck is no longer the limiting factor, large particles with up to 35 primary particles per aggregate are observed.As the shear rate increases from case 3 to 9, the aggregates are only formed up to smaller sizes compared with case 3.
Note that the segmentation shown in Fig. 5 is not used in the simulations.In a deterministic population balance framework, such a representation would probably also have been chosen to investigate size and structure evolution.In the case of the deterministic population balance approach, equations would have been needed to model the aggregation and its effect on structural changes.However, the information needed to mechanistically formulate this equation (i.e. the explicit structure the aggregates) would not be available.In the present approach, the structure is a result of the simulation and no structural assumptions must be made.The full morphological information is available at all times.Fig. 6 shows randomly selected particles for cases 1, 3, 6 and 9 at the final simulation time.As already obser ved from the histogram plots in Fig. 5, the particles are much smaller and more uniform for the small growth rates.For larger growth rates (cases 3 and 9), the final system comprises non-aggregated particles as well as large structures.For increasing shear rates, the particle size tends to be smaller.The largest aggregates observed in cases 1-9 are given in Fig. 7. Obviously, the maximum observed size increases with increasing growth rate and decreases with increasing shear rate.From case 3 in Fig. 6 and Fig. 7, it can already be anticipated that the resulting structures are more open.For the final particle state, the aggregation is not yet limited by the collision efficiency.All particles thus stick together regardless of their orientation.As the collision efficiency decreases, only compact aggregates are allowed to be formed.This is clearly observed for cases 1 and 3. Also, the particles for case 9 seem to be more compact.However, as the visual interpretation also depends on the orientation of the particles, the results may be misleading.To quantify the structural differences, the compactness is used.Fig. 8 shows the averaged compactness for all aggregates with exactly 15 primary particles.Only the comparison for the same number of primary particles is meaningful.Fig. 7 Aggregates with the largest number of primary particles at the end of the simulation for the cases 1-9.
Fig. 8 Average compactness for aggregates with 15 primary particles per aggregate.To ensure a certain statistical relevance, data is only shown if at least 100 aggregates were available for averaging.
An averaging over the complete particle set would e.g.strongly pronounce the many un-aggregated initial particles with a compactness value of 1.To obtain a statistically relevant averaging, only the states are reported where at least 100 aggregates consisting of 15 primary particles were present.As the abscissa, the average number of primar y par ticles per aggregate of the whole system is used.As discussed above, the simulation is stopped when the average number of primary particles per aggregates reaches 10.In cases 1, 4 and 7, the shear rate is already too large to allow the formation of any aggregate with 15 primary particles.In all other cases, it takes a while until enough aggregates with 15 primar y particles have been formed.The plots reveal that the compactness of the particles increases with increasing shear rate (case sequences: 2→5→8, 3→6→9).If particles reach the same size (number of primary particles) in different shear environments, higher shear rates promote the formation of compact structures.Similarly, an increasing growth rate (case sequences: 2→3, 5 →6, 8→9) leads to more open structures as the particles tend to stick to each other regardless of their relative orientation.

Conclusion
The presented Monte Carlo framework facilitates the investigation of the structural and size evolution of particle aggregates.In contrast to alternative approaches, no predefined assumptions on the aggregate shape and structure need to be made.However, the actual aggregation event has to be considered in very high detail.Because the fully detailed geometric information is available for modelling the aggregation event, this poses new challenges on mechanistic modelling.The current assumptions can still be considered very restrictive in terms of predictive power.Simplifications have mainly been chosen for computational reasons.However, even with this simplified approach, a mechanistic modelling of the expected behaviour can be realized.Note that the contribution of this work is not the phenomenological results discussed in the previous section.They actually represent the desired and expected behaviour.Instead of reproducing this behaviour with a descriptive model, the present paper aims at providing a method to address these changes in a mechanistic way.The presented technique is a step towards predictive modelling of the intrinsic kinetics of aggregation, and allows the simulation of highly complex aggregate shapes.Using the current simplifications, the simula-tion of fully complex shapes of particle aggregates is possible with state-of-the-art desktop computers.The computation time was in the order of 4 h for each run.With more detailed mechanistic models for the aggregation event, as they are planned in future work, one might use a multiprocessor machine.These mechanistic models will then be validated with experimental results and compared with alternative modelling approaches.

Nomenclature
C collision frequency di,sgg,subst apparent size of particle i derived from the substitute system G growth rate Li,j axis j of the substitute system of particle i M dimensionless number characterising the aggregation efficiency Ni number concentration for the size class i r effective aggregation rate Δt inter-event time Vagg aggregate total volume (sum of primary particle volumes) V control volume of the Monte Carlo scheme Xi,agg,subst data matrix for the substitution system β effective aggregation rate kernel Ψ aggregation efficiency γ ˆ shear rate σ* apparent yield stress μ kinematic viscosity ω compactness

Fig. 1
Fig.1Full and substitute characterization for an aggregated and a primary particle.

Fig. 2
Fig. 2 Collision efficiency as a function of the dimensionless group M.

Fig. 3
Fig.3Size-dependent efficient aggregation rate kernel.Size is given as the volume equivalent sphere diameter of the substitute ellipsoid.For better readability, the plot shows the rate kernel evaluated at the centre of each size interval.In the simulation, the discrete values on the intervals are used.

Fig. 4
Fig. 4 Number concentration as a function of time for the different cases 1-9.

Fig. 5
Fig. 5 Two-dimensional histogram representation of the final particle ensemble.Each bar shows the frequency of occurrence of aggregates in a certain interval of compactness and primary particles per aggregate.

Fig. 6
Fig. 6 Selected, representative particle morphologies for the cases 1, 3, 7 and 9.The aggregates were sampled randomly from the active set.

Table 1
Parameters for the 9 cases simulated in this study