Modeling Polymer Microencapsulation Processes using CFD and Population Balance Models

1. Introduction

Microencapsulation is a process that involves enveloping solids, liquids, or gaseous materials with a thin polymeric shell to form microcapsule particles. These systems are classified based on the particle size of microparticles, microcapsules, or microspheres. Emulsification is commonly used in microencapsulation processes, where the core material is dispersed in an organic solvent containing the wall-forming polymer. This dispersion is then emulsified in a non-miscible liquid, leading to the formation of droplets surrounded by a polymeric shell. The organic solvent is subsequently evaporated, leaving behind microcapsules with the active ingredient encased within a thin polymeric membrane [1,2,3]. Along with the influence of drop sizes on mass and heat transfer, the droplet size distribution plays a crucial role in the design and scale-up of essential processing equipment such as chemical reactors, mixers, and separators [4]. The ability to control and optimize the droplet size distribution is essential for ensuring the efficient performance and scalability of these unit operations.

Microencapsulation processes often involve complex fluid dynamics, mass transfer, and interfacial phenomena that are challenging to fully understand and optimize using experimental methods alone [5]. The heterogeneous nature of multiphase systems, such as emulsions and suspensions, encountered in microencapsulation further complicates the process and requires a comprehensive understanding of the interactions between various components [6]. Additionally, scaling up microencapsulation processes from lab-scale to industrial-scale production poses significant challenges, as fluid dynamics and mass transfer characteristics can change significantly with varying operating conditions and geometries [7]. Computational Fluid Dynamics (CFD) modeling can play a crucial role in addressing these challenges by providing valuable insights into the underlying mechanisms, enabling virtual optimization and screening of process parameters, and facilitating the successful scale-up of microencapsulation technologies [8]. The development of robust CFD models can significantly improve the efficiency, effectiveness, and quality of microencapsulation processes, leading to enhanced product performance and broader applications in fields such as drug delivery [9,10], the petroleum industry [11], food processing [12,13] or cosmetics/personal care [14,15].

The population balance approach is widely used to model the size distribution of droplets, bubbles, or crystals that evolve and change due to various phenomena such as nucleation, growth, coalescence, and breakage within a flowing system. The population balance model is a balance equation that tracks the changes in the size distribution of these dispersed species, similar to mass, energy, and momentum balances. Early research proposed a solution for the population balance equation in a well-mixed batch system, utilizing the moment transform technique to convert the population balance equations into a set of ordinary differential equations, which enabled tracking of size distribution changes [16]. Later, the moment transform approach was employed to develop a numerical technique for modeling the growth and aggregation of particles in a suspension, such as calcium oxalate monohydrate crystals, facilitating the simulation of kidney stone formation [17].

In recent years, advanced CFD techniques have been utilized to simulate multiphase flow processes. Some researchers have attempted to couple the population balance model (PBM) with transport equations through CFD modeling of slurries, emulsions, and gas-liquid systems. For instance, CFD coupled with a PBM has been used to investigate the aggregation of solid particles in a slurry system [18]. In this work, the Eulerian approach for multiphase flow and a k-ε turbulence model were used to simulate the turbulent two-phase flow. Other authors employed a similar modeling approach to estimate droplet size evolution in an oil-in-water emulsion [19]. In both instances, the researchers reported the mean particle or droplet size evolution over time but did not consider the size distribution of the particles, focusing only on particle growth.

A discrete PBM method coupled with a Eulerian-Eulerian approach for multiphase flow and a $k-\epsilon$ turbulence model was validated for a gas-liquid (air in water) dispersion in a tank [20], using the moving reference frame (MRF) technique and Luo's model for bubble breakage and coalescence [4]. In another study, a discrete PBM method was solved. using a sequential solution approach—first simulating the turbulent flow field and then solving the population balance equations. Due to the extensive computational time required for calculating hydrodynamic variables coupled with population balance equations, the authors recommended parallelizing the model in future work to reduce computational costs [21].

