Describing the Drying and Solidification Behavior of Single Suspension Droplets Using a Novel Unresolved CFD-DEM Simulation Approach

Abstract

With increasing computational capacities and advances in numerical methods, a simulation-based process and product design approach for spray dried products has become the focus of many researchers. In this context, a novel unresolved CFD-DEM simulation approach for single suspension droplets is presented, which models the drying and solidification process in a hot gas environment. The solidification process is described by the formation of bonds between the primary particles in the suspension droplet, showing a plausible solidification behavior, which starts when a critical solid concentration is reached at the droplet surface. The drying conditions for the single-droplet simulations are determined from a large-scale spray dryer simulation for different droplet size classes. The resulting solid particles from the CFD-DEM simulation show a higher tendency to form hollow particles at high drying rates when the locking point is reached at earlier stages of the drying process. Using these results, suitable operating conditions for the spray dryer can be determined, which favor the formation of particles with a targeted morphology.

1. Introduction

Spray drying processes are widely used in industry to produce large amounts of granular material from liquid solutions or suspensions. Due to the increasing demand for materials with different granule properties, size distributions and particle morphologies, the use of computational methods for the design and optimization of spray drying processes has gained a significant scientific interest [1].

The spray drying process can be divided into three main stages: (1) the atomization step, in which a spray of small droplets is created from a liquid feed; (2) the drying and solidification process of the single droplets; and (3) the separation of the granular material from the drying gas. From a phenomenological perspective, the drying process is often divided into two distinct stages. During the first drying stage, the droplet shrinks due to the unhindered evaporation of the solvent material until the droplet solidifies and forms a wet particle [2]. This drying stage is also called constant rate period as the droplet surface decreases at a constant rate, which is also expressed by the so-called d²-law [3]. The second drying stage, also called the falling rate period, consists of the drying of the residual solvent in the wet particle, which is limited by the transport processes of the solvent vapor within the particle itself [2].

For over a hundred years, spray dryers have been operated with a focus on dairy and detergent applications [1]. More recently, the focus has also shifted towards the production of more specialized granular materials like catalysts that require a specific particle morphology to achieve optimal performance. The development and design of spray drying processes for such materials poses a great challenge [4]. On the one hand, the mechanisms of particle formation occurring at droplet scale are difficult to model due to their inherent complexity, and on the other hand, during the scale-up of laboratory and pilot scale experiments to production scale, changing flow-patterns and changing atomization behavior may severely impact the resulting product properties [5].

During the last decades, several methods for the numerical modeling of granular flows have been developed. Since most industrial processes involving granular materials include a secondary fluid phase (typically air or water), the effect of solid–fluid interactions must be accurately modeled. A promising approach for simulating granular-fluid systems is the coupled CFD-DEM method, with which the motion of the particles is resolved by means of the Discrete Element Method (DEM), while the Computational Fluid Dynamics (CFD) method is used to calculate the interstitial fluid flow [6]. Reviews on DEM and CFD-DEM methodologies have been published by Zhu et al. [7,8] and recently by El-Emam et al. [9].

The coupled CFD-DEM method became a crucial tool for understanding the physics of particulate systems and for optimizing industrial process design, as described in the works of Golshan et al. [10]. In this work, a novel methodology for predicting the average granule size and morphology in spray drying applications for suspensions is presented, as shown in Figure 1. On the one hand, this approach consists of an unresolved CFD-DEM simulation of a single evaporating spray droplet containing a suspension of solid primary particles, and on the other hand, of a CFD simulation of the spray dryer. The latter is used to evaluate the average drying conditions that the single droplets may experience within the spray dryer.

The CFD-DEM simulation used for droplet-scale simulation is of an unresolved type, meaning that the flow around each individual particle is not resolved since the particle sizes are much smaller than the computational grid. This unresolved approach drastically reduces the computational time in comparison with resolved CFD-DEM, albeit at the cost of losing access to small-scale flow structures.

The information obtained from the macroscopic spray-drying simulation is used to set the time-dependent boundary conditions for the microscopic droplet-scale simulation. The microscopic simulation resolves the mass, momentum and energy exchange between a single spray droplet and the surrounding gas, modeled as two distinct fluid phases using a Volume of Fluid (VOF) method. Furthermore, the momentum and energy exchanges between the single droplet and the suspended solid primary particles, which are modeled as DEM particles, are resolved. The interface between the liquid droplet and the gas deforms due to the balance between drag and surface tension forces, while it recedes due to the solvent evaporation. The primary solid particles recirculate inside the droplet, are repulsed by the liquid–gas interface and agglomerate to form a crust or a full solid granule. Due to the high number of primary particles inside the droplets, a coarse-graining approach is used. A detailed description of the model employed for the microscopic simulation is presented in Section 2.

The internal circulation of the primary particle depends on many factors such as the droplet size, the surface tension and the viscosity of the liquid–solid mixture, as well as the evaporation rate. The constant rearrangement of the solid suspension might affect the formation of coherent solid structures inside the droplet such as shells at the surface. However, this behavior depends also on the material’s properties, such as the concentration and nature of the binder present in the slurry. These complex and interwoven dynamics, which determine the size and the morphology of the dried granules in spray drying processes, have been investigated in this work using unresolved CFD-DEM simulations at a droplet scale.

In the literature, there are different approaches that deal with the detailed modeling of the evaporation process of single droplets. Breinlinger et al. [11,12] used coupled CFD-DEM simulations for this purpose and derived a capillary force model that is based on the Young–Laplace and Young equations. The model neglects the local curvature of the solid–liquid contact line. For their simulations, a coarse-graining of the particles from 1 µm to 20 µm was used to reduce the computational cost. From the simulation results, the authors showed the influence of the cohesive and capillary forces on the resulting granule morphologies. As a first indicator for the validity of the model, the simulation results were compared with results for dried particles from the literature.

Nishiura et al. [13] combined the VOF and DEM methods, along with the Constrained Interpolation Profile (CIP) method, introduced by Yabe et al. [14]. The CIP method, a semi-Lagrangian approach, calculates the motion of individual particles within a gas–liquid multiphase flow. The authors employed the liquid immersion technique to observe the internal structure of granules and investigated the impact of drying conditions on granule morphologies, considering the drying rate distribution around the droplet under varying air flow conditions. They found that increased friction and vertical capillary forces in the top part of the droplet led to the formation of particles with depressions. In another study by Miyazaki et al. [15], the coupled DEM-CIP simulations were used to examine the binder segregation during the drying process. They observed that adding binder material resulted in hindered particle movement due to increased viscosity and reduced surface tension, leading to the formation of more hollow particles. Washino et al. [16] employed the CIP method in their study to solve the liquid–gas two-phase flow and used the VOF method to account for the multiphase nature of the system. They introduced a novel approach, known as the continuous capillary force (CCF) model, to consider the stresses between suspended particles induced by capillary forces. Instead of explicitly tracking the contact line, they considered a transition region of finite thickness. The interaction force between the fluid and solid was modeled using the immersed boundary method. The authors demonstrated that the CCF model allows the evaluation of the capillary force at low computational costs, even in the case of an inexact determination of the location of the contact line.

