A thick wall concept for robust treatment of contacts in DEM simulation of highly polydisperse particulate systems

Modelling particulate systems with the Discrete Element Method (DEM) is an established practice, both in the representation and analysis of natural phenomena and in scale-up and optimization of industrial processes. Since the method allows tracking individual particles, each element can possess geometrical, physical, mechanical or chemical surface properties different from those of the other particles. One example is a polydisperse particulate system, i.e., characterized by a size distribution, opposed to the idealized monodisperse case. In conventional DEM, a softer particle stiffness is commonly adopted to reduce the computational time. It might happen that artificially soft particles, when colliding against a wall boundary, exhibit such large, unrealistic overlap that they “pass through” the wall and exit the domain. In the case of highly polydisperse systems, this often occurs when fine particles are pushed against the wall by coarse particles with masses several orders of magnitude larger. In the manuscript, a novel method is proposed, named thick wall, to allow the particles in contact with the walls to experience relatively large overlaps without ending up ejected out the domain. In particular, a careful way to calculate the particle-wall overlap and force unit vector can accommodate normal displacements larger than the maximum usually allowed, i.e., typically the particle radius, thereby preventing particles from being expelled from the domain. First, critical velocities for which single particles and pairs of fine/coarse particle escape the domain are analytically characterized using the linear and the Hertz models. The thick wall concept is then introduced and its effect on the maximum critical velocity is demonstrated with both contact models. Finally, application to pharmaceutical powder composed of carrier (coarse) and active pharmaceutical ingredient (API) (fine) particles in a shaken capsule prove this to be an example of vulnerability to the phenomenon of fine particle ejection and to significantly benefit from the thick wall modification.

1 Introduction

Particulate material processing often involves dealing with some degree of polydispersity, e.g., the particle size inherently follows a distribution, or the process requires particles of different sizes to be mixed for improved efficiency. Despite the importance of powder sampling and preparation, sieving and other classification processes are always imperfect. In binary or multi-solid mixtures small-size particles arise naturally, which may then exhibit segregation and final product inhomogeneity. Also, fines may be generated by abrasion of coarse brittle materials. If the small particles are sufficiently fine, they are subjected to adhesive forces and are often found attached to coarser particles, potentially detaching as a result of collisions with the containing walls or by the shearing action of a fluid stream. Examples include formation of dust during bulk solids handling (Schulz et al., 2022; Schulz et al., 2023), the deaggregation of active pharmaceutical ingredient (API) powder from coarser carrier particles in dry powder inhalation devices (Sharma and Setia, 2019; Chaurasiya and Zhao, 2020; Spahn et al., 2022), the production of active material with a layer of conducting material (e.g., carbon black) during dry mixing in battery electrodes manufacturing (Lischka et al., 2024), the use of solid flow-aids to improve powder flowability or to track powder mixing (Fulchini et al., 2017; Karde et al., 2023; Khala et al., 2023) and the improvement of powder bed homogeneity in additive manufacturing of ceramic materials by using bimodal particle size distribution (Sofia et al., 2018; Schmidt and Peukert, 2022).

The Discrete Element Method (DEM), also in combination with computational fluid dynamics (CFD-DEM), is commonly utilized to model the behavior of granular materials, due to its versatility in dealing with very different flow regimes, packing density and dispersion of particle properties (Golshan et al., 2020; Kieckhefen et al., 2020). DEM has been successfully used to study the segregation of particles differing in size and density in fluidized beds (Di Renzo et al., 2011; Peng et al., 2016), the mixing of electrically charged API and carrier particles for inhalation in a vibrating container (Yang et al., 2015), the mixing of binary particle beds in a rotating drum (Huang and Nakagawa, 2023), the size segregation during powder spreading for additive manufacturing using laser powder bed fusion (Shaheen et al., 2021), the detachment of API from carrier particles in dry powder inhalers (Cui and Sommerfeld, 2015; Tong et al., 2017; Ariane et al., 2018; Alfano et al., 2021a; Alfano et al., 2022a; Alfano et al., 2022b).

Particles trajectories are solved for each solid element, including during interparticle contacts and impacts with walls, by considering the compression and rebounding stages. Indeed, contact forces are generally modelled using the soft-sphere model, which assumes the particles to keep their spherical shape and bases the force computation upon the calculation of the particle overlap. Linear spring-dashpot or Hertzian-type spring with some dissipation mechanisms is adopted in the direction normal to the contact plane, with similar mechanical models plus surface sliding in the tangential direction (Zhu and Yu, 2006; Zhu et al., 2007). In a collision, the material parameters like the normal spring stiffness or (Young’s) modulus of elasticity primarily determine the particle maximum deformation and the contact duration, the latter being responsible also for the maximum integration time-step (Kruggel-Emden et al., 2010). The realistic interparticle or particle-wall contact time is orders of magnitude smaller than the typical characteristic processing times. So, very commonly, softer than realistic stiffness constants are set, even two or more orders of magnitude smaller. Systematic ways to reduce the stiffness have been proposed to minimize its influence on the intended simulated quantities (Hærvig et al., 2017; Washino et al., 2018; Washino et al., 2024; He et al., 2021). For cohesionless particles, the differences in the contact dynamics does not affect, if not marginally, the kinematics (Thornton et al., 2011), so for moderate to granular flows with little inertial effects the evolution of the particles trajectories is unaffected (see e.g., in Kuo et al., 2002; Grima and Wypych, 2011). For cohesive particles, special treatments have been proposed, such as in (Washino et al., 2018; Washino et al., 2024). On the other hand, the larger integration time-step and reduced simulation time greatly benefit the feasibility of the intended DEM study. In addition to particle stiffness, the number of particles is another parameter significantly affecting the computational time required in DEM simulations. Depending on the algorithm used, the computational time may have a quadratic (worst case) or quasi-linear (best-case) dependence on the number of particles. Strongly polydisperse systems are often characterized by a high number of fine particles, thus significantly impacting computation times.

Given the particle properties, a reduced material stiffness can be computed that allows for the quickest simulations of the particle flow without incurring in excessive particle overlap. The maximum overlap depends on the impact velocity, so conditions may arise leading to localized high overlaps. In monodisperse particle systems, the critical velocity leading to overlaps comparable to the particle radius are very rarely encountered. However, simulations of polydisperse systems are much more vulnerable in this regard. It will be shown in the following subsections that in the case of particle of significantly different sizes, the expected interparticle or particle-wall overlap becomes overly exaggerated at much lower impact velocities, up to the point that particles can cross the wall boundaries and leave the geometric domain. An obvious solution would be to use a larger value of the material stiffness, so that the maximum overlap decreases to acceptable values. The required increase can be rather significant, with the consequence that the contact duration decreases correspondingly, and the required smaller time-step leads to much higher simulation times.

The objective of this work is to characterize the conditions of size difference, material properties and operating conditions leading to the above-referenced unrealistic overlap and to propose and test a robust calculation method to overcome such limits, without sacrificing the advantageous increased time-step. The first part, characterization of the conditions for unrealistic overlaps, will help understand and predict what particulate systems (e.g., highly polydisperse, or cohesive fine-on-coarse particle systems) may exhibit vulnerability to the phenomenon. So, the conditions leading to the unfavorable phenomenon are discussed and characterized first analytically for mono- and poly-disperse systems, in terms of critical impact velocity (the maximum velocity before particle crossing the wall occurs) as a function of size difference and material properties. Then, a modified and more robust calculation of the overlap is proposed, which aims to prevent particles from crossing the boundaries even when high overlaps are expected, yet without requiring changes in the integration time-step. Such a solution would prevent the onset of fine particles ejection out of the system, at essentially no additional computational cost. Expected limitations of the solution are discussed. Finally, an application of the DEM-CFD model to the flow of a fine-coarse pharmaceutical mixture in a capsule is presented and discussed, comparing the results of the conventional method with those of the proposed method.

2 DEM model of contacts with walls

We start by recalling the basic elements of the contact model in DEM simulations. The collision of two perfectly spherical particles and of one particle with a flat wall are considered for simplicity, although the subsequent treatment is amenable to other particle or wall shapes without significant changes. Also, for the sake of simplicity, only normal elastic contact will be considered in the mathematical derivation. This will include both linear spring and Hertzian contact models. For an analytical treatment of the elastic-dissipative contacts the interested reader is referred to Antypov and Elliott (2011).

