An integrated DEM–FEM approach to study breakage in packing of glass cartridges on a conveyor belt

Abstract

The use of glass for pharmaceutical new applications such as high-technology drugs, requires the strictest container inertness. A common theme of paramount importance in glass container integrity preservation is the detailed mechanism driving the sudden failure due the crack propagation. Using a combination of discrete element method (DEM) and finite element method (FEM), a stress map for glass cartridges packed into an accumulation table and transported by a conveyor belt at a fixed velocity is obtained under realistic conditions. The DEM calculation provides a full description of the dynamics of the cartridges, as approximated by an equivalent sphere, as well as the statistics of the multiple collisions. The FEM calculation exploits this input to provide the maximum principal stress of different pairs as a function of time. Our analysis shows that, during their transportation on the conveyor belt, the cartridges are subject to several shocks of varying intensities. Under these conditions, a crack may originate inside the cartridge in the area of maximal tensile stress, and propagate outward. Estimated stresses are found in good agreement with real systems.

Graphic abstract

1 Introduction

Frictional contacts between macroscopic objects may lead to deformation and breakage. The analysis of conditions for breakage in a large packing of objects with multiple contacts and collisions requires the evaluation of phenomena at different scales, from large scale force propagation (with strong fluctuations in space and time) down to local microscale stress accumulation. Here we analyse a case study, inspired by an industrial problem, which reveals several unexpected facets and illustrates the connection between material mechanics, physics and granular statistical physics.

The use of glass for pharmaceutical keeps increasing, especially for new applications such as high-technology drugs, requiring the strictest container inertness. Indeed, high-end biotech products keep entering the market, posing new technical and regulatory challenges for the development of injection systems [1]. Among the primary containers, syringes and cartridges are becoming the key players thanks to their superior performances and easy usability, replacing the reconstitution filling procedures for vials. Following the market trends, new manufacturing processes and advanced technologies have been developed, making available a new range of primary packaging aimed at increasing patient safety and reducing the risk of market recall. Regarding this, one of the most critical topics remains the mechanical resistance: glass breakages can lead to serious issues like container integrity failure, microbiological contamination of a sterile package and glass particulates generation. Almost always, the decrease of the mechanical resistance [2] is caused by strength-reducing damage introduced during forming and handling processes. For this reason, studies for risk assessment, based on physical and functional factors, become fundamental, starting from the early stages of the primary container manufacturing. The above problem is not confined to technological issues, but it also bears important consequences in terms of slow dynamics non-linearity in brittle materials [3], and in the problem of slow flows of granular materials drained by conveyor belts [4, 5], not dissimilar from the widespread problem of granular discharge [6]. From this perspective, a numerical modelling that simulates the critical steps of the container life-cycle can be a powerful tool to reduce the risk of those adverse events. Figure 1 shows the typical conditions that our calculation would aim at reproducing.

In this paper, we will discuss the glass breakage issue for cylindrical cartridges, that are one of the most commonly used containers in pharmaceutical processes.

Fig. 1

a A typical production line with the accumulation table for the cartridge; b three cartridges with characteristic cracks in different parts of the cartridge: body area (left), shoulder area (center), fire polish area (right)

Cartridges are directed through an accumulation table toward the filling line (see Fig. 1a). During this process, they experience a high packing condition that may lead to the onset of invisible micro-fractures (examples are given in Fig. 1b) that affect the glass integrity and their safe use in the filling process. The idea of the present study is to identify the critical conditions for this to occur. This will be attained by a three-step process.

• 1. Discrete element method (DEM) and finite element method (FEM) to justify and to define the equivalent sphere model together with its parameters (input data for step 2).

  Firstly, the optimal geometrical and dynamical conditions will be studied using a discrete element method (DEM) calculation [7] as provided by the Yade package [8]. To this aim, we have built a model for a cylindrical cartridge as shown in Figure S1 of Supplementary Information. Each cylindrical cartridge is composed by 632 identical small beads—a smaller number of small spheres would imply too a poor resolution and make the shape of the simulated cartridges unfaithful to the real ones. A full analysis of the dynamics of many such cartridges proves however a daunting task, especially under realistic conditions, because it would involve the bookkeeping of each of the spheres composing each single cartridge, with a corresponding huge computational cost. Therefore we provide evidence that under high packing conditions the most relevant features of the dynamics can be still captured by using a set of equivalent spheres, that were obtained by the same glass material and total mass of the cartridges. A visual representation of this is displayed in Fig. 2a. As cartridges in the filling line are highly packed, their dynamics is essentially driven by by translations in the plane perpendicular to gravity and rotations (for instance due to friction) around the axis of symmetry of the cylindrical cartridge, i.e. other movements such as inclinations/falling of a cartridge are negligible. Hence, the use of an equivalent sphere turns out to be fully justified at this stage, as it will be further discussed below. A large number of these equivalent spheres is then inserted into a computational box mimicking the accumulation table displayed in Fig. 1a, with an underline conveyor belt whose velocity can be tuned. In this way, different geometries matching those of the accumulation table can also be tested, and the final optimal configuration quickly identified. See Sect. 3.1

• 2. DEM simulations to obtain input data for step 3;

  In the second step, we will perform a direct DEM calculation of the velocity field of the cartridge population (as represented by their equivalent spheres), as well as the number of contact points for selected items, as a function of time, during the process of accumulation toward the filling line. The velocity field and other dynamical information derived from the DEM analysis will be then used as an input to perform a finite element method (FEM) analysis [9] on a small cluster of cartridges, using the Abaqus package [10]. See Sect. 3.4

