A multivariate regression parametric study on DEM input parameters of free-flowing and cohesive powders with experimental data-based validation

One of the key challenges in the implementation of discrete element method (DEM) to model powder’s flow is the appropriate selection of material parameters, where empirical approaches are mostly applied. The aim of this study is to develop an alternative systematic numerical approach that can efficiently and accurately predict the influence of different DEM parameters on various sought macroscopic responses, where, accordingly, model validation based on experimental data is applied. Therefore, design of experiment and multivariate regression analysis, using an optimized quadratic D-optimal design model and new analysis tools, i.e., adjusted response and Pareto graphs, are applied. A special focus is laid on the impact of six DEM microscopic input parameters (i.e., coefficients of static and rolling friction, coefficient of restitution, particle size, Young’s modulus and cohesion energy density) on five macroscopic output responses (i.e., angle of repose, porosity, mass flow rate, translational kinetic energy and computation time) using angle of repose tests applied to free-flowing and cohesive powders. The underlying analyses and tests show, for instance, the substantial impact of the rolling friction coefficient and the minor role of the static friction coefficient or the particle size on the angle of repose in cohesive powders. In addition, in both powders, the porosity parameter is highly influenced by the static and rolling friction coefficients.


Introduction
One of the big challenges in the pharmaceutical industry is to optimize the powder solid handling process performance. This significantly relies on the flowability of the granular material in different process stages including hopper discharge, powder feeding, blending, mixing and die filling. Flowability, which describes the behavior of powder, is dependent on a range of fundamental powder properties. The flow properties of powders can be affected by various essential particle properties, including mean particle size and distribution, particle shape, surface roughness and moisture content. Powder flow also depends on the interparticle interactions where, for instance, flowability decreases with increasing cohesiveness of the powder. According to Parker et al. [1], cohesion is defined as the bonding or joining of two particles of the same material. A bulk solid with poor flowability is considered to be cohesive, which is due to the relatively high interparticle forces compared to the particle weight [2]. The inherent strength of the interparticle forces of attraction (often called cohesive forces) is derived from a combination of van der Waals forces and electrostatic charge on the surface of the particles. These forces hold particles close to their neighbors and thus result in the formation of agglomerates. Capillary force (liquid bridge), van der Waals force, electrostatic force and the weight force are all considered as adhesive forces, mainly interparticle forces between the particles and the wall. Understanding the physical meaning of these properties and their impact on the final product design is crucial in the development and successful production of solid dosage forms.
During their operation, industrial machines dealing with powders might encounter diverse problems which, due to their complexity, are most of the time difficult to analyze by traditional means. Therefore, many studies are extensively applied to optimize their operations, to deeply understand their functionalities, and to predict possible defects. For modeling a granular material, several numerical approaches on different scales and with different degrees of complexity and computational costs can be found in the literature. The microscopic discrete element method (DEM) (Cundall and Strack [3]) accounts for the most popular schemes and provides an effective approach for the understanding of granular material's mechanical behavior and flowability on the microscale. In particular, DEM is a particle-based method used to model movement and interaction of bulk material [2]. Each particle is treated individually and its motion is described by Newton's equation of motion. A soft-sphere approach is used, where colliding particles are allowed to slightly overlap resulting in a repulsive force [4]. DEM has been proven to be an important and efficient tool in several fields such as powder metallurgy, chemical engineering, pharmaceuticals and agriculture. To mention some, DEM simulations have been applied in granular-related processes such as die filling [5], bulk compression [6], powder mixing [7], particles packing [8,9], agglomerates fragmentation [10], flow in screw conveyor [11][12][13] and outflow from a hopper. Apart from DEM, other approaches, such as continuum porous media mechanics, can be applied to compute granular material mechanical responses on larger (macroscopic) scales and capture important behaviors such as wave propagation or fracture [14][15][16]. However, these approaches assume the continuity of the matter and cannot explicitly capture the interactions between the particles. Therefore, the focus in this work would be merely on DEM approach.
One of the most challenging aspects in applying efficiently the DEM numerical approach is choosing the appropriate contact models and then finding the right values of the input particle parameters involved in any used model. As particles introduced in granular materials are considerably small in size, measuring the properties of every single particle is considered very costly, challenging and timeconsuming. Yet, some trials using atomic force microscopy have been given to measuring the individual particles' yield stress [17], Young's modulus (Y) [18] and interfacial energy [19]. The state-of-the-art approaches for calibration of DEM parameters in most cases are based on the trialand-error, where the track of the calibration process is lost in many cases by the human calibrators, especially when many factors are involved [20]. It is often the case that the DEM model is calibrated against bulk measurements and the input parameters are adjusted until the outputs of the model match with the experimental observations [21]. Therefore, it is of great interest for researchers to optimize the process of calibration. In the literature, different tools for DEM parameters calibration, some of which are (semi-) automatized, can be found. Al-Hashemi et al. [22] presented an extensive review related to the angle of repose (AoR), methods of measurements, appropriate applications and the influencing factors of different parameters. Coetzee [23] discussed different DEM calibration approaches covered in the literature over the last 25 years to assist future researchers to improve on the existing ones. Liu et al. [24] investigated the impact of particle aspect ratio on the DEM simulation output in studying powder flow from hoppers. The authors developed a modified description of the classic Beverloo equation. Zhou et al. [25] studied the influence of particle-particle and particle-wall coefficients of friction, i.e., static and sliding on the AoR of a sandpile test. They constructed a power-law relationship between measured AoR and the input parameters and showed that all friction coefficients followed a positive dependency, whereas particle size followed a significant negative dependency. Yan et al. [26] utilized a multilevel statistical analysis approach on hopper flow in DEM to study the influence of coefficients of friction and restitution, as well as Young's modulus. Using principal component analysis (PCA), the empirical power-law approach was found to agree with a semiempirical, which compared AoR to contact model parameter sensitivity matrices. The results in all cases showed that the coefficient of static friction was the critical parameter, which allows the gap to be bridged between bulk and interparticle interactions. In addition, it is shown that reducing the shear modulus (i.e., from 10 7 and 10 11 Pa) proved to have an insignificant effect on the powder flowability using free-flowing material which was also proven by Lommen et al. [27].
Boukouvala et al. [28] used a reduced-order approach to assess the importance of DEM parameters where PCA was coupled with surrogate models mapping which identified key areas of sensitivity. Major principal components enabled a reduced and computationally more efficient model to be engaged for future simulations. El Kassem et al. [29] calibrated the AoR of a pharmaceutical powder by varying coefficients of rolling friction and restitution using DoE, while determining the coefficient of static friction using the FT4 powder rheometer. Souihi et al. [30] used a statistical design of experiments (DoE) to provide a more efficient approach in investigating roller compaction. Orthogonal partial least squares (PLS) regression was used to analyze the results of a reduced central composite face-centered design in DoE. Wilkinson et al. [31] added that such approach using regression analysis has not been used so far for DEM simulations, but is highly recommended for developing more universal models for particle flow in a more efficient way. He proved that the Y has a negative impact on flow energy under specific conditions using a cohesive material. In this work, the effect of Y with the combination of varying the particle size in the same design model is studied which was not covered by [26] and [31]. To add up, Johnstone [32] used DoE to perform different calibration and validation methods based on experimental measurements. Benvenuti et al. [33] used artificial neural network to identify DEM simulation parameters. However, Rackl et al. [34] used a methodical calibration approach based on Latin hypercube sampling and Kriging, where they showed that different parameter combinations might lead to similar results. Wei et al. [35] studied the effect of three particle shapes on the AoR and bottom porosity distribution of a heap formed by natural piling using DEM. They found that the static friction has the major impact on the AoR and the porosity followed with the particle shape. Moreover, Li et al. [36] investigated the influence of three factors, namely moisture content, rice velocity and depth of the rice layer on the porosity of the flowing rice layer based on mass conservation and the poroelasticity theory. These results showed that the latter three factors have a significant influence on the porosity. Cheng et al. [37] developed a DEM calibration approach of granular soils based on the sequential quasi-Monte Carlo (SQMC) filter. They calibrated the micromechanical parameters of the contact laws against the stress-strain behavior of Toyoura sand in drained triaxial compression conditions at different confining pressures. Do et al. [38] developed two automated DEM calibration approaches using genetic and direct optimization algorithms.
Moreover, in addition to the calibration process, one of the existing challenges in DEM is the high computational cost needed to perform real industrial-scale simulations because of the extremely large number of particles being studied simultaneously at every time step. One of the most efficient approaches utilized to overcome this "obstacle" is to scale up the particle sizes and thus, reduce the total number of particles in the system according to coarse-graining method [39][40][41][42][43]. The most challenging issue in following such a method lies in adapting the material properties in such a way that compensates for the enlargement of the particles. In this study, the coarse-graining scheme proposed by Bierwisch et al. [42] was applied. The basis of this scheme is the preservation of the kinetic and gravitational potential energy densities as well as the volume fraction of particles as the original system with unscaled particles.
In this study, a novel simplified semiautomated parametric study and calibration approach is developed using an optimized multivariate regression analysis (MVRA). Initially, experimental powder characterization and AoR tests are performed. DoE method, using quadratic D-optimal design model, are applied for the sake of minimizing the number of simulations in an optimized way that retains high significance by recommending suitable possible combinations between DEM microscopic input parameters that should be tested together. The D-optimal design proved its efficiency by building a robust model by saving a lot of time and efforts compared to other models that require much more data points such as full factorial design. The proposed simulation runs from the DoE were performed in a parallel manner by implementing a code programming through the open-source software LIGGGHTS ® [44]. The studied microscopic input parameters are the coefficients of static and rolling friction, coefficient of restitution, particle size, Young's modulus and cohesion energy density. The effects of these input parameters on five macroscopic output responses (targets), i.e., AoR, porosity, mass flow rate, translational kinetic energy and computation time, are studied through MVRA, using Cornerstone software, to find out which parameters and in what sense do they influence these targets. This systematic approach was done for a free-flowing powder and a cohesive powder, namely SpheroLac 100 and Inhalac 251, respectively. This method gives a more robust understanding of the impact of the interparticle input parameters on the flowability of the powders at different ranges (low, center and high) by tuning them in a methodical way rather than trial-and-error. Since it is possible to have the same AoR with various DEM input parameter combinations [45,46], the porosity is measured along with the AoR in DEM simulations in order to determine a more reliable and practical value for the studied input parameters. Accordingly, a prediction model is built to validate the regression analysis and to calibrate our experimental AoR and porosity values by selecting one optimized parameter combination.

