Next Article in Journal
Digital Twin Modeling for Smart Injection Molding
Previous Article in Journal
Faster Evaluation of Dimensional Machine Performance in Additive Manufacturing by Using COMPAQT Parts
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Study on Powder Spreading Quality in Powder Bed Fusion Processes Using Discrete Element Method Simulation

by
Panagiotis Avrampos
* and
George-Christopher Vosniakos
*
Manufacturing Technology Laboratory, School of Mechanical Engineering, National Technical University of Athens, Heroon Polytehniou 9, 15773 Athens, Greece
*
Authors to whom correspondence should be addressed.
J. Manuf. Mater. Process. 2024, 8(3), 101; https://doi.org/10.3390/jmmp8030101
Submission received: 22 April 2024 / Revised: 10 May 2024 / Accepted: 11 May 2024 / Published: 16 May 2024

Abstract

:
Powder deposition is a very important aspect of PBF-based additive manufacturing processes. Discrete Element Method (DEM) is commonly utilized by researchers to examine the physically complex aspects of powder-spreading methods. This work focuses on vibration-assisted doctor blade powder recoating. The aim of this work is to use experiment-verified DEM simulations in combination with Taguchi Design of Experiments (DoE) to identify optimum spreading parameters based on robust layer quality criteria. The verification of the used powder model is performed via angle of repose and angle of avalanche simulation–experiment cross-checking. Then, four criteria, namely layer thickness deviation, surface coverage ratio, surface root-mean-square roughness and true packing density, are defined. It has been proven that the doctor blade’s translational speed plays the most important role in defining the quality of the deposited layer. The true packing density was found to be unaffected by the spreading parameters. The vertical vibration of the doctor blade recoater was found to have a beneficial effect on the quality of the deposited layer. Ultimately, a weighted mean quality criteria analysis is mapped out. Skewness and kurtosis were proven to function as effective indicators of layer quality, showing a linear relation to the weighted means of the defined quality criteria. The specific weights that optimize this linearity were identified.

1. Introduction

Selective Laser Sintering/Melting (SLS/SLM) is a Powder Bed Fusion-based (PBF) additive manufacturing (AM) process. During this process, a three-dimensional object is manufactured by sequentially sintering or melting powder layers via a focused laser beam at the predefined points. By the superposition of the successive layers, an object is formed [1,2].
Most commercial SLS/SLM machines feature a built-in powder deposition system (PDS) that is responsible for performing the task of spreading even, homogeneous layers of powder one above the other. The most widespread PDSs use a simple, compliant doctor blade which, as it translates parallel to the fabrication piston, evens out the top surface of the powder, creating a layer of the desired theoretical layer thickness. The recoater works in conjunction with a dual powder piston, one of which, the feeding piston, moves upwards to provide the necessary amount of powder, and the other one, the fabrication piston, moves downwards in order for the powder to be deposited on its surface, where it will then be sintered or melted. Alternate methods of powder deposition have been examined [2]. However, it has been found via Analytical Hierarchy Process (AHP) analysis that the most beneficial is the purely mechanical powder deposition, i.e., the one that utilizes a simple mechanical recoater, such as a doctor blade or a roller [2]. It has been proven via DEM simulations that the counter-rotating roller is more beneficial compared to the simple doctor blade since it promotes a higher packing density of the deposited layer [3]. Furthermore, researchers have used DEM simulations to examine the most important powder deposition aspects in terms of the translational speed of the recoater [4], blade shape [5,6], roller rotational speed [4], roller diameter [7], etc.
It is a fact that powder spreading greatly affects the finished part’s quality, as has been proven by many works, e.g., [8,9,10]. More specifically, Ziegelmeier et al. proved that low powder bed volume fraction is connected to the high porosity of the finished part, while low surface roughness of the un-sintered powder layer leads to poor surface quality of the finished part [11]. Despite that, not many works focus on the powder kinematics that govern the process. Some works utilize DEM simulations to reach conclusions about the powder spreading process; however, not many works systematically identify layer quality criteria, strictly define them and follow statistical analysis or DoE to reach safe, universal conclusions. Haeri et al. examined the performance of roller and doctor blade recoaters on the spreading of elongated particles with various aspect ratios based on the packing fraction and the surface roughness of the powder bed. However, they did not provide a robust statistical analysis to identify possible interactions of the spreading parameters. They also omitted to adequately define the volume fraction of the powder bed [3]. Tan et al. studied a composite powder made of elongated polymer and glass fiber. They calibrated the powder via angle of repose checks and only examined roller recoaters. However, they only examined layer thickness, spreading velocity and powder mixture composition without providing enough data on possible interactions between the examined parameters [12]. Wang et al. examined the possible connection between various recoater shapes to the surface coverage and packing density of the powder bed; however, they did not adequately examine recoating translational velocity. Instead, they focused more on the theoretical layer thickness, i.e., the gap between the bottom of the recoater and the fabrication surface [13]. In all these works, the common point is a lack of a robust statistical analysis that would provide a quantification of the qualitative results and observations provided by the DEM simulations.
To the authors’ knowledge, the present work is the only one that provides a systematic method of tackling the complex aspect of powder bed characterization and optimization via DEM simulations, by simultaneously focusing in five different aspects, namely the following: (i) thorough examination of minimizing computation time methods without hampering the validity of the simulations; (ii) calculation of particle size distribution as long as adhesive and cohesive forces and established contact model selection; (iii) accurate definition of six powder layer quality criteria (layer thickness deviation, surface coverage ratio, surface root mean square roughness, true packing density, skewness and kurtosis); (iv) Taguchi DoE with subsequent ANOVA for all the quality criteria to identify optimum parameter sets and develop regression equations of the criteria, as long as identifying possible interactions between the spreading parameters with possible ways to compensate for, e.g., higher spreading speed; (v) statistical analysis that proves skewness and kurtosis as equivalent powder bed quality criteria based on weighted means of the predefined four with least-squares analysis for best fit.
In this work, in order to evaluate mechanical powder spreading via a doctor blade, a Taguchi Design of Experiments (DoE) is implemented. The tests were carried out via Discrete Element Method simulations in the EDEM (version 2022.2) software package, which was kindly provided by AltairΤΜ. Firstly, the authors present the reasoning behind the selection of the particle size distribution and the particle contact model selection for computational time economy without damaging the reliability of the simulations. Secondly, the simulation geometry and bodies, as well as the simulation steps, are presented. Thirdly, the authors define the deposited layer quality criteria and present the ways in which they are calculated. Then, the Taguchi DoE is presented, with its selected L27 orthogonal array and the ANOVA and regression equations, as calculated by the statistical analysis. Afterward, the authors approach the quality criteria via their weighted means and connect them to surface skewness and kurtosis for a novel method of equivalent surface quality evaluation. Finally, conclusions are presented, along with suggestions for future work.

2. Materials and Methods

2.1. State-of-the-Art in Modelling of Powder Spreading

2.1.1. Critical Timestep and Simulation Time

In order for the DEM simulation to keep a low computational time while running in a stable manner, the critical time step Δtc is calculated using the Rayleigh equation [14],
Δ t c = π · r m i n β ρ G
where r m i n is the minimum particle radius in the simulation’s domain [15,16], ρ is the powder bulk density, G is the powder material’s shear modulus, and β can be approximated by Equation (2) [17,18], with ν being the powder material’s Poisson ratio.
β = 0.8766 + 0.163 ν
It is common practice to use lower G values for particle materials compared to the real-world value, granted that preliminary checks ensure small deviations in terms of results and powder behavior.
Zhang et al. (2020) selected a value of 3 GPa [7], outside the suggested range (1 MPa–1 GPa). However, this value is still two orders of magnitude below its actual value (150 GPa) [19,20,21]. After running shear cell tests to define the interparticle as well as the particle/wall coefficients of restitution, static and rolling friction, they ran cross-checking simulations/experiments of the angle of avalanche and angle of repose tests, finding deviations of +1.1% and +1.6%, respectively, validating the realistic behavior of the powder within the simulation. However, Chen et al. (2017) [22] examined how many orders of magnitude can decrease E without affecting the simulation’s results for an angle of avalanche test. If ν is constant, then E and G are linearly related: E = 2G∙(1 + 2ν).
Chen et al. proved that when E varied between the real-world value, E 0 , and the value, 0.001 E 0 , the simulation was in agreement with the real experiment in terms of avalanche angle, powder mixing mechanics and behavior. However, when E was set to any value below 0.0007 E 0 , the angle of avalanche was reduced by approximately 15% (from 33° to 28°), and, due to the reduction in shearing forces, the mixing rate was decreased [22]. However, drum powder mixing relies more on shear forces than doctor blade powder recoating. Fouda et al. (2019) studied doctor blade spreading of a mono-sized Ti6Al4V powder. They calibrated the simulation via cross-checking simulation and experimental results of an angle of repose test. They also studied the shear modulus range of 1 MPa to 100 MPa, i.e., 2.2 · 10 5 G 0 to 2.2 · 10 3 G 0 , with G 0 being the real-world value (42.5 GPa). Only minimal powder kinematics differences were observed; hence the smallest value was selected.

2.1.2. Particle Cohesion/Adhesion Models

“Cohesive” forces refer to those that are applied between particles of the same material, while “adhesive forces” refer to those that are applied between different materials, e.g., between the particles and the doctor blade, the deposition plate, or any other bodies of the simulation. The motion of very fine powders (D ≤ 100 μm) is not solely governed by gravity but also by attractive Van der Waals forces [23].
There are four different models that describe auto-adhesive particle behavior due to Van der Waals forces; the JKR (Johnson–Kendall–Roberts) model [24], the DMT (Derjaguin–Muller–Toporov) model [25], the BD (Bradley–Derjaguin) model [26] and the MD (Maugis–Dugdale) analytical model [27].

Hertz–Mindlin with JKR Model

The general relationship between normalized normal force and normalized relative particle approach in the JKR model is given by Equation (3). Figure 1 depicts the contact stages of two particles as described via the JKR model.
δ n δ t o = 3 F n J K R F p o J K R + 2 + 2 1 + F n J K R F p o J K R 1 2 3 2 / 3 F n J K R F p o J K R + 2 + 2 1 + F n J K R F p o J K R 1 2 1 / 3 ,   F n J K R F p o J K R 1
where F p o is the maximum tensile force required to break the contact among the two particles (i.e., the “pull-off force”) (see Equation (4)):
F p o J K R = 3 2 π R Γ
And δ t o is the relative approach at the tear-off moment (see Equation (5)):
δ t o = 3 F p o 2 16 R E 2 1 / 3
JKR model gives the Equations (6) and (7) for the normal force F n J K R and relative approach δ n , respectively, as functions of the contact patch radius, a .
F n J K R = 4 E a 3 3 R ( 8 π E Γ a 3 ) 1 / 2
δ n = a 2 R 2 π Γ a E 1 / 2

Tabor Parameter and Model Applicability

JKR and DMT initially appeared to be competitive with each other. However, they soon appeared to be limited to a range of solutions governed by the non-dimensional Tabor parameter μ (see Equation (8)) [23]. Grierson et al. (2005) [29] state that “Tabor’s parameter is physically equivalent to the ratio between the normal elastic deformation caused by adhesion (i.e., in the absence of applied load) and the spatial range of the adhesion forces themselves”.
μ = R Γ 2 Ε 2 z 0 3 1 / 3
where z 0 is the equilibrium separation in the Lennard–Jones potential, Γ is the work of adhesion (see Equation (9)):
Γ = γ 1 + γ 2 γ 12
where γ 1 and γ 2 are the surface energies of the two contacting bodies and γ 12 is the interface energy. For the same material contacting particles, γ 12 = 0 , and Equation (9) gives Γ = 2 γ , since γ 1 = γ 2 = γ .
Opinions on contact adhesion model applicability within the areas of the Tabor parameter vary. Thornton [23] suggests that, for μ < 0.1 , the DMT model is suitable, while for μ > 5 , the JKR theory is preferable. For the values that lay in between, Maugis (1992) developed the analytic Maugis–Dugdale (MD) approach [27]; however, its analytic nature renders it more complicated than DMT and JKR. Nevertheless, it has been proven that MD is accurate in all μ areas [30], while JKR and DMT are the MD approximations at high and low μ values, respectively. Greenwood (1997) states that, for μ 3 , JKR theory accurately approximates the actual radius of contact, the shapes and pressure distributions [26]. Greenwood also proves that DMT theory is “wrong in theory and in practice” (Greenwood, 2017) [26]. Instead, the Bradley–Derjaguin Equation (10) describes accurately how the tensile load varies with regard to the separation. Equation (10) is accurate for all values of μ for large positive separations and for negative separations for small values of μ [26].
T δ n = 2 π R Γ 4 3 z 0 δ n 2 1 3 z 0 δ n 8
MD model implementation in a DEM simulation would drastically increase the simulation time since it is analytic. In general, the JKR model is applied for large-diameter and more compliant materials (smaller E), while the DMT model is preferable for smaller bodies of more rigid materials (larger E) [30]. Table 1 summarizes the particle contact model range applicability.

