Material Evaluation and Dynamic Powder Deposition Modeling of PEEK/CF Composite for Laser Powder Bed Fusion Process

Polymeric composites such as Poly-ether-ether-ketone (PEEK)/carbon fiber (CF) have been widely utilized due to outstanding performances such as high specific strength and specific modulus. The PEEK/CF components via powder bed fusion additive manufacturing usually show brittle fracture behaviors induced by their poor interfacial affinity and inner voids. These defects are strongly associated with powder packing quality upon deposition. The particle dynamic model has been widely employed to study the interactions of particle motions. Powder property, bulk material property, and interfacial features of composite powders are key factors in the particle dynamic model. In this work, an efficient and systematic material evaluation is developed for composite powders to investigate their deposition mechanism. The discrete element method is utilized to simulate the dynamic behaviors of PEEK/CF composite powders. The powder properties, bulk material properties, and interfacial features of powders are calibrated and justified by experimental measurement, numerical simulation, and design of experiments. The particle dynamic model can explain the powder flow behaviors and interactions. The experimental and simulation AOR results show a maximal deviation of 4.89%. It reveals that the addition of short CF particles can assist the flow of PEEK powders and improve the packing quality of the composite powders. The results show an experimental improvement of 31.3% and 55.2% for PEEK/CF_30wt% and PEEK/CF_50wt%, with a simulated improvement of 27.4% and 50.2% for corresponding composite powders.


Introduction
Laser powder bed fusion (LPBF) technology, as an essential category of multiple additive manufacturing technologies, is able to produce polymeric composite materials economically and with obvious flexibility [1,2]. This technique is widely used in aerospace [3,4], automotive manufacturing [5,6], optical components [7][8][9], and other industrial fields. The technique of LPBF uses a high-energy intensity laser beam to selectively melt the powder, which can fabricate more complex and superior structures due to its smaller laser spot and thinner layer thickness [10][11][12]. Polymeric composites such as Poly-ether-ether-ketone (PEEK)/carbon fiber (CF, specific chopped fiber in this work) have attracted the attention of researchers as well as industrialists due to outstanding performances such as high specific strength and specific modulus. However, the laser-sintered PEEK/CF components are usually brittle because of their poor interfacial affinity and inner voids. These defects are strongly associated with the powder packing quality upon deposition. Thus, the investigation of the powder deposition process incorporating the detailed material properties is

Contact Model Theory
The discrete element method (DEM) is utilized to simulate the dynamic behaviors of PEEK/CF composite powders. In the DEM simulation, the particles are modeled as discrete entities to describe the dynamics of particle flow. The translation, rotation, and location of the particles can be characterized by Newton's laws of motion, which is described as follows [27,28]: where m i , ω i , v i , and I i are the mass, angular velocity, translational velocity, and moment of inertia, respectively. F c,i is the contact force, M c,i is the contact torque, and R i is the rotation matrix. The relationship between cohesive force and particle size is shown in Figure 1a [29,30]. It appears that the cohesive forces such as Van der Waals force, electrostatic force, and capillary force are greater than the gravity of particles when the particle size is less than 100 µm, which reduces the powder fluidity and causes inconvenience for the powderspreading process. In reality, the powder is preheated and dried before the powderspreading process. Thus, the capillary force can be ignored.

