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.


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 In this regard, researchers have made efforts on the formulation of th discrete element method (MC-DEM) as an attempt to implicitly introduc formability. One way is to enhance Hertz's elastic contact theory and rew equations such that not only local particle deformation but also the global taken into account. The global deformation is defined as the result of m imposed on a single particle by neighboring particles. Brodu et al. [22] su nique wherein the strain field acting on a single particle is coupled wi Hertz's equation to account for the global deformation. Brodu et al. [22], novel model by predicting the compression profiles of a packing of hydr pressed at low stresses. As an alternative, the stress field might be used in As a result, a multi-contact model that takes contact dependency into ac posed by Frenning [23]. The particle global deformation was related to the tensor in this case. Giannis et al. [24] introduced a stress-based multi-con takes anisotropic particle deformation into account and was validated for rials with in the elastic regime. Attempts have also been made to explicitl formable particles in the framework of the DEM in order to address the bas of Hertz's theory. A sophisticated model was presented by Rojek et al. [25 approach of the so-called deformable discrete element method (DDEM). conceptually similar to the method provided by Brodu et al. [22]. The ma that particle deformability is introduced explicitly with the DDEM. The pe tensor generates the isotropic particle deformation. As a result, the new de induces the formation of new contact points (not accessible with classical D issues of this method are that (a) only isotropic particle deformation is cons the high computing cost limits its use to 2D cases [25] or 3-D [26] cases wit ber of particles.
This work is divided into two parts: (a) using an extension of Giannis contact model, to consider plastic deformation of pharmaceutical particles, pare the results of the experiments with model calculations based on the M work.
The article's outline is as follows: Section 2.1 contains the basic equat sical formulation of the DEM; multi-contact modeling is being discussed i 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.

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 m i , .. a i , I i and . ω i are the mass, acceleration, moment of inertia and angular velocity for particle i, respectively; F nij , F tij , τ ij are the normal force, tangential force, and torque acting on particles i and j at contact points, respectively; g 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.

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 k n = 4 3 E * √ R * is the normal stiffness coefficient, with R * = R i R j R i +R j the effective radius and E * = is the effective Young's modulus. In this expression, ν and G represent the particles Poisson's ratio and shear modulus, respectively. The normal overlap is δ n , . δ n is the relative velocity in normal direction of interacting particles and γ n the viscoelastic damping constant for normal contact viscosity. The tangential force is [24]: where k t = 8G * √ R * is the tangential stiffness coefficient and the effective shear modulus. The tangential overlap is δ t , . δ t is the relative velocity in tangential direction of interacting particles and γ t the viscoelastic damping constant for tangential contact viscosity. force, a static situation requires the use of an elastic spring, i.e., a non-zero remaining tan gential force in static equilibrium due to activated Coulomb friction. By applying a torqu to the contacting surfaces, rolling friction can be controlled. The rolling friction constan 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 durin simulation. In a strict sense, it is assumed that particles are undergoing some kind o pseudo deformation, and this model is known as the truncated Hertz-Mindlin model [40 Moreover, these models include binary interactions between two particles, which implie that during contacts, particles are in touch through a single point.

Multi-Contact Adhesive Elastic-Plastic Model
The linear or non-linear elasticity theory, on the other hand, is only applicable t small deformations. However, when it comes to powder compaction, plasticity prevai (flattening in contact areas), necessitating the modeling of elastic-plastic spheres in con tact. As a result, Hertz theory must be extended to cases in which particles are deforme plastically. Persson and Frenning [40], for instance, presented an extension of classica Hertz theory to account for elastic-plastic contacts. In this example, a limiting contact pres sure was added, whereas adhesion was not added, such that plastic deformation begin after the contact region's maximum pressure is achieved.
We propose a novel adhesive elastic-plastic multi-contact model that combines con cepts from the Luding [17] and Edinburgh [28] adhesive elastic-plastic (hysteretic) contac models and the multi-contact model proposed by Giannis et al. [24], the pseudo-code o the algorithm utilized is briefly described in the Appendix A. For the first time, this mode is being used in this study to investigate the behavior of elastic-plastic medicinal mater als. When two particles collide with one another elastic and plastic deformation (linea and non-linear force-displacement curve) occur. A nonlinear contact model is propose that takes into consideration both elastic-plastic contact deformation and adhesion. Th adhesive plastic force is: The tangential overlap, δ t , between particles obtained by integrating surface relative tangential velocities during elastic deformation of the contact is given as [37,38]: where v t is the velocity component tangential to the contact surface and ∆t is the time-step. The tangential and normal forces are connected by Coulomb's law, F t ≤ µF n , in the event of sliding, there is dynamic friction. F t = µF n . The dynamic and static friction coefficients are assumed to be equal in this case., µ = µ d = µ s . 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), τ ij , 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.

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 materi-Pharmaceutics 2021, 13, 2194 6 of 16 als. 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 k 1 , the loading-unloading branch stiffness k 2 , the adhesion branch stiffness k c , the plastic overlap (deformation) δ 0 and the constant pull-off force F 0 . During initial loading the contact model follows the virgin loading path k 1 , until the maximum overlap is reached at δ max . The maximum overlap δ max is a contact-specific history-dependent parameter that is updated and saved. During unloading the contact will alter from virgin loading k 1 to unloading/reloading k 2 , which depends on δ max . At δ max , the force is decreasing from its value to zero at overlap δ 0 , which resembles the plastic contact deformation (remaining overlap). The plastic overlap is defined as: for cases where the limit is met, with k 1 = k 2 results in δ = 0 (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, k 2 → ∞ captures a perfectly plastic contact. Unloading below δ results in attractive adhesion forces until the minimum force is equal: Pharmaceutics 2021, 13, 2194 7 of 1 Figure 3. A non-linear hysteretic, adhesive force-displacement ( ) relation in normal direction. Th slope * of the unloading and reloading branch interpolates between and a maximum stiffnes .
Moreover, to address the fundamental assumption of the classical DEM, which treat each contact locally as a binary pair interaction, Giannis et al. [24] presented a nonloca model which takes into account the mutual influence between contacts. While this mode has been verified for cases in the elastic regime, in this study we will extend its applica bility to capture plasticity. And the overlap δ min 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 k 2 is made dependent on the maximum overlap δ max (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 ϕ f 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 k 1 and normal overlap δ n is used to calculate the force. While stiffness k 1 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 k n = 4 3 E * √ R * 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 F n = k 1 δ 3/2 n . By contrast, the unloading stiffness k 2 is load dependent as in Luding's model. Additionally, the new contact model was supplied with a non-linear adhesion F n = −k c δ 3/2 n . 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 (v), the contact area (A) 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:

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

Experimental Methods
Compaction experiments were performed applying a Styl'One evolut simulator (CS; Medel'Pharm, Beynost, France). This equipment can accura compaction process and allows for in-depth investigation of powder pr extract force/displacement profiles. In-die data were evaluated by applyi ANALIS (Medel'Pharm, Beynost, France). Generic profile was applied for reach compression stresses of approx. 30, and 180 MPa.

Numerical Methods
In a wide range of applications, including this study, the discrete e (DEM) is used to model and analyze granular materials. However, predi be correct if the input parameter values are carefully chosen. There are a n parameters that need to be tuned depending on the contact model used; th known as calibration. An in-depth review into calibration is provided by 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 A ij is the contact are between the interacting particles. The isotropic component of the stress is the pressure P ij = 1 3 tr(σ i ) + tr σ j , with tr(σ) = σ xx + σ yy + σ zz and σ i , σ j the stress tensors of particle i and j, 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.

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

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 Pharmaceutics 2021, 13, 2194 9 of 16 ANALIS (Medel'Pharm, Beynost, France). Generic profile was applied for compaction to reach compression stresses of approx. 30, and 180 MPa.

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.

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.
Pharmaceutics 2021, 13, 2194  In this section, the calibration for the new multi-contact model is show and as discussed in Section 4.1, the system under consideration is a cube w (0.8 mm) 3 (2nd RVE in Figure 5) along x-y-z directions. The system under   In this section, the calibration for the new multi-contact model is show and as discussed in Section 4.1, the system under consideration is a cube wi (0.8 mm) 3 (2nd RVE in Figure 5) along x-y-z directions. The system under contains 698 particles for the MCC-A material and 1193 particles for the M Figure 6. Converging analysis of the suggested RVEs.

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 R 2 value between the experimental and simulated data was calculated from: where y indicates the stress response of the experimental data,ŷ indicates the stress response of the simulated data and y indicates the mean of the stress response of the experimental data. The value of R 2 was used to evaluate the accuracy with which the simulated input parameters fit the experimental data; a successful fit was attained when R 2 was close to 1. Therefore, when R 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 ϕ f , a high and constant value was selected (low contact stiffness) to achieve a high contact stiffness; when necessary, the prefactor β was tuned accordingly.

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

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.

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

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

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 multicontact 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 pref-actor β 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 κ 2 and the cohesion stiffness κ c were the only parameters that were calibrated. The loading stiffness κ 1 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.