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 ( δ 0 = 0 ), the following linear ODE governs the dynamics for all the contact duration: d 2 δ d t 2 + ω 2 δ = 0 (1) in which ω = k / m is the oscillation pulsation, and m = ρ 4 3 π R 3 is the particle mass and ρ its density. As discussed in more detail in Di Maio and Di Renzo (2004), the solution is δ t = v 0 ω sin   ω t (2) and the collision duration is τ = π / ω . Since the contact is a perfect sinusoidal oscillation, the maximum overlap occurs at half of the oscillation duration, i.e., δ max = v 0 ω (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: δ * = δ R ;   t * = t τ (4)

The ODE can be re-expressed in the following terms: d 2 δ * d t * 2 + π 2 δ * = 0 (5) and the initial conditions are δ 0 * = 0 and d δ * d t * = π v 0 R ω .

A characteristic dimensionless number can be defined, which is named here “impact number” I n : I n = v 0 R ω (6)

The dimensionless solution of Eq. 5 is δ * t * = I n   sin   π t * (7) and δ max * = I n .

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 = 4 3 E R δ 3 2 (12)

The corresponding ODE governing the transient motion becomes non-linear d 2 δ d t 2 + E ρ π R 5 / 2 δ 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 δ max = 15 16 m R E v 0 2 2 5 (14)

By expressing the particle mass as function of the density and size, the maximum dimensionless overlap is easily obtained as δ max * = δ max R = 5 4 π ρ E v 0 2 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 0 c , L = E / π ρ (Eq. 11) appears.

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 ω 12 = k 1 / m 1 + m 2 , which for similar density solids leads to ω 12 = k 1 / 4 3 π ρ 1 R 1 3 + R 2 3  

The collision duration is τ 12 = π / ω 12 . The maximum allowed overlap is R 1 . So, the dimensionless variable will be considered as follows δ * = δ R 1 ;   t * = t τ 12 (17)

The dimensionless impact number is I n 12 = v 0   12 R 1 ω 12 (18) and, as before, the solution for the highest value of the dimensionless displacement is δ 12 , ⁡ max * = I n 12 .

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 δ 12 , ⁡ max = 15 16 m 1 + m 2 R 1 E 1 v 0 2 2 5 (23)

The dimensionless form of such overlap can be expresses as δ 12 , ⁡ max * = δ 12 , ⁡ max R 1 = 5 4 π ρ 1 E 1 1 + q 3 v 0 2 2 5 (24)

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).

The wall normal force is expressed by F = − k δ e n (27)

The (scalar) overlap and the unit vector for the direction are computed, respectively, by δ = R − d (28) e n = d d (29) in which 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 d changes direction with respect to the wall and so does the unit vector 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: δ = R + sign d ⋅ n d (30) e n = − sign d ⋅ n d d (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 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.

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 0 c , L = δ max * 3 4 k π ρ R (32) yielding a linear dependence on the value of δ max * . For the case of fine-on-coarse spheres in contact, the more general result is v 0 c 12 , L = δ max * 3 4 k 1 π ρ 1 R 1 1 + q 3 (33)

With the Hertz theory approach, the two results are slightly different, i.e., v 0 c , H = δ max * 5 4 4 5 E π ρ (34) and v 0 c 12 , H = δ max * 5 4 4 5 E 1 π ρ 1 1 + q 3 (35) yielding a superlinear dependence of the critical velocity on the maximum dimensionless overlap.

As one example, with a maximum dimensionless overlap δ 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.

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.

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.

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 d v → i d t = ∑ j = 1 N F → c o n t , i j + F → g , i + ∑ k = 1 N f F → f i c t , i k (36) I i d ω → p i d t = ∑ j = 1 N T → c o n t , i j + T → r (37) in which m i , v i , I i and ω p i 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, ∑ j = 1 N F → c o n t , i j and gravity, 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 ∑ k = 1 N f F → f i c t , i . In Eq. 37, the torques considered for the rotational motion are contact torque ( ∑ j = 1 N T → c o n t , i j ) and rolling resistance ( T → r ), calculated according to the Constant Directional Torque model, CDT (Ai et al., 2011): T → r = − μ r R e q F c o n t , i j n ω → p i − ω → p j ω → p i − ω → p j (38) where μ r is the rolling friction coefficient, R e q is the equivalent radius, F c o n t , i j n is the magnitude of the normal contact force, ω → p i and ω → p j 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 e l , J K R = 4 π γ E e q a 3 2 − 4 3 E e q R e q a 3 (39) δ n = a 2 R e q − 4 π γ a E e q (40) where a is the radius of the contact area, E e q the equivalent Young modulus, δ 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: F c o n t , i j n = F e l + η n δ n 1 4 v n (41) where η 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.

Two DEM simulations were conducted, with and without the thick wall effect (wit critical overlap set to δ 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.

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   μ m (as required to reach the patient’s lower airways when inhaled), while the carrier has a diameter of 100   μ 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.

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.

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.

It is useful to investigate the causes and compare with the critical velocity values for the contact model used (Hertz). The critical velocities v 0 c , 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 0 c 12 , 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 0 c 12 , H = 0.17   m / 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.