Contact Model Theory
The discrete element method (DEM) is utilized to simulate the dynamic behaviors of PEEK/CF composite powders. In the DEM simulation, the particles are modeled as discrete entities to describe the dynamics of particle flow. The translation, rotation, and location of the particles can be characterized by Newton's laws of motion, which is described as follows [27,28]: (2) where mi, ωi, vi, and Ii are the mass, angular velocity, translational velocity, and moment of inertia, respectively. Fc,i is the contact force, Mc,i is the contact torque, and Ri is the rotation matrix.
The relationship between cohesive force and particle size is shown in Figure 1a [29,30]. It appears that the cohesive forces such as Van der Waals force, electrostatic force, and capillary force are greater than the gravity of particles when the particle size is less than 100 µm, which reduces the powder fluidity and causes inconvenience for the powder-spreading process. In reality, the powder is preheated and dried before the powderspreading process. Thus, the capillary force can be ignored.
(a) (b) (c) Figure 1. Interaction between particles and particle modeling method: (a) the relationship between particle force and particle diameter [30]; (b) contact force diagram between particles; (c) multiplesphere method to model the CF particle. Figure 1. Interaction between particles and particle modeling method: (a) the relationship between particle force and particle diameter [30]; (b) contact force diagram between particles; (c) multiplesphere method to model the CF particle. (1) Hertz-Mindlin contact model The rationality of the Hertz-Mindlin contact model has been widely recognized in the particle-related research field [31,32]. In order to accurately characterize the stress-strain relationship between particles considering adhesion, this work adopts the Hertz-Mindlin and Johnson-Kendall-Roberts (JKR) model for polymers and their composite powders, which is implemented using EDEM software. Among them, the Hertz-Mindlin model is the foundation of the DEM method, and the parameter of surface energy in the JKR model can comprehensively characterize the bonding effect between particles caused by the van der Waals force and the electrostatic force [33]. Therefore, the Hertz-Mindlin and JKR models possess universality in the DEM simulation of other thermoplastic polymers and their composite powders.
In the Hertz-Mindlin model, the normal and tangential force components between particles are calculated based on the Hertzian contact theory and the Middlin-Deresiewicz method, respectively [34,35]. Both normal and tangential force between particles have damping components. The damping coefficient is related to the restitution coefficient [36]. The tangential friction force between particles follows Coulomb's law of friction, and the rolling friction between particles is computed based on the Contact-Directional-Constant-Torque model [27,37]. Figure 2 shows the schematic diagram of the Hertz-Mindlin model for adjacent particles i and j, where the normal contact force on particle i can be expressed as where F CN,E is the normal elastic force, and F CN,D is the normal damping force. The F CN,E is computed based on the Hertz contact theory, which is described as where α n is the normal contact displacement (overlap of adjacent particles). The Contact constant K n is expressed as where G is the shear modulus of the particle, r is the radius of the local curvature, and υ is the Poisson's ratio. (1) Hertz-Mindlin contact model The rationality of the Hertz-Mindlin contact model has been widely recognized in the particle-related research field [31,32]. In order to accurately characterize the stressstrain relationship between particles considering adhesion, this work adopts the Hertz-Mindlin and Johnson-Kendall-Roberts (JKR) model for polymers and their composite powders, which is implemented using EDEM software. Among them, the Hertz-Mindlin model is the foundation of the DEM method, and the parameter of surface energy in the JKR model can comprehensively characterize the bonding effect between particles caused by the van der Waals force and the electrostatic force [33]. Therefore, the Hertz-Mindlin and JKR models possess universality in the DEM simulation of other thermoplastic polymers and their composite powders.
In the Hertz-Mindlin model, the normal and tangential force components between particles are calculated based on the Hertzian contact theory and the Middlin-Deresiewicz method, respectively [34,35]. Both normal and tangential force between particles have damping components. The damping coefficient is related to the restitution coefficient [36]. The tangential friction force between particles follows Coulomb's law of friction, and the rolling friction between particles is computed based on the Contact-Directional-Constant-Torque model [27,37]. Figure 2 shows the schematic diagram of the Hertz-Mindlin model for adjacent particles i and j, where the normal contact force on particle i can be expressed as where FCN,E is the normal elastic force, and FCN,D is the normal damping force. The FCN,E is computed based on the Hertz contact theory, which is described as where αn is the normal contact displacement (overlap of adjacent particles). The Contact constant Kn is expressed as where G is the shear modulus of the particle, r is the radius of the local curvature, and υ is the Poisson's ratio. The normal damping force F CN,D between particles is expressed as where v re,n is the component of the relative velocity of particle i relative to particle j along the normal direction. k n is the normal contact stiffness, β * is the damping coefficient, and m * is the equivalent mass. Therein, the damping coefficient and normal contact stiffness can be expressed as where e is the restitution coefficient between particles. The equivalent mass m * and equivalent curvature radius r of particles are, respectively, described as where m i and m j represent the masses of adjacent particles i and j. r i and r j are the radiuses of adjacent particles i and j.
The tangential contact force F CT between adjacent particles can be expressed as where F CT,E is the elastic force along the tangential direction, and F CT,D is the damping force along the tangential direction. The tangential elastic force F CT,E can be expressed in incremental form as where F CT,E,(n−1) is the tangential elastic force of the previous analysis step; ∆F CT,E is the increment of the tangential elastic force in the current analysis step, which is described as where ∆α t is the tangential displacement increment, and k t is the tangential spring stiffness, expressed as where v re,t is the relative tangential velocity between particles, and ∆t is the step time. Meanwhile, the tangential damping force between adjacent particles can be expressed as It is worth noting that the tangential force between adjacent particles is related to the frictional force between particles (µ s F CN , µ s is the static friction coefficient). In the simulation process, the rolling friction cannot be ignored, which is achieved by applying a torque to the contact surface of the particles: where R i is the distance from the contact point to the mass center, µ r is the rolling friction coefficient, and ω i is the unit angular velocity vector at the contact point. (2) JKR theory The Hertz-Mindlin with JKR contact model can form an adhesive force contact model between particles, mainly considering the influence of van der Waals forces in the contact area. In this model, the normal elastic contact force between particles is computed based on the JKR theory [38], including the overlap of particles δ, interaction parameters, and surface energy γ, which can be described as 1 where E * is the equivalent Young's modulus, r * is the equivalent curvature radius, and α is the contact radius of particles. It is discovered that several parameters are involved in the aforementioned model, such as static friction coefficient, rolling friction coefficient, restitution coefficient, and surface energy. Therefore, it is essential to carry out the parameter calibration work through the material characterization.

