A new data-driven framework for prediction of molten pool evolution and lack of fusion defects in multi-track multi-layer laser powder bed fusion processes

This study proposes a new data-driven approach to predict the area fraction of lack of fusion and morphology of 3D-printed parts produced using the laser powder bed fusion process. Different processing parameters were considered to develop this methodology, such as laser power, laser speed, hatch spacing, and layer thickness. Melt pool geometrical parameters were numerically predicted using a mathematical model optimized using the Nelder-Mead optimization method and a machine learning model. The Gaussian process was later employed to further modify the layerwise deposition. The accuracy of the proposed methodology was then extended to predictions of the area fraction of lack of fusion and the morphological aspects of the melt pools in 3D-printed parts by comparing them with the experimental observations. The results revealed that this data-driven approach, which is independent of the thermophysical properties of the material, is capable of reasonably predicting the area fraction of lack-of-fusion formed during material deposition using the laser powder bed fusion process.


Introduction
Laser powder bed fusion (LPBF), a track-by-track followed by a layer-by-layer fabrication process, is a promising approach to producing customized parts with complex geometrical features [1] through computer-aided design (CAD) models [2]. LPBF is considered a suitable approach in aerospace and medical applications [3] and in producing mission-critical parts remotely [4] if a high level of quality assurance and process reliability [5] could be obtained. One main challenge in the LPBF process is to manufacture defect-free parts [6][7][8]. Mainly porosity forms due to lack-of-fusion, balling phenomenon, and keyhole effect [9].
Mechanical properties of the 3D-printed parts are highly affected by lack-of-fusion during processing [10,11]. Incomplete melting of the powder material forms lack-offusion porosity [12,13], which could be due to improper selection of the processing parameters, such as scanning speed, laser power, laser spot size, hatch spacing, and layer thickness [14]. The parts produced in high-quality require a great understanding of defect formation mechanisms and the approaches to avoid or minimize defects. Also, optimization of the processing parameters is mainly based on a significant number of trial and error cycles [15,16]. Therefore, considering all the complex mechanisms, parameters, and interactions of parameters, affecting the quality of 3D-printed parts is challenging and costly.
To improve the quality of 3D-printed parts, the processing parameters, affecting the melt pool dimensions and defects formation mechanisms are required to analyze [17][18][19][20][21][22]. Simonds et al. [23] quantified the relationship between the laser power absorption and melt pool geometry of the samples made of Ti-6Al-4V powder. Also, the effects of laser power and laser speed on the melt pool dimensions of the samples made out of steel alloy EN25 powder [24] and Inconel 718 [25] were evaluated. No clear relationship was observed between melt pool dimensions and build height [25]. A power-velocity process map was generated for controlling the melt pool dimensions in a study by Soylemez et al. [26]. Yadroitsev and Smurov [17] indicated that the processing parameters could have a stability zone, which forms defect-free tracks, and Beuth et al. [18] illuminated process mapping for several additive manufacturing (AM) systems. The volume fraction of lack-of-fusion porosity formed due to insufficient overlap of melt pools was predicted [19], and the processing parameters were optimized for a faster build rate [20]. In general, more than 130 processing parameters, affecting the final quality of the LPBF parts, have been recognized [27], which make the design challenging [28]. Also, obtaining reliable experimental data and optimizing the processing parameters are time-consuming, costly [29], and this knowledge is not mostly transferrable from one material to another [30].
One way to optimize the processing parameters for improving the quality of 3D-printed parts and process reliability is the application of in situ monitoring systems [31], which requires an efficient defect detection method, such as image processing. On the contrary, machine learning (ML) approaches can be employed to address the need of analyzing data [29,32,32] and develop predictive models to determine the optimal processing parameters in an efficient way [33]. Kamath [28], for example, applied regression trees and Gaussian process models to predict the melt pool depth. Tapia et al. [34] also developed a Gaussian process-based surrogate model to predict the melt pool depth. Meng et al. [35] developed a Gaussian process-based model to predict the remelted depth. Lee et al. [36] employed Bayesian ridge regression, kernel ridge regression, linear regression, nearest neighbors regression, random forest regression, and support vector machine models to predict the melt-pool geometries of a singletrack deposition using LPBF. Oh and Ki [37] developed two successive generators to predict the melt pool cross-section based on a convolutional neural network (CNN) for a singletrack deposition. The artificial neural networks (ANNs) model, a machine learning computational approach that consists of several processing elements, was also employed to predict the melt pool depth and peak temperature [38], melt pool width and height [39], and the width, depth, and height of the deposited track produced by direct metal deposition [40]. ML approaches help correlate the processing parameters with the output parameters (i.e., melt pool dimensions), leading to identifying the mechanisms, affecting the porosity in 3D-printed components.
Predicting the morphology of multi-track and multilayer samples is difficult. The studies reviewed earlier, excluding the study by Oh and Ki [37], which predicted the melt pool shape in a single-track deposition, only predicted the melt pool dimensions but not the shape of the melt pools. Predicting the shape of a melt pool can provide more intuitional results, as this type of data contains more information than the numerical values of melt pool geometries (i.e., width and depth). Moreover, predicting the morphology of multi-tracks and multi-layer 3D-printed parts is even more desirable because this is the final result of the printing process. However, so far, limited studies have developed a model to predict the morphology of multitrack and multi-layer 3D-printed parts. In light of this fact, this study proposes a new methodology to predict the melt pool cross-sections and consequently morphology of multi-track and multi-layer 3D-printed parts produced using the LPBF process at different processing parameters. The new methodology helps estimate the area fraction of lackof-fusion in these parts at various processing parameters. In light of this fact, this study aims to develop a new methodology to (1) estimate the area fraction of lack-offusion in these parts at various processing parameters and (2) predict the melt pool cross-sections and consequently morphology of multi-track and multi-layer 3D-printed parts produced using the LPBF process at different processing parameters. To accomplish this, a range of volumetric energy density was employed along with the melt pool boundary curves extracted from the micrographs taken from the samples to form a function optimized later with the method of least squares [41], Nelder-Mead [42], and random forest regression methods [43]. Eventually, the melt pool boundaries (considering bi-directional scanning melt pool layers with 90 rotation) and powder layers were further optimized using the Gaussian process regression [44,45], and the layers were added one-by-one mimicking the printed process.
In this study, the data preparation procedure was discussed in Section 2. The data was used for the final plot in Section 3. In Section 4, two ML models were evaluated for data prediction, and the area fraction of the lackof-fusion measured from the micrographs and predicted from computationally plotted images were compared together.

