Linking particle properties to layer characteristics: Discrete element modelling of cohesive fine powder spreading in additive manufacturing

Abstract Particle properties play a key role in determining the flowability of cohesive fine powders and consequently the quality of parts fabricated by additive manufacturing. However, how the characteristics of a spread powder layer are linked to particle properties remains not well understood, largely due to the limitations of available numerical models and characterisation approaches. This study thus established an efficient discrete element modelling framework to address these issues, combining GPU computing with a novel methodology for particle stiffness scaling. The validity of the DEM model was verified on both packing and spreading of cohesive fine powders. The effect of particle cohesion on the characteristics of a spread powder layer was systematically analysed for non-cohesive and weakly to strongly cohesive powders. The structure of spread powder layer was qualitatively illustrated and quantitatively evaluated by not only commonly used metrics, such as coordination number, packing density and surface profiles but also particle clustering and pore characteristics. While reducing particle cohesion leads to an enhanced powder flowability, different regimes were identified in the relationship between layer characteristics and particle cohesion. The results showed that powder spreadability has a complicated dependence on the strength of particle cohesion, where the underlying mechanisms can be understood in part via a dimensionless inertial length. These findings not only provide valuable metrics to quantitatively evaluate the quality of a spread powder layer, but also enable a better understanding of the physics underlying the spreading of cohesive fine powders.


Introduction
Additive manufacturing (AM) has attracted significant attentions across a broad range of industrial sectors due to its advantages of fast prototyping and near-net-shape production. In powder-bed-based AM processes, particles within selected areas are fused together by either a high energy source (laser, electron beam or plasma) or a binder, in a layer-upon-layer fashion to build up objects of complex shapes [1]. The quality of additively manufactured part is affected by many parameters, including both process parameters such as laser power, scanning speed, spot size, preheat temperature and feedstock properties such as powder surface chemistry and morphology [2][3][4][5][6]. In these processes, powder material is often deposited with the aid of a spreading blade or a rotating roller, in which particles are subject to localized shearing and compression. The quality of spread powder layer plays a key role in determining residual porosity and sintering strength and hence the structural and mechanical properties of the fabricated parts [7][8][9][10]. A poor spreading may lead to undesirable features, such as patchy coverage, size segregation and non-uniform density distribution. An in-depth understanding of the link between microscopic particle properties and characteristics of a spread powder layer is thus critical for the evaluation of powder spreadability and optimization of part quality.
Considerable efforts have been devoted to investigating powder spreading process using both experimental and numerical approaches. Standardized tests have been used to assess powder flowability for powder bed fusion, such as funnel flowmeter, static angle of repose, rotating drum, shear cell, FT4 powder rheometer and ball indentation [11][12][13][14][15]. However, stress conditions within these tests are different from that of powder spreading. A testing device of relevance to the process conditions is thus required to either probe powder behaviour during spreading or to measure powder spreadability [16][17][18][19][20]. Van den Eynde et al. [17] proposed a lab-scale powder spreader to mimic powder spreading in laser sintering, which allows the overall packing density of deposited powder layer to be measured on a balance. Sun et al. [18] developed a powder deposition system identical to that used in an EBM machine, where surface quality of spread powder layer is analysed by processing images taken by a digital camera. Ahmed et al. [19] used a cutter blade to manually spread a small heap of powder. The presence of empty patches is extracted from SEM images of a spread powder layer, where particles are immobilised by an adhesive spray. Particle-scale characterisation of a deposited powder layer and direct observation of spreading, however, is challenging due to the lack of transparency of the spreading rig and small particle size, making it difficult to track the motion of individual particle [21,22]. In situ observation using high-speed high-energy x-ray imaging has been reported but it is at present limited to dynamics of free surface flow in front of spreading blade [22,23]. Most importantly, there is a lack of quantitative understanding, as the effects of particle properties are difficult to isolate in powder feedstock. In view of these inherent limitations of experimental approaches, discrete element modelling (DEM) has shown great potential to understand powder spreading from the particle scale, as it allows interparticle forces of various sources to be explicitly considered and the resulting packing structure to be readily characterised. Application of DEM has revealed that layer homogeneity is affected by both powder properties (such as, particle cohesiveness [21,24,25], sliding friction [21], particle shape [26,27], size [24] and its distribution [28]) and process parameters (such as, spreading speed [26,29,30], blade clearance [21,31] and blade type [25,26,29,30,32,33]). DEM is also helpful to understand physics underlying powder spreading [31,34]. Nan et al. [31] attributed the presence of empty patches to the occurrence of transient jamming in front of the blade. Recently, Fouda and Bayly [34] applied DEM to benchmark the behaviour of non-cohesive spherical particles. The reduction of packing density in a deposited powder layer compared to the heap in front of the blade was attributed to three major mechanisms: shear-induced dilation, particle rearrangement in the blade gap and particle inertia. To date, few DEM studies have been dedicated to the role of particle cohesion, despite particle cohesion presents a strong impact on the spreading behaviour and subsequently uniformity of the deposited layer [21,24,25,28]. Chen et al. [21] showed that increasing particle cohesion leads to an increased dynamic angle of repose in front of spreading blade and a reduced continuity of powder flow during spreading. Meier et al. [24] found that increasing particle cohesion leads to a decreased packing density but an increased surface roughness. The same conclusion was also drawn by Wang et al. [25]. However, different findings were reported regarding the homogeneity of powder bed. Meier et al. [24] reported a decreased layer homogeneity with increasing cohesive forces while Wang et al. [24] found the layer homogeneity first increases and then decreases with particle cohesion. Despite being novel in understanding the spreading process, the layer characteristics and the underlying mechanisms of cohesive fine powder spreading remain not fully understood.
It should also be noted that most DEM studies mentioned above adopted a reduced particle stiffness that is several orders of magnitude smaller than the true physical value. This is because of the unrealistically small timestep required if using the original particle stiffness, making DEM simulations computationally unaffordable. A reduced particle stiffness is thus normally adopted to accelerate DEM calculations. For non-cohesive particles, reduced particle stiffness has little impact on bulk behaviour, if overlap between particles remains small. However, this is not the case for cohesive particles, particle cohesion needs to be scaled with particle stiffness, otherwise more kinetic energy would be dissipated due to an enlarged interparticle overlap. Several scaling laws have been proposed for this purpose, depending on the force models used in DEM modelling [35,36]. For example, in the work of Nan et al. [31], the experimentally measured surface energy (i.e. 9 mJ/m 2 ) was scaled with particle stiffness from 211 GPa to 2.1 GPa, giving a value of 1.4 mJ/m 2 for use in the spreading simulation. However, our recent study showed that these existing scaling laws cannot preserve the original particle behaviour [37]. In particular, the sliding and rolling resistances are largely underestimated, undermining the predictive capability of DEM and consequently confidence on the results. This is especially the case for powder spreading where enduring contacts dominate the spreading process. Not only is the validity of the existing scaling approaches questionable, but also the previous studies are quite limited in characterising the spread powder layer, primarily focusing on packing density and surface roughness. To the best of our knowledge, little work has been reported to systematically address pore characterises of a spread powder layer, although it directly linked to some key phenomena in additive manufacturing, such as laser absorption, particle sintering and the presence of defects in a fabricated product.
The aim of this study is thus to establish a discrete element modelling framework for the spreading of cohesive fine powders, where GPU (Graphic Processing Unit) computing [38][39][40] is combined with a novel stiffness scaling methodology [37] to enable a fast and reliable simulation. The characteristics of spread powder layer are qualitatively illustrated and quantitatively evaluated using a recently developed digital-based approach [41]. How the quality of spread powder layer varies with particle cohesion and the underlying mechanisms are systematically addressed. For the first time, particle clustering and pore characteristics of a spread powder layer are revealed. Since this study is mainly focused on the fundamental understanding of the role of particle cohesion in powder spreading, detailed comparisons between modelling and spreading experiments are thus not aimed here. In this work, particle cohesion is quantified via a dimensionless Bond number (Bo) which is defined as the ratio of maximum pull-off force between particles to particle gravity.