A mixture multiphase flow model in combination with a $k-\epsilon$ turbulence model, along with Lehr’s model for droplet breakage and coalescence [22], was also tested for liquid-liquid separation in hydro cyclones [23]. Several researchers have coupled CFD flow solutions with a Direct Quadrature Method of Moments (DQMOM) population balance model to study aggregation and breakage processes in one-dimensional domains. These studies mainly reported on computation times, numerical methods, and benchmarks of different CFD packages [24,25,26]. The effect of oil phase viscosity in droplet formation frequency in an oil-in-water microchannel emulsification system has also been studied by means of CFD [27]. In this study, the authors compare drop size distributions obtained through CFD modeling against experimental results; unfortunately, no relevant details about the simulation strategy adopted are facilitated.

Many studies have focused on modeling droplet breakage [28,29,30,31] and coalescence [32,33,34,35] phenomena, which are significant in the design of dispersion processes. However, predicting particle sizes in inhomogeneous flows with spatially and temporally varying velocity magnitudes and directions remains a challenge, especially at the microencapsulation particle length scales. Breakup and coalescence processes determine source terms in the population balance equation, while velocity vectors from the transport equations determine the convective term. An important computational issue is the solution method: sequential (first calculating the flow field and turbulence, then solving the population balance) versus simultaneous (solving all equations together from the start). The accuracy of predictions is influenced by the chosen solution method [27].

Several authors have combined CFD and PBM for particle size distribution prediction in emulsion processes. A study on laboratory-scale mixing tanks investigated the particle size distribution in water-in-oil emulsions using Conroe and Troika oil, examining the effects of different impeller speeds (300, 400, 600 rpm) on water droplet breakage and coalescence [36]. The study presents a computational model using Fluent 6.3.26. A $k-\epsilon$ turbulence model, a Euler-Euler multiphase model, and PBM were selected to simulate the mixing process. Validation through comparison with experimental data showed satisfactory similarities. Other studies extended this work by coupling the PBM model with a Euler-Lagrange multiphase model to see the effect of the droplet’s trajectories on the size distribution and [37,38]. Another study compared experimental and simulation data for a baffled mixing tank, using a Nexbase oil mixture and 3.5% NaCl brine at various agitator speeds [39]. By altering impeller speeds and introducing surfactant SPAN 80, the authors obtained experimental particle size distributions using particle vison and measurement (PVM) probes. Further research [19] highlighted the limitations of this methodology for particle size predictions in low-volume fraction and three-phase flow problems in double emulsion processes.

There are few references in the literature regarding the application of CFD+PBM for nanoparticle production processes (~10¹-10² nm), particularly involving high-shear agitators [40]. The existing studies primarily focus on supercritical anti-solvent [41] and flash nanoprecipitation processes [42,43]. For microcapsule production via continuous emulsion, CFD+PBM approaches have been employed to determine power consumption scale-up parameters [44]. However, no references were found addressing particle size distribution predictions in these processes, indicating a gap in current research. This underscores the need for further investigation to enhance process optimization and control in nanoparticle production.

In the present study, a CFD model of a polymer microencapsulation process, created using ANSYS FLUENT 2020 R2, is presented. A discrete PBM is used to predict particle number densities and cumulative size distributions obtained through single emulsion. Obtained values at different operating conditions were then compared against experimental data obtained from laboratory tests and errors are computed and analyzed. This approach aims to enhance the understanding and optimization of polymer microencapsulation processes, contributing to improved product and process performance, and broader applications in various fields.

2. Materials and Methods

The present study was conducted in two main phases: the first focused on the experimental polymer microencapsulation process, while the second involved the development and application of the CFD simulation model. The experimental phase focused on particle formation and size distribution, utilizing optical techniques for characterization. The numerical modeling phase employed CFD simulations to replicate the experimental conditions. An evaluation of the applicability of the CFD method for polymer encapsulation modeling was performed based on the results obtained.

2.1. Experimental Methodology

In the experimental setup, 1% w/v polyvinyl alcohol/water solution (PVA) was prepared as the continuous phase, and 20% w/v polyethyleneimine/dichloromethane solution (PEI) served as the dispersed phase. Table 1 presents the physical properties of the materials used.