In soft-sphere DEM simulations, the particles are modelled in their motion as rigid (i.e., non-deformable) spheres and walls as rigid flat surfaces. When the distance of the particle center from the wall is less than the particle radius, a contact is identified. The spherical shape is also used to compute the particle deformation during the collision, measuring it as the particle to surface or inter-particle overlap (Figure 1). According to the displacement-driven formulation, at each time instant, the elastic repulsive force is computed from the overlap, by means of geometrical and stiffness material parameters, followed by integration of Newton’s second law of dynamics for the new velocity and new position, ready for the next time integration. Such scheme allows tracking the instantaneous particle positions, deformation (overlap), velocity and forces during the collisions.

FIGURE 1

FIGURE 1. Particle in contact with a wall (A) or with another particle (B); comparison between actually deformed contacting shapes and DEM-modelled undeformed representation, showing the overlap concept.

As anticipated, we are concerned with the possibility that particles undergo excessive overlaps, with the somewhat unexpected result that they may definitely cross the wall boundaries. As will be shown for the particles with different sizes, this phenomenon has severe consequences, much more often in the case of particle-wall collisions (Figure 2) than in particle-particle collisions. There are two reasons for that: 1) the inter-particle contact overlap must reach the sum of the two radii ($R_1+R_2$) in order for the particles to cross each other and, no matter what the force value becomes, 2) the particles remain inside the domain. On the other hand, for particle-wall contacts, the particle to wall critical overlap is the particle radius ($R$) and, once it is overcome, the particle undergoes a reversed overlap, getting repulsed in the other direction, i.e., further across the wall (see the blue overlap and the wall force as red arrows in Figure 2), until escaping the simulation domain and then becoming unrecoverable. In the classical algorithm, this is essentially due to the calculation formula for the overlap and the reversal of the normal unit vector pointing from the particle center to the wall surface, whose versus change as soon as the particle center moves beyond the wall surface.

FIGURE 2

FIGURE 2. Collisions giving rise to very high modelled overlap, showing the wall contact force with a red arrow and the overlap with a blue line. (A) high-speed impact of a single particle against a wall shown with actual deformation (left) and undeformed DEM model (right); (B) wall impact of two particles with the coarse one (2) strongly pushing, through its inertia, the fine one (1) towards the wall. Under certain high-loading conditions (as depicted for the undeformed two spheres), the calculated overlap “reverses” (see accompanying text) and the wall force is directed towards the wall rather than against it. Deformations are exaggerated for illustration purposes.

In a typical implementation with the linear spring model, the elastic force is expressed by $\mathit{F}=-k\delta {\mathit{e}}_n$, in which $k$ is the normal spring stiffness constant, $\delta$ is the overlap and the unit vector is ${\mathit{e}}_n=\frac{\mathit{d}}{\left|\mathit{d}\right|}$, where $\mathit{d}$ is the distance vector pointing from the particle center to the closest point on the wall surface. The overlap is calculated, along the normal vector connecting the particle center and the wall closest point, by $\delta =R-\left|\mathit{d}\right|$ and $R$ is the particle radius (Figure 3).

FIGURE 3

FIGURE 3. Schematic representation of the variables involved in the calculation of the wall normal force (see accompanying text for the variables’ definition and role in the contact model formulas). (A) Common low overlap conditions; (B) very high deformation case, with wrong calculation of the overlap and unit vector; (C) reversed wall contact conditions equivalent to scheme (B).

For the purpose of this study, contact conditions can be distinguished in two cases: common, low-overlap contacts, in which the deformation is relatively low; the distance vector points toward the wall, the overlap represents the actual total deformation, and the resulting force is wall repulsive (Figure 3A); very high deformation case, in which the particle center overcomes the wall surface; the overlap does not represent the total deformation and the distance vector points towards the internal domain, with the consequence that the force is wall attractive (Figure 3B). The second case is equivalent to a reversed wall contact (Figure 3C): the unit vector changes direction, and so does the force, and the overlap corresponds to the deformation of a wall on the other side of the domain; since the process diverges the particle will end up ejected.

Therefore, the focus of what follows will be on particle-wall contacts and the relevant critical conditions. Conditions leading to the critical point where the undeformed overlap reaches the particle radius are analytically derived for the case of a single particle impacting a flat wall and for a coarse/fine two-particle system impacting a flat wall (Figure 2). Both cases are investigated assuming a linear spring model and the Hertzian contact theory.

2.1 Single particle contact: linear spring model

Key results for the current analysis are the evolution of the particle overlap and the collision duration, which will be discussed below for the different cases. As anticipated, the focus will be on the purely elastic force in the wall normal direction, neglecting the effects of other contact forces (e.g., plastic or viscoelastic components). Gravity force is neglected as well, since it is several orders of magnitude smaller than contact forces for the particle sizes considered in this study.

Assuming that a particle approaches the wall with an initial velocity $v_0$ and the analysis starts from the moment when it touches the wall (${\delta }_0=0$), the following linear ODE governs the dynamics for all the contact duration:

$\frac{d^2\delta }{d{t}^2}+{\omega }^2\delta =0(1)$

in which $\omega =\sqrt{k/m}$ is the oscillation pulsation, and $m=\rho \frac{4}{3}\pi R^3$ is the particle mass and $\rho$ its density. As discussed in more detail in Di Maio and Di Renzo (2004), the solution is

$\delta \left(t\right)=\frac{v_0}{\omega }\sin \left(\omega t\right)(2)$

and the collision duration is $\tau =\pi /\omega$. Since the contact is a perfect sinusoidal oscillation, the maximum overlap occurs at half of the oscillation duration, i.e.,

${\delta }_{\mathit{max}}=\frac{v_0}{\omega }(3)$

Before proceeding further, it is useful to adopt a dimensionless formulation, so that general relationships will be more easily identified. The following variable changes are set:

${\delta }^{*}=\frac{\delta }{R};t^{*}=\frac{t}{\tau }(4)$

The ODE can be re-expressed in the following terms:

$\frac{d^2{\delta }^{*}}{d{t^{*}}^2}+{\pi }^2{\delta }^{*}=0(5)$

and the initial conditions are ${\delta }_0^{*}=0$ and $\frac{d{\delta }^{*}}{d{t}^{*}}=\pi \frac{v_0}{R\omega }$.

A characteristic dimensionless number can be defined, which is named here “impact number” $In$:

$In=\frac{v_0}{R\omega }(6)$

The dimensionless solution of Eq. 5 is

${\delta }^{*}\left(t^{*}\right)=In\sin \left(\pi t^{*}\right)(7)$

and ${\delta }_{\mathit{max}}^{*}=In$.

2.1.1 Linear model—critical impact conditions

The critical condition for our analysis is attained when the overlap reaches the particle radius, i.e., the critical, maximum overlap allowed before particle crossing the wall occurs, using the conventional calculation method described above. At that stage, the normal vector pointing from the particle center to the wall closest point reverses direction, and the overlap starts becoming computed in the opposite direction (Figure 3). As a result of this change, the repulsive force becomes directed in the opposite direction and the particle is accelerated towards the other side of the wall, i.e., it gets ejected out of the domain. It is noteworthy that no reduction in the time-step or change in the integration algorithm can solve this specific issue, as it is inherent in the approximation of undeformed bodies (spherical particles and flat walls), and unrelated to numerical accuracy of the integration scheme. Only using stiffer materials can mitigate the occurrence of particle ejection, but it comes at the cost of longer simulation times.

Analytically, this critical condition is easy to determine, by simply setting ${\delta }_{\mathit{max}}=R$ or ${\delta }_{\mathit{max}}^{*}=1$. The combination of parameters leading to the critical condition is expressed by

$I{n}_c=\frac{v_{0c}}{R\omega }=1(8)$

The critical impact velocity $v_{0c}$ can be derived from Eq. 8 and expressed as a function of the geometrical, physical and mechanical properties of the particles:

$v_{0c,L}=R\omega =R\sqrt{\frac{k}{m}}=\sqrt{\frac{3}{4}\frac{k}{\pi \rho R}}(9)$

where the subscript $L$ denotes the linear model.

