A discrete element framework for the numerical analysis of particle bed-based additive manufacturing processes

This paper investigates the potential of the discrete element method to simulate the physics of particle bed-based additive manufacturing. This method naturally captures the discrete aspects of additive manufacturing processes, such as material addition. The proposed discrete element framework uses constitutive relations for loose powder, bonding kinematics and the thermo-mechanical behaviour of bonded particles. The mechanical bond interactions consist of beams that interconnect the particles. These beams are able to transfer forces as well as moments. The thermal conductive bond interactions assume an effective conductive area and density to account for the voids in the system. Simulated compression tests reveal that the macroscopic Young’s modulus and Poisson’s ratio of the bonded material are controlled by only two micro-scale parameters. Furthermore, a heat conducting rod of both powder and bonded material is simulated and compared to a continuum finite element simulation. The proposed discrete model is able to simulate a complete printing process, capturing the solid material behaviour accurately. A simulation of a printed sample shows various additive manufacturing aspects such as: the deposited powder layer, G-code input, heat source interaction, contact, bonding, thermal conduction and the accumulation of residual stresses and deformations.


Introduction
Over the past years, significant progress has been made in the development of various Additive Manufacturing (AM) techniques. As a result, a variety of materials can be printed with an ever increasing quality and precision. Nevertheless, many aspects of the process rely on the experience obtained in previous print jobs. Especially in the development of new products, a number of trial runs are often required before the desired product quality is achieved [1,2]. In principle, numerical simulations are able to provide a fundamental understanding of the AM processes and can be of use in predicting and correcting phenomena that may result in undesired product properties or shape [3].
The numerical analysis of AM processes poses various challenges from a simulation perspective. Generally, conventional modelling techniques, such as the Finite Element Method (FEM), require a continuous remeshing as material is added to the domain during analysis [1]. In addition, complex constitutive relations in combination with a coupled thermo-mechanical analysis are required to correctly represent the phase transition of materials [1,2]. Finally, the analysis spans a wide range of time-and length scales, which calls for innovative modelling techniques to reduce computational costs [1][2][3][4].
Regarding the remeshing in FEM, the addition of new material is typically done using quiet elements or inactive element activation. In the method of quiet elements, the complete product is discretised in advance, whereby the material properties are set to negligible values that are increased when the material is actually added [5].
Marc G. D. Geers and Joris J. C. Remmers contributed equally to this work.
Alternatively, in the inactive elements method, the elements are only algorithmically introduced when the corresponding material is added [5]. The biggest drawback of both methods is that since the domain is discretised beforehand, the placement of the elements is rigid and forces new material to be aligned with already placed material [6]. As a result, these methods are not able to accurately capture misalignment issues during printing, which is the case when the new material is not placed correctly on top of the previous layer, that has deformed due to thermal or chemical shrinkage.
To capture these effects, the Discrete Element Method (DEM) is a suitable simulation alternative. In most AM processes, the addition and manipulation of material is done locally and in a step-by-step fashion, which can be closely represented by a DEM method. The method excels in handling discrete variations and is not constrained to a specific grid or shape [7,8]. Moreover, the failure mechanism often observed in AM products is fracture due to residual stresses, which essentially is a discrete phenomenon that can be implemented in DEM in a natural fashion [9][10][11][12][13].
To accurately predict the properties and shape of an AM product it is important to include the local small-scale variations in the material. When using DEM, distinct interfaces between different phases, materials or layers are captured naturally when different powder size distributions or material properties are used. The evolution of microstructures, due to geometrical particle stacking for example, is accounted for naturally in DEM. This is in contrast with the aforementioned FEM simulations that homogenise the bulk material properties, which results in loss of resolution and loss of heterogeneity [1,3,6].
The discrete element method has proven to be suitable for simulating bulk material processes, such as granular flow or powder mechanics [14][15][16][17]. Deposition of the powder is a key feature in many bed based AM applications since the packing density of the deposited layer influences the resulting strength of a product. The thermodynamic behaviour of AM processes has also been simulated with DEM, to study heat transfer in a powder bed or to analyse a melt pool [18][19][20][21][22]. Combining the mechanical and thermodynamical processes results in a print process simulation, which has only been accomplished in a limited number of papers [23][24][25][26] that are discussed below.
Zohdi [23] presented a framework that focuses on particle movement, contact forces, near field interaction, laser absorption, beam interference, conductive heat transfer and thermal softening. The main focus of that paper is a droplet of particles being deposited and subsequently heated with a laser. Particles are bonded when a critical interpenetration distance is reached. Various numerical input parameters, such as the normal contact or bond compliance, are not directly linked to physical material properties.
Steuben et al. [24] presented a basic DEM framework for AM with contact, bonds, conduction, laser heating and shows the simulation of subsequent layers in a Selective Laser Sintering (SLS) process. To capture some essential trends, relatively simple constitutive models and boundary conditions have been applied. For example, the mechanical behaviour of bonded particles is determined by a normal bond stiffness, which is chosen equal to the contact stiffness. The normal stiffness is computed with the Hertz-Mindlin model and a non-physical parameter for the drag coefficient is used.
Lee et al. [25] included contact, conduction and laser interactions, which allows to simulate powder deposition, re-coating, laser heating and the holding stage. The paper primarily focuses on thermodynamic aspects, but is limited due to the uniform laser source, uniform particle size and the relatively low number of particles. The effect of the laser settings on the temperature distribution of a powder bed are studied, but the subsequent bonding of particles was not included.
Xin et al. [26] introduced the most elaborate framework to-date. This model includes bonding, necking, thermal expansion and softening, phase-change, heat conduction and laser heating. It is able to simulate the formation of subsequent layers during a SLS process, enabling computation of the material properties of the product in a second framework. However, the simulation includes only a simple thermodynamic model for both bonded and non-bonded particles, whereby the mechanical bond parameters depend on the size distribution, packing and constitutive laws.
In the present paper, we present a multi-physics, thermo-mechanically coupled discrete element framework for AM processes that is able to simulate a powderbed-based AM process from start to finish. This framework departs from the work of Xin et al. and takes into account the proper physical behaviour of the bulk powder, the solidification process and the constitutive behaviour of the sintered solid material. The input parameters are uniquely linked to the target material properties, i.e. a specific material can be unambiguously modelled in the powder state as well as the solid state by defining its corresponding input parameters. Consequently, the framework can be used to gain fundamental understanding of the process and to predict the resulting properties. This framework enables the simulation of the print process of a powder-bed based AM process, allowing for a detailed analysis from powder to product.
The key attribute of this paper is a thermo-mechanically coupled DEM framework that is capable of simulating the expected constitutive behaviour of powder, solidification and solid material. Simulating the constitutive behaviour of a solid material is not straightforward with DEM, as there are "intrinsic numerical" voids due to the discrete elements, which affect the mechanical properties. Research is done into simulating solids with DEM, for both mechanical [27][28][29][30][31][32] and thermodynamical problems [27,[33][34][35]. The framework assumes different models for non-bonded and bonded particles to comply with the material structure differences.
The paper is ordered as follows. In the next section, the proposed computational framework is elaborated. The constitutive equations for the mechanics and thermodynamics and the associated explicit integration methods are presented. Next, in Sect. 3, a convergence study is performed, to identify the fine-scale parameters to simulate a continuum-like solid. To validate the model, a compression and a heat conduction test are performed in Sect. 4, highlighting the accuracy of the model for both the powder and the solid behaviour. In Sect. 5, the capabilities of the proposed framework are demonstrated with a representative additive manufacturing use case. Finally, conclusions and recommendations are drawn in Sect. 6.
Throughout the paper, the following mathematical notations are used. Scalars are described by normal font characters. Vectors and quaternions are depicted by bold characters. Subscripts i and j are used to indicate which particle(s) the variable is associated with.

DEM framework
In the model, both the powder and the solid material are represented by spherical particles. The state of these particles is governed by coupled thermo-mechanical balance equations. The interactions between the particles in the powder phase and the particles representing a solid phase are described in Sects. 2.1 and 2.2, respectively.
In all cases, we consider a 3D domain Ω that is filled with spherical particles i with different radii R i , as shown in Fig. 1. Each particle has a position x i , a current rotation represented by a quaternion q i and a homogeneous temperature T i .

Powder phase
For particles in the powder phase the DEM framework includes contact mechanics and thermodynamic conduction. Contact is represented by two particles having a small overlap, see Fig. 2a.
A particle in contact with another particle generates a repelling contact force, see Fig. 2b. In the model, two particles are in contact when the overlapping distance u i−j is positive. The overlap between two particles i and j is defined as The unit vector normal to the contact surface, pointing from particle i to j, is defined as The contact force F c i is assumed to be linearly dependent on u i−j : where k is the contact stiffness and t i−j is the tangential displacement vector, depending on the relative velocity [26,36,37]. is derived via dimensional analysis and depends on (1) (2)  particle size distribution, packing and constitutive laws [26]. The moments acting on a particle in contact are generated by the tangential part of the force and its arm r i : The effective behaviour of a system of particles is controlled by the contact stiffness k in Eq. (3). By properly computing this parameter the properties of a powder can be accounted for. In this study, we use the Cundall contact law, which relates the stiffness to the elasticity modulus E and radii R i and R j by assuming a cylindrical rod of load-carrying material with length D ij and cross section A ij in between the particles [26]: It is assumed that the temperature T i is constant in a particle [38]. Convection and radiation are ignored, which means that heat can only transfer between two particles when they are in contact [24,25]. The conductive heat transfer c i scales with the aforementioned contact area A ij , which increases when particles are pressed together. Similar to the simulation of the contact force, a cylindrical rod of material, that effectively governs the thermal conduction, is adopted between the particles. The conductive heat transfer scales with the difference in temperature between the particles, which is positive for particle i if the temperature of the adjacent particle j is higher [24,25,38,39]: In this equation, is the thermal conductivity.

Solid phase
The mechanical behaviour of a solid is approximated by means of bonded particles and a beam bond model [10]. Each bond is represented by a massless cylindrical beam with a specific modulus and beam radius, see Fig. 3. This beam is able to transfer elongation and bending forces as well as bending and twisting moments. Different from a contact force, a bond force can be either tensile or compressive. This concept is similar to that of a lattice/network model with mass-points at the intersecting points and beams as the interconnects, albeit without a fixed grid [39,40]. It is assumed that the relative displacement of two particles is small, thus the deformation of the beams fits in the infinitesimal theory. Therefore, the forces and moments can be computed with the analytical solutions of the Euler-Bernoulli beam equations [27][28][29]. Figure 3 shows a schematic representation of a bond in undeformed and deformed configuration. When the bonding criteria are met, a bond is formed between two particles. The state at creation of the bond is considered to be the equilibrium state, the particles do not overlap for clarity of the graph. Figure 3a depicts two bonded particles in equilibrium with their respective rotations in the global domain Ω . Whenever a bond is deformed as depicted in Fig. 3b, a 2D plane can be found that confines the bending. The initial and current orientations indicate how the particles should rotate with respect to the bond vector to be in equilibrium. The beam tends to return to its equilibrium state and the mismatch with the deformed state results in forces and moments that drive the translations and rotations of the particles.
The local beam directions, a set of unit vectors, have to be chosen in a systematic way to exploit the analytical beam solutions. As bending of a cylindrical beam takes place in a 2D plane, the Y i vector that describes the bending, should be parallel to this plane. The X vector represents the bond vector, which is the vector e ij pointing from particle i to j, capturing the torsion in this plane. The X i vector corresponds to the current deformed state of the bond, which can be computed from the deformation rotation i . The Y i vector is chosen perpendicular to the X i and X vector following a right hand rule. Lastly, the Z i vector is set perpendicular to the X and Y i vector to create a right handed coordinate system. The analytical solutions, in the reference frame of the beam, for discrete element i are [28]: In Eq. (7), E is the microscopic elasticity modulus of the beam, A b is the beam cross-sectional area, D b the equilibrium bond length and I b the area moment of inertia of the beam. y i is the angle between the equilibrium and deformed states of particle i around the Y-axis. Equation (8) includes J b , which is the polar area moment of inertia of the beam, and G , which is the microscopic shear modulus. The shear modulus is computed with the microscopic elasticity modulus and the microscopic Poisson's ratio, which is taken equal to the macroscopic Poisson's ratio [28]. The areas and area moments of inertia are computed with the beam radius R b . The beam radius is linked to the radii of the corresponding particles by the bond beam ratio r : The complete deformation of the bond is described by elongation, bending and torsion. To compute the forces and moments of a deformed bond the following steps apply. The bending of the cylindrical beam, i.e. bond rotation, takes place in the 2D bending plane, see Fig. 4. The bond rotation quaternion q b is obtained with Eq. (34) by computing the angle b and rotation axis n b from the initial bond vector b 0 to the current bond vector X: The equilibrium orientation quaternion q eq i in the current configuration is obtained by adding the initial bonding orientation q 0 i and the bond rotation quaternions, according to Eq. (36): The deformation rotation quaternion i from the deformed state to the equilibrium state can be computed by subtracting the current orientation quaternion q i from the equilibrium orientation: The deformation rotation can be decomposed into angles around the X -and Y i -axis, torsion x i and bending y i , respectively, by the so-called swing-twist decomposition, see (37). Now, the forces and moments can be computed with Eqs. (7) and (8). The resulting forces and moments have to be translated back to the global coordinates of domain Ω.
As opposed to the powder state, the assumption that conduction only takes place through the overlapping crosssectional area between two particles falls short. The discrete elements should represent a volume without the remaining artificial voids introduced by the discrete element method [34]. An effective conductive heat transfer area and volume fraction should therefore be used in the computations. For a regularly stacked domain, the effective area and volume  A 2D view of a regular cubic packing, which can be exactly converted to a continuum-like solid for heat conduction fraction can be calculated exactly, such as for the cubic stacking shown in Fig. 5. However, for an irregular domain this is not as straightforward and an alternative approach needs to be used.
Studying the packing of cubic volumes with edges of size 1 that are filled with a single discrete element would result in each particle having 6 neighbours in three dimensions (i.e. the coordination number). These cubic volumes, with edges of twice the size of the radius of a particle, reflect a simple cubic packing. The effective volume fraction f e of the discrete element is computed by dividing the volume of the particle by its projected cubic volume: 4 The surface factor is calculated by: S e = 1 2 ∕ 1 2 2 = 4 . The surface factor is used to compute the continuum-like solid conductive area by multiplying the surface factor with the square of the radius of the particle.
A handful of shapes can completely fill a 3D domain: tetrahedron, cube (hexahedron), octahedron, dodecahedron and icosahedron. Each shape corresponds to a specific stacking which results in a distinct number of neighbours. The effective conductive area and volume fraction [27,33,35], which corresponds to a specific number of neighbors, can therefore be extrapolated from the regular packings, as shown in Fig. 6. This figure is used to determine what values to use for S e and f e for each discrete element, by looking at their number of neighbors.
Implementing the method above requires a modification in the computation of the conductive heat transfer and an adaptation to the mass of particle i. When calculating the conductive heat transfer in Eq. (6), the area A ij is replaced with the effective conductive area A e ij representing a solid: The mass m i is substituted by the effective mass m e i : To determine the effective surface area and volume fraction, the number of neighbours for each particle should be known at each iteration, which is determined during the interaction computations.

Thermo-mechanical coupling
When particles are sintered, mass moves from the particle volumes to the link in between [19,26], as shown in Fig. 7. This is referred to as neck growth. It effectively brings the mass centres of the particles closer together, reducing the bond length D b , while their radii are assumed to remain unchanged. A bond is then created through which forces and moments can be transferred. A pair of particles bond when the particles physically overlap and their temperatures exceed the bonding temperature T s : This bonding criterion can be extended, for example, by requiring a minimum contact pressure or a maximum relative velocity [24].
The particles and bonds experience both thermal expansion and shrinkage. This is implemented by calculating the change of a certain length with a linear coefficient of thermal expansion L times the temperature difference ΔT with respect to the initial temperature. Both radius R i of a particle and the equilibrium length of a bond D b are affected by thermal expansion: where R 0 i and D 0 b are the initial radius and initial equilibrium bond length.

Damping and gravity
To eliminate dynamic oscillations, damping is added. It takes the form of a force or moment for each interaction of each particle and scales with (angular) velocity, v i and i . The damping is controlled by a damping parameter 0 < ≤ 1 , which defines the percentage of elastic elastic energy that is conserved [41]. It can be used to model a physical viscous-dissipative material or to attain an equilibrium state faster: In Eq. (18) m eq is the reduced mass. I is a mass moment of inertia and k is a rotational stiffness in Eq. (19). The reduced mass is given by In the case of contact, the rotational contact stiffness and the mass moment of inertia of the corresponding particle are computed as In the case of a bond, the rotational bond stiffness and the mass moment of inertia, calculated with the reduced mass of the bond, are computed with Finally, each discrete element experiences a gravity force: where g 0 is the gravitational constant and e g the unit vector pointing in the direction of the gravity.

Explicit time integration
To move forward in time, Newton's second law of motion is solved for translations in Eq. (24) and rotations in Eq. (25) for each particle: where a i is the acceleration vector and i the angular acceleration vector of particle i. I i is the mass moment of inertia of particle i, which can be calculated with Eq. (21). The degrees of freedom (DOF) for each individual particle are explicitly updated every time increment. The translation and rotational DOF are updated with the Velocity-Verlet integration scheme [24,42], except for the current rotation, which is updated using quaternion operations. The temperature of a particle is updated using a discrete first-order integration scheme solving the thermal energy balance [18,25,26,33,38]: where c i is the specific heat capacity of particle i, m i its mass and h i the heat added to the particle by the heat source. For each translational DOF there is an equivalent rotational DOF, however the integration of rotations follows a slightly different scheme than the positions [43]. The change in rotation of the current iteration Δq i is computed from the half-step angular velocity

Implementation details
To simulate an AM printing process in a DEM framework, various algorithmic details need to be clarified.

Control and boundary conditions
The trajectory of the heat source as function of time and other print process instructions are given as G-code input. A heat source is modelled by assuming a Gaussian distribution of intensity at the laser spot position on the powder layer. Depending on the power, spot size and powder properties, the heat source intensity decreases over penetration depth, according to the Beer-Lambert law [4,18,38,44,45]. Heat is added to a particle according to (24) Δt and n i = where I h max is the peak intensity of the laser and I h 0 prescribes the normal laser intensity profile at the top surface as a function of x i and y i , which are the coordinates of particle i with respect to the heat source; is the attenuation coefficient and a the absorptance. The A prj i is the projected cross-sectional area of particle i and z i is the depth coordinate of the particle with respect to the heat source.
In this paper, the laser intensity on the surface is assumed to take a Gaussian bell shape. It is determined by the laser settings: laser power P L and laser radius R L . The peak intensity and surface intensity profile are described by The heat induced by the laser as a function of depth is shown in Fig. 8.
Although the focus of the model is on predicting the material transformation it is also possible to model particle deposition. However, in this paper, we will realise representative layers of powder using a deposition procedure, where the particles are dropped from a small height. Particles can be given random perturbations to their initial positions and a random velocity. The particles initially start at a fraction of the minimum radius and grow over a short period of time to their actual radius. The dropping and growing of particles increases randomisation and results in a densely packed powder bed.

Stresses
Stresses are computed using the Cauchy-Born hypothesis, which relates the Cauchy stress tensor to forces and displacements. This method is originally deduced from structured lattices, but performs adequately in most DEM simulations [36,[46][47][48][49]. The stress tensor of a selected volume V is computed, summing over all interaction pairs N p in its domain, according to the following equation: where F i are the interaction forces acting on particle i. Velocities and moments do not contribute to the stresses in this expression. As this paper does not focus on the dynamic behaviour of the powder and the velocities are small, the tangential force and moments for powder can be neglected.

Numerical optimization
To optimize the particle interaction computations, a cell-domain algorithm is implemented [38,42,50]. The domain is split in cells of equal size, which is larger than the largest particle diameter. At the start of each increment the particles are allocated to their corresponding cell, particle interactions are only checked and computed for particles in their specific cell and their neighbouring cells. This drastically reduces the dependency of the computation time on number of particles n from O(n 2 ) to O(n log n) [42,50].
Finally, to remove the effects of boundaries and ensure faster convergence, periodic boundaries are implemented. Periodic boundaries imply that the simulation repeats periodically across its boundaries. The boundaries in X-and Y-direction are made periodic. When the bonds between particles are created, the particles close to the boundaries at one side are virtually copied to the other side and bonding is computed accordingly. The deformation of a bond is tracked over time and bonding parameters are updated. The interactions of bonded particles are computed with these updated parameters, resulting in a periodically bonded simulation, where the boundaries are free to move.

Parameter identification for the solid phase
To obtain the expected continuum material behaviour at the macro scale the corresponding parameters of the beam bond model should be identified. To obtain these microscale parameters a compression test is performed. The macroscopic Poisson's ratio and Young's modulus E of the sample can be determined with the Cauchy stress and engineering strain , assuming that the sample is isotropic at the macroscopic level. The sample is expected to contain a comparable amount of interactions in all directions, resulting in global isotropy. The stresses are determined by summing the forces of the particles at the planes between the prescribed and free particles, or at the outer boundaries, and dividing by the corresponding area. The strains are obtained by calculating the current dimensions with the average displacement of these planes and comparing them to the initial dimensions. To compute the Poisson's ratio and Young's modulus the following relations are used: Before identifying the parameters, the convergence of the results is assessed for the Poisson's ratio and Young's modulus. For this purpose, the sample in Fig. 9 Figure 10 reveals that approximately 2000 particles are required for an adequately converged solution of the Poisson's ratio and Young's modulus. For a converged solution of 18,374 particles, the beam ratio can be varied from 0.1 to 1.0, which affects the macroscopic mechanical properties. The Poisson's ratio and Young's modulus are plotted as a function of the microscopic beam ratio in Fig. 11a. The Poisson's ratio is clearly governed by the beam ratio. The same analysis can be repeated by varying the microscopic modulus, ranging from 2 kPa to 200 GPa. The Poisson's ratio and Young's modulus are plotted as a function of the microscopic modulus for a beam ratio of 0.3 in Fig. 11b.   Fig. 9 Overview of the compression sample. The prescribed layer is shown in red, the constrained layer in grey and the free particles in blue. The shaded green area shows the dimensions of the domain that is approximated  The graph shows that the microscopic modulus minimally affects the macroscopic Poisson's ratio. This indicates that the macroscopic mechanical properties are each controlled by a single input parameter: first the Poisson's ratio through the beam ratio r and next the Young's modulus through the microscopic elasticity modulus E .

Validation powder and solid behaviour
The implemented models in Sect. 2

Compression test
The input parameters for a solid material were determined from Figs. 11a, b: a beam radius ratio of 0.3 for ≈ 0.36 and a microscopic elasticity modulus of 5.0 ⋅ 10 6 Pa for a macroscopic E ≈ 0.23 MPa. With these parameters and Table 1, a compression test is performed on a solid sample. The same compression simulation procedure is used as for the parameter identification, as described in Sect. 3. This compression test for the continuum-like solid results in a Poisson's ratio of approximately 0.357 and a Young's modulus of 2.37 ⋅ 10 5 Pa.
In Fig. 12, the evolution of the Poisson's ratio and Young's modulus during the compression test is shown, computed with Eqs. (32) and (31). At small strains, stresses are (close to) zero due to some initial rearranging of the particles, resulting in asymptotic behaviour for the Poisson's ratio and Young's modulus, i.e. the actual values are more reliable for larger strain values. The Young's modulus is also determined by fitting the slope of the compressive stress, computed with Eq. (30), versus the strain, shown by the green line.

Heat conduction test
To validate the thermodynamic solid approximation of the DEM framework the following test-case is used: a conductive rod with two different boundary temperatures at opposing sides and no heat conduction through the other boundaries. Figure 13 shows a rod of dimensions 0.04 × 0.01 × 0.01 m, which is bonded and initially at 100 • C. A temperature of 500 • C is imposed on the left boundary, while the right boundary is kept at 0 • C. The temperature profile along the length of the rod is compared to a converged finite element simulation [51] at different time intervals.  Multiple simulations have been carried out with different amounts of particles to discretise the rod: 225, 744, 1531, 3695, 6534, 13,007, 30,800 and 54,143 particles. Even for a small number of particles, the shape of the temperature profile is clearly visible. The results for the temperature profile converge adequately for just a few thousand particles. For fewer particles, the overshoot, when the temperature increases, and the undershoot, when the temperature decreases, become more significant. Figure 14 shows a section of the temperature profile at approximately a quarter of the length of the rod for different amounts of particles. Figure 14 reveals that the spread of the temperature in the rod is more significant for simulations with fewer particles, where each point represents a discrete element in the DEM simulation. The overshoot can also be observed: on average, for 3695 particles the difference with the FEM solution is larger than for 13,007 particles. With 13,007 particles, shown with all particles in Fig. 15, the maximum error during the simulation is around 5 • C, which is approximately 1 % of the total temperature range. Figure 15 reveals that the DEM temperature profiles follow the finite element simulations closely over time. At every time interval, the points are scattered in a small region around the FEM temperature profiles (dashed lines). Hence, the thermodynamic continuum approximation is adequate and can be used to simulate transient heat conduction problems in a solid.
To study the conduction of loose powder, the same method is used as for the finished product material, but without bonding the sample. The results are shown in Fig. 16. The heat conductivity of a powder bed of steel particles is experimentally measured to be around 0.6 W/mK, about 100 times lower than a solid sample [52][53][54]. This experimentally measured effective conductivity of a powder bed is used to simulate the FEM result.
Again the results are scattered around the FEM simulation in Fig. 16. However, this time the scattering is more significant, which is the result of the heterogeneities caused by the powder geometry. The FEM result assumes equal conductivity in all directions, while in DEM heat conduction only occurs in contact direction. Some particles have smaller  contact areas and therefore conduct less heat, which result in a large disparity in temperature along the width of the rod. As expected, the powder bed performs worse in terms of conducting the heat from one boundary to the other, which is why the temperature profile shows two steep slopes at either side and a central part where the temperature does not change much over the time.

AM use case
To illustrate the capabilities of the DEM framework, a typical AM use-case is presented. An H-shape is printed with outer dimensions of 3.3 × 3.3 × 2.2 mm with layers of 0.2 mm [55], see Fig. 17. A base with a height of 0.2 mm is printed and between a height of 1.2 to 1.6 mm the overhang is printed. The bodies at the left and right side of the H-shape have a width of 0.9 mm. A spot size of 0.3 mm is used with a power of 300 W and a printing speed of 20000 mm/min. The material is descretised by spherical discrete elements with a uniformly distributed radius between 35 to 50 micrometer. The material parameters from Table 1 are used, except for the Young's modulus ( E = 2.0 ⋅ 10 4 Pa), microscopic modulus ( E = 2.0 ⋅ 10 5 Pa) and thermal conductivity (20.0 W/mK). The shrinkage after initial bonding due to neckgrowth is set to 5%.
To start the simulation of the additive manufacturing process, the powder deposition of a layer should be captured in a representative way. This is done at the onset of each printing step with the deposition procedure of Sect. 2.6. The result of the deposition is presented in Fig. 18a. During the simulation, it can also be seen how the powder bed is sintered by the heat source moving according to the G-code input in Fig. 18b. An approximation can be made whether the input achieves the desired result. The spacing, direction and lay-out of scanning lines might have an influence on the printing process, which can be analysed by means of the figure. Figure 18b shows the laser interaction with the powder bed. The Gaussian bell shape can be recognised in Fig. 18c. The laser scanning line has a width of approximately 0.3 mm. The laser penetration depth is more difficult to asses during the printing process due to the previously printed layers. However, when the overhang is printed, excessive material addition is observed at the overhang of the H-shape. The laser penetrates the material with a depth of approximately 0.25 mm, while layers of only 0.2 mm were printed. With the DEM framework, various scanning settings can be studied to minimize this excessive material addition, for example by changing the spot size and power in certain parts of the print.
During the printing process, the material expands and shrinks due to thermal expansion. In addition, the bonded material shrinks due to the occurrence of neck growth. The expansion and shrinkage of the material result in internal stresses in the finished product and powder bed. These stresses emerge during the printing process, see Fig. 18d. When the product is removed from the base plate it typically warps due to the presence of residual stresses. The material takes a configuration in which the residual stresses are relaxed, which is usually not the desired product shape.
The H-shape shows distinct distortions after the simulation, see Fig. 19. The overall resulting deformations remain comparable to an experimentally printed sample [55]. The joints of the H-shape are deformed inwards, due to the shrinkage of the overhang during printing. The sidewalls are curved upwards, visible from both the front and side faces, due to thermal strains during printing, where a new layer is deposited on a shrinking prior layer. Some layers are observed that have detached, i.e. bonding between subsequent layers was not sufficient in these cases. Another part that does not have the intended shape is the overhang of the H-shape. Due to the laser penetration depth extra material is sintered at the bottom of the layer, known as the z-bonus.
The simulation shows that the settings used for printing do not result in the desired final product. To gain insight into the cause of the undesired printing phenomena different process aspects can be studied, such as the powder distribution, the scanning pattern, heat source interaction and deformation of the bonded particles. For this specific case, the majority of the deformations seem to originate from the bonding of particles and the subsequent shrinkage. With this framework, different scanning strategies can be assessed to verify if the effect of the shrinkage due to neckgrowth can be minimised.

Conclusion
In this paper, a discrete element framework is presented that is able to simulate a complete powder bed-based AM process. The framework describes and captures the thermomechanical properties of both the powder and solid phase in a physically adequate manner. Additionally, it includes a model to simulate the bonding. The input parameters for solid material can be identified by performing compression tests. The simulation of a printing process allows to assess the resulting properties and the full history of a product. Complex process parameters, such as residual stress formation or laser-powder interaction, can be scrutinised to gain fundamental understanding of the process.
The solid particle interactions are validated by simulating a compression test and a heat conductive rod. Forces, displacements and stresses of the compressed sample are analysed and a representative stress-strain curve is retrieved from the simulation. The mechanical properties of the DEM framework are extracted and shown to be in agreement with the expected material properties. Heat conduction in a rod is simulated with the DEM framework and compared with the converged FEM equivalent. The temperature profile of both the powder and solid material simulation are in close agreement with the FEM result at different time intervals and positions. The mechanical and thermodynamic solid models show promising results for simulation of a solid material with DEM.
A sample is printed to demonstrate the complete simulation of an AM process. The simulation is able to capture multiple aspects of the AM printing process such as: the deposited powder layer, G-code input, heat source interaction, contact, bonding, heat conduction, (residual) stresses and deformations. Deformation, shrinkage and warpage are observed, which are important for controlling the dimensional accuracy of products. Powder deposition is simulated by a procedure which reproduces a representative powder layer. The simulated laser interaction tracks the described G-code, which controls the melt pool dimensions. The particle interactions result in the final printed product, in which Fig. 18 Pictures of the DEM framework simulation of an H-shape during different parts of the printing process stresses and deformations can be adequately assessed. The proposed framework enables the simulation of the print process of powder-bed based AM processes, allowing for a detailed analysis from powder to product.

Rotations and quaternions
To keep track of the rotations of each individual discrete element quaternions are used. A quaternion is able describe a 3D rotation with just four parameters, similar to complex numbers with i , j and k imaginary, being: Quaternions contain a scalar or real part w and a vector or imaginary part x, y and z. The scalar part contains information about the rotation angle , while the vector part contains information about the rotation axis n . A quaternion can therefore be computed from an angle-axis representation as follows: A quaternion can also be computed between two vectors or from any other rotation representation. When representing rotations, a unit quaternion is used, which has a norm of 1. For unit quaternions the inverse q −1 is equal to the conjugate q * , which is simply computed by multiplying the imaginary (33) q = w + xi + yj + zk. part with −1 . Taking the inverse of a quaternion also inverses the rotation it represents, thus flipping its direction. A quaternion can represent a vector when w equals 0 and the imaginary part equals the vector. When using quaternions, a vector in quaternion representation p can be easily rotated by a quaternion q with the Hamilton product * : The addition and subtraction of multiple rotations is relatively simple and compact, but not commutative. The total rotation q 12 of rotation q 1 followed by rotation q 2 is computed by Any arbitrary rotation can be decomposed into two independent rotations around a specified axis: the first one around and the second one perpendicular to the specified axis. The decomposition of a rotation quaternion q r around an axis p a by computing the swing rotation quaternion q s after the twist rotation quaternion q t is called the swing-twist decomposition: First, the vector projection p t of the rotation axis p r of q r onto the specified axis of decomposition p a is needed for the construction of the twist: With the twist, the swing can be computed according to Eq. (37), where the twist q t is defined by: (36) q 12 = q 2 * q 1 .
(37) q r = q s * q t thus q s = q r * q −1 t .
(38) p t = p r ⋅ p a p a ⋅ p a p a Fig. 19 Comparison of the intended print in green and the resulting print with discrete elements