• 3. FEM simulations to assess the likelihood of crack initiation;

  This final step will provide a complete stress map of the cartridges, which allows the identification of the critical regions where microfractures can occur under that specific geometry. See Sect. 4.1.

The plan of the paper is as follows. In Sect. 2 we discuss the mapping of a cartridge into an equivalent sphere, and the theoretical background underlying the DEM and FEM analyses. In Sect. 3, we report all the obtained results, first from the DEM to identify the optimal geometry of the filling line, as well as the velocity field, and then from the FEM to study the stress map of a small cluster of cartridges subject to collisions as predicted by the DEM calculation (Sect. 4). Finally, Sect. 5 summarizes some key messages and presents some possible future perspectives.

Fig. 2

a Effective spheres representing the cartridge in DEM simulations; b the collision geometry of two cartridges

2 Theoretical background and methods

2.1 The Hertz contact model for full spheres, Reissner theory of thin shells and the cartridge behavior

Consider a cartridge travelling on a filling line and impacting a neighboring one as schematically depicted in Fig. 2b. The total deformation can be split into two main contributions: the mutual indentation depth and the deflection of the thin glass walls due to membrane and bending actions. If small displacements and strains along with linear elasticity hold, then the two terms can be calculated separately, and the system can be modeled by combining the corresponding stiffness in series. Furthermore, since the contact area is limited to a small zone, the two contacting bodies can be assimilated to two equivalent empty spheres, keeping the same mass of the cartridges. Under these assumptions, the maximum mutual indentation can be calculated following Hertz’ theory of non-adhesive elastic contact between two spheres [9]

$$\begin{aligned} F=k_Hd^{\frac{3}{2}}_H \end{aligned}$$

(1)

where $k_H$ is the Hertz stiffness

$$\begin{aligned} k_H=\frac{4}{3}\sqrt{R^*}E^* \end{aligned}$$

(2)

and:

$$\begin{aligned} \frac{1}{R^*}= & {} \frac{1}{R_1}+\frac{1}{R_2} \end{aligned}$$

(3)

$$\begin{aligned} \frac{1}{E^*}= & {} \frac{(1-\nu ^2_1)}{E_1}+\frac{(1-\nu ^2_2)}{E_2} \end{aligned}$$

(4)

The membrane and bending deformation is evaluated after Reissner theory [11, 12], and the maximum deflection of one shell is given by the linear relation

$$\begin{aligned} F=k_Rd_R \end{aligned}$$

(5)

where:

$$\begin{aligned} k_R=\frac{4}{\sqrt{3(1-\nu ^2)}}\frac{Et^2}{R} \end{aligned}$$

(6)

For a two-cartridge system accounting for both Hertz indentation and Reissner theory, the energy balance equation can be written as

$$\begin{aligned} \frac{1}{2}m_1v_1^2+\frac{1}{2}m_2v_2^2=\frac{1}{2}k_{R1}d_{R1}^2 +\frac{1}{2}k_{R2}d_{R2}^2+\frac{2}{5}k_Hd_H^{\frac{5}{2}} +\frac{1}{2}(m_1+m_2)v_{12}^2 \end{aligned}$$

(7)

The relative velocity ${\mathbf {v}}_{12}$ is obtained by imposing the linear momentum balance

$$\begin{aligned} {\mathbf {v}}_{12}=\frac{m_1{\mathbf {v}}_1+m_2{\mathbf {v}}_2}{m_1+m_2} \end{aligned}$$

(8)

Analogously to $E^*$ and $R^*$, we can define the reduced mass as:

$$\begin{aligned} \frac{1}{m^*}=\frac{1}{m_1}+\frac{1}{m_2} \end{aligned}$$

(9)

Equations (8) and (9) substituted in Eq. (7) lead to:

$$\begin{aligned} m^*v^2_{rel}=k_{R1}d_{R1}^2+k_{R2}d_{R2}^2+\frac{4}{5}k_Hd_H^\frac{5}{2} \end{aligned}$$

(10)

where ${\mathbf {v}}_{rel}$ is the relative velocity of the two objects:

$$\begin{aligned} |{\mathbf {v}}_{rel}|=|{\mathbf {v}}_{1}-{\mathbf {v}}_{2}| \end{aligned}$$

(11)

Considering (1) and (5), Eq. (10) can be written as:

$$\begin{aligned} m^*v_{rel}^2=\frac{F^2}{k_{R1}}+\frac{F^2}{k_{R2}} +\frac{4}{5}\frac{F^{\frac{5}{3}}}{k_{H}^{\frac{2}{3}}} \end{aligned}$$

(12)

which is an implicit expression for the maximum force transmitted, including the effect of both Hertzian and membrane-bending deformation for the two empty spheres. Equation (12) can be easily solved numerically for the maximum force F. As an alternative, assuming that Hertz contact can be expresses by a linear law [13]

$$\begin{aligned} k_{linear}=\frac{F}{d_H}=\frac{4}{3}\left( \frac{15}{16}\right) ^{\frac{1}{5}} R^{*\frac{2}{5}}E^{*\frac{4}{5}}m^{*\frac{1}{5}}v_{rel}^\frac{2}{5} \end{aligned}$$

(13)

an explicit expression for the contact force is obtained as:

$$\begin{aligned} F=\sqrt{k_{eq}m^*}|{\mathbf {v}}_{rel}| \end{aligned}$$

(14)

where the equivalent stiffness is simply given by:

$$\begin{aligned} \frac{1}{k_{eq}}=\frac{1}{k_{R1}}+\frac{1}{k_{R2}}+\frac{1}{k_{linear}} \end{aligned}$$

(15)

Figure 3 reports the contact force versus the Young modulus as computed using Hertz only theory, Reissner only theory, the approximate relation (15) and a finite element analysis of the two cartridges. See Sect. 3.1 for a discussion of the results.

Fig. 3

Maximum contact force as computed using the finite element method, Reissner only theory, Hertz only theory and a combination of Reissner and Hertz models

2.2 The discrete element method (DEM)

The discrete element method (DEM) [7] is a time-stepping algorithm that hinges upon the following calculation cycle:

• 1. Contacting forces are generated on the basis of the interactions between spheres

• 2. Spheres are displaced based on the forces previously computed and on the current positions

The second step is a straightforward application of Newtons’ Second Law that for a rigid body has a part related to forces

$$\begin{aligned} m_{i} \frac{d {\mathbf {v}}_i}{dt}= & {} \sum _{j=1}^N \left[ {\mathbf {F}}_{ij}^{(c)}+ {\mathbf {F}}_{ij}^{(nc)} \right] +{\mathbf {G}}{i} \end{aligned}$$

(16)

where ${\mathbf {F}}_{ij}^{(c)}$ and ${\mathbf {F}}_{ij}^{(nc)}$ are the contact and non-contact forces acting on the ith sphere due to all others, and ${\mathbf {G}}_{i}$ is the gravity, and a part related to torques

$$\begin{aligned} {\mathbf {I}}_{i} \frac{d \varvec{\omega }_{i}}{dt}= & {} {\mathbf {M}}_{i} \end{aligned}$$

(17)

where ${\mathbf {I}}_{i}$ is the moment of inertia tensor of the ith sphere, and $\varvec{\omega }_{i}$ and ${\mathbf {M}}_{i}$ are the angular velocity of the ith sphere and the total torque acting on the ith sphere, respectively. The Newton’s equations are integrated through the customary leapfrog scheme (also known as the Verlet scheme) [14].

The first step, on the other hand, is a more delicate one because depends upon the model at hand. It hinges upon the force-displacement principle describing the relationship between the relative movement of the two spheres and the contacting forces. The latter can always be decomposed in two parts: the normal direction that is always pointing from the center of one sphere to the other, and strongly depends on the contact (overlap) between the two bodies; and the tangential direction, whose relative movement is composed of a rotation and a translation components.

Once that both steps have been performed, time is advanced of one step, and the cycle is repeated again, until a prescribed number of time steps is reached.

Particular care has to be given to the boundary conditions, as the system is usually constrained by walls that may be constituted by different materials, and velocities cannot be assigned to walls.

All DEM simulations were carried out using an in-house adaptation of the Yade package [8].

2.3 The finite element method (FEM)

Many physical problems can be posed by a set of (partial) differential equations and suitable boundary conditions. The finite element method (FEM) [15] is a numerical approach to obtain an approximate solution to such boundary value problems.

In matrix notation, let ${\mathbf {A}}({\mathbf {u}})={\mathbf {0}}$ be a set of differential equations defined over a domain $\varOmega$ with boundary conditions ${\mathbf {B}}({\mathbf {u}})={\mathbf {0}}$ on the boundary $\varGamma$ of the domain. In the finite element approach, the set of governing differential equations needs to be transformed into an integral form:

$$\begin{aligned} \int _{\varOmega } {\mathbf {v}}^{T} {\mathbf {A}}\left( {\mathbf {u}}\right) \,d\varOmega +\int _{\varGamma } \overline{{\mathbf {v}}}^{T} {\mathbf {B}} \left( {\mathbf {u}}\right) \,d\varGamma =0 \end{aligned}$$

(18)

where ${\mathbf {v}}$ and $\overline{\mathbf {v}}$ are two sets of arbitrary test functions. FEM cuts the spatial domain into a number of simple geometric sub-domains $\varOmega _{e}$—the finite elements—which are assumed to be interconnected only at a discrete number of common points—the nodes—on their boundaries $\varGamma _{e}$. The values ${\mathbf {a}}_i$ of the unknown functions ${\mathbf {u}}$ at the nodal points are the basic unknowns of the problem. The complete approximate solution is obtained by interpolating the nodal values ${\mathbf {a}}_i$ by means of interpolation functions ${\mathbf {N}}_i$ (the so-called shape functions). Therefore, the finite element method seeks the solution in the following approximate form

$$\begin{aligned} {\mathbf {u}} \approx \overline{{\mathbf {u}}} = \sum _{i} {\mathbf {N}}_i {\mathbf {a}}_i = {\mathbf {N}}{\mathbf {a}} \end{aligned}$$

(19)

Obviously, the approximate solution (19) cannot satisfy both the differential equations and the boundary conditions in a general case.

As a consequence, in Eq. (18) instead of the arbitrary functions ${\mathbf {v}}$ and $\overline{\mathbf {v}}$, a set of approximate functions is considered

$$\begin{aligned} {\mathbf {v}}= & {} \sum _{j=1}^n {\mathbf {w}}_j \delta {\mathbf {a}}_j \end{aligned}$$