Given the particle and wall properties, Eq. 9 allows determining the maximum impact velocity that a wall can withstand without letting the particle trespass it. Note that the purely elastic collision is the most extreme case since there is the minimum critical normal impact velocity and the maximum overlap. By introducing the contributions of dissipative forces, the maximum overlap decreases and the critical velocity increases.

Values for the critical impact velocity are shown in Figure 4 as a function of different combinations of typical values of the relevant properties (radius, density, stiffness). A higher rigidity, with a larger stiffness constant value, leads to a higher critical velocity, with a slope 0.5 on the log-log plot (Figures 4A,B). The values for a light small particle, for example, $R=5\mu \mathrm{m}$ and $\rho =800\text{\hspace{0.17em}kg}/{\mathrm{m}}^3$, lays well above 100 m/s even for the smallest value of the stiffness constant ($k=500\mathrm{N}/\mathrm{m}$). Heavier particles, for example, with density typical of metals ($\rho >7000\text{\hspace{0.17em}kg}/{\mathrm{m}}^3$) show a critical velocity above 100 m/s for stiffness constant $k>2000\mathrm{N}/\mathrm{m}$. For particles of 1 mm of diameter ($R=500\mu \mathrm{m}$) the possibility to pass through the wall is observed at velocities ranging from below 10 m/s for heavy particles and spring stiffness lower than 2000 N/m, to above 100 m/s for lighter particles with a spring stiffness higher than 20,000 N/m (Figure 4B). The analysis in terms of the particle radius shows a decrease of the critical velocity with a slope −0.5 on the log-log plot. With a spring stiffness $k=1000\mathrm{N}/\mathrm{m}$, the critical velocities are in the range 10–200 m/s (Figure 4C). By increasing the spring stiffness to 10,000 N/m the critical velocity increases to $v_{0c}>100\mathrm{m}/\mathrm{s}$ for the heaviest particle considered with radius $R<30\mu \mathrm{m}$ or with radius up to $100\mu \mathrm{m}$, provided the density is similar to that of glass ($\rho \approx 2500\text{\hspace{0.17em}kg}/{\mathrm{m}}^3$) or lower (Figure 4D).

FIGURE 4

FIGURE 4. Critical velocity as a function of the elastic spring stiffness and particle density for radius $\mathit{R}=\mathrm{5}\mathrm{\mu }\mathrm{m}$ (A) and $\mathit{R}=\mathrm{500}\mathit{\mu }\mathrm{m}$ (B) and as a function of the particle radius and density for spring stiffness $\mathit{k}=\mathrm{1000}\mathrm{N}/\mathrm{m}$ (C) and $\mathit{k}=\mathrm{10000}\mathrm{N}/\mathrm{m}$ (D).

It is interesting to derive the relations for the critical impact velocity if the elastic stiffness is related to the real material properties, like the Young’s modulus of elasticity $E$. This parameter appears in the Hertz theory of the elastic contact between a spherical particle and a wall, which shows a non-linear dependence of the elastic repulsive force on the macroscopic displacement (i.e., the overlap). Considering the expression of the force as a function of the overlap in Hertz’s theory (see Eq. 12 below), an equivalently linear spring stiffness is $k_{LH}=\frac{4}{3}E\sqrt{R}\sqrt{\delta }$. If a reference value for the overlap equal to the particle radius—relevant to the very high deformation cases considered here—is adopted, a relation between the spring stiffness $k_{LH}$ and the Young’s modulus is obtained as

$k_{LH}=\frac{4}{3}ER(10)$

This can be substituted into Eq. 9 to give the following variant of the critical impact velocity expression

$v_{0c,L}=\sqrt{\frac{E}{\pi \rho }}(11)$

The interesting result is that the critical velocity becomes independent of the particle radius and changes only with the material modulus of elasticity and the particle density. As a curiosity, it is worth noting that $v_s=\sqrt{E/\rho }$ is the velocity of one-dimensional compressive waves in the solid, so the critical velocity happens to be related to this physically well-defined material parameter by a proportionality factor $1/\sqrt{\pi }$. The specific factor depends on the particular choice of the reference overlap adopted in deriving Eq. 10. However, the independence of the radius is obtained irrespective of this choice.

2.2 Single particle contact: Hertz model

The full consideration of mechanical normal stress-strain relation during a contact of homogeneous sphere with a flat wall was derived by Hertz and later extended to contacts between particles of different sizes and properties (Johnson, 1987). The well-known macroscopic elastic force-displacement relationship is

$F=\frac{4}{3}E\sqrt{R}{\delta }^{\frac{3}{2}}(12)$

The corresponding ODE governing the transient motion becomes non-linear

$\frac{d^2\delta }{d{t}^2}+\frac{E}{\rho \pi R^{5/2}}{\delta }^{\frac{3}{2}}=0(13)$

and cannot be easily solved analytically even for the purely elastic case. A few results are available (see e.g., Maw et al., 1976), such as the collision duration and the maximum displacement, which are sufficient for our purpose. In particular, the maximum displacement is

${\delta }_{\mathit{max}}={\left(\frac{15}{16}\frac{m}{\sqrt{R}E}{v}_0^2\right)}^{\frac{2}{5}}(14)$

By expressing the particle mass as function of the density and size, the maximum dimensionless overlap is easily obtained as

${\delta }_{\mathit{max}}^{*}=\frac{{\delta }_{\mathit{max}}}{R}={\left(\frac{5}{4}\pi \frac{\rho }{E}{v}_0^2\right)}^{\frac{2}{5}}(15)$

Interestingly, inside the parenthesis in Eq. 15, the square of the impact velocity $v_0$ is multiplied by a factor in which the square of the critical velocity using the linear model and the stiffness defined in terms of the Young’s modulus $v_{0c,L}=\sqrt{E/\pi \rho }$ (Eq. 11) appears.

2.2.1 Hertz model—critical impact conditions

Following the results of the Hertz theory of elastic contact applied to a collision, the critical condition, i.e., ${\delta }_{\mathit{max}}^{*}=1$, can be found from Eq. 15 as

$v_{0c,H}=\sqrt{\frac{4}{5}\frac{E}{\pi \rho }}(16)$

where the subscript $H$ denotes the Hertz model.

It is worth noting that the critical velocity is independent of the particle radius, but changes only with the material density and Young’s modulus.

Figure 5 shows the plots of the critical velocity according to the Hertz theory as a function of typical values of the Young’s modulus for a range of density values. Values above 200 m/s are easily obtained provided that the Young’s modulus is higher than 1 GPa, as typical of many solid materials. For stiffer systems, like steel ($E\approx 200\text{\hspace{0.17em}GPa}$), the critical velocity is well above 2000 m/s. On the other extreme, even for very deformable solids, for example, with $E=0.01\text{\hspace{0.17em}GPa}$, the critical velocity is always of the orders of tens of meters per second. All such velocity values can be considered rather high, such that impacts at velocities approaching the critical values can only be observed, if any, in fast flows.

FIGURE 5

FIGURE 5. Critical velocity according to the Hertz contact model as a function of the Young’s modulus of elasticity and particle density.

Also, it turns out slightly different from the value obtained with the linear model, even when the linear spring stiffness was related to the Young’s modulus (Eq. 10). In particular, the ratio of the critical velocities obtained with the two models is $\frac{v_{0c,H}}{v_{0c,L}}=\sqrt{\frac{4}{5}}\approx 0.89$.

3 DEM model of fine-coarse particle collision

In Section 2, the critical velocity leading to the particle center passing through the wall was calculated using the linear and the Hertz model. It is observed that the resulting velocities are generally high. In the present Section, cases in which relatively fine particles are pushed towards the wall by relatively coarser particles and the corresponding critical velocities are investigated. Like above, the analysis will be conducted separately using the linear model and the Hertz theory. The reference configuration is as depicted in Figure 2B, in which the particle close to the wall is denoted by 1 and the other one pushing it is denoted by 2.

3.1 Fine-coarse contact: linear spring model

For simplicity, it will be assumed that the two particles 1 (fine) and 2 (coarse) move integral with one another, as if they are kept together by a strong adhesive force. The analysis will be substantially valid also for the case in which the particle 1 moves parallel to the wall (e.g., rolling on it) and the particle 2 hits it perpendicularly, pushing it against the wall. The two particles will be assumed to have the same density, as the generalization to the case of different densities is trivial.