Parameter Classification of the Discrete Element Model
The first step of the DEM simulation is to calibrate numerous parameters through comprehensive material characterization. Figure 3 shows the parameters involved in the DEM model, which can be divided into three categories: (1) Solid-phase parameters: solid phase density, elastic modulus, and Poisson's ratio; (2) Contact parameters between particles: static friction coefficient, rolling friction coefficient, restitution coefficient, and surface energy; (3) Powder-phase parameters: particle morphology, particle size, and powder fluidity. (2) JKR theory The Hertz-Mindlin with JKR contact model can form an adhesive force contact model between particles, mainly considering the influence of van der Waals forces in the contact area. In this model, the normal elastic contact force between particles is computed based on the JKR theory [38], including the overlap of particles δ, interaction parameters, and surface energy γ, which can be described as where E * is the equivalent Young's modulus, r * is the equivalent curvature radius, and α is the contact radius of particles. It is discovered that several parameters are involved in the aforementioned model, such as static friction coefficient, rolling friction coefficient, restitution coefficient, and surface energy. Therefore, it is essential to carry out the parameter calibration work through the material characterization.

Parameter Classification of the Discrete Element Model
The first step of the DEM simulation is to calibrate numerous parameters through comprehensive material characterization. Figure 3 shows the parameters involved in the DEM model, which can be divided into three categories: (1) Solid-phase parameters: solid phase density, elastic modulus, and Poisson's ratio; (2) Contact parameters between particles: static friction coefficient, rolling friction coefficient, restitution coefficient, and surface energy; (3) Powder-phase parameters: particle morphology, particle size, and powder fluidity. For composite powders, the contact forms between particles can be further divided into three patterns of "particle-particle", "particle-reinforced phase", and "reinforced phase-reinforced phase". Specifically, the particle size and the static friction coefficient between particles are directly measured. The restitution coefficient is obtained via the For composite powders, the contact forms between particles can be further divided into three patterns of "particle-particle", "particle-reinforced phase", and "reinforced phase-reinforced phase". Specifically, the particle size and the static friction coefficient between particles are directly measured. The restitution coefficient is obtained via the FEM simulation. The rolling friction coefficient and surface energy are achieved using the "PB-BBD-GA" parameter-design method. This method extends the GA on the basis of traditional PB and BBD methods, avoiding the uncertainty of traditional methods in predicting the optimal parameter combination. The "PB-BBD-GA" method can quickly obtain the functional relationship between significant parameters and target values and determine the unknown parameter combination more accurately and efficiently.

Analysis Step and Particle Modeling
The DEM simulation is proceeding based on the time step. In the numerical computation, a differential equation needs to be transformed into a difference equation for the iterative solution. The time step is the time difference between two adjacent iterations. The time step is usually set to 20-30% of Rayleigh time, which makes the simulation stable and timesaving. The Rayleigh time is determined as [39] where R min is the minimum radius of spherical members and ρ is the density. In this work, the Rayleigh time is 6.4175 × 10 −9 s and the DEM results are recorded every 0.001 s. The particle model is established based on its morphology in the DEM simulation, and the spherical particles are computationally efficient compared with the non-spherical particles. Importantly, the contact model for spherical particles should include rolling resistance to model the rotational motion of actual non-spherical particles. In contrast, the rolling resistance can be omitted when utilizing the non-spherical particles [40]. In this work, the spherical particles are modeled and the rolling resistance represented by the rolling friction coefficient is determined last. Meanwhile, a multi-sphere method, by which multiple spheres with the same or different sizes are clumped to approximate the irregularly shaped particle (such as CF in this work), is employed in the model [41], as shown in Figure 1c.

Material Characterization
This work selects the materials of PEEK (ZYPEEK_330PF, Jilin Joinature Polymer Co., Ltd., Changchun, China) and PEEK/CF (CF from NanJing WeiDa Composite Material Co. Ltd., Nanjing, China) powders, which possess the features of poor powder fluidity and high difficulty in powder deposition process, as the representative research object. The adopted powders are the PEEK, PEEK/CF_30wt%, and PEEK/CF_50wt%, which are featured by the weight fraction of CF. The inherent properties of the compositions of PEEK and CF are listed in Table 1, which includes important input factors of the subsequently introduced powder dynamics model. As shown in Figure 3, a large number of constants need to be determined, which are described in detail as follows.
(1) Angle of repose Powder fluidity is a comprehensive characterization that deeply impacts the powder deposition process. The AOR can reflect the powder fluidity in a simple and convenient fashion. In this paper, the AOR is measured according to GB11986-89. As shown in Figure 4, a funnel is fixed above the horizontal plate, through which the particle sample falls, and the AOR of the cone is measured. Both experiments and the DEM model follow this standard. Figure 3, a large number of constants need to be determined, which are described in detail as follows.

