Mesoscopic-Scale Numerical Investigation Including the Influence of Process Parameters on LPBF Multi-Layer Multi-Path Formation

As a typical laser additive manufacturing technology, laser powder bed fusion (LPBF) has achieved demonstration applications in aerospace, biomedical and other fields. However, how to select process parameters quickly and reasonably is still the main concern of LPBF production. In order to quantitatively analyze the influence of different process parameters (laser power, scanning speed, hatch space and layer thickness) on the LPBF process, the multi-layer and multi-path forming process of LPBF was predicted based on the open-source discrete element method framework Yade and the open-source finite volume method framework OpenFOAM. Based on the design of experiments method, a four-factor three-level orthogonal test scheme was designed, and the porosity and surface roughness data of each calculation scheme were extracted. By analyzing the orthogonal test data, it was found that as the laser power increased, the porosity decreased, and as the scanning speed, hatch space, and layer thickness increased, the porosity increased. In addition, the influence of laser power and scanning speed on surface roughness showed a trend of decreasing first and then increasing, while the influence of scanning distance and layer thickness on surface roughness showed a monotonous increasing trend. The order of the influence of each process parameter on porosity was: scanning speed > laying thickness > laser power > hatch space, and the order of the influence of each process parameter on surface roughness was: hatch space > layer thickness > laser power > scanning speed. So the porosity of the part is most sensitive to scanning speed, and the surface roughness is the most sensitive to hatch space. The above conclusions are expected to provide process control basis for actual LPBF production of the 316L stainless steel alloy.

reasonably select process parameters so that the porosity [3] and surface roughness [4] are at a low level. Therefore, in-depth study of the in uence of process parameters on the LPBF process plays an important role in accelerating the industrial application of LPBF technology.
The core of the LPBF process is that the metal particles gradually melt after being heated by the laser, and then solidify to form a solidi ed track. This process is a typical high temperature and high transient physical process, where the highest temperature can reach 4000 • C [24], and the characteristic time of molten pool ow is on the order of microseconds [25]. Although the current experimental methods such as X-ray in-situ imaging technology and high-speed cameras can observe the LPBF process in real time [26], the experimental methods are dif cult to accurately obtain the three-dimensional information of the LPBF process, such as the evolution of the three-dimensional shape of the molten pool. As a typical quantitative research method, numerical simulation has become a powerful auxiliary method for analyzing LPBF process [27][28][29]. At present, the simulation research for the LPBF process is mainly based on the following three scales: (a) macroscopic-scale, that is, by establishing a simpli ed model, the transition between the three phases "powder bed-metal liquid-phase-metal solid-phase" is equivalently handled to obtain the temperature eld and stress eld distribution [30][31][32][33][34][35]; (b) Mesoscopic-scale, that is, the heating effect of the laser on the metal particles is directly calculated to describe the dynamic behavior of the mesoscopic molten pool [36][37][38][39][40][41]; (c) Microscopic-scale, that is, based on the thermal eld information of the LPBF process, the growth behavior of dendrites is described [42,43]. Yuan et al. [44] predicted the LPBF process at different scanning speeds based on the mesoscopic scale, and identi ed three mesoscopic molten pool states (unstable state, transition state, and steady state), and found that a stable molten pool with a small recessed area was bene cial to obtain a higher density part. The obtained simulation research results [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] provide a way to predict the temperature, molten pool ow, stress and other information of the LPBF process, but the current simulation research lacks a comprehensive analysis of the in uence of process parameters on the LPBF multi-layer and multi-path process based on the mesoscopic scale, which is also the focus of this paper. This paper was based on the open-source discrete element method (DEM) framework Yade and the open-source nite volume method (FVM) framework OpenFOAM to predict the multi-layer multi-path process of LPBF. In order to quantitatively analyze the effects of different process parameters (laser power, scanning speed, hatch space and layer thickness) on the porosity and surface roughness, based on the design of experiments (DOE) method, a four-factor threelevel orthogonal test scheme was designed, and the calculated results were compared with the experimental results. This paper is expected to provide a basis for process control in actual LPBF production.

