Modeling of Continuous Slug Flow Cooling Crystallization towards Pharmaceutical Applications

Abstract

The rising trend towards continuous production in the field of small-scale crystallization has generated many creative concepts for apparatuses for the production of active pharmaceutical ingredients. One of these promising apparatuses is the Slug Flow Crystallizer (SFC), which enables the adjustment of the particle size distribution and the achievement of high yields through its alternating slug flow. To realize and understand the crystallization inside the SFC, high experimental effort has been necessary until now. Therefore, a mechanistic model considering the hydrodynamics of slug flow, the energy and mass balances, and the crystallization phenomena of growth and agglomeration inside the apparatus was developed. Its purpose is to improve the understanding of the process, estimate the effects of operating parameters on target properties, and predict crystallization behavior for different substance systems with minimal experimental effort. Successful modeling was validated with experimental results for the substance system l-alanine/water. Furthermore, the robustness of the model was evaluated, and guidelines were presented, enabling the transfer of the model to new substance systems.

1. Introduction

The trend toward continuous production is becoming increasingly prevalent and researched in the area of small-scale crystallization (<100 mL min⁻¹), which can be used in the field of active pharmaceutical ingredients (APIs). Here, continuous crystallization is a promising strategy for process intensification, offering various advantages compared to the established batch processes. Continuous operation allows production with constant product properties and eliminates batch-to-batch variability. Because the equipment is smaller, it comes with lower facility costs and space requirements [1,2]. Additionally, spatial separation of the crystallization phenomena nucleation and growth is possible and can improve product quality control [3].

In API crystallization, the polymorphic identity, shape, mean size, and size distribution [3,4,5] are the important product properties. To achieve a narrow particle size distribution (PSD), it is necessary to ensure that all particles receive the same treatment. Therefore, a narrow residence time distribution (RTD) of the liquid and solid phases and the avoidance of particle settling are crucial [6,7]. Furthermore, the crystallizer must be designed to promote the desired phenomena and hinder the undesired [1,6]. Nucleation and agglomeration can either be desired or undesired mechanisms, depending on the product requirements and the setup [8,9]. It is important to suppress attrition, promote crystal growth as the preferred mechanism for reducing supersaturation, and ensure a residence time (RT) of several minutes to achieve the desired average product size. Different small-scale crystallizer concepts aiming to solve these challenges have been developed and are divided into mixed suspension mixed product removal (MSMPR) and tubular crystallizers [1,2,8,10].

The tubular crystallizer investigated in this publication is the Slug Flow Crystallizer (SFC). This concept uses a second, immiscible fluid phase to segment the mother liquor into slugs of characteristic length. Due to this segmentation, a plug-flow-like flow profile is achieved at laminar conditions. The interfaces trap the crystals in the slugs and no back-mixing of the liquid or solid phase occurs. Moreover, the interfaces and the wall’s no-slip condition induce so-called Taylor vortices, enhancing mixing, heat transfer, and suspension of the particles [8,11,12].

Important design and operating parameters that need to be selected for optimum SFC operation are the nucleation strategy, design of the mixer contacting the fluids (slug formation zone), choice of a second fluid, cooling strategy, tubing material, and length and inner diameter of the tubing. Besides the operating conditions, these parameters influence the performance of the SFC.

In the literature, inner tubing diameters $d_{\mathrm{i}}$ ranging from 2 mm [13,14] to 4.6 mm [15] are used. Small inner diameters are preferred because the surface and viscous forces are dominant over gravity and inertia forces in this microscale range, which is essential to maintaining a stable slug flow [16]. With increasing tubing diameter, slug coalescence or rupture becomes more probable, and the range of allowable operating conditions is reduced. While a smaller diameter allows a more robust operation, the tubing diameter is limited downwards by the risk of fouling and blockage and the required throughput.

All slugs should have a similar size and shape and must fill the whole cross-section [16]. To segment the mother liquor into slugs of the desired length, either gas (often air) or a second, immiscible liquid is used. Su and Gao [17] state that using a wall material with a high three-phase contact angle leads to convex slugs and the inhibition of migration to adjacent slugs in an air-separated system. An increased contact angle minimized secondary nucleation at the wall, resulting in a narrower PSD, corresponding to the dependency of heterogeneous nucleation on the contact angle. The observation by Su and Gao is further supported by various groups successfully applying air-liquid separated slug flow crystallization [3,16,18,19,20]. Termühlen et al. [20] showed that a convex slug shape is important to achieve narrow RTDs of the solid and liquid phases and to prevent back-mixing and wall film formation. In the absence of a wall film, the velocities of the liquid and gaseous phases are equal [16,20]. Therefore, the total volume flow rate, defined as the sum of the liquid and gaseous volume flow rate, can be used for the calculation of the RT, which is an important parameter in the SFC and in crystallization in general. Long RTs are desired to achieve sufficient growth [9]. The RT in the SFC can be increased by reducing the volume flow rate, increasing the tubing inner diameter, or increasing the tubing length. The slug stability limits the allowable diameter, and mostly the upstream process determines the volume flow rate. Therefore, the length is the only free design parameter that can be chosen to attain the desired RT. However, the pressure drop resulting from longer tubing is a challenge, so until now, a compromise has to be found between design parameters and operating parameters in order to provide a sufficiently long RT.

In the literature, reproducible crystallization in the SFC was proven [13,18,20,21]. However, the experimental data are limited to the investigated substance systems and operating conditions and cannot be extrapolated. To close this gap, process models are useful, generating new insights and increasing the transferability of the experimental results [5]. Generally, the evolution of a PSD over time during a crystallization process can be modeled using a population balance that incorporates all crystallization phenomena (growth, nucleation, agglomeration, and breakage) and the operating conditions [22,23,24]. With a PBE-based crystallization model, the user can predict the resulting PSD for an operating point, estimate kinetic parameters from experimental data, or optimize the process conditions to obtain the desired PSD [25,26]. The PBE calculates the change in number density distribution $n=dN/dL$, which is derived from the number of particles $N$ of a specific size $L$ per unit volume. Equation (1) depicts the general formulation of the PBE. The PBE can be formulated based on the particle volume $V$ or length $L$. The two formulations are equal and can be converted using a form factor [25,26]. For our application, a formulation based on the crystal size is chosen because it is convenient in the consideration of growth.

$$\frac{\partial n}{\partial t}+\frac{\partial \left(Gn\right)}{\partial L}+n\frac{\partial V}{V\partial t}+{\displaystyle \sum }_k\frac{{\dot{V}}_{\mathrm{i}}{n}_{\mathrm{i}}}{V}+D\left(L\right)-B\left(L\right)=0$$

(1)

The first summand $\frac{\partial n}{\partial t}$ computes the change in number density distribution with time, the second summand $\frac{\partial \left(Gn\right)}{\partial L}$ crystal growth, the third $n\frac{\partial V}{V\partial t}$ changes in crystallizer volume, the fourth ${\displaystyle \sum }_k\frac{{\dot{V}}_{\mathrm{i}}{n}_{\mathrm{i}}}{V}$ in- and outgoing streams, and the last death and birth by nucleation, agglomeration, and attrition $D\left(L\right)-B\left(L\right)$.

Figure 1a schematically shows the implementation of a PBE regarding only the crystallization phenomena. The crystallization behavior can be approximated with knowledge of supersaturation and operating conditions at a time t_j. Therewith, the PBE can be integrated to calculate the number density distribution and, thus, the PSD at the following time step t_j+1. With the change in PSD, the mass balance is evaluated. Combined with the energy balance, the supersaturation $S$, and following, the next time step t_j+2 of the PBE can be computed. Therefore, for calculating the PBE, expressions for all observed crystallization phenomena and mass and energy balances are required.

Although modeling has a long history for batch and continuous crystallizers, only a few applications in SFC have been found in the literature until now [4,18,28]. It was experimentally and theoretically shown that an SFC can be modeled as an independent, well-mixed batch crystallizer by observing one slug over the SFC’s length [3,29]. Therefore, no mass transfer occurs between adjacent slugs [18] and every slug receives the same treatment. This assumption of treating a slug as a batch crystallizer is valid if the RTDs of the solid and liquid phases are narrow and equal, which is the case in the SFC investigated in this publication [20]. The time scale of the batch crystallizer can be transformed to the length scale of the SFC using the slug’s velocity. If gas is used as a separating agent, its volume will expand due to the pressure drop decreasing along the flow path and will thus accelerate the liquid slugs [30,31,32,33]. This acceleration enhances heat transfer and suspension but reduces the RT. To the best of the author’s knowledge, all SFC modeling studies in the literature use a constant velocity, neglect the acceleration, and do not calculate the pressure drop.