(20)

$$\begin{aligned} \overline{\mathbf {v}}= & {} \sum _{j=1}^n \overline{\mathbf {w}}_j \delta {\mathbf {a}}_j \end{aligned}$$

(21)

where $\delta {\mathbf {a}}_j$ are arbitrary and n is the number of unknowns of the problem. In this way, the integral form (18) leads to a set of ordinary algebraic equations

$$\begin{aligned} \delta {\mathbf {a}}_{j}^{T}\left[ \int _{\varOmega } {\mathbf {w}}_j^{T} {\mathbf {A}}\left( {\mathbf {N}} {\mathbf {a}}\right) \,d\varOmega +\int _{\varGamma } \overline{\mathbf {w}}_j^{T} {\mathbf {B}} \left( {\mathbf {N}} {\mathbf {a}}\right) \,d\varGamma \right] =0 \qquad j=1,\ldots ,n \end{aligned}$$

(22)

from which the unknowns ${\mathbf {a}}_i$ can be determined.

Since ${\mathbf {A}}\left( {\mathbf {N}} {\mathbf {a}}\right)$ and ${\mathbf {B}}\left( {\mathbf {N}} {\mathbf {a}}\right)$ represent the errors due to the substitution of the approximate solution (19) into the differential equations and boundary conditions, Eq. (22) is a weighted integral of the residuals. According to the choice of the set of functions ${\mathbf {w}}_j$ and $\overline{\mathbf {w}}_j$, different weighted residuals methods are obtained. In the finite element approach, the shape functions are used as weighting functions (Bubnov–Galerking method)

$$\begin{aligned} {\mathbf {w}}_j = {\mathbf {N}}_j \end{aligned}$$

(23)

Finally, by virtue of the property of definite integrals, the approximate solution can be obtained after each element contribution is assembled:

$$\begin{aligned} \int _{\varOmega } \left( \bullet \right) \,d\varOmega + \int _{\varGamma } \left( \bullet \right) \,d\varGamma = \sum _{e=1}^m \left( \int _{\varOmega _e} \left( \bullet \right) \,d\varOmega _e + \int _{\varGamma _e} \left( \bullet \right) \,d\varGamma _e \right) \end{aligned}$$

(24)

where $(\bullet )$ is a short-hand notation for an arbitrary function, and where $\varOmega _e$ is the domain of each element, $\varGamma _e$ its part of boundary and m is the number of finite elements in the operated discretization.

3 Results

3.1 Justification of the use of equivalent spheres

As anticipated, one major difficulty of the DEM calculation is the very high computational cost required by the use of the full cartridge model shown in Fig. S1 of SI. We set up a computational box with geometries matching that of the real accumulation table for the cartridge, and inserted a large number of cartridges in a high packing conditions as displayed in Fig. S2(a) (details on the computational protocol will be given below). A full dynamic evolution of the system (see Fig. S2(b) for its evolution after few seconds) shows that cartridges are driven by the underlying belt by maintaining their upright position, in view of the high packing condition that they experience. Hence, all rotational degrees of freedom associated with the the tilt and possible fall for a cartridges are essentially ruled out under such conditions. The cartridge can then be replaced by a simpler object, such as a sphere, to infer the relative kinematics (i.e. velocities) and dynamics (i.e. forces) to first approximation. The accuracy of this approximation will be checked a posteriori by a comparison of the full evolution of a system of cartridges in one specific case.

In order to implement a mapping of a cartridge into an equivalent sphere, we need to ensure that the relative contact forces during the collisions of spheres as a function of the elapsed time will be in quantitative agreement with those originating from the cartridges. While simple in principle, care must be used as several possibilities arise. As presented in Sect. 2, a cartridge can be either considered as filled deformable sphere with a fully elastic behavior (Hertz theory) [16,17,18], or as thin shallow spherical shell (Reissner model) [11, 12, 19], or as combination of the two (combined Hertz–Reissner model) [20,21,22,23]. In order to identify the optimal model, a numerical analysis of the collision of two spheres, within the three different models described above, is compared with the collision of two cartridges obtained from a direct finite element (FE) numerical method. In all cases, the spheres were made of the same material as the cartridges and equivalent dimensions so to keep the same mass (see Table S1 in supplementary information) and the two spheres or cartridges were designed to collide at a relative velocity of $20\,{\text {cm}}\,{\text {s}}^{-1}$ as pictorially illustrated in Fig. 2b. In the cartridge case, one of the two cartridges stands in the upright position and the other is tilted at a fixed angle, with the two moving toward each other with the same relative velocity as in the case of the spheres.

The same collision was repeated at increasing values of the elastic modulus E, and the maximum contact force was calculated with the different methods. The result of this study is reported in Fig. 3 and shows that, to model two impacting cartridges, it is necessary to couple Hertz theory of contact with Reissner model for thin shells. Taken individually, Hertz and Reissner theories clearly overestimate the contact force, proving that for colliding cartridges neither the shell bending-membrane deformation nor local Hertzian deformation can be neglected. On the contrary, if both stiffnesses are considered in series, a good approximation is obtained with the equivalent spheres. A direct DEM calculation (not shown in Fig. 3) matches those of the combined Hertz–Reissner model. As a consequence, the dynamics of the system can be studied properly by means of the discrete element method, taking into account equivalent spheres. As anticipated, the DEM calculation considering the realistic shape of the cartridges would turn out to be quite prohibitive from the computational point of view, with approximately two order of magnitude larger computational times.