Materials
In this study, two pharmaceutical crystalline lactose monohydrate powders, purchased from MEGGLE Wasserburg GmbH & Co. KG (Germany), were selected. SpheroLac 100 was used to model the free-flowing powder, and InhaLac 251 was utilized to model the cohesive one. They are mainly characterized by different particle sizes, which in turn affect their flowability. SpheroLac 100 is mainly used in capsule filling, blends, premixes, sachets, triturations, whereas InhaLac 251 is a dry powder inhaler used in capsule-or blister-based formulations.

Experimental angle of repose (AoR)
The AoR is defined as the slope of a poured conical pile of loose uncompacted bulk solid material [18]. It is a characteristic related to the cohesiveness of the powder [47] and used as a common method to characterize powders according to their flowability [48][49][50][51]. It is a common method to calibrate DEM simulation parameters, where the angle at which a heap of powder settles depends on the strength of the interparticulate friction. This simple test case was chosen in order to perform a parametric study of the DEM microscopic input parameters. As depicted in Fig. 1, a funnel made of 1.4404 stainless steel was placed above a catching container made of borosilicate 3.3 glass. Initially, the funnel was loosely filled with 5 g of powder where the lower opening of the funnel was closed. Then, the funnel was opened to let the powder discharge under gravity forming a pile on the bottom of the catching container. The test for each powder was repeated five times, and an average of the obtained angles was considered as a value for the AoR for each powder.