The three studies for modeling SFC found in the literature follow different approaches with regard to the cooling concept and the crystallization phenomena considered. All studies neglect attrition as it was not observed in their experiments due to the low shear forces and the absence of stirrers in SFC [4,18,28]. Besenhard et al. [28] coupled a PBE with computational fluid dynamics (CFD) to model the crystallization of acetylsalicylic acid from ethanol inside an air-separated SFC. They used CFD to follow the particle trajectories inside a slug. For the PBE, they considered the crystallization phenomena of growth and agglomeration. Besides attrition, primary nucleation is neglected because seed crystals are used. Agglomeration is calculated using an empirical kernel depending on the particle sizes, growth rate, an agglomeration constant, and slug velocity. They considered the velocity because they state that a higher velocity increases the collision frequency and, therefore, the agglomeration rate. The agglomeration constant is fitted to one experiment and validated using two further ones. The first experiment can be described well with the fitted agglomeration rate. However, further experiments performed at lower slug velocities cannot be described accurately. Therefore, this agglomeration kernel cannot be applied for extrapolating to different operating conditions.

Rasche et al. [4] used seed crystals and considered only growth in their model. They focused on modeling the influence of different cooling concepts on the supersaturation profile and crystal growth. They compared calculations with four water baths in a row and multiple counter-current cooling jackets. They state that water baths lead to supersaturation peaks at each inlet, which can induce unwanted secondary nucleation. The use of cooling jackets, however, results in a smooth temperature profile, reducing the maximum supersaturation while providing a similar yield. Based on their modeling, they conclude that an operation with cooling jackets instead of baths is preferred. However, the experimental validation did not take place.

Mozdzierz et al. [18] investigated the crystallization of l-asparagine monohydrate from water in an air-separated SFC. They described a similar setup as Rasche et al. [4] using multiple cooling jackets in a row. They extend the apparatus with a nucleation zone where nuclei are formed continuously by ultrasonication. The model includes growth and nucleation, for the latter primary heterogeneous nucleation through ultrasonication and secondary nucleation are considered. Because they observed neither agglomeration nor attrition in their experiments, these phenomena are neglected in the model. They re-calculated 27 different experiments with varied cooling rates and ultrasonication amplitudes. Most calculated crystal sizes lay in the range of ±20% around the measured value. However, some computations lay outside the ±40% interval and show high deviations. In their experiments, the resulting particle size was underestimated. They attribute this to an overestimation of the influence of ultrasonication on primary nucleation. However, in some of their experimental PSD data, a non-negligible coarse fraction is observable, which is not predicted by the model. This could be caused by agglomeration, which is not incorporated in their model.

Up to now, no model has been published for slug flow crystallization that considers the hydrodynamics of slug flow, namely the influence of pressure drop [30,31,32,33], slug length, and RT decrease through gas expansion. Further, agglomeration cannot be predicted reliably with the existing empirical kernels. All aspects are crucial for the accurate modeling of crystallization and the prediction of product crystals inside the SFC. Hence, this publication aims to derive a model that considers the influences of the slug flow hydrodynamics on the crystallization behavior inside the apparatus. Mechanistic approaches to describe the crystallization phenomena are used to guarantee the transferability of prediction for different operating parameters. In this way, an optimal operating window for the production of high-quality product particles with minimum experimental effort can be determined. Afterward, the robustness of the model is evaluated using sensitivity analysis, and a guideline is presented, enabling the transfer of the model to new substance systems.

2. Overview of Experimental Data for Model Derivation

This section sums up previously conducted and published experiments [20,34] used for model derivation. The experimental procedures are outlined briefly. For a more detailed description of the experiments, the reader is referred to the literature [20,34].

All cooling crystallization experiments were conducted using the setup shown in Figure 2, and the model system l-alanine/water was used. l-alanine was purchased from Evonik Industries AG (Essen, Germany) with a purity ≥99.7%. Ultrapure water with a total organic carbon content ≤3 ppb, purified by a Milli-Q^® Advantage A10 apparatus (Merck KGaA, Darmstadt, Germany), is used as a solvent to prevent primary heterogeneous nucleation. Equation (2) shows the temperature-dependent solubility curve of l-alanine in water regressed by Wohlgemuth and Schembecker [27] and validated by literature data [35,36].

$$c^{*}\left(g_{\mathrm{Ala}}{g}_{\mathrm{sol}}^{-1}\right)=0.11238·exp\left(9.0849·{10}^{-3}·{\vartheta }^{*}\left(°C\right)\right)$$

(2)

In the nucleation zone, l-alanine seed crystals ($w_{\mathrm{solid}}=0.01 g_{\mathrm{solid}}{g}_{\mathrm{sol}}^{-1}$) within the sieve range of 160–200 µm were added in a saturated solution (${\vartheta }^{*}=50 °C, V_{\mathrm{vessel}}=500 \mathrm{mL}, n_{\mathrm{stirrer}}=450 \mathrm{rpm}$). After adjustment of the three-way valve under the vessel, the suspension is transported to the slug formation zone by means of a peristaltic pump (Ismatec Reglo Digital MS-4/12, $d_{\mathrm{i}}=2.29 \mathrm{mm}$, Pharmed). This enables crystal-friendly transport without abrasion for the seed fraction used (see Figure S1). Alternatively, cleaning intervals using water as a cleaning agent can be realized via the setting of the three-way valve.

In the slug formation zone, the suspension is divided into equal-sized slugs via the opposite feed of the segmentation medium (synthetic air, Messer Griesheim, Germany) in a polypropylene T-junction ($d_{\mathrm{i}}=3 \mathrm{mm}$). The gas flow rate is regulated by a needle valve (NV-001-HR, Bronkhorst) and a mass flow meter (El-Flow-Select, Bronkhorst, The Netherlands) based on the measured pressure drop along the entire SFC. The slug lengths were measured by image analysis of a video taken after slug formation [37].

In the growth zone, the slugs pass along a fluorinated-ethylene-propylene (FEP) tubing ($d_{\mathrm{i}, \mathrm{tubing}}=3.18 \mathrm{mm}$) with a length of up to 26.5 m. The tubing material was chosen for its hydrophobic material properties and allows for convex slugs, resulting in a narrow residence time distribution of the liquid and solid phases [20,38]. In order to determine the influence of the suspension on the target parameters, the residence time is varied by adjusting the volume flow rates (suspension and gas volume flow). Due to the expansion of the gas bubbles, the actual residence time in the apparatus is shorter than the set residence time. The actual residence time is measured by observing a single defined gas bubble along the tubing in steady-state operation.

The process medium (PM) temperature in the slugs is controlled by a tube-in-tube co-current cooling concept. The tempering medium (TM) enters the cooling jacket at the same temperature as the PM. This TM is cooled by heat loss to the environment and, therefore, cools the PM. The temperature profile of the PM can be adjusted using the TM flow rate. At high TM flow rates, the heat loss to the environment has limited influence, and the cooling rate is low, resulting in a high outlet temperature of the PM. Contrarily, at low TM flow rates, the heat loss has a strong influence, resulting in low outlet temperatures and high cooling rates. The PM temperature profile can be measured non-invasively by measuring the TM temperature [20]. The cooling rate is, therefore, adjusted by setting the final temperature of the TM. The measured slug lengths, pressure drops, and operating conditions are given in Tables S1 and S2 in the Supporting Information for all experiments used.

At the end of the crystallizer, the flow pattern and the goodness of suspension inside the slugs are monitored via video and analyzed using image analysis in order to evaluate the slug and gas bubble length and to identify the suspension state of particles inside the slugs [7,37,39]. After evaluating the video, the suspension state is output via dimensionless key figures (χ_50,H and χ_50,V) to describe the horizontal and vertical distribution of the particles. Afterward, the slugs were collected at the end of the crystallizer, the concentration in the liquid phase was determined gravimetrically, and the PSD was determined by means of a dynamic image sensor system (QICPIC, Sympatec, LIXELL dispersion system, 1024 × 1024 resolution, module M6). For evaluating the PSD, the median particle size $d_{50,3}$, the width of particle size distribution $d_{90-10,3}$, and the agglomeration degree $Ag$ are calculated with an in-house MATLAB script [40,41].

