A 2D hysteretic DEM model for arbitrarily shaped polygonal particles

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Controlling powder flow and dosing for small scale processes, such as 3D food printing, is challenging. This is especially the case when the used materials are widely different and have properties ranging from soft viscoelastic to hard elastic, spherical-and cube-shaped to irregular and non-cohesive to cohesive. For the development of e.g. multimaterial 3D food printing [1] by deposition through actuated miniaturized hoppers [2,3,4] a model is necessary to precisely determine equipment and process parameters. A suitable model should be efficient, be able to simulate arbitrarily shaped particles with a range of material properties and be applicable for miniaturized hoppers.
The Discrete Element Method (DEM), first described by Cundall and Strack [5], is a standard way of modeling the behavior of granular systems by describing the interaction between the distinct particles that comprise the system. Particles are modeled as discrete elements, with a specific non-deformable particle shape, upon which contact and non-contact forces act. Every time step, the equations of motion for each individual particle are solved and their position and velocity are updated accordingly. These values along with inter-particle forces are used to describe parameters of granular systems, such as particle flow, packing density, stress patterns, etc. By describing particles as discrete non-deformable shapes and using simple force laws the discrete element method is suited for simulating large numbers of particles in limited computational time. Despite the simplification, DEM models can capture the behavior of granular systems well.
A myriad of force models have been developed for describing normal contact forces. Most fall in one of these categories, namely: linear spring-dashpot, simplified Hertz-Mindlin and Dersiewicz, linear hysteretic, non-linear or case based models or a combination of these, [6,7] all describing the macroscopic viscoelastic-plastic response of two contacting particles. Tangential contact forces are modeled with friction models for sliding and rolling friction, normally based on stick-slip [8,9,10], and models for non-contact forces, such as gravity, van der Waals and electrostatic are implemented to account for these forces. These force models are historically developed for spherical particles to minimize mathematical and computational complexity. Recently, methods have been developed that include particle shape in DEM models. Shapes are described in a multitude of ways, including ellipsoid, composites of primitive shapes, superquadrics and polyhedrons [11,12,13]. To correctly process non-spherical particles the force models are adjusted to accommodate the shape description of the particles. Detailed reviews of the different force models and particle shape descriptions are described by Zhong et al. [14] and Lu et al. [15]. Often the restriction is made that the particle shape should be convex, this stems from the fact that most models only were designed to handle single contact points between particles. However, most food particles are highly irregular making these models not suitable for simulating food printing processes. Especially for hysteretic force models in combination with irregular concave particles inconsistencies become apparent.
Modeling highly irregular shapes in 3D is computational intensive. When using accurate particle shape descriptions in the form of irregular polyhedron, the computational time will become significant. DEM simulations in 2D have an asset value particularly in new applications and DEM analysis [16], with the benefit of reduced computational time. In a 2D simulation particle rotation and translation in the third dimension is neglected. If this effect is expected to have only a minor influence on result, modeling in 2D is preferred. Recent studies using a 2D convex polygon DEM model vary from describing particle flow in fluids [17,18] to contact treatment [19] and development of contact detection methods [20]. As far as we know, no model using irregular polygons (or polyhedrons) in combination with a hysteric force model accounts for the merging and splitting of intersection areas. Rather the potential discontinuity in force is solved by treating contact areas independently.
In this work, we propose an extension of a 2D hysteretic force model for arbitrarily shaped particles that will be modeled by irregular polygons. This extension specifically focuses on how the force memory, inherent in any hysteretic force model, should be processed in the case when specific transitions of the intersection areas occur that are only possible with concave arbitrarily shaped particles. A 2D model is chosen for computational reasons and is sufficient for qualitative studies on determining the influence of key parameters on several 3D print processes. Care is taken to make the model extendable towards 3D by splitting the geometric intersection and the force calculation. Hence, the improvement to the force model in 2D is also applicable to 3D hysteretic force models. Specific attention is given to the case where multiple intersections are present between two single polygons. The model is validated based on packing density of gravitational deposited specific shapes.
In this paper, we will give a short description on how the particles are modeled as arbitrarily shaped polygons and how the main parameters are defined. An extension of a hysteretic force model will be given that has a consistent force trajectory during the merging and/or splitting of contact areas that can occur with arbitrarily shaped particles. The used contact detection method will be described along with the algorithm to detect whether the particle contact is a new contact, a continuation, a merging or splitting contact. As proof of principle the deposition of three different particle packings is shown with differently shaped particles.