Discrete element method (DEM)
In DEM, grains, referred to in the following as particles, are treated as rigid bodies, which have translational and rotational degrees of freedom assigned to their center of mass. In this study, the focus is on the behavior of the free-flowing and the dry cohesive powder particles. In the DEM, particle interactions are modeled using the well-established Newton's equations of motion, contact laws and overlap relationships. The soft-particle approach, originally developed by Cundall and Strack [3], is followed in this work, where the particles are allowed to undergo small deformations, i.e., overlaps. Thereafter, these deformations are used to calculate the elastic and plastic deformations besides the frictional forces between the particles [4]. As two particles i and j come in contact with each other, different interactions, i.e., forces and torques due to gravity, deformations under collisions, as well as static and rolling frictions, might occur. Translational and rotational motion of particle i with mass m i and moment of inertia I i are represented by the following equations: where g , v i and i are the gravity vector, translational velocity and rotational velocity of particle i, respectively. The interaction between particle i and j at a defined time step leads to normal and tangential forces, i.e., F n ij and F t ij , respectively. R i is the vector between the center of particle i and the contact point, where the tangential force F t ij is applied.
r ij is the torque due to rolling friction between the two particles.

Hertz-Mindlin contact model
Hertz-Mindlin [52,53] contact law is the most common applied approach in modeling contacts between particles. It is used in calculating the contact forces introduced due to the slight considered overlap between particles. The frictional contact force between two granular particles i and j comprises normal F n ij and tangential F t ij nonlinear contact forces given as Each contact force consists of two terms, where the first term stands for the nonlinear elastic Hertz model in the normal direction (spring force) and the linear elastic Mindlin model in the tangential direction (shear force). The second term stands for the damping force (modeled as a dashpot), where a dissipative term is added for both normal and tangential directions to account for energy dissipation during collisions through inelastic deformation and friction. The normal and tangential forces between two spheres i and j during a collision are given by where k n and k t are the elastic constants for the normal and tangential contact, n ij is the normal contact overlap, given by | n ij | = R i + R j − l ij , where l ij is the distance between the centers of the two particles, t ij is the tangential contact overlap or displacement vector which is truncated to satisfy a frictional yield criterion, given by the integral of the tangential relative velocity through the collision time they are in con- and v t ij are the normal and tangential components of the relative velocity of the two particles at the contact, n and t are the viscoelastic damping constant for normal and tangential contact. Figure 2 illustrates the mechanical system of the interactions of particles in the Hertz-Mindlin contact model. Under sufficient tangential forces, particles will slip relative to each other or relative to surfaces they are in contact with. For noncohesive particles subjected to a constant normal force, the extent of slippage under tangential force is determined by where ss , represented by a shear slider, is the static friction coefficient between the two particles in the tangential direction. Thus, the tangential force described by a spring and a damper is limited to the Coulomb friction.
Consequently, based on Eqs. (4) and (5), the normal and tangential forces, F n ij and F t ij can be expressed in detail as where Y * is the equivalent Young's modulus of the two colliding particles. This is defined by 1 where v i and v j are Poisson's ratios, R * is the equivalent radius, defined by 1 , and m * is the equivalent mass, R * n ij are the normal and tangential contact stiffness with G * being the equivalent shear modulus, defined as Besides, is the damping ratio coefficient, which is a function of the coefficient of restitution, e , and given by = ln (e)∕ √ ln 2 (e) + 2 . The coefficient of restitution e is the ratio of relative velocity after collision to the relative velocity before collision. A value of e = 1 refers to a perfectly elastic collision, whereas e = 0 would mostly be an inelastic collision.

Rolling friction model
In addition to sliding friction, rolling friction is present as a contact-dependent force [55]. Rolling friction is a major parameter that serves as a correction factor by compensating for the use of spherical particle shapes in DEM [56], instead of representing their actual shapes in reality. In the literature, different rolling resistance models have been developed and reviewed [56]. In this work, the alternative elastic-plastic spring-dashpot (EPSD2) rolling resistance model [57], implemented in LIGGGHTS open-source code, is applied to calculate the torque r ij,t+Δt at an incremental time step as where r is the rolling coefficient. In this model, the damping torque is disabled for simplicity.

Cohesion model
The modified simplified Johnson-Kendall-Roberts cohesion model (SJKR), implemented in LIGGGHTS, is used to model the cohesion of particles. This model originates from the JKR model [58], where it represents the influence of various cohesive forces, such as Van der Waals within the contact zone. It can model high adhesive systems, such as fine and dry powders or wet materials. It adds an additional normal force component F coh to the Hertzian contact law tending to maintain a larger contact area (Fig. 3), i.e.
where C is the cohesion energy density in J m −3 and a is the particle contact area, given as where R i and R j are the radii of sphere i and j, respectively, and l ij is the distance between the particle centers.

Noncontact force
In addition to the contact forces, the gravitational acceleration g in the z-direction (with ẑ as the unit vector) yields the gravitational force of particle i, F grav i , which can be expressed as

Time step in DEM simulation
In the DEM simulations of particles flow in this work, the time step (Δt) has been selected to be less than 20% of the critical time step (Rayleigh time) T R in order to ensure the stability of the integration scheme. T R is calculated based on the Rayleigh wave propagation, which represents the elastic wave propagation along the surface of one particle to the adjacent contacting particle according to the following equation [59]: where ρ is the particle density, G is the shear modulus and is the Poisson's ratio.

