Modeling the Impact of Non-Ideal Mixing on Continuous Crystallization: A Non-Dimensional Approach

Research Article - Peer Reviewed Conference Proceeding ESCAPE 35 - European Symposium on Computer Aided Process Engineering Ghent, Belgium. 6-9 July 2025 Jan F.M. Van Impe, Grégoire Léonard, Satyajeet S. Bhonsale, Monika E. Polańska, Filip Logist (Eds.) [DOI_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 1 Modeling the Impact of Non-Ideal Mixing on Continuous Crystallization: A Non-Dimensional Approach Jan Trnkaa, František Štěpáneka* a University of Chemistry and Technology, Department of Chemical Engineering, Prague, Czech Republic * Corresponding Author: [email protected] ABSTRACT Mathematical modeling is essential for the effective control of many chemical engineering processes, including crystallization. However, most existing crystallization models used in industry and academia assume ideal mixing. As a result, the unclear effects of imperfect mixing on crystallization, reported in experimental studies, remain largely unexplained. In this work we aim to address this gap in understanding by examining antisolvent crystallization processes on a general theoretical level, using a novel dimensionless model. To address the impact of mixing on crystallization, we employ the Engulfment model coupled with a population balance, and we nondimensionalize the model equations. Using this model, we explore the dependence of the mean particle size on the homogenization rate, represented by the Damköhler number for crystallization. Moreover, we study the impact of mixing at various values of the model's kinetic parameters to simulate difference in properties of individual products. We show that we are able to explain the complex interaction between crystallization and mixing, proving our model can serve as a tool for achieving a better understanding of the processes involved. Finally, due to its efficiency and reduced number of parameters, the model is suitable for direct fitting to experimental data. Keywords: crystallization, modeling, mixing, continuous, non-dimensional INTRODUCTION Crystallization is a fundamental technique used for separating and/or purifying solids in various industries, including pharmaceuticals and food production. The process conditions during crystallization significantly affect the particle size and purity of the final product. Both particle size and shape influence various material properties, such as flowability, filterability, and dissolution behavior. Effective control of the crystallization process can enhance product quality and eliminate the need for additional steps like comminution or granulation. Unlike cooling or evaporative crystallization, reactive and antisolvent crystallization involve mixing two liquids of distinct compositions. Supersaturation is induced through mixing, either due to the product being synthesized (reactive) or the reduction of solubility (antisolvent). However, most models currently used in both industry and academia assume ideal mixing and fail to account for the system's dependence on mixing dynamics, thereby limiting our understanding of the process. Numerous experimental studies have addressed the impact of mixing on crystallization. The effect of increasing mixing intensity on particle size is inconsistent, varying across different substances and products [1]. The mean crystal size may increase, decrease, or even reach a maximum or minimum as the mixing intensity varies. Even using the same crystallizer type for the same product may lead to qualitatively different outcomes. Additionally, the rate of homogenization has also been shown to affect the width of the particle size distribution (PSD), particle shape, and/or the extent of agglomeration. Although these dependencies are well documented, a systematic approach to studying them has not been established. Research, both on the modeling and experimental sides, tends to focus on a single system at a time. In contrast to the current trend of developing complex models, we aim to use a relatively simple modeling framework, keeping the number of parameters low and enabling extensive parametric studies necessary for our theoretical research. The authors have already sufficiently explained the [LAPSE_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 2 mechanism behind the existence of a maximum in the dependency of particle size on mixing intensity for the semibatch case [2]. We now aim to extend our research using a novel dimensionless model for continuous crystallization. As in other areas of chemical engineering, nondimensionalization offers valuable insights by reducing the number of parameters and revealing characteristic system properties. It also provides suitable scaling of variables, improving both the precision and efficiency of numerical computations. MODELING Current Modeling Approaches The population balance equation (PBE) is the standard approach for crystallization modeling. Depending on the level of detail in describing the hydrodynamics involved, three main modeling approaches are commonly considered: ▪ The ideal mixing assumption is still widely used as it allows for neglecting mixing in the model equations, thus simplifying calculations. This approach assumes an infinitely fast homogenization rate and is therefore applicable to systems that do not exhibit a dependence on the agitation rate. ▪ Models based on compartmentalization represent a balanced approach to account for the complexity of mixing. These models divide the studied volume into several compartments and describe the material flux between them. Within each compartment, mixing is typically assumed to be infinitely fast. One of the most popular models in this category is the mechanistic micromixing model developed by Baldyga and Bourne, also known as the Engulfment model [3]. ▪ CFD-based approaches, where the PBE is integrated within the CFD framework, are perhaps the most rigorous. However, even CFD cannot fully capture mixing processes at the molecular level (i.e., micromixing). Since crystallization is a molecular process, these methods often require additional micromixing models. Alternatively, they assume that the homogenization process is limited only by macroscale mixing (i.e., fast micromixing). Despite the complexity and high computational cost, the improvement in accuracy over the Engulfment model appears to be limited. Model Description Coupling PBE With the Engulfment Model We have based our approach on the Engulfment model, as we believe it provides the most suitable choice for our theoretical study. In this Lagrangian model, the system volume is discretized into two well-mixed regions: (1) the mixed zone, enriched with the reference compound, and (2) the surrounding bulk fluid. Mixing is described as the expansion of the mixed zone through the engulfment of the bulk, leading to dilution of the compound in the mixed zone if no additional source is present. This approach was initially developed in reactor engineering but has since been applied to crystallization (e.g., [4]). Although the model was originally developed based on turbulent flow mixing mechanisms, we find its framework applicable even to non-turbulent flows. According to the Engulfment model [3], the concentration of a compound in the mixed zone evolves over time according to Eq. 1 as volume fraction of the mixed zone (𝑋) expands. The rate of homogenization is determined by a single constant, the mixing time 𝑡𝑚𝑖𝑥 (Eq. 2). 𝑑𝑐 𝑑𝑡 = 1 𝑋𝑑𝑋 𝑑𝑡 (𝑐𝑏−𝑐)+𝑟 (1) 𝑑𝑋 𝑑𝑡 =1 𝑡𝑚𝑖𝑥 𝑋 (2) Coupling the engulfment model with PBE results in Eq. 3, accounting only for primary nucleation and growth (neglecting agglomeration and breakage). The molar balance is then described in Eq. 4. These equations are valid for both reactive and antisolvent crystallization in batch or in continuous tubular crystallizers operating at steady state, for which time corresponds to the coordinate time. 𝜕𝑓 𝜕𝑡 +𝐺𝜕𝑓 𝜕𝐿 =1 𝑡𝑚𝑖𝑥 (𝑓𝑏−𝑓), 𝑓(0,𝑡)=𝐽 𝐺 (3) 𝑑𝑐 𝑑𝑡 = 1 𝑡𝑚𝑖𝑥 (𝑐𝑏−𝑐)−3𝑘𝑣𝜌𝑐𝑟 𝑀𝐺𝜙2+𝑟 (4) Model Nondimensionalization To study the effect of mixing on crystallization at a theoretical level, we have developed an efficient method for nondimensionalizing Equations 3 and 4. We begin by introducing the following dimensionless variables: 𝜏 = 𝑡 𝑡0, 𝜆= 𝐿 𝐿0, 𝜒= 𝑐∗ 𝑐0, 𝑆 = 𝑐 𝑐0𝜒, 𝑛=𝐿0𝑉0𝑓 (5) 𝐺 =𝐺 𝐺𝑟𝑒𝑓 , 𝐽󰆻=𝐽 𝐽𝑟𝑒𝑓 , 𝑟 =𝑡0 𝑐0𝑟, 𝑐0=𝑐𝑟𝑒𝑓 (6) Defining: 𝑉0=𝑘𝑣𝜌𝑐𝑟𝐿0 3 𝑀𝑐𝑟𝑐ref ,𝐿0=𝑡0𝐺𝑟𝑒𝑓,𝑡0=( 𝑐0𝑀𝑐𝑟 𝑘𝑣𝜌𝑐𝑟𝐽𝑟𝑒𝑓𝐺𝑟𝑒𝑓 3)1 4 (7) leads to significant reduction in the number of model parameters. The reference values 𝐽𝑟𝑒𝑓,𝐺𝑟𝑒𝑓,𝑐𝑟𝑒𝑓 can be set arbitrarily, although they do affect scaling. The appropriate value for 𝑐𝑟𝑒𝑓 is the maximal solubility while the suggested choice of 𝐽𝑟𝑒𝑓 and 𝐺𝑟𝑒𝑓 is discussed later. Assuming isothermal conditions and solubility as a sole function of the solvent volume fraction 𝜑, the final [LAPSE_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 3 form of our non-dimensional model is described by the following equations: 𝑑𝑆 𝑑𝜏 =𝐷𝑎𝑐𝑟𝑆(𝑆𝑏𝜒𝑏 𝑆𝜒 +(𝜑−𝜑𝑏) 𝜒𝑑𝜒 𝑑𝜑 −1)+1 𝜒(𝑟 −3𝐺   2) (8) 𝜕𝑛 𝜕𝜏 +𝐺 𝜕𝑛 𝜕𝜆 =𝐷𝑎𝑐𝑟(𝑛𝑏−𝑛), 𝑛(0,𝜏)=𝐽󰆻 𝐺  (9) Derivation of the non-dimensional model leads to the emergence of a dimensionless number, 𝐷𝑎𝑐𝑟, defined as the ratio of homogenization rate to crystallization rate (Damköhler number for crystallization): 𝐷𝑎𝑐𝑟 =𝑡0 𝑡𝑚𝑖𝑥 =𝑣𝑚𝑖𝑥 𝑣𝑐𝑟 (10) We further assume expressions for growth and crystallization rates according to Eq. 11. By setting the reference values equal to the rate coefficients, we conveniently reduce the number of kinetic constants from four to two, as shown in Eq. 12. 𝐺 =𝑘𝐺(𝑆−1)𝑔, 𝐽= 𝑘𝐽exp(−𝑗 ln2𝑆) (11) 𝐺 =(𝑆−1)𝑔, 𝐽󰆻=exp(−𝑗 ln2𝑆) (12) The constants 𝑗 and 𝑔 together with the solubility data remain the only product properties needed as input to the model. Application to Continuous Antisolvent Process In this study, we focus on isothermal antisolvent crystallization in a tubular device developed in our research group. The process is schematically illustrated in Figure 1. A solution of candesartan cilexetil, an active pharmaceutical ingredient (API), in acetone is injected perpendicularly to the length direction of the tube into a stream of water. Mixing of the two streams generates supersaturation, which leads to nucleation and subsequent growth of the product particles. The initial conditions of the mixed zone are given by the properties of the organic phase (𝑆(0)=1,𝜑(0)=0.85) while bulk is represented by pure water. The volume flow rate ratio of organic to inorganic phase is 1:9 (corresponds to 𝑋(0)=0.1). Setting the bulk variables 𝑆𝑏, 𝜑𝑏 and 𝑛𝑏 to zero and omitting the reaction term results in a simplified set of equations: 𝑑𝑆 𝑑𝜏 =𝐷𝑎𝑐𝑟𝑆(𝜑 𝜒𝑑𝜒 𝑑𝜑 −1)−3𝐺   2 𝜒 (13) 𝜕𝑛 𝜕𝜏 +𝐺 𝜕𝑛 𝜕𝜆 =−𝐷𝑎𝑐𝑟𝑛, 𝑛(0,𝜏)=𝐽󰆻 𝐺  (14) The solubility data of candesartan cilexetil used for our study were taken from literature [5]. The model of the continuous device is implemented in the Python environment. Eq. 13 is converted into a set of ordinary differential equations (ODEs) using 1D finite volume method with Koren flux limiter. The set of ODEs is then integrated numerically using the Runge–Kutta method (RK45). All simulation results are steady with respect to the coordinate time, assuming an infinitely long tube. Figure 1. Schematic representation of the setup and spatial discretization of the continuous ASP device according to the Engulfment model (the shade of blue color represents local concentration of API). RESULTS The Mixing Impact on the Mean Particle Size As mentioned, our model takes three product-specific inputs: parameters 𝑗 and 𝑔 and the solubility function 𝜒(𝜑). In our study, we vary the kinetic parameters in order to address the influence of the product properties on the results of our simulations. The effects of changing the solubility curve are not presented in this study. As 𝑡0 is a constant for a given substance, increasing 𝐷𝑎𝑐𝑟 has the meaning of increasing the rate of homogenization. We proceed with a sensitivity analysis of the effect of kinetic constants 𝑗 and 𝑔 on the mean particle size (𝜆43) and its dependence on 𝐷𝑎𝑐𝑟. The Influence of 𝑔 on 𝜆43(𝐷𝑎𝑐𝑟) The results of varying 𝑔 along with 𝐷𝑎𝑐𝑟 at constant value of 𝑗 are depicted in Figure 2. As shown, all the reported scenarios for the dependency of the mean particle size on mixing intensity are covered by our model. In agreement with the experimental studies, the particle size according to our model may decrease, increase, reach a maximum or minimum or remain constant with change in the mixing intensity. All the parametric curves are shaped similarly. At very low 𝐷𝑎𝑐𝑟, the particle size decreases while increasing the mixing rate. This happens as the result of monotonous increase in maximal supersaturation with 𝐷𝑎𝑐𝑟 demonstrated on Figure 3. At higher supersaturation, more nuclei are formed due to enhanced nucleation rate, resulting in decrease in particle size. The crystal size reaches minimum at 𝐷𝑎𝑐𝑟 =󰇗 10−1. Interestingly, the parametric curves switch their order shortly before reaching the minimum. This happens as the maximal supersaturation reaches the value of two (the power function argument of the growth rate reaches one). For more intense mixing, 𝜆43 grows with 𝐷𝑎𝑐𝑟. This most likely happens because the increase in nucleation [LAPSE_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 4 rate diminishes at high supersaturation while the growth rate keeps accelerating due to the nature of the respective equations (𝐽󰆻 is limited unlike 𝐺 ). However, the increase of 𝜆43 is discontinued quite abruptly at values of 𝐷𝑎𝑐𝑟 unique for every curve. The value of the local maximum and its location increases with higher growth rate exponent. To understand the sudden decrease, let us describe the kinetics of crystallization according to our model in the phase space of the non-dimensional concentration (𝑆𝜒) and the solvent volume fraction (Figure 4). The initial state marks the composition of the organic phase. The yellow line represents the conditions close to perfect mixing, where the mixing and crystallization events are separate. The transition from the initial state to the mixed state along the straight line is caused by dilution of the island by engulfment of the bulk. Any deviation from the yellow line represents induced crystallization while mixing. The maximal possible supersaturation is reached at the yellow line at 𝜑𝑆𝑚𝑎𝑥 =󰇗 0.42. For 𝜑 >𝜑𝑆𝑚𝑎𝑥 the supersaturation always grows due to mixing while for 𝜑<𝜑𝑆𝑚𝑎𝑥 mixing causes decrease in 𝑆. If the mixing rate is too fast relative to the crystallization rate, crystallization is not induced before reaching 𝜑𝑆𝑚𝑎𝑥, causing the crystallization to begin at significantly lower supersaturation, reducing the final size of the crystals at high values of 𝐷𝑎𝑐𝑟. This behavior is encoded in Eq. 13 as the term 𝜑 𝜒𝑑𝜒 𝑑𝜑 −1 is positive at 𝜑 >𝜑𝑆𝑚𝑎𝑥 and negative at 𝜑<𝜑𝑆𝑚𝑎𝑥. The value of 𝜑𝑆𝑚𝑎𝑥 is therefore determined only by the shape of the solubility function. The Influence of 𝑗 on 𝜆43(𝐷𝑎𝑐𝑟) Let us now consider the scenario of varying 𝑗 along with 𝐷𝑎𝑐𝑟 at constant 𝑔. Overall, increasing 𝑗 hinters the nucleation rate and thus promotes growth, causing a general increase in particle sizes. Increasing 𝑗 also leads to later onset of nucleation induced at higher supersaturation, further favoring growth over nucleation. The results of the simulations are shown in Figure 5. For low values of 𝑗, the shape of the curves remains unaltered compared to the results in Figure 2. However, at values roughly from 10 to 80, the system undergoes a qualitative change. In this process, the local maximum vanishes and the decrease in 𝜆43 associated with approaching perfect mixing is turned into an increase. Previously, for 𝑗=1, the reason for the decrease in particle size was sudden drop in supersaturation due to fast mixing, resulting in less pronounced growth. However, nucleation at 𝑗 =100 is about hundred times slower. At these conditions, lower supersaturation produces significantly smaller number of particles, resulting in seemingly paradoxical increase in particle size. Figure 2. Dependence of the mean non-dimensional particle size on the Damköhler number for crystallization with varying the growth rate exponent. Figure 3. Evolution of the supersaturation over time (maximal supersaturation increases with 𝐷𝑎𝑐𝑟). The x-axis is scaled by 𝐷𝑎𝑐𝑟 for better comparison. Figure 4. Representation of the evolution of mixing and crystallization in the phase space of the non-dimensional concentration and the solvent volume fraction. [LAPSE_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 5 Figure 5. Dependence of the mean non-dimensional particle size on the Damköhler number for crystallization with varying nucleation constant. Parametric Fitting As the developed non-dimensional model has only two kinetic parameters in contrast to four in the dimensional one, it is easier to fit the model to experimental data. Moreover, using the non-dimensional model requires no knowledge about the product properties other than its solubility behavior. To make use of these advantages, we further present how to use the model for parametric fitting. Let 𝑳 denote the vector of 𝑁 measured crystal sizes at different volume flow rates 𝑽󰇗 (the mixing rate in our system is assumed to be dependent only on 𝑉󰇗). As evident from Eq. 15, dividing 𝑳 by one of its elements gives the same results for both dimensional and non-dimensional data: 𝑳 𝐿[0] =𝐿0𝝀 𝐿0𝜆[0] =𝝀 𝜆[0] (15) We use this identity to define the objective function as follows: 𝐹 =∑(𝐿[𝑖] 𝐿[0]−𝜆[𝑖] 𝜆[0])2 𝑁 𝑖=0 (16) The remaining problem to be solved is finding the link between the flow rates and the corresponding Damköhler numbers. Based on the research done on similar micromixers [6], we expect the mixing time to be inversely proportional to the volume flow rate to the power of 1.5. Thus: 𝐷𝑎𝑐𝑟 =𝑡0 𝑡𝑚𝑖𝑥(𝑉󰇗)=𝑡0 𝐾𝑉󰇗1.5 =𝑐𝑉󰇗1.5 (17) where the coefficient 𝑐 is unknown. After applying the same strategy as for the crystal sizes, we get: 𝑫𝒂𝒄𝒓 𝐷𝑎𝑐𝑟[0] =𝑐𝑽󰇗1.5 𝑐𝑉[0] 󰇗1.5 =𝑽󰇗1.5 𝑉[0] 󰇗1.5 (18) The points at which 𝝀 are to be evaluated from simulations are therefore: 𝑫𝒂𝒄𝒓 =𝐷𝑎𝑐𝑟[0] 𝑽󰇗1.5 𝑉[0] 󰇗1.5 (19) where 𝐷𝑎𝑐𝑟[0] is unknown and therefore it is another parameter to be optimized. Overall, finding the kinetic parameters presents an optimization problem: minimize 𝑔,𝑗,𝐷𝑎𝑐𝑟[0] 𝐹 (20) We have used experimental data measured by our group to test the use of this method and validate our model. The best fit of the model to the experimental data is presented in Figure 6. In the experiments, we have used two distinct mixing units, the T-junction and the FDmiX mixer. The latter device increases the homogenization rate by passively introducing flow oscillations at otherwise laminar conditions. We have used a combination of direct grid search and the genetic algorithm to solve the optimization problem. The model qualitatively describes the measured trends quite well. The value of 𝐷𝑎𝑐𝑟[0] for the FDmiX was found at higher values than for the T-junction as expected. As mixing is very slow in the T-junction, the particle size drops quite rapidly before it levels out as explained in describing Figure 2. However, aggregation may also contribute to the larger size of particles measured at low mixing rates. On the other hand, the FDmiX seems to operate at mixing rates close to the local maximum, which is very well predicted by our model. Ongoing CFD analysis of the flow in the T-junction has revealed that mixing is not finished before reaching the outlet, which our model is not accounting for. Expanding the model to address this issue may further improve the results in the future. Figure 6. The best fit of the model to our experimental data for two different mixers. [LAPSE_DoNotChange] Syst Control Trans 4:XXXX-YYYY (2025) 6 CONCLUSION In this work, we present a simple, yet efficient crystallization model developed for studying the interaction between crystallization and mixing in batch and continuous tubular crystallizers. We have found an efficient way to nondimensionalize the model equations, significantly reducing the number of parameters. As a result, a new non-dimensional number has emerged during the process – the Damköhler number for crystallization – representing the ratio of homogenization and crystallization rates. We use the model to study the process of continuous antisolvent crystallization by means of parametric sensitivity analysis. We were able to undercover the complex interaction between nucleation, growth and mixing and increase the understanding of the processes involved. In addition to its suitability for theoretical research, our model is also convenient for fitting the kinetic parameters to experimental data as the reduced dimension of the search space streamlines solving of the optimization problem. Despite the simplicity of the model, we have shown it is able to account for all the reported scenarios of mixing impact on particle size and to fit measured data sufficiently well. ACKNOWLEDGEMENTS The financial support is from the project OP JAK INTER-MICRO (registration number CZ.02.01.01/00/22_008/0004597), the Specific University Research (MSMT), co-funded by European Union. LIST OF SYMBOLS 𝑐∗,𝜒 solubility [mol m−3], [1] 𝑡, 𝜏 time [s], [1] 𝑓, 𝑛 population density [m−4], [1] 𝐿, 𝜆 crystal size [m], [1] 𝜙2,  2 second moment [m−1], [1] 𝐽, 𝐽󰆻 nucleation rate [m−3], [1] 𝐺,𝐺  growth rate [m s−1], [1] 𝑟, 𝑟 reaction rate [mol m−3 s−1], [1] 𝑐 concentration [mol m−3] 𝑆 supersaturation [1] 𝑋 island volume fraction [1] 𝜑 solvent volume fraction [1] 𝑀 molar mass [kg mol−1] 𝑘𝑣 volume shape factor [1] 𝜌𝑐𝑟 crystal density [kg m−3] 𝑘𝐺 growth rate coeff. [m s−1] 𝑘𝐽 nucleation rate coeff. [m−3] 𝑔 growth rate exponent [1] 𝑗 nucleation constant [1] 𝐷𝑎𝑐𝑟 Damköhler number [1] REFERENCES 1. Qu Y, Cheng J, Mao ZS, Yang C. A perspective review on mixing effect for modeling and simulation of reactive and antisolvent crystallization processes. Reaction Chemistry & Engineering 6:183-196 (2021) https://doi.org/10.1039/D0RE00223B 2. Trnka J, Maggioni GM, Štěpánek F. Development and Application of a Simplified Non-Ideal Mixing Model for Semi-Batch Crystallization. Computer Aided Chemical Engineering 53:193-198 (2024) https://doi.org/10.1016/B978-0-443-288241.50033-8 3. Baldyga J, Bourne JR. Simplification of micromixing calculations. I. Derivation and application of new model. The Chemical Engineering Journal 42:83-92 (1989) https://doi.org/10.1016/0300-9467(89)85002-6 4. Ståhl M, Rasmuson ÅC. Towards predictive simulation of single feed semibatch reaction crystallization. The Chemical Engineering Journal 64:1559-1576 (2009) https://doi.org/10.1016/j.ces.2008.12.001 5. Cui, P.; Yin, Q.; Gong, J.; Wang, Y.; Hao, H.; Xie, C.; Bao, Y.; Zhang, M.; Hou, B.; Wang, J., Thermodynamic analysis and correlation of solubility of candesartan cilexetil in aqueous solvent mixtures. Fluid Phase Equilibria 337:354362 (2013) https://doi.org/10.1016/j.fluid.2012.09.027 6. Lindenberg C, Mazzotti M. Experimental characterization and multi-scale modeling of mixing in static mixers. Part 2. Effect of viscosity and scale-up. Chemical engineering science 64:42864294 (2009) https://doi.org/10.1016/j.ces.2009.06.067 © 2025 by the authors. Licensed to PSEcommunity.org and PSE Press. This is an open access article under the creative commons CC-BY-SA licensing terms. Credit must be given to creator and adaptations must be shared under the same terms. See https://creativecommons.org/licenses/by-sa/4.0/