DEM model for arbitrarily shaped polygons
The developed DEM model uses irregular polygons for the description of the particle shape. The 2D material can be thought of as 3D material with unity length in the direction perpendicular to plane of the shape description. The used force models are extended to be applicable to arbitrarily shaped polygons. Initially, the model will be described based on only single contact points between particles. Afterwards, the force model for multiple contact points between particle pairs will be discussed.
Simple arbitrarily shaped polygons are the discrete elements of the model. Simple polygons are here defined as: linear, closed and nonself intersecting without interior holes. Each polygon is defined by an ordered set of vertices {p 1 , …, p n }, where n is the vertex number, that represents the shape of a particle. The vertices are ordered in a counter-clockwise direction. Particles are defined by their shape, motion parameters, material properties and some bookkeeping information. Definitions for the shape and motion parameters are given in Fig. 1 for the i th particle. The position of a particle r i , is defined by its center of mass. The velocity, acceleration and attitude in a space-fixed frame of reference are v i , a i and θ i respectively. Here, the attitude in 2D is the angle with respect to the x-axis and is positive in the counter-clockwise direction.
A DEM simulation consists of tracking the state of a granular system by solving the equations of motion for the set of discrete particles. The position and velocity are determined by solving the differential equation where F i is the force experienced by the particle and m i the mass of particle i. The attitude, angular velocity and angular acceleration are similarly determined by where θ i , ω i and α i are the attitude, angular velocity and acceleration, respectively. M i is the torque on the particle with respect to its center of mass and I i the area moment of inertia. Force F i can be written as the summation of all forces acting on the particle, where F ij c represents the particle-particle contact forces between the i th and j th particle, F ij c the inter-particle non-contact forces imposed by particle k on particle i, F i damp an additional damping force and F i g the gravitational force. The torque M i , is the summation of all torques acting on the particle, where M ij slide , M ij roll and M i damp are the torques due to sliding, rolling and damping respectively and the cross product term is the torque caused by the normal contact forces acting on the particle. Where c ij is the position of the center of mass of the intersection area with respect to the particle center. How to define a single contact point and direction of the normal for non-spherical particles is still an debated issue in literature [14]. In our model, the center of mass of the intersection area is used as the contact point. This gives a continuous evolution of the trajectory of the contact point without discontinuities for a single intersection area. For the direction of the normal, a weighted intersection normal is used and will be detailed below.
The magnitude and direction of F ij c and M i are determined by the developed contact force model that relates the intersection area A ij of particle i and j and properties of these particles to force and torque. The direction of torque is always perpendicular to the particle plane and pointing into or out of the plane. The force direction is always parallel to the plane. The main parameters for polygon-polygon particle intersection area are shown in Fig. 2. The polygon-polygon particle intersection area is defined, similar to a polygonal particle, as a counter-clockwise polygon with a center of mass denoted by r ij and an intersection area of A ij . A geometrical intersection line is defined by For spherical particles the normal of the geometrical intersection line can be used as the direction of the normal force. For arbitrarily shaped particles this is not the case. However, the direction of the geometric intersection line running from begin point r ij b to end point r ij e is important since it defines the direction of the normal of the intersection n ij , pointing in the correct direction, meaning inward or outward. For arbitrarily shaped polygons it cannot be guaranteed that n ij will point outwards if it is simply defined as r ij − r i . To guarantee the direction of n ij , r ij b is defined as the first intersection point encountered in a counterclockwise direction when one starts traversing the vertices of particle i from outside the intersection. Subsequently, the other intersection point becomes r ij e . The mass distribution of the intersection is taken into account by using a weighted intersection normal, as amongst others, described by Matuttis and Chen [13]. For the weighted intersection normal, two intersection lines are defined, one running from the center of mass of the intersection to the begin point, given by and one running from the center of mass towards the end point The normals of these intersection lines are weighed to give the weighted intersection normal The weighted intersection tangential becomes by rotating the weighted normal intersection 90 ∘ in the counterclockwise direction. The weighted intersection normal and tangential are used as force direction. Besides the force direction several contact velocities are defined. The relative contact velocity is defined as where v i , v j , ω i and ω j are the translational velocity of particle i and j and the rotational velocity of those particles respectively. c ij and c ji are the positions of the center of mass of the intersection with respect to the particle centers. The normal and tangential components of the velocity become and respectively. For arbitrarily shaped particles it is important that an objective relative rolling contact velocity is used. An objective rolling velocity can be defined in several ways. Here the definition of Zhao [21] is used where the relative rolling velocity is