Under these assumptions, the inertia of both particles is to be repulsed by the containing walls, while the maximum allowed penetration equals the radius of particle 1, if such particle is to remain inside the boundary. The equation governing the motion of the colliding body is the same as in Section 2.1, but with different parameters. The oscillation pulsation becomes ${\omega }_{12}=\sqrt{k_1/\left(m_1+m_2\right)}$ , which for similar density solids leads to ${\omega }_{12}=\sqrt{k_1/\left(\frac{4}{3}\pi {\rho }_1\left(R_1^3+R_2^3\right)\right)}$

The collision duration is ${\tau }_{12}=\pi /{\omega }_{12}$. The maximum allowed overlap is $R_1$. So, the dimensionless variable will be considered as follows

${\delta }^{*}=\frac{\delta }{R_1};t^{*}=\frac{t}{{\tau }_{12}}(17)$

The dimensionless impact number is

$I{n}_{12}=\frac{v_{012}}{R_1{\omega }_{12}}(18)$

and, as before, the solution for the highest value of the dimensionless displacement is ${\delta }_{12,\mathit{max}}^{*}=I{n}_{12}$.

3.1.1 Linear model—critical impact conditions

The critical condition for a fine-coarse contact is achieved when the highest overlap reaches the radius of particle 1, i.e., when the dimensionless overlap is $I{n}_{12}=1$. Similar to the case of a single particle, the critical velocity can be computed by

$I{n}_{12c}=\frac{v_{0c12}}{R_1{\omega }_{12}}=1(19)$

yielding

$v_{0c12,L}=R_1\sqrt{\frac{3}{4}\frac{k_1}{\pi {\rho }_1\left(R_1^3+R_2^3\right)}}=\sqrt{\frac{3}{4}\frac{k_1}{\pi {\rho }_1{R}_1\left(1+q^3\right)}}(20)$

in which the particle size ratio $q=R_2/R_1$ has been introduced. The expression in Eq. 20 allows examining the critical velocity for any combination of the four parameters, $k_1$, $R_1$, ${\rho }_1$ and $q$.

It is interesting to compare the critical velocity value for the fine-coarse contact with that of a single particle contact, as their ratio depends only on the particle size ratio

$\frac{v_{0c12,L}}{v_{0c,L}}=\frac{1}{\sqrt{1+q^3}}(21)$

Figure 6 shows the plot of the critical velocity ratio of the fine-coarse contact over single particle contact as a function of the particle size ratio $q$, where 1 denotes the particle touching the wall and 2 the particle pushing the other one. It is observed that for low particle size ratio, the ratio of velocities is essentially unity. Only when $q$ approaches 1 does the velocity for the two-(equal-)size critical velocity become noticeably smaller than the single particle case. This is obvious, as the sum of masses of two equal particles leads to a different inertial behavior of the contact bodies. Indeed, for $q<1$ the pushing particle is smaller than the particle in contact and the critical velocity of the latter is essentially unaffected by the (less massive) pushing one. For coarser pushing particles ($q>1$), the decrease in critical velocity is much more evident. As the size ratio increases, the ratio of critical velocity follows a—3/2 sloped asymptote. This means that at a size ratio $q=5$ the critical velocity decreases by about one order of magnitude with respect to the single particle case. If the size ratio is 10, the critical velocity drops down to 3% of the value without the pushing particle and for $q=20$ it drops by about two orders of magnitude.

FIGURE 6

FIGURE 6. Critical impact velocity ratio (fine-coarse case over single particle contact), calculated according to the linear contact model as a function of the particle size ratio.

In some of the applications discussed in Section 1, the particle size ratio can reach up to 50, which corresponds to a decrease of the critical velocity to less than 0.3% of the single particle critical value, or even more.

Similar to the case for a single particle, let us use the definition of the contact stiffness in terms of the Young’s modulus (Eq. 10). The corresponding variant of the critical impact velocity is

$v_{0c12,LH}=\sqrt{\frac{E_1}{\pi {\rho }_1\left(1+q^3\right)}}(22)$

which yields the same decrease of the critical velocity ratio as that plotted in Figure 6. As noted earlier, the result in terms of the mechanical parameter $E_1$ is independent of the absolute particle radius.

3.2 Fine-coarse contact: Hertz model

The same approach used in the linear model is adopted here, by substituting the mass of the single particle with the total mass, $m_1+m_2$, and the maximum overlap as the radius of the particle 1.

The analytical expression of the maximum overlap is

${\delta }_{12,\mathit{max}}={\left(\frac{15}{16}\frac{m_1+m_2}{\sqrt{R_1}{E}_1}{v}_0^2\right)}^{\frac{2}{5}}(23)$

The dimensionless form of such overlap can be expresses as

${\delta }_{12,\mathit{max}}^{*}=\frac{{\delta }_{12,\mathit{max}}}{R_1}={\left[\frac{5}{4}\frac{\pi {\rho }_1}{E_1}\left(1+q^3\right)v_0^2\right]}^{\frac{2}{5}}(24)$

3.2.1 Hertz model: critical impact conditions

The critical velocity at which the overlap reaches the radius of the particle contacting the wall (1) as a function of the relevant variables is obtained by setting ${\delta }_{12,\mathit{max}}^{*}=1$, yielding

$v_{0c12,H}=\sqrt{\frac{4}{5}\frac{E_1}{\pi {\rho }_1\left(1+q^3\right)}}(25)$

By comparing the obtained value with that of a single particle, one obtains

$\frac{v_{0c12,H}}{v_{0c,H}}=\frac{1}{\sqrt{1+q^3}}(26)$

exactly like in the linear model case. This indicates that a decrease in the critical velocity by one or two orders of magnitude for size ratio $q=5$ or 20, respectively, is observed irrespective of the model considered.

4 A thick wall concept for particle-wall contacts

Under the conditions discussing above for the fine-coarse contact, the probability for the finer particles to be pushed out of the domain by the coarse particles is rather high, even in flows at moderate velocities. This evidently leads to reconsider the importance of a robust treatment of the wall contacts, if the objective is to avoid the risk of having fine particles ejected by coarse particles.

4.1 Robust solution of the normal contact

A simple modification of the overlap calculation is proposed, by conceptually adding a sort of solid “thickness” to the wall. In Section 2, we briefly discussed the conventional implementation of the particle-wall normal contact force. The key steps are the expression of the wall force exerted on the particle as a function of the overlap and the unit vector and the calculation of these last two variables. They are discussed here in relation to Figure 7 where a high-deformation contact is schematically shown with the variables involved in the conventional calculation of the force (Figure 7A) and in the calculation proposed according to the thick wall concept (Figure 7B).

FIGURE 7

FIGURE 7. Schematic representation of the high deformation particle-wall contact with the particle center beyond the wall boundary. (A) Conventional calculation of the overlap and the unit vector—see Eqs. 28, 29—for the contact force, which leads to a force directed towards the wall (the particle will be pushed out of the internal domain, across the wall); (B) robust calculation of the overlap and the unit vector according to the thick wall concept; the particle experiences a high overlap, corresponding to the very high deformation and the force direction is correctly pointing towards the internal domain (the particle will eventually rebound inside the domain).

The wall normal force is expressed by

$\mathit{F}=-k\delta {\mathit{e}}_n(27)$

The (scalar) overlap and the unit vector for the direction are computed, respectively, by

$\left.\left.\delta =R-\right|\mathit{d}\right|(28)$

${\mathit{e}}_n=\frac{\mathit{d}}{\left|\mathit{d}\right|}(29)$

in which $\mathit{d}$ represents the distance vector pointing from the particle center to the closest point on the wall surface.

As introduced in Section 2.1, when the particle center (with its undeformed spherical shape) goes beyond the wall surface, the distance vector $\mathit{d}$ changes direction with respect to the wall and so does the unit vector ${\mathit{e}}_n$ defining the force direction. In addition, the value of the overlap computed by Eq. 28 is smaller than the actual deformation, which should be bigger rather than smaller than the radius $R$ (see Figure 7B). Therefore, we propose to modify the two calculations as follows:

$\left.\delta =R+\text{sign}\left.\left(\mathit{d}\cdot \mathit{n}\right)\right|\mathit{d}\right|(30)$

${\mathit{e}}_n=-\text{sign}\left(\mathit{d}\cdot \mathit{n}\right)\frac{\mathit{d}}{\left|\mathit{d}\right|}(31)$

In which the sign function outputs +1 or −1 depending on the sign of the scalar product between the distance vector and the surface normal vector $\mathit{n}$ (see Figure 7B). Note that the surface normal is assumed to point toward the internal domain, otherwise the sign before the sign function must be the opposite in both formulas.

The above simple modification allows a far more robust consideration of the particle wall contact under extreme deformation cases. The particles can go with their center beyond the wall surface or even completely beyond it, without causing unexpected consequences. The increased overlap will determine a higher repulsive force, eventually safeguarding the presence of the particle inside the domain.

The robust calculation of the overlap discussed above shows no particular limit, e.g., a maximum deformation. The details in the implementation of DEM computational models can introduce limitations associated with the wall contact detection algorithm. Indeed, all codes implement some sort of distance cut-off to detect contacts, which means that if the particle is pressed against the wall so that its center is too far (inside) from the surface, then the contact with the wall will no longer be identified and considered active. This is expected to be a very extreme case, unlikely even for coarse particles pushing fine ones with significant particle size difference.

4.2 Critical velocity corresponding to higher maximum overlap

A simple way to investigate how the critical velocity changes in an expanded overlap range is to adapt the previous calculation by expressing it as a function of the maximum dimensionless overlap.

Using the linear model, the more general definition of the critical velocity becomes

$v_{0c,L}={\delta }_{\mathit{max}}^{*}\sqrt{\frac{3}{4}\frac{k}{\pi \rho R}}(32)$

yielding a linear dependence on the value of ${\delta }_{\mathit{max}}^{*}$. For the case of fine-on-coarse spheres in contact, the more general result is

$v_{0c12,L}={\delta }_{\mathit{max}}^{*}\sqrt{\frac{3}{4}\frac{k_1}{\pi {\rho }_1{R}_1\left(1+q^3\right)}}(33)$

With the Hertz theory approach, the two results are slightly different, i.e.,

$v_{0c,H}={{\delta }_{\mathit{max}}^{*}}^{\frac{5}{4}}\sqrt{\frac{4}{5}\frac{E}{\pi \rho }}(34)$

and

$v_{0c12,H}={{\delta }_{\mathit{max}}^{*}}^{\frac{5}{4}}\sqrt{\frac{4}{5}\frac{E_1}{\pi {\rho }_1\left(1+q^3\right)}}(35)$

yielding a superlinear dependence of the critical velocity on the maximum dimensionless overlap.

As one example, with a maximum dimensionless overlap ${\delta }_{\mathit{max}}^{*}=2$, i.e., assuming that the maximum overlap can extend up to one particle diameter, with the linear model, the critical velocity would turn out doubled with respect to the conventional case. With the Hertz model, the increase in maximum impact velocity before particle ejection would reach a value nearly 2.4 times higher than the original one.

The advantage of the proposed formulation is that the higher possible impact velocity is achieved without the need to use stiffer solids, which would cause the necessary integration time step to decrease and the total computational burden for the same simulation time to increase correspondingly.

4.3 Discussion of known limitations

At the end of Section 4.1 potential causes for overlap limits have been identified, mainly associated with the contact detection algorithm. Here, other limitations associated with geometrical considerations are briefly discussed.

The proposed solution is expected to work for truly thick walls. If thin walls are considered, the approach may introduce artificial and possibly unwanted thickness. It is the case of walls separating two regions of the system both populated with particles (Figure 8A). For example, the vortex finder in cyclones separates an internal region for the gas outflow and an annular region where the gas-solid mixture enters tangentially, so special care must be exercised.

FIGURE 8

FIGURE 8. Problematic cases for application of the thick wall concept when actual solid walls are thin. Forces shown as red arrows and overlap as blue segments. (A) Thin laminar wall separating two regions populated with particles, for which the wall thickness generates fictitious wall effects outside the actual thickness of the wall. (B) Sharp wedge-shaped wall, in which fictitious wall regions appear that exert action on the surrounding particles.

Other similar problematic cases include walls with wedges, pointed tips or other types of very sharp edges, where the thickness may act as corner smoothing or generate artifacts, possibly changing the particle behavior in the close proximity of these regions (see Figure 8B). The extent of the issue depends directly on the “depth” of the thick wall, i.e., on the maximum overlap allowed for particles going beyond the wall surfaces. In both cases shown in Figure 8, the forces shown, represented by red arrows, show that the particles are attracted to the closest wall surface (actually being repelled by the wall surface on the opposite side). So, because of the combined action of the two thick boundaries inside the real wall, the particles will likely end up trapped inside the wall at an equilibrium surface halfway between the two wall surfaces.

Another limitation is associated with the validity of the mechanical relations under such large overlaps. Hertz theory for the contact of an elastic sphere with a flat wall is valid in the limit of small deformations, a condition that evidently cannot be enforced here, where deformations are very significant. Unfortunately, it is far from easy to estimate the error introduced by extrapolating outside the original validity of the theory. On the other hand, considering the other assumptions commonly accepted in DEM simulations, such as perfect sphericity of the particles, highly softened elastic moduli, ideally smooth surfaces, homogeneous materials and so on, we speculate that in general the proposed thick wall concept remains within the same limits of acceptability, without further deteriorating the overall level of approximation.

5 DEM application to pharmaceutical particle motion in capsule

To demonstrate the usefulness of the thick wall concept, we show an application in transient modelling of the motion of carrier and active pharmaceutical ingredient (API) particles of pharmaceutical interest in an inhaler capsule during a simulated inhalation cycle. The computational model and simulation conditions details are available in (Alfano et al., 2022c; Alfano et al., 2023) and a brief summary is reported below.

5.1 Computational model and tool

As anticipated earlier, the Discrete Element Method (DEM) is used to simulate the particle-particle-structure system inside a rotating and vibrating capsule, representing the characteristic evolution during an inhalation process. Each particle is tracked by solving Newton’s equations of motion for the translational and rotational components:

$m_i\frac{{d\overrightarrow{v}}_i}{dt}={\displaystyle \sum _{j=1}^N}{\overrightarrow{F}}_{cont,ij}+{\overrightarrow{F}}_{g,i}+{\displaystyle \sum _{k=1}^{N_f}}{\overrightarrow{F}}_{fict,ik}(36)$

$I_i\frac{d{\overrightarrow{\omega }}_{pi}}{dt}={\displaystyle \sum _{j=1}^N}{\overrightarrow{T}}_{cont,ij}+{\overrightarrow{T}}_r(37)$

in which $m_i$, $v_i,I_i$ and ${\omega }_{pi}$ are the i-th particle mass, velocity, moment of inertia and angular velocity, respectively. In Eq. 36, the summation of external actions includes contact forces, ${\displaystyle \sum _{j=1}^N}{\overrightarrow{F}}_{cont,ij}$ and gravity, ${\overrightarrow{F}}_{g,i}$. In addition, to account for the capsule motion in the frame of reference of the capsule, fictitious forces arise including centrifugal, Euler and Coriolis forces, indicated by ${\displaystyle \sum _{k=1}^{N_f}}{\overrightarrow{F}}_{fict,i}$. In Eq. 37, the torques considered for the rotational motion are contact torque (${\displaystyle \sum _{j=1}^N}{\overrightarrow{T}}_{cont,ij}$) and rolling resistance (${\overrightarrow{T}}_r$), calculated according to the Constant Directional Torque model, CDT (Ai et al., 2011):

${\overrightarrow{T}}_r=-{\mu }_r{R}_{eq}\left|F_{cont,ij}^{\left(n\right)}\right|\frac{{\overrightarrow{\omega }}_{pi}-{\overrightarrow{\omega }}_{pj}}{\left|{\overrightarrow{\omega }}_{pi}-{\overrightarrow{\omega }}_{pj}\right|}(38)$

where ${\mu }_r$ is the rolling friction coefficient, $R_{eq}$ is the equivalent radius, $\left|F_{cont,ij}^{\left(n\right)}\right|$ is the magnitude of the normal contact force, ${\overrightarrow{\omega }}_{pi}$ and ${\overrightarrow{\omega }}_{pj}$ are the angular velocities of the two contacting objects $i$ and $j$.

