Numerical investigation of fiber orientations and homogeneity in powder bed fusion of fiber/polymer composites

ABSTRACT The packing characteristics of fiber/polymer powder in powder bed fusion additive manufacturing exhibit a high correlation with the mechanical behaviours of printed composite parts such as homogeneity and anisotropy. A discrete element model has been developed to investigate the packing characteristics of glass fiber/polyamide 12 (PA12) powder, which include fiber orientations, fiber homogeneity, and packing density. The predicted probability distributions of fiber orientations in the powder bed are comparable with those measured in glass fiber–reinforced PA12 composites printed via multi jet fusion. Three types of fibers with different length distributions are adopted to study the effects of the fiber length distribution on their packing characteristics. The simulation results reveal that a large average fiber length is beneficial to fiber alignment in the powder spreading direction but lowers the fiber homogeneity and packing density of the powder bed. Furthermore, varying the fiber length can provide an effective way to regulate fiber orientations in the powder packing process, which would help achieve satisfactory anisotropic mechanical properties for composite parts.


Introduction
Fiber-reinforced polymer composites have found wide applications in the automotive, aerospace, and biomedical industries because of their high corrosion resistance and strength-to-weight ratios.To meet the growing need for parts with complex geometries and customized functionalities, the additive manufacturing (AM) of polymeric materials has garnered tremendous interest among research scholars.AM technologies involving polymeric materials can be categorized into several branches including powder bed fusion (PBF) (Jing et al. 2017;Arai et al. 2018), material extrusion (Goh et al. 2018), vat photo polymerization (Yang et al. 2020), and sheet lamination (Klosterman et al. 1998).Compared with other AM technologies involving polymeric materials, PBF, including selective laser sintering and multijet fusion (MJF), possesses advantages such as being applicable to a wide variety of materials, not requiring support structures, and the capability of printing parts with intricate internal geometries (Yuan et al. 2019).To broaden the range of applications of fiberreinforced polymer composites processed by PBF, several existing concerns must be addressed, such as fiber orientations, fiber homogeneity, and the minimization of porosity in printed composite parts.
The fiber orientations of fiber-reinforced polymer composites have an essential impact on the anisotropic properties of the composite parts fabricated by both traditional manufacturing methods and novel AM technologies.Significant improvements in the mechanical strength and stiffness of composite parts have been observed in the dominant direction of fiber alignment (Parandoush and Lin 2017;Mortazavian and Fatemi 2015).In conventional manufacturing processes, a variety of approaches, including shear flow in narrow channels (Trebbin et al. 2013), magnetic fields (Kimura, Umehara, and Kimura 2010;Song et al. 2014), and electric fields (Yang et al. 2017), have been adopted to regulate the fiber orientations in polymer composites.AM technologies can provide effective control of the fiber orientations in successive layers, thereby tailoring the mechanical properties of the printed parts for various applications (Tekinalp et al. 2014;Yunus et al. 2016).Fiber-reinforced polymer composites printed via PBF also exhibit anisotropic mechanical properties due to fiber alignment along the powder spreading direction (Türk et al. 2017;Goodridge et al. 2011;Yan et al. 2011;Badini et al. 2020).The PBF process can be divided into three sub-steps: powder packing, energy input, and melting and coalescence (Chatham, Long, and Williams 2019).Although the melting and coalescence of polymer powder can influence the fiber orientation distribution due to the polymer melt flow and the shrinkage of the powder bed, powder packing is one of the main processes to determine the fiber orientations in printed composites.In the powder packing process, fibers tend to align along the spreading direction of the blade or roller spreader because of the particle shear flow developed in the powder pile (Nan and Ghadiri 2019;Chen et al. 2021).
Fiber homogeneity is crucial for enabling homogeneous mechanical properties in composite parts.Since mechanical failures are typically initiated from the weakest regions in the parts, inhomogeneous mechanical properties of such parts would lower their overall mechanical strength.An inhomogeneous fiber distribution among different composite parts may also result in a significant deviation between their mechanical strength (Tekinalp et al. 2014).To reduce fiber inhomogeneity, a few studies have focused on various methods for preparing composite powders (Goodridge et al. 2011;Yan et al. 2011;Christ et al. 2015).Owing to the distinct difference in the shapes, sizes, and material properties between polymer particles and fibers, particle segregation may occur during the powder packing process, thereby inducing an inhomogeneous fiber distribution within the powder bed.Presently, the research on fiber homogeneity in fiber/polymer powder beds is lacking.
In terms of the cost and challenges of studying the flow dynamics of powder particles in experiments, numerical modelling provides an effective route to understanding the underlying physical mechanisms during the particle flow and optimising the processing parameters to improve the packing quality of the powder bed.The discrete element method (DEM), which was proposed by Cundall and Strack for analyzing rock mechanics (Cundall and Strack 1979), has been extensively applied to simulate powder packing processes in AM.The DEM is capable of predicting the packing characteristics of polymer and metal powders such as packing density and surface roughness (Parteli and Pöschel 2016;Haeri et al. 2017;Nan et al. 2018;He, Hassanpour, and Bayly 2020).It has also been employed to study fiber alignment in the packing process of fiber/polymer powder (Chen et al. 2021).To print high-quality customized composite parts, a systematic study of the packing characteristics of the powder bed, including packing density, fiber orientations, and fiber homogeneity, is required.
In this work, a DEM model has been employed to investigate the packing characteristics of glass fiber/ polyamide 12 (PA12) powder in PBF, including fiber orientations, fiber homogeneity, and packing density.The predicted probability distributions of the fiber orientations in the powder bed are compared with those measured in glass fiber-reinforced PA12 composites printed via MJF.The effects of the average fiber length on the packing characteristics of the glass fiber/PA12 powder are investigated.Further studies of the effects of the fiber size on fiber orientations are performed, which provide guidance for tailoring the anisotropic mechanical properties of printed composites.