Normal force
Contact forces are determined from the intersection and particle parameters, based on the developed area dependent hysteretic force model. For equations only pertaining to two particles the particle indexes will be dropped from the terms in the equations. Hysteretic force models for spherical particles were developed, first by Walton Braun and later extended by a.o. Thornton, Vu-Quoc and Luding [22,23,8,24], to include plasticity as opposed to a velocity-based damping term as done in other models. A simple linear hysteretic force model is used instead of more involved non-linear models since it has been proven several times that more complicated models do not necessarily yield better results when compared to experimental data [6]. Forces in the developed model are based on intersection area instead of penetration depth making the resulting force shape dependent. A schematic of the developed force model is shown in Fig. 3. The spring stiffness coefficients k l , k r and k d are for the initial loading, unloading/ reloading and detaching respectively. For the case where both particles do not have the same coefficients the harmonic mean of the coefficients is used. The force trajectory can be split into three parts, the loading force given by k l A, the unloading/reloading force given by k r (A − A 0 ) and the detaching force given by −k d A. The loading trajectory is only followed if the force is greater than the previous maximum force F max during this particle collision, otherwise the unloading/reloading  trajectory is followed. The detaching trajectory is followed if the unloading force would be smaller than the detaching force. The magnitude of the hysteretic normal force can be written piecewise as.
where A is the intersection area of the specific particle-particle interaction, A 0 is the value for the intersection area below which the force becomes negative (attractive) and is calculated with A max is the maximum intersection area during a single particle-particle collision and should be kept in memory for the entire collision and updated accordingly when A > A max . For arbitrarily shaped particles with multiple possible intersection areas between two particles, A max can become ill-defined; this issue will be addressed in Section 2.4. The hysteretic force model dissipates energy only when traversing the transition points A min or A max . Within the reloading trajectory no energy dissipation is present leading to continuous vibration of densely packed systems. To alleviate this problem often a velocity-based damping term is added. Again for multiple intersection areas this will lead to unphysical behavior. To include damping, a term based on the derivative of the overlap intersection is used, much in line with using the derivative of the indentation depth for spherical particles. The normal force thus becomes where γ is a damping coefficient.

Tangential force
The particle experiences a tangential force due to friction and tangential contact velocity. A sliding-sticking friction model is implemented, based on the work of Luding [24], where the normal force is coupled to the tangential force via Amonton-Coulombs law. Rolling friction is implemented in a similar fashion.
Upon particle contact, a tangential spring ξ t is created, with t the time of the current time step and an initial length of zero, and will be kept in memory and updated during the complete collision. The tangential spring of the previous time step will be projected on the new tangential line with where ξ t−Δt is the tangential spring of the previous time step and Δt is the size of the time step. The projected static friction force is calculated, as well as the static and dynamic friction limit respectively, where ζ is a damping coefficient and μ s and μ d are the static and dynamic friction coefficients. The actual tangential force is dependent on whether the system is in a state of static or dynamic friction. Whether the frictional state is Static or Dynamic in the current time step is determined by the case rules: where state t is the frictional state of the current time step and state t−Δt of the previous one. Graphically, the different cases are shown in Fig. 4. Based on the relevant case the tangential spring is updated and the tangential force is calculated. Static Case: For the static friction case the tangential spring is incremented with the tangential displacement of the center of mass of the contact area and the tangential force becomes the with Eq. (17) calculated static friction force: Dynamic Case: For the dynamic friction case the tangential spring is made consistent with Coulombs condition and the tangential force is equal to the dynamic friction limit: Updating the tangential spring with Eq. (23) results in the projected static friction force in the next time step to be F s ≈ F t . The torque generated by the frictional force due to sliding becomes M ij slide = F t c ij .  Background damping is incorporated to include air resistance and overall damping of the complete system. The background damping force is, where γ b is the background damping constant. In a similar fashion rotational movement is dampened with, where ζ b is the rotational background damping constant and R i is the radius of particle i based on an area equivalent circular particle. Note, that these values should be kept minimal to avoid artificial overdamping of the system.