In the reviewed cooling crystallization experiments, the influences of the slug length, the RT, the supersaturation, and the particle suspension and mixing inside slugs on the product quality (PSD, yield) were evaluated. Since all process parameters are interdependent, it is not possible to change each parameter separately while keeping the other parameters constant, e.g., a change in flow velocity directly affects RT, particle suspension, slug length, and cooling rate. The complex interactions of design and operating parameters and their influences on product quality attributes for the setup used in the reviewed experiments are illustrated in Figure 3.

3. Model Derivation for the Slug Flow Crystallizer

The complex relationships, shown in Figure 3, are now supposed to be incorporated into the model so that the design and operation of the SFC can be optimized and the product quality attributes can be predicted. For the model, only the growth zone is considered due to the addition of seed crystals in the nucleation zone. Since, inside the growth zone, neither nucleation nor attrition was observed in the setup and experiments reviewed [20,34], the evolution of the particle number density can be calculated considering only growth and agglomeration phenomena. No in- and outgoing streams must be considered because one slug is considered a batch crystallizer, and no mass transfer or exchange of fluid elements is assumed to take place between adjacent slugs. Furthermore, the volume of the slug is assumed constant. Therefore, the PBE (Equation (1)) simplifies to Equation (3).

$$\frac{\partial n}{\partial t}+\frac{\partial \left(Gn\right)}{\partial L}+D\left(L\right)-B\left(L\right)=0$$

(3)

The resulting PBE is a partial differential equation consisting of two partial differentiations, one with respect to time $\frac{\partial n}{\partial t}$, and one to crystal size $\frac{\partial \left(Gn\right)}{\partial L}$. The PBE is formulated according to the time scale. However, the temperature profile of the SFC is measured over the length scale of the SFC. Therefore, the PBE must be converted to the length scale to combine it with the energy balance. In the SFC, slugs accelerate due to gas expansion, which must be considered in the conversion. For a non-constant velocity, the time scale $t$ can be converted to the length scale $z$ using the velocity $v$ and its spatial derivative $\partial v/\partial z$ [42]. Equations (4)–(7) depict the conversion of an arbitrary function $f$. Using the conversion introduced in Equation (7), the PBE can be calculated over the SFC’s length. For a constant velocity, its spatial derivative is zero, and Equation (7) is simplified.

$$\frac{\partial f}{\partial z}=\frac{\partial f}{\partial t}·\frac{\partial t}{\partial z}$$

(4)

$$t=\frac{z}{v}$$

(5)

$$\frac{\partial t}{\partial z}=\frac{\partial }{\partial z}\left(\frac{z}{v\left(z\right)}\right)=\frac{v-z\frac{\partial v}{\partial z}}{v^2}$$

(6)

$$\frac{\partial f}{\partial z}=\frac{v-z\frac{\partial v}{\partial z}}{v^2}·\frac{\partial f}{\partial t}$$

(7)

In the following, the population balance (Equation (3)) is set up, and the steps and influences shown in Figure 1b are taken into account. Prior to solving the population balance, it is essential to analyze and describe the hydrodynamics of slug flow. This is crucial in order to accurately determine and comprehend the location of each slug/batch within the tubing. Afterward, expressions for mass balance and supersaturation, energy balances, and crystallization phenomena (growth and agglomeration) are derived.

3.1. Solving the Population Balance Equation

The reviewed crystallization experiments were conducted in the SFC using seed crystals, and the seed crystal data were given to the model as an initial condition. Furthermore, the PBE is solved by discretization. The investigated particle size domain $L$ is geometrically discretized into $N$ classes. The discretization scheme is selected such that the volume of one particle in a class $V_{\mathrm{i}}$ is doubled compared to the volume of one particle in the prior class $V_{\mathrm{i}-1} \left(V_{\mathrm{i}}=2V_{\mathrm{i}-1}\right)$. Therefore, the characteristic particle size in interval $i$ is $L_{\mathrm{i}}=\sqrt[3]{2}{L}_{\mathrm{i}-1}$, resulting in non-constant differences between the classes $\Delta L_{\mathrm{i}}=L_{\mathrm{i}}-L_{\mathrm{i}-1}$. This classification is chosen because it enables us to cover a great size range. Furthermore, it corresponds to the particle size classification of the dynamic image analysis sensor QICPIC (Sympatec GmbH, Clausthal-Zellerfeld, Germany), the cam sizer used to determine the PSD in the reviewed experiments. Additionally, it is advantageous to consider agglomeration because two particles of class $i$ can agglomerate to one particle of class $i+1$ by conserving the volume. [22]

3.2. Modeling Slug Flow Hydrodynamics

The slug flow hydrodynamics must be modeled to evaluate the slug’s acceleration and the actual RT, which is reduced compared to the hydrodynamic RT calculated for the inlet conditions and, therefore, the following time available for crystal growth. The acceleration results from the gas expansion, which depends on solvent evaporation, pressure drop, and cooling. The gas expansion $V_2/V_1$ can be approximated using the ideal gas law considering changes in the pressure $p$, amount of substance $n$, and temperature $T$. Their influence is evaluated by investigating a typical setup of the SFC with operating conditions given in Table S3.

Synthetic air, consisting only of nitrogen and oxygen supplied by a gas cylinder, is used as a segmentation medium. Because the air is water-free, a certain amount will evaporate along the length of the SFC, limited by the temperature-dependent vapor pressure. For a worst-case estimation, slugs consisting of water only are calculated. The influence of l-alanine on evaporation can be neglected because the activity coefficients of l-alanine/water solutions are near unity [43]. Furthermore, the influence of crystallization on the gas bubbles is neglected. The volume increase in gas bubbles is calculated assuming that the air at the outlet is saturated. Therefore, Raoult’s law can be applied to find the outlet mole fraction of water $y_2$ using Equation (8) and following the change in the amount of substance $n_2/n_1$ using Equation (9). The gas expansion through temperature change and pressure drop can be calculated directly from experimental data. The determined volume change of each individual effect is given in

Table 1. Gas expansion $V_2/V_1$ of the individual effects calculated with the ideal gas law for the operating conditions given in Table S3.

	Evaporation	Cooling	Pressure Drop
V2/V1	1.04	0.94	2.00

.

$$y_2=\frac{p_{0,2}^{LV}}{p_2}$$

(8)

$$\frac{n_2}{n_1}=\frac{1}{1-y_2}$$

(9)

Evaporation results in a small volume increase, cooling in a small decrease. Pressure drops in the range of 1 bar, however, result in a doubling of the slug volume, clearly dominating gas expansion. Following this, the pressure drop is the main factor influencing the gas expansion. It depends on the hydrodynamic pressure drop of both phases and the pressure drop at the interfaces [44]. Therefore, the number of interfaces in the SFC is required, which can be calculated with the number and length of a unit cell consisting of a gas bubble and a liquid slug.

After the pressure drop, the slug length plays a significant role in affecting mixing, particle suspension, and heat transfer [37,44,45]. Several factors affect the length of the slug, including the proportion of liquid to total volume flow rate (liquid hold-up), the velocity and properties of the phases, the method of phase supply, the bubble detachment mechanism, as well as the size, geometry, material, and three-phase contact angle of the mixer [37,46,47]. Many parameters significantly influence the slug formation and the slug’s length and, therefore, a correlation for calculating the slug length must be derived for a similar setup and operating conditions. Etminan et al. [48] reviewed various slug flow correlations for flow patterns, flow regime transitions, slug lengths, and pressure drops. The correlation given in Equation (10) can be used to estimate the slug length $L_{slug}$ for the SFC investigated in this publication since the operating conditions lie in the range investigated by Qian and Lawal [47] (Table S4).

$$\frac{L_{\mathrm{slug}}}{d_{\mathrm{SFC}}}=c_1·{\left(1-{\epsilon }_{\mathrm{L}}\right)}^{c_2}·{\epsilon }_{\mathrm{L}}^{c_3}·R{e}^{c_4}·C{a}^{c_5}$$