3.2 Results for the dynamics of the spheres from DEM: choice of the optimal set-up

Fig. 4

Top view of our computational box with highlighted the tagged spheres used to probe the dynamics: triangular center (green), rectangular center (cyan), upper vertex (yellow), lower vertex (blue), first line (magenta). Also shown are the amplitude A and the direction of the velocity $V_{\mathrm{belt}}=20\,{\text {cm}}\,{\text {s}}^{-1}$ for the underlying belt (color figure online)

The discrete element method will provide indications for the optimal set-up, the velocity field of the equivalent spheres as a function of time, the evolution of the number of neighboring spheres in contact with a given one, identifying the most challenging scenarios to be considered for a detailed modelling by means of finite element method. We first consider the dynamics of 840 equivalent spheres, whose calculation can be carried out under realistic values of the glass elastic modulus ($E = 45\,\mathrm{GPa}$) within a reasonable computational time. The aim of this study is to identify the optimal set-up geometry in terms of possible contact forces and therefore stress peaks for the cartridges.

Figure 4 depicts the geometry of the computational box that represents the first part of the production line shown in Fig. 1a. In both cases (spheres and cartridges), all entities are inserted into a planar accumulation table confined by vertical walls and whose lower plane is moving with constant belt velocity $V_{\mathrm{belt}}$ driving them toward a constricted section ending into a small opening of amplitude A. The length of the tilted wall and the angle univocally define the geometry of the opening. In the real accumulation table, an external agent is acting beyond the opening preventing a clogging of the cartridges. In order to mimic this effect, the amplitude A is subject to a harmonic oscillation of period T. A number of spheres in different and strategic positions within the computational box were tagged and monitored during the dynamics, as shown in Fig. 4. We denoted as ‘triangular center’ (green), ‘rectangular center’ (cyan), ‘upper vertex’ (yellow), ‘lower vertex’ (blue) four such spheres and labelled them with different colours. In addition, we further tagged as ‘first line’ (magenta) all those spheres lying close to the constricted section.

Fig. 5

Selection of the optimal amplitude geometry with the force averaged over the first line (average force) as a function of time for different values of the parameters. a $V_{\mathrm{belt}}=20\,{\text {cm}}\,{\text {s}}^{-1}$, $T=1\,{\text {s}}$ and different amplitudes A; b $V_{\mathrm{belt}}=30\,{\text {cm}}\,{\text {s}}^{-1}$, $A=5\,{\text {cm}}$ and different periods T; c $T=2\,{\text {s}}$, $A=5\,{\text {cm}}$ and different values of $V_{\mathrm{belt}}$

We explored four different values of the amplitude $(A = 1, 2, 5, 7.5\,{\text {cm}})$, three different values of the belt velocity ($V_{\mathrm{belt}} = 10, 20, 30\,{\text {cm}}\,{\text {s}}^{-1})$, and three different periods $(T = 0.5, 1, 2\,{\text {s}})$ aiming at finding the configuration minimizing the average value of forces acting on the spheres. Here ’average’ means average over all considered spheres (coinciding with the force acting on that sphere in the case of a single sphere). The evolution of these forces was then monitored for all tagged spheres shown in Fig. 4. Figure 5a reports the result of the average force versus time for different values of the amplitude A ranging from 1 to $7.5\,{\text {cm}}$ for the ’first line’ spheres. The smallest possible case ($A = 1\,{\text {cm}}$) shows a significant increase of the average force at longer times presumably due to a clogging effect, since the spheres have diameters of about $1.6\,{\text {cm}}$. On the other hand, large values of the amplitude ($A = 5; 7.5\,{\text {cm}}$) give rise to very high spikes in the average force. This suggests the choice $A = 2\,\mathrm{cm}$ as optimal for the amplitude. Similar arguments, hinging upon Fig. 5b, c lead to the final choice of $T = 2\,{\text {s}}$, and $V_{\mathrm{belt}} = 20\,{\text {cm}}\,{\text {s}}^{-1}$. As mentioned, all forces in Fig. 5 were monitored in the first line spheres, but we have also computed the forces acting on all other tagged spheres with similar results. Interestingly, our final choice of $V_{\mathrm{belt}}$ turns out to be lower than the value $30\,{\text {cm}}\,{\text {s}}^{-1}$ used in the real production line. Our result then suggests that a lower value of $V_{\mathrm{belt}}$ would help in preserve the integrity of the cartridges, as one could expect on physical basis. These optimal values will then be used in all following calculations. A movie representing the dynamics under this condition is reported in Supplementary Information (Movie 2).

3.3 Universality of the coordination number in DEM calculations

Fig. 6

Coordination number as a function of time at increasing values of the Young modulus E for a all spheres; b ‘first line’ spheres; c ‘rectangular center’ sphere. In all cases $V_{\mathrm{belt}}=10\,{\text {cm}}\,{\text {s}}^{-1}$, $T=2\,{\text {s}}$, $A=2\,{\text {cm}}$

In close analogy with granular materials [7], we define the ‘coordination number’ as the number of neighboring spheres in contact with a given sphere. Figure 6 reports the results of this calculation as a function of the time for different choices of the tagged spheres and for different values of the Young’s modulus E. In Fig. 6a, b the coordination number is averaged over all spheres and first line spheres, respectively. Figure 6c shows the result for the ’rectangular center’ sphere that clearly displays an alternations between the two packing structures. Note that the initial configuration is a square packing, as visible from the snapshot of the initial conformation (see Fig. 4).