Multiple contact points
For convex particles only one overlap between two particles can exist. Concave particles can however have multiple contact points as shown in Fig. 5. Due to this, single contact force models cannot be used as will be shown.
As already mentioned, the additional damping is based on the rate of change of the contact area with a properly adjusted damping coefficient. Damping based on velocity is not correct for concave particles since the damping would scale with the number of contact points. The incorrectness of this can be easily visualized by taking a single overlap area and splitting this overlap in two by introducing an infinitesimal small gap while keeping all else constant. When handled as two separate contact points the damping will suddenly be accounted for twice increasing the damping twofold. A damping coefficient coupled to the rate of change of the contact area as shown in Eq. (15) avoids this. Since damping should only be present during the loading phase, one cannot have an attractive damping force, the damping force becomes which is input for the normal force expression in Eq. (15). The normal contact force law as detailed in Eq. (14) is also not suitable for multiple contacts. The force law is correct when multiple contacts stay separate but it neglects the possibility of contacts merging or splitting. A trivial case of this is shown in Fig. 5 where two separate contacts of one particle with another merge into a single contact and then split into two separate contacts again. The forces of the left and right intersection area before the merge and after the split are simply numbered in the subscripted and denoted F ij, 1 n and F ij, 2 n . With the contact force law of Eq. (14) this will result in a discontinuity in the force as shown in Fig. 6. The loading trajectory of the left and right contact follow the same curve until the intersection areas merge into a single intersection area. The red and blue line represent this loading trajectory and are plotted with the red line on top of the blue line. The yellow line represents the trajectory of the merged contact and the green and purple line represent the unloading and detaching trajectory of the left and right intersection after the split. As can be seen when two separate contacts merge during loading the forces are additive as should be and no special precaution is needed. Upon unloading, the force is not divided correctly over the two contact intersections as can be verified by adding the contact forces of both intersections. The force on the particle just before the split is 0.3 N while the force just after the split on the left intersection is 0.25 N and on the right intersection 0.75 N, resulting in a total force after the split of 1.0 N. The right intersection even has more force than before the split, so a force discontinuity occurs, resulting in a sudden addition of kinetic energy and thus velocity of the particle. This sudden addition is best visualized with the energy balance where the total energy of the system is with where E total , E pot , E kin and ΔE el−pl are the total energy, potential energy, kinetic energy and elastic/plastic energy respectively. h is the height of the particle above the plane where it would have zero potential energy, g is the acceleration due to gravity and s is the traveled path length. The energy of both the elastic and plastic deformation are captured in one term to easily verify that there is conservation of energy. The energy balance of the case shown in Fig. 6 is visualized in Fig. 7. At the transition point where a single intersection area is split into two separate intersection areas, a kink in the graph can be seen, at 3.4 ms. A sudden increase in kinetic energy and corresponding decrease in elastic/plastic energy can be seen correlated to the discontinuity in the force trajectory. This discontinuity is caused by an incorrect value of A max . A max must be recalculated with the information from the previous time step and  split accordingly. A max,1 and A max,2 should be calculated as if no merging ever happened. For an infinitesimal small time step between one intersection and an intersection split in two, the following must hold, Due to the linear nature of the used force law the maximum intersection area for intersection areas A 1 and A 2 scales with the weighted proportion of A s via, As is easily seen, this equation does not only hold for an intersection area split into two parts, but also in multiple parts. So calculation of the force after a splitting event is correct for any number of contact points between particles. The approximation is valid for small time steps and made for easier implementation so that no coupling of variables in this time step exists. The same trivial case as shown in Fig. 6 is presented in Fig. 8 with the force trajectory according to Eq. (33). As can be seen in Fig. 8, the summation of the force after the split of contact area is the same as just before the split. The force before the split is 0.3 N, the same as earlier, while the force on the left intersection is 0.075 N and on the right 0.225 N. So the addition of the left and right forces is 0.3 N and thus no discontinuity in force is present. The resulting energy balance shown in Fig. 9 shows a smooth energy evolution without any kinks as one would expect. The absence of the sudden increase in kinetic energy leads to a lower final value of the kinetic energy overall. Important to note is that correct handling of a splitting intersection area does not only affect the normal contact force but also the torque, resulting in quite different behavior on a particle level.

Algorithms for arbitrarily shaped polygons
For a DEM model with arbitrarily shaped polygons similar algorithms can be used as in conventional DEM with some exceptions. Algorithms for the determination of the intersection area, center of mass and inertia should be tailored for polygons. Algorithms for these properties are known and often implemented in geometry libraries. In our model the CGAL library [25] (Computational Geometry Algorithms Library) is used for calculating the intersection area of polygons with arbitrary precision arithmetic using the GMP [26] (GNU Multiple Precision Arithmetic Library) and MPFR library [27] (GNU Multiple Precision Floating-Point Reliable Library). Contact detection algorithms are also known, but the efficiencies of these algorithms are mostly given for spherical particles. For polygonal particles, efficiency not only scales with the number of particles but also with the number of vertices of the polygon. Especially the last step in any algorithm, checking whether a possible intersection is indeed an intersection is costly, so an algorithm with as little false positive intersection candidates is necessary. For any hysteretic force model details of the previous time step need to be in memory. For the described hysteretic model with multiple contact points the contact area also needs to be kept in memory for determining the merging and splitting of intersection areas.

