Particle-Scale Modeling to Understand Liquid Distribution in Twin-Screw Wet Granulation

Abstract

Experimental characterization of solid-liquid mixing for a high shear wet granulation process in a twin-screw granulator (TSG) is very challenging. This is due to the opacity of the multiphase system and high-speed processing. In this study, discrete element method (DEM) based simulations are performed for a short quasi-two-dimensional simulation domain, incorporating models for liquid bridge formation, rupture, and the effect of the bridges on inter-particular forces. Based on the knowledge gained from these simulations, the kneading section of a twin-screw wet granulation process was simulated. The time evolution of particle flow and liquid distribution between particles, leading to the formation of agglomerates, was analyzed. The study showed that agglomeration is a rather delayed process that takes place once the free liquid on the particle surface is well distributed.

1. Introduction

An accurate description of aggregation and breakage rate expressions is of great interest for the high shear wet granulation research and development. Specifically, this interest is fueled by the need for rational process design, i.e., to produce granules with the desired quality and maximize yield. Mechanistic modeling is seen as a powerful tool to develop this understanding regarding the wet granulation process and equipment design. Many mechanistic expressions, i.e., so-called rate kernels, have also been developed in this regard, describing aggregation and breakage rates [1]. However, these expressions often require the specification of different material properties, such as wetness, porosity, and yield strength. The mesoscale frameworks like population balance equations (PBEs) typically require closures for (i) collision rates, (ii) particle collision velocities, or data on other particle-scale phenomena. Additionally, material fluxes (i.e., the axial and radial velocities) need a suitable description for successful dynamic simulation of the mixing process. Unfortunately, the particle-scale information for these closures is known to be dependent on process parameters, equipment geometry, or material properties that cannot be directly computed. Hence, a common strategy is to calibrate these unknown parameters required in the closures using experimental data [2]. Despite calibration, such closure models have limited applicability due to the lack of predictive capabilities outside of their experimental domain [3]. This is unsatisfactory, especially if one aims for process optimization or change the granulator’s geometry to improve the performance.

In contrast, the discrete element method (DEM) is a micro-scale framework that tracks individual particles, or granules, as they move through space and time. This highly-accurate physical modeling tool can provide details on the flow of particles or granules. This information, e.g., collision rates and velocity profiles, can then be fed into the rate kernels, incorporating effects due to current process parameters, the equipment geometry, and material properties. Based on the level of detail, this can further be classified as particle-scale and granule-scale simulations, finding different applications.

The “particle-scale” simulations can be used for extensive investigation of particle-scale phenomena such as mixing of entities since particles in a granular assembly are treated as a separate entity and tracked. Previous attempts of particle-scale simulations include the work of Lian et al. [4] and McCarthy [5]. However, this also makes such simulations considerably more expensive than the granule-scale DEM-based simulations as a large number of entities (i.e., particles) have to be tracked.

DEM-based simulations that track individual granules (i.e., ensembles of particles) are frequently used [3,6,7]. This type of simulation is referred to as “granule-scale” since phenomena on the scale of the granules, typically with a size of 100 to 1000 µm, can be directly simulated. The advantage of these simulations is that a comparably small number of entities (i.e., granules) have to be tracked, and that granules are allowed to break or aggregate. This opens the door to investigate the effect of the actual geometry of a twin-screw granulator (TSG). Also, variations in granule properties, such as porosity and liquid content, are often accounted for in this type of DEM simulations using empirical correlations for, e.g., the coefficient of restitution or Young’s modulus [8]. Closures for mechanisms such as aggregation, breakage, and consolidation, which cause changes in particle size and other properties, need to be provided in granule-scale DEM simulations. Therefore, often a multi-scale approach is suggested where DEM data is collected and used within a population balance model (PBM) by either one-directional or bi-directional coupling [7]. Most of the multi-scale studies, such as Gantt et al. [9] used one-directional coupling of DEM with PBM to evaluate mechanistic coalescence kernels. Goldschmidt et al. [10] used DEM simulations to solve a PBM, replacing small particles with larger ones as they successfully coalesce. Bouffard et al. [11] used DEM to evaluate a spatial transport in a compartmental PBM. Additionally, Reinhold and Briesen [12] developed a coupled PBM-DEM model for wet granulation in which DEM simulations were used to evaluate a mechanistic aggregation rate kernel. Dhenge et al. [13] first used DEM to investigate surface velocity of powder in a TSG with conveying screws. Recently, Zheng et al. [14] used a GPU-enhanced DEM simulation for investigating the transport of dry particle systems through a TSG. Barrasso et al. [3] implemented a bi-directional coupling between PBM and DEM to evaluate collision frequencies and liquid distribution as a proof-of-concept. The model showed sensitivities to the screw configuration, process parameters such as screw speed, liquid-to-solid ratio, and material properties such as binder viscosity and pore saturation. However, this bi-directional coupling of a PBM-DEM model for wet granulation depends on the assumption that liquid is distributed evenly with respect to particle volume, which is not the case as per experimental studies [15]. Recently, Tamrakar et al. [16] integrated the complex interactions such as powder wetting behavior, capillary and viscous liquid bridge formation as well as binder dissolution within DEM to simulate two binder addition approaches—wet and dry, during high shear wet granulation (HSWG). However, particle level effects, like breakage or attrition, were omitted from the DEM model. In fact, achieving a homogeneous distribution of liquid, which largely controls the constitutive mechanisms of the twin-screw wet granulation, is the primary requirement for driving the granulation process towards a required state [17]. Therefore, investigating the drivers for particle and liquid mixing in the TSG is key to our understanding of the twin-screw wet granulation process.

