COMSOL 6.4 - Turbulent Aggregation of Nanoparticles

Nanoparticles have become widely used in a range of applications. For example, in pharmaceuticals, particles on the nanoscale, such as nanocrystals, are often used as an effective means of drug delivery. Meanwhile in material sciences, nanoparticle coating is an increasingly common way to improve or change the properties of materials. Understanding the factors affecting the size of these particles is an important step to develop production processes that optimize the desired particle properties.

Model Definition

The term nanoparticle is used to describe particles on the nanoscale. This example models the turbulence-induced aggregation of nanoparticles in the size range of a few hundred nanometers in a beaker fitted with a square blade impeller to predict the size distribution of volume-equivalent spherical particle diameters. The model is set up with a 2D geometry describing the aggregation and movement of particles in a turbulent flow in the mixing plane (Figure 1). The entire domain is modeled as a rotating domain and the impeller is set to a rotational speed of 350 rpm.

Figure 1: Model geometry consisting of beaker with impeller.

The aggregation and flow is modeled with the Precipitation in Fluid Flow interface which couples a Size-Based Population Balance interface with a Single-Phase Fluid Flow interface. A k-ε model is used to account for the turbulent flow resulting from the mixing. In this example we exclude nucleation, growth and dissolution. Furthermore, the particles are also assumed small enough that breakage is negligible. The size distribution is modeled using a discretization of 20 size intervals, or bins, and a logarithmic interval distribution. For aggregation, logarithmic distributions are recommended. The aggregation contribution is modeled as (Ref. 1)

(1)

where L is the characteristic length and K is the aggregation kernel. This model assumes spherical particles in which case the characteristic length, L, is the diameter. The aggregation kernel is written with both a perikinetic and orthokinetic contribution. They describe the collisions due to Brownian motion and convective transport by the fluid respectively. They are given by (Ref. 2 and Ref. 3)

(2)

and

(3)

where Lj and Lk are the sizes of the two colliding particles. Furthermore, ν is the kinematic viscosity (m2/s), ε is the turbulent dissipation rate (m2/s3), ρ is the fluid density (kg/m3), kb is the Boltzmann constant (m2·kg/s2/K) and T is the temperature (K). The turbulent dissipation rate are solved for by the fluid flow interface. All particles are modeled as spherical where the size, L, is the diameter such that the system, in absence of nucleation, growth and dissolution, conserves the total particle volume. Assuming that the collision contributions (Equation 2 and Equation 3) are independent of each other and that the contributions are additive we can write

(4)

Here, W is the stability ratio. In this model a joint stability ratio is applied to both collision contributions. The stability ratio can range from 1 (no stabilization) to infinity (complete stabilization) and the specific value will depend on the system. Experiments are often required to find a suitable stability ratio. For example, if ions are adsorbed on the particle surface, the stability ratio will increase due to repelling electrostatic forces. Or, if the solution is supersaturated, it is possible for the crystal particles to form strong bonds using the solute, in the form of agglomerates, which gives a lower stability ratio. As was done by Schwarzer and Peukert in their model of barium sulfate nanoparticle aggregation (Ref. 2), we therefore test multiple stabilization ratios. Three weak stabilization scenarios are tested; W = 1, W = 10, and W = 100. A stabilization of W = 1 is equal to having no stabilization and is a common approximation for systems with negligible electrostatic forces.

At time t = 0, there is an arbitrary initial distribution of particles according to (modified from Ref. 4)

(5)

if L > Loffset and 0 otherwise. The constants N0, Lmean, and Loffset control the initial particle concentration, mean size, and starting point, respectively. The distribution can be seen in Figure 2.

Figure 2: Initial size distribution.

Results and Discussion

The velocity field can be seen in Figure 3. The largest velocity magnitude is located at the tips of the impeller. This is also where we will see the most orthokinetic aggregation.

Figure 3: Velocity field in the mixing plane.

Figure 4 shows the concentration of particles at 0.2, 0.4, 0.8, and 1.2 seconds after mixing has begun with w = 1. A more pronounced decrease in particle concentration is seen near the impeller where the turbulent dissipation rate is higher and the orthokinetic aggregation contribution is larger.