Contact searching
A computational costly part of a DEM simulation is the contact detection algorithm. Naively checking whether spherical particles overlap would entail taking one particle and checking if it overlaps with any other in the simulation space by checking if R i + R j < ‖r i − r j ‖, where R is the radius of a particle. This method would result in a complexity of O n 2 À Á for n particles and is not feasible. Better methods are proposed in the form of neighborhood algorithms such as Verlet-Neighbor List (VL), Linked Cell (LC) and Linked Linear List (LLL), also called Sweep and Prune, where the complexities are from O n 5=3 À Á , O n ð Þ and O n log n ð Þ ð Þrespectively. The neighborhood algorithms yield potential contact candidates that in a final iteration need to be checked if they indeed overlap.
The LLL algorithm defines bounding boxes (BB) around the particles and keeps the boundaries in an ordered list. Each time step the ordered list is updated with an insertion sort algorithm where changes in order are linked to the creation or deletion of potential contact candidates. Graphically, this is shown in Fig. 10. For implementation details see   Powder Technology 378 (2021) 327-338 [28]. Although the algorithm is not the most optimal one for spherical particles, for polygonal particles it is more efficient due to the minimized boundary boxes. For polygonal particles the final iteration where the actual overlap is checked, is computational heavy and scales with the number of polygon vertices, while for spherical particles the actual overlap check is O 1 ð Þ. The minimized boundary boxes will lead to the smallest number of potential contact candidates and thus the LLL algorithm becomes efficient to use.

Merging and splitting detection
For arbitrarily shaped particles determining whether a contact is a newly formed contact or the evolution of a contact already present in the previous time step is not trivial. Additionally, distinctions must be made between new contacts, merging contacts, splitting contacts and normal evolution of contact. An efficient algorithm for polygonal particle is described which is independent of time step and only assumes a sorted memory structure for the polygon vertices.
Each time step all polygon intersections are determined and stored in memory. A polygon intersection as shown in Fig. 2 will always consist of at least three or more points, containing zero or more vertices of particle i, zero or more vertices of particle j and at least two intersection points that can, but do not have to, coincide with vertices of the particles. Two intersecting polygons result in one or more intersection polygons. Each intersection polygon will be stored as an ordered list of vertices P ij in a set. Additionally, the vertex numbers of particle i and j that make up the intersection polygon are stored in two separate sets, V i and V j respectively.
Each time step all polygon intersections are determined and stored in a tree-like structure with the intersection stored in a set

with a link to the previous intersection(s). A polygon intersection has a previous intersection if
or textually described: a polygon intersection has a previous intersection if the intersection of the both vertex number sets of the current time step with their counterpart from the previous time step are not empty for any intersection polygons between i and j of the previous time step. Links to previous intersections are made for all cases where the intersection set is not empty. As with Eqs. (33), (34) does not only hold for two contact point, but also holds for an arbitrary number of contact points between particles. This arrangement leads to a tree-like structure as shown in Fig. 11, which normally due to internal ordering can be traversed very efficiently for look-up queries. Since only the previous time step has to be kept in memory, the previous of the previous time step can be deleted every iteration to preserve memory. For a polygon intersection several arrangements are now possible.

New Intersection: If the parent (current polygon intersection) has no children (previous linked polygon intersection) the intersection is new.
Continuating Intersection: If the parent only has a single child and the child only has a single parent the polygon intersection is a normal continuation from the previous time step.
Merging Intersection: If the parent has two or more children, previous polygon intersections have merged into a single one and thus it is a merging intersection.
Splitting Intersection: If the parent has a single child, but the child has more than one parent the previous polygon intersection has split into two or more intersections and thus it is a splitting intersection.
Since the detection of the type of intersection does only rely on references between intersections and not on any information about the intersection itself this method can be used for any shape regardless of implementation.

Stability and reliability
The developed model is shown to be stable and reliable. Although the implementation of the model is not completely deterministic due to the use of sets and maps, in practice it behaves as such. This gives rise to highly reliable and reproducible results. During a simulation of particle packing, a relative smooth damping of the kinetic energy can be seen as shown in Fig. 12, after the initial peak of the conversion from potential energy to kinetic energy. This smooth transition shows that the simulation is stable with no sudden addition of energy if the time step is chosen appropriately. The apparent noise in the low energy region is only due to a few particles still settling.

Simulation parameters
The input parameters of the simulations are chosen to behave sugarlike. The shape of the particles is defined with a shape fitting algorithm and the initial particle configuration is based on a grid or random close packing algorithm.