The novel simulation framework presented in this study makes use of the possibilities offered by CFD and DEM to reproduce physical systems at different length and time scales. At the macroscale, the trajectories and temperature histories of droplets in the whole apparatus are determined using pure CFD simulations. The data are then used as input for the boundary conditions on the microscale in highly resolved coupled CFD-DEM simulations of individual droplet drying. This approach allows linking of the operating conditions of a spray drying process with the resulting particle properties, like size and morphology, which is crucial for operating and optimizing spray drying processes.

2. Materials and Methods

2.1. Model Structure

In this section, the models used for the microscopic simulations are described in more detail. The microscopic single-droplet simulation, on which this work is focusing, is a 4-way coupled CFD-DEM simulation of a single evaporating droplet containing a suspension of solid particles. Due to the presence of three phases, namely gas, liquid and solid, a VOF solver (cfdemSolverInterDyMPC from CFDEM^®couplling 5.3.0) has been coupled to the Aspherix^® DEM solver (version 5.3.0). CFDEM^®couplling is an extension of the open-source software package OpenFOAM^® [17] designed for coupled CFD-DEM simulations. In this work, the OpenFOAM^® version 8 is used.

cfdemSolverInterDyMPC is a coupled CFD-DEM solver based on InterDyMFoam, which is a VOF solver from OpenFOAM^® for numerically solving the Navier–Stokes equations in cases of two incompressible, isothermal, and immiscible fluids. In contrast to InterDyMFoam, cfdemSolverInterDyMPC also considers a discrete solid phase with the continuity, momentum, and energy equations, as expressed in Equations (1)–(4). Detailed information on the numerical procedure of a coupled CFD-DEM simulation can be found in [18].

A cubic simulation domain, in which the center corresponds to the center of mass of the droplet at every time point, has been selected, as shown in Figure 2a. The hot gas enters the computational domain from the bottom and may leave the system at the top or sides of the cubic domain. The inlet velocity and temperature in the microscopic simulation correspond to the average relative velocity between the droplets and the gas, and to the average gas temperature at the droplets, respectively. The latter information is obtained by performing and post-processing a macroscopic simulation of the entire spray dryer [19,20].

The gas and liquid phase are both incompressible and have constant thermophysical properties. Furthermore, both phases consist of a single species, namely nitrogen for the gas and pure water for the liquid. Here it should be remarked that the vapor concentration has not been explicitly modeled. Nonetheless, the overall fluid mass is conserved, since the liquid lost by the droplet through evaporation is converted into an equal mass of the surrounding gas.

2.2. The CFD Model

The governing equations for the two-phase system are the following [21,22]:

(1)$\frac{\partial \left({\epsilon }_f{\alpha }_{\mathcal{l}}\right)}{\partial t}+\nabla ·\left({\epsilon }_f{\alpha }_{\mathcal{l}}\mathit{u}\right)=-\frac{\dot{m}}{{\rho }_{\mathcal{l}}}$

(2)$\frac{\partial \left({\epsilon }_f{\alpha }_g\right)}{\partial t}+\nabla ·\left({\epsilon }_f{\alpha }_g\mathit{u}\right)=\frac{\dot{m}}{{\rho }_g}$

(3)$\frac{\partial \left({\epsilon }_f\mathrm{ }\rho \mathrm{ }\mathit{u}\right)}{\partial t}+\nabla ·\left({\epsilon }_f\mathrm{ }\rho \mathrm{ }\mathit{u}\otimes \mathit{u}\right)=-{\epsilon }_f\nabla p\mathrm{ }+\nabla ·\left({\epsilon }_f\mathrm{ }\mathit{\nu }\mathrm{ }\nabla \mathit{u}\right)+{\epsilon }_{f\mathrm{ }}\rho \mathrm{ }\mathit{g}+{\epsilon }_f{\mathit{f}}_{\sigma }-{\mathit{R}}_{pf}$

(4)$\frac{\partial \left(c_p\rho T\right)}{\partial t}+\nabla ·\left(c_p\rho T\mathit{u}\right)=\nabla ·\left(k\nabla T\right)-\dot{m}\Delta h_v$

Equations (1) and (2) are the continuity equations for the two fluid phases. ${\epsilon }_f$ is the fluid volume fraction, and ${\alpha }_{\mathcal{l}}$ and ${\alpha }_{\mathrm{g}}$ are the liquid and gas volume fractions, respectively. u is the velocity vector, $\dot{m}$ is the volume specific mass transfer between the two phases, and ${\rho }_{\mathcal{l}}$ and ${\rho }_{\mathrm{g}}$ are the density of the liquid and gas, respectively. Equations (3) and (4) are the Navier–Stokes and the energy equations. Here, the fluid density, kinematic viscosity, specific heat capacity and thermal conductivity of the fluid mixture are defined as follows:

$\rho ={\alpha }_g{\rho }_g+{\alpha }_{\mathcal{l}}{\rho }_{\mathcal{l}},\nu ={\alpha }_g{\nu }_g+{\alpha }_{\mathcal{l}}{\nu }_{\mathcal{l}},c_p={\alpha }_g{c}_{p,g}+{\alpha }_{\mathcal{l}}{c}_{p,\mathcal{l}},k={\alpha }_g{k}_g+{\alpha }_{\mathcal{l}}{k}_{\mathcal{l}}.$

The other terms in the Navier–Stokes equation are the pressure p, the gravitational acceleration vector g, the surface tension force ${\mathit{f}}_{\sigma }$ and the tensor representing the momentum exchange between the fluid and the particulate phase ${\mathit{R}}_{pf}$. Finally, the temperature T and the latent heat of vaporization $\mathsf{\Delta }{h}_v$ can be found in Equation (4).

The mass transfer $\dot{\mathrm{m}}$ between the two fluid phases is due to the evaporation of the liquid phase of the droplet. This term has been modeled by assuming a linear relation between the fluid temperature and the evaporation rate, as shown in the formula below:

(5)$\dot{m}=C_{ev}\mathrm{ }\max \left(T-T_{ref},\mathrm{ }0\right)\mathrm{ }\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{ }\mathrm{ }\mathrm{ }{\alpha }_{\mathcal{l}}\in [0,0.5[,$

(6)$\dot{m}=0\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{ }\mathrm{ }{\alpha }_{\mathcal{l}}\in \mathrm{ }[0.5,\mathrm{ }1],$

where $T_{ref}$ is a reference temperature, below which the evaporation does not occur, and $C_{ev}$ is the evaporation constant. As can be seen from Equation (6), the evaporation in the droplet bulk is prevented by setting the mass transfer term equal to zero for ${\alpha }_{\mathcal{l}}\ge 0.5$, where ${\alpha }_{\mathcal{l}}=0.5$ denotes the interface between the gas and liquid, as shown in Figure 2b. A verification of the CFD model in terms of mass and energy conservation is presented in Appendix A.

The shear force exerted by the gas on the droplet surface results in an internal recirculation of the primary particles. The momentum exchange between the droplets and fluid is calculated by applying the drag correlation of Koch and Hill [23]. At the droplet surface, the primary particles experience a capillary force which pushes them towards the center of the droplet. This interaction has been modeled based on the approach from Breinlinger et al. [11]. The internal circulation and the receding movement of the droplet surface determine the collisions between the primary particles, which are resolved by the DEM part of the model.

2.3. The DEM Model

The DEM was introduced by Cundall and Strack [24]. A description of the fundamentals is provided in this section. Further details on the contact physics and implementation issues are available in the literature [6,25,26].

Since the DEM is a Lagrangian method, the equations of motion are solved for all particles in the computational domain to track their trajectories. The force and momentum balances for a generic particle with index $i$ can be derived from Newton’s second law and read as follows:

(7)$m_i{\ddot{\mathit{x}}}_i={\mathit{F}}_{i,n}+{\mathit{F}}_{i,t}+{\mathit{F}}_{i,f}+{\mathit{F}}_{i,b},$

(8)$I_i·{\dot{\mathit{\omega }}}_i={\mathit{T}}_{i,t}+{\mathit{T}}_{i,r},$

where ${\mathit{F}}_{i,n}$ is the normal particle–particle contact force vector and ${\mathit{F}}_{i,t}$ is the tangential particle–particle contact force vector. ${\mathit{F}}_{i,f}$ is the force caused by the interphase interaction of the particles with the surrounding fluid, including the drag, pressure gradient, buoyancy, and capillary effects. Other body forces like gravity, electrostatic or magnetic forces are summarized by ${\mathit{F}}_{i,b}$. In Equation (8) for the rotation of a particle, ${\mathit{T}}_{i,t}$ is the torque caused by the tangential force acting on the particle in the point of contact, while ${\mathit{T}}_{i,r}$ is an additional torque on the particle that can be used to model non-sphericity by means of a so-called rolling friction model [27]. The moment of inertia, and linear and angular accelerations are denoted as $I_i$, ${\ddot{\mathit{x}}}_i$ and ${\dot{\mathit{\omega }}}_i$, respectively.

In the soft-particle DEM approach used in this study, the particles can overlap with each other or with the walls. The simplest contact model is a linear spring-dashpot system, which models the normal and tangential contact forces as follows [24]:

(9)${\mathit{F}}_n={-k}_n{\mathit{\delta }}_n+c_n{\dot{\mathit{\delta }}}_n,$

(10)${\mathit{F}}_{\mathit{t}}=\min \left[\left|c_t{\dot{\mathit{\delta }}}_t+k_t{\int }_{t_{c,0}}^{t_c}{\dot{\mathit{\delta }}}_{t}dt\right|,\mu {\mathit{F}}_n\right].$

Here, $k_n$, $k_t$ and $c_n$, $c_t$ are the normal and the tangential spring and damping coefficients, respectively. ${\mathit{\delta }}_n$ and ${\dot{\mathit{\delta }}}_n$ are the normal overlap and the relative normal velocity of the particles, and $\dot{{\mathit{\delta }}_t}$ is the relative tangential velocity of the particles. The tangential overlap ${\mathit{\delta }}_t$ is truncated to fulfill the Coulomb friction criterion, where $\mu$ is the coefficient of friction. The integral term represents an incremental spring that stores energy from the relative tangential motion, representing the elastic tangential deformation of the particle surface during contact. The dashpot part accounts for the energy dissipation of the tangential contact.

To describe the drying and solidification process of a single suspension droplet, additional force models are included to account for the capillary force acting on the primary particles in the vicinity of the droplet surface as well as the solidification mechanisms occurring via the formation of solid bonds between two particles. While the capillary force on particles at the droplet surface is considered as a body force, the solid bond model has the form of a contact model. Both models are introduced in the following paragraphs.

2.3.1. The Capillary Force Model

The capillary force ${\mathit{F}}_{\mathit{c}}$ acts on the primary particles due to the receding liquid surface. According to Breinlinger et al. [11], a simple capillary force model can be obtained by integrating the surface tension along the line of contact. This leads to the following expression for the capillary force:

(11)${\mathit{F}}_{\mathit{c}}=2\pi r\sigma \sin \left(\omega \right)\sin \left(\omega +\mathsf{\Theta }\right)\mathit{n},$

where $r$ is the radius of the sphere, $\sigma$ and $\mathsf{\Theta }$ are the surface tension and contact angle of the liquid, $\omega$ is the wetting angle, and $\mathit{n}$ is the unit vector normal to the droplet–gas interface and directed from the droplet to the gas, as shown in Figure 3. The wetting angle can be expressed as $\omega =-\arccos \left(-\frac{h_c}{r}\right)$, where $h_c$ is the distance of the DEM particle center from the droplet–gas interface. If the particle center is outside the droplet, then $h_c>0$, otherwise $h_c\le 0$.

The original model from Breinlinger et al. [11] has been modified by the authors as follows:

(12)${\mathit{F}}_{\mathit{c}}=\frac{2}{\beta }\pi r\sigma \sin \left(\omega \right)\sin \left(\omega +\mathsf{\Theta }\right)\mathit{n} \mathrm{f}\mathrm{o}\mathrm{r}-\beta r < h_c<\beta r,$

where $\beta$ is the interface scale parameter. If $\beta =1$, then Equation (12) corresponds to Equation (11), which means that the modified Breinlinger model corresponds to the original one. For $\beta >1$, the capillary force starts acting before the particle intersects the interface and its magnitude is scaled by the factor $\beta$.

In the absence of other external forces, the change in energy that a particle experiences while crossing the interface can be calculated as the integral of the capillary force along the particle trajectory. It can be mathematically demonstrated that the work performed by the capillary force on the particle does not depend on interface scale parameter; see Equation (13):

(13)$E_c={\int }_{-r}^r\left|{\mathit{F}}_{\mathit{c}}\right|dh = {\int }_{-\beta r}^{\beta r}\frac{\left|{\mathit{F}}_{\mathit{c}}\right|}{\beta }dh=\frac{8}{3}\pi \sigma r^2\cos \left(\mathsf{\Theta }\right).$

Hereafter, the formulation in Equation (12) is considered as it allows stable CFD-DEM simulations with larger CFD-DEM coupling times due to lower oscillations of the capillary force. The implementation of the capillary force model in CFDEM^®coupling has been verified by measuring the change in kinetic energy of a solid particle crossing the liquid–gas interface for different contact angles, as described in Appendix A.

2.3.2. The Particle Bond Model

The DEM particles inside the droplet can form bonds with each other according to the model proposed by Potyondy et al. [28]. The bond formation and breakage are schematized in Figure 4. If two initially unbound DEM particles come in contact and their distance drops below a certain threshold, then a bond is formed. By introducing the relative distance at which bonds form, hereafter indicated as $d_{b,min}$, the condition for bond creation reads as follows:

(14)$\frac{d}{r_1+r_2}<d_{b,min},$

where $d$ is the distance between the particle centers, and $r_1$ and $r_2$ are the particle radii.

The bond, modeled as a cylindrical beam, reacts to normal, tangential, and torsional deformations with elastic and damping forces and torques, which depend not only on the size of the bond but also on material properties like, for example, the Young’s modulus and the Poisson ratio of the bond. The bond can break either because the stress, either normal or tangential, surpasses a critical value or because the distance between the particles is larger than a critical bond distance. The microscopic simulations performed for this work have been performed considering the critical relative distance as a breakage criterion, which means that the bond between two particles breaks if the following condition is met:

(15)$\frac{d}{r_1+r_2}>d_{b,max},$

where $d_{b,max}$ is the maximum relative distance between bonded particles. In the following section, the simulations for different values of $d_{b,min}$ and $d_{b,max}$ are presented, showing the important impact of this parameter on the resulting granule morphology.

3. Results

3.1. Model Analyis

In this section, an analysis of the single-droplet model behavior is first presented to show the plausibility of the modeling approach. In a sensitivity analysis regarding the bond parameters, the impact on the resulting particle morphology is investigated and finally, the influence of different drying conditions on the particle formation processes is shown.

Simulation Setup

The simulation setup for the model analysis is shown in Figure 5. The remaining model parameters are listed in Table A1. The first property to be investigated is the evaporation rate of the single suspension droplet. The first drying stage of a droplet can be described by the d²-law, which states that the droplets surface decreases linearly in time as follows:

(16)$d^2=d_o^2-\beta t,$

with the initial droplet diameter $d_0$ and the evaporation coefficient $\beta$ [29].

The d²-laws for two different gas inlet temperatures are shown in Figure 6. Here, the evaporation constant $C_{ev}$ was adjusted to match the experimental findings of two single-droplet experiments in an acoustic levitator presented by Buchholz et al. [30]. The initial droplet sizes for the experiments are 1.15 mm for the experiment at 60 °C and 1.54 for the experiment at 140 °C. The resulting evaporation coefficients have values of 8300 and 31,800 µm²/s, respectively.

As illustrated in Figure 6, the model follows the d²-law reasonably well. Thus, the chosen approach to model the evaporation process within the CFD-DEM framework is suitable for describing the evaporation process in a plausible way.

A moving droplet in a surrounding gas environment experiences a drag force which leads to an internal circulation of the fluid bulk [31]. An analytical solution for a moving liquid droplet in a gas, neglecting a possible evaporation process, is given by the similarity with the Hadamard–Rybczynski solution [32]. The resulting velocity magnitude of the single-droplet simulation is compared with this theoretical solution at the center point of the droplet, shown as a black dot in Figure 7. The velocity magnitude of the center point in the simulation equals 0.0957 mm/s, which is in good agreement with the value of 0.0729 mm/s from the Hadamard–Rybczynski solution. The deviations may be explained by the numerical error from the simulation mesh and the additional deformation of the droplet, which is neglected in the analytical Hadamard–Rybczynski equation. In the unresolved CFD-DEM model, the minimum size of the numerical grid is limited by both the simulation time and the size of the primary particles in the suspension droplets. The cells should be roughly an order of magnitude larger than the largest primary droplets for the hypothesis of an unresolved CFD-DEM simulation to hold.

As mentioned above, the simulated droplet shows a deformation due to the drag and pressure forces that act from the surrounding gas flow on the droplet surface. The surface tension acts in return against this deformation such that a constant deformation is shown from the simulation results, when the evaporation is negligibly small or switched off. An empirical finding by Hsiang & Faeth [33] shows a correlation of the droplet deformation with respect to the Weber number $We=\frac{\rho u^{2}d}{\sigma }$ as follows:

(17)$\mathsf{\Psi }=1+0.19\sqrt{We}.$

Here, $\mathsf{\Psi }$ is the droplet aspect ratio, i.e., the ratio of the length of the major axis to the minor axis of the ellipsoid body of the deformed droplet. For the CFD-DEM simulations, the Weber number has been adjusted by adjusting the velocity of the entering gas accordingly. The comparison of the simulation results with this empirical correlation of the aspect ratio is given in Figure 8, showing a good agreement of the simulations with the empirical findings of Hsiang & Faeth [33].

As described in Section 2.3, the solidification of the suspension droplet is captured in the model by allowing the formation of bonds between the primary particles in the droplet when the particles come into contact. These bonds may break if the relative distance between the bonded particles becomes too large due to counteracting forces on the particles. In two exemplary simulations, the solidification behavior of the droplets is investigated over time for different settings of the bond parameters $d_{b,min}$ and $d_{b,max}$, as shown in Figure 9. For case (a), the bond parameters are $d_{b,min}=1.8$ and $d_{b,max}=2.4$, and the inlet gas temperature is 80 °C. For case (b), $d_{b,min}=1.8$ and $d_{b,max}=3.5$, with an inlet gas temperature of 170 °C. A bigger value of $d_{b,min}$ leads to an earlier bond formation at a smaller overlap of the two primary particles. A bigger value for $d_{b,max}$ results in stronger bonds, meaning the primary particles need a larger relative distance from each other before the bond breaks.

As can be seen in Figure 9, weaker bonds (i.e., low $d_{b,max}$) lead to the formation of an overall spherical and dense agglomerate, while stronger bonds (i.e., high $d_{b,max}$) lead to an earlier start of the solidification of the suspension droplet (1.7 s vs. 4 s). The formed particle also shows a hollow structure as a strong crust forms at the droplet surface. Regarding the overall solidification behavior, the bonds tend to form in the vicinity of the droplet surface when a critical solid concentration is locally reached, leading to more overlap between the DEM particles. This behavior mimics the actual particle formation behavior very well when binder material is present. As the binder does not evaporate during the drying process, it accumulates in the surface region when the droplet shrinks, leading to stronger binding forces between the primary particles and, thus, also to a tendency for early crust formation at the droplet surface.

The solidification behavior of the model can also be quantified by the evaluation of the Mean Coordination Number (MCN), which is defined as the mean number of bonds per primary particle in the suspension droplet. The temporal change in the MCN is shown for both simulation cases in Figure 10.

In the first case (blue line), the MCN shows a rather steady increase after 1.5 s and a sharp increase after 4.5 s, where the latter event indicates the reaching of the locking point. At the locking point, the droplet forms the dense spherical particle structure seen in the far right of Figure 9a. For the second case (green line), the MCN sharply increases already after 0.5 s. Here, the primary particles start to quickly form bonds at the droplet surface, which leads to the formation of a strong shell at the droplet surface and, thus, to the formation of a particle with a hollow core structure.

3.2. Particle Morphologies

The impact of the bond parameters on the resulting particle morphology is further investigated for different pairs of bond parameters. The morphology map of the resulting solidified suspension droplets is shown in Figure 11. The model parameters are listed in Table A1. Except for the bond parameter pairs, the remaining simulation parameters and the drying conditions are identical for all simulations in this parametric study.

As shown in Figure 11, particles with a hollow core are always formed at $d_{b,max}=4$, namely for the highest value of the bond strength. In the case of a medium bond strength, namely $d_{b,max}=3$, the hollow particles form only when the bonds have the tendency to form earlier, i.e., for $d_{b,min}=1.85$ and $d_{b,min}=1.9$. The shell obtained in the last two cases is thicker than the one for $d_{b,max}=4$. This can be phenomenologically explained by the lower bond strength between the primary particles ($d_{b,max}=3$), which leads to more breakage events until the crust reaches a certain thickness at which it can sustain the capillary forces imposed by the liquid–gas interface. Furthermore, the hollow particles with a thicker crust also take longer to form a stable crust and are, therefore, smaller than the hollow particles in the right column, which form a stable crust earlier. A dense spherical particle forms only in the top-left case, where the bonds between the primary particles only form when the primary particles are densely packed within the droplet. For the remaining cases, the bonds form earlier or are slightly stronger, such that irregular particle structures form, which have a higher porosity than a densely packed spherical particle.

3.3. Drying Rate

As shown before in the literature, the drying rate of the droplets has a strong impact on the resulting particle morphologies [34]. The drying rate influences the speed at which the droplet surface shrinks, which may lead to a fast accumulation of primary particles at the droplet surfaces, i.e., a critical solid concentration is reached locally, and therefore, to an earlier start of the solidification process at the droplet surface.

To investigate the impact of the drying rate on the particle morphology, the evaporation of zeolite suspension droplets in a Niro Minor pilot-scale spray dryer was simulated. This application scenario has been investigated experimentally and numerically by Buchholz et al. [20]. In their spray drying experiments, it could be seen that the resulting zeolite particles showed both dense and hollow structures. In numerical investigations, a pure CFD simulation was performed to evaluate the drying conditions for the differently sized droplets and particles within the spray dryer [19,20]. For this purpose, the authors set up an extended CFD model, which also considered the different drying kinetics of the droplets before and after the solidification. The evaluation of the mean drying conditions for the different droplet size classes in terms of the mean values and the standard deviations for the evaporation coefficient $\beta$ according to the respective $d^2$-laws are shown in Figure 12. The evaporation coefficients of the smaller droplets show larger standard deviations as the smaller droplets may move on very different trajectories through the spray dryer, leading to the strong deviations in the actual drying conditions of the single droplets [19].

From the diagram in Figure 12, three cases for the evaporation coefficients are defined for the single-droplet simulations to account for the mean and extreme points of the drying rates, namely $\beta =$ 3000, 5000 and 7000 μm²/s. The simulations are conducted for two cases of initial droplet size, $d_0=25\mathrm{\mu }\mathrm{m}$ and $d_0=40\mathrm{\mu }\mathrm{m}$. The main model parameters are listed in Table A1, and the resulting particle morphologies are shown in Figure 13.

For both droplet sizes, the resulting particles show denser structures at lower drying rates and hollow structures for higher drying rates, i.e., larger values of the evaporation coefficient $\beta$. This is explained by the faster-receding droplet surface at higher evaporation rates, which leads to an accumulation of DEM particles in this region. Due to the high solid concentration, these particles then have a higher tendency to overlap and form bonds.

For smaller droplets, the resulting particles show a denser structure for lower drying rates and a hollow structure for higher drying rates. The larger particles show more irregular structures, such as mushroom-cap morphology, at the lowest value of the drying rate. The mushroom-cap morphology may be explained by primary particles that may sediment or accumulate at the lower part of the droplets when the drying process takes longer. The overall particle formation behavior shown by the simulation, especially the tendency to form hollow structures at higher drying rates, shows the suitability of the applied modeling approach to describe the resulting particle morphologies of spray dried granules.

4. Conclusions

In this work, a modeling approach based on coupled CFD-DEM simulations is presented to capture the evaporation of a suspension droplet its solidification process. The model, implemented in CFDEM^®coupling and Aspherix^® (version 5.3.0), has been thoroughly verified (see Appendix A) and validated using experimental (see Figure 6), analytical (see Figure 7), and empirical results (see Figure 8). The time and size-dependent boundary conditions for the droplet-scale CFD-DEM simulations have been obtained from CFD simulations of lab-scale spray drying experiments [20,21].

The CFD-DEM simulations of the evaporating droplet show a strong dependency between the final granule morphology and the parameters of the bond model, which can be phenomenologically linked to the properties and concentration of a binder material in the spray suspension. It was shown that strong and early-forming bonds between the primary particles lead to the formation of hollow granules with thin crusts, while intermediate and weak bonds lead to hollow granules with thicker crusts and densely packed spheres, respectively. Furthermore, the simulation framework has been proved a valuable tool for relating the drying conditions that the droplets experience in a spray dryer to the average final granule morphology. It was shown that the drying rate influences both the internal circulation of the primary particles inside the droplet, as well as the center-directed capillary force from the droplet surface, which in turn leads to different particle morphologies. The main findings can be summarized as follows:

• (1). Simulations and experiments on different scales are combined to describe and investigate the particle formation mechanisms in spray dryers, both on a droplet and apparatus scale.

• (2). The CFD-DEM model is thoroughly analyzed for its physical plausibility, and the impact of the main model input parameters has been determined. The internal flow field and the droplet deformation due to the surrounding gas flow show comparable results to analytical and empirical results from the literature.

• (3). The droplet simulations illustrate the complex behavior of particle–particle–liquid movement and interactions, which influence the particle formation process.

• (4). The results from the droplet-scale simulation show that the morphology of spray dried particles is influenced both by the effective drying conditions and the material properties.

In future work, the parameters of the bond model could be calibrated by performing controlled experiments in an acoustic levitator [30,35] or by drying single falling droplets [36] and by comparing the resulting granule morphology with the simulation results. By using different types and/or amounts of binder, the authors could aim to find a correlation with the parameters of the bond model and, potentially, derive an empirical closure that will simplify future applications of the CFD-DEM simulation framework.

Author Contributions

Conceptualization, M.B., D.W., R.T., C.G. and S.H.; methodology, M.B., D.W., R.T., C.G. and S.H.; software, R.T. and C.G.; validation, M.B.; formal analysis, M.B.; investigation, M.B., D.W. and R.T.; data curation, M.B.; writing—original draft preparation, M.B., D.W. and R.T.; writing—review and editing, M.B., D.W., R.T., C.G. and S.H.; visualization, M.B. and D.W.; supervision, C.G. and S.H.; project administration, C.G. and S.H.; funding acquisition, C.G. and S.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Acknowledgments

Parts of this work are reprinted from [19] with permission from Cuvillier Verlag, Göttingen.

Conflicts of Interest

Authors Moritz Buchholz and Dominik Weis were employed by the company “BASF SE“. Authors Riccardo Togni and Christoph Goniva were employed by the company “DCS Computing GmbH“. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The authors declare that this study received funding from BASF SE. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Schematic representation of the (a) macroscopic CFD and (b) microscopic single-droplet CFD-DEM simulations.

Figure 2. (a) Calculation domain of CFD-DEM simulation and interactions between primary particles, the liquid and surrounding gas to describe the drying process of a single suspension droplet in a hot gas environment. (b) Schematic representation of the liquid volume fraction field [Forumla omitted. See PDF.] inside the droplet (dark blue), outside the droplet (white) and at the interface (light blue).

Figure 3. Schematic illustration of the relevant properties determining the capillary force, according to Breinlinger et al. [11].

Figure 4. Schematic illustration of the bond creation at a certain overlap of two particles and the breakage mechanism between two particles.

Figure 5. Simulation setup for the analysis of the modeling approach.

Figure 6. d2-laws of the single-droplet simulation in comparison with two single-droplet drying experiments at 60 °C and 140 °C [30]. Simulation 1 and 2 have been performed with evaporation coefficients of 8300 and 31,800 µm2/s, respectively.

Figure 7. Simulation results for the velocity profile of a single droplet moving in a gas environment. The vectors are shown in the center cut plane through the droplet. In the simulations, the droplet has a diameter of 2 mm. The inlet gas velocity for the simulation scenario is [Forumla omitted. See PDF.] and the gas temperature is [Forumla omitted. See PDF.].

Figure 8. Comparison of the aspect ratio between the single-droplet simulation and the empirical correlation by Hsiang & Faeth [33]. In the simulations, the droplet has a diameter of 2 mm. The inlet gas velocity for the simulation scenarios is [Forumla omitted. See PDF.], and the gas temperature is [Forumla omitted. See PDF.].

Figure 9. Time series of simulation results for two simulations with different bond parameters and gas inlet temperatures: (a) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] at 80 °C gas inlet temperature, (b) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] at 170 °C gas inlet temperature. The images show a vertical cut through the center of the droplet. The particles that are further away from the center plane are colored in darker grey. The formed bonds are colored in white.