As shown in
(1) Angle of repose Powder fluidity is a comprehensive characterization that deeply impacts the powder deposition process. The AOR can reflect the powder fluidity in a simple and convenient fashion. In this paper, the AOR is measured according to GB11986-89. As shown in Figure  4, a funnel is fixed above the horizontal plate, through which the particle sample falls, and the AOR of the cone is measured. Both experiments and the DEM model follow this standard. (2) Particle morphology and particle size distribution Particle morphology and particle size distribution have important impacts on powder fluidity and computational efficiency. The scanning electron microscope (SEM, TM4000PLUS, HITACHI, Tokyo, Japan) equipment is utilized in the high vacuum circumstance and at a voltage of 5 kV to acquire the SEM images that characterized the morphological features of PEEK powders and PEEK/CF composite powders. The objective area of the SEM image is limited. Therefore, the Malvern laser particle sizer (Mastersizer 2000, the Malvern Panalytical Ltd, Marvin City, UK) is utilized to measure the accurate particle size distribution of PEEK and CF.
(3) Bulk density The measurement method for the bulk density of PEEK and PEEK/CF powders is shown in Figure 5, which refers to GB/T 16913. . Prepare a 100 mL graduated cylinder and place the funnel at a certain height directly above the cylinder. The powder falls freely from the funnel mouth and fills the cylinder. Thereafter, screed the surface of the graduated cylinder, place it on the electronic scale to measure its mass, and obtain the bulk density of the powders. (2) Particle morphology and particle size distribution Particle morphology and particle size distribution have important impacts on powder fluidity and computational efficiency. The scanning electron microscope (SEM, TM4000PLUS, HITACHI, Tokyo, Japan) equipment is utilized in the high vacuum circumstance and at a voltage of 5 kV to acquire the SEM images that characterized the morphological features of PEEK powders and PEEK/CF composite powders. The objective area of the SEM image is limited. Therefore, the Malvern laser particle sizer (Mastersizer 2000, the Malvern Panalytical Ltd, Marvin City, UK) is utilized to measure the accurate particle size distribution of PEEK and CF.
(3) Bulk density The measurement method for the bulk density of PEEK and PEEK/CF powders is shown in Figure 5, which refers to GB/T 16913. . Prepare a 100 mL graduated cylinder and place the funnel at a certain height directly above the cylinder. The powder falls freely from the funnel mouth and fills the cylinder. Thereafter, screed the surface of the graduated cylinder, place it on the electronic scale to measure its mass, and obtain the bulk density of the powders.

(4) Static friction coefficient
The inclined plane test can be used to measure the static friction coefficient. The friction coefficient between different materials can be achieved by replacing the material of the inclined plane or sample. At first, the inclined plane is placed horizontally, and the sample is placed at the fixed position of the inclined plane. Then, lift the slope slowly, and stop the inclined plane when the sample starts to slide. The position of the inclined plane at this moment is recorded through a camera, and the tilt angle is read. The static friction coefficient is the tangent value of the inclination angle. It is worth noting that when measuring the static friction coefficient between powders, it is necessary to evenly sprinkle some corresponding powder samples on the inclined surface of the same material, which makes the results more accurate. The restitution coefficient is the ratio of the separation velocity and the approaching velocity of the two objects along the normal direction of the contact before and after the collision, which is only related to the material of the collision objects [44]. Generally, it can be experimentally measured. Fix the plate, and a ball falls freely above the plate. The coefficient of restitution can be calculated by recording the fall height and bounce height of the ball (the restitution coefficient between different materials is achieved by changing the materials of the plate and the free-falling ball). However, measuring the restitution coefficient associated with pure CF is a struggle due to the difficulty in producing the products in pure CF in reality. Thus, it is appropriate to calculate the restitution coefficient using the simulation. The process of free-falling and bouncing for a small ball is achieved via FEM simulation (as shown in Figure 6) [45]. The restitution coefficient is computed as where v app is the speed of the ball approaching the plate, v off is the speed of the ball away from the plate, and e denotes the restitution coefficient. (4) Static friction coefficient The inclined plane test can be used to measure the static friction coeffi tion coefficient between different materials can be achieved by replacing t the inclined plane or sample. At first, the inclined plane is placed horizon sample is placed at the fixed position of the inclined plane. Then, lift the slo stop the inclined plane when the sample starts to slide. The position of the at this moment is recorded through a camera, and the tilt angle is read. The coefficient is the tangent value of the inclination angle. It is worth noting th  Experimentally measuring the rolling friction coefficient and surface energy is a struggle and can even be impossible to achieve. Thus, the DOE including the PB, BBD, and GA is employed to determine them.
PB design is a two-level (−1 and +1) experimental design method, which attempts to screen and determine the factors that significantly impact the results with the least number of trials. In this work, AOR is used as the response value to screen eight factors corresponding to the powders of PEEK/CF_30wt%. These factors with their initial parameter spaces are listed in Table 2. When determining the parameter space, the ratio of high and low levels of the two parameters is set to be identical in order to prevent significant fluctuation between the investigated factors [46]. The detailed design and results of the PB design are discussed in Section 4.2. The BBD design is utilized to further investigate the optimal points of the key factors based on the PB results. The computed results are analyzed afterward by the response surface method. The corresponding variables are coded as follows: where Xi is the coded value of the independent factor, x0 is the real value at the center point, xi is the real value, and ∆xi is the value of step-change.
The BBD model can be described by the quadratic equation: where Y is the predicted response value (AOR in this work), Xi and Xj are the coded independent variables, β0 is the offset term, βi is the linear term, βii is the squared term, and βij is the interaction term. GA is developed by imitating the mechanism of natural biological evolution. It is an efficient, parallel, global searching method, which can automatically acquire and Experimentally measuring the rolling friction coefficient and surface energy is a struggle and can even be impossible to achieve. Thus, the DOE including the PB, BBD, and GA is employed to determine them.
PB design is a two-level (−1 and +1) experimental design method, which attempts to screen and determine the factors that significantly impact the results with the least number of trials. In this work, AOR is used as the response value to screen eight factors corresponding to the powders of PEEK/CF_30wt%. These factors with their initial parameter spaces are listed in Table 2. When determining the parameter space, the ratio of high and low levels of the two parameters is set to be identical in order to prevent significant fluctuation between the investigated factors [46]. The detailed design and results of the PB design are discussed in Section 4.2. The BBD design is utilized to further investigate the optimal points of the key factors based on the PB results. The computed results are analyzed afterward by the response surface method. The corresponding variables are coded as follows: where X i is the coded value of the independent factor, x 0 is the real value at the center point, x i is the real value, and ∆x i is the value of step-change.
The BBD model can be described by the quadratic equation: where Y is the predicted response value (AOR in this work), X i and X j are the coded independent variables, β 0 is the offset term, β i is the linear term, β ii is the squared term, and β ij is the interaction term.
GA is developed by imitating the mechanism of natural biological evolution. It is an efficient, parallel, global searching method, which can automatically acquire and accumulate knowledge about the searching space, and the searching process is adaptively controlled to obtain the optimal solution. The entire GA is achieved based on the Matlab GA toolbox in this work.