In the contact forces, the Hertzian elastic component is complemented by an adhesion/cohesion model, named the Johnson-Kendall-Roberts (JKR) model (Johnson et al., 1971), including attractive forces for negative overlap. The model equations read:

$F_{el,JKR}=4\sqrt{\pi \gamma E_{eq}}{a}^{\frac{3}{2}}-\frac{4}{3}\frac{E_{eq}}{R_{eq}}{a}^3(39)$

${\delta }_n=\frac{a^2}{R_{eq}}-\sqrt{\frac{4\pi \gamma a}{E_{eq}}}(40)$

where $a$ is the radius of the contact area, $E_{eq}$ the equivalent Young modulus, ${\delta }_n$ the normal overlap, and γ the cohesion surface energy.

In addition to the elastic and cohesive forces, a viscous/damping force contribution is considered to account for a constant coefficient of restitution. The total normal contact force is the sum of the JKR contribution and the viscous damping term:

$\left|F_{cont,ij}^{\left(n\right)}\right|=F_{el}+{\eta }_n{\delta }_n^{\frac{1}{4}}{v}_n(41)$

where ${\eta }_n$ is the damping coefficient and $v_n$ is the normal component of the impact velocity. The tangential contact force is also taken into account according to a modified Mindlin-Deresiewicz no-slip solution, as reported by (Di Renzo and Di Maio, 2005).

As for the fictitious forces generated by the moving (i.e., non-inertial) frame of reference, two features have been implemented: 1), the gravity vector was made to rotate according to the angular velocity assumed for the capsule and 2) formulations for three force contributions, named the centrifugal force, Euler force (related to angular acceleration), and Coriolis force. The reader is referred to our previous work (Alfano et al., 2022c) for the details of the corresponding expressions.

The DEM implementation available in the MFiX code, version 18.1.5 (Garg et al., 2012), as modified by our group as detailed in Alfano et al. (2021b); Alfano et al. (2022c); Alfano et al. (2023), is adopted. The original open-source MFiX-DEM code has been extended to include the CDT rolling friction model, the JKR elastic-cohesive model including attractive forces also for negative overlaps upon detachment, all the fictitious forces discussed above and the implementation of the proposed thick wall concept, by means of Eqs. 30, 31.

5.2 Simulated capsule motion

Two DEM simulations were conducted, with and without the thick wall effect (wit critical overlap set to ${\delta }_{\mathit{max}}^{*}=2$), each lasting 10 ms, considering the capsule of a Dry Powder Inhaler (DPI). The capsule is depicted in Figure 9 with its sizes and its resting position in a DPI. It represents the same system studied in our previous works (Alfano et al., 2022c; Alfano et al., 2023). The capsule rotates with a constant angular acceleration of 887 rad/s², starting from rest. At the end of the simulation (10 ms), the capsule rotates with an angular velocity of approximately 5,000 rpm. Additionally, the capsule undergoes oscillatory motion parallel to the capsule axis, with a frequency of 30 Hz and an amplitude of 0.8 mm.

FIGURE 9

FIGURE 9. Side (A) and front (B) views of the capsule, along with its dimensions. The initial position inside the inhaler is shown in (C).

Within the capsule, 1 mg of a typical DPI formulation is present, composed by the carrier-API pair lactose-salbutamol. The powder configuration is illustrated in Figure 10. Salbutamol API particles adhere to a much larger carrier particle (lactose) to enhance the flowability properties of the formulation. The weight fraction of the active ingredient is approximately 1%. The active ingredient has a diameter of $5\mu \mathrm{m}$ (as required to reach the patient’s lower airways when inhaled), while the carrier has a diameter of $100\mu \mathrm{m}$. The mechanical and physical properties of the particles are detailed in Table 1. This simulation is particularly challenging due to the large size ratio between the two solid phases, i.e. 20 to 1. Additionally, the very small diameter of the API particles forces the use of a very small time-step in the DEM simulation, of the order of nanoseconds. The collision duration for fine particles is, in fact, 82 ns. The DEM time-step is set to 1/10 of the collision time.

FIGURE 10

FIGURE 10. API-carrier configuration. (A) single coated carrier particle; (B) all particles inside the capsule. Total mass: 1 mg; API is about 1% (w/w).

TABLE 1

TABLE 1. DEM simulation parameters.

Without the implementation of special treatments for the wall contacts, at various points during the simulation, some API particles end up ejected from the capsule walls. This is shown in a snapshot in Figure 11A, which provides a close-up of a moment when the phenomenon occurs for the first time (at about 0.1 ms). Fine violet particles are clearly visible below, outside the capsule. A deeper investigation of the local dynamics reveals that the responsibility is to be attributed to the mechanism shown in Figure 2B, for which a carrier particle either hits an API particle already touching the wall or it is the same carrier particle with sticking API fine particles that pushes some of them actually through the walls and out of the capsule. The process of particle ejection from the bottom part of the capsule is significant in the very first few instants. Then, it goes on at a smaller rate throughout the entire simulation time.

FIGURE 11

FIGURE 11. Close up of an instant in time (about 0.1 ms) when fine API particles (violet) get compressed between coarse carrier particles (grey) and the capsule walls. (A) Simulation without thick wall, (B) simulation with thick wall (the critical overlap is approximately 2 times the fine particle diameter).

In Figure 11B, the particle configuration is shown for the simulation with the thick wall, at the same moment in time and with the exact same angle. It can be observed that in this case, despite still having the identical simulation setup, particle configuration and integration time-step, the API particles do not exit the capsule geometry.

An analysis of the number of particles ejected outside the capsule over time is reported for both simulations in Figure 12. Despite its simplicity, it can be observed that the use of thick wall proves very effective in significantly restricting the phenomenon, with only 4 particles found outside the geometry, compared to more than 150 particles in the simulation without thick wall in such a short simulation time.

FIGURE 12

FIGURE 12. Number of ejected particles during the simulation in the absence or presence of the thick wall implementation.

It is useful to investigate the causes and compare with the critical velocity values for the contact model used (Hertz). The critical velocities $v_{0c,H}$ for the API and carrier particles are found to be 6.5 m/s, and 5.8 m/s, respectively. The resulting values are similar because both materials share the same Young’s modulus, and there is only a small difference in density. The API particle, being less dense, exhibits a slightly higher critical velocity.

In scenarios where the carrier compresses the API particle, with a size ratio of 20, a dramatic decrease is observed, as the critical velocity $v_{0c12,H}=$ 0.073 m/s is obtained. In this case, the velocity is 89 times lower than for the single API particle case.

Recalculations are done under thick wall conditions. For the API particle, the critical velocity increases to 15.5 m/s and for the carrier particle to 13.9 m/s. In the fine-coarse case, the critical velocity for API ejection becomes $v_{0c12,H}=0.17\mathrm{m}/\mathrm{s}$.

The simulation was conducted using the JKR model, which is derived from the Hertz model regarding the elastic contribution. In comparison with the linear model, values of the elastic constants can be estimated by imposing equality between the critical velocities obtained analytically (Eq. 22; Eq. 25), resulting in $k_1=500$ N/m, $k_2=10000$ N/m, respectively for fine and coarse particles. With these values, roughly the same percentage of particles expelled from the system is expected. Evaluations on critical velocity can thus be useful for estimating the parameters of the linear model.

Figure 13 shows the velocity magnitude distribution for carrier and API particles at the end of the simulation. It is observed that the majority of particles, both API and carrier, exhibit velocities lower than 1 m/s. However, the maximum velocity recorded for an API particle at this moment in time is 5.8 m/s (not shown), a value close to the critical velocity. The maximum velocity for a carrier particle is about 1.7 m/s. The data presented thus suggest that the escape of API particles is more likely due to compression between the carrier and the wall, although the possibility that some of the API particles are expelled due to having velocities greater than the critical velocity cannot be excluded. What can be concluded from the number of escaped particles is that the use of thick wall implementation almost completely eliminates the problem of particle ejection, by simply leading to an increase of the maximum critical velocity, without affecting the integration time-step and the corresponding required simulation time.

FIGURE 13