Model description
In the present study, all the simulations were conducted using an inhouse GPU-based DEM code (HiPPS). Its predictive capabilities have been demonstrated in many particulate systems, such as particle packing [37,38], compaction [38,42] and particle-fluid systems [39,40]. The DEM approach is thus only briefly outlined here. More details of force models can be found in our previous studies [37].
For a fine particle of radius R i , mass m i and moment of inertia I i , the translational and rotational motions are governed by the Newton's second law of motion. The resulting governing equations can be written as, where v i and ω i are the translational and rotational velocities, respectively; R ij is the vector pointing from the centre of the particle to the contact point with a magnitude of particle radius R i . F cn,ij , F dn,ij and F v,ik are the normal contact force, normal damping force and van der Waals force, respectively. F ct,ij and F dt,ij are the frictional force and the damping force in the tangential direction. The torques acting on the particle consists of a torque due to the tangential forces (i.e. the first term on the right side of Eq. (2)) and a rolling resistance torque due to the asymmetric distribution of the contact pressure, given as, where μ r is the rolling friction coefficient.
Normal contact between two particles is assumed to be viscoelastic while tangential contact force is governed by the Mindlin and Deresiewicz theory [43]. The normal and tangential force models are written as, The normal contact force F cn,ij (i.e. the first term on the right side of Eq. (4)) acts in the normal direction n pointing from the centre of particle i to that of the particle j, which is a function of the normal overlap, δ n , between two contacted particles.
. Y is particle Young's modulus and is the effective radius and m * is the effective mass. The normal damping force F dn,ij (i. e. the second term on the right side of Eq. (4)) is responsible for energy dissipation due to inelastic collisions, in which S n = 2E * ̅̅̅̅̅̅̅̅̅ ̅ R * δ n √ and β n is the normal damping coefficient. μ t is the coulomb friction coefficient. δ t is the particle displacement in the tangential direction and δ t = δ t /|δ t |. v t = v t /|v t | is the unit vector of the relative velocity in the tangential direction. δ t,max = μ t δ n (2− ν)/(2− 2ν) and β t is the tangential damping The van der Waals force between two spherical particles follows the Hamaker theory [44], given by, where A is the Hamaker constant. s is the interparticle separation distance and calculated as s = . A newly proposed stiffness scaling methodology [37] is adopted here to accelerate the spreading simulations of cohesive fine powders. Although it is a general practice to use reduced particle stiffness to accelerate DEM simulations; however, our recent study [37] showed that, for cohesive particles, there are many cases where the previously proposed scaling laws fail to preserve original particle behaviour, especially for contact-dominated systems, leading to under-predicted sliding and rolling resistances and a poorly resolved non-contact van der Waals interaction [37]. Different from the existing scaling approaches [36,45], the new scaling methodology consists of three key components: an established scaling law for contact adhesion (i.e. Eq. (7)), a modified normal contact force for the calculation of sliding and rolling resistances (i.e. Eq. (8)) and a new prediction-correction scheme for the estimation of non-contact van der Waals force (i.e. Eq. (9)), as summarised below, The superscripts R indicates a reduced particle property while O refers to the original particle property. F M cn is the modified normal contact force, F R v and F O v are the van der Waals forces calculated by the reduced and the original Hamaker constant, respectively. This modified calculation of sliding and rolling resistances assumes that transient response of normal contact force has little influence on the bulk behaviour of a contact-dominated system, as it dies out quickly due to viscous damping. ∆t DEM is the timestep used in the DEM calculation. v i0,n is the predicted velocity without accounting for the van der Waals interaction when two particles are not in contact while v i,n is the corrected velocity when this non-contact van der Waals interaction is considered. Detailed calculation of the corrected velocity and discussion on the applicability of this approach can be found in [37].