Powder Deposition Process
The reliability of the complete model is supposed to be experimentally verified. The comprehensive indicator AOR can be used as the target, which reflects the overall performance of the studied powder. In addition, this work investigates the effect of CF on the surface properties of the deposited powder layer through powder deposition experiments (similar to the scraper form of the actual LPBF process) and corresponding simulation.
The surface roughness is the key factor in the powder deposition process. It has a great influence on the sintering process and the subsequent recoating process. Thus, it is necessary to investigate powder deposition through the powder dynamics model. By comparing with the corresponding experiment, the powder dynamics model can be further verified by this practical application.
The powder deposition is experimentally conducted based on a piece of slotted glass, which is shown in Figure 7a. The preparation steps are as follows: (1) prepare a piece of glass with a groove of 0.8-1 mm; (2) load enough powder manually in front of the groove; (3) use another piece of glass to scrape the loaded powders slowly from left to right with a tilt angle of around 45 • ; (4) the surface of the powders in the groove is smooth to the naked eye, as shown in Figure 7b; (5) the surface morphology of the target zones (Figure 7c) is measured using a microscope (J6000, Keyence, Shanghai, China). The observed area is located at the center of the powder layer.
In the DEM model, set the same conditions as the experiment. When the simulation of the powder deposition process is completed, obtain the surface particle contour curve at the center of the powder layer (as shown in Figure 7d). The fluctuation of the contour curve reflects the flatness of the powder surface.
accumulate knowledge about the searching space, and the searching process is adaptiv controlled to obtain the optimal solution. The entire GA is achieved based on the Mat GA toolbox in this work.

Powder Deposition Process
The reliability of the complete model is supposed to be experimentally verified. T comprehensive indicator AOR can be used as the target, which reflects the overall perf mance of the studied powder. In addition, this work investigates the effect of CF on surface properties of the deposited powder layer through powder deposition experime (similar to the scraper form of the actual LPBF process) and corresponding simulation The surface roughness is the key factor in the powder deposition process. It ha great influence on the sintering process and the subsequent recoating process. Thus, i necessary to investigate powder deposition through the powder dynamics model. comparing with the corresponding experiment, the powder dynamics model can be f ther verified by this practical application.
The powder deposition is experimentally conducted based on a piece of slotted gla which is shown in Figure 7a. The preparation steps are as follows: (1) prepare a piece glass with a groove of 0.8-1 mm; (2) load enough powder manually in front of the groo (3) use another piece of glass to scrape the loaded powders slowly from left to right w a tilt angle of around 45°; (4) the surface of the powders in the groove is smooth to naked eye, as shown in Figure 7b; (5) the surface morphology of the target zones (Figu 7c) is measured using a microscope (J6000, Keyence, Shanghai, China). The observed a is located at the center of the powder layer.
In the DEM model, set the same conditions as the experiment. When the simulat of the powder deposition process is completed, obtain the surface particle contour cur at the center of the powder layer (as shown in Figure 7d). The fluctuation of the conto curve reflects the flatness of the powder surface.

Results and Discussion
This section firstly describes the material properties of PEEK and PEEK/CF in de including particle morphology, particle size, and AOR, determines the contact parame between particles (such as rolling friction coefficient, restitution coefficient, and sur energy), and experimentally verifies the reliability of the determined parameters. ondly, the powder deposition process is investigated through the experiment and DEM simulation. The influence of CF particles on the surface performance of the pow deposition layer is analyzed.

Material Evaluation
The particle size distribution is obtained as input information in DEM, and it sig cantly influences computational efficiency. The SEM pictures of pure PEEK PEEK/CF_30wt% are given in Figure 8. It shows that the powder particles of PEEK almost ellipsoidal. The diameter of CF is uniform, but the length distribution is obvio dispersed. To specifically characterize the particle size information, the powders of P and CF are detected through the laser particle sizer, and the results are given in Tab and Figure 9.