Mathematical Model
(1) Modeling the spreading powder process The premise of describing the LPBF process based on the mesoscopic scale is to obtain the particle distribution of the powder bed. As a numerical calculation method for solving and analyzing the motion law and mechanical characteristics of complex discrete systems, DEM can be used to characterize collision and friction behaviors in particle systems [45]. The LPBF spreading powder process was predicted here based on DEM [39], and the calculation process was divided into the following four steps: (a) The metal particles were regarded as elastic spheres with contact friction effect, and the particle size distribution curve of the powder was obtained by means of experimental methods or literature data; (b) A particle cluster was generated in a certain space above the substrate and spread loosely on the substrate under the action of gravity; (c) The roller was moved to spread the powder, during which the contact friction between the roller and the particles and between the particles were considered; (d) After spreading powder, the sphere center coordinates and radius of the powder bed particles were derived.

(2) Laser heat source model
The way the laser energy is applied is very important to describe the LPBF process based on the mesoscopic scale. In order to ensure the reasonable application of laser energy and high calculation ef ciency, the multi-re ection absorption effect was ignored and the "metal phase-gas phase" interface was traced rstly in real time, then the elements directly acted by the laser were determined according to the laser spot center and radius. Using the elements directly acted by the laser as the starting points, the elements within a certain distance along the direction of gravity were found and marked as the elements indirectly acted by the laser, and nally the laser energy was applied in the form of a body heat source. Eqs. (1)-(5) is the laser heat source model established herein, in which, f metal represents the proportion of laser energy in the horizontal local area of a certain element acted by the laser, and Eqs. (4) and (5) is the volume of uid (VOF) model [46]. It should be noted that Eqs. (1)-(3) are original.
where, Q laser is the laser body energy density, W/m 3 ; f metal is the laser energy distribution coefcient; q laser is the laser surface energy density, W/m 2 ; ∆z is the equivalent size of the element, such as the side length of a cube, m; η is the laser absorption rate of the metal; P laser is the laser power, W; R is the spot radius, m; x, y are the horizontal coordinates of the element's center, m; x 0 , y 0 are the horizontal coordinates of the spot's center, m; v laser is the scanning rate, m/s; t is the time, s; V metal is the metal-phase volume of the element, m 3 ; V sum is the sum of the metalphase volumes of the elements acted by the laser under the same local horizontal area, m 3 ; α 1 , α 2 are the volume fractions of the metal-phase and gas-phase of the element, respectively; u is the velocity, m/s.  12) is the mass conservation equation. In order to ensure the ef ciency of numerical solution, two assumptions were made: the ow of liquid metal and gas were regarded as laminar ow of incompressible Newtonian uid, and the mass loss caused by metal gasi cation was not considered.
where, ρ, ρ 1 , ρ 2 are the mixed, metal and gas densities, respectively, kg/m 3 ; T is the temperature, K; c e , c 1 , c 2 are the equivalent [41], metal and gas speci c heat capacities, respectively, J/(kg · K); k is the mixed thermal conductivity, W/(m · K); h con is the interface convection heat transfer coef cient, W/(m 2 · K); T con is the external convection temperature, K; σ s is the Stefan-Boltzmann constant, W/(m 2 · K 4 ); ε is the emissivity; T rad is the outside radiation temperature, K; L gas is the metal vaporization latent heat, J/kg; m is the metal molecular mass, kg; k B is the Boltzmann constant, J/K; P 0 is the standard atmospheric pressure, Pa; T gas is the metal vaporization temperature, K; p is the pressure, Pa; τ is the stress tensor; µ is the hybrid viscosity, Pa·s; g is the acceleration of gravity, m/s 2 ; K C is the drag coef cient of the mushy zone [47], 1/s; f liquid is the liquid-phase ratio; C K is a self-de ned smaller value, such as 10 −6 ; σ is the surface tension coef cient, N/m; dσ dT is the coef cient of surface tension change with temperature, N/(m·K); κ is the interface curvature, 1/m; n is the normal vector of the interface.

Figure 1:
Calculation ow of LPBF multi-layer multi-path process: (a1-a3) spreading the rst layer of powder bed (particle size distribution → generation a particle cluster → complete spreading powder); (b1-b3) forming the rst layer of powder bed (geometry model → CFD simulation → STL model of the rst formed layer); (c1-c3) spreading the second layer of powder bed (particle size distribution → generation a particle cluster → complete spreading powder); (d1-d3) forming the second layer of powder bed (geometry model → CFD simulation → STL model of the second formed layer)

