Oblong True Shape Tablets as an Industrial Scale Digital Twin

1. Introduction and Motivation

1.1. DEM contributions

Many anhancements of the Discrete Emlement Method (DEM) are already achieved over the years. Starting with two-dimensional simulations of some hundred particles in the early 90th, nowadays threedimensional simulations are standard. Most of them are using spherical particles.

A leapfrog enhancement of DEM was possible with CUDA programming environment for the Graphics Processing Unit (GPU) [1]. This allowed to enter the million range using regular workstations [2]. Suggestions to the authors PhD mentor M. Kuna in 2010 resulted in a particle shape definition for polyhedra [3], adapted to GPU computation by E. Siegmann lateron [4].

The acceptance of digital twin tools depend on their capabilities to reflect real world problems. It is important to develop simulation tools which are capable to model realistic content. Here, in our current study, we focus on the shape of pharmeceutical tablets.

1.2. Tablet Coating

Unit Operation: Pharmaceutical tablet coating is a process used to apply a coating layer to a tablet to improve its appearance, taste, stability, and ease of swallowing. Tablet coating also helps to mask the unpleasant taste and odor of certain drugs and protects the tablet from moisture and light.

The tablet coating process involves several steps, which can vary depending on the type of coating and equipment used. Here is a general overview of the process:

Preparation: Before coating, tablets must be prepared by cleaning and polishing to remove any impurities or residues from the manufacturing process. Mixing the coating solution: The coating solution is prepared by mixing the coating material (usually a polymer) with solvents, plasticizers, and other additives. The solution is mixed until it is homogeneous and has the desired viscosity and density.

Coating application: The coating solution is then applied to the tablets using a coating machine. There are several types of coating machines available, including pan coaters, fluidized bed coaters, and spray coaters.

Drying: After the coating is applied, the tablets are dried to remove the solvents and other volatile components from the coating solution. This step is critical to ensure the uniformity and integrity of the coating.

Polishing: After drying, the tablets may be polished to improve their appearance and remove any rough edges. Overall, pharmaceutical tablet coating is a complex process that requires careful control of various parameters, including the type of coating material, coating thickness, and drying time, to ensure the quality and consistency of the coated tablets.

1.3. DEM - specialized for Pharmaceutical Processing

In the last decade DEM simulations have been widely used to analyze flows of granular media with the main emphasis on flows where the interaction between particles and a fluidic phase can be neglected. Extensions for extensively used industrial processes which are involve liquid spray modeled by ray-threading technique, for example film coating, dry impregnation and resulting wet imbibition, were invented and added successfully already by the author together with Rutgers Catalyst Consortium research group [5]. Here, liquids are handled implicit by altering friction, particle mass or by mimicking of liquid bridge force effects; Only particle-particle contacts are considered for mechanical interaction. The describes implementation works in-place, so whilst simulation runtime, properties (liquid content, particle mass, friction coefficients or cohesive forces) can change and alter the granular flow. Another approach to simulate spray is a kind of a droplet fog, shown in [6].

In contrast to well established simulation methods like FEM and CFD, DEM does not allow refinded regions for more detailed resolution. Often narrow slits or tiny gaps dictate the minimum particle size which has to be used all over the simulation volume. Moreover, users whom want to incorporate complex shaped particles have to choose between the number of particles and an acceptable simulation time which is lowered in most cases; Even though, a promising extrapolation approach can speedup many simulations significantly [7]. But in realistic pharmaceutical unit operations the number of tablets or capsules often surpass the million range in terms of particle numbers. This makes it even harder to incorporate geometrical details. Many promissing approaches for true shape particle modeling are covered by [8,9,10,11].

Some approximations toward tablets were carried out with superquadrics [12], [13]. But all of the flexible shape featuring superquadric objects have only smooth rounded edges. In the work of Yu et al., flat cylindical sharply edged objects are introduced [14]. The first true shape modeling approach for bi-convex tablets in 2D were found in the work of Song et al. [10]. Wassgren et al. demonstrated a lab-scale simulation study of bi-convex tablets [15]. Based on Song et al., the author and H. Kureck developed a 3D model and described the implementation in [16], which is tailored for GPU computing and runs similar fast than the simple sphere approach!