In our present study to investigate particle and liquid mixing in the mixing zone of a TSG, our focus was first to track individual particles in a short quasi-two-dimensional simulation domain. Our simulations incorporate detailed models for liquid-bridge formation, rupture, as well as the effect of the bridges on inter-particular forces. Thereafter, the same simulation approach was applied to the mixing zone of a TSG, and the effect of process parameters on particle and liquid mixing rate was investigated.

2. Particle Scale Modeling Approach

2.1. Particle Flow Model

Computer simulations using the DEM were performed using the software package LIGGGHTS [18]. In a dense particle system forces acting on a particle can be decomposed into a body force $F_i^b$ on particle i, as well as into a contact $F_{ij}^{con}$ and cohesion force $F_{ij}^{coh}$ between particles i and j. The interaction of the granular flow with the surrounding air has been neglected in our work due to the high particle fill ratio and the large particle-to-air density ratio. The so-called soft-sphere method, in which collisions are modeled as enduring contacts, was used. The hypothetical overlap between particles, occurring over a few time steps, is used to model the particles’ deformation. Contact forces are a function of these overlaps. The forces are expressed with the use of a spring, dash-pot and slider model, which separate forces into normal and tangential forces, as shown in Figure 1. Based on these forces, the trajectories of individual particles were computed using Newton’s equation of translational and rotational motion.

The contact force $F_{ij}^{con}$ acting on particle i, from particle j can be decomposed into a normal ${\mathbf{F}}_i^{c,n}$ and a tangential component ${\mathbf{F}}_i^{c,t}$, which are modelled as:

$${\mathbf{F}}_i^{c,n}=k_n\delta {\mathbf{n}}_{ij}-{\eta }_n{\mathbf{v}}_{ij}^n$$

(1)

$${\mathbf{F}}_i^{c,n}=-k_t\delta {\mathbf{u}}_{ij}^t-{\eta }_t{\mathbf{v}}_{ij}^t$$

(2)

where ${\mathbf{u}}_{ij}^t$ is the tangential overlap and $\delta$ is the normal overlap between particles i and j:

$$\delta =r_i+r_j-|{\mathbf{r}}_i^{*}+{\mathbf{r}}_j^{*}|$$

(3)

By Newton’s third law, particle j experiences the same contact force as particle i but in the opposite direction. The characteristic contact time is $t_c$ given by

$$t_c=\pi /\omega$$

(4)

where

$$\omega =\sqrt{k_n/m_{\mathrm{eff}}-{\eta }_n^2/\left(4m_{\mathrm{eff}}^2\right)}$$

(5)

This contact time sets the time step $\Delta t$ for the integration of Newton’s equation of motion. The restitution coefficient $e_n$, which is the ratio of relative velocities of the particle after and before collision, is approximated as:

$$e_n=\exp \left(\frac{-\pi {\eta }_n}{m_{\mathrm{eff}}\sqrt{4k_n/m_{\mathrm{eff}}-{\eta }_n^2/m_{\mathrm{eff}}^2}}\right)$$

(6)

to calculate the damping ${\eta }_n$.

2.2. Liquid Bridge Model

