Structural fluctuations in thin cohesive particle layers in powder-based additive manufacturing

Producing dense and homogeneous powder layers with smooth free surface is challenging in additive manufacturing, as interparticle cohesion can strongly affect the powder packing structure and therefore influence the quality of the end product. We use the Discrete Element Method to simulate the spreading process of spherical powders and examine how cohesion influences the characteristics of the packing structure with a focus on the fluctuation of the local morphology. As cohesion increases, the overall packing density decreases, and the free surface roughness increases, which is calculated from digitized surface height distributions. Local structural fluctuations for both quantities are examined through the local packing anisotropy on the particle scale, obtained from Vorono\"{\i} tessellation. The distributions of these particle-level metrics quantify the increasingly heterogeneous packing structure with clustering and changing surface morphology.


Introduction
Powder-based additive manufacturing techniques, like powder bed fusion, have garnered considerable interest [1,2,3] for their ability to facilitate rapid prototyping and the production of highly customizable parts.These methods enable efficient manufacturing by minimizing the need for material removal and extensive support structures, which in turn reduces production time and material waste.However, the quality and efficiency of powder-based techniques are far from ideal.Non-uniform powder packing during spreading is one of the major issues that limit the range of available powder materials and impair printing quality.Various types of structural defects in the deposited powder layer have been observed, which strongly correlate to defects in sintered parts [4,5,6].Since commonly used particle sizes are far below 100 µm, cohesion between particles can impair the spreading and deteriorate the quality of the powder layer through reduced powder flowability and cohesion-induced powder clustering.Understanding the influence of cohesion on spreading requires detailed measurement of the packing structure under various levels of cohesion, which is expensive and difficult to obtain experimentally [7,8,9,10].Characterizing a thin particle layer is also challenging, as most of the existing metrics are meant for bulk characterization.
One prominent method utilized to understand and design the powder spreading process is the Discrete Element Method (DEM), which is a particle-based simulation technique that computes particle trajectories from the interaction forces.Various DEM-based studies have investigated powder spreading with the aim of improving the quality of the powder layer [11,12,13,14,15,16,17,18]. Simulations can reflect behaviors of real powders [19,20] during spreading as they can be calibrated by experiments of powder flowability [21,13,22,19], which allows detailed studies of the influence of process parameters and powder properties on spreading.For example, Parteli and Pöschel [23] showed in simulations that a fast spreading process increases the surface roughness of a cohesive powder layer for roller spreading.Nasato et al. [24] found that small frequency and amplitude of a vibrating recoater lead to low powder bed porosity.Non-spherical powders with realistic particle shapes were also considered when investigating how the recoating velocity influences the bed porosity [25,12].Shaheen et al. [13] showed that powder layer defects are more likely to occur with higher particle rolling and sliding friction.
In these studies, the prerequisite of establishing the relation between process and material parameters and the layer quality is a detailed and informative characterization of the packing structure, which can be challenging for cohesive particles due to effects like clustering.While the global packing density is informative and widely used, it does not contain information of how the particles are spatially arranged.Therefore, the spatial fluctuation of the packing structure is also important, especially for highly cohesive powders where the packing tends to be heterogeneous [3].To this end, local density is often calculated using binning and coarse-graining where the averaging length must be chosen [18,13].Metrics based on the Voronoï cell volume can also be used, which does not require hand-picking an averaging length scale.For example, Phua et al. used Voronoï-tessellation of particles in 3D to calculate the average packing fraction of powder layers [26].However, examining the global distribution of the Voronoï still does not offer the complete picture of how density fluctuates.Here, we adopt a Voronoï-tessellation based method [27] to quantify local structural anisotropy, which is an inherent property of non-crystalline packing of particles and is associated with critical mechanical properties in disordered packings, such as jamming [27], plasticity [28], and shear band formation [29,30,31].This method does not require choosing a density threshold to identify voids and it yields a meaningful distribution of local anisotropy in a deposited powder layer, based on which the hetereogeneity of the packing can be quantified.
The surface roughness of the deposited powder also plays a crucial role in determining the functionality and aesthetics of the final product.Achieving the desired surface finish is essential for optimizing performance and ensuring consistent product quality.Surface roughness, similar to density, is also influenced by the interplay between process parameters [25,16,32,33] and material properties such as cohesion, particle size distribution [14] and particle shape [12,34].In particular, cohesion strongly influences the surface roughness during spreading.The powder bed surface roughness increases with cohesion due to powder agglomeration and particle removal caused by particle-to-blade cohesion during spreading [18,20].In DEM simulations, the surface roughness is typically evaluated by measuring the local surface height determined by the maximum vertical coordinate of the powder bed and monitoring its spatial variation.Using this variation as a metric of uniformity, surface roughness is calculated as the mean deviation of surface height from the powder bed average height [14,18].Experimentally, the surface height can be determined using optical 3D digital microscopy [32] or high-speed laser profilometry [33].Surface roughness can be quantified either using planar profile measurements in two dimensions [12,25] or areal measurements in three dimensions [18,35].
While metrics like the standard deviation of the global surface height distribution is informative, it does not offer a complete description of the surface profile.In this study, we evaluate the skewness and the kurtosis of the height distribution calculated using an efficient digitization method [35].These characteristics offer further insight into the presence of local outliers in surface roughness and the extent to which they deviate from the mean surface plane.We also address the problem that for a given set of surface height values, the distribution cannot well describe the local fluctuations because a spatial rearrangement of the height values does not change the distribution.This is similar to the aforementioned problem that the global packing density cannot sufficiently describe the heterogeneity of the packing.To this end, we quantify the spatial fluctuations of the free surface height of the powder layer through a coarse-graining approach to calculate the squared local spatial gradient of the Voronoï cell-averaged height.This quantity again yields a meaningful distribution that can be described by a single parameter, quantifying the height fluctuations.