Simulation conditions and data sampling
Simulations were conducted in a rectangular box with dimensions 40d p × 600d p , d p is the particle diameter. The model configurations are illustrated in Fig. 1, which represents a generic powder spreading system, comprising a rigid spreading blade and a flat base substrate. Periodic boundary conditions are applied to the front and rear directions of the rectangular box (i.e. horizontal y direction as shown in Fig. 1) to eliminate the so-called wall effects. Idealized mono-sized spherical particles are used here, as spherical shape are typical for particles produced by gas or plasma atomization [24,46]. For simplicity, a typical  viscoelastic and adhesive contact force model is used here, which can be further improved by introducing more realistic force models, for example, considering plastic deformation [47] and particle sintering [48].
One complete simulation consists of three stages: die filling, heaping, and spreading. At the stage of die filling, a rigid vertical gate is placed 40d p in front of the blade, forming a confined space with a base area of 40d p × 40d p . First, particles were randomly generated within this confined space without initial overlaps and velocities. Secondly, particles were settled under the effect of gravity until forming a stable packing. At the heaping stage, the front gate was removed, which allows a stable heap to be formed in front of the blade. Finally, for the stage of spreading, the blade was lifted up to create a fixed blade clearance (i.e. 1.5d p ) and was translated at a constant velocity in the positive x direction, allowing particles to be deposited on the base substrate.
Particle properties are those typical of Ti-6Al-4V powder, following that used by Meier et al. [24]. The modelling parameters are summarised in Table 1. It is assumed that the rigid blade, front gate, and base substrate have the same properties as that of particles. A minimum cut-off distance s min is introduced in the calculation of van der Waals interactions, same as that used in the work of Parteli et al. [49]. A time step of 2 × 10 − 8 s is used for particles of reduced particle stiffness. Particle cohesion is specified via a dimensionless Bond number (Bo). In the case of van der Waals force, Bond number is defined as, which is the ratio of the maximum pull-off force to particle gravity. A wide range of particle cohesions are examined here by varying Bo from 0 to 400, which can be justified by the values found in the literature. For example, Bo is around 255 in the work of Nan et al. [31] for 45 µm stainless steel particles and cohesive force is two orders of magnitude larger than particle gravity for a plasma-atomized Ti-6Al-4V powder with a median particle diameter of 17 µm [24]. A total length of 500d p is spread for all the cases considered here. The first 100d p of a spread powder bed and 10d p immediately behind the spreading blade are excluded from data analysis. Data sampling was conducted within a length of 390d p , as shown in Fig. 1. A recently developed digital-based approach [41] is adopted to characterise the spread powder layer in term of packing density, surface profile and pore characteristics. The working space is discretised into voxels, with those voxels located within solid particles labelled as solid voxel. To aid data sampling, a Cartesian 2D sampling grid along the base surface is used to sample local properties. As the actual height of powder layer varies with both particle and process parameters and under certain circumstances it may be larger than the blade clearance. A gap-based packing density is thus defined here, which differs from the conventional definition but enables a consistent comparison among different cases and can quantify  the mass per unit area being deposited on the build plate. The local packing density at each grid node is calculated as a ratio of particle volume to the volume of a sampling domain with a height same as the blade clearance (i.e. h b ), (11) where s Ω represents size of the local sampling domain Ω. The local sampling domain has a parallelepiped shape, with a square base and a dimension of 3.0d p . The local packing density is simply calculated as a ratio of the number of solid voxels to the total number of voxels within the local sampling domain. The surface height at each grid node is determined as the maximum vertical coordinate (i.e. z-coordinate) of solid voxels found within the local sampling domain. It is of practical interest to characterise the region which can potentially lead to failure. Therefore, two types of pores are characterised for the spread powder layer, namely the so-called density pore and chamber pore [41]. The density pore is identified by thresholding the spatial distribution of local packing density, which allows us to identify those less populated areas in a spread powder layer. Here, as an example, the density threshold is set to 0.1. The concept of so-called chamber pore is adapted here to quantify the size of empty patches. The size of density pore is defined as its equivalent circular area diameter (i.e. the diameter of a circle of same area as the pore region) while the size of chamber pore is defined as diameter of the largest inscribe circle. Details of the calculation approach can be found in the previous study [41].