When the granulation liquid, denoted as liquid phase in what follows, is added to dry powder particles in the TSG, the surface of the particles is wetted, and the liquid on the particle surface is distributed between the particles in contact. Liquid distribution in the moving dense particle system progresses due to (i) liquid transfer facilitated by the particle-particle collisions and (ii) the convective transport of the particle motion by throughput force of incoming material and drag force of the rotating screws. Depending on the ratio between solid and liquid volume, different wetting states are achieved [19]. The cohesive forces between the wetted particles strongly depend on the exact shape of the liquid bridge and the pressure difference between the gas and liquid phase. The stability of a liquid bridge between two particles and its rupture thus mechanistically dictates the macroscopic granulation status. For simulating cohesive forces, the capillary force model proposed by Mikami et al. [20] was used. In the case of bridge formation, model C proposed by Mohan et al. [21] was used for tracking the liquid distribution. In this model, when wet particles touch each other, a certain fraction of the liquid on the particle surface ${\varphi }_{tr}$ was exchanged at a finite rate $Q_{tr}$ calculated as

$$Q_{tr}{}_{,i}=-\frac{c_f}{t_{\mathrm{ref}}}(L_{p,0}{\varphi }_{tr}-V_b)$$

(7)

where $c_f$ is a dimensionless filling rate coefficient and $t_{\mathrm{ref}}$ is the reference liquid bridge filling time calculated as

$$t_{\mathrm{ref}}=\frac{r_{\mathrm{eff}}{\mu }_l}{{\sigma }_l},\mathrm{with}{r}_{\mathrm{eff}}=\frac{2r_i{r}_j}{{(r_i+r_j)}^2}$$

(8)

and,

$${\varphi }_{tr}=1-\sqrt{1-\frac{r_j^2}{{(r_i+r_j)}^2}}$$

(9)

For a more detailed discussion regarding the liquid transfer model see the study by Mohan et al. [21].

Liquid Loading, Bridge Volume Fraction & Liquid Bridge Coordination Number

In a TSG, liquid is added to a region of the granulator which is later distributed by the shearing motion of the particles. The amount of liquid load on the particles in the liquid addition region was calculated as

$$Q_{lod}=\frac{Q_i}{\frac{\pi }{6}{n}_{p,liq}{d}_p^3}{t}_{\exp }$$

(10)

At the steady state condition, the liquid bridge volume is given by

$$V_b=\frac{1}{2}(L_p{}_i+L_p{}_j){\varphi }_{tr}$$

(11)

The global-average fraction of liquid in liquid bridges relative to the total amount of liquid in the system, i.e., the liquid bridge fraction, was calculated as:

$$V_{bf}=\frac{{\sum }_{i=1}^{n_p}{V}_b{}_i}{{\sum }_{i=1}^{n_p}(V_b{}_i+L_p{}_i)}$$

(12)

As the liquid tracking model assesses the per-particle liquid content, the number of liquid bridges and particle-particle contacts at each time step was also tracked. This information can be used to compute $n_{b,i}$ i.e., the average number of liquid bridges per particle as:

$${\overline{Z}}_b=\frac{2{\sum }_{i=1}^{n_p}{n}_{b,i}}{n_p}$$

(13)

2.3. Simulation Set-Up and Input Parameters

2.3.1. Simple Periodic Simulation Box

In order to investigate the solid-liquid mixing under well-controlled flow conditions, an assembly of about 1500 mono-disperse particles (the exact number varied based on the fill ratio) of diameter $d_p$ and density ${\rho }_p$ were placed in a periodic box of equal length and width ($H/d_p=20$) at particle volume fractions ranging from $0.3$ to $0.5$ (see Figure 2 for a typical setup). A homogeneous shear flow, with the shear gradient pointing in the x-direction, was applied based on Lees–Edwards boundary conditions [22]. The y-direction is the spanwise direction. The reference data for the rheology in this set-up is available without wall effects because each point in the simulation domain is statistically identical. The dimensionless shear rate (Equation (14)), which is defined as the ratio of a characteristic deformation time of the particles and a characteristic shearing time, was used to characterize the shearing rate:

$${\gamma }^{*}=\gamma d_p^{3/2}/\sqrt{k_n/{\rho }_p}$$

(14)

For soft particles, ${\gamma }^{*}$ $>{10}^{-1}$ [23], thus the simulations were carried out for a fixed ${\gamma }^{*}$ of 1.