Industrial flow simulations at large scale, which require specialized CPU cluster hardware can be found in [12] and [17] but are often available to very restricted user community only. The shown examples in this paper are running on workstaion hardware on a single GPU.

1.4. pharmAssist DEM

In order to be able to simulate large granular systems, a high-performance DEM code was developed based on CUDA technology [18,19]. pharmAssist can simulate millions of particles on inexpansive GPU hardware as well as on any available Pro-Series Nvidia GPUs. It is a family of code packages of the so called "assist’s".

1.5. Numerical Precision

In order to serve the common double precesion capability of current GPU generations by one code base, the code is built using a ’guarded type definition’ of type real which unfolds to either type float or double together with all CUDA toolkit definitions of data types make_eal?() accordingly (?=$\{2,3,4\}$).

Commodity GPUs, widely used due to budget restrictions, have high single precision performance and can be used by the application of Froude theory as successfully shown in [20] and [21].

1.6. Visualization

DEM particle simulations in the million-particle range are still a challenge today but already common in some research groups. Vizualisation capabilities have to keep pace. Therefore, all available shapes (spheres, ellipsoids, polyhedra, capsules and tablets) have a shader-based rendering in our vizualisation tool. Our Kivy-based Pre- and Postprocessor has been tested to remain interactive up to several millions of particles on single a graphics card (Vizualisation presented in this paper was performed on a Nvidia RTX2080Ti using the VizAssist tool. Kivy is a crossplatform User Interface and supports also mobile devices running IOS or Android.

2. Implementation

2.1. Numerical Model

The wide spreaded ’soft’ particle approach is used in our work where spherical particles exert forces due to their mass m and geometrical overlap $\delta$ at the contact points c. The total force F and torque M acting on a particle

$$F=\sum _c{f}^c+mg$$

(1)

$$M=\sum _{c}r\times f^c$$

(2)

is the summation of all contact forces $f^c$ with other particles or boundaries and the gravitational force. By consideration of Newtons second law at time t, we obtain the acceleration $a\left(t\right)$ which can be used for an explicit time step integration scheme with a time step $dt$ in order to obtain the new velocities $v(t+dt)$ and corresponding positions $x(t+dt)$ for each particle in the system. Here, the classical Verlet integration scheme is employed for positions [22]

$$a\left(t\right)=F/m$$

(3)

$$x(t+dt)=2x\left(t\right)-x(t-dt)+a\left(t\right)d{t}^2.$$

(4)

Velocities are obtained by central difference approximation

$$\begin{array}{c}\hfill v(t+dt)=\frac{1}{2dt}\left(x(t+dt)-x(t-dt)\right)\\ \hfill +a\left(t\right)dt.\end{array}$$

(5)

The DEM accounts for translational as well as for rotational degrees of freedom. Both are handled separately by the application of the Verlet algorithm (3 - 5) and an Euler integration (7) in the latter case. Total torque M and moment of inertia I are used to obtain the angular acceleration $\dot{\omega }\left(t\right)$

$$\dot{\omega }\left(t\right)=\frac{M}{I}$$

(6)

and the new angular velocity at time $t+dt$

$$\omega (t+dt)=\omega \left(t\right)+\dot{\omega }\left(t\right)dt.$$

(7)

In order to obtain the new particle orientation $q(t+dt)$ an integration of the angular velocity is performed. Here, $q(t+dt)$ is a quaternion representation of particle orientation, see [23]. Details of the more complex implementation of a corrector-predictor method used in pharmAssist, can be found in [24].

The force laws in normal and tangential direction of our im[plementation follow the guidance of S. Luding, see [25]:

$$f_n=k_n{\delta }_n-D_n{\dot{\delta }}_n-f_{ch}$$

(8)

$$f_t=-min\left(k_t{\delta }_t-D_t{\dot{\delta }}_t,\mu \left|f_n\right|\right)·sgn\left({\delta }_t\right)$$

(9)

where $\mu$ is the static Coulomb friction coefficient, $k_n$, $k_t$ are spring constants and $D_n$, $D_t$ account for viscous damping in normal and tangential direction respectively. The cohesive force $f_{ch}$ is only active, if cohesive forces are considered, e.g., in the case of liquid bridges. The tangential overlap

$${\delta }_t={\delta }_t\left(t\right)+{\delta }_t(t-dt)$$

(10)

allows for restoring forces if the contact had existed already in the previous time step.Therefore, the evolution of the tangential overlap ${\delta }_t$ can be written also as

$${\delta }_t\left(t\right)={\int }_{t0}^t{\dot{\delta }}_t\left(\tau \right)d\tau .$$

(11)

In the algorithm, where the time of the lasting contact $[t-t_0]$ is discretized by k time steps dt

$$\begin{array}{c}\hfill {\delta }_t\left(t\right)=\sum _k{\dot{\delta }}_{t_k}\left(kdt\right)dt,\\ \hfill k=1,\left(1\right)\cdots [t-t_0]/dt.\end{array}$$

(12)

The reason for keeping the spring history by integration over contact time is to achieve tangential elasticity. Without tracking the evolution of ${\delta }_t$ it is impossible to reach a finite force where the tangential velocity ${\dot{\delta }}_t=0$ [26]P.13ff,[27],Ch.4ff. Tangential elasticity is a sufficient condition in order to simulate quasi-static and static granular flows. There, a particle assembly is moving with less granular temperature or it is resting, in both cases with respect to gravity. Quasi-static granular flows capture the most physical processes like mixing, compaction, filling and storage of granular matter. Without tangential elasticity a particle assembly cannot withstand gravitational forces by traction. It would be unstable until interlocking is taking place, which can –moreover– not occur when the particles have a pure spherical shape.

Sliding friction, a rolling resistance and twist spin about the contact normal vector is implemented in analogy. Collisions between particles and STL walls are handled in a similar fashion, comparable to the particle-particle collisions case [23].

2.2. Parallelization and Implementation

The pharmAssist code employes the common cell-space logic, a regular subdivision of the simulation volume into cubic cells [28]. This allows a fast detection of the adjacent neighbors, because the search for particles is restricted to a volume within a cell. While the search time for selecting two entities from a set of N entities is of $\mathcal{O}\left(N^2\right)$, the cell space logic has a linear growth characteristic of $\mathcal{O}\left(N\right)$ operations only. This is very important for DEM simulations when particle numbers exceed a few thousand. The cell-space implementation is tolerant in terms of particle size. Particles with an arbitrary size distributions within a limit of about 1:5 can be simulated without remarkable performance loss. Classical particle in cell (PIC) codes run a huge for-loop over the cells around a particle, collecting neighboring forces and integrating the resulting acceleration to obtain new positions and velocities.

A tree algorithm for STL handling follows the reference implementation with modified aligned bounding boxes for static, dynamic and free floating arbitrary shaped STL geometries (coating drums, railway sleepers, dummy objects etc) [29].

Method: Our software redefines the concept of parallelization of DEM with strong regards to the GPU capabilities. The light weight threads technology offered by programming models like Nvidia CUDA or AMD HIP allow for the simultaneous execution of up to hundreds of millions of such threads [19,30]; Morover, it is often required in terms of utilization of the potent hardware which fosters Teraflop performance.

Domains: Where domain decomposition of the cell-space is the common parallelization strategy on multi-CPU systems, where the simulation volume is divided into equal sized partitions for each of the CPU cores. The number of threads is equal to the number of domains. Such a coarse approach does not work on GPUs. Therefore, software cannot be simply ported to GPU, it has to be rewritten in total.

Thread-per-cell: The first good fitting approach is the split of the simulation volume into subsets of cells of the cell-space. Here, each cell is assigned to its neighborhood of 3x3x3 cells. Then, each of the cell is assigned to one thread and is ’looking’ for surrounding cells and contacting particles therein. So the number of threads is equal to the number of cells, typically reduced to the number of particles because N particles can only occupy a maximum of N cells (the particle center location defines the belonging). But on polydisperse systems number of particles per cell increases and performance is wasted in looping and search for more much than 3x3x3 particles in the latter case.

Thread-per-particle: This approach already allows a fast execution of the contact force computation but the varying number of contacts leads to unpredictable execution times for each of the thread-per-particle assigned threads (the one incorporating the most contacts delayes the execution of a whole threadblock). Typical contact numbers ranging from 0 – 8 contacts for each particle.

Thread-per-contact: Here, the situation is similar to the previous approach but a search step, which builds up a list of true contacts, is placed at the beginning of each of the time-steps. After the detection of contact pairs the contact list is prepared for seperate execution in a one-thread-per-contact implementation. As a result of an average of 5 – 8 contacts per particle (and many more blind contacts, when the bounding sphere of a non-spherical particle declares a potential contact) and the overhead of the blinds, twenty million up to hundred million threads are running in parallel for a simulation with a moderate size of only one million to five million particles.

2.3. Oblong Tablet Shape

Due to the manufacturing process in tablet press machines, most of tablets have a smooth cylindric section alongside the compaction axis, the so called band. The oblong shape of a tablet is characterized by a circumferending band with opposite parallel sections at the long sides, see Figure 1.

Some tablets land on their band and stand on the plane which one can observe in Figure 1. Here, two out of around thirty tablets stand upright. Such rather simple effects can only be reflected in particle simulations, when geometry features of the particles are modeled accordingly. In Figure 2 a similar simulation case of some tablets are dropped onto a rigid plane under gravity shows comparable behavior (25 out of 455 tablets land upright).

The geometry definition we use in the present work is an extension of the simple but smart idea of Song et al. [10], outlined in Figure 3. We extend this approach by adding a further L-parameter, namely $L_4$. Where $L_1$, $L_2$ and $L_3$ are corresponding to band hight, total hight and band width of the tablet, in the oblong shape extension $L_4$ represents the band lenght or total length of the oblong tablet. A fix set of the shape parameters $\{L_1,L_2,L_3,L_4\}$ defines an oblong tablet shape. Thereafter, such a tablet can be scaled using one size parameter $R_0$ only, see Figure 3.

The main idea of oblong particle shape definition is to use the bi-convex Song-tablet in Figure 3 as countour shape object (CSO). A contour shape is created by moving around such a geometrical defined CSO. Very similar to spherocylinders, which form by shifting a perfect shpere along a line, see Figure 4. Details of spherocylinder collision are given in [31]. This shape serves to model pharmaceutical capsules. In short, the collision between two capsules is handled by the collision of the two spheres at the position where both capsule axes are having its shortest distance.

When we simply replace the sphere in Figure 4 with a bi-convex tablet, shifting along the capsule axis, the oblong tablet shape is created, see Figure 4. Whereby the spherocylinder forms round caps at both ends the oblong tablet instead has ends of exactly the shape of a half bi-convex tablet.

All pairwise contact conditions for the oblong shaped tablets remain the same than for pairwise bi-convex tablet collisions, explicitly given in [16].

The only difference is the handling of band contacts, which are treated as capsule collisions using their cylindrical part. This modeling enables stapling of tablets at band-band contacts, see Figure 5 and simulates the (unwanted) effect of tablet climbing out of the box, see Figure 6.

3. Basic Application Example

3.1. Experimental Device

We use a rotating drum for our simulations, sketched in Figure 7, which is widely used in unit operations in pharmaceutical processing or food industries such as tablet coating or pan coating candy production. Drum coaters have been investigated to great extent for non-spherical particles by Hlosta et al. [32]. Their findings regarding fill levels depending on the asphericity of particles helped us to set up our simulation case.

In this synthetic case study we use an STL-based geometry of an empty drum without baffles or scrapers. The drum is filled up to 40% filling height with 320.000 oblong tablets, shown in Figure 7. The total mass of the batch is 160kg.

3.2. Simulation Setup

In order to study the influence of drum speed to the contact pair statistic, four different simulations have been ran. This simulations ranges from 6rpm to 18rpm and were carried out as shown in Figure 8.

For all cases from 6 rpm to 18 rpm, the simulation ran for 10 seconds. So our cases made 1, 1.5, 2 and 3 revolutions resp. The tablet size was 8.6 mm radius ($R_0$ in Figure 3). The shape defining set for a unit tablet with a length of 1 is given by $\{L_1=0.15479,L_2=0.23972,L_3=0.4,L_4=1.0\}$. Using this dimensions we obtain a tablet mass of 500mg.

Simulation parameters are given in Table 1.

4. Results

4.1. Contact tracking

In order to get a data base for statistical analysis of frequency of the occurance of contacts, we have markers in the simulation code which are written for each contact pair and which are summarized for each contact type. Whilst the simulation run, every 0.1s the histograms data are stored.

4.2. Contact pair statistics

Simulations allow to access data not available by experimental means at all. Having complex shaped tablets available, we decided to make a classification on different parts (cap, band rim) of an oblong tablet, depicted in Figure 5. The resulting six contact pair combinations are counted at time step intervalls of 0.01s. We labeled the pairs band-band, band-cap, cap-cap, rim-band, rim-cap and rim-rim. The histogram data are computed in-processing on the GPU and have no significant influence of the total performance.

All histograms in Figure 9 and Figure 10 are arranged as stacked barcharts, so one can compare the evolution of different contact pairs.

Different evolutions of contact pair variability can be observed up to four seconds of drum operation, cauesed by the different rpm of the drum in each case. Higher rpm values leads to higher fluctuations after the drum starts to operate. From four seconds onwards the evolution is steady with no remarkable visual differences between the cases running 6 rpm to 12 rpm. One can observe that cap-cap pairs dominate; followed by rim-cap pairs. The probably most damage causing rim-rim contacts are very low in all cases (around 2%). The fastes case has not reached steady pair distribution yet. The black, red and green pairs show an ongoing shift within each other. In contrast, all pairs with rim involved (upper three colors), are already steady after three seconds of simulation, see right diagram in Figure 10.

In order to reach better comparison, we plot data after the 10 seconds together in one diagram, see Figure 11. Now some interesting findings occur when we compare how the situation evolves with increasing drum speed in terms of their frequencies (comparing heights of the same colored sections). The absolute %-values for each of the pair frequencies are provided in Table A2.

Tendencies of the contact pair variability are summarized in Table 2. Further enhancement of the statistics might involve contact forces multiplied by contact duration or collision frequency, in order to describe wear depending on contact pairs etc.

4.3. Performance

The presented simulations (320.000 oblong tablets) ran for 10 seconds with a speed of 155 steps/s. Wall clock time to reach the 10 seconds was 1.8h, SP@RTX2080Ti). ( 121 steps/s DP @GV100). A comparable case with 1 Million oblong tablets ran with 41 steps/s DP@GV100, see also Table A1.

5. Industrial Scale Application

Coating processes as a pharmaceutical unit operation for bare tablet cores typically ran up to several hours. Whilst processing, the tablet bed is countinuously mixed in radial direction caused by circulation of tablets between the free surface of the bed and the cylinder wall. In order to enable axial mixing along the drum axis, multiple baffles are placed around the inner wall of the drum. The so improved mobility of the tablets decreases unwanted twinning, an effect where freshly coated wet tablets leaving the spray zone and stick together with another tablet, breaking up lateron again and remain coating film imperfections. For the mitigation of wetting effects, hot air streams through the tablet bed, shortening the drying phase. Therefore, the walls of such drum coaters are completely perforated like screen door windows. Together with the flow through property the perforation leads to a nearly frictionless transport of tablets by the moving drum wall, reflected by a high friction value (table ). Neither airflow, heating and spray are considered in the presented simulation case.

The used device geometry of the industrial simulation case is approximated to real devices which oparate batch sizes in the hundreds of kg range. The shape of the baffles is usually patented and therefore also simply visually adapted in the spirit of having smoothed internals guiding gently the tablet flow, see Figure 12.

Whilst mixing, the contact type pair statistic is recorded in order to evaluate the situation after steady state is reached.

In Figure 14 the pair-wise interaction frequency is depicted as it evolves over simulation time. The dominant pair is the cap-cap contact followed by band-band and rim-cap contact pairs. Less frequent are band-cap and rim-band pairs, rim-rim pairs are least frequent but can be most harmful with regard to shape degradation like edge imperfactions and notches.

Each granular assembly is characterized by its number of contacts between the granules. The so called coordination number C. In our case, oblong tablets are the individual objects, which interact by inter-tablet forces throughout these contacts. Where monodisperse spherical granular packings are known to have coordination numbers of about C = 5-6, the elongated tablet shape has higher coordination numbers ranging between C = 7-8, shown in Figure 15.

Bounding sphere contacts are potential contacts which are evaluated in the contact force computation part of the simulation code. Over 50 % of the contact calculations effort are wasted but are unkown to rasult in a true contact upfront the computational task. In the discussed case the total inter-tablet contacts are larger than 18 millions and the blind contacts caused by the bounding sphere detection reaches 50 million total contacts to compute in each individual timestep.

In order to get a more detailed analyzis of the acting inter-tablet forces, we give a summary of the probability density split into the six contact pairs, shown in Figure 16.

Each of the pairwise contact sets have their own forces and their distribution. This six differend probability density functions derived from histogram data by running averages) show different contributions. Contact pairs incorporating the tablet band tend to have a higher level of contributing forces, their curves are arranged higher than those for rim contact pairs. The extracted maximum normal force magnitudes at different times of the simulation are given in Table 4.