Validation of stiffness scaling
The validity of using reduced particle stiffness in DEM modelling is examined in both particle packing and spreading. Only two levels of particle cohesion (i.e. Bo = 5 and 50) are tested here, due to a high computational demand associated with the case of original particle stiffness. Fig. 2 shows the global packing densities as a function of stiffness scaling ratio for packed beds formed at the die filling stage. The global packing density is calculated as a ratio of solid volume to the volume of space occupied by the particles, where L x and L y are the box dimension in the horizontal directions (i.e. 40d p ), respectively. z min and z max are the minimum and maximum coordinates of a packed bed in the vertical direction. To reduce variability and surface effect, particles on the free surface of the packed bed are excluded from the calculation. The lower packing density seen with a higher Bo is qualitatively consistent with previous findings, as the interparticle cohesion resists relative motion between particles and hence leads to a loose packing [50,51]. For each Bo, the global packing density remains almost independent of scaling ratio, the variability is within 0.2% of that of the original particle stiffness in both cases, confirming the capability of the present model to achieve a stiffness-independent packing of cohesive fine powders. The validity of the stiffness scaling approach is further examined in the case of powder spreading. Fig. 3 compares local packing density distributions of the spread powder layer where different scaling ratios are adopted for the calculation. Scaling can be seen to have little influence on the packing density distribution. As shown in Fig. 4, the relative differences between different median local packing densities are < 3.6% for Bo = 5 and <5.3% for Bo = 50. The interquartile ranges are also relatively constant, the relative difference between different scaling is < 4.5% for a Bo = 5 and <13.7% at Bo = 50. This shows that the properties of the spread bed are well maintained when the particle stiffness is reduced by three orders of magnitude. To be computationally efficient, a scaling ratio of 0.001 is adopted in the present study. Fig. 5 shows the resulting packing structure of spread powder layers affected by particle cohesion. For illustration purposes, only a section of the layer is shown. Particles are coloured by coordination number (CN) (i.e. number of particle-to-particle contacts), where two particles are considered in contact when the distance between two centroids is smaller than a critical value (i.e. 1.005d p ) [50,51]. Spatial distribution of packing density is also examined, as it reflects the uniformity of mass being deposited on the substrate. Local packing density is sampled using a 2D grid, with a grid size of 1.0d p and a local sampling size of 3.0d p at each grid node.