A gradient of liquid content on the particles was imposed in the middle zone of 2$d_p$ width in the periodic box, both in positive and negative y-directions. The liquid content of the particles in this zone was fixed, while keeping other particles dry (i.e., liquid content and bridge volume equals zero). The liquid content on the particles was defined based on a dimensionless reference film thickness $ϵ$ = $h_0$/$r_{\mathrm{eff}}$ where, $h_0$ is proportional to the reference liquid content on the particles via $L_{p,0}$ = $d_p\pi h_0$. A typical snapshot of such a simulation at the initial state (i.e., time $t_0$) is presented in Figure 2. The material properties applied in these simulations are based on measurements on lactose particles by Perkins et al. [24] and using water as a granulation liquid. The simulation parameters are listed in

Table 1. Summary of input parameters used in particle-scale modeling.

Quantity	Symbol	Value	Unit
Particle diameter	$d_p$	1.00E-03	[m]
Young’s modulus	G	3.45E+9	[N/m2]
Initial particle velocity	$v_x$, $v_y$	1, 0.1	[m/s]
Coefficient of restitution	$e_n$	0.9	[–]
Coefficient of friction	µ	0.1	[–]
Poisson ratio	$\nu$	0.33	[–]
Film thickness	$h_0$/$r_{\mathrm{eff}}$	1.00E-02	[–]
Dimensionless filling rate coefficient	$c_f$	1	[–]

. The spring stiffnesses and damping coefficients in normal and tangential direction were calculated from these parameters considering the LIGGGHTS User Manual for the chosen contact model [25]. The timestep size during simulations was selected to correspond to the timestep of 20% of the Rayleigh time criterion, $d{t}_r$ calculated as:

$$d{t}_r=\pi r\frac{\sqrt{{\rho }_p/G}}{0.1631\nu +0.8766}$$

(15)

where ${\rho }_p$ is particle density, G is the shear modulus and $\nu$ is Poisson’s ratio.

2.3.2. Mixing Zone of a TSG

To investigate the mixing of granulation liquid and formulation powder in the mixing zone of a TSG, DEM simulations were performed using the approach discussed in Section 2 and knowledge gained from the base case of the simple periodic simulation box (Section 2.3.1). A 3D CAD model for the kneading discs and barrel were created in SALOME and the triangulated surface mesh (STL) was created using Gmesh. The geometry was based on the screw design which is used in the TSG of the ConsiGma-25 unit (GEA Pharma Systems, Collette™, Wommelgem, Belgium) (Figure 3). The screws were rotated at a fixed speed. Approximately 46,000 mono-disperse particles of diameter $d_p$ and density ${\rho }_p$ were placed in the free space between the fixed-mesh of barrel and the rotating-mesh of two kneading discs. A periodic interface was created for the axial boundaries of the mixing section to allow particles that exit the section to enter from the other side. This periodic interface avoids to have walls as axial boundaries that would induce potential negative effects (e.g., wall layering).

3. Results and Discussion

3.1. Solid-Liquid Mixing in the Simple Periodic Simulation Box

The solid-liquid mixing in a dense particle system was first investigated by the DEM simulation for a sheared assembly of particles discussed in Section 2.3.1. Shearing of particles resulted in (i) convective transport of liquid by particle motion (observed as liquid distribution), and (ii) the transfer of liquid from the particle surface to the liquid bridges (and vice versa) which is considered as conductive transport. The latter transport mechanism also results in the formation of agglomerates (Figure 4).

To investigate the effect of changes in the bulk density and fill ratio of the powder bed inside a TSG barrel, the volume fraction of the particles was changed in the simulations. To investigate the effect of liquid addition rate during granulation, the initial liquid loading on the particle in the wetting zone was changed in the simulations. The effect of liquid addition method at one point or two points, which basically affects the width of the addition zone, was investigated by changing the width of the wetting zone. The response of all these investigations on mixing and granulation was recorded in terms of change in the average number of liquid bridges per particle ${\overline{Z}}_b$ and the fractional amount of liquid in bridges $V_{bf}$. An increase in ${\overline{Z}}_b$ indicates an increasing aggregation level. Normally, the trends for ${\overline{Z}}_b$ and $V_{bf}$ should follow the same pattern. However, a change in their relation indicates the dominance of convective or conductive transport during wet granulation. Moreover, an increase in the screw speed leads to an increase in the shear level, which was investigated by increasing the applied shear to the particle assembly in the shear box.