2.2. Simulation Setup

2.2.1. The Simulation Software Platform

The simulation platform used was the EDEMTM, version 2022.2, provided by AltairTM. EDEMTM is a Discrete Element Method simulation platform that provides the capability of handling particles of various sizes, shapes and materials. EDEMTM provides the capability of importing physical bodies designed in 3D CAD environments, such as SolidWorksTM, and populating certain virtual geometrical shapes with particles of predefined physical properties and particle size distribution. EDEMTM offers a wide range of contact models from which the users can select the most suitable model combinations for the simulation they desire to run. These models include friction, electrostatic or Van der Waals adhesion and cohesion, air drag, restitution and others, covering a wide range of physical phenomena. Finally, EDEMTM offers the capability of NVIDIATM CUDA computing. Graphics card programming drastically decreases the simulation time for problems that are compatible with parallel computing.

2.2.2. Powder Modelling

The powder selected for the experiments was the SA-ZL-20 spherical alumina (Al2O3) powder provided by the manufacturing company Sinoenergy GroupΤΜ (Beijing, China). The powder specifications can be seen in Table 2.

Material Properties Modelling

The following values were used in the simulation: Poisson’s ratio ( ν ) of alumina equal to 0.3 [7], the solid’s bulk density ( ρ ) equal to 3820 kg/m3 [7], the shear modulus ( G ) equal to 1∙107 Pa. Furthermore, the alumina–alumina and alumina–SAE 304 stainless steel restitution (cr), static friction (csf) and rolling friction (crf) coefficients were used [7], and the values were as follows: c r , a l a l = 0.5 , c s f , a l a l = 0.34 , c r f , a l a l = 0.05 , c r , a l 304 = 0.52 , c s f , a l 304 = 0.2 and c r f , a l 304 = 0.05 .
The minimum particle radius is set at r m i n = 2.5   μ m , in order to not allow very small particles to enforce a very small timestep and to prevent the powder from displaying an extremely cohesive behavior. Similarly, the maximum particle diameter is set at D m a x = 75   μ m , 25% smaller than the theoretical layer thickness h l , t h (100 μm), in order to prevent quality deterioration of the surface [8]. Finally, as explained in Section 2.1.1, the shear modulus is set at a low value (10 MPa) after checking that this value minimizes the simulation time without affecting the validity of the results (see Section 3.2).

Particle Size Distribution

The powder selected was spherical alumina; hence, a uniform spherical shape was selected for all particles. Spreadsheet calculations were carried out in order to examine which kind of distribution best fits the three PSD numbers (D10, D50 and D90) provided by the powder manufacturer.
In many studies, a lognormal distribution is used to approximate the particle size distribution of a given powder; however, no method of defining the lognormal distribution parameters based on the aforementioned powder PSD specifications has been documented.
Let Z be a standard normal variable; then, the variable X = e μ + σ Z follows the lognormal distribution with parameters μ and   σ , μ , + , σ > 0 are real numbers. The mean and the standard deviation of the variable X and the cumulative distribution function are given by Equations (11), (12) and (13), respectively; e r f is the error function.
μ Χ = e 2 μ + σ 2
σ Χ = μ Χ 2 e σ 2 1
F Χ x = 1 2 1 + e r f l n x μ σ 2
The following equations characterize the particular powder employed:
D 90 = 47.5 D 50 = 21 D 10 = 8.2 F X 47.5 = 0.9 F X 21 = 0.5 F X 8.2 = 0.1   (14)
  (15)
  (16)
Equations (15) and (16) are used to calculate μ and σ , using Equation (13): μ = 3.0445 and σ = 0.7338. If Equation (14) is used for validation, it yields FX(47.5) = 0.867 ≈ 0.9. If Equation (13) is used to calculate D 90 , it gives x = 53.78 ≅ 47.5 with a deviation of 13.2%, which is deemed acceptable. Thus, using Equations (11) and (12), the values of μΧ and σX result as μΧ = 27.488 and σΧ = 23.216, and the values of their normalized counterparts (by the particle diameter Dave = 50 μm): 0.549758 and 0.464323. Furthermore, the user sets the maximum and minimum particle size caps, namely 0.1 and 0.5, as a number to be multiplied by the average particle size (“Physical Radius”).

Contact Model Selection

The model requires a value for “contact radius”, which has to be greater than the physical radius, and it is only used in the case where cohesive and adhesive attractive forces are to be taken into consideration.
There is a multitude of mechanisms by which microscopic particles may adhere to each other. In the case of relatively strong bonds, it may be solid, cemented or glued by a viscous liquid. Weaker bonds may be provided by pendular liquid bridges, Van der Waals forces, electrostatics or electro-magnetic fields [23]. In this work, however, there is no presence of viscous liquid or moisture that could create liquid bonds. Electrostatic charge is ignored since the minimum particle diameter is at 5 microns, which is not at the nanoscale, where the forces due to electrostatic charges begin to surpass gravitational forces. Even though some electrostatic charge can be developed, either via contact between the insulator (aluminum oxide) and the substrate (stainless steel plate) or via contact of particles of alumina with each other, granted they are of dissimilar size [31], they would be negligible due to the rather large average size of the powder sample used. Furthermore, even in the case of a fluidized bed of alumina particles, the charge of alumina dust is relatively neutral [32]. Hence, during powder deposition, where particle excitation is weaker, even more insignificant charges will develop. Similarly, no electro-magnetic fields exist, so only Van der Waals forces are to be considered. Furthermore, no plastic deformation was taken into consideration in this work since the materials are of very high hardness, and the low-speed collisions result in no permanent deformations.
For contact model selection, μ is calculated for a-alumina particles as follows:
The surface energy of amorphous alumina (a-Al2O3) is γ a A l 2 O 3 = 0.97 ± 0.04   J m 2 [33]. However, this value refers to nanoparticles of diameters between 2 and 5 nm and is valid for surface areas greater than 370 m2/g. The diagram of Figure 2 shows how the surface area of spherical alumina particles changes with the change in the particle’s diameter. The surface energy depends on particle size, the general consensus being that it decreases with the decrease in particle size [34]. This is justified since the number of next neighbors of surface atoms reduces with decreasing particle size [34]. Tepesch et al. (2000) calculated the surface energy of stoichiometric (1-Al-terminated) amorphous alumina (a-Al2O3) as γ a A l 2 O 3 = 2.13   J m 2 [35]. This surface is the most stable expression of a-alumina. Indeed, this value comes in accordance with the rule of thumb that as the particle size increases, the surface energy also increases. Hence, for two contacting particles of a-alumina, Γ 4.26   J m 2 .
R * and E * are calculated by the Equations (17) and (18)
R = R 1 R 2 R 1 + R 2
E = E 2 1 ν 2
Finally, the equilibrium separation in the Lennard–Jones potential for a-Al2O3 is z 0 a A l 2 O 3 = 5.1   Å = 5.1 · 10 10   m [36]. Figure 3 shows how μ varies with the diameter of the two contacting particles. G is adjusted to 1.29 · 10 11   P a , so that, combined with a v of 0.3, E is the maximum within the limit provided in [21], i.e., 413 GPa.
The Tabor parameter μ increases with the increase in the particle size (see Equation (8)). If the first among the two particles in contact has a diameter of 5 μm, then μ ranges between 1.492 and 1.840 for the second particle’s diameter of 5 and 75 μm, respectively. If the first particle’s diameter is 75 μm, however, μ ranges between 1.840 and 3.681 for the second particle’s diameter of 5 and 75 μm, respectively.
Judging by Figure 3, it can be deduced that, for the real-world material, almost all contacts of the larger particles (around 60–75 μm diameter) could be modeled via JKR or at least hybrid JKR-BD, as Greenwood suggests for 2 μ 3 values [26]. However, for smaller particles (around 5–15 μm diameter), the BD model (see Equation (10)) or even the analytic MD approach would be better. The MD approach would render the simulation much slower; however, the BD and hybrid JKR-BD models would be beneficial in areas where the JKR seems to deviate significantly from the analytical solution of the adhesive elastic contact problem.
In the trial runs that were performed in this study, in order to determine which shear modulus value should be used in the simulations, the shear modulus and, thus, Young’s modulus was varied from its actual value to a value 4 orders of magnitude smaller. The Tabor parameter’s value increases with the decrease in the relative modulus of elasticity ( E ) taking it to larger values by 3 orders of magnitude. Hence, the selection of the “Hertz–Mindlin with JKR v2” model in the EDEM software is perfectly justified.

Parametrization of Cohesive Forces

“Contact radius” (see Section 2.2.2, Sub-Section: “Contact Model Selection”) will, from this point on, be called R c , E D E M , and is given by the Equation (19).
R c ,   E D E M R p a r t i c l e + | δ t o |
where R p a r t i c l e refers to the “Physical Radius”, and δ t o is the relative approach of the two particles in contact (both of R p a r t i c l e radius) at which the contact breaks, with “to” standing for “tear-off” (see Equation (5)).
F p o is the maximum tensile force required to break the contact among the two particles (i.e., the “pull-off force”), given by Equations (4) and (20) for the JKR and Bradley (or DMT) model, respectively [23]
F p o B r a d l e y = 2 π R Γ
Finally, Equation (5), after replacing by Equations (4), (17) and (18), takes the form of Equation (21):
δ t o = 16.655 R 1 R 2 Γ 2 1 ν 2 2 R 1 + R 2 E 2 3
So, Table 3 is created by taking R 1 = R 2 = R a v e = 25   μ m .
It remains, to the authors’ knowledge, untested whether the contact radius setting has to be chosen based on the real-world G value or based on the G value that is set within the simulation environment in order to increase the timestep. Table 3 compares three different cases. By examining μ for each case, it is noticeable that the higher the particle’s stiffness, the smaller its relative approach at bond breakage. Indeed, when the material becomes less compliant and elastic, μ decreases by 3 orders of magnitude, from 1768.56 to 3.22, while the respective approach at bond breakage also decreases from around 1.5 μm to around 3 nm.
In order for adhesive/cohesive forces to facilitate, two particles or a particle and a body need to come into physical contact first. This check happens via the physical radii of the two bodies. However, the contact radius creates a zone around each particle where a negative approach is possible. As the two particles collide, they simultaneously demonstrate their compliance/elasticity and their cohesion. Upon collision, they behave like springs, initially with a positive approach, where they decelerate as they deform more and more. However, this collision is less “violent”, and they decelerate slower than they would if the cohesive forces had not been present since they are now attracted to each other due to Van der Waals forces. Eventually, the approach will stop increasing, and it will decrease again; however, when it becomes zero, the particles will not separate, but the “bottleneck effect” will develop since the particles will remain in contact due to the cohesive forces, and further work is required in order to create two new, separated surfaces. Figure 1 explains the stages of contact between two particles. The larger EDEM’s contact radius is, the larger the negative approach it will allow, hence the larger the “zone of influence” each particle will create around it, within which attracting adhesive forces will act. Hence, the larger the EDEM’s contact radius, the more cohesive the behavior the powder will exhibit. The real-world material is much stiffer and less cohesive compared to the one in the simulations. For the simulations, the value of the contact radius is 26 μm, which would be valid in the case where the relative approach at bond breakage would be 1 μm. By setting the difference between the physical and contact radii of the average particle at 1 μm, the simulation is easily scalable, while the increased cohesion makes this a study “on the safe side”, since, if the suggested system can counter the simulation material’s cohesion, the real-world material will be much easier to spread evenly.
It is of some interest to examine the fluctuation of the particles’ behavior in the simulation with regard to their real-world behavior (see Table 4).
Table 4 shows that the real-world material is stiffer and less compliant, which is also reflected in the calculated μ. The most suitable model for the real-world particles would be the BD, while the simulation implements a JKR pull-off force, which is 25% lower. Despite their higher rigidity, which is translated to shorter collisions of smaller deformations and sharper force profiles during contact with real-world particles, these demand a 25% larger pull-off force in order to escape from cohesive forces by their neighboring particles, which slightly compensates for the less cohesive behavior they demonstrate and is also visible by examining the relative approach at bond breakage.