FIGURE 13. Velocity magnitude distribution for (A) API and (B) carrier particle at the end of the simulation (approx. 10 ms).

6 Conclusion

The problem of fine particles possibly subjected to ejection out of the domain in DEM simulation of highly polydisperse systems is addressed. Critical velocities leading to particle passing through walls are characterized analytically, showing the role played by the most relevant geometrical, physical and mechanical parameters, using both the linear model (with the spring stiffness) and the Hertz theory (with the Young’s modulus). Conditions are identified for which coarse particles can push the fine ones out of the domain as a result of excessive overlaps. This is not due to time-step limitations or approximated numerical algorithms, but solely to the difference in the colliding masses and ultimately to their size difference. A possible worsening factor is the—rather frequent—use of softer than real stiffnesses. A simple modification to the implementation of the contact overlap and force unit vector is proposed, conceptually named thick wall, according to which contacts with larger overlap are allowed, like up to the particle diameter, or even more, without sacrificing the integration time-step. Pharmaceutical powders for inhalation formulations are often composed by a coarse carrier coated by fine API particles. The shaking motion of a capsule for Dry Powder Inhalers is simulated by DEM. It is shown that the ejection issue is rather frequently encountered and fine particles continuously leave the capsule domain during shaking. By implementing the thick wall modification (maximum overlap twice the fine particle diameter), the problem nearly completely vanishes. It is noteworthy that the simulation time is totally unaffected by the change.

Limitations of the thick wall concept associated with the excessive overlap and cases in which thin or sharp wall shapes may lead to geometrical and physical artifacts are also discussed. Such limitations are inherent in the concept and cannot be easily overcome. Current basic DEM assumptions like calculations based on undeformed shapes naturally lead to sacrificing accuracy under all conditions. In the future, detailed consideration of deformable, non-spherical bodies interacting with deformable walls may lead to more physically based evaluation of the contact forces. However, the far larger additional cost of computing local deformations and stresses may be justified only with very high accuracy specifications. With the current computational power, this is only likely in systems with few, quasi-static particles. Similarly, additional elements of the contact could be introduced in future studies to increase the realism of the model, such as irregular particle shapes, dissipative phenomena occurring at high deformations and even particle breakage models.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

FA: Data curation, Formal Analysis, Investigation, Methodology, Software, Visualization, Writing–original draft, Writing–review and editing. GI: Data curation, Formal Analysis, Investigation, Methodology, Software, Writing–original draft, Writing–review and editing. FD: Conceptualization, Funding acquisition, Resources, Supervision, Writing–original draft, Writing–review and editing, Methodology, Software. AD: Conceptualization, Funding acquisition, Resources, Supervision, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work has benefitted from the support of the following projects, which are gratefully acknowledged: ICSC National Research Centre for High Performance Computing, Big Data and Quantum Computing, funded by European Union—Next Generation EU (grant number CN00000013, CUP: H23C22000360005), through Piano Nazionale di Ripresa e Resilienza, M4C2, Investimento 1.4 and PRIN PNRR “TRIBOTWIN” project (CUP: H53D23008560001) funded by European Union—Next Generation EU, component M4C2, investimento 1.1.

Conflict of interest

The 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 author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ai, J., Chen, J. F., Rotter, J. M., and Ooi, J. Y. (2011). Assessment of rolling resistance models in discrete element simulations. Powder Technol. 206, 269–282. doi:10.1016/j.powtec.2010.09.030

CrossRef Full Text | Google Scholar

Alfano, F. O., Benassi, A., Gaspari, R., Di Renzo, A., and Di Maio, F. P. (2021a). Full-scale DEM simulation of coupled fluid and dry-coated particle flow in swirl-based dry powder inhalers. Ind. Eng. Chem. Res. 60, 15310–15326. doi:10.1021/acs.iecr.1c02864

CrossRef Full Text | Google Scholar

Alfano, F. O., Di Maio, F. P., and Di Renzo, A. (2022a). Deagglomeration of selected high-load API-carrier particles in swirl-based dry powder inhalers. Powder Technol. 408, 117800. doi:10.1016/J.POWTEC.2022.117800

CrossRef Full Text | Google Scholar

Alfano, F. O., Di Renzo, A., and Di Maio, F. P. (2023). Discrete element method evaluation of triboelectric charging due to powder handling in the capsule of a DPI. Pharmaceutics 15, 1762. doi:10.3390/pharmaceutics15061762

PubMed Abstract | CrossRef Full Text | Google Scholar

Alfano, F. O., Di Renzo, A., Di Maio, F. P., and Ghadiri, M. (2021b). Computational analysis of triboelectrification due to aerodynamic powder dispersion. Powder Technol. 382, 491–504. doi:10.1016/J.POWTEC.2021.01.011

CrossRef Full Text | Google Scholar

Alfano, F. O., Di Renzo, A., Gaspari, R., Benassi, A., and Di Maio, F. P. (2022b). Modelling deaggregation due to normal carrier–wall collision in dry powder inhalers. Processes 10, 1661. doi:10.3390/pr10081661

CrossRef Full Text | Google Scholar

Alfano, F. O., Sommerfeld, M., Di Maio, F. P., and Di Renzo, A. (2022c). DEM analysis of powder deaggregation and discharge from the capsule of a carrier-based Dry Powder Inhaler. Adv. Powder Technol. 33, 103853. doi:10.1016/j.apt.2022.103853

CrossRef Full Text | Google Scholar

Antypov, D., and Elliott, J. A. (2011). On an analytical solution for the damped Hertzian spring. EPL Europhys. Lett. 94, 50004. doi:10.1209/0295-5075/94/50004

CrossRef Full Text | Google Scholar

Ariane, M., Sommerfeld, M., and Alexiadis, A. (2018). Wall collision and drug-carrier detachment in dry powder inhalers: using DEM to devise a sub-scale model for CFD calculations. Powder Technol. 334, 65–75. doi:10.1016/j.powtec.2018.04.051

CrossRef Full Text | Google Scholar

Chaurasiya, B., and Zhao, Y.-Y. (2020). Dry powder for pulmonary delivery: a comprehensive review. Pharmaceutics 13, 31. doi:10.3390/pharmaceutics13010031

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, Y., and Sommerfeld, M. (2015). Forces on micron-sized particles randomly distributed on the surface of larger particles and possibility of detachment. Int. J. Multiph. Flow 72, 39–52. doi:10.1016/J.IJMULTIPHASEFLOW.2015.01.006

CrossRef Full Text | Google Scholar

Di Maio, F. P., and Di Renzo, A. (2004). Analytical solution for the problem of frictional-elastic collisions of spherical particles using the linear model. Chem. Eng. Sci. 59, 3461–3475. doi:10.1016/j.ces.2004.05.014

CrossRef Full Text | Google Scholar

Di Renzo, A., Cello, F., and Di Maio, F. P. (2011). Simulation of the layer inversion phenomenon in binary liquid--fluidized beds by DEM–CFD with a drag law for polydisperse systems. Chem. Eng. Sci. 66, 2945–2958. doi:10.1016/J.CES.2011.03.035

CrossRef Full Text | Google Scholar

Di Renzo, A., and Di Maio, F. P. (2005). An improved integral non-linear model for the contact of particles in distinct element simulations. Chem. Eng. Sci. 60, 1303–1312. doi:10.1016/j.ces.2004.10.004

CrossRef Full Text | Google Scholar

Fulchini, F., Zafar, U., Hare, C., Ghadiri, M., Tantawy, H., Ahmadian, H., et al. (2017). Relationship between surface area coverage of flow-aids and flowability of cohesive particles. Powder Technol. 322, 417–427. doi:10.1016/j.powtec.2017.09.013

CrossRef Full Text | Google Scholar

Garg, R., Li, T., and Pannala, S. (2012). Documentation of open-source MFIX–DEM software for gas-solids flows. Available at: https://mfix.netl.doe.gov/doc/mfix-archive/mfix_current_documentation/dem_doc_2012-1.pdf.

Google Scholar

Golshan, S., Sotudeh-Gharebagh, R., Zarghami, R., Mostoufi, N., Blais, B., and Kuipers, J. A. M. (2020). Review and implementation of CFD-DEM applied to chemical process systems. Chem. Eng. Sci. 221, 115646. doi:10.1016/J.CES.2020.115646