(10)

Using the dimensionless numbers liquid hold-up ${\epsilon }_{\mathrm{L}}=Q_{\mathrm{L}}/Q_{\mathrm{tot}}$, Reynolds number $Re=\left({\rho }_{\mathrm{L}}{v}_{\mathrm{slug}}{d}_{\mathrm{SFC}}\right)/{\eta }_{\mathrm{L}}$, and Capillary number $Ca=\left({\eta }_{\mathrm{L}}{v}_{\mathrm{slug}}\right)/{\sigma }_{\mathrm{L}}$, the correlation considers the influence of operating conditions (phase ratio and slug velocity), tubing geometry (tubing inner diameter), and properties of the liquid phase (density, viscosity, and surface tension). Those are weighted by Qian and Lawal with the correlation parameters $c_i$ given in

Table 2. Comparison of the correlation parameters by Qian and Lawal [47] and the fit to experimental data reviewed in this publication.

	${\mathit{c}}_1$ /-	${\mathit{c}}_2$ /-	${\mathit{c}}_3$ /-	${\mathit{c}}_4$ /-	${\mathit{c}}_5$ /-
Qian and Lawal [47]	1.637	−0.893	−0.050	−0.075	−0.0687
Fit to exp. data	1.969	−1.102	−0.035	−0.176	−0.0605

.

To validate the applicability of the correlation on the investigated SFC, experimental data measured and published by Termühlen et al. [37,39] and Kufner et al. [49] are reviewed. Experiments measuring the mean slug length were executed with suspension and water at varying flow rates, liquid hold-ups, and solid mass fractions (Table S1). The slug length is approximated using the measured volume flow rate at the inlet. Figure 4a shows the parity plot for the measurements and calculations according to the correlation by Qian and Lawal [47] (orange triangles). Besides two outliners, the slug length is overestimated by their correlation. Most calculated slug lengths lie outside the ±20% confidence interval with a mean absolute error of 26.63%. Therefore, the correlation parameters $c_i$ are fitted to the experimental data to reduce the mean error of the slug length (Figure 4a, blue squares). The fitted data spread around the parity, and the mean absolute error is reduced to 7.30%. With the fitted correlation, the slug length can be predicted reliably. The fitted correlation parameters are given in Table 2, too.

Based on the slug length calculation, the gas bubble lengths $L_{\mathrm{G}}$ and the length of a unit cell ($L_{\mathrm{UC}}=L_{\mathrm{slug}}+L_{\mathrm{G}}$) can also be determined, which is needed in the next step to identify the pressure drop in the apparatus. According to Kashid and Agar [44], the pressure drop in slug flow can be calculated using Equation (11) as the sum of the single-phase hydrodynamic pressure drops of the liquid $\mathsf{\Delta }{p}_{\mathrm{hyd},\mathrm{L}}$ and the gaseous $\mathsf{\Delta }{p}_{\mathrm{hyd},\mathrm{G}}$, which are weighted by the liquid hold-up ${\epsilon }_{\mathrm{L}}$ and the gaseous hold-up $\left(1-{\epsilon }_{\mathrm{L}}\right)$, respectively, and the multiphase term $\mathsf{\Delta }{p}_{mp}$ resulting from the interactions of the two phases. The single-phase hydrodynamic pressure drops $\mathsf{\Delta }{p}_{\mathrm{hyd},\mathrm{i}}$ can be calculated using the Hagen–Poiseuille equation depending on the slug velocity $v_{\mathrm{slug}}$, the phase’s viscosities ${\eta }_{\mathrm{i}}$, and the tubing’s inner diameter $d_{\mathrm{SFC}}$ [44,48]. The multiphase pressure drop $\Delta p_{\mathrm{mp}}$ can be calculated depending on the surface tension ${\sigma }_L$, the inner tubing diameter $d_{\mathrm{SFC}}$, and the three-phase contact angle $\theta$ [44].

In the case of crystallization, the influence of the particles on the pressure drop of the liquid must be considered, too. This is performed by estimating the viscosity of the suspension ${\eta }_{\mathrm{susp}}$. For low particle volume fractions (≤10%), viscosity is a linear function of the volume fraction [50,51,52,53] valid for the evaluated experiments. Mueller et al. [50] showed that the suspensions viscosity ${\eta }_{\mathrm{susp}}$ can be calculated with Equation (12) derived by Einstein [54,55] depending on the solution’s dynamic viscosity ${\eta }_{\mathrm{sol}}$ and the particle volume fraction φ.

$$\begin{gathered} \mathsf{\Delta }{p}_{\mathrm{tot}}=\mathsf{\Delta }{p}_{\mathrm{hyd}}+\mathsf{\Delta }{p}_{\mathrm{mp}}={\epsilon }_{\mathrm{L}}·\mathsf{\Delta }{p}_{\mathrm{hyd},\mathrm{L}}+\left(1-{\epsilon }_{\mathrm{L}}\right)·\mathsf{\Delta }{p}_{\mathrm{hyd},\mathrm{G}}+\mathsf{\Delta }{p}_{\mathrm{mp}} \\ ={\epsilon }_{\mathrm{L}}·\frac{32·{\eta }_{\mathrm{susp}}·v_{\mathrm{slug}}}{d_{\mathrm{SFC}}^2}+\left(1-{\epsilon }_{\mathrm{L}}\right)·\frac{32·{\eta }_{\mathrm{G}}·v_{\mathrm{slug}}}{d_{\mathrm{SFC}}^2}+N_{\mathrm{UC}}·\frac{4\sigma }{d_{\mathrm{SFC}}}\cos \left(\theta \right) \end{gathered}$$

(11)

$${\eta }_{\mathrm{susp}}={\eta }_{\mathrm{sol}}·\left(1+2.5\mathsf{\phi }\right)$$

(12)

For calculating $\Delta p_{\mathrm{mp}}$, the pressure drop along one interface must be multiplied by the number of interfaces. The number of interfaces $N_{\mathrm{UC}}$ can be estimated using the length of a unit cell $L_{\mathrm{UC}}$ consisting of a liquid slug $L_{\mathrm{slug}}$ and a gas bubble $L_{\mathrm{G}}$ (unit cell visualized in Figure 5) and the fraction of the length of the SFC, according to Equation (13). Each unit cell has two interfaces.

$$N_{\mathrm{UC}}=2\frac{L_{\mathrm{tubing}}}{L_{\mathrm{UC}}}-1$$

(13)

In calculating the number of unit cells, the gas expansion is neglected because the expansion can only be computed with the knowledge of the pressure drop. The expansion reduces the number of interfaces and, therefore, the pressure drop. In contrast, the slug velocity would increase, resulting in a higher hydrodynamic pressure drop. For the first estimation, both effects are neglected. Figure 4b shows the parity plot for the total pressure drop $\mathsf{\Delta }{p}_{\mathrm{tot}}$. The operating conditions are given in Table S2. Most of the data lie inside the ±20% confidential interval, and the mean absolute prediction error is 8.0%. However, some data points lie outside the confidence interval. Remarkably, especially experiments with a high ${\epsilon }_{\mathrm{L}}$ are underestimated, while experiments with low ${\epsilon }_{\mathrm{L}}$ are predicted well. Since these experiments have a similar number of interfaces but a lower ${\epsilon }_{\mathrm{L}}$, the underestimation of $\mathsf{\Delta }{p}_{\mathrm{tot}}$ can be attributed to underestimating the suspension hydrodynamic pressure drop. The reliable prediction of most data points validates the assumption of neglecting the two counteracting effects of gas expansion on pressure drop. With knowledge of slug length and pressure drop, gas expansion and the actual RT can be calculated at an arbitrary position inside the SFC using the ideal gas equation (Equation (S1)).

Figure 4c shows the parity plot for the hydrodynamic and actual RTs (calculated for the inlet conditions) [39,49]. All experiments used for the validation of the expressions for the actual RT calculation cover the whole operating range of the SFC. The liquid hold-up ${\epsilon }_L$ varied between 0.33 and 0.66, total flow rate $Q_{\mathrm{tot}}$ between 20 and 60 mL min⁻¹, and the SFC tubing length $L_{\mathrm{tubing}}$ between 7.5 and 37 m. The experimental conditions and measured and calculated RTs are given in the appendix in Table S5. Figure 4c shows that the hydrodynamic RT overpredicts the measured RT and some points lie outside the confidence interval. Using the correlations for the slug length and pressure drop described before, the actual RT can be predicted well for all operating conditions. The calculated actual RTs spread around the parity line and all points lie inside the ±20% confidential interval.

