Calibration of DEM for Cohesive Particles in the SLS Powder Spreading Process

In this paper, a new DEM calibration procedure based on two different types of procedures to compare simulation with experiments is proposed. The aim is to find the values of the interfacial adhesive surface energy and the coefficient of rolling friction between the particles to be used in the simulation. The approach adopted is the so-called Bulk Calibration method. The experimental values of the angle of repose and unconfined yield strength, found with a static testing method and by shear testing, respectively, are compared, respectively, with the angle of repose, found in a simulation reproducing the experimental procedure, and the unconfined yield strength, obtained from an idealized uniaxial testing procedure. The simulated DEM particles are spheres equipped with the Hertz Mindlin with JKR contact model. The results suggest that a bulk calibration approach is not able to provide results that are consistent with two simple bulk property evaluations and, therefore, direct ways to estimate the surface energy based on the evaluation of interparticle forces, for example, should preferably be adopted.


Introduction
Interest in additive manufacturing (AM) is growing in many applications, among which are biomedical, aerospace, buildings [1], due to its capability to create complex parts directly from raw materials with minimal waste and the possibility to use many types of materials such as metals, ceramics, polymers, composites [2]. Moreover, powders can be reused, making AM more affordable compared with traditional manufacturing methods [3], while other advantages are the time and cost reduction it provides [4]. Among the AM processes, Selective Laser Sintering (SLS) is a process by which artifacts are created by alternating spreading and laser sintering of layers of particles, done on the basis of a CAD file [5]. With respect to other manufacturing processes, SLS has the great advantage of being relatively inexpensive, but it has also some disadvantages, such as of being a slow process [2] and of being able to produce artifacts with size limitations, poor surface finishing [6] and low mechanical properties [2]. In particular, the latter two disadvantages are mostly related to a bad distribution of the powder in the formation step of the powder layer. Hence, the importance of focusing on the powder distribution process during the layer formation simulated both experimentally and numerically.
To understand why some irregularities of the layer are obtained, simulations using the discrete element method (DEM) can be a powerful tool [7]. DEM is a computational method by which the granular flow is simulated at the individual particle level. Often, a soft sphere approach [8] is used, in which the particles are modeled as discrete entities and their motions and trajectories are described by their translational and rotational motion, governed by the Newton's second law [9]. To calculate the forces between the particles, as well as between the particles and the objects around the particles, DEM requires appropriate contact models. Many of these models are reported in the literature, and differ for the filled with bulk material. Karkala et al. [28] simulated both the dynamic yield strength test and shear cell experiment to calibrate the following DEM parameters: surface energy, coefficients of sliding and rolling friction. Coetzee [29] calibrated the contact stiffness using a confined compression test and the coefficients of sliding and rolling friction using either the individual or combined results from a draw down, rotating drum and annular ring shear test. Chen et al. [30] evaluated the sliding friction coefficient and the interfacial adhesive surface energy between the particles by comparing the macroscopic flow profile of the powder during layering between the numerical and experimental results. Angus et al. [31] combined sliding test experiments and DEM simulations of both homogeneous simple shear and the FT4 shear cell to evaluate the coefficients of rolling and sliding friction; they figured out that one single simulation type could not be sufficient to calibrate the DEM parameters. In fact, homogeneous simple shear simulations returned multiple pairs of both the coefficients. Simulations of the full FT4 shear cell with different pairs obtained with the previous simulation procedure were necessary to find the correct pair.
In order to elucidate the powder spreading process in SLS, a new experimental apparatus has been designed and built at the University of Salerno. It allows for the reproduction of different spreading conditions by changing the process conditions determined by the kind of the spreading tool and its movement parameters, such as the gap below the spreading tool, its shape, and in case of rollers the rotational velocity. Therefore, the tool is able to experimentally simulate the spreading of the powder in different configurations of the SLS apparatus, in order to understand the mechanisms that determine the formation of a layer and optimize the process conditions; however, the proper modeling of the process is necessary to gain an insight into the process.
In this work, the calibration of a DEM to be used in future for the modeling of SLS spreading is carried out. The contact model taken into consideration is the Hertz-Mindlin and JKR models to describe the contact mechanics. The use of these models is supported by the fact that the particles in the SLS spreading process do not experience high stresses and, therefore, the intrinsic limit of this model not accounting for plastic deformations at contact points may not be significant. The bulk calibration approach is used to evaluate the interfacial adhesive surface energy and the coefficient of rolling friction between the particles, while the main objective of the calibration procedure will be to find a value of the interfacial adhesive surface energy and of the coefficient of rolling friction between the particles. Two different independent calibration procedures to estimate these two parameters are proposed: one is based on the match of the static angle of repose estimated with the formation of a powder heap, as obtained on both experiments and model simulations; the other relies on the comparison between the unconfined yield strength experimentally measured with a shear testing procedure and the unconfined yield strength calculated from an idealized uniaxial testing model simulation. Each of the two calibration methods provides an independent set (couples) of values of the model interfacial adhesive surface energy and of the coefficient of rolling friction between the particles able to satisfactorily describe the experiments with the model. The scope will be to verify whether it is possible to use these two independent sets to find at least a single couple of values of the model surface energy and of the coefficient of rolling friction between the particles to be used in both model simulation procedures to satisfactorily describe all the experimental results.