Two key messages can be learned from this analysis. Firstly, the average coordination number is nearly independent of the Young’s modulus E, thus suggesting that the way the spheres pack at the accumulation table does not depend on E. This means that the coordination number can be safely computed at low E, with much smaller computational endeavour. The second point is related to the structural rearrangement signalled by the oscillation between 4 and 6 coordination pertaining to a square or hexagonal packing, respectively. Note the oscillating period that is identical to the driving oscillation of the opening amplitude ($2\,{\text {s}}$), as one could expect from the outset.

3.4 Scaling of the velocity field in DEM calculations

Fig. 7

Velocity as a function of time for the “rectangular center” sphere a Young modulus range $10^7\,{\text {Pa}} \le E \le 10^9\,\text {Pa}$; b Young modulus range $11.25 \times 10^9\,{\text {Pa}} \le E \le 45 \times 10^9\,{\text {Pa}}$

Unlike the coordination number, contact forces are clearly strongly dependent on the Young’s modulus. We can infer this dependence by performing the same calculation—i.e. same initial conditions and same parameters other than the Young’s modulus—at different values of the Young’s modulus E, again as a function of time. Rather than reporting the magnitude of the contact forces, it proves more useful to report the magnitude of the velocity of the spheres. This is because the velocity field will be eventually input to the FEM analysis, so it is important to have a clear picture of its behaviour as a function of time.

We have explicitly monitored the velocities of all tagged spheres. Here we only consider the ‘rectangular center’ sphere because it best represents the ‘bulk behaviour’ of the systems. The other values can be found in Figures S3 and S4 of Supplementary Information. In view of the large difference in the values of forces depending on Young’s modulus, we have found expedient to gather plots in groups depending on the range of E values. Figure S3a of Supplementary Information reports the results for values $E=10^7,10^8,10^9\,{\text {Pa}}$ whereas Figure S3b of Supplementary Information refers to values of the Young’s modulus E from $11.25 \times 10^9\,{\text {Pa}}$ to $45 \times 10^9\,{\text {Pa}}$, the highest of which corresponding to the realistic value in the experimental set-up.

As visible in Fig. 7a, b, the ‘rectangular center’ sphere will reach the highest velocities after $\approx 50\,{\text {s}}$, with velocities ranging from $\approx 100\,{\text {mm}}\,{\text {s}}^{-1}$ for the lowest Young’s modulus value of $E=10^7\,{\text {Pa}}$ to $\approx 250\,{\text {mm}}\,{\text {s}}^{-1}$ for the highest (and realistic) value of $E=45 \times 10^9\,{\text {Pa}}$. These values will be used later on in the FEM analysis.

3.5 Velocities versus coordination number from DEM calculations

Fig. 8

a Coordination number and b velocity of the ’triangular center’ sphere as a function of time

Fig. 9