DEM model for glass fiber/polymer powder
A particle-based DEM model is developed to simulate the packing process of glass fiber/PA12 powder in PBF.The model is built on the open-source DEM platform YADE, which is written in the C++ programming language and allows for the flexible implementation of new algorithms and models (Šmilauer et al. 2015).Herein, the particle contact model based on the Hertz-Mindlin theory is presented.The features of glass fibers and PA12 powder particles including the size distributions and shapes are represented methodically in the model.A detailed mathematical model for describing multi-sphere PA12 particles and sphero-cylindrical glass fibers can be found in our previous work (Tan et al. 2021).

Discrete element method
In the DEM, the translation and rotation of a particle i are governed by Newton's second law (Cundall and Strack 1979;Tan et al. 2021), which are described by where m i is the mass of the particle, v i is its velocity, t is the time, F n is the normal component of the contact force, F t is the tangential component of the contact force, g is gravitational acceleration, F v is the van der Waals force of the particle, I i is its moment of inertia, ω i is its angular velocity, and M is its total moment.The contact force is described by a spring-dashpot model based on the Hertz-Mindlin theory (Tsuji, Tanaka, and Ishida 1992;Thornton, Cummins, and Cleary 2013).Meanwhile, the van der Waals force is calculated according to the Derjaguin-Muller-Toporov model (Derjaguin, Muller, and Toporov 1975).The contact force and van der Waals force between two interacting particles i and j are expressed as where k n and k t are the normal and tangential stiffness coefficients, respectively; u n and u t are the normal and tangential displacement, respectively; η is the damping coefficient; v r n and v r t are the normal and tangential components of the relative velocity, respectively; n and n t are unit vectors in the normal and tangential directions, respectively; µ is the coefficient of friction; R* is the effective radius; γ is the surface energy density.

Glass fiber/PA12 powder
Polymer particles typically possess complex geometric shapes, while glass fibers display rod-like shapes, as shown in Figure 1(a).In this model, PA12 particles and glass fibers are approximated using multi-sphere particles and sphero-cylinders, respectively, as shown in Figure 1(b).Figure 1(c) presents several representative PA12 shapes captured by two-dimensional scanning electron microscope (SEM) images and their representative multi-sphere particles used in the simulation.Each multi-sphere particle is confined within a sphere with a minimum diameter D, which follows the particle size distribution of PA12 powder obtained from experiments (Tan et al. 2021).Notably, the PA12 powder particles display a median size of 52.6 μm.
Figure 2 shows the optical microscopic images and length distributions of three types of glass fibers used in the simulation.The glass fibers display a constant radius of 8 µm, and their lengths are analyzed using the ImageJ software.To perform accurate measurements, the three types of glass fibers were spread uniformly on the platform of the microscope to minimize overlaps, as shown in Figure 2. The average lengths L a of the three types of glass fibers are measured to be 124, 227, and 267 µm.
The material properties of glass fibers and PA12 powder particles in the simulation are listed in Table 1.For contact interactions between a glass fiber i and a PA12 particle j, the effective parameters, i.e. the coefficient of restitution λ e , the coefficient of friction μ e , and the surface energy density γ e , are expressed as To expedite the computational process, the Young's moduli of PA12 powder and glass fibers are set to be 500 times smaller than their actual values, and accordingly, the surface energy density is also scaled down by the same amount.A computational time step of 10 −8 s, which is determined by the criteria in the reported work (Tan et al. 2021), is used for the simulations.