Materials
The polymeric powder used in the experiments to be simulated with the DEM is Domo's Sinterline unfilled polyamide 6 3400 HT110 Natural (PA6), purposely designed for the SLS process. The material properties derived from the literature or from the technical data sheet of the powder, and used in the DEM, are reported in Table 1. Table 1. Material properties of PA6 powder to be used in the simulation: d 10 , d 50 and d 90 are the particle diameters corresponding to the 10th, 50th and 90th percentile of the cumulative particle size, respectively; ρ s is the density of the solid; ρ b is the bulk density of the powder. Value with reference was taken from the literature, the other values were taken from the technical data sheet of the powder. To characterize the shape distributions of the particles to be reproduced with the DEM, SEM images have been taken by using the TM3030Plus Tabletop Microscope (Hitachi, JP). The images are reported in Figure 1. They show that the particles are clearly non-spherical and inhomogeneous. for the SLS process. The material properties derived from the literature or from the technical data sheet of the powder, and used in the DEM, are reported in Table 1. To characterize the shape distributions of the particles to be reproduced with the DEM, SEM images have been taken by using the TM3030Plus Tabletop Microscope (Hitachi, JP). The images are reported in Figure 1. They show that the particles are clearly nonspherical and inhomogeneous.

Direct Measurement of Particle Properties to Be Used in DEM
Some properties, which must be known to properly simulate the powder, were taken from the technical data sheet of the powder, such as the Young's modulus of the particle Processes 2021, 9, 1715 5 of 20 material, or were found in the literature such as the Poisson's ratio of the particle material and the coefficient of restitution between the particles, and these are reported in Table 2. A procedure similar to that described by Nan et al. [1] was followed by Lupo [33] to evaluate the coefficient of sliding friction between the particles, the coefficient of sliding friction with an aluminium surface, and the heap angle. Both the coefficients of sliding friction are reported in Table 2, and the heap angle was 36.96 • ± 2.5 • . Table 2. Properties inserted in the EDEM ™ model for the interaction between the particles (PA6) and between the particles and the geometries made of aluminium: e P−P and e P−G are the coefficients of restitution between the particles and between the particles and the geometry, respectively; µ S, P−P and µ S, P−G are the coefficient of sliding friction between the particles and between the particles and the geometry, respectively; µ R, P−P and µ R, P−G are the coefficient of rolling friction between the particles and between the particles and the geometry, respectively; ν is the Poisson's ratio; E is the Young's modulus. Values with references were taken from the literature, while the others were taken from the technical data sheet of the powder or directly measured as described in the text. Regarding µ S, P−G there are two rows: the first reports the values of µ S, P−G used in the unconfined yield strength (UYS) calibration; the second reports the value of µ S, P−G used in the static angle of repose (AOR) calibration, and which has been measured as described in the text.

Material
PA6 Aluminium e P−P (-) A Schulze ring shear tester RST-01.01, equipped with an S-cell (internal volume of the cell of 203.58 cm 3 ; external and internal annulus diameter of the lid of 118 and 62 mm, respectively) was used to shear test the powder and to derive the unconfined yield strength in condition of absence of consolidation, i.e., the major principal stress, σ 1 , equal to zero. Both the standard procedure recommended for the ring shear tester RST-01.01 and the ASTM-standard (ASTM D6773, 2002) were followed. Before the shear tests, the powder has been placed in the oven for 1 h at 100 • C to remove the humidity. Tests at consolidation stresses between 50 and 90 Pa were carried out and 4-5 shear points were registered to obtain the corresponding yield locus. A software application, developed in the LabVIEW environment (National Instruments), was used to acquire, visualize, and record the main data, which were measured during the shear experiments. Moreover, the RSV 95 software version 2.1.0.1 for the RST-01.01 Schulze shear tester was used to derive the unconfined yield strength of the powder from the powder yield loci.