Results and Discussion
This section firstly describes the material properties of PEEK and PEEK/CF in detail, including particle morphology, particle size, and AOR, determines the contact parameters between particles (such as rolling friction coefficient, restitution coefficient, and surface energy), and experimentally verifies the reliability of the determined parameters. Secondly, the powder deposition process is investigated through the experiment and the DEM simulation. The influence of CF particles on the surface performance of the powder deposition layer is analyzed.

Material Evaluation
The particle size distribution is obtained as input information in DEM, and it significantly influences computational efficiency. The SEM pictures of pure PEEK and PEEK/CF_30wt% are given in Figure 8. It shows that the powder particles of PEEK are almost ellipsoidal. The diameter of CF is uniform, but the length distribution is obviously dispersed. To specifically characterize the particle size information, the powders of PEEK and CF are detected through the laser particle sizer, and the results are given in Table 3 and Figure 9.

Results and Discussion
This section firstly describes the material properties of PEEK and PEEK/CF in detail, including particle morphology, particle size, and AOR, determines the contact parameters between particles (such as rolling friction coefficient, restitution coefficient, and surface energy), and experimentally verifies the reliability of the determined parameters. Secondly, the powder deposition process is investigated through the experiment and the DEM simulation. The influence of CF particles on the surface performance of the powder deposition layer is analyzed.

Material Evaluation
The particle size distribution is obtained as input information in DEM, and it significantly influences computational efficiency. The SEM pictures of pure PEEK and PEEK/CF_30wt% are given in Figure 8. It shows that the powder particles of PEEK are almost ellipsoidal. The diameter of CF is uniform, but the length distribution is obviously dispersed. To specifically characterize the particle size information, the powders of PEEK and CF are detected through the laser particle sizer, and the results are given in Table 3 and Figure 9.     It shows that the average particle size of PEEK is larger than that of CF. The majority of PEEK particles are in the range of 30-60 µm, while the lengths of numerous CF particles are scoped by 5-20 µm.
The restitution coefficient and the static friction coefficient are directly acquired according to the finite element simulation and inclined plane test, respectively. The corresponding results are summarized in Table 4. It appears that the PEEK-related static friction coefficients reach higher values compared with the CF-CF, which is caused by the irregular shapes and unsmooth surface property of the PEEK particles. The glass-related static friction coefficients are minimal due to the glossy surface. As for the restitution coefficient, the PEEK particles possess good elasticity, inducing the PEEK-related restitution coefficients to present larger values than the CF-related restitution coefficients. The bulk density results of PEEK and PEEK/CF powders are shown in Table 5. It is It shows that the average particle size of PEEK is larger than that of CF. The majority of PEEK particles are in the range of 30-60 µm, while the lengths of numerous CF particles are scoped by 5-20 µm.
The restitution coefficient and the static friction coefficient are directly acquired according to the finite element simulation and inclined plane test, respectively. The corresponding results are summarized in Table 4. It appears that the PEEK-related static friction coefficients reach higher values compared with the CF-CF, which is caused by the irregular shapes and unsmooth surface property of the PEEK particles. The glass-related static friction coefficients are minimal due to the glossy surface. As for the restitution coefficient, the PEEK particles possess good elasticity, inducing the PEEK-related restitution coefficients to present larger values than the CF-related restitution coefficients. The bulk density results of PEEK and PEEK/CF powders are shown in Table 5. It is observed that the bulk density of PEEK/CF is higher than that of the PEEK powder, and the bulk density of PEEK/CF increases with the increase in the CF mass fraction. Because the solid-phase density of CF is slightly higher than that of the PEEK, most CF particles have smaller sizes, filling the pores between PEEK particles.

Results of DOE
The rolling friction coefficient and surface energy cannot be directly measured through experimental methods, so the "PB-BBD-GA" method is used to acquire the optimal combination of rolling friction coefficient and surface energy. Therein, the PB design aims to determine the significant parameters of the related rolling friction coefficient and surface energy. The detailed design and response results of the PB test are listed in Table 6. The significance analysis of the model and corresponding factors is given in Table 7.