Figure 4: Particle concentration at different times for W = 1.

Figure 5 shows the average particle concentration and average size across the beaker. The average size can be seen to increase from 200 to 388, 229, and 203 nm while the particle concentration decreases from 6 × 1016 to 3.1 × 1015, 3.6 × 1016, and 5.7 × 1016 1/m3 for W = 1, 50, and 100, respectively.

Figure 5: Average particle size and number of particles.

Finally, Figure 6 shows the average size distribution at different times and stability ratios. The distribution is given at t = 0, 3, 6, and 10 seconds. For no stabilization, W = 1, the distribution almost immediately collapses into a smaller number of larger particles. This is also what was seen in Figure 5. The distribution for W = 10 flattens out at a slower rate. At this stability ratio the change in distribution profile is still pronounced. The combinations of sizes that can form a particles are fewer at smaller sizes, so for small particles we see that the death term (sink term) dominates. At the highest modeled stability ratio, W = 100, aggregation still takes place but the overall distribution profile does not change much over the modeled time span.

Figure 6: Size distribution at different times.

References

1. M. Ölander, Numerical Simulations for Battery Recycling, master’s thesis, KTH, Royal Institute of Technology, 2023.

2. H.C. Schwarzer and W. Peukert, “Combined experimental/numerical study on the precipitation of nanoparticles,” AICHE J., vol. 50, no. 12, pp. 3234–3247, 2004.

3. A. Raponi, S. Romano, G. Battaglia, A. Buffo, M. Vannie, A. Cipollina, and D. Marchisio, “Computational Modeling of Magnesium Hydroxide Precipitation and Kinetics Parameters Identification,” Cryst. Growth Des., vol. 23, no. 12, pp. 4748-4759, 2023.

4. W. Zhang, T. Przybycien, J.M. Breuer, and E. von Lieres, “Solving crystallization/precipitation population balance models in CADET, Part II: Size-based Smoluchowski coagulation and fragmentation equations in batch and continuous modes,” Comput. Chem. Eng., vol. 192, 2025.

Chemical_Reaction_Engineering_Module/Mixing_and_Separation/turbulent_aggregation

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select .

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file turbulent_aggregation_parameters.txt.

Definitions

Initial distribution

In the toolbar, click

In the window for , type Initial distribution in the text field.

In the text field, type initdist.

Locate the section. In the text field, type (L>L_offset)*3*(L-L_offset)^2*(N0)/((L_mean-L_offset)^3)*exp(-(L-L_offset)^3/(L_mean-L_offset)^3).

In the text field, type L.

Locate the section. In the text field, type 1/m^4.

In the table, enter the following settings:


Argument	Unit
L	m

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	L	0	400[nm]	0	m

Click

Geometry 1

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type D_beaker/2.

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file turbulent_aggregation_impeller_geom.txt.

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the object only.

Click

Click the

button in the toolbar.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

Component 1 (comp1)

Rotating Domain 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the f text field, type rpm_impeller.

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Wall 1

In the window, under > ε (spf) click .

In the window for , click to expand the section.

From the list, choose .

Wall 2

In the toolbar, click

and choose .

Select Boundaries 1–8 only.

Pressure Point Constraint 1

In the toolbar, click

and choose .

Select Point 6 only.

Size-Based Population Balance (pbsb)

In the window, under click .

In the window for , locate the section.

In the ρp text field, type rho_p.

Locate the section. From the list, choose .

In the I text field, type 20.

In the L0 text field, type L_offset.

In the LI text field, type 2000[nm].

Locate the section. Select the checkbox.

In the K text field, type 1/W*(sqrt(pi/(15*8))*(pop.Lj+pop.Lk)^3*sqrt(max(ep,eps)/spf.nu)+2*k_B_const*T/(3*spf.nu*spf.rho)*(pop.Lj+pop.Lk)*(1/pop.Lj+1/pop.Lk)).

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the n0 text field, type initdist(pop.L).

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size 1

In the window, under > click .

Select Boundaries 1–8 only.

Size 2