3.3. Mass Balance and Supersaturation

The modeling of crystallization phenomena is connected by solving the mass balance, where the actual and saturated concentration must be known to identify the supersaturation at every time step. The saturation concentration of component i (here: l-alanine) $c_i^{*}$ is calculated according to Equation (2) for the respective temperature along the SFC tubing determined by the derived energy balance. The actual concentration of component i $c_i$ is calculated using the mass balance given in Equation (14) as the ratio of the dissolved mass of component i $m_{\mathrm{i},\mathrm{sol}}$ to the solution mass $m_{\mathrm{sol}}$. Both masses are time-dependent and decreased by crystallization. The dissolved mass is calculated as the difference in the total mass of component i and crystal mass $m_{\mathrm{c}}$. The solution mass $m_{\mathrm{sol}}$ is the difference in the initial solution mass $m_{sol,0}$ and the crystallized mass ($m_{\mathrm{c}}\left(t\right)-m_{\mathrm{seed}}$). The crystal mass in one slug $m_{\mathrm{c}}$ can be calculated with Equation (15) depending on the slug volume $V_{\mathrm{slug}}$, crystal volume shape factor $k_{\mathrm{V}}$, discrete particle size $L_{\mathrm{i}}$, and discrete particle number distribution $n_{\mathrm{i}}$. Furthermore, the crystal mass derivative is required to calculate the crystallization heat (Equation (16)). The crystal number distribution is the only time-dependent quantity, and its derivative is known from the PBE.

$$c_{\mathrm{i}}=\frac{m_{\mathrm{i},\mathrm{sol}}\left(t\right)}{m_{\mathrm{sol}}\left(t\right)}=\frac{m_{\mathrm{i},\mathrm{sol},0}-\left(m_{\mathrm{c}}\left(t\right)-m_{\mathrm{seed}}\right)}{m_{\mathrm{sol},0}-\left(m_{\mathrm{c}}\left(t\right)-m_{\mathrm{seed}}\right)}$$

(14)

$$m_{\mathrm{c}}=V_{\mathrm{slug}}·{\rho }_{\mathrm{c}}·k_{\mathrm{V}}·{\displaystyle \sum }_{i=1}^N{L}_{\mathrm{i}}^3·n_{\mathrm{i}}$$

(15)

$$\frac{\partial m_{\mathrm{c}}}{\partial t}=V_{\mathrm{slug}}·{\rho }_{\mathrm{c}}·k_{\mathrm{V}}·{\displaystyle \sum }_{i=1}^N{L}_{\mathrm{i}}^3·\frac{\partial n_{\mathrm{i}}}{\partial t}$$

(16)

Afterwards, the supersaturation can be calculated according to Equation (17) depending on the actual concentration $c_{\mathrm{i}}$ and the saturation concentration $c_{\mathrm{i}}^{*}$, or expressed depending on the relative supersaturation $\sigma$.

$$S=\frac{c_{\mathrm{i}}}{c_{\mathrm{i}}^{*}}=\mathsf{\sigma }+1$$

(17)

3.4. Energy Balances

To calculate the saturation concentration at every time step in the PBE, the energy balances within the outer jacket for the TM and for the unit cell within the tubing for the PM have to be solved. In order to set up the energy balances, the time-dependent, batch-wise processes inside a unit cell traveling through the SFC must be coupled to the steady-state, length-dependent processes occurring in the TM. Figure 5a visualizes the SFC tube with one unit cell consisting of a gas bubble and a liquid slug, the TM, and its interaction with the environment. The TM is modeled by segmenting the tube into infinitely small elements and calculating the energy balance in these elements. One segment is shown by the red dotted rectangles in Figure 5a. The TM is modeled respective to the spatial length scale resulting from this segmentation. Its energy balance is determined by the convective heat flows in and out of the segment ${\dot{Q}}_{\mathrm{conv},\mathrm{i}}$, and heat transfer rates with the unit cell ${\dot{Q}}_{\mathrm{UC},\mathrm{TM}}$ and the environment ${\dot{Q}}_{\mathrm{TM},\mathrm{ENV}}$. The steady-state energy balance of the TM can be calculated by Equation (18) [20].

In the energy balance, the oscillatory behavior of the heat transfer rate from the unit cell to the TM must be considered [56]. At one time, the TM at position $z_{\mathrm{i}}$ is in contact with the PM, and at the next time with the segmenting gas, resulting in oscillatory heat transfer rates. With Equation (19), the heat transfer from the unit cell ${\dot{Q}}_{\mathrm{UC},\mathrm{TM}}$ is calculated as the weighted sum of the individual heat transfer rates to incorporate this oscillatory behavior. Heat transfer from the PM to the TM (${\dot{Q}}_{\mathrm{PM},\mathrm{TM}}$) is weighted with the liquid hold-up ${\epsilon }_{\mathrm{L}}$, and heat transfer from the gaseous phase to the TM (${\dot{Q}}_{\mathrm{G},\mathrm{TM}}$) with the gaseous hold-up $(1-{\epsilon }_{\mathrm{L}})$. Detailed expressions for ${\dot{Q}}_{\mathrm{conv},\mathrm{i}}$,${\dot{Q}}_{\mathrm{UC},\mathrm{TM}}$,${\dot{Q}}_{\mathrm{TM},\mathrm{ENV}}$ are given in the Supporting Information Equations (S2)–(S4).

The temperature of the fluid entering the segment is known from the previous segment. The temperature of the fluid leaving the segment is unknown. Therefore, the difference in the heat flow entering and leaving the segment is approximated by the first coefficient of the Taylor series, using the first spatial derivative of the TM’s temperature $T_{\mathrm{TM}}$. The spatial temperature profile of the TM can be calculated by solving Equation (20), where $c_{\mathrm{p},\mathrm{TM}}$ represents the temperature-dependent heat capacity and ${\dot{m}}_{\mathrm{TM}}$ the mass flow of the TM.

$$\frac{d{Q}_{TM}}{dt}=0={\dot{Q}}_{conv,in}-{\dot{Q}}_{conv,out}+{\dot{Q}}_{UC,TM}-{\dot{Q}}_{TM,ENV}$$

(18)

$${\dot{Q}}_{UC,TM}={\epsilon }_L·{\dot{Q}}_{PM,TM}+(1-{\epsilon }_L)·{\dot{Q}}_{G,TM}$$

(19)

$$\frac{\partial T_{\mathrm{TM}}}{\partial z}=\frac{{\dot{Q}}_{\mathrm{UC},\mathrm{TM}}-{\dot{Q}}_{\mathrm{TM},\mathrm{ENV}}}{{\dot{m}}_{\mathrm{TM}}·\left(T_{\mathrm{TM}}\frac{\partial c_{\mathrm{p},\mathrm{TM}}}{\partial T_{\mathrm{TM}}}+c_{\mathrm{p},\mathrm{TM}}\right)·\Delta z}$$

(20)

Figure 5b shows the energy balance of the PM considering one unit cell. The interactions between the PM and the segmenting gas (${\dot{Q}}_{\mathrm{G},\mathrm{PM}})$, their heat flows to the TM (${\dot{Q}}_{\mathrm{G},\mathrm{TM}}$ and ${\dot{Q}}_{\mathrm{PM},\mathrm{TM}}$), and the heat of crystallization ${\dot{Q}}_{\mathrm{c}}$ have to be considered (Equations (S5) and (S6)). Furthermore, the solvent evaporation ${\dot{Q}}_{\mathrm{EVP}}$ has to be evaluated and can be calculated with Equation (S7). ${\dot{Q}}_{\mathrm{G},\mathrm{PM}}$ can be neglected because both phases enter the SFC at the same temperature. In order to identify the most important magnitudes on energy balance, the heat flows are calculated for the cooling of the PM from 50 °C to 30 °C. The magnitudes of the individual contributions to the total heat flow are listed in