Rolling Friction Coefficient
Surface Energy Virtual Parameters  As shown in Table 7, both R2 and Adjusted-R2 are close to 1, which means that the model and the chosen response value of AOR are reasonable. The significance can be reflected by the p-value. The p-value less than 0.05 correlates with the significant factor. Therefore, the parameters of PEEK-PEEK rolling friction coefficient (X 0 ), PEEK-CF rolling friction coefficient (X 1 ), and PEEK-PEEK surface energy (X 5 ) have significant effects on the AOR. Among them, the influence of X 0 , X 1 , and X 5 are in descending order, which can be obtained through the F-value. This indicates that the effect of the rolling friction coefficient is greater than that of surface energy in this work. Since the ellipsoidal particles of PEEK are replaced by round particles, the rolling friction coefficient should be utilized to make up for the shape approximation in the powder dynamics simulation. Analogous treatment is also found in the literature of [47,48]. Regarding the parameter of rolling friction coefficient/surface energy, the general influence order of PEEK-PEEK, PEEK-CF, and CF-CF is in descending trend, which means that the interaction between CF particles is minimally effective in the composite powders PEEK/CF.
Since the significant parameters (X 0 , X 1 , X 5 ) are determined by the PB test, the remaining insignificant parameters can be reasonably fixed before the BBD design. The parameters of X 2 , X 3 , and X 4 are fixed as 0.01, X 6 and X 7 are set to 0.001. Sequentially, the detailed BBD design (correlated with the parameters of X 0 , X 1 , and X 5 ) and response-value of AOR are listed in Table 8. The variance analysis of the whole model and the corresponding factors are given in the following Table 9. Herein, the insignificant quadratic and cross-terms are eliminated for the convenience of analysis and subsequent calculation. The p-value of the model is extremely small, which shows that the entire model is reliable and the AOR can be reasonably predicted. Interestingly, the significance of cross-terms is greater than that of the quadratic term (eliminated here) in this numerical example. Based on the results of the BBD design, the quadratic regression model for the response value (AOR) and the investigated variables (X 0 , X 1 , X 5 ) can be established as: Therein, Equation (25) is in terms of the actual factors and can be utilized to predict the AOR. Equation (26) is in terms of the coded factors, which is necessary for the subsequent GA. In Equation (26), the term A corresponds to X 0 , B corresponds to X 1 , and C corresponds to X 5 . The upper (−1) and lower limits (1) of the value range of term A correlated with the upper (0.01) and lower limits (0.2) of factor X 0 , which is also correspondingly applied for the terms B and C.
When the predictive equation of AOR is acquired based on the BBD design, the corresponding optimal parameters of rolling friction coefficient and surface energy can be achieved by making the predicted AOR equal to the experimental value (34.23 • for the PEEK/CF_30wt%). Thus, the objective equation for the GA is expressed as Equation (27) is input into the MATLAB GA toolbox with the value range of [−1, 1] for all corresponding variables. The results are given in Table 10. So far, all the essential factors have been determined for the powder dynamics model.