Experimental material characterization
This section briefly describes the experimental methods used to find out some material parameters utilized in the DEM simulations. The parameters are divided into two categories: material parameters and interaction parameters [60]. The material parameters include the shape, size, density, elastic and shear modulus as well as Poisson's ratio, whereas the interaction properties are coefficients of restitution, static and rolling frictions and adhesion.

Particle size and shape characterization
Camsizer XT from Retsch Technology GmbH (Germany) was used to measure the volume moment mean (D [3,4]), which is also called De Brouckere Mean Diameter, and shape of particles via dynamic image analysis. The D [3,4] indicates around which central point (center of gravity) of the frequency the volume/mass distribution would rotate. The advantage of this method over other particle sizes calculation methods is that it does not require the number of particles in its calculation.
The Camsizer XT helps in approximating the particle shape, where the particles are described in two-dimensional images in terms of the minimum and maximum Feret diameters (Fe min , Fe max ). Fe min is the length of the minor axis and Fe max is the length of the major axis. The ratio of Fe min to Fe max is called the aspect ratio (ASR). In this, the aspect ratio with values between 0 and 1 reflects the elongation of a particle and deviation from a sphere where small values correspond to elongated particles and higher values correspond to spherical particles: According to Li et al. [61], a shape factor called circularity or roundness (R o ) is an additional factor that indicates to what extent a particle has a spherical shape. It is based on the projected area (A) of the particle and the overall perimeter of the projection (P) according to the following equation:

Density and porosity
The bulk density (BD) and the tapped density (TD) were analyzed via the jolting volumeter (Jel STAV II) from J. Engelsmann AG (Germany). A certain mass of powder was filled into a graduated cylinder, and the volume and mass were then recorded. The test was repeated five times to insure accuracy and repeatability where then an average value was considered. The bulk density represents the average density of the powder sample, where the bulk volume contains voids between individual particles. Therefore, the bulk density of a powder depends on both the material density of the particles and their spatial arrangement. The volume of a bulk solid depends on the size and the shape of the particles. Besides, the tapped density was attained after mechanically compacting the powder sample by up to 1250 taps.
The particle (true) density (ρ p ) was obtained using the Ultrapyc 1200e, Automatic Density Analyzer from Quantachrome Instruments (USA). It is a pycnometer that measures the true density of powder samples using an inert gas (helium) to measure its volume by applying Archimedes' principle of displacement and the technique of gas expansion (Boyle's law). Consequently, the powder bed porosity (ε) is calculated as described in the following equation:

Carr's index
Based on the BD and TD, powder flowability of bulk solid could be characterized by the Carr's index (CI), or Carr's compressibility index [18]. CI is an indication of the compressibility of a powder, where larger changes indicate poor flowability [62]. It is a measurement of the relative difference between tapped volume and untapped volume of the powder [47], i.e.

Cohesion and flow function coefficient
FT4 powder rheometer from Freeman Technology Ltd (UK) was used to find out the cohesion (C) and flowing factor ( ff c ) values. A commonly applied approach is the Mohr-Coulomb model that fits the Mohr stress circles to the yield locus that in turns identifies the major principal stress ( MPS ) and unconfined yield strength ( UYS ). In this, UYS is related to MPS via the material flow function coefficient ff c as Flow function ff c is commonly used to describe the flowability of a bulk solid. In addition, the cohesion value is identified by the shear cell test [63]. This value is the shear stress at which the normal stress is equal to zero for the best fit of the shear and applied normal stress plot [64]. Generally, a cohesive powder will have higher values for cohesion and UYS and consequently a low ff c . In particular, ff c values below 4 denote a very cohesive powder behavior (poor flowing), between 4 and 10 easy-flowing powder and above 10 free-flowing powder.

Static friction
The shear cell test in the FT4 powder rheometer was executed to find out an approximation of the angle of internal friction (AIF) between particles based on the linearized yield locus of the material. It is the angle between the axis of normal stress (abscissa) and the tangent to the yield locus [65]. It measures the shear stress for various normal stresses to describe the magnitude of the shear stress that powder can sustain [65]. From this angle, the interparticle static friction coefficient, denoted by µ s,pp , indicates the resistance between two contacting particles while being sheared [2]. This can be calculated as

Wall friction
The wall friction test, using the FT4 powder rheometer, provides a measurement of the sliding resistance between the powder and the surface of the process equipment [66]. This test is important for understanding the discharge behavior from a container, especially when using an auger. The measurement principle of a wall friction test is very similar to the shear cell test, where the powder is sheared against a material resembling process equipment wall and not against the powder. The process equipment wall used in this study is a rounded stainless steel disk with 0.05 µm surface roughness. The data from the test are represented as a plot of shear stress against normal stress, allowing to determine the wall friction angle (WFA). It is the arctan of the ratio of the wall shear stress to the wall normal stress [65]. The greater is the wall friction angle, the higher is the resistance between the powder and wall.
Similarly, the coefficient of static friction between a particle and wall (µ s,pw ) can be calculated from the wall friction angle as

Moisture content
According to Schulze [18], the holding forces or adhesive forces with the presence of liquid bridges (capillary force) are higher than any other holding forces in bulk materials, as long as the contact distance is very low. The content of moisture sheds the light on some characteristics such as cohesion or agglomeration tendency and can help to describe the flow behavior. The moisture content (MC) was experimentally determined by the halogen drying method using the Halogen Moisture Analyzer type HR73 from Mettler Toledo (USA). A sample amount was heated by a halogen lamp, and thus the moisture is extracted from the powder. Due to this drying process, the bulk loses some weight, and thus the moisture loss is calculated from the difference between the initial and final weights. This moisture loss corresponds to the moisture content of the powder in the initial state.