Table 3. The calculated values for the individual heat flow for a total energy balance of a unit cell over the SFC for cooling of the PM (l-alanine/water) from ${\vartheta }_{\mathrm{in}}$ = 50 °C to ${\vartheta }_{\mathrm{out}}$ = 30 °C, a total volume flow rate of $Q_{\mathrm{tot}}$ = 20 mL min⁻¹, and a liquid hold-up ${\epsilon }_{\mathrm{L}}$ of 0.5.

Effect	${\dot{\mathit{Q}}}_{\mathit{P}\mathit{M},\mathit{T}\mathit{M}}$	${\dot{\mathit{Q}}}_{\mathit{G},\mathit{T}\mathit{M}}$	${\dot{\mathit{Q}}}_{\mathit{c}}$	${\dot{\mathit{Q}}}_{\mathit{E}\mathit{V}\mathit{P}}$
Magnitude/W	13.32	0.0037	0.0038	0.0238

.

Table 3 shows that the energy balance of a unit cell is dominated by the cooling of the PM. The gaseous phase’s heat flow is magnitudes lower due to its lower density and heat capacity. The influence of crystallization is low due to the low crystallization enthalpy for the substance system l-alanine/water. Evaporation is negligible despite the high evaporation enthalpy because the evaporating mass is low. Therefore, the influences of evaporation and cooling of the gaseous phase can be neglected, simplifying the energy balance of a unit cell to Equation (21). The heat released by crystallization is considered despite its low influence because it can be computed at no further cost since the crystallizing mass is known from the mass balance.

Because the slugs are well mixed, it is assumed that the temperature in the slug is uniform, and the slug can be modeled as one isothermal unit. However, the temperature of the TM decreases along the slug’s length continuously (see Figure 5). Therefore, the slug is, at one time, in contact with TM at different temperatures, resulting in a decreasing heat transfer rate over the slug length. Therefore, the heat transfer must be integrated over the slug’s length at every time step (Equations (S8) and (S9)) and coupled with the length-dependent steady-state profile of the TM (Equations (S10) and (S11)). The temperature profile of PM can be calculated by Equation (22), including the slug’s acceleration along the tubing.

$$\frac{d{Q}_{PM}}{dt}={\dot{Q}}_c-{\dot{Q}}_{PM,TM}$$

(21)

$$\frac{\partial T_{\mathrm{PM}}}{\partial z}=\frac{1}{v}·\frac{\partial T_{\mathrm{PM}}}{\partial t}-\frac{z}{v^2}·\frac{\partial v}{\partial z}\frac{\partial T_{\mathrm{PM}}}{\partial z}$$

(22)

The energy balances derived above are validated by using experimental data measured by Termühlen et al. [20,39]. Figure 6 shows these initially calculated temperature profiles for TM and PM (orange and blue dashed lines) with experimental data (crosses) for an SFC tubing length of $L_{\mathrm{tubing}}$ = 7.5 m, a TM flow rate of $Q_{\mathrm{TM}}$ = 37 mL min⁻¹, and a PM flow rate of $Q_{\mathrm{PM}}$= 10 mL min⁻¹. The low TM flow rate results in a high influence of heat loss due to the low heat capacity in the experiments. Therefore, steep cooling and a low outlet temperature are achieved for a short tubing length.

It is shown that the initially calculated temperature profiles (dashed lines) deviate strongly from the experimental data. The modeled temperatures are higher than the measured ones. Obviously, the calculated heat transfer rates are too low. Since the only heat sink in the system is the heat loss from the TM to the environment, this heat transfer rate determines the temperature profile of the TM and PM. Therefore, the key to improving the model prediction is the better estimation of the heat loss to the environment. The transfer from the outer wall to the environment is the highest contribution to the heat transfer resistance (compare Table S6). To adjust the temperature profile to the experimental data, this heat transfer resistance must be reduced because it is the rate-limiting step.

The convective heat transfer coefficient from the wall to the environment is calculated for one horizontal tube in a resting environment. The tubing of the SFC, however, is vertically coiled (compare Figure S2). This vertical alignment influences the heat transfer by natural convection in two ways [57,58,59,60]. On the one hand, the lower coil locally increases the temperature of the environmental air, reducing its density. This heated air rises and is in contact with the coil above, resulting in an increased air temperature and lower heat transfer rate. On the other hand, the same effect increases convection compared to pure natural convection at a single horizontal pipe, increasing the heat transfer coefficient. These two effects counteract. The dominant factor depends on the operating conditions and the setup. In the case of the SFC used here, the tubing is vertically coiled, inducing a chimney effect inside the coil. This further increases the convection and the heat transfer coefficient. Therefore, the second effect dominates, resulting in an increased heat transfer rate to the environment. To account for this effect, the calculated heat transfer coefficient to the environment must be enhanced. Since all SFC configurations investigated are geometrically similar, the resulting enhancement factor is assumed to be a constant value. The enhancement factor was fitted to all available temperature profile data. The solid lines in Figure 6 depict the TM and PM temperature profiles calculated with the increased convective heat transfer rate (solid lines), which coincide well with the experimental data.

3.5. Modeling Crystal Growth

According to the PBE (Equation (3)), the partial differential equation has to be solved and it is discretized on the particle size grid introduced above. Therefore, the growth term in the PBE is solved using a high-resolution semi-discrete finite volume scheme developed by Koren [61] and adapted for crystallization by Qamar et al. [62,63,64,65]. This method was already successfully implemented in the modeling of l-alanine/water-cooling crystallization [27,66]. For a detailed description of this solution method, the reader is referred to the literature [61,62,63,64,65]. The expression of the growth rate $G$ is determined using the Burton, Cabrera, and Frank model [67], depending on the relative supersaturation $\sigma$ and the two model parameters $A_{\mathrm{BCF}}$ and $B_{\mathrm{BCF}}$ (Equation (23)). Furthermore, size-independent growth is assumed ($G\ne f\left(L\right)$).

$$G_{\mathrm{BCF}}=A_{\mathrm{BCF}}·{\sigma }^2·\tanh \frac{B_{\mathrm{BCF}}}{\sigma }$$

(23)

The growth rate constants $A_{\mathrm{BCF}}$ and $B_{\mathrm{BCF}}$ are derived from experimental data at operating conditions similar to the SFCs. Experimental data measured and published by Steenweg et al. [34] are evaluated. They performed batch cooling crystallization of l-alanine from an aqueous solution using seed crystals ($w_{\mathrm{seed}}=0.01{ \mathrm{g}}_{\mathrm{Ala}}{\mathrm{g}}_{\mathrm{sol}}^{-1}$) in a temperature range from $\vartheta =50 °\mathrm{C}$ to $30 °\mathrm{C}$ with a cooling rate of $\kappa =1.8 {\mathrm{K} \min }^{-1}$. These operating conditions match the SFC’s for the low cooling rate (compare Table S7). Figure 7a shows the calculated concentration profiles of two batch experiments using the fitted rate constants (solid lines) compared to rate constants from the literature [66] (dotted lines).

Using literature growth rate constants, higher concentrations compared to the experimental data are modeled. Therefore, the growth rate is underestimated. This coincides with the results of the calculation of the SFC. By fitting the rate constants, the experimental data can be reproduced well. Figure 7b shows the parity plot comparing the calculated concentration reduction in the SFC with the measured concentrations using the literature and the fitted constants for growth rate determination. Better model accuracy is achieved using the fitted growth rate, leading to a reduction in the mean absolute error to 10 %. Because the calculated values spread around the parity line, the remaining error is not systematic and can be attributed to the propagation of model and parameter uncertainties and experimental inaccuracies. This method confirms the finding that the growth constants must be determined individually for each operating point/operating range in order to obtain the most reliable result possible for modeling crystallization [27]. In the following, the growth rate constants fitted in this publication will be used in the SFC model.

3.6. Modeling Agglomeration

Agglomeration is modeled using birth $B\left(L\right)$ and death rates $D\left(L\right)$. In the case of a continuous PBE, birth and death are calculated by solving the integral terms given in Equations (S12) and (S13) in the Supporting Information [25]. Because the PBE is evaluated at discrete particle sizes only, the integrals must be approximated at these particle sizes. Therefore, it must be considered that the size of the formed agglomerates must lie on the grid points.