CrossRef Full Text | Google Scholar

Grima, A. P., and Wypych, P. W. (2011). Investigation into calibration of discrete element model parameters for scale-up and validation of particle-structure interactions under impact conditions. Powder Technol. 212, 198–209. doi:10.1016/j.powtec.2011.05.017

CrossRef Full Text | Google Scholar

Hærvig, J., Kleinhans, U., Wieland, C., Spliethoff, H., Jensen, A. L., Sørensen, K., et al. (2017). On the adhesive JKR contact and rolling models for reduced particle stiffness discrete element simulations. Powder Technol. 319, 472–482. doi:10.1016/J.POWTEC.2017.07.006

CrossRef Full Text | Google Scholar

He, Y., Hassanpour, A., Alizadeh Behjani, M., and Bayly, A. E. (2021). A novel stiffness scaling methodology for discrete element modelling of cohesive fine powders. Appl. Math. Model 90, 817–844. doi:10.1016/J.APM.2020.08.062

CrossRef Full Text | Google Scholar

Huang, C., and Nakagawa, M. (2023). Effects of rotation axis on mixing behavior of dissimilar particles in rotating drums. Powder Technol. 428, 118868. doi:10.1016/J.POWTEC.2023.118868

CrossRef Full Text | Google Scholar

Johnson, K. L. (1987). Contact mechanics. Cambridge: Cambridge University Press.

Google Scholar

Johnson, K. L., Kendall, K., and Roberts, A. D. (1971). Surface energy and the contact of elastic solids. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 324, 301–313.

Google Scholar

Karde, V., Khala, M. J., Hare, C., and Heng, J. Y. Y. (2023). Use of shear sensitive coloured guest component to track powder mixing in adhesive binary mixtures. Powder Technol. 414, 118068. doi:10.1016/J.POWTEC.2022.118068

CrossRef Full Text | Google Scholar

Khala, M. J., Hare, C., Karde, V., and Heng, J. Y. Y. (2023). A numerical analysis of the influence of material properties on dry powder coating performance. Chem. Eng. Res. Des. 193, 158–167. doi:10.1016/J.CHERD.2023.03.028

CrossRef Full Text | Google Scholar

Kieckhefen, P., Pietsch, S., Dosta, M., and Heinrich, S. (2020). Possibilities and limits of computational fluid dynamics–discrete element method simulations in process engineering: a review of recent advancements and future trends. Annu. Rev. Chem. Biomol. Eng. 11, 397–422. doi:10.1146/annurev-chembioeng-110519-075414

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruggel-Emden, H., Stepanek, F., and Munjiza, A. (2010). A study on adjusted contact force laws for accelerated large scale discrete element simulations. Particuology 8, 161–175. doi:10.1016/J.PARTIC.2009.07.006

CrossRef Full Text | Google Scholar

Kuo, H. P., Knight, P. C., Parker, D. J., Tsuji, Y., Adams, M. J., and Seville, J. P. K. (2002). The influence of DEM simulation parameters on the particle behaviour in a V-mixer. Chem. Eng. Sci. 57, 3621–3638. doi:10.1016/S0009-2509(02)00086-6

CrossRef Full Text | Google Scholar

Lischka, C., Gerl, S., Kappes, J., Chauhan, A., and Nirschl, H. (2024). Experimental & simulative assessment of mixing quality for dry Li-Ion cathode production in an Eirich intensive mixer. Powder Technol. 431, 119072. doi:10.1016/J.POWTEC.2023.119072

CrossRef Full Text | Google Scholar

Maw, N., Barber, J. R., and Fawcett, J. N. (1976). The oblique impact of elastic spheres. Wear 38, 101–114. doi:10.1016/0043-1648(76)90201-5

CrossRef Full Text | Google Scholar

Peng, Z., Joshi, J. B., Moghtaderi, B., Khan, Md. S., Evans, G. M., and Doroodchi, E. (2016). Segregation and dispersion of binary solids in liquid fluidised beds: a CFD-DEM study. Chem. Eng. Sci. 152, 65–83. doi:10.1016/j.ces.2016.05.032

CrossRef Full Text | Google Scholar

Schmidt, J., and Peukert, W. (2022). Dry powder coating in additive manufacturing. Front. Chem. Eng. 4, 995221. doi:10.3389/fceng.2022.995221

CrossRef Full Text | Google Scholar

Schulz, D., Reinecke, S. R., Woschny, N., Schmidt, E., and Kruggel-Emden, H. (2023). Development and evaluation of force balance based functions for dust detachment from bulk particles stressed by fluid flow. Powder Technol. 417, 118257. doi:10.1016/J.POWTEC.2023.118257

CrossRef Full Text | Google Scholar

Schulz, D., Woschny, N., Schmidt, E., and Kruggel-Emden, H. (2022). Modelling of the detachment of adhesive dust particles during bulk solid particle impact to enhance dust detachment functions. Powder Technol. 400, 117238. doi:10.1016/J.POWTEC.2022.117238

CrossRef Full Text | Google Scholar

Shaheen, M. Y., Thornton, A. R., Luding, S., and Weinhart, T. (2021). The influence of material and process parameters on powder spreading in additive manufacturing. Powder Technol. 383, 564–583. doi:10.1016/J.POWTEC.2021.01.058

CrossRef Full Text | Google Scholar

Sharma, R., and Setia, G. (2019). Mechanical dry particle coating on cohesive pharmaceutical powders for improving flowability - a review. Powder Technol. 356, 458–479. doi:10.1016/J.POWTEC.2019.08.009

CrossRef Full Text | Google Scholar

Sofia, D., Chirone, R., Lettieri, P., Barletta, D., and Poletto, M. (2018). Selective laser sintering of ceramic powders with bimodal particle size distribution. Chem. Eng. Res. Des. 136, 536–547. doi:10.1016/j.cherd.2018.06.008

CrossRef Full Text | Google Scholar

Spahn, J. E., Zhang, F., and Smyth, H. D. C. (2022). Mixing of dry powders for inhalation: a review. Int. J. Pharm. 619, 121736. doi:10.1016/J.IJPHARM.2022.121736

PubMed Abstract | CrossRef Full Text | Google Scholar

Thornton, C., Cummins, S. J., and Cleary, P. W. (2011). An investigation of the comparative behaviour of alternative contact force models during elastic collisions. Powder Technol. 210, 189–197. doi:10.1016/j.powtec.2011.01.013

CrossRef Full Text | Google Scholar

Tong, Z. B., Yang, R. Y., and Yu, A. B. (2017). CFD-DEM study of the aerosolisation mechanism of carrier-based formulations with high drug loadings. Powder Technol. 314, 620–626. doi:10.1016/j.powtec.2016.10.004

CrossRef Full Text | Google Scholar

Washino, K., Chan, E. L., and Tanaka, T. (2018). DEM with attraction forces using reduced particle stiffness. Powder Technol. 325, 202–208. doi:10.1016/j.powtec.2017.11.024

CrossRef Full Text | Google Scholar

Washino, K., Nakae, S., Yamagami, R., Chan, E. L., Tsuji, T., and Tanaka, T. (2024). Scaling of attraction force and rolling resistance in DEM with reduced particle stiffness. Chem. Eng. Res. Des. 203, 501–519. doi:10.1016/J.CHERD.2024.02.006

CrossRef Full Text | Google Scholar

Yang, J., Wu, C.-Y., and Adams, M. (2015). DEM analysis of the effect of electrostatic interaction on particle mixing for carrier-based dry powder inhaler formulations. Particuology 23, 25–30. doi:10.1016/j.partic.2014.12.007

CrossRef Full Text | Google Scholar

Zhu, H. P., and Yu, a. B. (2006). A theoretical analysis of the force models in discrete element method. Powder Technol. 161, 122–129. doi:10.1016/j.powtec.2005.09.006

CrossRef Full Text | Google Scholar

Zhu, H. P., Zhou, Z. Y., Yang, R. Y., and Yu, A. B. (2007). Discrete particle simulation of particulate systems: theoretical developments. Chem. Eng. Sci. 62, 3378–3396. doi:10.1016/J.CES.2006.12.089

CrossRef Full Text | Google Scholar

Nomenclature