DEM
DEM is a numerical procedure which allows to describe the mechanical behavior of assemblies of particles [9]. There are two main approaches that can be used to calculate contact forces: the hard-sphere approach and the soft sphere approach [8]. In the first case, each binary collision is modeled as an instantaneous process, and the properties of the particles after the collision are related to the properties of the particles before the collision by means of momentum and energy balances. Therefore, this method can be used only for those systems where multiple and simultaneous collisions are unlikely to occur, such as dilute systems. Moreover, this method cannot include long-range interparticle forces in the Processes 2021, 9, 1715 6 of 20 model [37]. On the contrary, the soft sphere approach deals with multi-particle collisions and the interparticle forces can be easily implemented, even if with higher computational time. In the soft sphere approach, the particles are modelled as individual entities and their movement is described by the translational and rotational motion, given by the Newton's second law: where m i , v i and ω i are the mass, translational velocity and angular velocity of particle i, respectively; t is the time, F c ij is the contact force due to interaction between the particle i and the nearby ones j or the wall; F nc ik represents the non-contact forces like van der Waals and electrostatic forces acting on particle i by particle k or the other possible forces; F f i is the particle-fluid interaction force; F t ij R is the torque imposed on particle i by particle j, I i is the moment of inertia and g is the acceleration due to the gravity. By integrating (2) and (3), new velocities and positions of the particles are obtained [9,38].
To implement a DEM simulation, a contact model has to be assumed with the corresponding parameters. A contact model describes how elements behave when they come into contact with each other, i.e., describes the forces deriving from particle collisions and contacts. As a result, particle interactions are considered in the calculations. There exist a number of force models which mostly allow particles to have deformation. The deformation is modeled as an overlap between particles. Among the most used contact models are the Hertz-Mindlin (no slip) and the Hertz-Mindlin with JKR (Johnson-Kendall-Roberts) cohesive model. The Hertz-Mindlin (no slip) is the default model used in the EDEM software (provided by DEM Solutions Ltd., Edinburgh, UK). It accounts for the effects of contact of two elastic spheres. In this model: the normal force component is based on Hertzian contact theory [11], the tangential force model is based on the work of Mindlin-Deresiewicz [13,14]. Hertz-Mindlin, along with the JKR cohesive model, is a cohesion contact model that also accounts for the influence of van der Waals forces within the contact zone and allows the user to model strongly adhesive systems, e.g., dry powders or wet materials.
According to the JKR model [15], the normal elastic contact force depends on the overlap δ and the interaction parameter, i.e., the interfacial adhesive surface energy Γ, as follows: where a is the contact radius, while E * and R * are the equivalent Young's modulus and the equivalent radius, respectively, which are defined as: where E i , ν i , R i , and E j , ν j , R j , are the Young's modulus, the Poisson ratio and the radius of each sphere in contact, respectively. The interfacial adhesive surface energy, Γ, is defined as follows: Processes 2021, 9, 1715 7 of 20 where γ 1 and γ 2 are the surface energies of the two spheres and γ 1,2 is the interface surface energy. For the special case where two spheres of the same material come into contact, the interface surface energy is zero (γ 1,2 = 0), thus the interfacial surface energy becomes Γ = 2 γ, where γ 1 = γ 2 = γ. For Γ = 0, F JKR turns into a Hertz-Mindlin normal force. The EDEM ™ software has been used to simulate both the formation of the powder pile for the measurement of the static angle of repose and the uniaxial compression test of a powder bed, with the aim of calibrating the interfacial adhesive surface energy between the particles of the Hertz-Mindlin model with the JKR model. Both the EDEM ™ 2018 and the EDEM ™ 2019 software packages, provided by DEM Solutions, Edinburgh, UK, have been used.

Model Tuning for SLS Applications and Calibration Procedures
By using the Hertz-Mindlin model with JKR cohesive model, elastic deformation and adhesion between the contacting particles are taken into account. Therefore, it is assumed that no plastic deformation between the particles in the contact point is involved for the powder considered. A model that does not account for the plastic deformation at contact points would be inadequate to completely describe the effects of powder consolidation. In fact, consolidation is mainly affected by two contributions: one is the increased number of contacts between the particles, the other is the increased strength of the binary forces between particles due to local plasticization of the contact point and the consequent increase of the surface of contact between particles. Several experimental results on the powder flow properties measured at increasing temperature have proven that the effect of temperature in the change of powder flow properties can only be explained by introducing the plastic deformation of the contact that is affected by temperature, which changes the particle material yield strength [39,40]. A purely elastic model, as the one considered in this paper, can only rely on the first mechanism to describe powder consolidation and, in general, may not be adequate for real powders. Nevertheless, in SLS application, powder consolidation stresses are very low and, for this special case, interparticle forces might be adequately accounted for by a purely elastic model. This also means that in the model simulation, the value of compression force is only relevant to realize the correct number of interparticle contacts. The relatively independent value of the measured unconfined yield strength of the material with the consolidation stress supports the hypothesis of purely elastic deformation of the particles at the contact points. Since the number of these contacts is directly related to bed porosity, the compression step is carried at values of the compression stress that are able to provide values of the bed porosity similar to the experimentally found in SLS processes. In order to compare simulation results with experimental results obtained with the shear tester in the powder, the experimental result obtained with this test was extrapolated to zero consolidation, which is more relevant to the SLS process. The Hertz-Mindlin model has been used also for the interaction between particles and the geometries; this means that no adhesive forces between the contacting bodies have been considered.
Spherical particles have been used to simulate the powder rather than more complicated shapes in order to reduce the simulation time. From the SEM images ( Figure 1) it appears evident that the polymeric powder is constituted by particles with different and irregular shapes. Initially, in order to replicate the actual shape of the particles, it has been tried to model them as clumped spheres [41,42] that is each particle can be approximated by a number of overlapping spheres of different sizes. However, by using this approach, simulations were excessively long, even if the particles were modelled as two or three clumped spheres. This is due to two different reasons: the first is the increase of the particles involved to simulate the powder; the second reason is that particles, which were already small, were modeled with particles even smaller and, consequently, the simulation time greatly increased. In fact, the simulation time step is chosen as a fraction of the Rayleigh time step, i.e., the natural oscillation period of the mass-spring system made by the particle Processes 2021, 9, 1715 8 of 20 and the elastic contact [43]. The Rayleigh time (∆t R ) depends on the particle diameter according to the following relationship: where R is the particle radius, ρ p is the particle density, G is the shear modulus and ν is the Poisson's ratio of the particle. Therefore, the smaller the spheres used to model the particles, the slower the simulation.
The particle size distribution of the powder has been approximated to a purely Gaussian distribution, centered in d 50 = 55 µm, with a width 3 σ, such that: By doing so, particles smaller than d 10 and bigger than d 90 have been cut out. Indeed, if particles smaller than d 10 were considered, the time step would have been too small and the simulation time would have been too long. This approximation seems to be reasonable, as demonstrated ( Figure 2) by the shape of the frequency particle size distribution (PSD) and of the cumulative one, based on the equivalent circle particle diameter (CE), as obtained with the Morphologi G3 (Malvern Panalytical Ltd., Malvern, UK). CE is the diameter of a circle with the same area as the 2D image of the particle. by a number of overlapping spheres of different sizes. However, by using this approach, simulations were excessively long, even if the particles were modelled as two or three clumped spheres. This is due to two different reasons: the first is the increase of the particles involved to simulate the powder; the second reason is that particles, which were already small, were modeled with particles even smaller and, consequently, the simulation time greatly increased. In fact, the simulation time step is chosen as a fraction of the Rayleigh time step, i.e., the natural oscillation period of the mass-spring system made by the particle and the elastic contact [43]. The Rayleigh time (∆t ) depends on the particle diameter according to the following relationship: where R is the particle radius, ρ is the particle density, G is the shear modulus and ν is the Poisson's ratio of the particle. Therefore, the smaller the spheres used to model the particles, the slower the simulation.
The particle size distribution of the powder has been approximated to a purely Gaussian distribution, centered in d = 55 μm, with a width 3 σ, such that: By doing so, particles smaller than d and bigger than d have been cut out. Indeed, if particles smaller than d were considered, the time step would have been too small and the simulation time would have been too long. This approximation seems to be reasonable, as demonstrated ( Figure 2) by the shape of the frequency particle size distribution (PSD) and of the cumulative one, based on the equivalent circle particle diameter (CE), as obtained with the Morphologi G3 (Malvern Panalytical Ltd., Malvern, UK). CE is the diameter of a circle with the same area as the 2D image of the particle.  Simulation of the static angle of repose resembles the experimental procedure. Inside the hopper, a dynamic virtual factory creates particles in a random way (Figure 3). The particles fall through the hole of the hopper and then they deposit on a layer of particles of the same powder; indeed, this layer is created immediately before the fall of the particles. By doing so, the experimental situation of the angle of repose test where fresh powder deposits on a previously deposited powder is replicated, while when a clearly visible heap is obtained, the creation of the particles is interrupted within the factory. The angle of repose is measured once the particles are all steady, while the hopper is a truncated cone. The radii of the bases are 100 and 400 µm, while the height of the cone is 1 mm. The properties of the powder and of the geometries, which have been used, are those reported in Tables 1 and 2. Processes 2021, 9, x FOR PEER REVIEW 9 of 21 Figure 2. (a) Frequency particle size distribution (PSD) and (b) cumulative PSD on a volume basis as function of the equivalent circle particle diameter (CE).