Packing structure
Particle cohesion is seen to have a significant impact on the resulting packing structure. For non-cohesive particles (i.e. Bo = 0), both particle clustering and large empty patches are visible in the powder layer. Local packing densities can reach as high as 0.3, indicating a closely packed structure within these clusters, although it is significantly lower than that of a random close packing (0.62-0.64) due to a dominant monolayered packing structure. In contrast, no significant variation of packing density can be observed in the spreading direction for cohesive particles. A slight increase of particle cohesion (i.e. 0 < Bo < 5.0) results in a significant reduction in the degree of particle clustering. The uniformity of the spread layer continues to improve with further increase of Bo up to 50, as evidenced by a further reduction in the size of empty patch and a more homogeneous density distribution. However, a sparser layer starts to emerge when Bo is beyond 200. It is found that, within this Bo range, particle assembly in front of the blade largely maintain its shape integrity during spreading due to very strong interparticle attractive forces. The deposition process can thus be viewed as particles being stripped off the base of a solid-like assembly by a moving cohesive plate. Consequently, particles being deposited on the substrate are reduced considerably with increasing cohesion. This is different from the previous work [31] where large variation in the size of empty patches is still observable when Bo is larger than 200. For all the cases considered here, coordination numbers are quite small. This can be related to the dominant mono-layered structure because of a small blade gap (i.e. 1.5d p ). Moreover, a large proportion of particles has zero particle-to-particle contact, especially at Bo = 50. And chain-like structure starts to show when Bo is larger than 1.0 due to the increase of interparticle attractive forces, and consequently leading to a reduction in the coordination number. Fig. 6 shows the frequency distribution of coordination number in the spread layer. For all the considered cases, coordination number is smaller than 5. For non-cohesive particles, almost 45% of the particles have zero particle-to-particle contact. Particles with a higher CN accounts for a smaller proportion. However, an appreciable ratio of particles with a CN of 4 is still obtained in the spread layer of non-cohesive particles. As shown in Fig. 6(a), for Bo < 1.0, the increase of particle cohesion decreases the ratio of free particles (i.e. no contact with other particles) while increases the ratio of particles with one contact. The particles having more than two contacts are less affected. Within this Bo range, particles gain a high level of mobility after being released from the blade gap. The increase of interparticle cohesion makes particles more likely to stick together to form doublet, and thus leading to a decrease in the ratio of free particles and an increase of particles with one contact.
A more complicated behaviour is observed with a higher level of particle cohesion. As shown in Fig. 6(b), the ratio of free particles increases with Bo up to 50 and then drops to a local minimum at Bo = 200, beyond which it increases again. However, an opposite trend is observed for particles having contacts. For Bo between 1.0 and 50, this behaviour reflects the presence of two competing mechanisms: collision-induced particle clustering and particle-wall cohesion. Within this Bo range, particle-wall cohesion takes over as the dominant mechanism enhancing layer uniformity while suppressing the form of large clusters. Consequently, the ratio of free particles increases while the ratio of particles having more than one contacts starts to decrease. However, when Bo > 50, particles have very limited mobility after being discharged from the blade. The interparticle cohesion becomes so strong that particles contacted with other particles in front of the blade are more likely to stay in contact after being released from the blade gap. Because of the small blade clearance (1.5d p ) used in the present study, particles are more likely to form doublet, thus resulting in a significant ratio of particles with a single contact. When Bo is further increased (i.e. Bo > 200), a static assembly begins to form in front of the blade due to the large interparticle cohesive force. Particles are less likely to be escaped from the blade clearance, thus leading to an increased ratio of free particles and a decrease in the ratio of particles with contacts.
To view vertical structure of the spread layer, Fig. 7 shows frequency distribution of particle centres in the vertical direction. All the distributions show a significant peak at 0.5d p , suggesting a dominant monolayered packing structure. For Bo < 1.0, more than 99% of the particles are touching the base substrate and thus are not shown here. With increasing Bo, particles located within the range of between 0.5d p and 1.2d p start to increase. It is interesting to see that, because of strong interparticle cohesion, there is a small number of particles (< 5%) located above 1.0d p when Bo > 50. Fig. 8(a) shows statistical distribution of local packing density in the spread layer as a function of Bo. The median packing density decreases slightly with Bo up to 200, suggesting a limited influence with this Bo range. However, this is followed by a sharp drop when Bo > 200, as less particles being discharged from the blade gap due to the formation of a static assembly. To measure the dispersion of local packing density, Fig. 8(b) shows the ratio of IQR to the median as a function of Bond number, which has a similar role as coefficient of variation (or relative standard deviation) but is more robust to the presence of outliers. A large IQR\Median ratio indicates a high level of data dispersion. The results in Fig. 8(b) indicates the presence of four regimes regarding the variation of IQR\Median ratio, suggesting there are different dominant mechanisms of fine powder spreading. For the cases considered here, a