Figure 10. Change in the MCN over time for two values of the bond parameter [Forumla omitted. See PDF.], leading to a dense spherical particle in the case of [Forumla omitted. See PDF.] for a gas inlet temperature of 80 °C (blue line), and to a hollow particle in the case of [Forumla omitted. See PDF.] for a gas inlet temperature of 170 °C (green line).

Figure 11. Morphology map of solidified suspension droplets for different pairs of bond parameters. The images show a vertical cut through the center of the droplet. The particles that are further away from the center plane are colored in darker grey. The formed bonds are colored in white.

Figure 12. Mean evaporation coefficients [Forumla omitted. See PDF.] and their respective standard deviations according to the [Forumla omitted. See PDF.]-law for different droplet size classes of a CFD simulation of a Niro Minor spray dryer experiment [20].

Figure 13. Morphology map of solidified suspension droplets for different drying rates, indicated by the evaporation coefficient [Forumla omitted. See PDF.]. The images show a vertical cut through the center of the droplet. The particles that are further away from the center plane are colored in darker grey.

Appendix A. Model Verification

The CFD model (Equations (1)–(4)) has been verified by performing a simulation of a suspended droplet of 1 mm radius undergoing evaporation over 40 s; see Figure A1. The inlet velocity is 1 cm/s and the inlet gas temperature is 360 K. The initial temperature of the system, i.e., the droplet and surrounding air, is 300 K. The mass and energy budgets for this simulation are represented in Figure A2. As can be seen in Figure A2a, the decrease in liquid mass (magenta line) perfectly balances the gas outflow (green line), while all other terms are negligible. The energy budget, displayed in Figure A2b, shows that the system is cooling down over time (see negative value of the internal energy variation shown by the blue line). This occurs because the thermal energy sink due to the evaporative cooling (see green line) is higher than the heat entering the system from the bottom inflow (see magenta line).