Numerical Setup
We employ DEM to obtain particle-scale information on powder layers created by a spreading process, using MercuryDPM [36].The simulation setup is shown in Figure 1.The powder is spread by a blade tool, moving at constant velocity v T along the spreading direction, x [37,38].We simulate a small slice of the powder bed of length 10 mm in the x-direction and width 1 mm in the lateral y direction where periodic boundary conditions are applied.For the subsequent analysis, we consider the range 0 ≤ x ≤ 7 mm].It is assumed that the substrate is flat and the coefficient of friction between the wall and the particles is equal to that of the particle-particle interaction.A log-normal particle size distribution is considered with mean particle diameter D 50 = 37 µm, D 10 = 24 µm, and D 90 = 56 µm.The particles are initially generated in front of the spreader tool at (x, y, z) ∈ [0.5, 2.5] mm ×[0, 1] mm ×[0, h] as shown in Figure 1(a), filling a total bulk particle volume of 0.75 mm 3 which is sufficient to create a powder layer of 10 mm length, 1 mm width, and tool gap H = 100 µm, where the tool gap is defined as the gap between the base of the blade and the substrate, as shown in Figure 1(b).The spreading process starts at a constant velocity v T = 10 mm/s with initially all particles at rest.It ends when the blade arrives at the end after 1.2 s, and the simulation ends at time 1.5 s once the system is relaxed again, i.e., when the kinetic energy is sufficiently low.

Contact models 2.2.1. Hertz-Mindlin visco-elastic contact model
The visco-elastic Hertz-Mindlin contact model (no-slip solution) [39,40] is employed to calculate the normal and tangential elastic contact forces between particles, respectively.The normal force for the Hertz visco-elastic model is given as where ξ = R i + R j − |⃗ r i − ⃗ r j | is the compression of two interacting particles i, j of radii R i and R j at positions ⃗ r i and ⃗ r j and ⃗ e n = (⃗ r i − ⃗ r j )/|⃗ r i − ⃗ r j | is the normal unit vector, A n = 5 × 10 −6 s is the normal dissipative parameter, calculated as in [41], considering a coefficient of restitution of 0.4 for the characteristic blade velocity 10 mm/s and with the effective radius R * .The effective elastic modulus, depends on the elastic moduli and the Poisson ratio of the material of particles i and j.
We model the tangential viscoelastic forces following the no-slip solution of Mindlin [42] for the elastic part and Parteli and Pöschel [23] for the tangential dissipative constant A t ≈ 2A n E * , which are capped by the static friction force between two particles.The tangential force is given by with the friction coefficient, µ, the effective shear modulus which for particles of identical material simplifies to G * = 4G 2−ν , and the tangential relative displacement of the particles, ds.