2.2.3. Experimental Validation

In order to ensure that the modeled powder behaves realistically in terms of flowability and kinematics in any possible situation, a cross-check between simulation and experiment is necessary. It has been proven in the literature that a simulation–experimental cross-check for the combination of angle of avalanche and angle of repose tests is enough to ensure whether the powder model is adequate and whether its rheological characteristics approach the real powder [7,8]. This is explained since “the avalanche angle experiment suggests the powder to gravitational, cohesive and shear forces inflicted by the wall of the cylinder till a dynamic equilibrium state is reached (if possible), whereas the repose angle experiment only examines the powder’s cohesion in respect with the gravitational forces, letting the powder reach a static equilibrium state.” (Avrampos et al., (2022)) [8].
The powder material properties, particle–particle and particle–body interaction coefficients, and contact model implemented in this work were the same as in the work of Zhang et al. (2020) [7] (see Section 2.2.2. Zhang et al. (2020)); in their work, they initially ran a Jenike shear test (ASTM D6773 [37]) to define the static friction coefficients. The rolling friction coefficients were set to a low value due to the brittleness of alumina, while the restitution coefficients were referenced by the literature [38]. The interparticle and particle–wall interfacial surface energies were calibrated according to the cross-check for the angle of repose test.
They found the experimental values for the angle of repose and angle of avalanche to be 25.8° and 36.7°, respectively, while the simulation values were 26.2° and 37.1°, respectively, showing a deviation of only +1.6% and +1.1%, respectively; hence, proving the realistic behavior of the powder model that was used [7].

2.2.4. Simulation Geometry and Kinematics

The powder spreading simulation features 6 physical bodies and a virtual shape necessary to use as a powder factory. The 6 physical bodies are the substrate onto which the powder is deposited; the doctor blade, which performs the powder spreading; and 4 metal plates that serve as the left, right, back and front border of the powder layer. These geometrical borders create a square sample area with a side of 1 mm. The height of the back and front border plates is equal to the theoretical layer thickness, i.e., 100 μm, while that of the right and left border plates is equal to 1 mm (see Figure 4a).
Above the sample square, there is a virtual cube with a 1 mm side, which functions as the powder factory (see Figure 4b). The simulation populates the area inside this virtual cube with spherical powder particles of a diameter that follows the lognormal distribution that was described in previous sections. This area is simultaneously populated for all particles at the start of the simulation, ensuring that no particles are in physical contact with each other at the moment of their generation.
Then, the particles free-fall vertically onto the sample square under the influence of gravity, creating a slope in front of the doctor blade that performs the deposition (see Figure 5).
The smaller particles fall slower compared to larger ones due to the effect of air drag. The simulation enables natural gravity and Schiller–Naumann drag model [39,40]; the air’s characteristics are as follows: ρ a i r = 1.225 k g / m 3 ,     ν a i r = 1.81 · 10 5 P a · s ,   u a i r = 0 i ^ + 0 j ^ + 0 k , ^   s c a l e = 1 , where ν a i r is the kinematic viscosity of the air; u a i r is the velocity of the air, indicating that the simulation features no air flow; and “scale = 1” means that the amount of drag applied to each particle is the full amount with regard to its size and not a fraction of it.
The material of the physical geometry is stainless steel 304; its properties are as follows [41]: ρSt.304 = 8000 kg⁄m3, νSt.304 = 0.275, GSt.304 = 62.3 GPa, ΕSt.304 = 193 GPa.
The powder material properties and the powder–powder and powder–body interaction coefficients were set as described in Section 2.2.2, Sub-section “Material Properties Modelling”.
After the slope was created and the particles were all “frozen”, i.e., at a status of negligible kinetic energy, the doctor blade performed the spreading process. The spreading consists of a superposition of a linear translational motion with a constant speed along the axis of deposition (x-axis) and a vertical (along z-axis) vibrational motion, i.e., a sinusoidal oscillation of constant frequency and amplitude throughout the deposition process. In this example, the frequency was set to f v i b = 2000   H z and the amplitude to A v i b = 5   μ m . Similarly, the spreading speed was set to u t r = 0.01   m / s .

2.3. Surface Evaluation and Quality Criteria

The result of each powder-spreading simulation is a square layer of (1 mm × 1 mm), with a theoretical thickness of 0.1 mm. Then, a file containing the coordinates of the centers of the spherical particles, the particles’ diameter and the particles’ ID is exported. Particle ID is a unique number for each particle assigned to them at the moment of their generation. The same file contains the timestamp of the moment, which corresponds to the extracted data, and the total mass of the particles that comprise the layer. Table 5 shows the format of the exported data.
Figure 6 is a top view of the finished deposited powder layer. The particle coloring happens either via their centers’ Z-coordinate (see Figure 6a) or via their diameter (see Figure 6b). Particles that are positioned higher or have a larger diameter are assigned red colors, and the ones that are positioned lower or have a smaller diameter are assigned blue colors. The exported file (see Table 5) includes the data of every particle that constitutes the layer. To calculate the actual layer thickness, only the points that constitute the top of the exported layer are mathematically defined via the following methodology.
Initially, a grid is defined within the 1 mm × 1 mm square within which the powder is located. The grid has x-step Δx = 1 μm and y-step Δy = 1 μm. The grid size is 5 times smaller than the smallest particle’s diameter (5 μm). An even finer grid was proved to not increase the results’ accuracy; contrariwise, it increases the computational time and even causes an overflow in some cases. Hence, the square has 1001 x-values and 1001 y-values, the respective lines intersecting each other, developing a grid of n n o d e s = 1,002,001 “nodes”. The start of the coordinate system is positioned at the center of the square. From each node, a vertical linear “ray” is cast. This vertical line intersects some particles of the layer. To calculate the layer points, a script in C language, which uses the exported data (see Table 5) as input, compares the z-coordinates of all the intersection points of each “ray” with the particles and keeps the one with the highest z-coordinate. If a “ray” has no intersections with any particles, then the layer point that refers to this node is the node itself. The code is explained via the flowchart in Figure 7.
The point cloud of the layer’s top surface is then used to quantify the layer’s quality by four quality criteria, namely the layer thickness deviation (LTD), surface coverage ratio (SCR), root-mean-square surface roughness (Sq-RMS) and true packing density (PD or PDtr), which will be defined in Section 2.3.1, Section 2.3.2, Section 2.3.3 and Section 2.3.4.

2.3.1. Actual to Theoretical Layer Thickness Deviation

It is necessary to examine how closely the actual layer thickness approaches the theoretical layer thickness, which corresponds to the vertical distance between the substrate’s surface and the plane that the lower edge of the doctor blade creates as it moves horizontally along the x-axis. The vertical displacement of the fabrication piston is commonly assumed to be equal to the layer height, which is inaccurate. The actual layer thickness is slightly smaller due to interparticle cohesive and particle-recoater adhesive phenomena, where the impact varies with basic spreading parameters such as the recoater’s translational speed or the layer height.
Based on the calculated point cloud of the layer’s top surface, the actual layer thickness deviation is given by Equation (22):
L T D = h l h l , t h = 1 n n o d e s i = 1 n x j = 1 n y z L A Y E R i , j h l , t h
The actual layer thickness corresponds to the mean plane of the powder layer surface, i.e., the plane z = h l . In order to calculate surficial parameters (see Section 2.3.3), a “capital-z” ( Z ) height variable is defined, which corresponds to the plane z = h l , instead of the reference plane z = 0 ; see Equation (23) as follows:
Z i , j = z L A Y E R i , j h l

2.3.2. Surface Coverage Ratio

By examination of the layer’s top view, uncovered areas might be spotted, i.e., areas where the substrate is visible (see Figure 8). In optimum deposition, the gaps are filled with smaller-sized particles, completely eradicating such defects.
These uncovered areas might be caused by a larger particle/particle cluster being dragged along the surface or being pushed at such a high speed that its inertia enables it to keep moving for a longer than the desirable distance after losing contact with the recoater, causing a gap behind it. To create high-quality layers, it is of paramount importance to minimize the uncovered nodes, i.e., the nodes for which the c-code (see Figure 7) has returned no “ray”–particle intersections; hence, it has zero particles above it. The surface coverage ratio (SCR) is calculated by dividing the number of uncovered nodes by the total number of nodes within the sample square (Equation (24)).
S C R % = n u n c o v e r e d n n o d e s · 100 %
The maximum SCR value is 100% by definition, corresponding to perfect coverage. It was estimated by examining the 27 layer samples developed for the Taguchi DoE (see Table 6) that layers with SCR below 98.5% are considered to be unacceptable for PBF processes since the surface quality deteriorates drastically below that SCR value.

2.3.3. Root-Mean-Square Surface Roughness

The importance of the deposited layer’s surface roughness to the finished product’s quality has already been stressed in the literature [8]. Deep valleys or high peaks are detrimental to the quality since they can be inherited by the next layers, affecting the surface quality or even the part’s dimensional accuracy. It is common practice to extract powder layer profiles via methods such as white light scanning [42] to evaluate the two-dimensional surface roughness. This method’s disadvantage is the possible omission of profiles where significant defects are present, falsely leading to a higher layer quality estimation. To address this issue, the areal method is implemented in this work to calculate the surficial root-mean-square (RMS) roughness.
S q (see Equation (25)) is among the most widely used parameters. It corresponds to the standard deviation of the height distribution, generates reliable statistics and enables stable results since it is not significantly influenced by scratches, contamination, and measurement noise.
S q = 1 A A   Z 2 x , y d x d y
where A is the sample area (in this work A = 1   μ m 2 ), and Z x , y is given by Equation (23). As noted in Section 2.3.1, in order to calculate any surficial parameters, the height Z x , y that denotes the top surface of the layer refers to the mean plane of the surface, z = h l (see Equation (23)). Hence, from now on, the term Z refers to z L A Y E R h l . Since it is not possible to formulate an analytical equation to describe the top surface of the powder layer, Equation (25) must be solved numerically.
The integral in Equation (25) is the volume between the surface Z 2 x , y and the plane Z = 0 . The calculation happens via a column of points with coordinates P i , j x i , j ,   y i , j , Z i , j 2 . This volume is divided into many infinitesimal volumes shaped like truncated right rectangular prisms with square bases of Δ x = Δ y = 1 μm. More specifically, they are prisms with square bases; however, the four vertical edges of each prism are of various lengths (or heights) (see Figure 9). These volumes are calculated via the trapezoid rule.
The volume secluded between the nodes ( N i , j , N i + 1 , j , N i + 1 , j + 1 , N i , j + 1 ) and the points ( P i , j , P i + 1 , j , P i + 1 , j + 1 , P i , j + 1 ) is given by Equation (26):
Δ V Z 2 = Δ x · Δ y 4 · Z i + 1 , j + 1 2 + Z i , j + 1 2 + Z i + 1 , j 2 + Z i , j 2
Then, the integral V Z 2 = A   Z 2 x , y d x d y is calculated by summation of every infinitesimal prism.
If a node is located at one of the 4 corners of the sample square, the only infinitesimal prism it is part of is the one located at the corner of the sample square. However, if a node is located on one of the four edges of the square sample, it is part of two different infinitesimal prisms; hence, it contributes twice to the total volume. Finally, if a node is interior to the sample square, i.e., it does not belong to either an edge or a corner of the sample square, then it belongs to four different infinitesimal prisms, contributing to the total volume four times. Hence, S q is given by
S q = 1 A V Z 2
Similarly, the arithmetical mean height ( S a ), skewness (Ssk) and kurtosis (Sku) are calculated by:
S a = 1 A A   Z x , y d x d y
S s k = 1 S q 3 1 A A   Z 3 x , y d x d y
S k u = 1 S q 4 1 A A   Z 4 x , y d x d y
Skewness (Ssk-unitless) represents the degree of symmetry of the surface heights about the mean plane. Skewness’ sign indicates the peak or valley predominance of a surface. Negative skewness values indicate an upwards-deviated distribution, i.e., a surface with deep valleys and smoother, more even plateaus. Positive skewness indicates a downwards-deviated distribution, i.e., a surface with high peaks with smoother and less steep valleys.
Kurtosis (Sku-unitless) is used to evaluate sharpness in the height distribution. More specifically, Sku > 3 indicates the presence of inordinately high peaks/deep valleys or lack thereof (Sku < 3) making up the texture. For normally distributed surface heights (i.e., bell curve), Ssk is 0, and Sku is 3. Surfaces described as gradually varying, free of extreme peaks or valley features, tend to have Sku < 3 [43].
Skewness is useful for monitoring the effect of wear on a surface (e.g., evaluating the abrasion and oil sump of lubricants for slide planes), while kurtosis is used for spotting the presence of either peak or valley defects that may occur on a surface [44,45].
Other useful surficial parameters, like the peak-to-valley roughness (Sz), the maximum peak height (Sp) and the maximum peak depth (Sv), mostly identify the magnitude of the largest defects of the surface and fail to examine the consistency and evenness of the surface, which the aforementioned metrics examine more accurately.