Figure A1. CFD simulation of an evaporating droplet. Vertical cross-section of the domain colored by (a) the liquid volume fraction [Forumla omitted. See PDF.] and (b) the temperature. The arrows denote the velocity field, while the white line represents the droplet interface for [Forumla omitted. See PDF.].

Figure A2. (a) Mass and (b) energy budgets from the CFD simulation of an evaporating droplet.

The capillary force model (see Equation (11)) has been verified by measuring the velocity jump of a solid particle crossing the fluid–gas interface for different contact angles $\mathsf{\Theta }$. In the absence of other external forces like gravity, the change in kinetic energy that a solid particle experiences when crossing the fluid–gas interface can be calculated as the integral of the capillary force along the particle trajectory, see Equation (13). Hence, the theoretical final velocity can be calculated from the initial kinetic energy and the work performed by the capillary force as follows:(A1)$$v_f=\sqrt{v_0^2+\frac{16\pi \sigma r^2\cos \mathsf{\Theta }}{3m}}$$

Figure A3 shows the particle velocity as a function of time for different values of the contact angle. The CFD-DEM simulations have been performed in the setup shown in Figure A4. The particle, which is initially immersed in water and has a velocity of 1 m/s, begins to cross the interface to air after 1 μs. Depending on the value of $\mathsf{\Theta }$, different behaviors are observed. For very hydrophilic solid materials (e.g., $\mathsf{\Theta }=0$ and $30°$), the particle’s velocity rapidly increases as the interface is crossed. For a less hydrophilic material ($\mathsf{\Theta }=60°$), the particle first decelerates and then accelerates again as its center crosses the interface. For a neutral material ($\mathsf{\Theta }=90°$), the particle does not cross the interface and bounces back. The final velocities obtained from the DEM simulation correspond to the theoretical values obtained from Equation (A1). A similar behavior can be seen in Figure A3b for an initial particle velocity of 5 m/s.