Particle shape
The parameters of the simulation are based on sugar particles of Südsucker S1 crystalline sugar. Sugar was sieved with a Retsch AS 200 sieve shaker with a sieve stack from 400 to 100 μm with a 100 μm interval. Only the sieve fraction that passed the 400 μm sieve but not the 300 μm one was used. The particle size distribution was measured with a Fig. 11. The four possible different cases shown for the polygon intersections between two particles. The top two rows show the particle arrangement with the top being the last time step. The bottom two rows show a schematic of the tree data structure. In the c++ code each time step has its own memory map that is implemented as a multi-map inside a map where the keys are particle references i and j and the value is the intersection data. Malvern mastersizer 2000 particle size analyzer. The particle distribution of all sieve fractions is shown in Fig. 13. The shape of the particles was determined with optical microscopy in conjunction with a shape fitting algorithm. A representative image of the 300-400 sugar sieve fraction is shown in Fig. 14 at a magnification of 5×, imaged with a Carl Zeiss Stemi 2000-C Stereo Microscope. The shape of the particles is defined with a list of coordinates arranged in a counter-clockwise direction with as origin the center of mass of the particle based on the binary image. In pre-processing steps the image is made binary with global thresholding, all border particles and particles touching other particles are removed, holes are filled and small fines are removed. The shape fitting algorithm will fit the image of a particle with a given number of vertices by minimizing the value of the heuristic function. The heuristic function used is the absolute ratio of the difference in area between the image of the sugar particle and the fitted polygon area described by the list of coordinates over the total image area. An additional constraint is added penalizing angles sharper than 22.5°for stability reasons. The value that is minimized calculated by the heuristic function is, where A fit is the area of the fitted polygon, A sugar is the area of the binary image of the particle, A im is the total image area, defined as the product of the image height and image width, and v θ is a value for penalizing sharp angles. A sugar is defined as the summation of the pixels of the binary image created with adaptive thresholding where all holes in the blob (Binary Large OBject) are filled. The image height and width are defined as the minimum size needed to encompass the blob. The value v θ is used for penalizing angles sharper than a set minimum angle between three vertices. The angle constraint is necessary to avoid that the algorithm collapses several subsequent points onto a single straight line in a non-sequential order. Above the minimum angle, v θ will be zero, for angles of zero radians v θ will be one, any value in between is linearly interpolated. The minimum angle value of 22.5°is dependent on the particle shape and number of fitted vertices and chosen pragmatically, based on a series of trial runs. Due to the imaging process a bias is introduced in the particle size distribution since particles will tend to rest on the largest side. An example of the results of fitting a sugar particle is shown in Fig. 15, where the fit is made for three numbers of vertices. As can be seen with a fit of 16 vertices the shape of the particle is welldescribed. What should be noted is that due to the particle imaging method a bias is introduced. Particles have a tendency to rest in the most stable position which in the method means on their largest flat side. This results in viewing particles only in this orientation and thus skews the particle size distribution towards larger particles.

Particle stiffness
For determining the spring stiffness coefficients for the initial loading and unloading/reloading of sugar, nano-indentation was performed Fig. 13. The particle size distribution of unsieved S1 sugar and of the three different sieve fractions. The discrepancy between the measured mass-median diameter D 50 value and the sieve size is because both techniques use a different particle diameter. Fig. 14. A color inverted image of a selection of sugar crystals from the 300-400 sieve fraction of S1 Südzucker sugar.  on single sugar crystals. Nano indentation was performed on an MTS nano indenter xp with a 10 μm diameter diamond spherical tip. The best indentation result is shown in Fig. 16, where no slippage or breakage of the sugar crystal occurred. The conventional load-displacement curve [29] is adjusted to use an intersection area instead of an indentation depth, in order to be used in the simulation. The calculated coefficients, k l and k r , are 5.8 × 10 6 and 9.2 × 10 6 for loading and un/re-loading respectively based on a body revolution based intersection area.

Initial particle configuration and parameters
As base input for a simulation a set of particles is analyzed for the shape description of the particles. For each simulation, this set is used multiple times to generate the required number of particles. For a new simulation all particles are placed in a grid or with a random close packing algorithm for circular particles [30] and given a small random velocity to break any specific order introduced by the algorithm. All particle parameters can be set individual to mimic the inhomogeneity of most food powders but in this case are kept all equal. The values for the parameters are shown in Table 1. The parameters are chosen such that the particles behave sugar-like. Further calibration is required for direct experimental comparison but that is outside the scope of this paper.

Results
Simulations were performed to validate the model against literature based on ellipse packing properties, to quantify the severity of the merging and splitting of the intersection areas, and to qualify the influence of particle shape on multi-material 3D food printing processes.