Experimental Verification
The calibrated DEM model for the PEEK/CF composite powders should be experimentally verified. As a comprehensive index of powder properties, AOR can be utilized to conduct the experimental verification work. The AOR experiment is carried out based on the standard GB11986-89. The experimental AOR results of PEEK, PEEK/CF_30wt%, and PEEK/CF_50wt% are shown in Figure 10a, with the detailed data listed in Table 11. The simulation results are shown in Figure 10b, and the detailed data are given in Table 12. It appears that the experimental and simulation results show a maximal deviation of 4.89%, which can indicate the high reliability and effectiveness of the developed powder dynamics model. Therein, Equation (25) is in terms of the actual factors and can be utilized to pred the AOR. Equation (26) is in terms of the coded factors, which is necessary for the sub quent GA. In Equation (26), the term A corresponds to X0, B corresponds to X1, and corresponds to X5. The upper (−1) and lower limits (1) of the value range of term A cor lated with the upper (0.01) and lower limits (0.2) of factor X0, which is also correspon ingly applied for the terms B and C.
When the predictive equation of AOR is acquired based on the BBD design, the c responding optimal parameters of rolling friction coefficient and surface energy can achieved by making the predicted AOR equal to the experimental value (34.23° for t PEEK/CF_30wt%). Thus, the objective equation for the GA is expressed as (27) is input into the MATLAB GA toolbox with the value range of [−1, for all corresponding variables. The results are given in Table 10. So far, all the essent factors have been determined for the powder dynamics model.

Experimental Verification
The calibrated DEM model for the PEEK/CF composite powders should be expe mentally verified. As a comprehensive index of powder properties, AOR can be utiliz to conduct the experimental verification work. The AOR experiment is carried out bas on the standard GB11986-89. The experimental AOR results of PEEK, PEEK/CF_30wt and PEEK/CF_50wt% are shown in Figure 10a, with the detailed data listed in Table  The simulation results are shown in Figure 10b, and the detailed data are given in Ta 12. It appears that the experimental and simulation results show a maximal deviation 4.89%, which can indicate the high reliability and effectiveness of the developed powd dynamics model.  It also reveals that the pure PEEK possesses the largest AOR, which implies the m imum of powder fluidity. With the addition of CF, the AOR of the composite powd decreases and the powder fluidity rises. Importantly, the AOR decreases with the increa in the weight fraction of CF, exhibiting improved powder fluidity. Since the powder fl idity is greatly affected by the cohesive effect between particles of PEEK-PEEK, PEEK-C and CF-CF. The influence of PEEK-PEEK, PEEK-CF, and CF-CF interfacial reactions the composite AOR is in descending order (discussed in Section 4.2). The PEEK-relat cohesion plays a dominant role in the PEEK/CF powder packing. Thus, the increased v ume fraction of CF leads to a decrease in total cohesive interaction between particles, sulting in improved powder fluidity.

Evaluation of the Powder Deposition
The investigation of the powder deposition process is necessary to achieve the e pected behaviors of the laser-sintered parts in PEEK/CF composites. The surface prop ties and recoating difficulty of powder layers can be evaluated in a qualitative fashi using the powder deposition process. The fluctuation of the powder surface is observ by an optical microscope (VHX-6000, KEYENCE (CHINA) CO., LTD., Shanghai, Chin Figure 11 shows the micro-pictures of the surface topography and undulation of PEE PEEK/CF_30wt%, and PEEK/CF_50wt%. It appears that the image quality of the surfa   It also reveals that the pure PEEK possesses the largest AOR, which implies the minimum of powder fluidity. With the addition of CF, the AOR of the composite powder decreases and the powder fluidity rises. Importantly, the AOR decreases with the increase in the weight fraction of CF, exhibiting improved powder fluidity. Since the powder fluidity is greatly affected by the cohesive effect between particles of PEEK-PEEK, PEEK-CF, and CF-CF. The influence of PEEK-PEEK, PEEK-CF, and CF-CF interfacial reactions on the composite AOR is in descending order (discussed in Section 4.2). The PEEK-related cohesion plays a dominant role in the PEEK/CF powder packing. Thus, the increased volume fraction of CF leads to a decrease in total cohesive interaction between particles, resulting in improved powder fluidity.

Evaluation of the Powder Deposition
The investigation of the powder deposition process is necessary to achieve the expected behaviors of the laser-sintered parts in PEEK/CF composites. The surface properties and recoating difficulty of powder layers can be evaluated in a qualitative fashion using the powder deposition process. The fluctuation of the powder surface is observed by an optical microscope (VHX-6000, KEYENCE (CHINA) CO., LTD., Shanghai, China). Figure 11 shows the micro-pictures of the surface topography and undulation of PEEK, PEEK/CF_30wt%, and PEEK/CF_50wt%. It appears that the image quality of the surface morphology of PEEK is worse than the other two composite powders. This indicates a violent surface fluctuation of the PEEK powders, which is caused by the inhomogeneous particle distribution-induced porosity. With the addition of the CF, the surface quality is improved and also meliorates with the increase in the weight fraction of CF, which is caused by the particle distribution of the adopted CF. As shown in Figure 9b, numerous CF particles exhibit the size range of 5-20 µm. These small and uneven CF particles exactly fill the inter-particle cavity among PEEK powders, and the surface of its powder bed is observed smoothly.  Figure 12 shows the cross-sections of the simulated powder deposition results. It also shows that drastic undulation occurs on the PEEK surface. The surface fluctuation of PEEK/CF_30wt% is slightly mitigated and improved a lot for the PEEK/CF_50wt%, which is consistent with the experimental results. The standard deviation of the surface contour curve of the powder deposition layer is used to quantitatively characterize the fluctuation of the layer surface, which is similar to the RMS roughness. The fluctuation of the layer surface decreases with the decrease in the standard deviation, and the corresponding powder layer becomes smoother.  Figure 12 shows the cross-sections of the simulated powder deposition results. It also shows that drastic undulation occurs on the PEEK surface. The surface fluctuation of PEEK/CF_30wt% is slightly mitigated and improved a lot for the PEEK/CF_50wt%, which is consistent with the experimental results. The standard deviation of the surface contour curve of the powder deposition layer is used to quantitatively characterize the fluctuation of the layer surface, which is similar to the RMS roughness. The fluctuation of the layer surface decreases with the decrease in the standard deviation, and the corresponding powder layer becomes smoother. The experimental and simulation results of the standard deviation of the surface contour curve are shown in Figure 13 and Table 13. According to the experimental results, the surface smoothness of PEEK/CF_30wt% is increased by 31.3% compared to the PEEK powder, while the surface smoothness of PEEK/CF_50wt% powder is improved by 55.2% compared to PEEK powder. According to the simulation results, compared to the PEEK powder, the surface smoothness of the PEEK/CF_30wt% is increased by 27.4%, and the PEEK/CF_50wt% powder is improved by 50.2%. The similar experimental and simulation results further verify the reliability of the DEM model.   The experimental and simulation results of the standard deviation of the surface contour curve are shown in Figure 13 and Table 13. According to the experimental results, the surface smoothness of PEEK/CF_30wt% is increased by 31.3% compared to the PEEK powder, while the surface smoothness of PEEK/CF_50wt% powder is improved by 55.2% compared to PEEK powder. According to the simulation results, compared to the PEEK powder, the surface smoothness of the PEEK/CF_30wt% is increased by 27.4%, and the PEEK/CF_50wt% powder is improved by 50.2%. The similar experimental and simulation results further verify the reliability of the DEM model. The experimental and simulation results of the standard deviation of the surface contour curve are shown in Figure 13 and Table 13. According to the experimental results, the surface smoothness of PEEK/CF_30wt% is increased by 31.3% compared to the PEEK powder, while the surface smoothness of PEEK/CF_50wt% powder is improved by 55.2% compared to PEEK powder. According to the simulation results, compared to the PEEK powder, the surface smoothness of the PEEK/CF_30wt% is increased by 27.4%, and the PEEK/CF_50wt% powder is improved by 50.2%. The similar experimental and simulation results further verify the reliability of the DEM model.    In summary, the method of systematic material evaluation possesses universal applicability for different types of composite powders, even multi-component composites.