Coarse-graining method
To find a trade-off between computational time and accuracy, it is reasonable to limit the number of particles in the DEM simulation. According to Bierwisch et al. [42], coarse-graining (CG) method is a good approach to save significant computational time but DEM simulations with good accuracy are performed in comparison with real tests. In the coarse-graining technique, the original particles with radius R are substituted by larger grains with radius R ′ . The scaled-up radius R ′ can be represented as a multiplication of the original grain radius R and a scaling factor S f . Moreover, CG method is based on an appropriate adjustment of interaction laws and equations of motion according to scaling rules (20) s,pw = tan(WFA).
for the material parameters in such a way that the energy densities, i.e., gravitational potential and kinetic energy, are preserved in the scaled-up system as the original one. As the potential density does not depend on the particle's radius while maintaining a constant volume fraction and particle's density, the particle's density is ought to be fixed upon applying the CG method. According to the coarse-graining proposed by Bierwisch et al. [42], the volume fraction is not affected and, therefore, the potential energy of the scaled system is comparable to the original system. Consequently, the mass of each particle introduced in the system is scaled up by a factor of S f 3 according to Eq. (21), but as the total number of particles is decreased also by the same factor, the total mass of the particles is conserved: Bierwisch et al. [42] also proved that the energy dissipation per volume and time is preserved if the coefficient of restitution remains the same after applying the coarse-graining method. Thus, the kinetic energy density is preserved because the scaling does not affect particle velocities. In connection with the applied test in this study, i.e., AoR, it is shown that coarse graining has no considerable effect on the resulting angles. In addition, it was observed that the shape of the heap is changing with increasing the particles sizes, where a noticeable tail of the ground heap and a more procumbent peak could be seen for larger particles. Thus, although the CG method leads to highly acceptable results in terms of bulk properties, it is hard to get exactly the same shape of the heap in the DEM simulation after applying the CG method.

Design of experiment
Design of experiments (DoE) is a systematic technique of planning experiments that is used to determine the cause-effect relationship between factors affecting a process and its responses [67]. It can help in identifying the best combinations of parameters that yield to the desired results and, thus, optimizing the output [32]. The design space of the experiments, in our case simulations, is identified according to the design model used by efficiently selecting the suitable factors combinations to be tested. The DoE method is used to investigate the effects and interactions of DEM input parameters on the output responses, which in turn will facilitate the calibration process. Quadratic D-optimal design model was used to obtain the largest amount of information from the smallest number of simulations using Cornerstone 6.1.1.1 software from camLine GmbH (Germany). The mathematics behind D-optimal design is closely related to the Least Squares method to find the concrete model equation on the basis of the design and the response data.

Multivariate regression analysis
Multivariate regression analysis (MVRA) is a statistical technique applied to analyze the variation of more than one outcome variable and thus estimates a regression model. This technique was utilized to observe and analyze the variation and impact of different DEM input parameters on the output parameters. Cornerstone software was used, where the inserted data values are standardized; that is, the original variables are scaled to have a mean of 0 and a variance of 1 as each attribute has a different unit of measurement.
The goodness-of-fit of the statistical model was described by the adjusted R 2 (coefficient of determination).

Fixed parameters through the studies
In the DEM simulations, the Hertz-Mindlin contact model and the rolling friction model (EPSD2) were applied for both the free-flowing and the dry cohesive studies. The adhesion contact model (SJKR) was added to the cohesive study. The values of Young's modulus (Y) and the Poisson's ratio (ν) for the funnel and the catching container were selected from the literature [69]. The boundary conditions such as restitution and friction between particle and walls (funnel and glass in this case) were kept constant to reduce the number of simulations in the DoE. Table 1 lists the input parameters in the applied two DEM AoR studies.

Parameters varied in the design of experiments
In the free-flowing study, five particle properties were varied as part of the DoE which were as follows: particle's Young's modulus (Y), particle size (d), coefficient of restitution (e pp ), coefficient of static friction (µ s,pp ) and coefficient of rolling friction (µ r,pp ). In addition to these parameters, cohesion energy density (C pp ) was varied in the dry cohesion study. The latter four parameters are related to the particle-to-particle interactions.
In the calibration process, the variation study of the parameters aimed to give the best fit between the DEM simulation and the corresponding reference data. To this end, the µ s,pp was determined using the FT4 powder rheometer, but because of reducing the Y parameter due to computational limitations, this coefficient was calibrated to get the same bulk solid behavior [45,70]. e pp and µ r,pp were also tuned because they are difficult to be determined for microscale powders. When using spherical particle shape, the rolling resistance is considered to be less dominant above the static friction as the contact area is smaller. The range value of the cohesion energy density was roughly estimated by performing presimulation runs or so-called initial screening investigations. Table 2 lists the values of the six particle properties that were varied in our studies. Three values, which constitute a low, center and high levels for each parameter, were Three-level design was used where each factor was considered at three levels. For the free-flowing powder, 40 simulation runs were executed whereas 53 runs for the dry cohesive one.  1 ν 0.21 Poisson's ratio particle ν 0.3 Coefficient of restitution particle-wall e pw 0.2 Coefficient of rolling friction particle-wall µ r,pw 0.2 Coefficients of restitution and friction wall-wall 0.2 Acceleration due to gravity g m/s 2 9.81 Cohesion energy density particle-wall 2 C pw J/m 3 1000 Cohesion energy density wall-wall 2 C ww J/m 3 100