Results and discussion
In the simulation, a counter-clockwise rotating roller is used to spread glass fiber/PA12 powder particles on a substrate to form a powder bed, as shown in Figure 3.To represent the previously printed layer, the substrate is composed of clumped spherical polymer particles.The domain of interest is constrained in a periodic cuboid whose width and length are equal to those of the substrate.The roller with a radius R = 10 mm has a translational velocity v and an angular velocity ω = v/R.The packing parameters used for the simulations are listed in Table 2.The fiber weight fraction in the initial powder pile is 20%.To eliminate the boundary effect, a wider powder bed is adopted for longer fibers.The width W of the powder beds for the packing of the three types of glass fiber/ polymer powder with L a = 124, 227, and 267 µm in the simulations is 1.0, 1.7, and 1.7 mm, respectively.

Fiber orientations and homogeneity
A snapshot of the packing process of fiber/polymer powder particles is shown in Figure 4.The powder layer thickness and the spreading velocity are 180 µm and 0.1 m/s, respectively.The x-z cross-section of the powder pile is illustrated in Figure 4(a).Notably, the particles in contact with the roller move at higher speeds than particles near the substrate, forming a velocity gradient between the roller and the substrate.According to the velocity vector in Figure 4(b), the powder particles tend to be dragged by the roller and predominantly move along the roller spreading direction, i.e. the x-direction.The x-component of the average particle velocity (as a function of z) at different regions (A1-A5) of the powder pile is shown in Figure 4(c), and it can be observed that the xcomponent of the average particle velocity increases along the building direction of the powder layers, i.e. the z-direction.The existence of a gradient of the x-component of the average particle velocity indicates the formation of a particle shear flow between the roller and the substrate, which results in the alignment of the fibers along the powder spreading direction (Guo et al. 2012), as shown by the schematic in Figure 4(d).
A glass fiber/polymer powder bed is displayed in Figure 5(a) (the PA12 powder particles and the glass fibers are delineated in blue and gray, respectively).After the polymer powder has been removed, the distribution of the glass fibers is presented in Figure 5(b).A schematic describing the fiber orientation angles, α and β, in a three-dimensional (3D) space is shown in Figure 5(c), and the probability distribution of the fiber orientations with respect to the angles α and β is shown in Figure 5(d).Two peaks can be observed at the points (α, β) = (90°, 0°) and (−90°, 0°), which indicates that the fibers tend to lie within the powder-spreading plane and align along the spreading direction.The probability distribution of the fiber orientations in the powder bed in terms of the angles α or β is shown in Figure 5(e).The probability values for −90°< α < −80°and 80°< α < 90°when −90°< β < 90°are 25.1% and 26.1%, respectively, while few fibers are oriented with the angle α being between −20°a nd 20°.Meanwhile, according to the fiber orientation probability distribution with respect to the angle β, the probability values for |β| < 10°and |β| < 30°when −90°< α < 90°are 26.3% and 61.0%, respectively.
The distribution of the fiber weight fraction w in different regions of the powder bed (Figure 5(a)) along the powder spreading direction is shown in Figure 5(f).The powder bed is divided into eight regions and the length of each region Δx is 0.5 mm.The initial fiber weight fraction w 0 is 20%.The deviation of the fiber weight fraction from the initial value can be attributed to particle segregation during the powder packing process.