2.3.4. True Packing Density

The packing density, or compaction ratio, is directly connected to the finished parts’ density and microhardness [46].
In many works, “theoretical” packing density is calculated by dividing the mass of a layer by the theoretical volume of the layer, which is calculated by multiplying its border dimensions along the three axes, x, y and z [47,48]. However, the actual layer height differs from the theoretical layer height (see Section 2.3.1). Theoretical packing density is an indicator of how well the powder layer is “filled” but fails to show how well the particles themselves have been arranged in order to increase their coordination numbers and minimize the air gaps between them. In this work, the true packing density is calculated by using the actual instead of the theoretical layer thickness (see Equations (31) and (32)).
P D k g m 3 = m L A Y E R V L A Y E R = 1 V Z i = 1 n p a r t i c l e s m p a r t i c l e , i
P D % = P D k g m 3 ρ p o w d e r ,   B U L K · 100 %

3. Results

3.1. Design of Experiments

The Taguchi DoE is a statistical tool used to reach conclusions when it comes to the effect of various factors on the result of an experiment with regard to some quality criterion that is set beforehand. The number of factors and the number of each factor’s levels leads to the selection of an orthogonal array (OA) suitable to accommodate all the necessary combinations of the factors’ levels in order to reach solid outcomes [49].
Four variables are examined:
  • Translational speed of the doctor blade along the deposition ( x ) axis ( u t r ), examined at three levels: 0.01 m/s, 0.05 m/s, 0.1 m/s. These levels were selected according to the current industrial SLS/SLM machines’ powder deposition speed standards. These were validated in both the prototype SLS/SLM powder deposition system designed [1,42] as well as via the industrial SLM machine of the Laboratory of Manufacturing Engineering of the School of Mechanical Engineering of NTUA (Z-rapid iSLM280, ZRapid, Suzhou, China).
  • Vibrational frequency of the doctor blade along the vertical ( z ) axis ( f v i b ), examined at three levels: 500 Hz, 1000 Hz, 2000 Hz. The frequencies selected were relatively low compared to, say, ultrasound frequencies to maintain a reasonably high simulation timestep. Furthermore, the simulations seek to examine whether a low-frequency-vibrating doctor blade benefits powder spreading compared to vibration-less deposition, which is the common powder deposition method encountered in the industry. Preliminary testing of vibration-less doctor blade powder deposition examined the deposited layer’s quality for a flat doctor blade at two different translational speed levels, i.e., 0.08 m/s, typical of commercial machines (see Figure 10a) and 0.01 m/s (see Figure 10b). The cross-check between vibration-less and oscillating recoater deposition proved that vibration enhances the homogeneity and evenness of the deposited layer.
  • Vibrational amplitude of the doctor blade along the vertical ( z ) axis ( A v i b ), examined at three levels: 1 μm, 2.5 μm, 5 μm. The amplitude values were selected in such a way that they are equal or less to the minimum particle diameter (5 μm), thus making it impossible for powder particles to escape behind the doctor blade via the gap between the doctor blade and the back border plate.
  • Angle of relief of the doctor blade ( θ r e l ), examined at three levels: 0°, 5°, 10°. The geometry of the doctor blade is a quasi-cubic 3D shape with a 1 mm side (see Figure 11). The feature differentiating it from a cube is its angle of relief, beginning after a 100 μm horizontal flat area, implemented to enhance the vibration’s powder compaction effect. Had this horizontal area not existed in the 5- and 10-degree doctor blades, negligible powder compaction would have been achieved via the blade’s edge. Furthermore, this would render the simulation non-realistic since no sharpened blades are used for powder deposition. The size of this feature (100 μm) was selected to be 25 μm larger than the largest particle diameter (75 μm).
The blade design resembles the relatively thick, commercially available doctor blades. Figure 11 shows that, as the angle increases, the secondary contact surface between the blade and the particles decreases. The primary contact surface is the left vertical doctor blade’s surface, which forces the powder slope to move toward the direction of u t r . The secondary contact surface is the area that applies vertical pressure via vibration onto the powder layer, i.e., the 100 μm wide strip and the fraction of the inclined surface behind this strip that interacts with the particles via the vibrational and adhesive forces. When the relief angle equals 0°, the secondary contact area is equal to the entire 1 mm2 of the doctor blade bottom area, while when the relief angle increases, the secondary contact area decreases to 0.1 mm2 plus the decreasing inclined contact zone, as described above.
The goal is to optimize the response of all quality criteria simultaneously. This is not a trivial feat since each response is individually optimized at different levels. This work attempts to find global optimum levels of the spreading variables for all of the quality criteria to achieve the best possible outcome. Initially, each of the four responses is examined separately.
Table 6 presents the setup of the powder spreading simulations in an L27 Taguchi orthogonal array. The L27 provides the capability of examining three interactions, which is enough for this case study. Examining more interactions would overcomplicate the design and lead to a number of trials comparable to the full factorial design. The quality criteria are described in Section 2.3.1, Section 2.3.2, Section 2.3.3 and Section 2.3.4, namely LTD, SCR, Sq-RMS and PD. The quality loss function is the lower, the better for LTD and Sq-RMS, and the higher, the better for SCR and PD. Sa, Ssk and Sku are also presented since they provide useful observations. The regression equations are calculated using the MinitabTM (version 17) software. All regressions used four continuous predictors, i.e., the four variables of the Taguchi DoE; according to these, four linear terms and three interaction terms were used, as shown in Table A1, Table A2, Table A3, Table A4 and Table A5 of Appendix A.

3.2. Shear Modulus Sensitivity