Figure A3. Particle velocity as a function of time obtained from the CFD-DEM simulation (solid lines), where the particle has an initial velocity of (a) 1 m/s and (b) 5 m/s. The dashed lines denote the particle’s final velocity derived analytically from Equation (A1).

Figure A4. CFD-DEM simulation of a solid particle (radius equal to 10 μm) crossing the interface from water (red) to air (blue) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]; (a) before and (b) after crossing the interface.

Appendix B. Simulation Parameters

Table A1

Overview of model parameters for different simulation cases.

Parameter	Model Analysis(Section 3.1)	Particle Morphologies(Section 3.2)	Drying Rate(Section 3.3)
Boundary/initial conditions
$u_{in}$ [m/s]	0.1	0.01	0.01
$T_{in}$ [K]	360	353	343
$u_0$ [µm]	2000	200	25/40
Material properties
${\rho }_w$ [kg/m3]	1000	1000	1000
$c_{p,w}$ [J/(kg K)]	4186	4186	4186
$k_w$ [W/(m K)]	0.6089	0.6089	0.6089
${\sigma }_w$ [N/m]	1 × 10−5	0.07	0.07
${\rho }_{air}$ [kg/m3]	1	1	1
$c_{p,air}$ [J/(kg K)]	1852	1852	1852
$k_{air}$ [W/(m K)]	0.04353	0.04353	0.04353
$\mathsf{\Delta }{h}_v$ [kJ/kg]	2256	2256	2256
$d_p$ [µm]	10	4	1
$E$ [Pa]	5 × 106	5 × 105	1 × 107
$\nu$ [-]	0.45	0.45	0.45
Numerical Parameters
$\mathsf{\Delta }{t}_{CFD}$ [s]	1 × 10−3	5 × 10−6	2.5 × 10−7
$\mathsf{\Delta }{t}_{DEM}$ [s]	1 × 10−5	2.5 × 10−7	2.5 × 10−9
$N_p$ [-]	2000	11,400	1420/5820
$C_{evap}$ [kg/(K s)]	−1 × 10−10	−1 × 10−14	−5 × 10−16/−2 × 10−15/−4 × 10−15