The observed maximum forces in Table 4 correspond to the decreasing force levels for the contact type pairs, extracted from the force distribution curves at 10s in Figure 16. However, wear and notches at contacts depends on acting force, the corresponding contact area as well as on the penetration depht, conclusions to real damage properties is an open question.

6. Summary

We have investigated the influence of mixing speed to various contact pairs forming whilst the continous flow of tablets in a rotating drum. To show capabilities of getting run-time data out of a simulation, we plotted histograms for contact pair statistics. The Results show no groundbraking effects but give already an idea how a digital twin might be used to analyze specific unit operations in tablet processing.

Future development work of the pharmAssist code will focus onto the modeling of other commonly used tablet shapes. Enhancements of analytical methods might extend to wear and tear analyzis, contact duration, trace of tablet cycling, occurence on free surface, flip rate (which side looks to the spray gun, how often), fill-level effects etc.

Appendix A.1

Table A1. Performance data for an 1 million oblong tablet simulation using a 40k STL geometry at seddled particle density and fully employed contact list.

	steps per second	/(step time for a full DEM cycle)
GPU type	SP	DP	SP | DP / total Power
RTX280Ti	30	12	243 | 155 /250W
P4000	17	2	101 | 60 /105W
GV100:	88	41	245 | 250 /250W
A100:	121	59	236 | 240 /400W