DEM studied responses
In preparation of the mesh for the simulations, all the faces that have no contact with the particles were deleted in Solid-Works in order to optimize the computation time. Figure 4 shows the modified mesh file of the funnel in ParaView filled with spherical particles as an initial state before removing the stopper and allowing the particles to flow out. Five macroscopic output responses, namely AoR, porosity (ε), mass flow rate ( ṁ ), overall translational kinetic energy (TKE) and the computation time (CT), on the two powders were studied. The porosity was measured in the funnel due to the simplicity in defining a common volume since no heap (zero AoR) was formed in some simulations. The ṁ at the discharge of the funnel was studied, where the mass values were measured every 25 ms. The overall translational kinetic energy (TKE) of the system was calculated over time to check how the dissipation of energy varies according to different simulation input combinations. The average of the TKE, after opening the funnel, was calculated and then normalized (NA-TKE) with the total discharge time of the particles. Table 3 reports the DoE followed and its simulation results of five output responses under the variation in five input parameters for the free-flowing powder. Table 4 shows the DoE used and the corresponding results of five output responses under the variation in six input parameters for the cohesive powder. Table 5 presents the results of the particle size and shape, as well as the bulk properties of the two powders. Based on these material characterization results, the spherical form was selected to represent the particle's shape in LIGGGHTS because both powders have high ASR and R o values. The values of the ρ p were used in the DEM setup. The low values of MC in the table prove that both powders are almost dry. Table 6 lists the flowability values according to Carr's index as well as the shear properties of the two studied powders. The results illustrate how SpheroLac 100 is considered a free-flowing powder because it has low CI, low C and high ff c . Conversely, InhaLac 251 has high CI, high C and low ff c values and thus it is considered as a cohesive powder. The values of µ s,pp were just set as reference values, whereas the values of µ s,pw were used in the simulation setup.

Angle of repose experiment
The AoR for each experiment was calculated by taking the average value between the left and right angles, θ 1 and θ 2 , respectively. The obtained results of the AoR for both powders were set as reference values in the validation of the regression analysis in the prediction model and thus in the calibration approach. Figure 5 shows the heap formed by the AoR test of the free-flowing powder (SpheroLac 100) in the glass-catching container where the average angle was 33°.