3.1.1. Effect of Change in Volume Fraction of Particles

The filling ratio inside the TSG varies significantly due to change in the screw speed, material throughput, and flow restriction by the kneading elements. In order to understand the effect of change in the fill ratio using simulations, the volume fraction of the particles in the simple periodic shear box was varied in the range of $0.3$ to $0.5$. Increasing the fill ratio resulted in a higher number of particles in the wetting zone (Figure 2). This is directly reflected by more liquid bridges for the high fill ratio at the beginning of the simulation (Figure 5a). However, as the shearing progressed ${\overline{Z}}_b$ increased for conditions with lower volume fractions ($\varphi =0.3,0.4$) whereas for the high volume fraction ($\varphi =0.5$) it first reduced and then equilibrated to gain some increase. The volume fraction of liquid in bridges $V_{bf}$ also showed the same patterns where the liquid bridge fraction was higher for higher $\varphi$, and it always remained so compared to the lower $\varphi$ despite shearing. However, while $V_{bf}$ equilibrated as shearing progressed for higher volume fraction ($\varphi =0.4,0.5$), it kept on increasing for low volume fraction ($\varphi =0.3$). Additionally, at a high fill fraction, while ${\overline{Z}}_b$ was lowering by initial shearing, $V_{bf}$ increased. This indicates shear elongation, thus more liquid in the capillary bridges followed by breakage of these bridges.

The results indicate that with a high fill ratio, the (bridge and agglomerate) breakage mechanisms dominate over agglomeration, which results in redistribution of liquid. However, at a low fill ratio, agglomeration is the principal mechanism of progress in granulation. In order to keep both agglomeration and breakage as the active mechanism of granulation during simulation, hereafter $\varphi$ has been fixed to intermediate level ($\varphi$ = 0.4) unless stated otherwise.

3.1.2. Effect of Change in Liquid Loading on Particles

In a TSG, the amount of liquid loaded on the first wetted particles in the wetting zone changes when liquid addition rate or the fill fraction is changed. The shear mixing distributes this differently liquid loaded particle to proceed to a certain granulation end-point. This was simulated by increasing the liquid loading of the particles in the wetting zone $Q_{lod}$ (Equation (10)). The effect of an increase in liquid loading resulted in an increase in both ${\overline{Z}}_b$ and $V_{bf}$ (Figure 6). At a low $Q_{lod}$ (3.75 $\times {10}^{-4}$), after an initial increase, ${\overline{Z}}_b$ equilibrates at a certain level and later started to increase (Figure 6a). However, at higher levels of $Q_{lod}$ (7.5 $\times {10}^{-4}$, 1.5 $\times {10}^{-3}$) ${\overline{Z}}_b$ kept on increasing. Interestingly, although $V_{bf}$ had already achieved equilibration (Figure 6b), the ${\overline{Z}}_b$ was increasing. The ${\overline{Z}}_b$ increased further at high shear and this elevation was more at higher $Q_{lod}$.

This indicates that when more liquid is available, the liquid on the particle surface gets rapidly transferred to the liquid bridges. At low liquid loading, in the attempt to form more bridges by changing rearrangement, some liquid is returned back to the particle surface. However, in particles with higher loading, the initial liquid bridges rearranged themselves to form more bridges without giving away the liquid to the particle surface.

3.1.3. Effect of Change in Liquid Addition Zone Width

The width of the wetting zone in a TSG and its effect on the effective distribution of liquid inside the barrel is an important aspect of solid-liquid mixing inside the granulator. This becomes especially crucial in the case of TSG because granulators from different manufacturers come with one or two liquid addition ports, which should influence the wetting zone. Varying the width of the wetting zone in the granulator while keeping the liquid flux constant allowed investigation of this effect.

Despite a drastic increase in the wetting zone width $W{Z}_{width}$ from 2$d_p$ to 8$d_p$, the pattern of increase in the liquid bridge fraction remained essentially the same (Figure 7). However, as the shearing time increases, in the case of a narrower wetting zone, the liquid from the bridge is given back to the particle surface, indicating a dominance of convective transport of liquid (Figure 7a). When the wetting zone was wide, there was a slight increase in ${\overline{Z}}_b$ indicating an increase in the volume of liquid in bridges. This also translated in $V_{bf}$ but the effect was delayed. While ${\overline{Z}}_b$ increased initially, $V_{bf}$ showed a very limited increase for a narrower wetting zone (Figure 7b). However, the ${\overline{Z}}_b$ increased when the $V_{bf}$ was already equilibrated for all the three different wetting zone widths. This increase was more pronounced for the larger wetting zone.