Experiment and surrogate models
To develop a predictive model, a set of experiments were conducted. After taking the optical micrographs, the numerical data required to predict the melt pools crosssections and morphology of the parts were extracted from the micrographs. These data were used to train two separate models (i.e., a mathematical melt pool boundary curve model and the random forest model), which were later used to numerically predict the melt pool shapes.

Experiment and micrographs
LPBF, a promising AM process due to its flexibility in feedstock and printed shapes [46][47][48][49], was employed to manufacture 33 cubes (cube size = 2 2 2 3 ) made of stainless steel, SS316L with different processing parameters (i.e., laser power, laser speed, hatch space, and layer thickness) listed in Table 1. The ranges of the processing parameters listed in Table 1 were selected in a way that causes the formation of lack-of-fusion defects in a wide range (i.e., low to high) in terms of total area fraction during manufacturing. Creating a dataset, covering a wide range of area fractions is needed for developing the model more accurately. The samples were then scanned to make electronic images with the size of 2815 pixels 2112 pixels, which were later used to train two separate models.
A bi-directional scanning strategy (rotated by 90 ) was utilized during the LPBF-based 3D printing, as shown in Fig. 1. If the scan direction is parallel to the -axis, the cross-sections (transverse cross-section) of all the melt pools on the layer can be observed on the -plane whereas the longitudinal profile of just a melt pool can be observed on the -plane, which looks like a long strip. These observations are shown in Fig. 1 Fig. 2, the top layer, where the melt pools are observed in full depth, was aligned with the -axis in each micrograph. For the rest of the layers aligned withaxis, the scanning order was from left to right. Therefore, the right corner of each melt pool overlapped with its rightside adjacent one except for the right-most melt pool in each layer. With this in mind, the geometrical information of the melt pools associated with different processing parameters is extracted, which can be fed as input to the surrogate models.