Effects of the fiber length on packing characteristics
Particle features, such as their shapes and sizes, are of vital importance in determining the packing characteristics of the powder bed.Glass fibers commonly possess cylindrical shapes characterised by their radius and length.An insightful understanding of the effects of the particle feature on the packing characteristics of fiber/polymer powder is essential to achieving a high-quality powder bed.Hence, the effects of the fiber length on the fiber orientations, fiber homogeneity, and packing density of the powder bed are investigated.The fiber orientations are assessed by three orientation parameters, which represent the preference of fiber alignment along the three orthogonal directions (x, y, and z) in the 3D space (Nguyen Thi et al. 2015): where N is the number of fibers.Fiber homogeneity is evaluated based on the normalised deviation of the fiber weight fraction S: where w i is the fiber weight fraction in a subdivided region of the powder bed, w 0 is the initial fiber weight fraction before the powder packing process, and n is the number of subdivided regions in the powder bed.
A small normalised deviation S indicates good fiber homogeneity in the powder bed.Three types of glass fiber/polymer powder with different average fiber lengths (see Figure 2) are used in the simulation, and their respective powder beds are shown in Figure 6.The fiber weight fraction, powder layer thickness, and spreading velocity are 20%, 180 µm, and 0.1 m/s, respectively.As shown in the regions delineated by the red dashed ellipses in Figure 6(b, c), more small voids occur in the powder beds with larger average fiber lengths (L a = 227 and 267 µm).Furthermore, inhomogeneous fiber distributions can be observed in all the powder beds.
Figure 7(a) shows orientation parameters for the powder beds with the three types of fibers.The orientation parameter a x is much larger than the parameters a y and a z , indicating a preference for fiber alignment along the powder spreading direction.The increasing parameter a x with the increase in the average fiber length means that long fibers are useful for aligning fibers along the powder spreading direction.The normalised deviation of the fiber weight fraction S and the packing density w in powder beds with different average fiber lengths are shown in Figure 7(b).An increase in the average fiber length leads to a decrease in packing density, and the normalised deviation of the fiber weight fraction in the powder beds with L a = 227 µm and L a = 267 µm is larger than that in the powder bed with L a = 124 µm.This phenomenon can be attributed to the poor flow ability of the glass fiber/PA12 powder with longer fibers, which reduces particle flow stability during their packing process.Additionally, long fibers display a tendency to tangle together, resulting in void formation and Optical microscopic images of cross-sections (x-y, x-z, and y-z planes) of the glass fiber-reinforced PA12 composites printed by the MJF test bed from HP Inc. are presented in Figure 8.A counter-rotating roller packing system for the glass fiber/PA12 powder was employed in the MJF process, in which the spreading velocity v and the fiber weight fraction w 0 were 0.2 m/s and 20%, respectively.The substrate was lowered down by a distance of 90 µm for each layer.Notably, large pores are discovered in the printed composites with average fiber lengths L a = 227 and 267 µm, which can be attributed to the lower packing density of their respective powder beds (refer to the predicted packing density in 7(b)).Furthermore, the inhomogeneous fiber distribution in these powder beds may lead to the formation of more pores in the printed composites.In the regions containing densely distributed fibers (refer to Figure 6), pores are more likely to be formed during the powder melting and coalescence processes as a result of insufficient polymer melt filling the voids between the fibers.
A comparison of the probability distributions of the fiber orientation angle β obtained from the simulated powder beds and actual MJF-printed composite parts is shown in Figure 9.The ImageJ software was used to conduct a statistical analysis of the fiber orientation  angle in the x-y cross-section of the composites.Due to the symmetry displayed by the fiber orientation angle β with respect to the spreading direction, the probability distribution of |β| is considered.The probability values of |β| < 10°in the printed composite parts in which L a = 124 µm, 227 µm, and 267 µm were measured to be 23.2%, 25.6%, and 26.7%, respectively.Typically, powder beds undergo shrinkage during the printing process as a result of their porous nature, and the actual thickness of a powder layer is proximately two times larger than the lowering distance (90 µm) of the substrate at each layer (Chen et al. 2020;Bidare et al. 2017).Hence, in the simulation, a powder layer thickness of 200 µm is selected, while the spreading velocity and the fiber weight fraction are identical to the values adopted in the experiment, i.e. v = 0.2 m/s and w 0 = 20%.The probability values for |β| < 10°when L a = 124 µm, 227 µm, and 267 µm are predicted to be 25.6%, 26.0%, and 29.1%, respectively.Notably, both experimental and simulation results indicate the tendency for fibers to align along the spreading direction.The predicted fiber orientation probability is slightly higher than the experimental result, which can be attributed to the powder melting and coalescence processes that lead to the reorientation of the fibers.3.3.Tailoring fiber orientations by varying fiber sizes Fiber-reinforced polymer composites printed via PBF usually possess the highest mechanical strength along the powder spreading direction due to favourable fiber alignment, while their mechanical strength in the building direction is the lowest due to undesirable fiber alignment and low interlayer strength.In many applications, the composite products subjected to loading in different directions are required to display isotropic mechanical strength, and an enhancement of their mechanical strength along the building direction is necessary.In Section 3.2, short fibers (L a = 124 µm) are found to display superior fiber alignment along the building direction (the z-direction) to the other two types of fibers (L a = 227 and 267 µm).Hence, changing the fiber length is a potential strategy to regulate fiber orientations along different directions, which allows the anisotropic mechanical properties of fiber-reinforced polymer composites to be tailored for various applications.
In this study, the fiber length is further decreased and the effects of the fiber length on the fiber orientations in the powder bed are investigated.Glass fibers of uniform lengths are used in the simulation with the powder spreading velocity, the powder layer thickness, and the fiber weight fraction being 0.1 m/s, 200 μm, and 20%, respectively.Powder beds corresponding to different fiber lengths (L a = 24, 32, 64, and 100 µm) are shown in Figure 10.The probability distributions with respect to the fiber orientation angle β and the orientation parameters (a x , a y , and a z ) of the powder beds with different fiber lengths are shown in Figure 11.In Figure 11(a), when the fiber length is reduced, the probability of the fibers being aligned along a direction within −40°< β < 40°is shown to decrease, which indicates that shorter fibers exhibit a more uniform distribution of the fiber orientations.
Figure 11(b) shows that a decrease in the fiber length leads to a decrease in the orientation parameter a x and an increase in the two orientation parameters a y and a z .Notably, the differences among the three orientation parameters are found to gradually decrease as the fiber length is reduced.As the fiber length is reduced from 100 µm to 24 µm, the orientation parameters a y and a z increase from 0.23 to 0.31 (increasing 34.7%) and from 0.11 to 0.29 (increasing 163.6%), respectively.In particular, the significant increase in fiber alignment along the z-direction (building direction) may enhance the mechanical strength of the composites along that direction.
The effects of the fiber radius r f on fiber orientations are also investigated, as shown in Figure 12.When the fiber radius is reduced from 8 μm to 4 μm, no significant changes in orientation parameters are observed, in contrast to varying the fiber length.In particular, the orientation parameter a y remained nearly constant when the fiber radius is reduced at different fiber lengths.Therefore, varying the fiber radius may not be an effective method to manipulate the orientation behaviour of fibers during the powder packing process.

