Modeling of High-Density Compaction of Pharmaceutical Tablets Using Multi-Contact Discrete Element Method.

The purpose of this work is to simulate the powder compaction of pharmaceutical materials at the microscopic scale in order to better understand the interplay of mechanical forces between particles, and to predict their compression profiles by controlling the microstructure. For this task, the new framework of multi-contact discrete element method (MC-DEM) was applied. In contrast to the conventional discrete element method (DEM), MC-DEM interactions between multiple contacts on the same particle are now explicitly taken into account. A new adhesive elastic-plastic multi-contact model invoking neighboring contact interaction was introduced and implemented. The uniaxial compaction of two microcrystalline cellulose grades (Avicel® PH 200 (FMC BioPolymer, Philadelphia, PA, USA) and Pharmacel® 102 (DFE Pharma, Nörten-Hardenberg, Germany) subjected to high confining conditions was studied. The objectives of these simulations were: (1) to investigate the micromechanical behavior; (2) to predict the macroscopic behavior; and (3) to develop a methodology for the calibration of the model parameters needed for the MC-DEM simulations. A two-stage calibration strategy was followed: first, the model parameters were directly measured at the micro-scale (particle level) and second, a meso-scale calibration was established between MC-DEM parameters and compression profiles of the pharmaceutical powders. The new MC-DEM framework could capture the main compressibility characteristics of pharmaceutical materials and could successfully provide predictions on compression profiles at high relative densities.

1. Introduction

The ability to predict the bulk behavior of granular materials is of great importance for many industrial applications (i.e., tableting, metal forming) when deformation has to be handled in a controlled manner. Pharmaceutical powders are a branch of granular materials and while undergoing high loaded compression under confined conditions, allow the formation of compact granules, especially tablets, a process known as compaction or tableting. Pharmaceutical powder compaction is a crucial production process for pharmaceutical manufacturing due to the prevalence of tablets as solid dosage forms. Understanding the compaction behavior is of practical importance to improve the efficiency of product development and the manufacturing performance [1,2].

Powder compaction is frequently modeled using either continuous or discrete numerical techniques, or a combination of both. On the one hand, the finite element method (FEM) is a continuum approach that allows for the representation of deformation at a larger scale when combined with a suitable constitutive law, such as the Drucker–Prager cap (DPC) [3] or modified DPC [4], but does not reveal the physics of the system at the particle level. For this reason, one way to model the mechanical response of discrete elements at a particle level is the application of the multiple particle finite element method (MPFEM) [5,6,7] where each individual particle is being meshed with finite elements. Another way is to combine FEM and discrete element method (DEM), referred as meshed discrete element method (MDEM) [8,9], an approach that uses the contact detection algorithm from DEM and applies it in the FEM context. The capability to consider particles as deformable bodies is the key benefit of these approaches. The critical problem is their computational complexity, which prevents them from being used in large-scale industrial operations.

On the contrary, granular materials are mathematically described as a collection of particles by using the soft discrete element method (DEM) [10], and the bulk behavior of the granular materials is determined based on the interactions between pairs of particles. Newton’s second law of motion is used to define the trajectories of each particle, and particle deformation is proportional to the overlap between particles that are in contact. When an overlap is detected, an appropriate contact model is called, which associates the overlap with the force experienced at the contact point.

The DEM is the only currently available technique that can provide insight at the particle level and has been used in a wide variety of applications, such as ceramics [11], pharmaceutical, and food industries [12,13]. Given that the DEM is an excellent tool for studying these applications, modeling of confined powder compression with DEM remains a challenge. Based on Hertz’s theory, particles are considered to be rigid bodies to which deformation occurs locally and concentrated at the contact points. However, the basic assumptions of Hertz’s theory (classical DEM) is limited on capturing small deformations [14] and is rendered invalid for large strains, which typically occur under high relative densities and high loads. As a result, Hertz’s contact theory must be extended or modified to account for the fact that a larger (flattening) contact area results in a higher contact stiffness. To address this issue, many existing classical DEM contact models have added contact stiffness depending on plastic deformation and flattening in the contact areas. Ganrer et al. [15] proposed another adhesive elastic-plastic contact model to simulate powder compaction. The authors calibrated the contact stiffness between-particles (p-p) and particles-walls (p-w) using simulations of mono-sized particles, to reach high relative densities, at the macro-scale, and to predict the compression profiles. In a similar way, Y. Gao et al. [16]. applied the Luding’s elastioplastic [17] contact model to model powder compaction. Coarse-grained particles were used in this method and contact stiffness between particles (and p-w) was calibrated for one material under uni-axial compaction and the calibrated values were then used to predict the compression profiles of mixtures of additional materials.

In fact, many of the existing classical DEM contact models regard each contact of the same particle as independent of neighboring contacts, which is reasonable for loose powders but unrealistic for high relative densities. Three important characteristics of the mechanical behavior of elastic-plastic contacts are induced by a packing of simultaneously deforming particles undergoing uni-axial compaction (Figure 1): (a) the die filling and rearranging of the particles in the first phase; (b) the deformation is initially elastic, at this phase the contact areas between the particles are small with each contact independent of its neighbors; (c) the material yields entering the plastic zone and, because the pores are now almost closed the spatial confinement creates high contact pressures that allows for an additional degree of resistance that becomes more significant with time. Fischmeister and Arzt [18] refer to this as “geometrical hardening”. At this last phase, the compact is pressed to (nearly) full density [19], and mechanically behaves like a porous solid. The development of contact pressure originates from changes in the crystal lattices or intermolecular interactions. It has been suggested that contact dependency arises at relative densities of 0.7 and higher [14,20,21]. In practice, a relative density higher than 0.8 is required to produce commercial pharmaceutical tablets with adequate mechanical strength. As a result, multiple-contact DEM contact models are necessary.

Figure 1

The three phases of a packing of simultaneously deforming particles.

In this regard, researchers have made efforts on the formulation of the multi-contact discrete element method (MC-DEM) as an attempt to implicitly introduce particles’ deformability. One way is to enhance Hertz’s elastic contact theory and rewrite its classical equations such that not only local particle deformation but also the global deformation is taken into account. The global deformation is defined as the result of multiple contacts imposed on a single particle by neighboring particles. Brodu et al. [22] suggested a technique wherein the strain field acting on a single particle is coupled with the classical Hertz’s equation to account for the global deformation. Brodu et al. [22], validated their novel model by predicting the compression profiles of a packing of hydrogel balls compressed at low stresses. As an alternative, the stress field might be used in the equations. As a result, a multi-contact model that takes contact dependency into account was proposed by Frenning [23]. The particle global deformation was related to the isotropic stress tensor in this case. Giannis et al. [24] introduced a stress-based multi-contact model that takes anisotropic particle deformation into account and was validated for relevant materials with in the elastic regime. Attempts have also been made to explicitly introduce deformable particles in the framework of the DEM in order to address the basic assumptions of Hertz’s theory. A sophisticated model was presented by Rojek et al. [25] proposed the approach of the so-called deformable discrete element method (DDEM). This method is conceptually similar to the method provided by Brodu et al. [22]. The main difference is that particle deformability is introduced explicitly with the DDEM. The per particle stress tensor generates the isotropic particle deformation. As a result, the new deformable shape induces the formation of new contact points (not accessible with classical DEM). The main issues of this method are that (a) only isotropic particle deformation is considered, and (b) the high computing cost limits its use to 2D cases [25] or 3-D [26] cases with a small number of particles.

This work is divided into two parts: (a) using an extension of Giannis et al. [24] multi-contact model, to consider plastic deformation of pharmaceutical particles, and (b) to compare the results of the experiments with model calculations based on the MC-DEM framework.

The article’s outline is as follows: Section 2.1 contains the basic equations of the classical formulation of the DEM; multi-contact modeling is being discussed in Section 2.2 of this document; the materials used in this study are given in Section 3 of this article; Section 4 presents the calibration strategy and numerical results; and the last section, concludes with some final thoughts.

2. DEM’s Theoretical Background

Particle deformations are reproduced in soft-particle DEM by overlaps between interacting particles. When an overlap is detected, a contact law is used to compute the contact forces (force–displacement) between two particles. The underlying assumption is that particle contacts are independent of one another, and therefore contact forces are resolved locally. Newton’s equations of motion are used in this method to determine the connection between particle motion and forces acting on each particle. The equations of a particle’s translational and rotational motion are: where , , and are the mass, acceleration, moment of inertia and angular velocity for particle , respectively; , , are the normal force, tangential force, and torque acting on particles i and j at contact points, respectively; is the acceleration due to gravity.

Different contact models can be used to express force–displacement laws at contact points. This study does not go into great depth on the various contact models and their related equations. Rojek [27] and Thakur [28] summarize the many contact models that are used in discrete particle simulations. In their works, O’Sullivan [29] and Thornton [30] go into deep details on contact models.

2.1. Classical Hertz–Mindlin Contact Model

The linear spring-dashpot model [31], in which the spring stiffness is considered to be constant, is the simplest contact configuration. To enhance the linear contact model, the Hertz theory [32] (classical) is used to calculate the force-displacement relation for contacting particles (e.g., nonlinear spring-dashpot model). In this case, the normal stiffness varies depending on the degree of overlap. The Hertz–Mindlin [33,34,35,36] contact model is another contact model for representing the force-displacement relation. This nonlinear model combines accuracy and simplicity of the Hertz theory in the normal direction with the Mindlin model in the tangential direction. This model includes a contact force as well as a viscous contact damping force at contact points. In the normal (n) and tangential (t) axes, these forces were calculated using elastic springs and dashpots (Figure 2). The normal repulsive contact force is: where is the normal stiffness coefficient, with the effective radius and is the effective Young’s modulus. In this expression, and represent the particles Poisson’s ratio and shear modulus, respectively. The normal overlap is , is the relative velocity in normal direction of interacting particles and the viscoelastic damping constant for normal contact viscosity. The tangential force is [24]: where is the tangential stiffness coefficient and the effective shear modulus. The tangential overlap is , is the relative velocity in tangential direction of interacting particles and the viscoelastic damping constant for tangential contact viscosity.

Figure 2

Contact force model illustrating particle interaction with normal and tangential forces.

The tangential overlap, , between particles obtained by integrating surface relative tangential velocities during elastic deformation of the contact is given as [37,38]: where is the velocity component tangential to the contact surface and is the time-step.

The tangential and normal forces are connected by Coulomb’s law, , in the event of sliding, there is dynamic friction. . The dynamic and static friction coefficients are assumed to be equal in this case., . In order to allow for a restoring force, a static situation requires the use of an elastic spring, i.e., a non-zero remaining tangential force in static equilibrium due to activated Coulomb friction. By applying a torque to the contacting surfaces, rolling friction can be controlled. The rolling friction constant directional torque (CDT), , used for this study, is given by [39]:

The particles in these models are assumed to be spherical and do not deform during simulation. In a strict sense, it is assumed that particles are undergoing some kind of pseudo deformation, and this model is known as the truncated Hertz–Mindlin model [40]. Moreover, these models include binary interactions between two particles, which implies that during contacts, particles are in touch through a single point.

2.2. Multi-Contact Adhesive Elastic-Plastic Model

The linear or non-linear elasticity theory, on the other hand, is only applicable to small deformations. However, when it comes to powder compaction, plasticity prevails (flattening in contact areas), necessitating the modeling of elastic-plastic spheres in contact. As a result, Hertz theory must be extended to cases in which particles are deformed plastically. Persson and Frenning [40], for instance, presented an extension of classical Hertz theory to account for elastic-plastic contacts. In this example, a limiting contact pressure was added, whereas adhesion was not added, such that plastic deformation begins after the contact region’s maximum pressure is achieved.

We propose a novel adhesive elastic-plastic multi-contact model that combines concepts from the Luding [17] and Edinburgh [28] adhesive elastic-plastic (hysteretic) contact models and the multi-contact model proposed by Giannis et al. [24], the pseudo-code of the algorithm utilized is briefly described in the Appendix A. For the first time, this model is being used in this study to investigate the behavior of elastic-plastic medicinal materials. When two particles collide with one another elastic and plastic deformation (linear and non-linear force-displacement curve) occur. A nonlinear contact model is proposed that takes into consideration both elastic-plastic contact deformation and adhesion. The adhesive plastic force is:

Figure 3 illustrates the force-displacement curve. The term “displacement” refers to the overlap of particles. The loading, unloading, re-loading, and adhesive branches e is defined by the loading branch stiffness , the loading-unloading branch stiffness , the adhesion branch stiffness , the plastic overlap (deformation) and the constant pull-off force . During initial loading the contact model follows the virgin loading path , until the maximum overlap is reached at . The maximum overlap is a contact-specific history-dependent parameter that is updated and saved. During unloading the contact will alter from virgin loading to unloading/reloading , which depends on . At , the force is decreasing from its value to zero at overlap , which resembles the plastic contact deformation (remaining overlap). The plastic overlap is defined as: for cases where the limit is met, with results in (no remaining overlap) yields to a special case of non-linear elasticity. Hence, the non-linear elastic Hertz–Mindlin contact model is included as a special case. On the other hand, captures a perfectly plastic contact. Unloading below results in attractive adhesion forces until the minimum force is equal:

Figure 3

A non-linear hysteretic, adhesive force-displacement () relation in normal direction. The slope of the unloading and reloading branch interpolates between and a maximum stiffness .

And the overlap is:

Attractive forces emerge as the unloading process continues.

In order to account for the fact that a larger contact surface leads to a higher contact stiffness, the coefficient is made dependent on the maximum overlap (history dependent parameter):

The behavior of the unloading slope described by Equation (11) is similar to that assumed by Luding [10,17], with the exception that nonlinear behavior is addressed here. Likewise, the limit of plastic flow overlap is given: where is the dimensionless plasticity depth, defined in relation to the reduced radius. The original contact model proposed by Luding can be characterized as a piecewise linear hysteretic model [17]. For the virgin loading, the contact normal stiffness and normal overlap is used to calculate the force. While stiffness is not a physical parameter according to Luding’s contact model, in this work it is depending on the Young’s modulus. The new stiffness is is identical to the one of Hertz theory (Section 2.1) and similar to the one of Edinburgh contact model. Furthermore, in contrast with Luding’s and in line with the Edinburgh [28] contact model, a non-linear force-displacement relation is proposed . By contrast, the unloading stiffness is load dependent as in Luding’s model. Additionally, the new contact model was supplied with a non-linear adhesion .

Moreover, to address the fundamental assumption of the classical DEM, which treats each contact locally as a binary pair interaction, Giannis et al. [24] presented a nonlocal model which takes into account the mutual influence between contacts. While this model has been verified for cases in the elastic regime, in this study we will extend its applicability to capture plasticity.

The main idea of the on how to account for multi contact effect is shown in Figure 4. Multiple contacts acting on a particle have been taken into account by using the trace of the average stress tensor coupled with the Poisson’s ratio (), the contact area () between interacting particles, and a material-dependent prefactor (). More information may be found here [24]. The new multi-contact law formulation yields to this equation:

Figure 4

A multi-contact modification of the classical DEM.

The first term in Equation (13) was defined before in Section 2.1 of this article (Equation (2)). The second term carries the information from neighboring particles operating on the particle. In this expression, is a dimensionless empirical prefactor that allows for particle geometry changes to be taken into consideration indirectly, is the Poisson’s ratio, and is the contact are between the interacting particles. The isotropic component of the stress is the pressure , with and the stress tensors of particle and , respectively. By combining Equations (6) and (13), the multi-contact model can be extended from linear to plastic deformation, yielding a new equation:

In the next sections, this new equation (Equation (14)) will be investigated and verified for the modeling cases of uni-axial compaction of pharmaceutical materials.

3. Materials and Methods

3.1. Materials

Two microcrystalline cellulose grades Avicel® PH 200 (FMC BioPolymer, Philadelphia, PA, USA) and Pharmacel® 102 (DFE Pharma, Nörten-Hardenberg, Germany) were studied in depth. Henceforth, the powders shall be referred to by the abbreviations MCC-A and MCC-P. A certain proportion of the MCC particles, particularly the larger particle sizes, are rounded agglomerates. Taking this into account, and due to the simplicity of this approach in simulation, spherical shapes are used in DEM simulations. The powder characteristics particle size distribution (PSD) and true density, are given in Table 1 and are available in the literature [41].

Table 1

Powder characteristics: PSD and densities [41].

Material	x10 (Q3) (μm)	x50 (Q3) (μm)	x90 (Q3) (μm)	Span (-)	True Density (kg m−3)
MCC-A	82.9	224.6	379.3	1.32	1541.1
MCC-P	28.3	86.5	173.8	1.68	1533.7

3.2. Experimental Methods

Compaction experiments were performed applying a Styl’One evolution compaction simulator (CS; Medel’Pharm, Beynost, France). This equipment can accurately control the compaction process and allows for in-depth investigation of powder properties and to extract force/displacement profiles. In-die data were evaluated by applying the software ANALIS (Medel’Pharm, Beynost, France). Generic profile was applied for compaction to reach compression stresses of approx. 30, and 180 MPa.

3.3. Numerical Methods

In a wide range of applications, including this study, the discrete element method (DEM) is used to model and analyze granular materials. However, predictions can only be correct if the input parameter values are carefully chosen. There are a number of input parameters that need to be tuned depending on the contact model used; this procedure is known as calibration. An in-depth review into calibration is provided by [42,43,44,45]. If Luding’s [17] original contact model is to be used, a total of 19 input parameters must be predefined or calibrated. A comprehensive experimental determination of these parameters would be highly time-consuming and labor-intensive to accomplish. Fortunately, not all parameters have the same effect on the simulation output. As a consequence, only the parameters that are most significant to the validity of the simulation results are considered. Material parameters for the single particles and used in this study were obtained from the study of Cabiscol et al. [41]. The calibration technique consists of bibliographical sources, experiments, and their replication by DEM simulations. DEM simulations of nano-indentation experiments were used to determine the fitting parameter that expresses the identical experimental results, in terms of force and displacement; therefore the Young’s modulus (E) of a single particle may be determined. The ring shear tester and the Jenike wall test were used for a direct determination of the sliding friction between particles (μs(pp)) and between particles and walls (μs(pw)). Later, the tumbling drum test was used to determine the final frictional and rotational DEM related parameters μs(pw), μr(pp), and μr(pw). This test functioned as a calibration test for the rolling friction as well as a second and final iteration for μs(pw), starting from the values obtained from the ring shear tester and the Jenike wall test. The complete calibration method is described in depth in the Cabiscol et al. study [41]. The DEM input parameters for the material properties are summarized in Table 1. The calibration of the input parameters of the new multi-contact model will be discussed in Section 4.1.1.

4. Results and Discussion

4.1. Determination of a Representative Volume Element (RVE)

Full-scale DEM simulations need a significant amount of computational power; to bypass this limitation, a representative volume element is used. The literature has several definitions of the representative volume element (RVE), the most notable of which are discussed in [46,47,48,49]. Although there is no single and exact definition of the RVE, the basic concept is that the RVE should be large enough to retain the microstructural information while being small enough in relation to the macroscopic structural dimensions to eliminate fluctuations. In this study, the RVE concept was used to speed up the simulations. The method used to determine an RVE was similar to that described by Wiącek et al. [50]. The basic ideas are as follows: (a) set an initial packing contained in a small domain (the smallest possible); (b) gradually expand the domain’s dimensions while maintaining the same particle size distribution (PSD) and packing density; (c) carry out numerical simulations to obtain the force-displacement curve; and (d) analyze the results to determine if they are converging.

The simulations presented here involve a series of uni-axial compaction tests in cubes of the following ranging sizes: (0.6 mm)3 (1stRVE), (0.8 mm)3 (2nd RVE), (1.0 mm)3 (3nd RVE), and (1.4 mm)3 (4th RVE) (see Figure 5). To eliminate the wall boundary effect, periodic boundaries were used along the X and Y axes. The cubes contain a top and a bottom plate. A series of uniaxial compression simulations were performed using our new multi-contact DEM model. After calibration, as will be discussed later in Section 4.1.1, and based on the data shown in Figure 6 we can confirm the existence of an RVE since the results are converging. However, due to the small size of the sample, the results of the first RVE underestimated its macroscopic stress–strain response. There is excellent agreement between the results of the following RVEs, as well as with the experimental data. It is therefore decided to use the second RVE for practical reasons (less computational time) as follows.

Figure 5

Determination of a representative volume element (RVE).

Figure 6

Converging analysis of the suggested RVEs.

4.1.1. Calibration Method for the Input Parameters of the Multi-Contact Model

In this section, the calibration for the new multi-contact model is shown. In this case, and as discussed in Section 4.1, the system under consideration is a cube with dimensions (0.8 mm)3 (2nd RVE in Figure 5) along x-y-z directions. The system under consideration contains 698 particles for the MCC-A material and 1193 particles for the MCC-P material with a particle size distribution for both materials given in Section 3.1 and Table 1. The particles are initially randomly positioned in a cubic system with periodic boundary constraints in order to minimize wall effects. After initial deposition, the particles are allowed to grow. Growth is terminated as soon as the desired packing density of approximately 59% is reached and is in line with experimental results reported in the literature [41,51]. The sample was then compressed uni-axially along the z-axis to a maximum target strain of 57% for MCC-A and 53% for MCC-P, then it was decompressed. A strain-driven simulation was used to achieve the maximum desired stress of 29 MPa for the MCC-A material and 25 MPa for the MCC-P material. The calibrated material parameters presented in Table 2 (Section 3.2) were used here. However, the input parameters for the multi-contact model were obtained using an iterative process to determine the optimum parameters that better fits the experimental results. A parameter optimization method was used based on a series of simulations, similar to the one presented by Gao et al. [16]. Given the experimental macroscopic stress and strain response the 2 value between the experimental and simulated data was calculated from: where indicates the stress response of the experimental data, indicates the stress response of the simulated data and indicates the mean of the stress response of the experimental data. The value of 2 was used to evaluate the accuracy with which the simulated input parameters fit the experimental data; a successful fit was attained when 2 was close to 1. Therefore, when 2 exceeded 0.95, the iterations needed for calibration were terminated. Figure 7 shows that experimental and simulated results are in excellent agreement, indicating that the calibration was successful. Table 3 summarizes the input parameters that were calibrated. For the dimensionless plasticity depth , a high and constant value was selected (low contact stiffness) to achieve a high contact stiffness; when necessary, the prefactor β was tuned accordingly.

Table 2

The input parameters for single particles and walls [41].

Property	Symbol	Units	MCC-A	MCC-P
Young’s modulus—particle (p)	E	Nm−2	2.58 × 108	1.34 × 109
Young’s modulus—wall (w)	E	Nm−2	7.62 × 1010	7.62 × 1010
Poisson’s ratio—particle	ν	-	0.30	0.30
Poisson’s ratio—wall	ν	-	0.31	0.31
Coefficient of restitution particle	COR(p-p)	-	0.352	0.346
Coefficient of restitutio—wall	COR(p-w)	-	0.352	0.346
Coefficient of sliding fric—(p-p)	μs(pp)	-	0.561	0.548
Coefficient of sliding f—(p-w)	μs(pw)	-	0.707	0.715
Coefficient of rollin—(p-p)	μr(pp)	-	0.3	0.3
Coefficient of rol—(p-w)	μr(pp)	-	0.01	0.01
Density	ρ	kg/m3	1541.1	1533.7

Figure 7

Calibration for the: (a) MCC-A material under uni-axial compaction at maximum target stress of 29 MPa; (b) MCC-P material under uni-axial compaction at maximum target stress of 25 MPa.

Table 3

Multi-contact model input parameters.

Property	Symbol	Units	MCC-A	MCC-P
Unloading stiffness	k2/k1	-	120	120
Adhesion stiffness ratio	Kc/k1	-	0.5	0.5
Dimensionless plasticity depth	φf	-	0.99	0.99
Prefactor of the MC-dem	β	-	1.3	1.5

4.1.2. Verification for Uni-Axial Compaction for MCC-A

In this section, the simulation results for compaction of the MCC-A material are shown. The system is identical to the one presented in Section 4.1.1. The sample was compressed uni-axially along the z-direction to a maximum target strain of 71%, and then decompressed. The target stress for this case is 180 MPa. The calibrated material input parameters given in Section 3.2 (Table 2) and, for the multi-contact DEM model given in Section 4.1.1 (Table 3), were used. When the calibrated prefactor β = 1.3 (Section 4.1.1 (Table 3)) was used, an excellent agreement between experimental and simulated results was achieved, as shown in Figure 8b. A value of prefactor β = 0.0 indicates that the multi-contact effect is not included and when β = 0.0 conventional DEM underestimates the macroscopic stress–strain response, as shown in Figure 8a. It is also clear from comparing Figure 8a,b the multi-contact effect predominates for strains higher than 0.2.

Figure 8

Verification for the MCC-A material under uni-axial compaction at maximum target stress of 180 MPa: (a) without the multi-contact effect (prefactor β = 0.0); (b) with the multi-contact effect (prefactor β = 1.3).

4.1.3. Verification for Uni-Axial Compaction for MCC-P

In this section, the simulation results for compaction of the MCC-P material are shown. The system is identical to the one presented in Section 4.1.1 for the MCC-P material. The sample was compressed uni-axially along the z-direction to a maximum target strain of 69%, and then decompressed. The target stress for this case is 185 MPa. The calibrated material input parameters given in Section 3.2 (Table 2) and, for the multi-contact DEM model given in Section 4.1.1 (Table 3), were used. When the calibrated prefactor β = 1.5 (Section 4.1.1 (Table 3)) was used, an excellent agreement between experimental and simulated results was achieved, as shown in Figure 9b. As expected with β = 0.0 conventional DEM underestimates, the macroscopic stress–strain response is shown in Figure 9a. It is also clear from comparing Figure 9a,b that the multi-contact effect predominates for strains higher than 0.2, around the same point as that seen for a the MCC-A material in Section 4.1.2.

Figure 9

Verification for the MCC-P material under uni-axial compaction at maximum target stress of 185 MPa: (a) without the multi-contact effect (prefactor β = 0.0); (b) with the multi-contact effect (prefactor β = 1.3).

5. Conclusions

In this study, by employing our new elastic-plastic multi-contact DEM model, the compaction profiles (stress–strain) of microcrystalline cellulose grades Avicel® PH 200 (FMC BioPolymer, Philadelphia, PA, USA) and Pharmacel® 102 (DFE Pharma, Nörten-Hardenberg, Germany) were successfully predicted. It was also shown that the multi-contact effect predominates for strains higher than 0.2. A calibration strategy to calibrate the input parameters, prefactor β, for the multi-contact model was presented here. The prefactor β was calibrated at low relative densities (low macroscopic stress) and subsequently used for high relative densities (high macroscopic stress). The new multi-contact model requires a separate calibration for each material as prefactor β is a material-dependent parameter. However, more research is needed to determine if this is also true for a mixture of other relevant materials.

In terms of the input parameters for the multi-contact model prefactor β, the unloading stiffness and the cohesion stiffness were the only parameters that were calibrated. The loading stiffness was related to the Young’s modulus of the material. The material input parameters were calibrated separately. The concept of a representative volume element (RVE) was used to speed up simulations. In comparison to alternative approaches that use coarse-grained particles, the RVE was preferred because the particle size distribution (PSD) can be maintained while using the RVE.

In future work, we will aim to calibrate the prefactor β as an intrinsic material parameter. The aim is to conduct uni-axial compaction simulation in a series of relevant materials, then create a comprehensive database, and finally, with the assistance of artificial intelligence (e.g., neural network), generalize the results. Furthermore, the results presented here are based on the assumption of perfect spherical particles; in a future attempt, the real shape of the particles should be addressed.