The procedure (Figure 1) began by placing 90 ml of PVA in a beaker equipped with a high shear IKA S18N-19G dispersing tool [45]. 10 ml of PEI were introduced from the top using a syringe, with the agitator impeller initially set to the lowest speed to induce turbulence and prevent PEI from sticking to the beaker bottom. Once the mixing domain stabilized, the impeller speed was adjusted to the test velocity. The dispersed phase flow rate was set to 0.33 ml/min, taking approximately half an hour to complete the PEI dispersion into PVA, marking the completion of the microencapsulation process. For this study, the impeller speed was set to 5000, 7500, 10000, 12000, and 15000 rpm.

After the microencapsulation process, the PEI/PVA dispersion obtained for each tested agitator velocity was washed in a Sorvall RC5B plus refrigerated centrifuge (Thermo Scientific, Waltham MA, USA) for three cycles. A few drops extracted from each processed sample were placed on aluminum specimen mounts, stored in a box, and subjected to -80°C for 24 hours before being dried in a freeze dryer using a Telstar LyoQuest apparatus [46].

Scanning electron microscopy (SEM) and dynamic scattering light (DSL) techniques were used to obtain qualitative and quantitative data on the morphology and size of the prepared particles. The DSL technique specifically measured the obtained particle sizes at different impeller speeds. The data was then analyzed, plotted, and used to validate the CFD model, as discussed in the following sections. For SEM, a Focus Ion Beam Zeiss Neon 40 instrument (Carl Zeiss, Oberkochen, Germany) was used. For DLS, a NanoBrook 90Plus Zeta Instruments (Brookhaven Instruments, Holtsville NY, USA) was used. Figure 2 shows examples of the SEM images obtained for each tested agitator speed.

2.2. CFD Simulation Strategy

A geometrical model of the experimental setup for the simple emulsion microencapsulation process was constructed using Ansys DesignModeler. The rotor and stator sections of the agitator were modeled as separate volumes connected through a numerical interface, enabling the use of a MRF approach for the rotational frame during the CFD simulation setup.

A tetrahedral mesh was generated using Ansys Meshing, with a global average element size of 0.4 mm for the entire domain. To ensure mesh independence, the surface meshes of the agitator rotor and stator walls were successively refined. Steady-state MRF CFD simulations were conducted using water as the operating fluid, with the agitator rotation speed set to 15000 rpm. The rotor blades' torque was recorded for each case and used as the mesh convergence criteria. A local average element size of 0.1 mm was selected for the agitator surface mesh, meeting the mesh independence requirements. The final mesh, shown in Figure 3, consisted of approximately 6400 kEl. Details on the mesh construction process and independency analysis can be found in [47].

Ansys Fluent 2020 R2 [48] was used to generate the CFD model, set up for a turbulent, steady-state flow. A MRF approach was utilized, with rotor blades' rotational speed set between 5,000 and 15,000 rpm to replicate experimental conditions. The model employed a Realizable k-ε model to account for turbulence, and a Eulerian implicit approach was selected for the multiphase modeling. Additionally, a discrete PBM with five particle size bins was implemented, using Luo’s model as the breakage kernel [4] and a turbulence aggregation kernel [48]. Table 2 summarizes the CFD model numerical setup.

To validate the CFD+PBM simulation strategy, a model was constructed to replicate the study by Roudsari et al. [29], which experimentally and numerically investigated a 15% v/v Conroe oil-water emulsion system at rotational speeds of 300, 400, and 600 rpm. As shown in Figure 4, there is a good agreement between the experimental results reported by Roudsari et al. and the outcomes of the CFD model, demonstrating the model's accuracy in simulating the emulsion system and its reliability for predicting droplet size distributions in emulsification processes.