Conclusions
A DEM model has been developed to study the packing characteristics of glass fiber/PA12 powder in PBF.The numerical model is validated through comparing its predicted fiber orientation probability distributions with experimental results.The effects of the average fiber length on the packing characteristics of glass fiber/PA12 powder are investigated.The simulation results suggest that a large average fiber length enhances fiber alignment along the powder spreading direction but lowers fiber homogeneity and the packing density of the powder bed.Varying the fiber length could provide plausible procedures to manipulate fiber orientations, thereby effectively tailoring the anisotropic mechanical properties of printed fiberreinforced polymer composites.well as cash and in-kind contribution from the industry partner, HP Inc.   Kun Zhou is a Professor in the School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore.

Figure 1 .
Figure 1.(a) SEM image of actual glass fiber/PA12 powder particles, (b) glass fiber/PA12 powder particles in the simulation, and (c) a sphero-cylinder for representing a glass fiber of length L and multi-sphere particles for representing PA12 powder particles confined within a sphere of diameter D.

Figure 2 .
Figure 2. Optical microscopic images (left) and length distributions (right) for three types of glass fibers with different average lengths: (a) L a = 124 µm, (b) L a = 227 µm, and (c) L a = 267 µm.

Figure 3 .
Figure 3. Roller-type packing system consisting of a counter-clockwise rotating roller and a substrate.

Figure 5 .
Figure 5. (a) Glass fiber/PA12 powder bed and (b) glass fiber powder bed after the PA12 powder has been removed; (c) the schematic describing fiber orientation angles in a 3D space; (d) the probability distribution of the fiber orientations in the powder bed as a multivariable function of α and β; (e) the probability distribution of the fiber orientations with respect to the angle α or β; (f) the distribution of the fiber weight fraction in different regions of the powder bed along the spreading direction.

Figure 6 .
Figure 6.Powder beds with different average fiber lengths: (a) L a = 124 µm, (b) L a = 227 µm, and (c) L a = 267 µm.The red dashed ellipses indicate regions in the powder beds where voids exist.The yellow dashed ellipses and yellow solid ellipses indicate dense and sparse fiber distributions, respectively.

Figure 7 .
Figure 7. (a) Orientation parameters and (b) the normalised deviation of the fiber weight fraction and the packing density in the powder bed with different average fiber lengths.

Figure 9 .
Figure 9.Comparison of the probability distributions with respect to the fiber orientation angle β for different average fiber lengths in simulations and experiments.

Pengfei
Tan is a Research Fellow in HP-NTU Digital Manufacturing Corporate Lab, School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore.Xiaojiang Liu is a Research Fellow in HP-NTU Digital Manufacturing Corporate Lab, School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore.

Figure 12 .
Figure 12.Orientation parameters versus the fiber length at which the fiber radius r f = 4 and 8 μm.

Figure 11 .
Figure 11.(a) Probability distributions with respect to the fiber orientation angle β in powder beds with different fiber lengths and (b) the values of orientation parameters at different fiber lengths.