Validation based on ellipse packing properties
The packing density of granular systems is of interest in almost all powder handling. As such, there is a large body of literature regarding packing densities. Packing densities of ellipses of different aspect ratios and contact friction parameters have been described by Guises et al. [31] and Dong et al. [32] against which the model described in this paper is validated. The paper from Guises et al. [31] models the packing fraction based on a combined FEM/DEM model, while the paper from Dong et al. [32] models the packing fraction based on a DEM model implementing an ODDS (Orientation Discretization Database Solution).
A set of simulations is performed where 288 mono-disperse elliptical particles are allowed to settle in a rectangular container with a width of 15 mm and a height of 20 mm. The elliptical particles are described as polygons each consisting of 20 vertices. The aspect ratio α of the particles is varied from one to five, where α = a/b, with a and b as the major and minor axes of the ellipse respectively. The areas of the ellipses are kept constant at 0.785 mm 2 irrespective of aspect ratio, for this the axes are rescaled with a ¼ ffiffiffi ffi . Friction parameters are varied to create a frictionless system, where μ s = μ d = 0 and one with friction where μ s = μ d = 0.5. The wall friction is kept equal to the particle-particle friction. The ellipses are generated in a loose grid pattern 10 mm above the bottom of the container with enough random initial velocity to destroy any order due to the initial configuration. The final configuration after the particles have completely settled is shown in Fig. 17.
As expected, frictionless particles align themselves, locally creating crystalline structures, whereas more disorder is observed with the higher friction coefficient. The packing fraction for all cases are determined by dividing the actual filled area over the total area. The total area used is a rectangle of 13 × 14 mm 2 (width × height), situated 1 mm from the sides and 1 mm from the bottom, this to avoid any Table 1 All the particle-shape independent parameters that are set at particle level. influences from the walls on the calculated packing fraction. The calculated packing fractions with respect to the ellipse aspect ratio and the values reported in literature are shown in Fig. 18. As can be seen, the model captures both the magnitude as well as the general trend of the packing fraction relatively well compared to literature. The trend shows that the packing fraction is a complex interaction between particle shape and friction parameters. The frictionless case has a very good correlation with the simulation result of Guises with only the extremes showing deviation. The frictionless case of Dong is included for completeness, but as noted in his paper: "Note that in our simulations we avoid to set μ s = 0, as the contact force model is established under the normal condition and will probably be invalid under such an extreme condition" [32](p. 507). For the friction case with μ s = 0.5 there is a good correlation within the error margin of the simulations. The reported optimum packing fraction at an aspect ratio of 1.6 is also captured within the error margin. As shown, the model is in good agreement with values reported in literature.

Merging and splitting of intersection areas
To verify the correct behavior of the model for merging and splitting intersection areas a simulation of 2,336 sugar particles each consisting of 20 vertices deposited in a closed container was simulated in 20,000 time steps with a resolution of 5 μs. As base input for the simulation 73 sugar particles of the 300-400 sieve fraction were used for the shape description of the particles. For the initial configuration a grid pattern was used for particle placement. Three time steps of the simulation are shown in Fig. 19, a video of the deposition can be found in the supplementary material. What can be seen in the enlargement of the last time step is that already some ordering takes place due to the particle shape. Flat sides have a tendency to align with neighboring flat sides and thus create force chains. What can also be seen is that a lot of multiple contact points exist between two particles.
With the existence of a large amount of multiple contact points the occurrence of merging and splitting of intersection areas should be taken into account as stated earlier. To quantify the effect a set of five simulations were done similar to the one shown in Fig. 19, with the particle shape fitted with 4, 8, 12, 16 and 20 vertices. The location, initial velocity and all other parameters were the same between the simulations with different vertex numbers. The average coordination number, number of particles in contact, of the particle packing after it settled around 0.07 s for the different shape fitting is shown in Fig. 20, along with two enlarged sections of the particle packings for the shape fitted with 4 and 20 vertices. The small peak and subsequent valley in the coordination number is caused by the last layer of particles hitting the underlying layer and the whole system bouncing back up a bit. During the deposition and settling of the packing new contact points will originate giving an indication how many new force trajectories, akin to the one shown in Fig. 3, will be created. Within the force trajectory a split of the intersection area can occur as shown in Fig. 8. The cumulative number of the creation of new contacts and splitting intersection areas is shown in Fig. 21. For quantifying the effect of a splitting event the ratio of splitting intersection areas over newly created contact points at a specific time step is taken. This metric gives an indication what the probability is that a force trajectory will go through this event. This is the probability that there is a sudden increase in kinetic energy during collision for noncorrected hysteretic models. The probability for this event occurrence with respect to the number of vertices fitted at time 0.07 s is shown in Table 2. As expected shapes fitted with four vertices have no splitting events since the shape is convex and no splitting event can occur. The higher the number of fitted vertices the higher the probability of a splitting event to occur although after 12 vertices the increase levels off. This is in accordance with expectations, after the rough shape is fitted adding more detail by adding more vertices will not result in significant different collision behavior and thus also not in different splitting behavior. This is specifically true in the case of sugar particles which have a relative cubic shape and can be fitted well with a small number of vertices.