a Snapshot of the configuration after $9.3\,{\text {s}}$; b velocity of the ‘triangular center’ sphere as a function of time; c snapshot of the configuration after $12.2\,{\text {s}}$; d average velocity of the ‘first line’ sphere as a function of time (here ‘middle’ means that we do not include extreme values in the averages; In both b and d the right vertical axis also reports the corresponding values of the coordination numbers

Another important outcome of the DEM analysis that will be exploited in the following FEM calculation is associated with the two following questions. How is the velocity field distributed among the different spheres? And is the resulting velocity deriving from single or multiple contacts of the spheres?

We now address these issues by contrasting the resulting velocities with the corresponding coordination numbers for the highest value $E = 45 \times 10^9\,{\text {Pa}}$ of the Young’s modulus, since this is the actual value expected in the accumulation table. As a preliminary calculation, we performed a full span of the entire dynamics of all spheres, as well as of the ’first line’ spheres. This is done in Fig. 8 that displays the coordination number (a) and the corresponding average velocities (b) for the ’triangular center’ spheres. This is a rather useful information as it allows to associate the velocity of a given group of spheres with the corresponding coordination numbers at the same time. As before, the characteristic periodicity of both quantities is associated with the period of $2\,{\text {s}}$ of the amplitude oscillation, with the coordination number relaxing from a high value $\approx 6$ corresponding to hexagonal packing, to a nearly vanishing value when the exiting amplitude A is at its largest extension. The velocity also oscillates accordingly around $20\,{\text {mm}}\,{\text {s}}^{-1}$ until a marked increase is eventually reached when the ’triangular center’ spheres reach the outlet after $\approx 33\,{\text {s}}$.

Having identified the most interesting time intervals of each set of spheres, we then zoomed into them to extract specific patterns. This is done in Fig. 9a, b for the ’triangular center’ sphere, and in Fig. 9c, d for the ‘first line’ spheres. In first case, the ’triangle center’ sphere is the closest one to the outlet (labelled in green), and Fig. 9a shows a snapshot of the configuration after $9.3\,{\text {s}}$ of simulation, whereas Fig. 9b displays the corresponding dynamical evolution of the velocity (left axis) and of the coordination number (right axis). After an initial lag of roughly $2\,{\text {s}}$, where the sphere is nearly immobile into the initial hexagonal conformation, the sphere starts to move with bursts of activities having $2\,{\text {s}}$ intervals, with velocities and coordination number oscillating from $\approx 80\,{\text {mm}}\,{\text {s}}^{-1}$ and 5–6 respectively, to nearly vanishing values. The two quantities appear to be out of phase, however, as one could expect on a physical ground, with the highest velocity achieved when the coordination number is low and conversely lowest velocities associated with high coordination. After approximately $12\,{\text {s}}$ the ‘triangular center’ sphere exits from the amplitude and its dynamic is no longer of interest. This is also visible from the snapshot of the simulation (Fig. 9c) after $12.2\,{\text {s}}$ where the (green) ‘triangular center’ sphere is no longer visible. In this case, the interest switches to the ‘first lines’ spheres, and Fig. 9d reports the corresponding average velocity and coordination number. This set of spheres is particularly interesting because it was initially at the boundary of the constricted section with square packing of spheres on the rectangular part of the computational box, and hexagonal packing on the triangular side. In this case, there is an initial slow dynamics with maximum velocities of $30\,{\text {mm}}\,{\text {s}}^{-1}$ and coordination numbers fluctuating from 1 to 5. Here it is important to notice that the tagged spheres that were originally aligned become shuffled in the sea of all other spheres already after $9\,{\text {s}}$, as it can be clearly seen from the snapshot of Fig. 9a (taken after $9.2\,{\text {s}}$ of simulations).

3.6 The case of the cartridges from FEM analysis

By examining the dynamics of the sphere population obtained by means of the discrete element method, it was possible to detect the typical scenario in terms of positions, relative velocities, and coordination numbers. These results allow the setting up of a detailed Finite Element Model calculation that will provide a complete stress map of the cartridges, and will then allow to provide a risk assessment for the single cartridge.

All FEM calculation have been carried out using the Abaqus Code Package [10]. The cartridge mesh is made of three-dimensional, linear tetrahedral elements with three degrees of freedom per node, and the base belt is meshed by using rigid elements. This choice has been dictated by the specific geometry of the cartridge that would required a too dense mesh to describe the neck and the bottom of the cartridge with hexahedrons. We have explicitly tested in few specific cases, that this choice is sufficient for our purposes and that it is not necessary to use quadratic or higher order tetrahedrons. A friction coefficient of 0.2 is considered for the contact among cartridges, that means that the contact force has both a normal and a tangential component.

Fig. 10

Cartridge positions in the following cases a maximum relative velocity; b average relative velocity; c case of no springs; d case where no-tension springs are considered

4 Discussion of FEM results

Fig. 11

a Normal contact force versus time for pair 1–2; b tangential contact force for pair 1–2; c normal contact force versus time for pair 1–3; d tangential contact force for pair 1–3

Consider a single cartridge (labelled as 1 in Fig. 10) surrounded by additional 6 other cartridges (labelled as 2–7 in Fig. 10) arranged into a hexagonal conformation that, as we saw from DEM calculation, is the typical conformation where each cartridge is expected to experience the highest stress. We consider the characteristic scenario where two neighbouring cartridges (labelled as 2 and 3) collide with 1 with relative velocities $V_2$ and $V_3$ along a specific direction, and the remaining neighbouring cartridges (labelled 4–7) are in contact with 1. Based on previously reported DEM analysis, this configuration is one of the most representative among the innumerable possible situations that could take place on the accumulation table. We have then considered the two different cases of maximum relative velocities (Fig. 10a) where $V_2=180\,{\text {mm}}\,{\text {s}}^{-1}$ and $V_3=120\,{\text {mm}} \, {\text {s}}^{-1}$), and the case of average relative velocity (Fig. 10b) where $V_2=50\,{\text {mm}}\, {\text {s}}^{-1}$ and $V_3=25\,{\text {mm}}\,{\text {s}}^{-1}$) where, in both cases, the colliding angles have been obtained by the velocity vectors of the DEM analysis at selected positions. For each of the two above cases we have further considered two different conditions. In the first one, all the 7 considered cartridges are uncoupled from the others (Fig. 10c), whereas in the second one cartridges 4–7 are coupled through massless springs to neighbouring ’ghost’ cartridges (Fig. 10d). This latter case is meant to implicitly capture the non-linear effects provided by the set of remaining cartridges at the simplest possible level. Hence, in the first case the confinement provided by the adjacent cartridges is neglected (corresponding to a group of cartridges slightly isolated from the rest of the cartridges), whereas in the second case the external cartridges (4–7) are restrained by using no-tension springs mimicking the confinement provided by the additional cartridges. While we have performed the analysis for both the maximum (Fig. 10a) and the average (Fig. 10b) cases, only the most interesting case relative to the maximum velocity will be reported in the following.

Our interest hinges on the evaluation of the contact forces to estimate the risk of rupture for the cartridges under the worst possible conditions. Both the normal and the tangential contact forces are presented, highlighting their dependency on the initial position, relative velocity and mechanical properties of the cartridges. We focus our attention on pairs 1–2 and 1–3 that are those experiencing the collisions. Figure 11a shows the normal and tangential components of the contact forces for the 1–2 and 1–3 pairs as a function of time, starting from a conformation where both pairs are not in contact. Results are reported for both the case with and without external springs and forces are monitored on cartridge 1 which collides with cartridges 2 and 3 both at the bottom and at the top.