DEM Calibration Using the Static Angle of Repose
Simulation of the static angle of repose resembles the experimental procedure. Inside the hopper, a dynamic virtual factory creates particles in a random way (Figure 3). The particles fall through the hole of the hopper and then they deposit on a layer of particles of the same powder; indeed, this layer is created immediately before the fall of the particles. By doing so, the experimental situation of the angle of repose test where fresh powder deposits on a previously deposited powder is replicated, while when a clearly visible heap is obtained, the creation of the particles is interrupted within the factory. The angle of repose is measured once the particles are all steady, while the hopper is a truncated cone. The radii of the bases are 100 and 400 μm, while the height of the cone is 1 mm. The properties of the powder and of the geometries, which have been used, are those reported in Tables 1 and 2. In the case of the simulation of the angle of repose, 40 CPUs have been used and the simulation time is approximately 8 h.

DEM Calibration Using the Unconfined Yield Strength
There are two methods to measure the unconfined yield strength of a cohesive powder subject to compression [44]. The simplest is to use a uniaxial tester in which the powder is inserted in a cylindrical mold. The powder is compressed with a piston lid to a defined consolidation stress, then when the compression force is released, the mold is open in order to free the agglomerated powder sample. If the compression force is sufficiently high, the compressed powder sample is agglomerated and is able to keep the shape of the cylindrical mold. The piston lid is then lowered again and the force necessary to break the cylindrical agglomerate is measured. The ratio between this force and the agglomerate cross section is the unconfined yield strength of the material, which is a function of the normal stress in the compression step. Unfortunately, in this procedure, it is very difficult to avoid the effect of the wall shear stress on the sample, reducing the effective compression force at the bottom of the sample [45] unless complex devices are used. In the case of the simulation of the angle of repose, 40 CPUs have been used and the simulation time is approximately 8 h.