Particle shape influence
The effect of simulating particles with their actual shape is easily seen in simulation. Two cases are shown here that simulate processes during multi-material 3D food printing, namely deposition through a miniaturized hopper and powder-bed compaction. Loading of a miniaturized hopper with sugar from the sieve fraction 300-400 μm with a nozzle diameter of 2 mm will block the flow of powder after some initial flow. The loading is simulated by placing the end result of the simulation shown in Fig. 19 above a hopper and removing the bottom plate. When simulated, this effect is reproduced as shown in Fig. 22, a video of the hopper loading supplied in the supplementary material. What can be seen is that the forces near the nozzle form a bit of an arch-like pattern but that the particles do not form a very clear arch pattern. What should also be noted is that force lines are developed where particles align with their neighbors.
The second case is powder-bed compaction. The creation of a powder bed is simulated by the deposition of 200 arbitrarily shaped particles in a square container. The generated shapes simulate the effect of particles with high concavity. Each particle has a unique set of parameters drawn from a Gaussian distribution to mimic the inhomogeneity of food powders; the parameters loosely resemble crystalline sugar. Particles are placed in two groups tightly packed with a random close packing algorithm for circles, with no initial overlap between particles. To guarantee no influence from the initial configuration each particle is   19. Three time steps of the deposition of sugar-like material from the 300-400 sieve fraction, with the particle shape fitted with 20 vertices into a 30 mm wide box. The bottom row shows an enlargement of a section of particles for better viewing of the particle shape and interactions. Colors indicate the initial height of the particle at the start of the simulation. given a random velocity vector. The velocity is chosen such that in the first 10% of the simulations on average the particles move a complete particle diameter in a random direction. The deposition is shown in Fig. 23. The effect of particle shape can be readily seen by comparison of the deposition of different particle shapes and comparing the properties of the final powder bed. A comparison is made between arbitrarily shaped polygonal, square and circular particles. The circular particles are modeled as polygons with 20 vertices but this does not have an appreciable effect on the behavior [33]. The area of the particles are matched, all other parameters (velocity, starting location, material parameters) are kept exactly the same. The starting state and final deposition after 0.2 s is shown in Fig. 24. The particle packing density calculated as percentage is 91.1, 93.9 and 90.0% for circular, square and polygonal particles respectively. The packing density is calculated in the area below the dashed line shown in Fig. 24. As can be seen the arbitrarily shaped polygons are packed less dense which is not surprising since this system is more likely to form voids. The square-shaped particles are packed slightly denser than the circular particles. This is caused by the tendency of particle sides to align with each other resulting in a decrease of voids.

Conclusion
A DEM model for arbitrarily shaped particles has been described where the forces are linked to the intersection area instead of the penetration depth. It has been shown that force models need to take into account the possibility that two particles can have multiple contact points between them. If this possibility is not taken into account discontinuities in the force will be present if a single contact area splits into two or more areas. Furthermore, it is shown that velocity-based damping is not appropriate for DEM models based on the intersection area.
A force model is presented based on the intersection area where contact areas can split or merge without a discontinuity in the overall force on the particle. An improvement to the linear hysteretic force model is made, where upon the splitting of a contact area the maximum area for each separate contact is recalculated as a weighted proportion of the non-split area. Furthermore, it has been shown that for an area based force model for arbitrarily shaped particles, damping should relate to the rate of change of the intersection area.

Table 2
The fractions of the cumulative number of intersection areas that underwent a splitting interaction over the number of newly formed contact points for shapes fitted with different number of vertices.  Colors only indicate an arbitrary particle number for distinction between particles. Fig. 24. The top three images show the starting condition of the system with all the same parameters with the exception of the shape. The bottom three images show the state after deposition has stopped. Colors only indicate an arbitrary particle number for distinction between particles.
Our proposed model is in good agreement with in literature reported packing ratios for ellipses of different aspect ratios. Simulated deposition of sugar shaped particles shows that for a dense packing the fraction of splitting interaction that occur per each force trajectory can be up to 7%. The influence of particle shape is shown to be significant due to particle alignment that cause shape dependent properties to emerge.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.