The birth $B_{\mathrm{Agg}}\left(L_{\mathrm{i}}\right)$ (Equation (24)) and death rate $D_{\mathrm{Agg}}\left(L_{\mathrm{i}}\right)$ (Equation (25)) according to Hounslow et al. [22] were used. For particles to leave their class $i$ by agglomeration, the size of the formed agglomerate must be the size of class $i+1$ or higher. Agglomeration processes not leaving the class cannot be accounted for in this discretization scheme and are neglected.

$$B_{\mathrm{Agg}}\left(L_{\mathrm{i}}\right)=\frac{1}{\Delta L_{\mathrm{i}}}\left(\frac{1}{2}{\beta }_{\mathrm{i}-1,\mathrm{i}-1}·n_{\mathrm{i}-1}^2·\Delta L_{\mathrm{i}-1}^2+n_{\mathrm{i}-1}·\Delta L_{\mathrm{i}-1}·{\displaystyle \sum }_{j=1}^{i-2}{2}^{\mathrm{j}-\mathrm{i}+1}{\beta }_{\mathrm{i}-1,\mathrm{j}}·n_{\mathrm{j}}·\Delta L_{\mathrm{j}}\right)$$

(24)

$$D_{Agg}\left(L_i\right)=\frac{1}{\Delta L_i}\left(n_i·\Delta L_i·{{\displaystyle \sum }}_{j=1}^{i-1}{\beta }_{i,j}·n_j·\Delta L_j+n_i·\Delta L_i·{{\displaystyle \sum }}_{j=1}^{\infty }{\beta }_{i,j}·n_j·\Delta L_j\right)$$

(25)

The birth and death rates depend on an agglomeration kernel $\beta$. In this publication, a mechanistic approach for the agglomeration kernel was applied, describing the occurring physical phenomena mathematically, which should result in a higher transferability to different operating conditions. Therefore, the approach of Mumtaz and Hounslow [68,69], which was improved by Hounslow and Pitt [70,71], was applied. The kernel ${\beta }_{\mathrm{i},\mathrm{j}}$ is calculated with Equation (26) using a collision frequency ${\beta }_{coll}$ and an agglomeration efficiency $\mathit{\Psi }$. The collision frequency determines the number of collisions between two particles of class $i$ and $j$ and depends on the operating conditions and flow regime. The efficiency determines how many collisions lead to agglomeration instead of rupture. The SFC is operated in a laminar flow. In the investigated operating conditions, the particles are not homogeneously suspended over the whole slug, resulting in an increased collision frequency and agglomeration tendency [20]. Therefore, collision is described considering collisions in a laminar flow field ${\beta }_{\mathrm{coll},\mathrm{LS}}$ and through differential settling ${\beta }_{\mathrm{coll},\mathrm{DS}}$. Because both processes are independent of each other, the total collision frequency is calculated using Equation (27) as the sum of the individual ones [72].

$${\beta }_{\mathrm{i},\mathrm{j}}=\mathit{\Psi }·{\beta }_{\mathrm{coll}}$$

(26)

$${\beta }_{\mathrm{coll}}={\beta }_{\mathrm{coll},\mathrm{LS}}+{\beta }_{\mathrm{coll},\mathrm{DS}}=\frac{1}{6}·\dot{\gamma }·{\left(d_1+d_2\right)}^3+\frac{\pi }{4}·{\left(d_1+d_2\right)}^2·\left(v_{\sin \mathrm{k},1}+v_{\sin \mathrm{k},2}\right)$$

(27)

The collision frequency ${\beta }_{\mathrm{coll},\mathrm{LS}}$ is calculated depending on the mean shear rate $\dot{\gamma }$ and the diameters of the two colliding spheres $d_1$ and $d_2$. The collision frequency through differential settling ${\beta }_{coll,DS}$ depends on the particle sizes $d$ and their settling velocities $v_{\sin \mathrm{k}}$. The agglomeration efficiency $\mathit{\Psi }$ is approximated by the approach of Mumtaz and Hounslow [68,69], who assume that two particles collide in a flow field at an arbitrary angle to each other. The particles adhere due to particle–particle interactions, and in a supersaturated state, a crystalline bridge can grow. The shear forces acting on the particles define a characteristic time during which growth can take place. A stable agglomerate forms if the bridge is strong enough to withstand the shear forces. If not, the particle disaggregates. To describe this, a characteristic Mumtaz number $M$ relating the strength of the bridge and the hydrodynamic disruptive force induced by shear on an aggregate is formulated (Equation (28)). $d_{\mathrm{eq}}$ is the equivalent diameter of the aggregate calculated as the geometric mean $d_{\mathrm{eq}}=\sqrt{d_1·d_2}$. [70,71]

$$M=\frac{\mathrm{strength} \mathrm{of} \mathrm{the} \mathrm{bridge}}{\mathrm{hydrodynamic} \mathrm{disruptive} \mathrm{force}}=\frac{G·L^{*}·{\sigma }_{\mathsf{\gamma }}}{d_{\mathrm{eq}}^2·\eta ·{\dot{\gamma }}^2}$$

(28)

The Mumtaz number depends on the operating conditions (shear rate $\dot{\gamma }$ and growth rate $G$), the particle sizes ($d_1$,$d_2$), and a substance constant $L^{*}{\sigma }_{\mathsf{\gamma }}$. This substance constant is unknown only and must be determined by fitting it to experimental PSD data. Further, a failure criterion must be defined to determine whether a single collision leads to agglomeration. This criterion decides, depending on the yield strength, whether the crystalline bridge is strong enough or the particles disaggregate so that the agglomeration efficiency $\mathit{\Psi }$ is 1 or 0. For more detailed information, the reader is referred to the literature [68,69,70,73]. This approach includes the assumption of homogeneous particle suspension inside the whole slug, leading to the underestimation of collision frequency if this is not the case. To take this into account, the effective suspension volume (Figure 8) is determined and used in the calculation for agglomeration. Therefore, the particle number density $n$, which refers to the total slug volume $V_{slug}$, is replaced by the effective particle number density $n_{\mathrm{eff}}$, which refers to the effective suspension volume $V_{susp}$ according to Equation (29).

$$n_{\mathrm{eff}}=n·\frac{V_{\mathrm{slug}}}{V_{\mathrm{susp}}}$$

(29)

The effective suspension volume ($V_{\mathrm{susp}}=L_{\mathrm{susp}}·A_{\mathrm{susp}}$), Equations (S14)–(S16), is calculated using the centroid of suspension, which is a dimensionless quantity describing the height and width of the suspension. This centroid of suspension can be measured following an approach by Termühlen et al. [7] by analyzing a video of the slugs recorded directly after slug formation. The parity plots for the PSD results of the modeled experiments using the mechanistic kernel described above and the effective number density for agglomeration calculation are shown as orange triangles in Figure 9 by comparing the calculated and experimental measured characteristic sizes $d_{50}$ and $d_{90}$.

Five experiments are calculated well, and the $d_{50}$-values lie near the parity line and the $d_{90}$ values inside the ±20% confidential interval. However, three experiments cannot be calculated well, whereby both characteristics are underestimated. This trend can be observed as well when investigating the solution space of the mechanistic kernel parameter (Figure S3). Comparing these three experiments, it stands out that the five well-described experiments share the same parent population of seed crystals, leading to different growth and agglomeration behavior (Table S9). Therefore, the combined fit resulting in one parameter set cannot be used to evaluate all experiments nor extrapolate to different operating conditions and seed crystal populations.

Even though the combined fit of the mechanistic kernel cannot describe all experiments well, a trend in the results of the mechanistic kernel can be observed. Figure 10 shows the results of the separately fitted mechanistic kernel parameter against the corresponding effective suspension volume. The kernel parameter $L^{*}{\sigma }_{\mathsf{\gamma }}$ can be approximated with a linear function in a log-log diagram. Therefore, the kernel parameter can be correlated to the effective suspension volume using a power-law approach (Equation (30)).

$$L^{*}{\sigma }_{\mathsf{\gamma }}=0.085·{\left(\frac{V_{\mathrm{susp}}}{V_{\mathrm{slug}}}\right)}^{1.731}$$

(30)