DEM Calibration Using the Unconfined Yield Strength
There are two methods to measure the unconfined yield strength of a cohesive powder subject to compression [44]. The simplest is to use a uniaxial tester in which the powder is inserted in a cylindrical mold. The powder is compressed with a piston lid to a defined consolidation stress, then when the compression force is released, the mold is open in order to free the agglomerated powder sample. If the compression force is sufficiently high, the compressed powder sample is agglomerated and is able to keep the shape of the cylindrical mold. The piston lid is then lowered again and the force necessary to break the cylindrical agglomerate is measured. The ratio between this force and the agglomerate cross section is the unconfined yield strength of the material, which is a function of the normal stress in the compression step. Unfortunately, in this procedure, it is very difficult to avoid the effect of the wall shear stress on the sample, reducing the effective compression force at the bottom of the sample [45] unless complex devices are used. For this reason, a second method is more often followed, and the unconfined yield strength is indirectly measured from shear testing experiments [46]. It is assumed that the major principal stress obtained from the Mohr circle describing the stress state during shear testing consolidation coincides with the consolidation stress in the uniaxial test and that the yield condition in the uniaxial tester is represented by the Mohr circle passing through the origin of the τ − σ plane and tangent to the static yield locus obtained with the shear tester at the same consolidation condition. This latter procedure was used to calculate the flow function of the material, i.e., the unconfined yield strength as a function of the major principal stress at consolidation. Simulation of the shear testing procedure is quite cumbersome with the DEM procedure, due to the length of the test, the number of particles to be used and complexity of the shear tester geometry. Conversely, by using a virtual experiment in DEM it is possible to easily eliminate wall friction, and therefore, an idealized uniaxial testing procedure was used to calibrate unknown DEM parameters and compare the obtained unconfined yield strength, f c (sim), with the value indirectly measured in shear testing experiments, f c (exp), (Figure 4). For this reason, a second method is more often followed, and the unconfined yield strength is indirectly measured from shear testing experiments [46]. It is assumed that the major principal stress obtained from the Mohr circle describing the stress state during shear testing consolidation coincides with the consolidation stress in the uniaxial test and that the yield condition in the uniaxial tester is represented by the Mohr circle passing through the origin of the − plane and tangent to the static yield locus obtained with the shear tester at the same consolidation condition. This latter procedure was used to calculate the flow function of the material, i.e., the unconfined yield strength as a function of the major principal stress at consolidation. Simulation of the shear testing procedure is quite cumbersome with the DEM procedure, due to the length of the test, the number of particles to be used and complexity of the shear tester geometry. Conversely, by using a virtual experiment in DEM it is possible to easily eliminate wall friction, and therefore, an idealized uniaxial testing procedure was used to calibrate unknown DEM parameters and compare the obtained unconfined yield strength, fc(sim), with the value indirectly measured in shear testing experiments, fc(exp), (Figure 4).   The geometry is the one represented in Figure 5a. The mold is made by a flat base and two hemicylindrical walls that are kept fixed on the base during the compression step and removed during the yield step. In order to eliminate wall shear stress, the coefficient of sliding friction and the coefficient of rolling friction between the particles and the wall are put to a very low value (10 −4 ) in the compression step. This idealized approach makes the size of the sample less significant. Thus, in order to minimize the number of particles used in the experiment, the diameter of the mold was set to 0.5 mm and the mold height to around 1.5 mm, that is three times the mold diameter. The simulation consists of four steps ( of sliding friction and the coefficient of rolling friction between the particles and the wall are put to a very low value (10 −4 ) in the compression step. This idealized approach makes the size of the sample less significant. Thus, in order to minimize the number of particles used in the experiment, the diameter of the mold was set to 0.5 mm and the mold height to around 1.5 mm, that is three times the mold diameter. The simulation consists of four steps ( Figure 5): (a) generation of the powder bed inside the mold; (b) compaction of the bed with a circular piston lid to achieve the desired consolidation stress in correspondence of a powder bed porosity of about 0.55; (c) raising of the lid to zero out the force exerted by the lid on the powder bed; (d) unconfined compression of the powder bed by means of the lid. The powder bed has been created by a static factory within the mold (Figure 5a). This factory allows to create particles with a random arrangement instantly. The total mass of the powder bed used is 1.4 × 10 −7 kg. Properties of the powder or regarding the interaction between the particles and the geometries, which have been used in the simulation, are reported in Tables 1 and 2.
In order to avoid the effect of powder weight on the stress acting at the bottom of the sample, simulations have been carried out by setting the acceleration due to the gravity equal to 0 m/s 2 .
The general rule for the implementation of the model is that the time step should be between 10 and 40% of critical Rayleigh time step [43], as suggested in the EDEM ™ 2019.2 documentation. In this work, time steps of 20 and 40% of the critical time step have been used to verify the reliability of the measured unconfined yield strength. Once the bed is created in the cylindrical factory, a piston lid (Figure 5b) with the same diameter of the powder bed compacts the bed going down at a scrolling velocity, V, of 0.005 ms −1 . To guarantee roughly the same initial void fraction of the polymeric powder in all the simulations and to evaluate the unconfined yield strength in correspondence of a major principal stress close to 0 Pa, the lid is stopped when the consolidation stress obtained during the compaction is almost the lowest major principal stress that the powder has experienced during the shear tests, i.e., 988 Pa. At this point, the lid is moved up to its initial position at the same velocity as the compaction step (Figure 5c). During this stage, the bed relaxes because of the elastic behavior of the particles. Then, the lid is moved down again at V = 0.001 ms −1 (Figure 5d) and the hemicylindrical walls of the cylinder are removed sideways in opposite directions at V = 0.05 ms −1 . The absence of adhesive forces between the particle and the wall ensures that this movement of the mold does not produce any traction stress on the powder that would affect its state of consolidation. During this step, the coefficient of rolling friction and of sliding friction between the particles and the lid The powder bed has been created by a static factory within the mold (Figure 5a). This factory allows to create particles with a random arrangement instantly. The total mass of the powder bed used is 1.4 × 10 −7 kg. Properties of the powder or regarding the interaction between the particles and the geometries, which have been used in the simulation, are reported in Tables 1 and 2.
In order to avoid the effect of powder weight on the stress acting at the bottom of the sample, simulations have been carried out by setting the acceleration due to the gravity equal to 0 m/s 2 .
The general rule for the implementation of the model is that the time step should be between 10 and 40% of critical Rayleigh time step [43], as suggested in the EDEM ™ 2019.2 documentation. In this work, time steps of 20 and 40% of the critical time step have been used to verify the reliability of the measured unconfined yield strength. Once the bed is created in the cylindrical factory, a piston lid (Figure 5b) with the same diameter of the powder bed compacts the bed going down at a scrolling velocity, V, of 0.005 ms −1 . To guarantee roughly the same initial void fraction of the polymeric powder in all the simulations and to evaluate the unconfined yield strength in correspondence of a major principal stress close to 0 Pa, the lid is stopped when the consolidation stress obtained during the compaction is almost the lowest major principal stress that the powder has experienced during the shear tests, i.e., 988 Pa. At this point, the lid is moved up to its initial position at the same velocity as the compaction step (Figure 5c). During this stage, the bed relaxes because of the elastic behavior of the particles. Then, the lid is moved down again at V = 0.001 ms −1 (Figure 5d) and the hemicylindrical walls of the cylinder are removed sideways in opposite directions at V = 0.05 ms −1 . The absence of adhesive forces between the particle and the wall ensures that this movement of the mold does not produce any traction stress on the powder that would affect its state of consolidation. During this step, the coefficient of rolling friction and of sliding friction between the particles and the lid and between the particles and the base is increased from 1 × 10 −4 to 1 to avoid particles from sliding on the wall of the geometries during the compression step. In fact, in previous simulations, where both the coefficients were maintained at low values, i.e., 1 × 10 −4 , the whole powder bed slipped away during the compression, preventing its deformation.
For the simulation, four CPUs have been used and the simulation time has been about 2 days.
The unconfined compression stage is carried out by imposing a downward movement to the lid and by registering, step by step, the total force exerted by the lid on the powder bed. In a real uniaxial experiment, the maximum of this force divided by the cross section of the bed should correspond to the unconfined yield strength. In fact, it corresponds with the fracture of the bed and the loss of its strength. However, in the numerical calculations, the force measured on the lid follows the general trend reported in Figure 6. The figure reports the force relative to the whole simulation in two panes. Figure 6a reports the force as a function of time, while Figure 6b reports the force as a function of the lid position. The first peak (point A) corresponds to the initial compaction of the bed before the unconsolidated compression. It appears that during the unconfined compression, the force does not reach a maximum value, but the lid continues to move for a long lid displacement without significant change; in fact, the force reaches a plateau value. The kinds of bed deformation observed are reported in Figure 7, and it appears clearly that in the simulations, the bed starts shearing without breaking (Figure 7b). Only after long time and at very large bed deformations does the force start increasing again and raises above the plateau value, due to the limited space remaining available for the bed compression. In the simulation, particles cannot fall due to the lack of gravity, but particles might, in principle, detach. This happens only at very low values of Γ (Figure 7a). In order to compare simulation results with the experiment, the yield of the powder is assumed to happen at the attainment of the plateau value of the compression force, that has been used to calculate the unconfined yield strength. In particular, an average value of the force among those obtained before the sudden increase has been used as the plateau value.   . Total normal force exerted by the lid on the powder bed as function of (a) time and (b) lid position. These graphs have been obtained for Γ = 1 × 10 −3 J/m 2 and μ , = 0.8 with a time step of 20% of the critical time step. The first peak (point A) corresponds to the initial compaction of the bed in the consolidation step, the black hyphenated line corresponds to the plateau value of the force assumed valid for the calculation of the unconfined yield strength. Measured values of the unconfined yield strength (f c ) as function of the coefficient of rolling friction (µ R, P−P ) provided a fixed value of Γ have been reported in a plot. Fitting these values of f c with a regression line, the couple of the parameters (µ R, P−P , Γ) which gives the experimental f c is the desired one.

Calibration According to the Static Angle of Repose
The static angle of repose has been measured with the setup and by following the procedure described in Section 2.1.3 using six repetitions was 36.96 • ± 2.5 • , where the deviation reported is the standard deviation of results. Different simulations have been carried out to reproduce the formation of the heap and to evaluate the corresponding static angle of repose by choosing a time step of 40% of critical Rayleigh time step. Indeed, this time step guarantees the same value of the angle, which can be obtained with a lower percentage value, e.g., 20%. Five different values of Γ have been chosen in the range between 10 −5 and 10 −2 J/m 2 . For each value of Γ, different values of the coefficient of rolling friction between the particles (µ R, P−P ) have been used ranging between 10 −5 and 0.8. The aim is to find the pairs of the parameters for which the experimental static angle of repose, α(exp), is matched by the one obtained by the simulations, α(sim), (Figure 8).
Results are reported in Figure 9, where the value of angle of repose resulting from simulation (α) is plotted as a function of the interfacial adhesive surface energy and the coefficient of rolling friction between the particles. Since the value of α is estimated with a visual procedure, the simulation result is subject to a certain variation represented by error bars. Points corresponding to the same Γ are plotted using the same symbol. In the figure, the experimental value of the angle of repose is reported as a continuous red horizontal line. The two parallel red hyphenated lines represent the standard deviation of the measured value obtained by experiment repetition. Figure 9a also reports interpolating polynomial for each set of data at similar Γ. Intersections of these polynomials with the horizontal line representing the experimental α help to find couples (µ R, P−P , Γ) that are able to predict the correct value of the angle of repose. To verify the reliability of the results, new simulations have been carried out by setting the parameters for the four couples (µ R, P−P , Γ) found with the procedure above. Resulting values of α are reported in Table 3, which shows that the procedure is fairly adequate. viation reported is the standard deviation of results. Different simulations have been carried out to reproduce the formation of the heap and to evaluate the corresponding static angle of repose by choosing a time step of 40% of critical Rayleigh time step. Indeed, this time step guarantees the same value of the angle, which can be obtained with a lower percentage value, e.g., 20%. Five different values of Γ have been chosen in the range between 10 −5 and 10 −2 J/m 2 . For each value of Γ, different values of the coefficient of rolling friction between the particles (μ , ) have been used ranging between 10 −5 and 0.8. The aim is to find the pairs of the parameters for which the experimental static angle of repose, α(exp), is matched by the one obtained by the simulations, α(sim), (Figure 8).
is fixed is fixed (exp)? Figure 8. Schematic of the static angle of repose calibration methodology in a flow chart. Γ and Γ P−G are the interfacial adhesive surface energy between the particles and between the particles and the geometry, respectively. Table 3. Comparison between the predicted α (α from the polynomial) and the one obtained from the simulations (α sim ) for the reported pairs of the parameters (Γ and µ R,P−P ) to verify the reliability of the prediction. σ is the standard deviation of α sim . polynomial for each set of data at similar Γ. Intersections of these polynomials with the horizontal line representing the experimental α help to find couples (μ , , Γ) that are able to predict the correct value of the angle of repose. To verify the reliability of the results, new simulations have been carried out by setting the parameters for the four couples (μ , , Γ) found with the procedure above. Resulting values of α are reported in Table 3, which shows that the procedure is fairly adequate.  Figure 10 reports the couple of the parameters (µ R, P−P , Γ) represented as red circles on a plane. The best interpolating function for points in this plot is a logarithmic function: To verify the reliability of the equation, three other pairs of the parameters from the curve in the same range of Γ and µ R, P−P, reported as blue rhumbs in Figure 10, have been chosen. By using these pairs, the value obtained of α from the simulation is very close to the experimental one (Table 4). This proves that all the pairs of the logarithmic curve are all possible candidates. Table 4. Pairs of parameters (Γ and µ R,P−P ) from the logarithmic curve for which the experimental α and that obtained from the simulation (α sim ) are matched. σ is the standard deviation of α from the simulation.

Γ (J/m 2 )
µ R,P−P (-) To verify the reliability of the equation, three other pairs of the parameters from the curve in the same range of Γ and μ , , reported as blue rhumbs in Figure 10, have been chosen. By using these pairs, the value obtained of α from the simulation is very close to the experimental one (Table 4). This proves that all the pairs of the logarithmic curve are all possible candidates.   In the case of Γ = 1 × 10 −2 J/m 2 , represented in Figure 9b, a regression line between the points representing α at µ R, P−P = 10 −4 and at µ R, P−P = 10 −5 has been considered to have a best estimation of the coefficient of rolling friction in case of Γ P−P = 1 × 10 −2 J/m 2 . From the line, a value of µ R, P−P = 5.92 × 0 −5 has been obtained. Repeating the simulation with this value of µ R, P−P , the resulting static angle of repose is 42.6 ± 1.6, which is somewhat higher than the experimental angle, compared to the other values of α obtained from the simulations with the other pairs of parameters of Tables 3 and 4. Therefore, it is better to consider values of Γ between 1 × 10 −5 and 5 × 10 −3 J/m 2 as possible candidates.

Calibration According to the Unconfined Yield Strength
Four different tests at different maximum loads (0.4, 0.5, 0.6, and 0.7 kg) have been carried out with the shear tester. To verify the reproducibility of the measurements, each test has been repeated three times. A value of f c has been measured for each value of σ 1 .
The average values of f c , obtained from the repetition, have been reported in Figure 11, where the flow factors (FF = σ 1 /f c ), which limit the Jenike flowability classes [47], are reported as hyphenated lines. Data show that at ambient temperature and for σ 1 > 1600 Pa the powder is free flowing (FF > 10), while for 975 < σ 1 < 1600 Pa, the polyamide 6 is easy flowing (4 < FF < 10). The good flowability of this polyamide is not surprising, being that the powder specifically engineered for the SLS process. Significantly, the variation of the unconfined yield strength is not very much affected by the compression stress. The unconfined yield strength has been extrapolated at σ 1 = 0 Pa from the regression line (red hyphenated line in Figure 11). Its value is about 141 Pa.
With the aim of finding new pairs, values of Γ and of the coefficient of rolling friction have been sought comparing the experimental values of the unconfined yield strength and the model results of uniaxial idealized test. Time steps of 20 and 40% of Rayleigh time step were adopted. In these simulations, Γ has been set between 10 −5 and 10 −2 J/m 2 and µ R, P−P has been changed between 0.1 and 1. Results are reported in Figure 12. The error bars correspond to different values obtained by measuring f c with a time step of 20% and 40%. The bars are quite short and it proves that the measured value does not change much if the time step is varied in the aforementioned range. In the tested range, it was possible to find only a single couple (µ R, P−P , Γ) for which the experimental f c is matched by f c from the simulation, and it is provided by the intersection of the regression line of simulation results corresponding to Γ = 1 × 10 −2 J/m 2 , with the red line corresponding to the experimental value of f c . The intersection occurs at µ R, P−P = 0.326 and, correspondingly, the f c obtained from the simulations is 141 ± 7 Pa.
carried out with the shear tester. To verify the reproducibility of the measurements, each test has been repeated three times. A value of f has been measured for each value of σ .
The average values of f , obtained from the repetition, have been reported in Figure 11, where the flow factors (FF = σ /f ), which limit the Jenike flowability classes [47], are reported as hyphenated lines. Data show that at ambient temperature and for σ > 1600 Pa the powder is free flowing (FF > 10), while for 975 < σ < 1600 Pa, the polyamide 6 is easy flowing (4 < FF < 10). The good flowability of this polyamide is not surprising, being that the powder specifically engineered for the SLS process. Significantly, the variation of the unconfined yield strength is not very much affected by the compression stress. The unconfined yield strength has been extrapolated at σ = 0 Pa from the regression line (red hyphenated line in Figure 11). Its value is about 141 Pa.  experimental value of f . The intersection occurs at μ , = 0.326 and, correspondingly, the f obtained from the simulations is 141 ± 7 Pa.
In the case of 10 −5 J/m 2 , for μ , < 0.8 it was not possible to find any simulation value for f because the particles separated before the unconfined compression due to the very low values of resulting forces keeping them together. This case is the one reported in Figure 7a. From Figure 12 it appears that the higher the Γ, the higher is the resulting f , fixed the value of μ , . In fact, the higher the value of Γ, the higher the van der Waals attractive force between the particles and, correspondingly, the higher f becomes. The comparison between Figures 9 and 12 shows that the only value of the interfacial adhesive surface energy for which both the experimental α and f are matched by the ones from the simulations is Γ = 1 × 10 −2 J/m 2 . However, this happens with different coefficients of the rolling friction namely with a value of μ , close to 5 × 10 −5 , calibrating results using the angle of repose simulation, and with a value of μ , = 0.326, in the case of the unconfined compression test simulation. Therefore, it would seem that using DEM as we have done, adopting the Hertz-Mindlin contact model with the JKR model for cohesion, it is not possible to calibrate the model so that it can contemporarily reproduce the correct values of the unconfined yield strength and of the static angle of repose. In the case of 10 −5 J/m 2 , for µ R, P−P < 0.8 it was not possible to find any simulation value for f c because the particles separated before the unconfined compression due to the very low values of resulting forces keeping them together. This case is the one reported in Figure 7a. From Figure 12 it appears that the higher the Γ, the higher is the resulting f c , fixed the value of µ R, P−P . In fact, the higher the value of Γ, the higher the van der Waals attractive force between the particles and, correspondingly, the higher f c becomes.
The comparison between Figures 9 and 12 shows that the only value of the interfacial adhesive surface energy for which both the experimental α and f c are matched by the ones from the simulations is Γ = 1 × 10 −2 J/m 2 . However, this happens with different coefficients of the rolling friction namely with a value of µ R, P−P close to 5 × 10 −5 , calibrating results using the angle of repose simulation, and with a value of µ R, P−P = 0.326, in the case of the unconfined compression test simulation. Therefore, it would seem that using DEM as we have done, adopting the Hertz-Mindlin contact model with the JKR model for cohesion, it is not possible to calibrate the model so that it can contemporarily reproduce the correct values of the unconfined yield strength and of the static angle of repose. This disappointing result may be related with the other circumstance that we have observed above, regarding the lack of an apparent yield with a maximum of the compression force in the uniaxial unconfined compression test, suggesting that the model used is not able to properly describe the powder shear mechanisms during the unconfined yield. It may be that alternative contact models, or the use of different particle shapes, should be tried to verify if the DEM model is able to describe with the same calibration parameters the unconfined yield and the correct static angle of repose.
A different method was used by Lupo [33], based on the approach developed by Ruggi et al., by extrapolating the tensile strength of the powder from a linearized yield locus and by calculating from this the theoretical pull off force the Rumpf equation. This approach provides a value for the surface energy of Γ = 2 × 10 −3 J/m 2 , for which a value of µ R, P−P = 1.9 × 10 −1 should be applied.

Conclusions
In this paper, a calibration procedure, based on the comparison between the results obtained from two different comparisons between simulation and experimental results, has been proposed to determine the interfacial adhesive surface energy between the particles of the Hertz-Mindlin model with the JKR model as well as the coefficient of rolling friction between the particles. These comparisons involve the evaluation of the static angle of repose from the formation of the powder heap and of the unconfined yield strength from the uniaxial compression test. The results of the simulations show that the two calibration procedures are not able to provide a single set of calibrated parameters. More direct ways to estimate the surface energy based on the evaluation of interparticle forces are suggested, rather than using a bulk calibration approach.