The emulsion microencapsulation simulations were executed on a workstation equipped with an Intel® Xeon® W-2265 CPU operating at 3.50 GHz. This processor features 12 cores with hyperthreading technology, resulting in a total of 24 logical processors available for parallel computations. The workstation also provides 128 GB of RAM to meet the memory requirements of the simulations. All cases were parallelized across 12 processors, achieving an average wall-clock time of 3.7 seconds per iteration with peak RAM usage reaching approximately 10.5 GB. Local bin fraction and flow velocity monitors were added to check for convergence. Simulations took up to 25000 iterations/26 hours to converge.

3. Results and Discussion

For all studied cases, the numerical results obtained for particle sizes and cumulative size distributions at different agitation speeds were compared against the SEM and DLS results obtained from the single emulsion experiments.

3.1. Experimental Particle Size Distributions

DLS was used to determine the experimental size distribution. For each agitator rotational speed tested, three samples were processed. Table 3 summarizes the obtained average particle sizes for each case.

Figure 5 shows an example of the size distributions obtained. The results indicate that for rotational speeds up to 10,000 rpm, the average particle diameter remains around 10³ nm. Beyond this speed, the average particle size drops almost by half with a 25% increase in agitation speed. It is known that the emulsification process exhibits two regimes based on droplet and turbulent flow interactions: turbulent inertial and turbulent viscous. In the turbulent inertial regime, drops are larger than the smallest turbulent eddies in the continuous phase, while in the turbulent viscous regime, drops are smaller than these eddies. Lower impeller speeds result in larger droplets due to reduced turbulent energy to break down the dispersed phase, whereas higher speeds produce smaller droplets as greater turbulent kinetic energy allows smaller eddies to fragment the drops [49].

3.2. CFD + PBM Model Performance

Due to the computational expense of solving PBM equations simultaneously with multiphase flow and turbulence equations, only 5 bins were used to compute particle number densities. Although this limitation can introduce discretization errors, literature suggests that 5 bins are sufficient for satisfactory results [27]. To validate the model, simulated droplet size distributions at different impeller speeds were compared to experimental data from DLS analysis.

3.2.1. Particle Number Densities

Figure 6 presents the particle number density contours from the CFD simulations for three different particle sizes and agitation speeds in a radial plane intersecting the agitator impeller. The results show higher particle number densities for larger particles at 5000 rpm, and the formation of smaller particles significantly increases with higher agitator rotational speeds. This is consistent with the experimental observations and the DLS data shown in the prior section.

The predicted number density contours for the smallest bin size (6.0 x 10^-6 m) at 10000 rpm and 15000 rpm display similar particle densities across most of the plane, whereas at 5000 rpm, these smaller particles are confined to the high-shear agitator area. This indicates that achieving turbulent viscous regimes throughout the entire flow domain is crucial for preventing the coalescence of small droplets formed near the agitator when they move to regions with lower turbulence. By visually comparing the contours for the smallest bin size, it becomes evident that the mixing system reaches a turbulent inertial regime at an agitation speed of 10000 rpm, ensuring consistent particle size distribution across the domain, thus optimizing the emulsification process.

3.2.2. Particle Size Distribution

Figure 7 shows the bin fraction for various particle sizes at different rotational speeds, comparing CFD simulations (blue bars) with experimental data (orange lines). At low agitation speeds (5000/7500 rpm), the experimental results show larger fractions for larger particle sizes, indicating that lower agitation speeds favor the formation of larger particles, consistent with a turbulent viscous emulsification regime. For these low-speed cases, The CFD simulations underestimates these larger particles, probably signaling that the turbulence model is overestimating the turbulent kinetic energy (thus obtaining smaller eddies) under those flow conditions.

As the rotational speed increases (10000 rpm or above), the experimental and CFD results exhibit a better agreement across all particle sizes, reflecting a transition towards a turbulent inertial regime. However, the simulations tend to slightly underestimate the fraction of the smallest particles. The cumulative distribution of particle diameters for various rotational speeds presented in Figure 8 highlight the aforementioned. The CFD model fails to capture the size distribution trend at lower speeds (5000 rpm and 7500 rpm), likely making it not suitable for turbulent viscous emulsification conditions