Non-linear cohesive model
To simulate particle cohesion, we incorporated adhesive forces described by the Johnson-Kendall-Roberts model [43] (JKR) and attractive forces using a model for non-bonded van der Waals interactions [8].The JKR adhesive force is computed as where γ is the surface energy density and a is the contact radius related to deformation, calculated using The maximum interaction distance at which the contact breaks under tension is given by The non-bonded van der Waals attractive force [8,44,45] reads where D min = 1.65 Å is a parameter introduced to avoid a singularity [8], D max is the maximum interaction distance of the van der Waals interaction, which is set as 1 µm [8] and A H is the Hamaker constant which relates to the surface energy density via

Material Parameters
The powder spreading process is simulated considering a metallic Ti-6Al-4V powder.The material and simulation parameters can be found in Table 1.According to the experimental measurement of the angle of repose of 41 • and matched with the simulation results of Meier et al. [19], the surface energy of Ti-6Al-4V is 0.1 mJ/m 2 .To study the effect of particle cohesion on the powder quality, we simulate the powder spreading process for varying surface energy γ from 0 to 0.5 mJ/m 2 in steps of 0.05 mJ/m 2 .In general, for cohesive bonds, the surface energy γ quantifies the energy associated with disrupting a bond between cohesive particles to create surface.Thus, varying γ is a meaningful representation of varying cohesion intensity between neighboring particles.We introduce the Bond number 11) to characterize the ratio between interparticle cohesion and gravity, where D 50 = 37 µm.Note that Bo = 54.5 corresponds to the surface energy of the Ti-6Al-4V powder with γ = 0.1 mJ/m 2 [19].For the given particle and material parameters, the values of γ ∈ {0, 0.05, 0. surface energy (γ) mJ/m 2 0 − 0.5 surface energy of Ti-6Al-4V mJ/m 2 0.1 interaction distance (D max ) µm 1

Local density characterization of the powder layer
The quality of the produced powder layer is closely related to the packing density of the particles prior to sintering [20].The density of the layer can be quantified by the ratio of the volume occupied by particles and the total volume, ϕ = V solid /V total .For sufficiently small V total , ϕ can be considered as a local variable.Low packing fraction values indicate loose structures that are prone to defects in the final product.In general, a high packing fraction is desirable for high product quality.
The packing density is calculated locally for subsections of each layer to provide information about the spatial variability of voids throughout the layer.To this end, the local packing fraction is calculated for horizontal strips across the spreading distance x, of fixed width equal to 1 mm.The strip size is chosen to be sufficiently large so that it contains a representative number of particles and voids for the calculation of the packing fraction.The density calculations are performed using YADE, where a dedicated algorithm exists for density calculations [46].Alternative, high-resolution techniques have been proposed to calculate the density of granular packings based on the exact partial intersection volume between spheres and mesh elements [47].
Powder spreading leads to unstructured, inhomogeneous layers of material with spatially varying packing characteristics.This is the motivation behind calculating the packing density locally, for subsections of the layer along the spreading direction, aiming to explore the degree of density inhomogeneity within each layer.Figure 2 shows values of the packing fraction for various values of cohesion as a function of Bo.Evidently, cohesive materials lead to loose packings of the powder material, which is in agreement with previous observations [20].For increasing cohesion from Bo = 0 to Bo = 272.3, the mean packing fraction reduces by nearly 60%, from ϕ ≈ 0.60 to ϕ ≈ 0.25, while the scattering of the values also increases slightly.The reduction of the average packing fraction is gradual for increasing Bond number (within the studied range of values), i.e., no sudden transitions are observed between layers made of powders with similar Bond numbers.The data points in Figure 2 are colored according to their distance, x, from the starting point of spreading, where a clear trend is not observed, indicating that the degree of scattering does not correlate with the spreading distance, x.

Heterogeneous packing of cohesive particles
Figure 3(a) shows an example of the deposited layer of highly cohesive particles (Bo = 190.6).We observe a heterogeneous structure comprising regions of dense and loose packing.At the particle scale, the density of neighbors surrounding each particle can be highly anisotropic, which could have important implications for subsequent processes such as heat transfer and phase change.Although such spatial fluctuations can be reflected by bin-averaged density, as done in Figure 2, and the degree of fluctuation of densities at different locations depends on a manually chosen bin sizes.To avoid the need to manually specifying sampling length scales, we adapt a method of characterizing the structural heterogeneity using a particle-level measurement based on the packing anisotropy from the Voronoï tessellation [27].

Anisotropy Vector and Divergence
The measure of the local structural anisotropy was developed by Rieser et al. from the observation that the center of a particle deviates from the centroid of its Voronoï cell in a disordered packing [28].Any two particles with a shared Voronoï cell face are defined as neighbors; from this, a Delaunay triangulation is generated by connecting groups of three mutual neighbors into triangles.Figure 3(a) includes a schematic illustration of the Voronoï tessellation calculated based on the projections of the particle positions on the xy-plane, which is plotted below the particle packing with the Delaunay triangle k and the anisotropy vector ⃗ C pointing from particle centers to corresponding Voronoï cell centroids.For a triangle representing a densely occupied area (overpacked, like the one depicted), all the ⃗ C vectors point outward.For a triangle representing a void (underpacked), the vectors point inward.We quantify the extent to which the vectors ⃗ C at the three vertices of a Delauney triangle point inward or outward by calculating the divergence of the vectors of a triangle k with area A k .This is calculated based on the concept of constant strain triangle in finite element analysis [48].The local structural anisotropy, Q k , calculated from the divergence, defined as where Ā is the average of all A k within the packing.By construction, Q k is dimensionless with a mean near zero.It is sensitive to the local structural anisotropy and has a geometrical significance: positive (negative) values correspond to overpacked (underpacked) regions.The challenging aspect here is to extend the 2D calculation to a quasi-2D thin free-surface layer with a height of two to three particle diameters, which is typical in powder spreading.Since the deposited layer is thin and many regions only contain a single layer of particles, as shown in Figure 3, we use the 2D projections of the particles on the xy plane for calculating the anisotropy.This simplification could lead to contributions of highly positive Q k , especially in non-cohesive packing, due to vertically aligned particles in a quasi-2D layer.To see the influence of such scenarios, we scale Q k with the ratio of the projected area of overlapping particles, A p , on the xy plane and the sum of the real area of the particles, A r = πr 2 i .The scaled Q ′ k is given as follows: For highly dense packings of vertically aligned particles, A p /A r < 1, and the anisotropy is reduced.For dilute packings, A p /A r = 1, and thus Q ′ k coincides with the original definition without correction.The distribution of Q k is a strong structural indicator that is associated with important mechanical properties of a disordered packing, including jamming and shear band formation [27,29,30,31].In dense and homogeneous regions, Q k fluctuates randomly, and a peak in the Q k distribution around Q k = 0 is typically observed [27].In heterogeneous regions with high anisotropy, the highly positive and negative Q k values show up together on the tails of the distribution, making them deviate from Gaussian.Figure 4(a) shows the distribution of Q k for different Bo.The majority of the Q k resides in the region around zero.It can be fitted to a Gaussian distribution using values between Qk −0.5 to Qk +0.5 for each data set (solid curves), where Qk is the mean of the distribution.For lower Bo, we observe a consistent slope of the distribution throughout the range Q k < 0, which is an indication of a homogeneous structure throughout the packing.In contrast, a transition of the slope of the distribution at Q k = −1 is clearly distinct for higher Bo, suggesting the coexistence of dense homogeneous regions and dilute heterogeneous regions for highly cohesive materials: for Q k < −1, the distribution deviates from Gaussian and becomes exponential-like for higher Bo.This exponential tail corresponds to the existence of highly underpacked sites distributed sparsely in the packing, as seen in Figure 3(c); for −1 < Q k < 1, the distribution is narrower with increasing Bo, indicating the existence of locally homogeneous packing.This variation in structure is also evident in cohesive systems shown in the experimental studies by Xiao et al. [29].Such a variation is not observed for non-cohesive experimental particle systems in Harrington et al. [31] for disordered particle packings.

Divergence Fields and Distributions
To quantify the difference in packing heterogeneity for different Bo, we show the standard deviation and the skewness of the Q k distribution in Figure 4(b).The standard deviation reflects the portion of highly anisotropic sites (triangles) in a packing.The skewness roughly compares the degree of anisotropy of loosely packed sites to densely packed sites.For higher cohesion, particles can sustain more voids during spreading and encounter higher local anisotropy, leading to a more heterogeneous overall packing structure.As a result, the standard deviation increases with Bo, which reflects the difference seen in Figure 3 in a quantitative way.The skewness decreases with Bo, which reflects the growing tail at the negative end of the Q k distribution.This corresponds to the fact that the void sites not only grow larger in number, but also have larger sizes at higher Bo.We compared the standard deviations and the skewness of the distributions with the anisotropy calculated using Equation 12 and Equation 13 for Q ′ k , respectively.The divergence Q k from the condition in Equation 12gives slightly higher values for both the standard deviation and the skewness but qualitatively shows the same behavior as Q ′ k .

Digitized free surface height characterization
Measuring the surface roughness of the powder layer is of interest, as it is related to the roughness of the final product, while distinct rough features are areas prone to become the source of defects.As discussed in the previous sections, packing density and structural packing anisotropy are integral elements in assessing the quality of the finished part.However, they do not provide information on the irregularity of the surface texture of the powder layer.To this end, it is useful to characterize the surface roughness of the produced layers corresponding to different cohesion values.Various aspects of roughness can be characterized via the calculation of independent quantitative indices that provide diverse morphological information on the layer's topography.
Previous work focused on the two-dimensional characterization of surface roughness features of powder layers, using indices that correspond to planar rough profiles, usually taken as representative of the real rough profile [23,24].Here, the three-dimensional surface profile of each powder layer was reconstructed for each layer after spreading is completed.The particles located near the top of the powder layer were identified, and points on their surface were calculated using a regular sampling grid [35], giving a surface height profile, z s , similar to Meier et al. [18].Note that if the surface of the substrate is directly exposed at a sampling point, a value of zero is recorded.Figure 5 shows the surface heights of the 11 studied powder layers of different cohesion.Interestingly, the surface height reaches a maximum value of up to 150 µm for larger Bond numbers, which is larger than the gap height of 100 µm, shown in Figure 1(b).This typically occurs for fine cohesive powders, due to decreased flowability of the powders when the cohesion effects dominate gravity and inertia, resulting in the formation of agglomerates with internal cavities and irregular surface profiles [17,49].The surface height of the cohesionless powder layer (Bo = 0) does not exceed the gap height.Figure 6 shows three example distributions of the measured surface height, which shows that the distribution widens as Bo increases.However, the shapes of the height distributions are rather complicated, and require many parameters to describe as listed in the following subsection.
A spike at z s = 0 exists for all three Bo, which corresponds to the exposed substrate surface.

Roughness characterization using height distributions
The current state-of-the-art for characterizing rough surfaces, as outlined in ISO 25178 [50], calculates roughness indices based on the surfaces of real, three-dimensional texture profiles.Using distributions of z s , the surface roughness is characterized in terms of arithmetic mean height (S a ), root mean square height (S q ), skewness (S sk ) and kurtosis (S ku ).The height deviation of the surface roughness from a mean surface height of the entire layer is used, which is z m .We show these roughness parameters in Figure 7 for increasing Bond number values, where the points are colored according to their distance x from the start of spreading, and along the spreading direction.A clear correlation was not found between any of the surface roughness parameters and their distance from the initial spreading position.We next discuss the significance of the roughness parameters individually.The arithmetic mean height is calculated as: where z = z s − z m is the height of a point on the layer surface, measured from the plane of mean surface height, z s is the measured free surface height, x and y the horizontal coordinates of the point along and transversely the spreading direction, and A the area occupied by the layer.It becomes evident in Figure 7(a) that for powders of increasing cohesion, the arithmetic mean height increases almost linearly with the Bond number up to values of Bo = 217.9.This indicates that cohesive powders lead to higher deviations from the mean height, and to rougher surface texture profiles.Also, the scatter of measurements increases slightly for larger Bond numbers, which points to the conclusion that more cohesive powders feature more heterogeneous profiles, with taller peaks and deeper valleys.To further validate this trend, the root mean square height is calculated as: It can be seen in Figure 7(b) that the root mean square height presents the same general trend as the scattering of the arithmetic mean, where more cohesive powders form layers with more heterogeneous height distributions.This is in agreement with findings from the literature [18,20].These two measures of the average height of the powder surface texture are informative regarding the extent of the roughness, but provide no information about their morphology.To this end, the skewness and kurtosis of the surface height profiles are examined.The height skewness is calculated as:  The solid horizontal lines note the global values of the surface roughness parameters for each layer.The sample points are colored according to their distance x from the starting point of spreading.The dashed line for S sk = 0 marks the threshold between profiles where most rough features appear above the mean plane (S sk < 0) and below it (S sk > 0).The dashed line for S ku = 3 marks the threshold between rough profiles with rounded peaks (S ku < 3) and with sharp ones (S ku > 3).
Skewness is a measure of the asymmetry of the layer height distribution around the mean plane.Negative skewness values (S sk < 0) indicate that the height distribution is skewed above the mean height plane, with a few deep valleys, zero skewness values (S sk = 0) correspond to a symmetric surface, where peaks and valleys occupy the same amount of surface in average, while positive values (S sk > 0) indicate that the height distribution is skewed below the mean height plane, with a few tall peaks.Figure 7(c) shows a monotonically increasing trend of skewness with increasing Bond number values, where powders with lower cohesion (Bo < 190.6) demonstrate negative skewness, with average cohesion (Bo ≈ 190.6) nearly zero skewness and with higher cohesion (Bo > 190.6) positive skewness values.The height kurtosis is calculated as: Like skewness, kurtosis describes a particular morphological aspect of the surface height distribution.Skewness is a metric of whether most of the rough profile is positioned above or below the mean height plane.Kurtosis provides information on the average shape of the surface texture asperities, and can be seen as a probability density sharpness of the rough features.Low kurtosis values (S ku < 3) indicate platykurtic surface texture profiles of well-rounded asperities presenting short tails, zero kurtosis values (S ku ≈ 3) correspond to mesokurtic profiles of Gaussian-like asperities characterized by medium-sized tails, while high kurtosis values (S ku > 3) correspond to leptokurtic surface profiles, with sharp, spike-like asperity characteristics presenting long tails.
Figure 7(d) shows the kurtosis values for the various powder layers of varying cohesion, where the parameter shows a non-monotonous, mostly declining trend for increasing Bond number.It becomes evident that the powder layer corresponding to zero cohesion (Bo = 0) features high kurtosis values (S ku > 3), while layers made of cohesive powders feature lower kurtosis values (S ku < 3).For the higher end of the studied cohesion levels (Bo > 217.9) kurtosis shows a mild increasing trend, which is however characterized by a high degree of scatter, making a further interpretation challenging.
Combining the observations of all roughness parameters for the studied powder layers of varying Bond number, it can be that increasing cohesion leads to powder layers characterized by increased roughness, where the layer lies mostly below its average height, and presents a few, rounded peaks.For less cohesive powders, the corresponding layers are characterized by less pronounced rough features of a sharper nature.These observations can possibly be explained by considering that cohesive particles tend to agglomerate into larger clusters, which appear to be more rounded at the scale of the full powder layer, compared to cohesionless particles which pack without demonstrating clustering, and thus it is more probable for them to have individual particles deposited on the surface of the layer, which macroscopically resemble sharp peaks.

Spatial fluctuation of the free surface height
While examining the digitized height distribution is informative, it does not contain information on the spatial arrangement of the height profile.For a given set of digitized height values, a permutation of their spatial arrangement does not change the distribution.This problem is analogous to the problem where a global packing density does not offer information on the homogeneity of the packing.Therefore, we again use the projection-based Voronoï and Delaunay tessellations as in Section 4 to address this issue.For a sphere packing, the digitized free surface height values for each sphere are spatially correlated as they can be fully described by the center coordinates and the radius of the sphere.To reduce this correlation, a spatial coarse graining at the length scale of a particle's diameter is required.We average the free surface height value, z s , in each Voronoï cell, defining a cell-averaged height, z v , at the center of each corresponding sphere, which is also a vertex in the Delaunay triangulation.The calculated distributions of z v for different cohesion are shown in Figure 8(a), which are colored by the corresponding Bo.Results show that the distribution widens with increasing Bo, which agrees with results in Figure 7.The variance of the z v distribution can be calculated as σ 2  zv , but as mentioned earlier, it does not contain information of the spatial height fluctuation.To demonstrate the permutation problem, a sketch is made in Figure 8(b) where the Delaunay triangles are drawn in black, and the height of each vertical stick from a vertex represents z v .For simplicity, we show idealized scenarios with only two height values, which can be organized into scenario A where the short surfaces (blue) and tall surfaces (red) are spatially segregated, and scenario B where they are mixed, with both cases having the same height distribution.The degree of "mixing" between taller and shorter surfaces needs to be quantified for a more complete description of the free surface profile.This can be described by how different the values are for the three vertices in an triangle, and this difference is small for most triangles in A and large for most triangles in B. To quantify this, we use the square of the first spatial derivative, |∇z v | 2 , and an integration of this quantity over the entire domain, which is essentially the Dirichlet Energy [51,52,53] which quantifies the degree of variation of a function in a given domain, with the function being the height, z v , that varies on the 2D domain A on the xy plane.In a lattice triangulation, this quantity can be digitized as where Λ k is the normalized Dirichlet energy for a single triangle k, and l, m, n are the vertices of k, and α lm is the angle facing the edge connected by l and m.The distribution of Λ k of all analyzed triangles for each Bo is shown in Figure 8(c).For each Bo, the distribution is a straight line on a log-lin scale suggesting an exponential distribution, P (Λ k ) = λe −λΛ k .Unlike the distribution of the digitized surface height with complicated shapes and spikes at z s = 0 (Figure 6), the exponential distribution can be conveniently described by a single parameter, λ, which sets the rate of decay for P (Λ k ).It can be seen from Figure 8(c) that the more cohesive cases have faster decays with higher values near zero.
To quantify the variation and the decay, we plot the total Dirichlet Energy for each Bo in Figure 8(d) and the fitted distribution parameter λ as an inset, both decreasing with Bo.Note that with the normalization by σ 2 zv , these two quantities truly reflect the blending of taller and shorter surfaces, not the spread of surface heights.These results quantitatively show that at the length scale set by particle size, higher cohesion results in less local height fluctuation, despite having a higher spread in height values.This is because low cohesion particles pack densely and homogeneously, and the surface height fluctuates at the particle scale, which is similar to scenario A in Figure 8(a).On the other hand, high cohesion particles form dilute and heterogeneous packings with clustering that is more similar to scenario B. In this sense, the local surface height fluctuations and the local packing anisotropy in a layer should be closely related, which is subject to future studies.

Conclusions
This work quantifies the structural features with a focus on density and surface roughness in powder layers in DEM simulations using realistic cohesive interaction forces.We first used a more traditional approach by calculating global values to show the general trend of decreasing density and increasing surface roughness as cohesion increases.The global structural features was calculated by digitization of the simulated spheres at a very fine scale and then samples globally by binning for density and by examining the distribution for the surface height profile.The increase in the surface roughness was then further interpreted by examining higher moments of the height distribution, including the skewness and the kurtosis, both show a gradual evolution for layers with neighboring Bond numbers.In particular, for Bo = 0 the skewness S sk < 0 and kurtosis S ku > 3, indicating that most rough features appear above the mean height plane and have sharp peaks, while for Bo = 272.3we observe the inverse trend, i.e. the skewness S sk > 0 and kurtosis S ku < 3, indicating that most rough features appear below the mean height plane and have more rounded peaks.
To highlight the increasing heterogeneity of the density and surface profile, we also developed Voronoï-based metrics that quantifies the spatial fluctuations of these quantities of interest.For density fluctuation, the divergence of the Voronoï anisotropy vector, Q k , was adopted for the thin deposited particle layers as a geometrical measure of their structural heterogeneity.The transition in the slope of Q k distribution at Q k = −1 displays a signature for the coexistence of dense regions with homogeneous structures as well as dilute regions with highly anisotropic structures, which is typical for cohesive materials.With increasing cohesion, both the standard deviation and skewness of the Q k distributions exhibit a consistent, monotonic change, indicating increasing structural heterogeneity of the deposited layer.
We quantified the fluctuation of the free surface height using the Voronoï cell-averaged height.Instead of focusing on the global distribution of this height, which contains no information on the spatial arrangement of the height values, we calculate the local squared spatial gradient as a measure of how well the taller and shorter surfaces are mixed.The distribution of the squared gradient is exponential which can be quantified by a single parameter.When normalized by the variance of the surface height, both the fitted distribution parameter and the total sum of the squared gradient show a decrease with increasing Bo.This quantitatively demonstrates that higher cohesion leads to reduced local height fluctuation despite the height values having a wider spread, which is possibly because that the packing density heterogeneity results in significant fluctuations of the free surface at a larger length scale.In contrast, at lower cohesion levels, particles densely and homogeneously pack, resulting in more surface height fluctuation at the particle scale.
The additional sets of metrics for the spatial fluctuation of density and surface height, combined with the global metrics, offer a more complete description of the packing structure than the traditionally used bulk-averaged values.This set of parameters can not only serve as a quantification of the quality of spreading but can also be used as a structural basis for modeling and analysis of subsequent processes, such as heat transfer and binder infiltration, as the heterogeneity of the packing structure on the particle level is important in these processes.For thicker layers and non-spherical particles, the 2D projection-based Voronoï calculation could lose its validity, but the same concept can be extended using real 3D set Voronoï analysis [54,26].

Figure 1 :
Figure 1: Numerical setup for powder spreading on a planar substrate during spreading showing (a) initial configuration and spreading setup for (b) cohesionless and (c) cohesive powders.

Figure 2 :
Figure 2: Local packing fraction ϕ as a function of Bond number Bo.The horizontal lines note the global packing fraction value for each layer.The sample points are colored according to their distance x from the starting point of spreading.

Figure 3 :
Figure 3: Characterizing the structural anisotropy of the deposited layer.(a) Granular packing of powder layer for Bo = 190.6.The points represent the projections of particle centers on the xy-plane.The corresponding 2D Voronoï tessellations are also shown on the plane.Also shown in (a), are a schematic representation of the particle packing with superimposed Voronoï tessellation (blue) and Delaunay triangles (green).The vectors Cp (red) point from particle centers to the centroids of Voronoï cells.The calculated Q k for (b) Bo = 0 and (c) Bo = 272.3,where the Delaunay triangles are colored by the corresponding Q k values.

Figure 4 :
Figure 4: (a) Probability density of normalized divergence of center-to-centroid vectors for the quasi-2D packing of powder deposited for different Bo.Q k > 0 regions are more densely packed than their surroundings; hence, we call these regions overpacked.Q k < 0 regions are more loosely packed than their surroundings and are, therefore, labeled underpacked.The solid curves are Gaussian fits Qk − 0.5 to Qk + 0.5.(b) Standard deviations (red circles) and skewness (blue squares) vs. Bo.The standard deviations and skewness of the distributions of Q k (solid) and Q ′ k

Figure 3 (
Figure 3(b) and (c) show the Q k map for Bo = 0 and Bo = 272.3,respectively, where the triangles are colored according to the corresponding values of Q k .For Bo = 0, the triangles share similar areas, and the Q k value fluctuates between positive and negative randomly in space.For Bo = 272.3,large triangles corresponding to underpacked regions exist, making the nearby Q k values highly positive or negative, indicating strong anisotropy.Note that the Q k value of each individual triangle is determined by its immediate neighborhood, rather than the overall packing density.The dense and homogeneous regions in both Bo = 0 and Bo = 272.3have low Q k values and only anisotropic regions show extreme values, which mostly exist in Bo = 272.3.

Figure 5 :
Figure 5: Height of powder layer surfaces for various Bond numbers from Bo = 0 to Bo = 272.3.

Figure 7 :
Figure 7: Surface roughness parameters (a) arithmetic mean height Sa (b) root mean square height Sq (c) skewness S sk and (e) kurtosis S ku shown as a function of Bo.The solid horizontal lines note the global values of the surface roughness parameters for each layer.The sample points are colored according to their distance x from the starting point of spreading.The dashed line for S sk = 0 marks the threshold between profiles where most rough features appear above the mean plane (S sk < 0) and below it (S sk > 0).The dashed line for S ku = 3 marks the threshold between rough profiles with rounded peaks (S ku < 3) and with sharp ones (S ku > 3).

Figure 8 :
Figure 8: Quantifying mixing of surface heights.(a) Distributions of Voronoï cell-averaged surface height for different Bo.(b) Illustration of poor mixing (left) and well mixing (right) that generate the same surface height distribution.(c) Distributions of the Dirichlet energy of individual triangles for different Bo.(d) The total Dirichlet energy for each Bo.Inset shows the fitted exponential distribution constant for each Bo.