Particle clustering
Particle clustering or agglomeration is a pronounced feature of cohesive particulate systems. However, the degree of particle clustering has not been addressed before in a spread powder bed, despite its important practical consequences in additive manufacturing. For example, its presence has an impact on the absorption of laser energy, heat transfer and mass transfer during melting. Fig. 9 illustrates the presence of particle clusters in the spread layer, with each coloured by its orientation relative to the spreading direction. For clarity, particles without contacting with other particles (i.e. single particle) are not shown here. Cluster size is seen to decrease significantly with Bo up to Bo = 5.0 and then increases slightly with further increase of Bo. Here, the cluster size is defined as the number of particles consisting it. The shape of the cluster is also affected, shifting from large agglomerates to chain-like structures. Chain-like structure is seen to dominant when Bo > 5 and most clusters are found in the form of doublet. Cluster orientation also shows a dependence on particle cohesion. Clusters for a Bo between 5 and 200 shows a preferred orientation perpendicular to the spreading direction while clusters in other cases are more randomly orientated.
To quantitatively evaluate the extent of particle clustering, particlebased frequency distribution of cluster size is presented in Fig. 10. Most clusters are found in the form of doublet, partly due to the small blade clearance adopted here (i.e. 1.5d p ). For non-cohesive particles, the cluster size can reach as large as 136. For weakly cohesive particles (i.e. Bo ≤ 1), as shown in Fig. 10(a), the ratio increases with Bo for small clusters (i.e. < 4) while drops for larger clusters (i.e. > 8). This clearly shows that the presence of particle cohesion can suppress the formation of large agglomerates. For moderately cohesive particles (i.e. Bo ≤ 50), clustering is further suppressed with increasing particle cohesion. More than 95% of the particles form a cluster with a size smaller than 9. However, the clustering ratio increases again when Bo ≥ 50. It is found that, within this Bo range, particles tend to be adhered to the substrate immediately after being discharged from the blade gap, as the increase of particle cohesion leads to higher particle-wall adhesion. At the same time, the interparticle attractive force is significantly enhanced, making particles more likely to maintain their contact status after being deposited on the substrate and thus promoting the formation of chainlike structured cluster. Because of the small blade clearance, chain-like structure is still the dominant form. For larger clusters (i.e. > 2), a further decrease of the clustering ratio can be observed for strongly cohesive particles (i.e. Bo ≥ 200). This is because interparticle attractive force becomes so strong that a static particle assembly formed in front of the blade, hindering the release of particles, and consequently reducing the formation of large clusters. Fig. 11 shows the particle-based orientation distribution of clusters in the spread powder layer. Here, clusters are classified into three groups: small cluster with a size of 2, medium cluster with a size of 3 or 4 and large clusters with a size larger than 5, to examine the correlation between cluster size and its orientation. As shown in Fig. 11(a), clusters of non-cohesive particles have no orientation preference, irrespective of the cluster sizes. For weakly cohesive particles (i.e. Bo ≤ 1.0), small clusters (doublets) tend to align in the spreading direction while no preferential alignment is observed for medium and large clusters. The increase of particle cohesion slightly increases the clustering ratio of doublets while not affecting much their orientations. In contrast, as shown in Fig. 11(b), clusters are more likely to be aligned perpendicular to the spreading direction when particle cohesion is further enhanced. It is especially the case when Bo ≥ 50. This again can be explained as the loss of mobility after particles being discharged and the increase of interparticle attractive forces. Consequently, small clusters and clusters that are aligned with the blade have a higher survival rate when escaping from the blade gap. Large clusters and clusters with other orientations are more likely to be disintegrated by transient jamming that occurs in front of the blade. However, for strongly cohesive particles (i.e. Bo > 200), large interparticle attractive force leads to the formation of a rigid assembly, and thus reducing particles being discharged and increasing the survival rates of clusters with other orientations, as particles are essentially being stripped off the base of the assembly. That is why the orientation of clusters is slightly more random at Bo = 400 compared to that of Bo = 200. Since large clusters are more difficult to survive compared to small clusters, the clusters with random orientations are dominated by the form of doublet. Fig. 12 shows the effect of particle cohesion on the surface profile of a spread powder layer, where surface height is sampled with the aid of a 2D grid, where the grid size is set to 1.0d p and a local sampling size of 1.0d p is used. For weakly cohesive particles (i.e. Bo ≤ 1.0), more than 80% of the spread powder layer has a surface height between 0.75d p and 1.0d p , indicating a mono-layered packing structure. In contrast, multiple layered structure starts to emerge when Bo ≥ 5. The value of zero corresponds to those empty patches observed in Fig. 5. Its ratio decreases with particle cohesion up to Bo = 50, beyond which it increases again. For cohesive particles, the surface height can be larger than 1.25d p but is still limited by the blade clearance (i.e. 1.5d p ). The ratio for 0.75d p < z ij ≤ 1.0d p decreases with particle cohesion while the ratio for z ij > 1.0d p first increases up to Bo = 50 and then decreases.