Numerical Solution
The open-source DEM framework Yade was used to calculate the spreading powder process, and the open-source FVM framework OpenFOAM was used to predict the LPBF multi-layer multi-path process. Fig. 1 shows the calculation ow of LPBF multi-layer multi-path process, which mainly includes the following steps: (a) Import the particle size distribution data into Yade, generate a particle cluster on the substrate, then move the roller to complete the spreading powder work of the rst layer, and export the particles' center coordinates and radius data of the rst layer of powder bed; (b) Import the particle information of the rst layer of powder bed into the 3D modeling software to obtain the geometric model, develop the LPBF process solver based on OpenFOAM, predict the forming process of the rst layer of powder bed, and derive the STL model of the rst formed layer; (c) Import the STL model obtained in step (b) into Yade, calculate the spreading powder process for the second layer, and derive the particles' center coordinates and radius data of the second layer of powder bed; (d) Import the particle information of the second layer of powder bed into the 3D modeling software to obtain the geometric model, predict the forming process of the second layer based on OpenFOAM, and derive the STL model of the second formed layer; (e) Repeat steps (c, d) to predict the forming process of the subsequent layer. It should be noted that the molten pool had been cooled down to the room temperature when anther spreading powder process started.

Calculation Parameters
The metal powder used here was 316L stainless steel, and its alloy composition (mass percentage) is: Fe 65.395%-Cr 17.0%-Ni 12.0%-Mo 2.5%-Mn 2.0%-Si 1.0%-P 0.045%-C 0.03%-S 0.03%. Tab. 1 shows the physical properties of 316L stainless steel and air obtained with the software JMatPro v7.0. The particle size distribution of the stainless steel powder satis ed a Gaussian distribution (the central value is 25 µm, and the variance is 5 µm), and the particle size was arti cially limited between 15 µm and 35 µm. When simulating the spreading powder process based on Yade, the set material parameters included density (7270 kg/m 3 ), contact friction angle (0.1), Young's modulus (195 GPa) and Poisson's ratio (0.3).
The focus of this paper is to analyze the in uence of process parameters on the LPBF multi-layer multi-path process based on the mesoscopic scale, and the considered process parameters include laser power, scanning speed, hatch space and layer thickness. In order to analyze the sensitivity of each process parameter, the orthogonal experiment method for the design of experiments was adopted, and Tab. 2 is the four-factor and three-level parameter table of process parameters, and Tab. 3 is the corresponding orthogonal test schemes. The reason for choosing the parameter values in Tab. 2 is that these parameter values are commonly used in actual LPBF schemes. Fig. 2 is a schematic diagram of the multi-layer multi-path forming process of each calculation scheme. Three layers of powder bed were formed together, and the laser scanned three paths for each layer of the powder bed. In order to control the size of the calculation domain as much as possible, the geometric dimensions of the calculation domain of each calculation scheme were X-direction (650 µm), Y-direction (2 * hatch space + 100µm) and Z-direction (layer thickness * current forming layer number + 100µm). In addition, the laser spot diameter was 70 µm, the mesh size was 2.5 µm, the calculation time step was 4 × 10 −7 s, the laser absorption rate was 0.35, the ambient temperature was 300 K, and the computing resource used was con gured as Intel Xeon Gold 5120 CPU (dual CPU, 56 threads, 96 GB memory).