EDEMTM tutorials suggest an initial value of G around 1∙107 Pa since the actual value drastically decreases the timestep (see Section 2.1.1). The suggested G value ranges between 1∙106 Pa and 1∙109 Pa.
In Table 7, both the critical timestep and the total simulation time of the three repetitions of Trial #9 of the Taguchi DoE are presented, along with the surface quality criteria of the deposited layer.
Three orders of magnitude G reduction caused the simulation duration to drop tenfold, from 168 h to 15 h, while only |LTD| and Sq were minimally affected (increased by 4.4% and 6%, respectively). By a further G decrease by another order of magnitude, again only |LTD| and Sq increased by a further 14.4% and 4%, respectively; the percentages are calculated as a fraction of their real-world G simulation (repetition #3) values.
To keep the duration of the 27 Taguchi trials within acceptable limits, the smallest G value was selected. The results deviate substantially only for |LTD| (18.8%), but not at such a level to justify running simulations that last 12 ( 25 / 2.5 ) or 67 ( 168 / 2.5 ) times more. The study remains “on the safe side”, because, by decreasing G, the quality of the layer seems to deteriorate. If the surface quality is considered acceptable by the 10 MPa G value, then the layer produced by the real-world G value will be of higher quality, also rendering it acceptable.

3.3. ANOVA Results and Discussion

Figure 12 and Figure 13 present the results corresponding to means and signal-to-noise-ratio, respectively. The former depicts the average response for each combination of control factor levels, while the latter is the ratio of the mean (signal) to the standard deviation (noise). Hence, the means plot shows the average performance, while the SNR plot shows the variance from the design. Additionally, Table A1, Table A2 and Table A3 depict the ANOVA results for the LTD, SCR and RMS, respectively.
Figure 12 and Figure 13 deduce that none of the four factors has a statistically important impact on the true packing density, which ranges from 66.8% to 68.4%, while the mean PD value varies between 67% and 67.7%. Indeed, the powder’s compressing mechanism relies mostly on gravity, since even the vibrating blade only momentarily pushes the layer downwards and then moves upwards again.
The blade vibrates between zlow = hl,th and zhigh = (hl,th + Avib); hence, it does not compress the layer below the actual layer thickness level. Even though the theoretical packing density varies in an inversely proportional manner compared to the layer thickness deviation, the true packing density will not be examined as a quality criterion since it was proven that the four factors do not substantially affect it. The fact that the true PD remains unchanged might be caused by the small particle size, since the cohesive forces are large enough to prevent the gravitational forces from increasing the coordination number of each individual particle during the recoater’s compressive motion.

3.3.1. Layer Thickness Deviation

The “Contribution” column of the ANOVA for the |LTD| (Table A1) shows that |LTD| is primarily affected by utr (87.49%) and much less by Avib (5.70%). fvib and θrel contribute negligibly, at 0.36% and 0.12%, respectively, similar to the interactions between the factors. The only interactions examined are the interactions between utr and the other factors, since it is obviously the most impactful factor. Clearly, the smaller utr and the larger Avib are, the more the actual layer thickness approaches the theoretical layer thickness.
Figure 14 shows that utr-Avib and utr-fvib interaction plots have approximately parallel lines; hence, they demonstrate a low degree of interaction, only traceable at higher speeds. Both couples show no interaction at low- and middle-speed levels; however, at higher speeds, it becomes more beneficial to maintain a low frequency and high amplitude. The interaction between utr and θrel reveals that, at low and high speeds, sharper blades with large angles of relief perform better, while, at middle speeds, the optimum performance belongs to the flat-bottomed blade.

3.3.2. Surface Coverage Ratio

Figure 12 and Figure 13 show that SCR, like |LTD|, is primarily influenced by the utr (78.44% contribution, Table A2). Again, an increase in utr causes a decrease in SCR, deteriorating the layer quality.
Table A2 proves that θrel, fvib and Avib play a small role (1.60%, 1.02% and 0.19%, respectively) in terms of SCR. Figure 15 proves that the pair utrrel has negligible interaction; the flat-bottomed blades show better results at all speeds. However, for both the pairs utr-fvib and utr-Avib, there only is no interaction at low and middle speeds, and SCR is positively affected by higher frequency amplitude. Contrariwise, at greater utr values, this trend gets reversed in terms of Avib; smaller amplitudes give better SCR at maximum speed. In terms of fvib, middle-frequency levels drastically deteriorate SCR at maximum utr, leading to the conclusion that the natural frequency of the system is located at approximately 1000 Hz.
The regression’s error is higher than 10%, probably because of the high inherent complexity of the process. However, this does not harm the credibility of the results, as shown in the calculated deviation between the simulation and the expected regression value (see Table 8).

3.3.3. RMS Surface Roughness

Figure 12 and Figure 13 for Sq-RMS show that the smaller utr is, the smaller Sq is; hence, the better the surface quality. Sq is affected more by θrel compared to previous criteria. Indeed, flatter blades tend to create smoother surfaces due to the increased contact zone between them and the surface. Again, high frequency and amplitude levels show better results in terms of Sq.
Table A3 shows that utr is again the prevalent, affecting factor (65.79%); however, θrel and the fvib are also important (5.61% and 2.89% contribution, respectively). Despite the weak influence of Avib (0.32%), its interaction with utr is quite significant (6.37%).
Interestingly, as seen in Figure 16c, even though at low utr, the effect of different θrel is negligible on Sq; as the speed increases, the worse the results that are achieved via sharper blades. Furthermore, while between low and medium speeds, there seems to be no utr-fvib or utr-Avib interaction, this changes between medium and high-speed levels. While at low and medium speeds, it is beneficial to apply high Avib, at high utr, it is beneficial to opt for a lower Avib. This is possibly due to the high amount of kinetic energy fed towards the particles, which surpasses a threshold, leading to an unstable powder bed behavior increasing the surficial defects.
The error surpasses 10% due to the process complexity; however, this does not harm the credibility of the experiment and the validity of the results and conclusions (see Table 8).
In conclusion, if it is absolutely necessary to vastly increase the spreading speed (e.g., for production reasons), the negative outcome in terms of Sq-RMS and SCR can be partially compensated by opting for flat-bottomed blades and reducing Avib while maintaining high fvib. This leads to a slight |LTD| increase, which is easily countered via adjusting the downwards motion of the fabrication piston, i.e., increasing hl,th, while poor SCR and Sq-RMS signify important surficial defects, possibly inheritable to the next layers.

4. Discussion

Process Parameter Optimization

The optimum parameter sets for all three quality criteria is by ¾ identical (referring to utr, fvib and Avib) at u t r = 0.01   m / s , f v i b = 2000   H z , A v i b = 5   μ m with the observation that flat-bottomed blades promote optimum SCR and Sq-RMS, while sharp blades promote higher |LTD|, compromising the other two (see Table 9).
Table 6 reveals that all trials have negative Ssk, as expected, since the top hemispheres of the particles create surfaces with smoother, rounder plateaus and deeper, sharper valleys, with an upper-side deviated distribution. Furthermore, almost every trial features Sku > 3 (apart from #23 and #24), which is also expected, indicating the presence of inordinately deep valleys since the valleys reach the substrate at uncovered nodes.
After arranging the lines of Table 6 by ascending (i) Ssk and (ii) Sku order, respectively, the three quality criteria are plotted versus Ssk and Sku, respectively, in Figure 17 and Figure 18. Each set of points is approximated linearly via the least-squares method, and the coefficients of determination (R2) can be seen on the plot.
The skewness plot reveals that, as the layer’s quality deteriorates (an increase in |LTD| and Sq-RMS; decrease in SCR), Ssk increases, approaching zero from the negative side. Similarly, the kurtosis plot reveals that the lower the layer quality, the smaller Sku becomes, approaching zero from the positive side. Figure 17 and Figure 18 provide the least-squares linear approximation of the three quality criteria versus Ssk and Sku. Therefore, Ssk and Sku serve as equivalent indicators of surface quality.
A weighted means (WM) of |LTD|, (100-SCR) and Sq-RMS is used for optimization after they are normalized in the range 0–1 (see Equations (33) and (34)).
W M S s k = w L T D n S s k · L T D n + w 100 S C R n S s k · 100 S C R n + w S q n S s k · S q n
W M S k u = w L T D n S k u · L T D n + w 100 S C R n S k u · 100 S C R n + w S q n S k u · S q n
These WMs and the experiment are connected via plotting the WMs vs. the respective Ssk and Sku values because of the observation that the smaller (negative) Ssk and the larger (positive) Sku values become, the higher the layer quality.
A simple algorithm combines the three weights, given that their sum must always be equal to 1, aiming to develop a line via the least-squares method, which gives the best fit (highest coefficient of determination (R2) value) among all the possible lines created by various weight combinations. The algorithm has a 0.05 weight-increment step. The weight combinations that lead to the best-fit lines are the following:
w L T D n S s k = 0.35 w 100 S C R n S s k = 0.05 w S q n S s k = 0.60   W M S s k = 0.6449 · S s k + 0.9549   R 2 = 0.8869
w L T D n S k u = 0.40 w 100 S C R n S k u = 0.20 w S q n S k u = 0.40     W M S k u = 0.1497 · S k u + 1.1224 R 2 = 0.9499
The line depicting WMSsk versus Ssk is an increasing function of Ssk (see Figure 19). Since the higher the (negative) Ssk value, the higher the layer quality, it is concluded that the smaller the quantity 0.35 · L T D n + 0.05 · 100 S C R n + 0.60 · S q n , the better the layer quality, according to Ssk.
Similarly, the line depicting WMSku versus Sku is a decreasing function of Sku (see Figure 20). Since the higher the (positive) Sku value, the higher the layer quality, it is concluded that the smaller the quantity 0.40 · L T D n + 0.20 · 100 S C R n + 0.40 · S q n , the better the layer quality, according to Sku.
Both Ssk and Sku lead to similar observations with regard to WMSsk and WMSku that need to be minimized in order to optimize the layer’s quality. The weights calculated are similar to the qualitative importance that an expert in AM would assign to them. SCR is the least important since it detects the extent to which large surface defects exist within the layer. Such defects can be caused by dragging an agglomerate along the surface, creating a groove that leaves the substrate visible, drastically decreasing the SCR. However, such defects are easily avoidable by implementing proper sieving of the powder prior to recoating. Despite the possibility of such defects being catastrophic for the surface quality, since they are usually inherited from layer to layer and cause serious defects in the geometrical accuracy or even the feature geometry of the finished part, the easiness in their avoidance renders them quite uncommon. That is why it is assigned the lowest weight, at 0.05 and 0.2, for Ssk and Sku evaluation, respectively.
Secondly, |LTD| shows how lower the mean line of the layer’s surface is compared to the plane of the theoretical layer thickness. It is important since it is directly connected to the layer’s PDth. Despite being important, it can be countered by appropriately adjusting (increasing) the hl,th by increasing the vertical downwards displacement of the fabrication piston so as to achieve the desired hl. Hence, it is assigned a medium weight, at 0.35–0.4 for Ssk and Sku evaluation, respectively.
Sq-RMS quantifies layer evenness and homogeneity, showing whether the finished surface will have a high or low roughness value. High Sq-RMS values cannot lead to serious inherited defects throughout the layers; however, they can affect the surface quality of the finished part, making it necessary to be post-processed in order for it to become ready for use. Furthermore, Sq-RMS is the most consistent quality criterion, as it will identify, without fail, large areas where the layer’s surface is below or above the mean line, identifying intense spreading fluctuations during the powder–recoater interaction. Hence, it is assigned the highest weight, at 0.6 and 0.4, for Ssk and Sku evaluation, respectively.
In conclusion, the empirical approach of the method agrees with the calculated analytical results after the statistical analysis. Ssk and Sku can be used solely in conjunction with the weighted means of the three quality criteria in order to examine the general quality of the layer. Examining Ssk or Sku is equivalent to examining a weighted mean of the three normalized quality criteria, |LTD|n, (100-SCR)n and (Sq-RMS)n. The weights are in agreement with the ones that would have been chosen by experts in the field. Table 9 is completed via the regression–validation simulations. Both Ssk and Sku indicate that the optimum parameter level combination for the SCR and RMS values leads to a layer quality higher than that corresponding to the parameter values optimizing |LTD|. Hence, the optimum parameter level combination is: u t r = 0.01   m / s , f v i b = 2000   H z , A v i b = 5   μ m , θ r e l = 0 ° .
The vibration’s positive effect on the spreading outcome is proven by comparing the optimum level combination with a vibration-less trial where optimum values for utr and θrel are chosen. Metrics/visual evaluation are shown in Table 10 and Figure 21, respectively.
The optimum level combination determined by Taguchi analysis is superior to its equivalent vibration-less trial by every quality criterion, including Ssk and Sku. A top-view examination of the layers shows multiple surface defects in the vibration-less layer, as the SCR indicates. These are located behind larger particles, revealing the detrimental effect of large particle dragging. Vibration alleviates large particle and agglomerate dragging and partially breaks down agglomerates, decreasing the defective areas.
In the vibration-less trial, the height of the particles displays higher randomness, with many particle centers being over hl,th, as shown in the legend. Contrariwise, the particles in the optimum vibration-assisted trial have a uniform height, with the color showing much less variation, while the top of most particles is just contacting z = hl,th plane (see Figure 21b). The positive effect of vibration, combined with the fact that the layer’s PDth increases with applied vibration, solely via the decrease in the |LTD| and not via increased particle packing, which is proven to be negligible, come in agreement with the findings in [50], who proved these via DEM simulations for both a vibrating roller and a vibrating doctor blade recoater.
Since Ssk and Sku serve as general layer quality indicators, ANOVA for these is provided in Table A4 and Table A5, respectively, while Equations (35) and (36) provide the regression Equations.
S s k = 0.587 0.90 u t r + 0.000295 f v i b + 0.1455 A v i b + 0.0112 θ r e l 0.00318 u t r f v i b 1.222 u t r A v i b 0.00901 A v i b θ r e l
S k u = 3.498 9.95 u t r + 0.000656 f v i b + 0.4665 A v i b + 0.0582 θ r e l 0.00613 u t r f v i b 3.07 u t r A v i b 0.0273 A v i b θ r e l
The regression equations above were calculated in order to minimize the error of the ANOVA. More specifically, the interactions were selected to maximize the sum of the contribution column of the ANOVA table, hence minimizing the error of the regression equation.
By using Equations (35) and (36), it is confirmed that the optimum parameter combination is the one shown in Table 10. The regression optimum values are as follows: S s k o p t , r e g = 1.771 , S k u o p t , r e g = 6.766 .
The simulation-calculated optimum values are as follows: S s k o p t , s i m = 1.543   S k u o p t , s i m = 6.530 .
The deviation between the simulation and the regression-calculated optimum values is −12.9% and −3.5% for skewness and kurtosis, respectively. Kurtosis’ value is within the 8.77%, while skewness slightly exceeds the 9.26% error predicted by ANOVA. However, this is not considered to be a very important error, even though it is indicative of the fact that kurtosis might be a more reliable indicator of layer quality, which was also revealed by the R2 values of the best-fit lines.

5. Conclusions

This work utilizes DEM simulations in conjunction with robust statistical analysis in the form of Taguchi DoE with subsequent ANOVA to evaluate the effect of spreading parameters on the deposited layer quality for PBF processes. Pertinent individual contributions are itemized as follows.
  • A method is presented for connecting D10, D50 and D90 of a powder to specific lognormal distribution parameters to accurately define the particle size distribution.
  • Different ways of maximizing the timestep were exploited to minimize the computational time without affecting the validity of the results.
  • The powder cohesion is parametrized to achieve realistic behavior of the powder.
  • An appropriate contact model is selected based on powder characteristics.
The quality of the spread layer is proven to deteriorate with the spreading speed increase. Furthermore, oscillation of the doctor blade recoater promotes even flattening of the deposited layer and lower RMS roughness values, along with higher surface coverage ratios.
Different weighted means of the three quality criteria, namely LTD, SCR and RMS, have been proven to be connected to the layer’s surface areal skewness and kurtosis. Furthermore, skewness and kurtosis were proven to serve as general quality criteria for a deposited surface. The lower the (negative) skewness and the higher the (positive) kurtosis values, the higher the quality of the layer.
As a future step, Taguchi experiments using a prototype and a commercial powder deposition system are aimed at evaluating the layers via 3D white light scanning in order to cross-check the simulation results with actual experimental results on a larger scale.

Author Contributions

Conceptualization, Methodology, Formal Analysis, Investigation, Writing—Original Draft Preparation: P.A.; Writing—Review and Editing, Supervision, Project Administration: G.-C.V. All authors have read and agreed to the published version of the manuscript.

Funding

P. Avrampos gratefully acknowledges financial support from a 4-year PhD Scholarship from the National Technical University of Athens (PA2020-2024).

Data Availability Statement

The datasets presented in this article are not readily available, firstly because the data are part of an ongoing study and secondly because of the large file size. Requests to access the datasets should be directed to Prof. G.-C. Vosniakos.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

Table A1. ANOVA for |LTD|.
Table A1. ANOVA for |LTD|.
Source DoF SeqSS Contribution Adj SS Adj MS F-Value p- Value
Regression71182.5994.68%1182.59168.94148.330.000
Utr (m/s)11092.8287.49%65.7865.78318.820.000
fvib(Hz)14.510.36%14.2814.2844.090.058
Avib (μm)171.255.70%37.6737.66810.780.004
θrel (deg)11.560.12%0.650.6530.190.670
Utr (m/s)*fvib (Hz)19.770.78%9.779.7722.800.111
Utr (m/s)*Avib (μm)12.660.21%2.662.6590.760.394
Utr (m/s)*θrel(deg)10.010.00%0.010.0140.000.950
Error1966.425.32%66.423.496
Total261249.01100.00%
Table A2. ANOVA for SCR.
Table A2. ANOVA for SCR.
Source DoF Seq SS Contribution Adj SS Adj MS F-Valuep- Value
Regression7264.64484.48%264.644378.06314.770.000
Utr (m/s)1245.72078.44%13.689136.8855.350.032
fvib(Hz)10.32011.02%0.00510.005120.020.889
Avib (μm)10.05810.19%0.62100.620982.430.136
θrel (deg)10.50001.60%0.00210.002130.010.928
Utr (m/s)*fvib (Hz)10.09220.29%0.09220.092180.360.556
Utr (m/s)*Avib (μm)10.62592.00%0.62590.625912.450.134
Utr (m/s)*θrel(deg)10.29610.95%0.29610.296111.160.296
Error1948.63015.52%48.6300.25595
Total26313.274100.00%
Table A3. ANOVA for Sq-RMS.
Table A3. ANOVA for Sq-RMS.
Source DoF Seq SS Contribution Adj SS Adj MS F-Value p- Value
Regression7166.03284.33%166.032237.18814.600.000
Utr (m/s)1129.53065.79%0.5140.51360.320.580
fvib(Hz)15.6962.89%1.85318.5291.140.299
Avib (μm)10.6320.32%11.329113.2876.970.016
θrel (deg)111.0455.61%0.0500.05040.030.862
Utr (m/s)*fvib (Hz)10.0000.00%0.0000.00000.000.996
Utr (m/s)*Avib (μm)112.5396.37%12.539125.3907.720.012
Utr (m/s)*θrel(deg)16.5903.35%6.59065.9024.060.058
Error1930.86015.67%30.86016.242
Total26196.892100.00%
Table A4. ANOVA for |Ssk| in powder spreading.
Table A4. ANOVA for |Ssk| in powder spreading.
Source DoF Seq SS Contribution Adj SS Adj MS F-Value p- Value
Regression7314.07790.74%314.0770.44868126.610.000
Utr (m/s)1238.47268.90%0.004010.0040080.240.631
fvib(Hz)10.165784.79%0.294220.29421717.450.001
Avib (μm)10.091492.64%0.338030.33803320.050.000
θrel (deg)10.091882.65%0.014370.0143710.850.367
Utr (m/s)*fvib (Hz)10.159004.59%0.143020.1430168.480.009
Utr (m/s)*Avib (μm)10.148794.30%0.148790.1487868.820.008
Utr (m/s)*θrel(deg)10.099112.86%0.099110.0991105.880.025
Error190.320409.26%0.320400.016863
Total26346.117100.00%
Table A5. ANOVA for Sku in powder spreading.
Table A5. ANOVA for Sku in powder spreading.
Source DoF Seq SS Contribution Adj SS Adj MS F-Value p- Value
Regression7301.81791.23%301.81743.11728.250.000
Utr (m/s)1243.77473.69%0.48600.48603.180.090
fvib(Hz)111.3293.42%14.52414.5249.520.006
Avib (μm)120.3686.16%34.73734.73722.760.000
θrel (deg)10.16240.49%0.38520.38522.520.129
Utr (m/s)*fvib (Hz)10.62661.89%0.53320.53323.490.077
Utr (m/s)*Avib (μm)10.93912.84%0.93910.93916.150.023
Avib (μm)*θrel(deg)10.90642.74%0.90640.90645.940.025
Error1928.9998.77%28.9990.1526
Total26330.816100.00%

Appendix B

Table A6. Table of variables and abbreviations.
Table A6. Table of variables and abbreviations.
Context Description/Definition Symbol/ Abbreviation Unit Equations Tables
GeneralAdditive manufacturingAM---
Selective Laser Sintering/Selective Laser MeltingSLS/SLM---
Powder Bed FusionPBF---
Powder deposition systemPDS---
Analytic Hierarchy ProcessAHP---
Discrete Element Method/ModellingDEM---
Design of ExperimentsDoE---
Analysis of VarianceANOVA--(Table A1, Table A2, Table A3, Table A4 and Table A5)
Particle size distributionPSD--(Table 2)
Hertz–Mindlin with JKR contact modelJKR--(Table 1), (Table 4)
Derjaguin–Muller–Toporov contact modelDMT--(Table 1)
Bradley–Derjaguin contact modelBD-(10)(Table 1), (Table 4)
Maugis–Dugdale analytical contact modelMD--(Table 1)
Layer thickness deviationLTDμm(22)(Table 6, Table 7, Table 8, Table 9 and Table 10)
Surface coverage ratioSCR(%)(24)(Table 6, Table 7, Table 8, Table 9 and Table 10)
Root-mean-square surficial areal roughnessSq-RMSμm(25), (27), (29), (30)(Table 6, Table 7, Table 8, Table 9 and Table 10)
Arithmetical Areal Surficial Mean HeightSaμm(28)(Table 6)
Areal Surficial SkewnessSsk-(29), (35)(Table 6, Table 9 and Table 10)
Areal Surficial KurtosisSku-(30), (36)(Table 6, Table 9 and Table 10)
True packing densityPDtr (PD)(%)(31), (32)(Table 6 and Table 7)
Theoretical packing densityPDth(%)--
LayerReal layer height/thicknesshlμm(22), (23)-
Theoretical layer height/thicknesshl,thμm(22)-
Material propertiesShear modulusGPa(1)(Table 3 and Table 7)
Young’s modulusEPa(18), (21)(Table 3 and Table 4)
Poisson ratiov-(2), (18), (21)(Table 3 and Table 7)
Powder bulk densityρkg/m3(1), (32)(Table 2 and Table 7)
Coefficient of restitution of materials A–Bcr,A-B---
Coefficient of rolling friction of materials A–Bcrf,A-B---
Coefficient of static friction of materials A–Bcsf,A-B---
Equilibrium separation in Lennard–Jones potential z 0 Å (8), (10)(Table 3)
Particle size and distributionParticle diameter distribution percentile (X%)DXμm(14)–(16)(Table 2)
Lognormal distribution’s mean μ Χ μm(11), (12)-
Lognormal distribution’s standard deviation σ Χ μm(12)-
Cumulative lognormal distribution function F Χ -(13)–(16)-
Particle contact modelEffective Young’s modulus of contacting particlesE*Pa(5)–(8), (18)(Table 3)
Effective radius at the point of interparticle contactR*m(4)–(8), (10), (17), (20)(Table 3 and Table 4)
Work of adhesionΓJ/m2(4), (6)–(10), (20), (21)(Table 3 and Table 4)
Contact patch radiusαμm(6), (7)-
Relative particle approachδnμm(3), (7), (10)-
Relative particle approach at tear-off pointδtoμm(3), (5), (19), (21)(Table 3 and Table 4)
Normal force of particles in contact (JKR model) F n J K R N(3), (6)-
Pull-off force of particles in contact (JKR model) F p o J K R N(3), (4)-
Pull-off forceFpoN(5), (20)(Table 4)
Tabor’s parameterμ-(8)(Table 1, Table 3 and Table 4)
Tensile loadTN(10)-
SimulationMinimum particle radius in simulationrminμm(1)(Table 7)
Rayleigh critical timestepΔtcnsec(1)(Table 7)
Physical radiusRparticleμm(19)-
Contact radiusRc,EDEMμm(19)-
Average particle radius (lognormal normalization)Raveμm-(Table 3)
Taguchi-DoETranslational deposition speed of the doctor bladeutrm/sec(35), (36)(Table 6, Table 8, Table 9 and Table 10)
Vertical vibrational frequency of the doctor bladefvibHz(35), (36)(Table 6, Table 8, Table 9 and Table 10)
Vertical vibrational amplitude of the doctor bladeAvibμm(35), (36)(Table 6, Table 8, Table 9 and Table 10)
Angle of relief of the doctor bladeθreldeg.(35), (36)(Table 6, Table 8, Table 9 and Table 10)

References

  1. Avrampos, P.; Vosniakos, G.-C. A prototype powder deposition system for an open Selective Laser Sintering machine. Procedia Manuf. 2020, 51, 755–762. [Google Scholar] [CrossRef]
  2. Avrampos, P. Optimized Development of a Prototype Selective Laser Sintering Powder Recoating System via Analytic Hierarchy Process. Int. J. Exp. Des. Process Optim. 2023, in press. [Google Scholar] [CrossRef]
  3. Haeri, S.; Wang, Y.; Ghita, O.; Sun, J. Discrete element simulation and experimental study of powder spreading process in additive manufacturing. Powder Technol. 2017, 306, 45–54. [Google Scholar] [CrossRef]
  4. Budding, A.; Vaneker, T.H.J. New Strategies for Powder Compaction in Powder-based Rapid Prototyping Techniques. Procedia CIRP 2013, 6, 527–532. [Google Scholar] [CrossRef]
  5. Haeri, S. Optimisation of blade type spreaders for powder bed preparation in Additive Manufacturing using DEM simulations. Powder Technol. 2017, 321, 94–104. [Google Scholar] [CrossRef]
  6. Shamsdini, S.; Ghoncheh, M.H.; Mohammadi, M. Effect of recoater-blade type on the mechanical properties and microstructure of additively manufactured maraging steels. Mater. Sci. Eng. A 2021, 812, 141104. [Google Scholar] [CrossRef]
  7. Zhang, J.; Tan, Y.; Bao, T.; Xu, Y.; Xiao, X.; Jiang, S. Discrete Element Simulation of the Effect of Roller-Spreading Parameters on Powder-Bed Density in Additive Manufacturing. Materials 2020, 13, 2285. [Google Scholar] [CrossRef]
  8. Avrampos, P.; Vosniakos, G.-C. A review of powder deposition in additive manufacturing by powder bed fusion. J. Manuf. Process. 2022, 74, 332–352. [Google Scholar] [CrossRef]
  9. Spierings, A.B.; Herres, N.; Levy, G. Influence of the particle size distribution on surface quality and mechanical properties in AM steel parts. Rapid Prototyp. J. 2011, 17, 195–202. [Google Scholar] [CrossRef]
  10. Leicht, A.; Fischer, M.; Klement, U.; Nyborg, L.; Hryha, E. Increasing the Productivity of Laser Powder Bed Fusion for Stainless Steel 316L through Increased Layer Thickness. J. Mater. Eng. Perform. 2021, 30, 575–584. [Google Scholar] [CrossRef]
  11. Ziegelmeier, S.; Christou, P.; Wöllecke, F.; Tuck, C.; Goodridge, R.; Hague, R.; Krampe, E.; Wintermantel, E. An experimental study into the effects of bulk and flow behaviour of laser sintering polymer powders on resulting part properties. J. Mater. Process. Technol. 2015, 215, 239–250. [Google Scholar] [CrossRef]
  12. Tan, P.; Shen, F.; Tey, W.S.; Zhou, K. A numerical study on the packing quality of fibre/polymer composite powder for powder bed fusion additive manufacturing. Virtual Phys. Prototyp. 2021, 16 (Suppl. S1), S1–S18. [Google Scholar] [CrossRef]
  13. Wang, L.; Yu, A.; Li, E.; Shen, H.; Zhou, Z. Effects of spreader geometry on powder spreading process in powder bed additive manufacturing. Powder Technol. 2021, 384, 211–222. [Google Scholar] [CrossRef]
  14. Burns, S.J.; Piiroinen, P.T.; Hanley, K.J. Critical time step for DEM simulations of dynamic systems using a Hertzian contact model. Int. J. Numer. Methods Eng. 2019, 119, 432–451. [Google Scholar] [CrossRef]
  15. Kremmer, M.; Favier, J.F. A method for representing boundaries in discrete element modelling—Part II: Kinematics. Int. J. Numer. Methods Eng. 2001, 51, 1423–1436. [Google Scholar] [CrossRef]
  16. Kafui, K.D.; Thornton, C.; Adams, M.J. Discrete particle-continuum fluid modelling of gas–solid fluidised beds. Chem. Eng. Sci. 2002, 57, 2395–2410. [Google Scholar] [CrossRef]
  17. Thornton, C.; Randall, C.W. Applications of Theoretical Contact Mechanics to Solid Particle System Simulation. In Micromechanics of Granular Materials, Proceedings of the U.S./Japan Seminar on the Micromechanics of Granular Materials, Sendai, Japan, 26–30 October 1987; Satake, M., Jenkins, J.T., Eds.; Elsevier: Amsterdam, The Netherlands, 1988; Volume 20, pp. 133–142. [Google Scholar] [CrossRef]
  18. Li, Y.; Xu, Y.; Thornton, C. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles. Powder Technol. 2005, 160, 219–228. [Google Scholar] [CrossRef]
  19. Bauccio, M. ASM Engineered Materials Reference Book, 2nd ed.; ASM International: Novelty, OH, USA, 1994. [Google Scholar]
  20. AZO Materials. Alumina as a Biomaterial (99.5% Alumina). Available online: https://www.azom.com/properties.aspx?ArticleID=105 (accessed on 15 November 2023).
  21. AZO Materials. Alumina—Aluminium Oxide—Al2O3—A Refractory Ceramic Oxide. Available online: https://www.azom.com/properties.aspx?ArticleID=52 (accessed on 28 November 2023).
  22. Chen, H.; Xiao, Y.G.; Liu, Y.L.; Shi, Y.S. Effect of Young’s modulus on DEM results regarding transverse mixing of particles within a rotating drum. Powder Technol. 2017, 318, 507–517. [Google Scholar] [CrossRef]
  23. Thornton, C. Granular Dynamics, Contact Mechanics and Particle System Simulations; Springer International Publishing: Cham, Switzerland, 2015; Volume 24. [Google Scholar] [CrossRef]
  24. Johnson, K.L.; Kendall, K.; Roberts, A.D. Surface energy and the contact of elastic solids. Proc. R. Soc. Lond. A Math. Phys. Sci. 1971, 324, 301–313. [Google Scholar] [CrossRef]
  25. Derjaguin, B.V.; Muller, V.M.; Toporov, Y. Effect of contact deformations on the adhesion of particles. J. Colloid Interface Sci. 1975, 53, 314–326. [Google Scholar] [CrossRef]
  26. Greenwood, J.A. Adhesion of Elastic Spheres. Proc. Math. Phys. Eng. Sci. 1997, 453, 1277–1297. [Google Scholar] [CrossRef]
  27. Maugis, D. Adhesion of spheres: The JKR-DMT transition using a dugdale model. J. Colloid Interface Sci. 1992, 150, 243–269. [Google Scholar] [CrossRef]
  28. Edem, A. The Hertz-Mindlin with JKR Version 2 Model. Available online: https://2022.help.altair.com/2022.2/EDEM/Creator/Physics/Base_Models/Hertz-Mindlin_with_JKR_V2.htm (accessed on 27 March 2024).
  29. Grierson, D.S.; Flater, E.E.; Carpic, R.W.K. Accounting for the JKR–DMT transition in adhesion and friction measurements with atomic force microscopy. J. Adhes. Sci. Technol. 2005, 19, 291–311. [Google Scholar] [CrossRef]
  30. Shi, X.; Zhao, Y.-P. Comparison of various adhesion contact theories and the influence of dimensionless load parameter. J. Adhes. Sci. Technol. 2004, 18, 55–68. [Google Scholar] [CrossRef]
  31. Mukherjee, R.; Gupta, V.; Naik, S.; Sarkar, S.; Sharma, V.; Peri, P.; Chaudhuri, B. Effects of particle size on the triboelectrification phenomenon in pharmaceutical excipients: Experiments and multi-scale modeling. Asian J. Pharm. Sci. 2016, 11, 603–617. [Google Scholar] [CrossRef]
  32. Ose, S. Dusting of Alumina and Electrostatic Charging. Part. Sci. Technol. 2003, 21, 237–245. [Google Scholar] [CrossRef]
  33. Tavakoli, A.H.; Maram, P.S.; Widgeon, S.J.; Rufner, J.; van Benthem, K.; Ushakov, S.; Sen, S.; Navrotsky, A. Amorphous Alumina Nanoparticles: Structure, Surface Energy, and Thermodynamic Phase Stability. J. Phys. Chem. C 2013, 117, 17123–17130. [Google Scholar] [CrossRef]
  34. Vollath, D.; Fischer, F.D.; Holec, D. Surface energy of nanoparticles—Influence of particle size and structure. Beilstein J. Nanotechnol. 2018, 9, 2265–2276. [Google Scholar] [CrossRef] [PubMed]
  35. Tepesch, P.D.; Quong, A.A. First-Principles Calculations of α-Alumina (0001) Surfaces Energies with and without Hydrogen. Phys. Status Solidi 2000, 217, 377–387. [Google Scholar] [CrossRef]
  36. Vilfan, I.; Deutsch, T.; Lançon, F.; Renaud, G. Erratum to: ‘Structure determination of the (3×3) reconstructed α-Al2O3(0001)’ [Surf. Sci. 505 (2002) L215]. Surf. Sci. 2003, 529, 281. [Google Scholar] [CrossRef]
  37. ASTM D6773-16; Standard Shear Test Method for Bulk Solids Using the Schulze Ring Shear Tester. ASTM International: West Conshohocken, PA, USA, 2008. [CrossRef]
  38. Tahir, A.; Rasche, S.; Könke, C. Discrete element model development of ZTA ceramic granular powder using micro computed tomography. Adv. Powder Technol. 2018, 29, 3471–3482. [Google Scholar] [CrossRef]
  39. Schiller, L.; Naumann, A. A Drag Coefficient Correlation. Z. Ver. Dtsch. Ingenieure 1935, 77, 318–320. [Google Scholar]
  40. Karunarathne, S.S.; Tokheim, L.-A. Comparison of the influence of drag models in CFD simulation of particle mixing and segregation in a rotating cylinder. In Proceedings of the 58th Conference on Simulation and Modelling (SIMS 58), Reykjavik, Iceland, 25–27 September 2017; pp. 151–156. [Google Scholar] [CrossRef]
  41. AZO Materials. Stainless Steel—Grade 304 (UNS S30400). 2023. Available online: https://www.azom.com/properties.aspx?ArticleID=965 (accessed on 27 December 2023).
  42. Pasalopoulos, S.; Avrampos, P.; Vosniakos, G.-C. Surface quality evaluation of non-sintered powder layers in Selective Laser Sintering by 3D scanning. Procedia Manuf. 2020, 51, 748–754. [Google Scholar] [CrossRef]
  43. Taylor, S.; Forrest, E.C.; Jared, B.H. Investigating Applicability of Surface Roughness Parameters in Describing the Metallic Additive Manufacturing Process. In Proceedings of the Solid Freeform Fabrication Symposium, Austin, TX, USA, 11–14 August 2019; Available online: https://www.osti.gov/biblio/1641526 (accessed on 9 May 2024).
  44. Michigan Metrology LLC. Ssk (Skewness) and Sku (Kurtosis). 2023. Available online: https://michmet.com/glossary-term/kurtosis/ (accessed on 22 December 2023).
  45. OLYMPUS-IMS-EVIDENT. Surface Roughness Measurement—Parameters. 2023. Available online: https://www.olympus-ims.com/en/metrology/surface-roughness-measurement-portal/parameters/ (accessed on 23 December 2023).
  46. Ziri, S.; Hor, A.; Mabru, C. Effect of powder size and processing parameters on surface, density and mechanical properties of 316L elaborated by Laser Powder Bed Fusion. In Proceedings of the ESAFORM 2021 24th International Conference on Material Forming, Liege, Belgium, 14–16 April 2021. [Google Scholar] [CrossRef]
  47. Fouda, Y.M.; Bayly, A.E. A DEM study of powder spreading in additive layer manufacturing. Granul. Matter 2020, 22, 10. [Google Scholar] [CrossRef]
  48. Wu, S.; Yang, Y.; Huang, Y.; Han, C.; Chen, J.; Xiao, Y.; Li, Y.; Wang, D. Study on powder particle behavior in powder spreading with discrete element method and its critical implications for binder jetting additive manufacturing processes. Virtual Phys. Prototyp. 2023, 18, e2158877. [Google Scholar] [CrossRef]
  49. Roy, R. A Primer on the Taguchi Method, 1st ed.; Society of Manufacturing Engineers: New York, NY, USA, 1990. Available online: https://catalogue.nla.gov.au/Record/1852247 (accessed on 21 April 2024).
  50. Nasato, D.S.; Briesen, H.; Pöschel, T. Influence of vibrating recoating mechanism for the deposition of powders in additive manufacturing: Discrete element simulations of polyamide 12. Addit. Manuf. 2021, 48, 102248. [Google Scholar] [CrossRef]
Figure 1. Visual depiction of the work of adhesion between the two surfaces and order of contact evolution [28].
Figure 1. Visual depiction of the work of adhesion between the two surfaces and order of contact evolution [28].
Jmmp 08 00101 g001
Figure 2. Surface area vs. particle diameter plot for spherical a-alumina particles.
Figure 2. Surface area vs. particle diameter plot for spherical a-alumina particles.
Jmmp 08 00101 g002
Figure 3. Tabor parameter curve vs. the second particle’s diameter for two extreme diameters of the first particle.
Figure 3. Tabor parameter curve vs. the second particle’s diameter for two extreme diameters of the first particle.
Jmmp 08 00101 g003
Figure 4. (a) The physical bodies of the powder spreading simulation; (b) semi-transparent cyan cube, which functions as a powder factory.
Figure 4. (a) The physical bodies of the powder spreading simulation; (b) semi-transparent cyan cube, which functions as a powder factory.
Jmmp 08 00101 g004
Figure 5. (Left) The initial moment of powder generation; (Middle) particles’ free-fall onto the substrate, over the sample square; (Right) particles’ free-fall complete. The slope has been created/deposition is ready to begin.
Figure 5. (Left) The initial moment of powder generation; (Middle) particles’ free-fall onto the substrate, over the sample square; (Right) particles’ free-fall complete. The slope has been created/deposition is ready to begin.
Jmmp 08 00101 g005
Figure 6. Sample top view of the deposited powder layer, with the coloring of the particles via (a) vertical position and (b) diameter.
Figure 6. Sample top view of the deposited powder layer, with the coloring of the particles via (a) vertical position and (b) diameter.
Jmmp 08 00101 g006
Figure 7. Flowchart of the top-surface layer calculation code.
Figure 7. Flowchart of the top-surface layer calculation code.
Jmmp 08 00101 g007
Figure 8. Black circles are uncovered areas where the substrate is visible, denoting coverage defects.
Figure 8. Black circles are uncovered areas where the substrate is visible, denoting coverage defects.
Jmmp 08 00101 g008
Figure 9. Infinitesimal prism visualization.
Figure 9. Infinitesimal prism visualization.
Jmmp 08 00101 g009
Figure 10. Powder deposition example (top and side views) with a non-vibrating doctor blade at deposition speed of: (a) 0.08 m/s; (b) 0.01 m/s.
Figure 10. Powder deposition example (top and side views) with a non-vibrating doctor blade at deposition speed of: (a) 0.08 m/s; (b) 0.01 m/s.
Jmmp 08 00101 g010
Figure 11. Side view of the doctor blade. Depiction of the angle of relief and the blade’s profile.
Figure 11. Side view of the doctor blade. Depiction of the angle of relief and the blade’s profile.
Jmmp 08 00101 g011
Figure 12. Means of the quality criteria: (a) LTD; (b) SCR; (c) RMS; (d) PD.
Figure 12. Means of the quality criteria: (a) LTD; (b) SCR; (c) RMS; (d) PD.
Jmmp 08 00101 g012
Figure 13. Signal-to-noise ratios: (a) LTD; (b) SCR; (c) RMS; (d) PD.
Figure 13. Signal-to-noise ratios: (a) LTD; (b) SCR; (c) RMS; (d) PD.
Jmmp 08 00101 g013
Figure 14. Interaction plots for means of |LTD|: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Figure 14. Interaction plots for means of |LTD|: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Jmmp 08 00101 g014
Figure 15. Interaction plots for means of SCR: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Figure 15. Interaction plots for means of SCR: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Jmmp 08 00101 g015
Figure 16. Interaction plots for means of Sq-RMS: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Figure 16. Interaction plots for means of Sq-RMS: (a) utr-fvib; (b) utr-Avib; (c) utrrel.
Jmmp 08 00101 g016
Figure 17. Quality criteria vs. Ssk.
Figure 17. Quality criteria vs. Ssk.
Jmmp 08 00101 g017
Figure 18. Quality criteria vs. Sku.
Figure 18. Quality criteria vs. Sku.
Jmmp 08 00101 g018
Figure 19. Best-fit line WMSsk versus Ssk.
Figure 19. Best-fit line WMSsk versus Ssk.
Jmmp 08 00101 g019
Figure 20. Best-fit line WMSku versus Sku.
Figure 20. Best-fit line WMSku versus Sku.
Jmmp 08 00101 g020
Figure 21. Layer quality comparison: (a) vibration-less (top and side); (b) vibration-assisted optimum (top and side). (Note: Magenta: back border plate; pink: front border plate).
Figure 21. Layer quality comparison: (a) vibration-less (top and side); (b) vibration-assisted optimum (top and side). (Note: Magenta: back border plate; pink: front border plate).
Jmmp 08 00101 g021
Table 1. Particle contact model applicability based on Tabor parameter range.
Table 1. Particle contact model applicability based on Tabor parameter range.
JKR
[24]
DMT
[25]
BD
[26]
MD
[27]
Thornton [23]μ > 5μ < 0.1Not mentioned0.1 ≤ μ ≤ 5
Shi [30]μ ≥ 5μ ≤ 0.1Not mentionedAny μ
Greenwood [26]μ ≥ 3Wrong in theory and in practiceAny μAny μ
Table 2. Specifications of the powder used in the simulations.
Table 2. Specifications of the powder used in the simulations.
Item Unit Typical Value
Particle sizeD10μm8.2
D50μm21
D90μm47.5
Specific Surface Area (S.S.A.)m2/g0.19
Electrical Conductivity (E.C.)μs/cm300
pH-8.5
Moisture%0.05
True Densityg/cm33.8
Spheroidization%96
Chemical
Composition
Al2O3%99.5
Fe2O3ppm300
Na2Oppm3500
Ion contentNa+ppm400
Table 3. Calculated relative approach and Tabor parameter for different shear modulus values.
Table 3. Calculated relative approach and Tabor parameter for different shear modulus values.
C A S E G
( P a )
ν E
( P a )
E *
( P a )
R a v e
( μ m )
R *
μ m
Γ J m 2 z 0
( Å )
δ t o
( μ m )
μ
S I M 1 1 × 107 0.3 3.2 × 1071.76 × 107 25 12.5 4.26 5.1 1.451 1768.56
S I M 2 1.5 × 108 0.3 4.8 × 1082.64 × 108 25 12.5 4.26 5.1 0.239 290.78
R W 3 1.29 × 1011 0.3 4.13 × 10112.27 × 1011 25 12.5 4.26 5.1 0.003 3.22
Table 4. The relation of parameter values between the real world and the simulation-defined material. In parentheses, it is stated how this parameter affects the powder’s behavior compared to its counterpart.
Table 4. The relation of parameter values between the real world and the simulation-defined material. In parentheses, it is stated how this parameter affects the powder’s behavior compared to its counterpart.
REAL-WORLD PARTICLES SIMULATION PARTICLES
Young’s modulus↑ (less elasticity)↓ (more elasticity)
Tabor’s parameter<5 (less elasticity)>100 (more elasticity)
Suitable model for F p o Bradley 2 π R Γ
(more cohesion)
JKR 3 / 2 π R Γ
(less cohesion)
Relative approach to bond breakage δ t o ( μ m ) 3 nm (less cohesion)1451 nm (more cohesion)
Table 5. Structure of the exported EDEM file for the particles constituting the layer of the simulation trials.
Table 5. Structure of the exported EDEM file for the particles constituting the layer of the simulation trials.
P a r t i c l e   I D D i a m e t e r X c e n t r e Y c e n t r e Z c e n t r e
11.39 × 10−5−0.0001−0.000170.000507
31.98 × 10−5−0.00038−0.00020.000531
20,6093.51 × 10−5−0.00013−0.000350.000531
Table 6. L27 Taguchi orthogonal array for powder spreading.
Table 6. L27 Taguchi orthogonal array for powder spreading.
# u t r m s f v i b
H z
A v i b
μ m
θ r e l
d e g .
L T D
( μ m )
L T D
( μ m )
S C R
( % )
S q
( μ m )
P D t r
%
S a
( μ m )
S s k S k u
10.015001.00−29.029.098.919.968.415.5−0.6723.889
20.015002.55−26.026.099.117.968.013.6−0.9064.852
30.015005.010−21.221.299.317.467.113.0−1.1795.572
40.0110001.05−26.726.799.118.167.613.8−0.9564.741
50.0110002.510−24.324.399.018.067.413.5−1.3515.626
60.0110005.00−22.222.299.316.168.012.1−1.5156.403
70.0120001.010−23.823.899.017.967.013.6−1.3385.570
80.0120002.50−23.023.099.117.067.812.8−1.5036.075
90.0120005.05−21.521.599.316.467.812.3−1.4316.232
100.055001.05−36.636.697.721.667.217.0−0.6853.347
110.055002.510−34.934.997.722.466.817.6−0.5773.362
120.055005.00−31.731.798.618.867.514.6−1.0424.358
130.0510001.010−36.936.997.920.466.816.1−0.8113.602
140.0510002.50−33.833.898.519.267.015.0−0.9263.993
150.0510005.05−32.132.198.419.967.115.5−0.8983.993
160.0520001.00−35.035.098.120.267.115.9−0.8693.698
170.0520002.55−32.432.498.619.267.114.9−0.9624.116
180.0520005.010−32.232.298.520.067.115.5−0.8153.956
190.105001.010−40.840.897.022.867.318.2−0.3993.051
200.105002.50−41.241.296.921.867.917.4−0.5973.059
210.105005.05−36.936.997.522.067.017.3−0.5963.310
220.1010001.00−41.441.497.021.767.617.3−0.6333.063
230.1010002.55−40.740.796.124.967.420.1−0.3232.732
240.1010005.010−37.037.094.928.766.922.9−0.2242.753
250.1020001.05−40.240.297.321.367.317.0−0.6223.169
260.1020002.510−41.141.197.322.567.518.0−0.3352.940
270.1020005.00−40.240.297.221.367.516.8−0.5583.604
Table 7. Calculated timesteps and surface quality criteria for various G values. The settings used were those of trial #9 of the Taguchi L27 array for powder spreading, i.e., U t r = 0.01   m / s , f v i b = 2000   H z , A v i b = 5   μ m and θ r e l = 5 ° .
Table 7. Calculated timesteps and surface quality criteria for various G values. The settings used were those of trial #9 of the Taguchi L27 array for powder spreading, i.e., U t r = 0.01   m / s , f v i b = 2000   H z , A v i b = 5   μ m and θ r e l = 5 ° .
REP
#
r m i n
μ m
ρ K g m 3 ν β G
P a
Δ t c
n s e c
t s i m
( h )
L T D
( μ m )
S C R
( % )
S q
( μ m )
P D
( %   o f   B D )
12.538200.30.9255107165.8612.521.5 (+18.8%)99.3
(−0.1%)
16.4
(+10%)
67.8
(+0.6%)
22.538200.30.92551.5 × 10842.8252518.9
(+4.4%)
99.4
(≡)
15.8
(+6%)
67.4
(≡)
32.538200.30.92551.5 × 10111.35416818.199.414.967.4
Table 8. Regression vs. simulation performance of optimum parameter sets for the three different quality criteria.
Table 8. Regression vs. simulation performance of optimum parameter sets for the three different quality criteria.
Optimum utr (m/s) fvib (Hz) Avib (μm) θrel (°) Regression Optimum Simulation Optimum Deviation
(sim-reg)/reg
ANOVA Regression Error
|LTD|0.01200051020.8 μm20.4 μm−1.9%5.32%
SCR0.0120005099.54%99.33%−0.2%15.52%
Sq-RMS0.0120005015.8 μm15.7 μm−0.6%15.67%
Table 9. Performance comparison of optimum |LTD|, SCR and Sq values (reg: “regression”; sim: “simulation”).
Table 9. Performance comparison of optimum |LTD|, SCR and Sq values (reg: “regression”; sim: “simulation”).
OPTIMUM utr
(m⁄s)
fvib
(Hz)
Avib
(μm)
θrel
(deg)
|LTD|reg
(|LTD|sim)
SCRreg
(SCRsim)
Sq-reg
(Sq-sim)
SsksimSkusim
|LTD|0.01200051020.8
(20.4)
99.51
(99.29)
15.9
(16.5)
−1.4126.304
SCR0.0120005021.4
(22.1)
99.54
(99.33)
15.8
(15.7)
−1.5436.530
Sq-RMS0.0120005021.4
(22.1)
99.54
(99.33)
15.8
(15.7)
−1.5436.530
Table 10. Vibration-less vs. vibration-assisted performance (reg: “regression”; sim: “simulation”).
Table 10. Vibration-less vs. vibration-assisted performance (reg: “regression”; sim: “simulation”).
utr
(m⁄s)
fvib
(Hz)
Avib
(μm)
θrel
(deg)
|LTD|reg
(|LTD|sim)
SCRreg
(SCRsim)
Sq-reg
(Sq-sim)
SsksimSkusim
NO VIB.0.01000-
(28.9)
-
(98.56)
-
(20.9)
−0.7163.962
OPTIMUM0.0120005021.4
(22.1)
99.54
(99.33)
15.8
(15.7)
−1.5436.530
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Avrampos, P.; Vosniakos, G.-C. A Study on Powder Spreading Quality in Powder Bed Fusion Processes Using Discrete Element Method Simulation. J. Manuf. Mater. Process. 2024, 8, 101. https://doi.org/10.3390/jmmp8030101

AMA Style

Avrampos P, Vosniakos G-C. A Study on Powder Spreading Quality in Powder Bed Fusion Processes Using Discrete Element Method Simulation. Journal of Manufacturing and Materials Processing. 2024; 8(3):101. https://doi.org/10.3390/jmmp8030101

Chicago/Turabian Style

Avrampos, Panagiotis, and George-Christopher Vosniakos. 2024. "A Study on Powder Spreading Quality in Powder Bed Fusion Processes Using Discrete Element Method Simulation" Journal of Manufacturing and Materials Processing 8, no. 3: 101. https://doi.org/10.3390/jmmp8030101

APA Style

Avrampos, P., & Vosniakos, G. -C. (2024). A Study on Powder Spreading Quality in Powder Bed Fusion Processes Using Discrete Element Method Simulation. Journal of Manufacturing and Materials Processing, 8(3), 101. https://doi.org/10.3390/jmmp8030101

Article Metrics

Back to TopTop