Input data
In AM, the printed part quality highly depends on the values of the processing parameter (e.g., laser power ( ), laser  [50][51][52]. To better evaluate the effects of processing parameters on the output, volumetric energy density ( ) was introduced as follows [53]: (1) Therefore, along with , , , and will act as descriptors for melt pool shape variations observed in the optical micrographs.

Melt pool boundary shape
The melt pool boundary shape was extracted from the topmost melt pool layer of each image manually. The extracted melt pool boundary shape data is incomplete because the right corner of the melt pool overlaps with the melt pool on its right side. To solve the issue, the melt pool boundaries were extended using Gaussian process regression [44]. All melt pools were moved to the center of the -plane, which means the centers of the melt pools are exactly the origin.

Melt pool width and depth
To predict the melt pool shape, the width and depth of the full-melt pools were measured from the micrographs (Fig. 2). In each micrograph, width ( ) and depth ( ) of all the melt pools on the top layer were measured, as shown in Fig. 2. The means of the width and depth values were then computed and denoted as and , respectively.

Standard deviations of geometries of the melt pools
In order to account for the melt pool overlapping (hatch space) effects during the layer-by-layer LPBF process, two parameters were measured from the melt pools within the layers other than the top layer: horizontal distance between the visible (left) corners of the adjacent melt pools ( ) (bottom inset of Fig. 2), and depth of the melt pools within the layers other than the top layer ( ) (bottom inset of Fig. 2).
In each optical micrograph, 50 melt pools were randomly selected within the layers other than the top layer, and the two parameters listed above were measured from these melt pools. The standard deviations of the two parameters and were then computed. The means of the horizontal distance between the corners ( ) and the mean of the depth of the melt pools ( ) were discarded, as and are estimates of and .

Melt pool boundary shape prediction
To predict the melt pool boundary shape, each point on the melt pool boundary curve extracted from the micrographs was converted into the polar coordinate system , where is the polar angle in radians and is the radial distance in micrometers. The values of the polar angles were restricted in the interval , as function , which is introduced later, may not have a period of 2 for , and the continuity of the lower half of the graph should be guaranteed. Considering the corresponding value of volumetric energy density ( ), a set of triplets 1 can be obtained, where is In this model, 1 , 2 , 3 , 1 , 2 , 3 , and are constants to be determined later with the method of least squares.
Considering as the sum of squares defined as follows: where is the vector, containing the parameters 1 , 2 , 3 , 1 , 2 , 3 , and . Nelder-Mead optimization method [42] was applied to minimize and to estimate the values of these parameters. Once all the values are found, then can be used to draw the melt pool boundary shape for the given value of volumetric energy density. The approach proposed in this study is indicated in a flowchart (Fig. 3).

Prediction of w, d , b , c
In each optical micrograph, four classes of parameters can be obtained, namely average width ( ), average depth ( ) of the full-depth melt pools, the standard deviations of the horizontal distance between the adjacent visible corners ( ), and the standard deviation of the depth ( ) of the melt pools within the layers other than the top layer. These four parameters are sample means and sample standard deviations, which are often used to estimate the corresponding true statistical values ( , , , ). Therefore, these data along with their corresponding input values ( , , , , and ) can be used to train the regression models. The regression model applied to predict these four parameters is similar to the algorithms introduced in [28,[34][35][36][38][39][40] and is not discussed here in this study. Then, the predicted values were used to draw the morphology of the 3D-printed part.
Seven different regression models were compared with the same dataset with a 5-fold cross-validation repeated 20 times. The average mean squared errors (MSE) [54] and average mean absolute percentage errors (MAPE) [54] were listed in Table 2   MAPE, random forest regression [55] was selected to predict these parameters, which are used for the prediction of the morphological aspects of the melt pools in the multitrack and multi-layer parts. Random forest regression is an ensemble of regression trees [56]. Each tree splits the data into smaller partitions iteratively by the splitting rule [57] with some random factors, and each terminal node contains a simple model to predict the data [58]. The mean of all the values predicted by the regression trees is taken as the final predicted value. The random forest regression model in this research contains 200 trees without an indication of the maximum depth. The predicted values for , , , and were denoted as , , , and , respectively.

Melt pool shape prediction
Predicting the morphology of multi-track and multi-layer 3D-printed parts involved multiple sequential steps, as indicated in Fig. 3. First, the melt pool shape predicted using the mathematical model was adjusted with the predicted values of and . Then, the predicted values and were used to generate noise data (variations), which was later employed to further optimize the melt pool boundaries. Finally, the morphology of the 3D printed part was plotted layer-by-layer using the adjusted melt pool shape and noise data.

Melt pool boundary shape adjustment
After comparing the values of width and depth extracted from the predicted melt pool boundary shape with the width ( ) and depth ( ) predicted with the random forest regression model, it was found that the error of the melt pool boundary shape model was much larger than that of the random forest regression model. One of the reasons for this discrepancy is that the melt pool shape model considered only the volumetric energy density ( ) as the input; however, all of the five parameters mentioned earlier were considered as the input parameters in the random forest regression model. To reduce the error of the model, the width and depth of the predicted melt pool boundary shape were required to adjust. To do this, the width and depth of the predicted melt pool boundary shape were set at and , respectively. If the input values are considered as 0 , 0 , 0 , 0 , and 0 , the melt pool width and depth can be predicted using the random forest regression model. Also, the melt pool boundary shape ( ) computed by can be presented by a collection of points: As shown in Fig. 4(a), 1 is quite arbitrary and different from the observed melt pool profiles. Furthermore, this shape might be different if different values are chosen for 1 , 2 , 3 , 1 , 2 , 3 , and . The interval was considered unchanged in both 1 and the training data for to guarantee the continuity and smoothness for the lower half of the graph of 1 . Even though the section under the dashed line ( -axis) is the primary section corresponding to a melt pool boundary curve, the whole loop is kept for future use.
The depth of the predicted shape was also calculated as: 0 2 . To relocate the graph along the z-axis, 1 was then converted into the form of Cartesian coordinates, and then all the points were shifted up along the -axis for units. The graph of 2 is also presented in Fig. 4 2 is the height of the whole curve 2 , then 2 intersects the -axis at two symmetrical points. This condition never fails throughout the whole testing process. If is the -coordinate of the intersection on the positive side of -axis, the graph of 2 then stretches or shrinks by the factor of with respect to the -axis.
Similarly, the graph of 3 (Fig. 4(c)) under the dashed line is the melt pool boundary shape with width of and depth of . Accordingly, 2 and 3 are represented by the following Eqs. 7 and 8, respectively.

Prediction of melt pool variations
For each alternative layer, two different types of curves are required to accurately mimic transverse and longitudinal pool profiles observed in Fig. 2. One curve represents the melt pool boundaries (solid curves in transverse crosssection YZ plane in Fig. 5), and the other one indicates the longitudinal profile of the melt pool (bound between solid and dashed curves in YZ plane in Fig. 5) corresponding a track rotated through 90 (bidirectional scanning strategy). The variations for these two types of curves were generated independently with the approach of the Gaussian process [44,45]. A Gaussian process is a stochastic process, consisting of a collection of random variables, such that any finite number of the random variable has a joint Gaussian distribution [44]. It is specified with a mean function and a covariance function , which are defined as: and the Gaussian process is written as follows: The covariance function of the Gaussian process for the variations on the horizontal wavy lines is defined as follows: As expressed in Eq. 12, is a squared-exponential function with length-scale and signal variance 2 . 1 and 2 can be any two -coordinates, and 1 2 is the the covariance of the values of the Gaussian process to be defined at 1 and 2 . Similarly, the covariance function of the Gaussian process for the variations on melt pool shape curves is defined as: . Similarly, and 2 are length-scale and signal variance, respectively. The length scales in Eq. 12 and in Eq. 13 were chosen to be 100 and 0.6, respectively, based on the experimental observations from and masurements. Also, in Eq. 12 and in Eq. 13 were calculated based on the values of and , respectively, which are discussed later in Section 3.4. Furthermore, the Gaussian processes for the variation on the horizontal wavy lines are defined as: where is the total number of horizontal wavy lines with variation in the plotted image. In Eq. 14, 1 were indicated as Gaussian processes with the mean function 0 and covariance function , and they were generated independently. 1 were added to straight horizontal lines to make these lines similar to the dashed lines shown in Fig. 5. Similarly, the Gaussian processes for the variation on the melt pool boundary curves are expressed as: where is the total number of melt pools plotted in the generated image. Again, the Gaussian processes 1 were generated independently. These variation data were added to the melt pool boundary curve 3 in the polar coordinate system. However, before adding the variations (Gaussian process) to 3 , 3 is required to be converted into the polar coordinate system. . This function, , can be represented as: Accordingly, the graphs of 3 and 4 are illustrated in Fig. 4(c). The only difference is that the former is in a Cartesian coordinate system and the latter is in the polar coordinate system. The variation obtained by the Gaussian process was then added to 4 to obtain 5 .
The graph of 5 was presented in Fig. 4(d). In Eq. 18, in 5 is constrained within the interval to restrict the generated graph in the lower half region, which is relevant for melt pool shape prediction. Once the noise data is added to the melt pool boundary curve 4 , the radial coordinate of the top section may switch to negative. Therefore, after adding the noise data, the top section is removed from the sum to avoid negative values on the radial coordinate.
Eventually, 5 was converted into a Cartesian coordinate system, and the center of the melt pool boundary curve was shifted from the origin to a different location as follows: 6 cos sin 5 .
As shown in Fig. 4(e) and (f), the shape of 5 is the same as the shape for 6 . The melt pool 5 was relocated because each melt pool should not only contain the unique variation data but also be located at a different place. Finally, 6 is the desired melt pool curve to generate the fingerprint of the melt pools within a given layer. 6 with the top section removed is shown in Fig. 4(f). 6 was first cut by a horizontal wavy line, and the part above this line was then erased. Only drawing the melt pool curve below a horizontal line is a step in plotting the final image (Section 3.3). Each point on this horizontal line was normally distributed with mean (the -coordinate of the center point) and some fixed variance value. Therefore, this line should fluctuate around . In this case, the section of the melt pool data should be kept completely below and above , which is why the section above the dashed lines was not completely removed in 1 through 6 .

Steps of image prediction
After estimating the melt pool shape using the trained data, the generated melt pools were plotted one by one to produce the morphology of an entire geometry of a multi-track and multi-layer 3D-printed component in the cross-section. The order of the plotting procedure was the same as the sequential steps followed during actual printing (Fig. 6). As the laser scan direction was rotated by 90 after printing each layer, the transverse (the top layer in Fig. 6(k)) and the longitudinal profile of the melt pool layers (the top layer in Fig. 6(m)) were plotted interchangeably.
Before plotting each layer, a powder layer was plotted on top of the previous layer in advance (the top layers in Fig. 6 The whole procedure of image prediction starting from the bottom layer through the top layer Fig. 6(a), (f), (h), (l), and (n)). The powder layer and/or lackof-fusion areas were indicated in black, while the melted areas were shown in white. The method of plotting the powder layer for any non-bottom layer was also the same, and the area between the top line of the current layer and the top line of the previous layer was presented in black. To create the bottom layer, as the powder was spread on the substrate, the area between the top line and the substrate was filled in black ( Fig. 6(a)). The substrate is a straight (flat) horizontal line with no variations, whereas the top line of any powder layer contains variations. The height of any powder layer is expected to be the layer thickness , which was considered a processing parameter.
Similar to any other lines drawn for image prediction, the top line is the sum of two components: the -coordinate of the flat horizontal line, which was denoted as and the noise data, which is a Gaussian process from the remaining ones of 1 (discussed in Section 3.2) with no repetition. Here, selection without repetition means that once a Gaussian process is selected to plot the image, this Gaussian process will never be selected again for further use. It is worth mentioning that was a single value and could also be considered as the expected mean value of the -coordinates for each layer. The value of changed at different layers, and the -coordinates of the substrate along with the values of for all layers constituted a finite arithmetic sequence with the common difference .
After plotting a powder layer, either the layer, representing the transverse profiles of the melt pools (the top layer in Fig. 6(k)) or the layer, representing the longitudinal profiles of the melt pools (the top layer in Fig. 6(m)) was plotted to make the powder layer melt. When the laser scan direction is along the -axis, the layer, representing the longitudinal profile was plotted, and if the scan direction is along the -axis, the layer, representing the transverse profiles of the melt pools was plotted. To plot either type of layer, the information of the powder layer z-value/height for the concerned layer number is required. This also acts as the top surface of the melt pools as indicated in Fig. 6 (e and o).
Plotting of the layer, representing the longitudinal profile of the melt pool is first discussed. In simple terms, two nearly parallel horizontal lines were drawn in black, and then the area between the lines was filled in white. The top of the corresponding powder layer was considered the top line. The bottom line was considered a new line with a variation. The expected mean value of the -coordinates of the bottom line was denoted as . The height of this layer is expected to be , which is the predicted value for the depth of complete melt pools. Therefore, . The variation on the bottom line was also selected from the remaining elements of 1 with no repetition.
The melt pools in a given deposited layer were plotted one by one from left to right, in the same order that they were printed. Therefore, the right corner of each melt pool overlapped with the one on its right side, which was magnified in Fig. 6(p). The variation on the melt pools was selected from the remaining 1 to avoid repetition, and the centers of the melt pools were located over the line . The -coordinates of the melt pool centers form a finite arithmetic sequence with the common difference in the hatch space value . In each melt pool, after the noise data and the central location were obtained, the variation in data was added to the original melt pool boundary curve 4 in the polar coordinate system. The melt pool was then shifted according to its central location. The steps followed in this process were the same as those followed through from 4 to 6 . It is worth mentioning that only the section of the melt pool boundary curve under the top line was plotted, which was marked in black. The interior area was bounded by the melt pool boundary curve, and the portion above the top line was then filled in white. The above process allows for differentiating the parameters associated with the melt pool, such as and .
Bi-directional scanning melt pool layers were plotted interchangeably. Once the last layer was plotted, the image prediction process is complete. The melt section within the substrate was disregarded, as it had little significance during multi-track multi-layer melting of subsequent layers. Similar to that in the optical micrographs, the entire part generated by the proposed approach shows full-depth layers, representing the transverse or longitudinal profiles of the melt pools, depending on the type of the last layer ( Fig. 6(p)). Therefore, with the help of predicted melt pool parameters and noise function, stitching of melt pool tracks in transverse and longitudinal layers can be achieved in a stepby-step process. Furthermore, the predicted melt pool shape parameters allowed visualization of the melt pools and melt pool overlap during subsequent layer formation along with a few non-melted regimes that often present lack-of-fusion related porosity (Fig. 6(p)).

Calculating A and B
As discussed earlier, in each micrograph image (Fig. 2), the horizontal distance between the left corners of adjacent melt pools ( ) and the depth of the melt pools ( ) in the non-top layers were measured, and their standard deviations ( and ) were computed (Section 2.3.3). Therefore, these two standard deviation values measured from the recorded images are expected to meet their corresponding computationally predicted values ( and ). As the sample standard deviation is a biased estimate of population standard deviation, sample variance was used. Hence, the following equations were obtained: where is the horizontal distance between the left corners of the and the 1 melt pools in the layer, is the depth of the melt pool in the layer (Fig. 7), and are the means of all the 's and 's, respectively.
is the total number of melt pools in the plotted image, which can be calculated as follows: (22) where is the number of layers appearing in the image, and is the number of melt pools in each layer. To find and indicted in Eqs. 12 and 13, and were replaced with the left side of Eqs. 20 and 21. Considering Eq. 20, in each non-top melt pool layer, the number of in each layer is one less than the number of melt pools in each layer. Thus, the total number of in the whole image is . For simple notation, let 1 and . Thus, Eq. 20 can be expressed as: For each and , can be written as: where stands for the horizontal distance between the visible corner and the center of the melt pool in the layer (Fig. 7). Also, the first term on the right-hand side of Eq. 23 can be written in the form of:  Fig. 7 Illustration of , , and in the non-top layer melt pools where 1 , 2 , and the in Eq. 26 denote the horizontal distances between the visible corner and the center of the 1 , 2 , and melt pools, respectively, in the same layer. Also, E 2 can be obtained as follows: Then, the expressions of E 2 1 , E 1 2 , and E 1 need to be obtained. To facilitate the calculation, some simplifications are required. The area enclosed by the melt pool boundary and the horizontal wavy line vary slightly in each track. Hence, in this region, the variation of any two points on the same curve should be very close, as the correlation of the two points is close to 1. Therefore, the covariance functions of the Gaussian processes for the variation on the melt pool boundary curve and the horizontal wavy line are close to the values of 2 and 2 , respectively. Based on this, the melt pool boundary curve was expressed by an equation of an ellipse with width 2 and height 2 2 (Fig. 8), where is the index of the melt pool and is a random variable, representing the variation on the melt pool boundary curve. The lower half of the ellipse is the major part of the melt pool, and the width and height of the lower half are ( 2 ) and where 0 is the center of the melt pool, and 0 2 . This will make sure the variation of the melt pool curve in the intersection region remain roughly the same. It was also assumed that in a given layer the bottom melt pool boundary in the longitudinal section can also represent a straight horizontal line of a melt pool bottom in the transverse section of the layer above this layer. As all the melt pools discussed in this section are located in the same layer, the -coordinates of the melt pool centers have the same value, which let us use a single variable 0 to represent this value. Considering as the vertical distance between the center of the melt pool 0 and the horizontal line around the region of the intersection (with subsequent layer) on the left side of the melt pool shown in Fig. 8, which can be expressed as 2 . Then, the horizontal distance between the intersection point and the center of the melt pool can be expressed using function , which is defined as follows: where can be the indices 1, 2, and . Only 's are dependent on each other. Also, covariance of 1 , 2 , and can be calculated as: and Cov 1 2 exp 2 and are also dependent. Any other pairs of the 's and 's are all independent. It is also required to make sure that the semi-major axis and the semi-minor axis of the ellipse are both positive. Thus: max 1 2 (32) and the intersection of the left side of the ellipse with the horizontal line is: Considering as the region that and can take values as follows: Therefore, E 2 1 can be expressed as: Similarly, considering as the joint probability density function of 1 , 1 , , and ( is an index of melt pool, which can be either 2 or ), which can be expressed as: Therefore, E 1 (42) Now, the left side of Eq. 20 can be expressed by and . The left side of Eq. 21 can also be expressed by: As the value of E 2 is the same for any and , this value is denoted by E 2 without indicating the indices. The depth of a melt pool in a non-top layer ( in Eq. 43 is the distance between the white triangle point on the lowest part of the melt pool boundary and the the white square point on the horizontal line right above it as indicated in Fig. 9, and the segment between these two points is a vertical line. The distance between these two points is approximately equal to the distance between the intersections of the vertical line through the center of the melt pool (the black circle in Fig. 9) with the melt pool boundary (the black triangle has partially been covered by the white triangle in Fig. 9) and with the horizontal line above (the black square has partially been covered by the white square in Fig. 9). Therefore, (44) where and are the distances between the melt pool center and the two intersection points described above, respectively, with 2 and  , and , which can be expressed as: if 0 otherwise (56) where is the pdf function of a multivariate normal distribution with mean and covariance matrix , where   (60) where x is the vector, containing all input variables, and k is a scaler presenting the dimension of x. and are the mean vector and covariance matrix of the random variables, respectively. A system of equations are formed using Eqs. 20 and 21 with two unknowns and . This system of equations can be solved with Newton-Raphson method,

Evaluation
To evaluate the above-described methodology, the area fraction of lack-of-fusion detected in the micrographs was compared with that in the computationally generated images. Due to the limited number of micrographs (33 images), leave-one-out cross-validation (LOOCV) [60] was applied during the evaluation process. For each round of the cross-validation, the training dataset was applied to train the model, which was developed mathematically, and the random forest model to predict melt pool boundary curves. Then, five input values from the test micrograph were used to plot the images with the trained predictive models throughout all the steps introduced in Sections 2 and 3. The micrographs and their corresponding generated images are presented in Fig. 10. It was found that the total area fractions of lack-of-fusion in the computationally generated images were roughly at the same levels as those observed that in the micrographs in consideration. However, a visual comparison may not be a proper reference to comment on the accuracy of the approach proposed in this study. To evaluate the accuracy of this methodology, the area fraction of lack-of-fusion was measured from the micrographs and computationally generated images, and these values were compared together. To do this, the area fraction of the lack-of-fusion was measured for each micrograph using the ImageJ software, and then the lack-of-fusion images were generated. In the plotting process, after a melt pool shape was filled in white, the boundaries distinguished by black were switched to white boundaries around the melt pools. Similarly, after the area between two horizontal wavy lines was filled in white, the color of these lines changed to white. The increase in the number of layers, , also resulted in an increase in the invisible background of the lack-of-fusion regions. An image, representing the lack-of-fusion extracted from the computationally generated cross-section image of a multitrack and multi-layer sample (see Fig. 11 (a)) generated using this methodology is depicted in Fig. 11 (b). For each cross-validation purpose, 2000 images were generated. White pixels on these images represent the melted area, and the black pixels indicate the unmelted/pore area. To find the total area fractions of lack-of-fusion, the percentage of the black pixels was computed for each lack-of-fusion image. These 2000 values along with their mean were then compared with the observed value . The standard deviation of the percentage values was also computed. The results of some cross-validation rounds are depicted in Fig. 12. The total area fractions of the lack-of-fusion of the corresponding micrographs are presented in these figures.
The difference between and was less than 1%. If the difference is presented in the form of standard deviation , it is expected that the number of standard deviations by which the total area fractions of lack-of-fusion in the micrographs is above or below the mean of the total area fractions of lack-of-fusion in the computationally generated images. This is called standard score and is expressed as [61]: The absolute value of the standard score is expected to be less than 2.576, which is the critical value for a 99% confidence interval [61]. The graph of standard scores versus the difference for all cross-validation rounds was illuminated in Fig. 13. Two conditions were distinguished in Fig. 13, namely the region of acceptable deviation from the mean of the total area fraction of lack-of-fusion and the  region of the acceptable standardized value of area fraction of lack-of-fusion. Circle marks represent the response in the desired region for both conditions, square marks represent that only one condition is met, and the triangle mark indicates that no condition is met. It was found that most of the marks are in the desired rectangular region as indicated in Fig. 13. Only one case did not meet the conditions; however, its standard score was very close to the threshold. By comparing the micrographs recorded experimentally and their corresponding computationally generated images Fig. 13 The standard score of area fraction of lack-of-fusion ( Fig. 10), the total area fractions of lack-of-fusion in both experiment and computationally predicted values (Fig. 12), and the standard scores (see Fig. 13), it can be concluded that the computational model proposed in this study can fairly predict the porosity formed during the experiment with different processing parameters.

Summary and future work
In this study, a new data-driven approach was developed to predict lack-of-fusion in a product printed using the LPBF approach. To develop this methodology, 33 cubes were printed at different processing conditions and scanned for training and test purposes. The micrographs were then used to develop a mathematical model for plotting melt pool shape using (the Nelder-Mead optimization method), and a machine learning model (i.e.,e random forest regression) was employed to predict , , , and . The predicted and were then used to adjust the melt pool shape created using the mathematical model, and and were later used along with the Gaussian process to produce noise data used to model the bi-directional melt pool and powder layers. These layers were then added one by one in the order they were printed to form the part. The model linked the processing parameters to the area fraction of lack-of-fusion formed during manufacturing by predicting one possible scenario that could occur in the morphology of a multi-track and multi-layer 3D-printed part. The current model was developed based on the data obtained from SS 316L and is independent of the thermophysical properties of the material used for 3D printing. Due to the independence of the model from the material properties, the model can be extrapolated, and the lack-of-fusion can be predicted for a wide spectrum of materials regardless of their thermophysical properties, if a larger dataset, covering a wide range of materials is employed.
The current model, similar to many other physics-based and data-driven-based models, also has some limitations. One main limitation of the proposed approach is that despite reasonably predicting the area fraction of lack-of-fusion, the model is incapable of accurately identifying the exact location of lack-of-fusion in the 3D-printed part. This means that the model succeeded in the volumetric mapping of lackof-fusion but not spatial mapping. Also, the accuracy of the model depends on the size of the dataset. The more data is provided, the more accurate the model would predict. These efforts are continued to advance the approach taken in the current study to address these shortcomings and will be presented in a future study. In a future study, the proposed model will be further improved to more accurately predict the morphology of a 3D-printed part and predict the spatial mapping of lack-of-fusion in addition to volumetric mapping.

Conclusions
In this study, a computational approach was developed to predict the multi-track and multi-layer morphology of the parts produced using LPBF processes considering laser power, laser speed, hatch space, and layer thickness as the processing parameters. This approach was then tested to predict the morphology of the samples produced at different processing parameters. The results revealed that the total area fraction of lack-of-fusion computationally predicted in this study is in a fairly good agreement with that estimated from the experimentally recorded micrographs. This approach provided us with multi-track and multilayer morphology of 3D-printed parts and reasonably well allowed us to measure the total area fractions of