Extraction of Porosity and Surface Roughness Data
Firstly, the simulation results were illustrated. Here, calculation Scheme 1 was taken as an example. Fig. 3 shows the temperature and molten pool distributions when forming each layer and each path. From the calculation results, it could be seen that, on one hand, the metal particles were rapidly melted by the laser heating, forming a molten pool, and even gasi cation; on the other hand, as the laser beam was in motion, the overheated molten pool continuously conducted heat to the surrounding metal and radiated heat to the external environment, making the molten pool in a quasi-steady state during the forming process. So it could be considered basically that the shape of the molten pool had always maintained a "droplet" shape. Fig. 4 shows the surface morphology of each formed layer. It could be clearly seen that there were bumps between adjacent solidi ed tracks. The reason for the bumps is that, on one hand, the surface of the molten pool formed a depression under the action of the gasi cation recoil, the so-called pinhole effect, so that the liquid metal in the molten pool was squeezed from the center of the laser spot to the surroundings; on the other hand, because the temperature of the molten pool surface presented an obvious distribution of "high spot center-lower surroundings," the liquid metal near the molten pool surface tended to ow from the spot center to the surroundings under the Marangoni effect. Under the combined action of the above two reasons, the height of the side of the molten pool was higher than that of the adjacent solidi ed track, and the side was too late to recover during the subsequent cooling and solidi cation process, and a bump was formed. It should be noted that in the actual LPBF process, spatter and powder bed erosion also have an important in uence on the forming process [48,49]. In summary, the simulation results showed that the LPBF multi-layer multi-path process were very complicated. In order to realize the sensitivity analysis of the process parameters, it is necessary to extract the key information from the simulation results to achieve the purpose of quantitative analysis. The porosity and surface roughness of the parts are the main quality indicators of concern in actual production, so the porosity and surface roughness were token here as the extraction information of the simulation results. Taking calculation Scheme 1 as an example, Fig. 5 shows the nal gas-solid interface obtained after forming three layers and three path, and the method to obtain this interface was to extract the isosurface of α 1 = 0.5. It could be seen from Fig. 5 that since the laser did not act on the front and rear ends of the X and Y directions, a large number of gas-solid interfaces appeared in these areas, so these areas needed to be avoided when extracting information on porosity and surface roughness. Fig. 6 is a partial volume area used for calculating porosity, and the top surface of the partial volume area needed to be lower than the upper surface of the formed zone in Fig. 5. The post-processing software ParaView v5.4.0 provided the element information summing operation for the local volume area, so the calculation method of porosity is as follows: (15) where, V i is the volume of the element, m 3 ; φ i is the gas volume ratio of the element; N is the number of elements; V total is the total volume of the statistical area, m 3 ; V gas is the total volume of the gas-phase in the statistical area, m 3 ; φ is the porosity of the statistical area. V total and V gas could be provided by ParaView, and the extracted porosity data of each calculation scheme can be seen in Tab. 4.  For the extraction of surface roughness data, a partial surface area (Fig. 7) was cut out from the upper surface of the interface in Fig. 5, and it was thought that the degree of uctuation of the partial surface area could re ect the surface roughness of the part. Fig. 8 shows the local area used to calculate the surface roughness under calculation Scheme 3. It could be clearly seen that due to the large hatch space at this time, there were metal particles that were not completely melted between the solidi ed tracks, and they would have a signi cant impact on the surface roughness, which was consistent with the actual situation. It should be noted that the size of the local area analyzed herein was limited to a few hundred microns, which is quite different from the actual surface size of the part, therefore, the so-called surface roughness data extracted here could only be regarded as a certain mapping of the actual surface roughness of the part, and the speci c value could not represent the actual surface roughness. ParaView provides the operation of extracting the node information of the partial surface area, so the calculation method of surface roughness is as follows: where, z i is the height coordinate of the node, m; m is the number of nodes; z average is the average height of the nodes in the local area, m; S rough is the surface roughness, m. z i could be provided by ParaView, and the extracted surface roughness data of each calculation scheme can be seen in Tab. 4.