Surface profile
To quantify the effect of particle cohesion on the surface profile, surface roughness (Ra) is presented in Fig. 13 as a function of Bo, which is calculated as arithmetical mean deviation of the surface height. A smaller value of Ra indicates a better surface homogeneity. Variation of surface roughness is found to be consistent with the distribution of local packing density (i.e. Fig. 8(b)), where four different regimes can also be observed. For weakly cohesive particles (i.e. Bo ≤ 1), surface homogeneity increases with particle cohesion, approaching a value around 0.12d p . With further increase of particle cohesion, a local minimum is seen at Bo = 50, beyond which Ra increases with Bo at a decreasing rate up to Bo = 200. For strongly cohesive particles (i.e. Bo > 200), Ra increases almost linearly with Bo. This is different from the observations of Wang et al. [25] where the surface roughness was found to increase with Hamaker constant.

Pore characterization
As shown in Fig. 14, for weakly cohesive powders (i.e. Bo < 1), density pore tends to have an irregular shape and possesses a rather random orientation. With increasing particle cohesion, the resulting pore becomes more elongated and aligns more perpendicular to the spreading direction. The number of pores and the corresponding pore size are also affected. The density pore size first decreases until Bo = 50 and then increases with particle cohesion. An opposite trend is observed for the number of density pores. At Bo = 400, a single large pore that covers most of the spread powder layer is obtained, corresponding to a very sparse particle distribution. Like that of the density pores, the size of chamber pore is also seen a decrease followed by an increase with Bo. For cohesive powders, chamber pores are distributed quite uniformly, with no clear localization observed in the spreading direction while large size variations of the chamber pore can be observed in the spread layer of non-cohesive particles.
A quantitative description of the role of particle cohesion can be revealed by the cumulative distribution of pore size, as shown in Fig. 15. For non-cohesive particles, the size of density pore can reach 24d p , which is significantly larger than those of cohesive particles. The slope of distribution first increases and then decreases with particle cohesion, peaking at Bo = 50. For the weakly cohesive powders (Bo < 1), a significant increase of the slope can be observed, suggesting a slight addition of particle cohesion can significantly reduce the pore size. Consistent with the distributions of density pore, the slope of the distribution of chamber pore also peaks at Bo = 50. For Bo < 200, most of the chamber pores are smaller than 4.0d p and the distributions almost follow a same trend till a pore size of 1.0d p , indicating a limited influence of particle cohesion on the ratio of small chamber pores.
To further explore the role of particle cohesion, both types of pores are classified into six size groups. As shown in Fig. 16(a), the increase of particle cohesion in the range of Bo ≤ 1 mainly reduces the ratio of large pores (>5d p ) while further increase of the particle cohesion leads to a change of the dominant small pore size, from 2d p to 3d p for Bo = 5, to 1d p -2d p for Bo = 50 and back to 2d p -3d p for Bo = 100. For Bo ≥ 200, the spread layer is again dominated by large pores (>5d p ). For the chamber pores, as shown in Fig. 16(b), the groups of large pores (>3d p ) follow a same trend that increasing particle cohesion leads to a decrease followed by an increase of the ratio. The very small chamber pores (<1.0d p ) remains relatively stable up to Bo = 50, beyond which pores smaller than 2d p are seen a sharp drop in the ratio. In contract, variation of the medium sized chamber pores (2d p -3d p ) is more sensitive to particle cohesion, indicating a dependence that is in line with other descriptors discussed before.