These results once again suggest that the liquid on the particle surface is quickly transferred to the liquid bridges. However, the liquid in bridges is later redistributed during the rearrangement of particles due to the shear to form more liquid bridges. The results also suggest that the liquid added to a broader area provides a clear benefit in terms of formation of more liquid bridges per particle, which is beneficial for granulation.

3.2. Solid-Liquid Mixing in the Mixing Zone of a TSG

To investigate the solid-liquid mixing in the mixing zone of a TSG, DEM simulations were performed using the set-up presented in Figure 3. About 46,000 particles were introduced, and the simulation was initialized for 0.5 sec after the wetting of particles (but not allowing any liquid transfer) exposed to a periodic interface in the z-direction to obtain representative randomness of particles (Figure 8a). After that, the amount of liquid on the particle surface and its transfer to liquid bridges by the shear mixing due to co-rotation of the two kneading discs was tracked in time. The shear mixing of particles resulted in a quick transfer of liquid from the wetted particle surface to neighboring dry particles in the barrel. This resulted in a rapid reduction in the liquid volume per particle over time from 0.15 to 2 $\times {10}^{-6}$ (see in Figure 8a–c and Figure 9a). The liquid on the particle surface was mostly transferred to the liquid bridges when the wetted powder was compressed between the closely inter-meshing kneading discs (see the right snapshots in Figure 8a, where particles with the maximum volume of liquid in bridge, i.e., particles connected with a higher number of liquid bridges are located in inter-meshing zone).

As granulation proceeds in TSG, both volumes of liquid on the surface of particles and the liquid in primary liquid bridges also reduced in this period (Figure 9b). This suggests that liquid transferred from particle surfaces was not contributing to agglomeration as a first step, but that conductive transport of liquid from the particle surface to liquid bridges between adjacent particles is the dominating first step to induce agglomeration. Hence, the primary liquid bridges were mostly present in over-wetted lumps, from which liquid was then convected under the shearing motion leading to the formation of newer liquid bridges. Further shear mixing by the kneading discs resulted in the equilibration of liquid on the particle surface and an increase in the number of liquid bridges. A simultaneous transfer of liquid from particle surface to liquid bridges led to oscillation in the mean liquid bridge volume. However, due to liquid mass conservation across particle surface and liquid brides, the liquid bridges contained a lower amount of liquid as time progressed and their number increased.

The results from this study suggest that the agglomeration is a delayed process, which is a follow-up event of equilibrated wetting of the particle. The convective mixing between solid-liquid followed by conductive transport of liquid from particle surface in the mixing zone leading to agglomeration confirms the importance of understanding the transport and mixing within a TSG. The results from this study are preliminary and require a more detailed study to understand the liquid distribution in other zones of the TSG, as well as the impact of screw geometry. Nevertheless, this modeling technique can be very useful for a detailed investigation of the mixing of solid particles and liquid. Since every particle in the simulation is tracked, the velocity and contact frequency information can be extracted to be used in refined PBE-based models for a more mechanistic simulation with local insight.

4. Conclusions

This study presented a theoretical analysis of the particle and liquid mixing based on particle-scale DEM simulations, in the context of HSWG using TSG. In the initial simulations, individual particles in a short quasi-two-dimensional simulation domain were tracked. Liquid bridge formation, rupture, as well as the effect of the bridges on inter-particular forces were simulated. Thereafter the same simulation approach was applied to the mixing zone of a TSG and the effect of process parameters on the particle and liquid mixing rate has been investigated. The results from this study demonstrate the potential of the particle-scale DEM method in modeling solid-liquid mixing and developing a detailed understanding of relevant phenomena. An extension of this study should be performed for detailed investigation of solid-liquid mixing along the complete length of a TSG. Since every particle in the simulation is tracked, the velocity and contact frequency information can be extracted to be used in PBE-based models for more mechanistic, still computationally realistic simulations. Moreover, there are multiple opportunities to increase the fidelity of the DEM simulations in the future, e.g., by considering viscous liquid bridge forces, or surface roughness effects.