Factors affecting angle of repose
• Free-flowing The average AoR for each simulation test was calculated. The values of AoR ranged between 0° and 39°, where in eight cases the AoR was equal to zero due to low µ s,pp as listed in Table 3. Figure 7 presents numerical examples of six AoR DEM simulations for the free-flowing powder formed in the catching container.
The model of the AoR for the free-flowing powder showed an adjusted R 2 = 0.954, indicating a very good fit to the data points. Figure 8a shows the adjusted response graph for the AoR. Each of the presented sections (columnwise) corresponds to an input parameter. In each section,  all the data points obtained from simulations are presented in addition to the effect drawn from these points. For a data point, the effects of the other factors included in the model were averaged; therefore, the data points were adjusted in such a way that the effect of the corresponding factor on the response is only presented (adjusted response graph). As the data points were recalculated, it can be seen that the actual values of the responses listed in Table 3 do not match with the adjusted values, yet the effect of each parameter can still be modeled. Finally, scattering points signify data points that cannot be described by the model due to some noise effects. From Fig. 8a, the following trends can be observed regarding AoR: Y and d S have no significant effects on the AoR, e pp and µ r,pp have a slight effect and µ s,pp has the major effect. The absence of d S as an effective factor to the AoR ensures the correct implementation of the coarse-graining scheme proposed by Bierwisch in our system. Moreover, our results contradict many studies which proved that the particle size has an impact on the AoR [22,25,47]. Figure 9 illustrates the Pareto graph of effects on the AoR of the free-flowing powder where the largest positive effects are shown on the left and the largest negative effects on the right, with factors of smaller estimated effects in the middle of the graph. It can be seen that µ s,pp has the highest positive effect (proportional) on the AoR, where, as µ s,pp increases, the AoR increases as well. The rest of the affecting factors, i.e., µ r,pp and e pp , interactions (e.g., µ s,pp * µ r,pp ), and quadratic terms (e.g., e 2 pp ) have lower effect compared to that of µ s,pp . To be stated, as e pp increases, the AoR decreases, which is previously confirmed by Wei et al. [35]. As e pp is the ratio of the relative velocity after collision to the relative velocity before collision of an object with a surface, the results can be explained by the fact that due to the high elasticity of the spherical particles, they can more likely bounce than to settle down along a free surface.
• Cohesive Figure 10 illustrates numerical examples of six AoR DEM simulations for the cohesive powder formed in the catching container.
The adjusted R 2 of the AoR of the cohesive powder model has a very good value of 0.915. Figure 11a shows the adjusted response graph for the AoR. It indicates that the introduction of cohesion in the system through the SJKR model added a new affecting parameter on the AoR, which is Y. Unlike the free-flowing case, it is proved that Y has a significant negative effect on AoR when C pp is involved. This ensures the importance of such a parametric study, where Y is not included in many other calibration processes  due to its computational cost impact, though it has influence. From Fig. 11a, it can be seen that the angle of repose increases with an increase in the sliding and rolling friction coefficients, which is proven in other studies as well [25,46,[71][72][73][74]. This positive correlation is similar in both freeflowing and cohesive powders. Figure 12 demonstrates the Pareto graph of effects on the AoR of the cohesive powder. It shows that µ r,pp has the  highest positive effect on the AoR, which is contrary to most of the studies that proved the superiority of µ s,pp in terms of effectivity. Therefore, it is very critical to always relate the effectiveness of the parameters to the circumference of the case study. Moreover, Y has a bigger impact on the AoR than µ s,pp .

Factors affecting porosity
• Free-flowing The regression analysis of the ε of the free-flowing powder resulted in a very good adjusted R 2 value of 0.981. The adjusted response graph of ε in Fig. 8b demonstrates the high significance positive correlation of µ s,pp on ε followed with µ r,pp , which was also proven by Wei et al. [35]. The other studied microscopic parameters indicated less significance, where the relative size of effects of these parameters is presented in Fig. 13.

• Cohesive
The adjusted R 2 of the ε of the cohesive powder was very good as well with a value of 0.983. Similarly to the freeflowing powder, µ s,pp had the highest positive influence on ε followed with µ r,pp as revealed in Fig. 11b. In addition, the Pareto graph of effect in Fig. 14 shows the superiority positive relative effect of µ s,pp on ε compared to the other parameters.

Factors affecting mass flow rate
• Free-flowing The adjusted R 2 of the ṁ model of the free-flowing powder was recorded to be 0.95. The adjusted response graph for the ṁ is shown in Fig. 8c. As observed, all the studied parameters affect the ṁ , but each in a different way which is also presented in the Pareto graph (Fig. 15). The Pareto graph shows that µ s,pp has the highest negative effect on the ṁ . This indicates that as µ s,pp increases, the particles' flow out of the funnel gets delayed due to the sliding friction Following this, the mass and ṁ of the six example simulations shown in Fig. 10 were plotted over time in one graph ( Fig. 16) to compare them with each other. The parameter combinations and the mean ṁ values of each run are presented in Table 3. Figure 16 shows that the slopes of mass graphs M_1 and M_9, corresponding to the first and ninth DEM simulations, have the highest ṁ while having the lowest values of µ s,pp (0.1). M_8 is showing the lowest slope ( ṁ ) with the assignment of having a high value of µ s,pp (0.85).

• Cohesive
The model fit of the ṁ in the cohesive study was very good with an adjusted R 2 = 0.957. Figure 11c presents the adjusted response graph for the ṁ , and the Pareto graph (Fig. 17) shows the weight value effect of the input parameters on the ṁ . It is apparently shown that ṁ has been highly negatively affected by µ s,pp , µ r,pp and d I in decreasing order of effectivity. A lower but significant positive effect has been witnessed for Y. This is due to the fact that particles with high stiffness tend to flow easier under the effect of gravity.
More clear visualization of how ṁ is varying with respect to the parameters can be seen in Fig. 18, which shows an example of the mass and ṁ of six simulation cases plotted over time. The varied input parameters and the mean ṁ values of each run are presented in Table 4.

Factors affecting kinetic energy
• Free-flowing The overall translational kinetic energy (TKE) was recorded every 1 ms starting from the moment of particle insertion inside the funnel till the discharge of the last particle inside the catching container. The results showed that µ s,pp has a negative impact on the TKE as illustrated in an example of three runs in Fig. 19.
From 0 to 1 s, the particles were inserted inside the funnel, and then there was 0.3 s settlement time for the particles before removing the stopper at 1.3 s. This resulted in having an almost identical character of variation of TKE from 0 to 1.3 s. The rapid increase in TKE within 0-0.1 s is explained by the falling (insertion) of the particles into the funnel.      Within 0.1-1 s, the distance between the insertion region and the first contact of the particles in the vertical direction decreased as the funnel was being filled with particles. Then, as the particles were settling down, the TKE started decreasing. In the time between 1 and 1.3 s, a rapid decrease in the TKE was observed as no more particles are inserted and most of the particles are got in rest (settled down). It can be seen that upon removing the stopper, the peak of E kin _7 was smaller than that of E kin _5 and then E kin _7 remained smaller and the total time for full discharge was longer compared to E kin _5. The only difference between these two runs is the value of the µ s,pp where it was 0.1 in run 5 and 0.85 in run 7. So, the high value of µ s,pp decreases the overall TKE of the system. On the other hand, E kin _2 has a smaller e pp than E kin _5 which showed a minimal difference in a small interval of time during the discharge of the particles. Thus, e pp has no significant influence on the overall TKE when µ s,pp is low. On the contrary, e pp has a significant influence on the overall TKE when µ s,pp and µ r,pp are both high as shown in Fig. 20.
Moreover, the particle size d S has a major influence on the overall TKE. Figure 21 plots the overall TKE graphs of simulation cases 1 and 24 where the d S in run 1 is 460 µm and the d S in run 24 is 690 µm.
The graph plot shows that both runs have similar overall TKE during the insertion time. After removing the stopper, the overall TKE of run 1 with lower d S (460 µm) is higher and reaches full discharge faster than run 24 with higher d S (690 µm). This can be explained by that the number of particles decreases with the increase in radii keeping the same mass; therefore, the number of collisions decreases and thus leads to a decrease in the TKE.
For a better understanding of the relative level of effectivity of these parameters, the adjusted response graph and Pareto chart are also designed for this response. The adjusted R 2 of the NA-TKE model of the free-flowing powder was 0.981. In Fig. 8d, the adjusted response graph for the NA-TKE is illustrated. The Pareto graph of effects on the NA-TKE of the free-flowing powder is shown in Fig. 22. It shows that the quadratic term of µ s,pp has the highest positive effect and µ s,pp has the highest negative effect on the NA-TKE, followed by µ r,pp . This result harmonizes perfectly with the definitions of these two factors, as they represent, more or less, the extent of resistance of a particle to remain in its position upon being exposed to a certain force. Therefore, as they increase, the particles are more hindered, and therefore the TKE decreases. And since the particles in our system do not tend to move in the rotational direction as in the translational direction, the µ r,pp is showing less effect than µ s,pp .

• Cohesive
The model fit of the NA-TKE of the cohesive powder was very good with an adjusted R 2 = 0.97. The adjusted response graph for the NA-TKE is shown in Fig. 11d. The Pareto graph of effects on the NA-TKE (Fig. 23) shows that Y has the highest positive effect and µ s,pp has the highest negative effect on the NA-TKE. As seen, a new factor, namely Y, compared to the free-flowing case, has been introduced as a significant parameter. The reason for that is the existence of the C pp in the system representing an additional normal contact force tending to stick the particles more to each other. As Y denotes the stiffness of the particles, exerting an  additional force on them along with less tendency to deform leads to some repulsion causing the higher velocity of the particles. This interactive effect of these two parameters can be seen also by the term C pp * Y having a positive impact on the NA-TKE. On the contrary, maintaining the same value of Y with increasing the C pp , decreases the NA-TKE, because the additional normal contact force ought to keep the particles in contact. It can be also observed that the particle size has a higher significant effect than in the free-flowing case.
In fact, this can be justified by the increase in contact area along with the increase in normal contact force, the thing that would cause an increase in the frictional force in order to maintain a constant µ s,pp . In Fig. 24, the overall behavior of the TKE of the system in the cohesive case generally looks similar to the freeflowing case. Some phenomena discussed earlier regarding the effects of some parameters can be seen.

Factors affecting computation time
• Free-flowing The CT model of the free-flowing powder had a very high adjusted R 2 value of 0.987. It is clearly shown in Table 3 how the CT is varying through the different simulation runs. The adjusted response graph for the CT is illustrated in Fig. 8e. In addition, Fig. 25 demonstrates the Pareto graph of effects on the CT of the free-flowing powder. It is indicated that Y has a big positive impact on CT, whereas d S has a smaller negative effect on it. In other words, as the Y becomes bigger, the CT increases and as d S increases, the CT decreases. Although e pp has much less effect on CT compared to Y and d S , yet it is important to notice that e pp has a negative effect which might be used in another applications needing relatively high computation time.

• Cohesive
The adjusted R 2 of the CT model of the cohesive powder was very high as well (0.991). Figure 11e illustrates the adjusted response graph for the CT, and Fig. 26 presents the Pareto graph of effects on the CT of the cohesive powder. Similarly, the results showed that Y and d I are the almost

Prediction model and experimental validation
The main aim of this section is to validate the regression models of the parametric studies obtained above and to validate the numerical model based on experimental data. Therefore, the reference values of the AoR and ε determined from the experiments will be our targets in the following prediction simulations.

Free-flowing
Based on the results obtained, a prediction model was built to estimate the AoR and ε upon using an optimized combination between them. The predicted AoR and ε for the free-flowing powder are calculated according to Eqs. (22) and (23), respectively, as follows:  The constants in both equations are presented in Table 8.
As the value 0.61, which represents the experimental value of µ s,pp of the free-flowing powder using the FT4 rheometer, is within the tested range of µ s,pp , it was fixed in the prediction model. In addition, according to the obtained results, Y was set to minimum and d S was set to maximum to reduce the CT. Using an integrated optimization feature algorithm in Cornerstone software, an optimized DEM parameter combination of µ r,pp and e pp , which aims to match our studied bulk responses (AoR and ε) with the experimental reference values with minimum deviation errors, is selected. Table 9 shows the combination of the DEM input parameters used in validating the prediction model and thus calibrating SpheroLac 100 (the free-flowing) with its experimental AoR and ε values. The values to be predicted were set as our experimental reference values.
The results in Table 9 validate the robustness of our regression model, and thus our powder was calibrated. On this occasion, Fig. 27 shows the numerical AoR of run 41 validated with the experimental AoR of the free-flowing powder.

Table 9
Executed simulation run for the validation of the process model of the AoR and the ε of the free-flowing powder can be seen in this study that this depends on the application and not only on the type of powder being used. The prediction model of the cohesive powder was validated by running a new simulation having µ s,pp equal to the experimental value. Based on selecting the optimized parameter combination, Table 10 lists the results of the prediction run simulation of the AoR and ε compared and validated with the experimental data.
The very low deviation errors confirm the validity of our study approach, and thus InhaLac 251 (the cohesive powder) was calibrated as well. Moreover, Fig. 28 shows the AoR of run 54 compared with the experimental AoR of the cohesive powder.

Summary and conclusions
In the presented work, a semiautomated calibration approach was introduced in which two different powders, i.e., freeflowing and cohesive, were systematically studied and their parameters numerically calibrated using a DEM-experimental data and following regression analysis. The work aimed, firstly, to perform a systematic statistical parametric study to understand the impact of six interparticle DEM microscopic input parameters on four macroscopic bulk responses and on the computational time. Secondly, numerical validation of the experimental angle of repose and porosity values using a statistical prediction model was performed. The examined microscopic parameters were Young's modulus, particle size, particle-particle coefficient of restitution, static and rolling friction coefficients and cohesion energy density. Interaction properties between particle and wall and between wall and wall were set constant, whereas particle-particle properties were extensively studied. Some of the highlights can be summarized as follows: • Applying powder characterization using measurement devices decreased the number of needed parameters to be calibrated, which improves the robustness of parameter tuning. • DoE proved to be an efficient tool in the calibration approach, whereas the regression analysis served in better understanding the impact of different DEM input parameters on the observed responses.   • The quadratic D-optimal design model ascertained its efficiency and reliability in saving computational time by setting up a small optimized number of simulations to be run. • For the free-flowing study, it was shown that the static friction coefficient has the most significant impact on almost all the output responses. The particle size and stiffness, however, had the highest impact on the computation costs. • Regarding the cohesive study, the results were in many cases similar to those of the free-flowing study except that the particle stiffness had a negative impact on the angle of repose, and the cohesion energy density had some impact on all the responses. It has been proven in this study that the involvement of cohesion increases the influence of the particle stiffness on some responses, e.g., the angle of repose. • Scaling up the particle size using the coarse-graining method proved to have no effect on the angle of repose in both studies. • The built statistical prediction model yielded a good quantitative agreement between the DEM simulation and the experimental result values.
Based on these results, future work will focus on additional numerical and experimental studies under different conditions to determine the validation scope of the selected parameters. One application would be running auger dosing experiments and studying the mass flow rate. This would allow to validate the DEM model using a real-scale machine, subjected to controlled initial and boundary conditions.