The correlation of the kernel parameter is now used to calculate all experiments combined. As a further improvement, the suspension correlation developed in our previous publication [49] is implemented to evaluate the effective suspension volume along the tubing since the particle properties also change along the crystallizer. The green diamonds in Figure 9 show the parity plots comparing calculated and measured characteristic sizes of the PSD $d_{50}$- and $d_{90}$-values. The calculation of the PSD using the correlated mechanistic kernel is improved compared to the combined fit of using one fitted substance property constant. Nearly all calculated $d_{50}$- and $d_{90}$-values lie inside the ±20% confidence interval. Using the correlation for the substance constant $L^{*}{\sigma }_{\mathsf{\gamma }}$ to obtain the mechanistic kernel and applying the particle suspension correlation to evaluate the effective suspension volume, all experiments can be evaluated independent of used seed crystals using one parameter set, and the transfer to various operating points is enabled. Furthermore, it is now possible with the model to display and evaluate target variables and states along the tubing. This is shown in Figure 11 as an example of the concentration and PSD curve along the tubing for an experiment with $Q_{tot}$ = 20 mL min⁻¹, ${\epsilon }_L$ = 0.5, $L_{tubing}$ = 26.5 m, and κ = 1.8 K min⁻¹ (experimental data according to experiment EPS in Table S8).

4. Model Robustness

A sensitivity analysis (SA) is utilized to identify the critical model parameters that affect the output of the SFC model to check its robustness. The model’s output contains uncertainty due to both the model formulation and the parameters used. The uncertainty can be categorized into three sources: Stochastic, subjective, and structural [74,75]. The SA focuses on the subjective uncertainty arising from incomplete knowledge of the model parameters. All the parameters and correlations used in the PBE model are derived from experimental data, which have a degree of uncertainty and propagate to the model output. The SA aims to determine the impact of perturbations in individual parameters on the model output and rank them based on their influence. To improve the predictive model’s performance, the uncertainty of the parameters with the strongest influence must be reduced. Further, these parameters should be determined precisely when transferring the model to different substance systems.

A global SA is used, which is based on the Morris method, determining the elementary effects of individual parameters by varying one parameter at a time and covering the entire parameter range [75,76]. In implementing the Morris method on the SFC model, all above-derived submodels, namely the slug flow hydrodynamics, energy and mass balances, as well as the crystallization phenomena crystal growth and agglomeration, are included. The influence of the model parameters on the relative yield $Y_{rel}$ (varying outlet concentration $c_{out}$ for a given start concentration) and characteristic values of PSD $d_{50}$ and $d_{90}$ are investigated.

Table 4. Model parameters of submodels, their mean values, and uncertainty ranges varied in the SA.

Submodel	Parameter	Mean Value	Unit	Uncertainty Range
Solubility	$k_{\mathrm{c},1}$	0.1124	${\mathrm{g}}_{\mathrm{Ala}}·{\mathrm{g}}_{\mathrm{sol}}^{-1}$	±5%
	$k_{\mathrm{c},2}$	9.085 × 10−3	$°{\mathrm{C}}^{-1}$	±5%
Slug flow hydrodynamics	$L_{\mathrm{slug}}$	calc	m	±20%
	$\Delta p$	calc	bar	±20%
Energy balance	${\lambda }_{\mathrm{FEP}}$	0.209	$\mathrm{W}·{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$	±20%
	${\alpha }_{\mathrm{ENV}}$	calc	$\mathrm{W}·{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$	±20%
Crystal growth	$A_{\mathrm{BCF}}$	1.160 × 10−4	$\mathrm{m}·{\mathrm{s}}^{-1}$	±20%
	$B_{\mathrm{BCF}}$	0.618	/	±20%
Agglomeration	$k_{\mathrm{Agg},1}$	0.085	$\mathrm{N}·{\mathrm{m}}^{-1}$	±20%
	$k_{\mathrm{Agg},2}$	1.731	/	±20%

lists the model parameters and their respective assumed uncertainty ranges.

The results of the SA are shown in Figure 12 displayed as elementary effects for the model parameters of the submodels for the respective target parameters (yield and characteristics of PSD). The elementary effects are only a qualitative measure and can only be used to rank the parameters according to their contribution. The absolute values of the effects have no relevance. Elementary effects lower than 0.5 are considered insignificant.

Figure 12 shows that the solubility curve’s two parameters $k_{\mathrm{c},1}$ and $k_{\mathrm{c},2}$, the pressure drop $\Delta p$ determining the RT, the convective heat transfer to the environment ${\alpha }_{\mathrm{ENV}}$, and the second parameter of the agglomeration correlation $k_{\mathrm{Agg},2}$ have a high impact on all target parameters. The first solubility factor has the most significant influence despite the lower uncertainty range of ±5% since it governs the driving force in crystallization. Furthermore, the SA highlights the importance of accurately modeling slug flow hydrodynamics, represented by the pressure drop $\Delta p$ and temperature profiles. Additionally, the low influence of agglomeration on the relative yield can be attributed to the reduced particle surface area and the need for more growth units to be incorporated to facilitate particle growth. In relation to the product PSD, the influence of agglomeration is high because it affects the particle size directly.

5. Guideline for Model Transfer to New Substance Systems

Based on experience, the model derivation, and the results of the SA, it is now possible to define the necessary steps for transferring the model to a new substance system (Figure 13).

The physical properties of the solution can be found in the literature or have to be determined experimentally and should be regressed to cover the operating range. Depending on these properties, the slug flow hydrodynamics and the energy balance can be predicted. The adjusted slug length correlation includes the properties of the solution and, therefore, is transferable to different substance systems. Additionally, the SA showed that uncertainties in the slug length have only a low impact on the product properties. Therefore, a certain degree of uncertainty in the slug length is tolerable. The energy balance calculation was adjusted to the SFC setup, and the temperature profile is mainly determined by the TM. Using a similar apparatus, the energy balance can be transferred to a different substance system without further adjustments.

The SA showed the importance of the solubility curve. Therefore, the curve should be regressed and validated based on precise experimental data. The growth rate constants can be determined from batch crystallization data at similar operating conditions. Derivation from batch data is recommended because probes can be sampled at will during the crystallization process, in contrast to the SFC. Further, the experimental effort required is substantially lower. The agglomeration behavior, however, must be studied by key experiments inside the continuous SFC because it cannot be transferred from different substance systems or apparatuses. Following the model-based approach, twofold crystallization experiments at two operating points covering the desired operating range must be conducted to identify the suspension volume and validate the suspension correlation [49], as well as to describe the agglomeration tendency for the respective substance system at the given suspension volumes. This is necessary to determine a correlation comparable to Equation (30). Further, this correlation can be validated with a further experiment. With these few experiments, the agglomeration behavior for different operating conditions can be mapped, significantly reducing the experimental effort.

6. Conclusions

In order to produce high-quality pharmaceutical products, it is crucial to be able to adjust specific target parameters, such as particle size distribution (PSD) or yield. To achieve this, the goal was to develop a transferable model that could predict the crystallization behavior inside the Slug Flow Crystallizer (SFC). This model considers the slug flow hydrodynamics, energy and mass balances, and the crystallization phenomena of growth and agglomeration. Its purpose is to improve understanding of the process, estimate the effects of operating parameters on target parameters, and predict crystallization behavior for different substance systems with minimal experimental effort.

The pressure drop inside the crystallizer mainly affects the reduction in residence time (RT) by gas expansion. A correlation for slug length was used, and with this value, the pressure drop and actual RT were determined across the entire operating range of the SFC. Furthermore, energy balances were conducted using heat transfer rates utilizing Nusselt correlations for convective heat transfer coefficients and were validated using experimental data. The Burton, Cabrera, and Frank model was used for crystal growth, while a mechanistic kernel describing the physical phenomena of agglomeration was connected with a correlation to predict the effective suspension volume, allowing for the reliable estimation of collision frequency and agglomeration rate. This combination enabled mechanistic modeling of crystallization behavior for different operating parameters independent of the seed crystals used. Furthermore, a sensitivity analysis was conducted to identify decisive model parameters and evaluate the influence of uncertainties on the relative yield and on the characteristic values of product PSD. In the end, a guideline to transfer the model to a new substance system is formulated, applying only a few key experiments.

The developed model helps to enable robust and long-term stable operation of the SFC for different substance systems so that the SFC can also be used in the process development phase of pharmaceutical substances, where only a few materials are available. For production purposes, it is necessary to set up model predictive control and ensure the traceability of the product, which can be achieved by this mechanistic model.