References

1. Tsotsas, E. Multiscale Approaches to Processes That Combine Drying with Particle Formation. Dry. Technol.; 2015; 33, pp. 1859-1871. [DOI: https://dx.doi.org/10.1080/07373937.2015.1047954]

2. Mezhericher, M.; Levy, A.; Borde, I. Theoretical Models of Single Droplet Drying Kinetics: A Review. Dry. Technol.; 2010; 28, pp. 278-293. [DOI: https://dx.doi.org/10.1080/07373930903530337]

3. Law, C.K. Recent advances in droplet vaporization and combustion. Prog. Energy Combust. Sci.; 1982; 8, pp. 171-201. [DOI: https://dx.doi.org/10.1016/0360-1285(82)90011-9]

4. Na, K.; Choi, M.; Ryoo, R. Recent advances in the synthesis of hierarchically nanoporous zeolites. Microporous Mesoporous Mater.; 2013; 166, pp. 3-19. [DOI: https://dx.doi.org/10.1016/j.micromeso.2012.03.054]

5. Sosnik, A.; Seremeta, K.P. Advantages and challenges of the spray-drying technology for the production of pure drug particles and drug-loaded polymeric carriers. Adv. Colloid Interface Sci.; 2015; 223, pp. 40-54. [DOI: https://dx.doi.org/10.1016/j.cis.2015.05.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26043877]

6. Kloss, C.; Goniva, C.; Hager, A.; Amberger, S.; Pirker, S. Models, Algorithms and Validation for Opensource DEM and CFD-DEM. Prog. Comput. Fluid Dyn. Int. J.; 2012; 12, 140. [DOI: https://dx.doi.org/10.1504/PCFD.2012.047457]

7. Zhu, H.P.; Zhou, Z.Y.; Yang, R.Y.; Yu, A.B. Discrete Particle Simulation of Particulate Systems: Theoretical Developments. Chem. Eng. Sci.; 2007; 62, pp. 3378-3396. [DOI: https://dx.doi.org/10.1016/j.ces.2006.12.089]

8. Zhu, H.P.; Zhou, Z.Y.; Yang, R.Y.; Yu, A.B. Discrete Particle Simulation of Particulate Systems: A Review of Major Applications and Findings. Chem. Eng. Sci.; 2008; 63, pp. 5728-5770. [DOI: https://dx.doi.org/10.1016/j.ces.2008.08.006]

9. El-Emam, M.A.; Zhou, L.; Shi, W.; Han, C.; Bai, L.; Agarwal, R. Theories and applications of CFD–DEM coupling approach for granular flow: A review. Arch. Comput. Methods Eng.; 2021; 28, pp. 4979-5020. [DOI: https://dx.doi.org/10.1007/s11831-021-09568-9]

10. Golshan, S.; Sotudeh-Gharebagh, R.; Zarghami, R.; Mostoufi, N.; Blais, B.; Kuipers, J.A.M. Review and Implementation of CFD-DEM Applied to Chemical Process Systems. Chem. Eng. Sci.; 2020; 221, 115646. [DOI: https://dx.doi.org/10.1016/j.ces.2020.115646]

11. Breinlinger, T.; Hashibon, A.; Kraft, T. Simulation of the Influence of Surface Tension on Granule Morphology during Spray Drying Using a Simple Capillary Force Model. Powder Technol.; 2015; 283, pp. 1-8. [DOI: https://dx.doi.org/10.1016/j.powtec.2015.05.009]

12. Breinlinger, T.; Hashibon, A.; Kraft, T. Simulation of the Spray Drying of Single Granules: The Correlation Between Microscopic Forces and Granule Morphology. J. Am. Ceram. Soc.; 2015; 98, pp. 1778-1786. [DOI: https://dx.doi.org/10.1111/jace.13557]

13. Nishiura, D.; Shimosaka, A.; Shirakawa, Y.; Hidaka, J. Simulation of Drying of Particulate Suspensions in Spray-Drying Granulation Process. J. Chem. Eng. Jpn.; 2010; 43, pp. 641-649. [DOI: https://dx.doi.org/10.1252/jcej.10we040]

14. Yabe, T.; Xiao, F.; Utsumi, T. The Constrained Interpolation Profile Method for Multiphase Analysis. J. Comput. Phys.; 2001; 169, pp. 556-593. [DOI: https://dx.doi.org/10.1006/jcph.2000.6625]

15. Miyazaki, S.; Nishiura, D.; Shimosaka, A.; Shirakawa, Y.; Hidaka, J. Revealing the Formation Mechanism of Granules by Drying Simulation of Slurry Droplet. Adv. Powder Technol.; 2011; 22, pp. 93-101. [DOI: https://dx.doi.org/10.1016/j.apt.2010.12.007]

16. Washino, K.; Tan, H.S.; Hounslow, M.J.; Salman, A.D. A New Capillary Force Model Implemented in Micro-Scale CFD–DEM Coupling for Wet Granulation. Chem. Eng. Sci.; 2013; 93, pp. 197-205. [DOI: https://dx.doi.org/10.1016/j.ces.2013.02.006]

17. OpenCFD Ltd. OpenFOAM—The Open Source CFD Toolbox. 2009; Available online: http://www.openfoam.com (accessed on 20 January 2024).

18. Xie, Y.; Chen, Y.; Fang, Z.; Zhou, H.; Wei, S.; Yang, L. The research of gas-solid fluidized bed bubbling behavior based on CFD-DEM coupled simulation. Chem. Eng. Res. Des.; 2023; 195, pp. 166-180. [DOI: https://dx.doi.org/10.1016/j.cherd.2023.05.053]

19. Buchholz, M. Towards a Multiscale Model-Based Design of Spray Drying Processes; Granulation; Number 23 in SPE-Schriftenreihe Cuvillier Verlag: Göttingen, Germany, 2023.

20. Buchholz, M.; Haus, J.; Pietsch-Braune, S.; Kleine Jäger, F.; Heinrich, S. CFD-Aided Population Balance Modeling of a Spray Drying Process. Adv. Powder Technol.; 2022; 33, 103636. [DOI: https://dx.doi.org/10.1016/j.apt.2022.103636]

21. Zhou, Z.Y.; Kuang, S.B.; Chu, K.W.; Yu, A.B. Discrete particle simulation of particle–fluid flow: Model formulations and their applicability. J. Fluid Mech.; 2010; 661, pp. 482-510. [DOI: https://dx.doi.org/10.1017/S002211201000306X]

22. Gopala, V.R.; van Wachem, B.G.M. Volume of fluid methods for immiscible-fluid and free-surface flows. J. Chem. Eng.; 2008; 141, pp. 204-221. [DOI: https://dx.doi.org/10.1016/j.cej.2007.12.035]

23. Koch, D.L.; Hill, R.J. Inertial Effects in Suspension and Porous-Media Flows. Annu. Rev. Fluid Mech.; 2001; 33, pp. 619-647. [DOI: https://dx.doi.org/10.1146/annurev.fluid.33.1.619]

24. Cundall, P.A.; Strack, O.D.L. A Discrete Numerical Model for Granular Assemblies. Géotechnique; 1979; 29, pp. 47-65. [DOI: https://dx.doi.org/10.1680/geot.1979.29.1.47]

25. Di Renzo, A.; Di Maio, F.P. Comparison of Contact-Force Models for the Simulation of Collisions in DEM-Based Granular Flow Codes. Chem. Eng. Sci.; 2004; 59, pp. 525-541. [DOI: https://dx.doi.org/10.1016/j.ces.2003.09.037]

26. Pöschel, T.; Schwager, T. Computational Granular Dynamics; Springer: Berlin/Heidelberg, Germany, 2005.

27. Liggghts^® Open Source Discrete Element Method Particle Simulation Code. Available online: https://www.cfdem.com/liggghts-open-source-discrete-element-method-particle-simulation-code (accessed on 20 January 2024).

28. Potyondy, D.O.; Cundall, P.A. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci.; 2004; 41, pp. 1329-1364. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2004.09.011]

29. Tuckermann, R.; Bauerecker, S.; Neidhart, B. Evaporation Rates of Alkanes and Alkanols from Acoustically Levitated Drops. Anal. Bioanal. Chem.; 2002; 372, pp. 122-127. [DOI: https://dx.doi.org/10.1007/s00216-001-1132-7]

30. Buchholz, M.; Haus, J.; Polt, F.; Pietsch, S.; Schönherr, M.; Kleine Jäger, F.; Heinrich, S. Dynamic Model Development Based on Experimental Investigations of Acoustically Levitated Suspension Droplets. Int. J. Heat Mass Transf.; 2021; 171, 121057. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121057]

31. Kallendorf, C.; Fath, A.; Oberlack, M.; Wang, Y. Exact Solutions to the Interfacial Surfactant Transport Equation on a Droplet in a Stokes Flow Regime. Phys. Fluids; 2015; 27, 082104. [DOI: https://dx.doi.org/10.1063/1.4928547]

32. Rybczynski, W. On the Translatory Motion of a Fluid Sphere in a Viscous Medium. Bull. Acad. Sci. Crac. Ser. A; 1911; 40, pp. 33-78.

33. Hsiang, L.-P.; Faeth, G.M. Near-Limit Drop Deformation and Secondary Breakup. Int. J. Multiph. Flow; 1992; 18, pp. 635-652. [DOI: https://dx.doi.org/10.1016/0301-9322(92)90036-G]

34. Walton, D.E.; Mumford, C.J. Spray Dried Products—Characterization of Particle Morphology. Chem. Eng. Res. Des.; 1999; 77, pp. 21-38. [DOI: https://dx.doi.org/10.1205/026387699525846]

35. Mondragon, R.; Julia, J.E.; Barba, A.; Jarque, J.C. Microstructure and mechanical properties of grains of silica nanofluids dried in an acoustic levitator. J. Eur. Ceram. Soc.; 2012; 32, pp. 4295-4304. [DOI: https://dx.doi.org/10.1016/j.jeurceramsoc.2012.07.017]

36. Vehring, R.; Foss, W.R.; Lechuga-Ballesteros, D. Particle formation in spray drying. J. Aerosol Sci.; 2007; 38, pp. 728-746. [DOI: https://dx.doi.org/10.1016/j.jaerosci.2007.04.005]