The particle size distribution obtained from simulations was analyzed using D10, D50, and D90 metrics, representing the particle diameters at the 10th, 50th, and 90th percentiles, respectively. Figure 9 illustrates these results, providing a clearer view of the trends observed in the cumulative size distributions shown in Figure 8.

Figure 10 presents the computed relative errors by bin size and agitation speed, revealing higher errors in particle size distribution at lower agitation speeds. It is worth noticing that the largest relative errors occur for smaller particle sizes. Despite minor discrepancies, the CFD simulations closely align with the experimental data at higher agitation speeds (>10000 rpm), validating the simulation strategy's accuracy in predicting particle size distributions under turbulent viscous operational conditions, with acceptable errors below 10%. However, the model fails to accurately capture the flow characteristics at lower agitation speeds (<7500 rpm) within the turbulent inertial regime, resulting in significantly higher mean errors, meaning that the simulation is unable to accurately capture the flow field and dynamics under those operating conditions. This behavior was also observed in the validation case results (Figure 4).

At lower agitation speeds, as the turbulence level decreases, the flow regime may transition to a turbulent viscous regime, where different mechanisms govern droplet breakup and coalescence. The CFD model may struggle to accurately capture these regime-dependent phenomena at lower flow velocities, likely due to the performance of the turbulence model used [50]. As turbulence intensity increases, the influence of interfacial forces, such as surface tension and viscous forces, becomes less significant compared to the turbulent forces acting on the droplets [51], leading to better agreement between experimental and simulation results at higher speeds.

These results highlight the necessity for further research to ensure suitable turbulence model selection for the CFD modeling process. A known issue of the Realizable k-ϵ model is its tendency to overpredict turbulent kinetic energy and underpredict energy dissipation in high mixing flow conditions [52], which may explain its higher inaccuracy at low agitation speeds. Exploring scale-adaptive turbulence models could potentially improve the accuracy of CFD simulations across varying flow regimes. Additionally, investigating alternative turbulence models that can better account for the transition between different flow regimes could enhance the model's predictive capabilities at lower agitation speeds.

4. Conclusions

Computational fluid dynamics (CFD) proves to be a valuable tool for gaining insight into complex processes such as polymer microencapsulation by emulsification. The combination of a population balance model with a Eulerian multiphase framework and a moving reference frame approach is effective in predicting number densities and particle size distributions for such processes.

The model successfully captures the intricate fluid dynamics involved, demonstrating good agreement between experimental and simulation results. Specifically, at higher agitation speeds (>10000 rpm), the CFD model accurately represents the underlying turbulent flow phenomena. The agreement between simulated and experimental data at these speeds indicates that the model effectively captures the turbulent inertial flow conditions, validating its capability to predict droplet behavior and size distribution under high-shear conditions.

However, discrepancies at lower speeds (<7500 rpm) reveal challenges in accurately capturing turbulent viscous regimes. In these low-speed cases, CFD simulations tend to underestimate larger particles, likely due to an overestimation of turbulent kinetic energy, resulting in smaller predicted eddies. This discrepancy at lower impeller speeds underscores the importance of comprehensive experimental validation across a wide range of operating conditions.

While the current CFD model demonstrates overall good predictive accuracy, its performance at lower agitation speeds needs improvement. Future work should focus on extensive experimental validation across a broader range of operating conditions, to better characterize the turbulent viscous and turbulent inertial regimes. This will help identify and address the limitations of the current CFD model, ensuring its reliability and accuracy for different flow conditions. The integration of more advanced turbulence models and improved computational techniques could further refine the predictions and provide deeper insights into the droplet dynamics and particle size distributions in polymer microencapsulation processes.

By addressing the identified discrepancies and enhancing the model's capabilities, CFD can become an even more powerful tool for optimizing and controlling microencapsulation processes, ultimately leading to more efficient and effective industrial applications. The continued development and validation of CFD models will play a crucial role in advancing our understanding of complex fluid dynamics and improving the design and scale-up of polymer microencapsulation processes.