A marked difference between bottom and top impacts is observed, and the influence of the external springs is more effective on the 1–3 pair. Consider pair 1–2 depicted in Fig. 11a for the normal and in Fig. 11b for the tangential components of the experienced force. The bottom impact forces result significantly smaller than the corresponding top forces, both for normal and tangential components. The springs play a role after $1\,\mathrm{ms}$, only on the bottom impact: the collision takes place a little earlier with a slightly higher normal and tangential force. For the 1–3 pair (Fig. 11c, d), during the initial phase of the impact, the springs do not have any influence on the contact, whereas after about $1.2\,\mathrm{ms}$ the dynamics of the system dramatically changes. The confinement provided by the outer cartridges increases the contact forces, leading to considerably larger normal forces up to $55\,\mathrm{N}$ for the bottom contact and $80\,\mathrm{N}$ for the top contact. For tangential forces values as high as $10\,\mathrm{N}$ are achieved for bottom contact and $14\,\mathrm{N}$ for top contact.

This effect can be ascribed to the combination of constraints provided by the base belt, that is driving the whole population, and by the neighboring cartridges (4–7), thus resulting into a more and more blocking in the case of no-tension springs. This phenomenon has been observed repeatedly for several cartridge pairs, revealing that secondary impacts play a fundamental role. Within the real accumulation table, this situation can be associated to the dangerous case where some elements get stacked into a narrow space while the following are going to impact against them. Movie 3 of Supplementary Information shows the FEM dynamics in the case where external springs are present.

4.1 Stress field, strain field and microfracture

Fig. 12

a Maximum principal stress vs time for pair 1–2; b maximum principal stress through thickness of cartridge 1 for pair 1–2; c maximum principal stress vs time for pair 1–3; d maximum principal stress through thickness of cartridge 1 for pair 1–3

An essential theme in container integrity preservation is understanding the mechanisms behind the sudden failure and finding possible strategies for the mitigation of its risk. To this aim, a possible risk assessment can be made through the evaluation of the forces and peak stresses developed in the cartridge during its path on the accumulation table. It is very interesting to investigate the variation of the maximum principal stress through the thickness of the glass wall, highlighting its correlation with the impact force. As presented in the Method section, the contact load is directly linked to the relative velocity, the mechanical parameters of the container and in general to the initial conditions. This aspect is exasperated under real conditions on the accumulation table, where the chaotic dynamics of the cartridges causes several, continuously varying shocks. In Fig. 12a, c the maximum principal stress is shown as a function of time for pairs 1–2 and 1–3. When considering the confinement springs, for the pair 1–3 the peak stress reaches $140\,{\text {MPa}}$ at the bottom of the cartridge number 1. At the instant of the maximum peak stress, the corresponding variation through the thickness is shown in Fig. 12b, d. The worst situation is recorded for pair 1–3 with no-tension spring where the maximum principal stress reaches $94\,{\text {MPa}}$ inside the cartridge at the contact zone. It is interesting to observe a shifting trend of the stress, with increasing thickness depth subjected to tension stresses with consequent danger for crack initiation and/or propagation. Figure 13 reports the contours of the maximum principal stress of cartridge 1 in correspondence of the peaks of normal contact forces for the various cases. The black arrow indicates the impacting cartridge. It is worth to underline that the interior and exterior part of the interested impact zone are in tension and compression respectively. The possible crack initiates from inside the cartridge and propagates toward outside. This is a very dangerous situation, because the containers can be damaged by microscopic scraps which, in the worst case, can lead to the sudden rupture of the cartridge and cause the stop of the filling line, with high costs for the company.

In terms of risk assessment, by considering a realistic scenario based on DEM kinematics outcomes used as input for the Finite Element modelling, it is possible to extrapolate what could happen in the real accumulation table when the dynamical system is extended to thousands of cartridges.

Fig. 13

a Pair 1–3 bottom, no springs; b pair 1–3 top, no springs; c pair 1–2 bottom, no springs; d pair 1–2 top, no springs; e pair 1–3 bottom, non-linear springs; f pair 1–3 top, non-linear springs; g pair 1–2 bottom, non-linear springs; h pair 1–2 top, non-linear springs

5 Conclusions

The interest in the present paper was triggered by practical questions with far reaching consequences, but eventually evolved into an interesting general problem on its own right. The filling procedure of glass cartridges used in pharmaceutical applications poses increasing challenges to guarantee the integrity of the cartridge and the patient safety during the injection system. This process was mimicked in our study via a sequential combination of DEM–FEM calculation that led to the complete stress field of the cartridges. We used discrete element methods (DEM) to study the full dynamics of a large number of equivalent spheres in an accumulation table of variable geometry. We have motivated the use of equivalent spheres by showing that their dynamics under these conditions is representative of the dynamics of the cartridges (see Movie 4 in SI for a visual representation) that would otherwise be out of the reach of present computational capabilities. In addition to the optimal accumulation table conditions, DEM calculation provided the distribution of the number of contacts, the intensity and directions of the corresponding velocities and positions as function of time within the non-linear elastic collision model. These results were then used as input for the finite element method (FEM) calculation that provided the contact forces as a function of time, the stress tensor field, as well as the critical cracking zones for realistic values of the glass composition. Here, the likelihood of a crack initiation is evaluated based on the maximum principal tensile stress. It would be extremely interesting to implement a line of experimental investigations where the numerical predictions reported in the present work could be tested in simplified case studies. Among the possible available techniques, a promising example is given by photoelastic force measurements, which have been used with success in experiments with granular materials [24]. Our findings pave the way toward a general approach to these class of systems, thus transcending the specific case reported here.