Table A1. Performance data for an 1 million oblong tablet simulation using a 40k STL geometry at seddled particle density and fully employed contact list.

	steps per second	/(step time for a full DEM cycle)
GPU type	SP	DP	SP | DP / total Power
RTX280Ti	30	12	243 | 155 /250W
P4000	17	2	101 | 60 /105W
GV100:	88	41	245 | 250 /250W
A100:	121	59	236 | 240 /400W

Table A2. Relative freqiencies of contact pairs in %.

drum rpm	band-band	band-cap	cap-cap	rim-band	rim-cap	rim-rim
06	15.7262	18.7596	28.7395	12.2372	20.8561	2.14923
09	16.1118	17.092	29.7451	12.2076	21.4282	2.16562
12	15.9853	16.7705	29.7367	12.3117	21.8448	2.23866
18	16.7701	15.0783	31.1261	12.1837	21.7753	2.18817

Table A2. Relative freqiencies of contact pairs in %.

drum rpm	band-band	band-cap	cap-cap	rim-band	rim-cap	rim-rim
06	15.7262	18.7596	28.7395	12.2372	20.8561	2.14923
09	16.1118	17.092	29.7451	12.2076	21.4282	2.16562
12	15.9853	16.7705	29.7367	12.3117	21.8448	2.23866
18	16.7701	15.0783	31.1261	12.1837	21.7753	2.18817