Mechanisms affecting powder spreading
In terms of the mechanisms underlying inhomogeneous particle distribution, it can be understood from the competing effects of particle mobility, particle-particle, and particle-wall cohesions. Animations of the spreading process can be found in the supplementary, in which particles are coloured by the velocity in the spreading direction. After being released from the blade clearance, particles gain a momentum to continue moving forward. The resistances to particle's motion are mainly composed of sliding and rolling resistances. Compared to rolling, sliding is the dominant mode of motion for a cohesive particle. A dimensionless inertial length, I μ , characterising the distance that a particle can travel before friction brings to a stop can be estimated as, where V s is the blade speed which is characteristic of the particle velocity immediately after being released from the blade gap. μ t is the sliding friction coefficient. F v is the pull-off force required to separate a particles from the base substrate. The μ t (mg + F v ) denominator is an estimate of the frictional force on the particle once it has left the blade gap. If particle-particle and particle-wall sliding friction coefficients are different, this term can be adjusted, and a modified inertial length can be defined. The dimensionless inertial length, I μ , is 8.49 for non-cohesive particles and 4.24 when Bo = 1.0. Within this Bo range (i.e. Bo ≤ 1), particles gain a high level of mobility after being discharged from the blade gap. They continue to move forward and collide with each other, dissipating kinetic energy and consequently leading to the formation of large clusters. Collision-induced energy dissipation can be thus attributed as the dominant mechanism responsible for the inhomogeneity of a spread powder layer, which is also one of the mechanisms giving rise to particle clustering in gas-solid fluidization [52]. For weakly cohesive particles (i. e. Bo ≤ 1), the increase of particle cohesion mainly reduces the cluster size and thus leading to an improved layer homogeneity.
When particle cohesion becomes more significant, I μ quickly drops to 1.41 for Bo = 5.0 and 0.16 for Bo = 50. Particles have very limited mobility within this range of Bo. Particles adhere to the substrate shortly after being released from the blade gap. Collision-induced particle clustering is thus significantly suppressed. Particle-wall adhesion takes over as the dominant mechanism which can enhance the uniformity of a spread layer. Further increase of particle cohesion (i.e. 50 < Bo ≤ 200), I μ drops to 0.04 at Bo = 200. Particles adhere to the base substrate immediately after being released from the blade gap. Chain-like cluster structure starts to dominate in the spread powder layer, which can be attributed to the increase of interparticle cohesion. However, when Bo > 200, It becomes very difficult for the particles to escape from the blade gap due to the formation of a static assembly in front of the blade and consequently less particles are deposited on the layer and particle distribution becomes sparser. In summary, the increase of particle cohesion leads to four competing mechanisms affecting the structure of a spread powder layer: 1) an inhomogeneous distribution for non-cohesive and weakly cohesive powders (i.e. Bo ≤ 1) due to collision-induced particle clustering, 2) an enhanced layer homogeneity for moderately cohesive powders (i.e. 1 < Bo ≤ 50) due to the increase of particle-wall adhesion, 3) a reduction in layer homogeneity when further increasing interparticle cohesion (i.e. 50 < Bo ≤ 200) due to the formation of chain-like clusters by strong interparticle attractive forces, and 4) particles being stripped off the base of a static assembly for strongly cohesive powders (i.e. Bo > 200).
It is worthwhile to mention that particle sliding will certainly be hindered if spreading over a rough surface formed by pre-melted/ sintered powder. Interlocking with the substrate may also occur, but the extent of which depends on surface morphology, particle size, shape, cohesiveness, and subjected stress. For SLM, the surface roughness can be in the order of 10 µm, spreading of small particles are thus expected to be affected more than that of larger particles. Moreover, weakly cohesive particles can be more affected than that of strongly cohesive particles, due to a larger particle mobility after being discharged from the blade gap. As the analysis performed in this paper is based on a small blade clearance, further study on the interaction between particle cohesion and blade clearance are worthwhile. Increasing the blade clearance allows more particles to be discharged, therefore, the empty patches observed in the present study may vanish. For weakly cohesive particles, particle rearrangement can be enhanced so that layer homogeneity and mean packing density can be improved. However, a more realistic contact force model should not have a significant influence on the results presented here, as cohesion force plays a dominant role in determining the sliding and rolling resistances for cohesive fine particles.

Conclusions
In the present study, a discrete element modelling framework was established, which combines a stiffness scaling methodology with GPU computing to enable fast and reliable particle-scale simulations of the spreading of cohesive fine powders in powder-bed-based additive manufacturing. Together with a digital-based approach for the characterisation of spread powder layer, a direct link between microscopic particle properties and layer characteristics was established. The role of particle cohesion in the spreading of spherical particles was systematically addressed, where the strength of particle cohesion is quantified by a dimensionless Bond number. The main findings are summarized as follows, • Particle stiffness can be safely reduced three orders of magnitude for both packing and spreading of cohesive fine powders, confirming the validity of the present approach in handling reduced particle stiffness. • Particle cohesion has a strong impact on the local structure of a spread powder layer. It prohibits the formation of large agglomerates for weakly cohesive powders while promotes the formation of chainlike structures for strongly cohesive powders. Increasing particle cohesion first decreases and then increases the cluster size, and  making these clusters aligned increasingly perpendicular to the spreading direction. • The median of local packing density remains relatively stable over a large range of Bo, which decreases slightly with particle cohesion up to Bo = 200 beyond which a sharp drop is observed. • The layer characteristics shows a consistent and complex dependence on particle cohesion, as captured by several metrics in terms of local packing density, surface roughness and pore size. For the cases considered here, an optimal level of particle cohesion is found around Bo = 50, showing the smallest variations of local packing density, surface roughness and pore sizes. • Four competing mechanisms are identified regarding the role of particle cohesion on powder spreading, which can be partly interpreted by a dimensionless inertial length.
Although this study presents an idealized spreading of cohesive fine powders, it reveals a more complex spreading behaviour than that might be anticipated. The modelling framework developed in this work represents a critical step towards understanding powder spreading in more realistic systems. Future work includes i) extension to more realistic boundary conditions, including the morphology of substrate and operating temperature, ii) quantitative model validation on the quality of spread powder layer and iii) the influence of material properties and their interaction with process conditions.

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