Sensitivity Analysis of Process Parameters (A) The in uence of different levels of process parameter on porosity
Based on the porosity and surface roughness data of each calculation scheme in Tab. 4, the in uence of process parameters on the LPBF multi-layer multi-path process could be quantitatively analyzed. For each process parameter, the in uence of different level parameters on porosity was analyzed. Taking laser power as an example, according to the data analysis strategy of the orthogonal experiment method, the level 1 parameter of laser power was used in calculation Schemes 1-3, so the porosity data of calculation Schemes 1-3 should be summed and averaged, and the results re ected the in uence of the level 1 parameter of laser power on porosity. Similarly, the in uence of level 2 and 3 parameters of laser power on porosity could be obtained (Fig. 9a). According to this strategy, the in uence of different level parameters of scanning speed (Fig. 9b), hatch space (Fig. 9c) and layer thickness (Fig. 9d) on porosity could be obtained. It can be seen from the data in Fig. 9 that the effects of laser power, scanning speed, hatch space and layer thickness on porosity all exhibited monotonicity. Among them, as the laser power increased, the porosity decreased, and as the scanning speed, hatch space and layer thickness increased, the porosity increased. Experimental research [50,51] also revealed the same law. The reason for this law is that the larger the laser power, the smaller the scanning speed, the hatch space and the layer thickness meant that the metal particles in a certain volume absorbed more laser energy, so that the metal particles could be more fully fused, and the gas in the molten pool had more suf cient escape time. It should be noted that since the highest laser power selected herein was only 250 W, the data in Fig. 9a did not re ect the entrainment phenomenon caused by the excessively high laser power, which easily causes the porosity to increase. With reference to the above data analysis strategy, Fig. 10 shows the in uences of different level parameters of laser power (Fig. 10a), scanning speed (Fig. 10b), hatch space (Fig. 10c) and layer thickness (Fig. 10d) on surface roughness. It can be seen from the data in Fig. 10 that the in uence of laser power and scanning speed on surface roughness showed a trend of rst decreasing and then increasing, while the in uence of hatch space and layer thickness on surface roughness showed a monotonous increasing trend. Experimental research [50] also revealed the same law. The reason is that, on one hand, too low laser power and too high scanning speed meant that the laser energy applied to the powder bed per unit time was too little, and it was dif cult to fully fuse metal particles to obtain a rough surface. Too high laser power and too low scanning speed meant that too much laser energy was applied to the powder bed per unit time, which caused the pinhole phenomenon and the Marangoni effect to be too signi cant, thereby obtaining an uneven solidi ed track surface. On the other hand, because the hatch space and layer thickness were not too small, as the hatch space and layer thickness gradually increased, the laser energy applied to a certain volume of the powder bed was reduced, and it was dif cult to fully fuse the metal particles to obtain a rough surface. Then, the degree of in uence of each process parameter on porosity and surface roughness was analyzed, that is, comparing each process parameter horizontally, and ranking the in uence of each process parameter on porosity and surface roughness. Taking the in uence of laser power on porosity as an example, according to the data analysis strategy in the orthogonal experiment method, the range (the difference between the maximum value and the minimum value) in Fig. 9a could represent the in uence of laser power on porosity. According to this strategy, the in uence of each process parameter on porosity (Fig. 11a) and surface roughness (Fig. 11b) could be obtained. It can be seen from the data in Fig. 11 that the order of the in uence of each process parameter on porosity is: scanningspeed > layerthickness > laserpower > hatchspace, and the order of the in uence of each process parameter on surface roughness is: hatchspace > layerthickness > laserpower > scanningspeed, that is, the porosity of the part is most sensitive to scanning speed, and the surface roughness of the part is the most sensitive to hatch space. The order of the in uence of process parameters obtained herein was not consistent with [23]. The reason is that the degree of in uence of each process parameter on porosity and surface roughness was closely related to the speci c values of the parameters, and the ranges of the parameters' values set herein were not completely consistent with the ranges of [23]. The above-mentioned data analysis conclusions help to provide a basis for process control for actual LPBF production of the 316L stainless steel alloy. OpenFOAM, the multi-layer multi-path forming process of LPBF was predicted herein. In order to quantitatively analyze the in uence of different process parameters (laser power, scanning speed, hatch space and layer thickness) on the LPBF process, based on the DOE method, a four-factor three-level orthogonal test scheme was designed, and the porosity and surface roughness data of each calculation scheme were extracted. (2) By analyzing the orthogonal test data, it was found that as the laser power increased, the porosity decreased, and as the scanning speed, hatch space, and layer thickness increased, the porosity increased. In addition, the in uence of laser power and scanning speed on surface roughness showed a trend of decreasing rst and then increasing, while the in uence of hatch space and layer thickness on surface roughness showed a monotonous increasing trend. (3) The order of the in uence of each process parameter on porosity is: scanning speed > layer thickness > laser power > hatch space, and the order of the in uence of each process parameter on surface roughness is: hatch space > layer thickness > laser power > scanning speed. The above conclusions are helpful to provide process control basis for actual LPBF production of the 316L stainless steel alloy.
Funding Statement: This work was supported by the Natural Science Foundation of Guangdong Province (No. 2019A1515012040).

Con icts of Interest:
The authors declare that they have no con icts of interest to report regarding the present study.