Use of Detailed Particle Melt Modeling to Calculate Effective Melt Properties for Powders

Free-form fabrication techniques reduce the cost of producing prototypes or small batch parts by creating parts directly from computer-aided design models. Selective laser melting (SLM) is a promising free-form fabrication process as it works with a wide variety of materials. SLM produces a solid object by selectively fusing successive layers of melt powder. A thin layer of powder is deposited on top of a piston. The surface of the powder is then scanned by a laser with a modulated power, fusing the powder to itself and the layer below where the cross section is intended to be solid and leaving it loose where it is not. When the scan of the layer is complete, a new layer of powder is deposited on top and the process repeats. After the build is complete, the loose powder is removed, leaving the final part [1]. SLM processing parameters (laser power and speed, scan pattern, preheat temperature, etc.) have a strong influence on the quality of the produced part. However, it is often difficult to determine the optimal processing parameters for a given material and geometry. Thus, experimentation and testing are often required when using new materials or geometries to determine the parameters needed to produce an acceptable part [2].

Accurate computational models of the SLM process have the potential reduce the experimentation and testing required in SLM. Bed-scale continuum models, in which the powdered material is treated as a continuous medium as opposed to a collection of individual particles, are particularly promising in their ability to handle large domains without incurring prohibitive computational expense. These models use effective material properties to capture the powdered nature of the medium. Bed-scale models typically use domain sizes on the order of tens of centimeters and simulation times on the order of minutes or hours.

H_p is the enthalpy of the powder, k_p is the thermal conductivity of the powder, f is the laser heat source, c_s is the specific heat of the solid powder, c_l is the specific heat of the liquid metal, L is the latent heat of the powder, ρ is the density of the powder, and g_p(T) is the metal temperature–melt fraction function. g(T) relations are typically calculated for bulk materials by measuring the temperature–enthalpy relationship for the material, H(T). Sharp changes in the slope of the H(T) curve are used to identify the temperatures at which phase change begins and ends [4], and linear relations used to estimate the melt fraction when the temperature is between the two [3]. This method assumes a homogeneous material and uniform melting within a sample. As powdered materials are not homogeneous, finite rate transport effects are critical in determining the rate of powder melting and thus the g_p(T) relation.

By simulating the melting of a powder at the particle scale, the g_p(T) relation can be calculated and applied in a bed-scale continuum model. Additionally, by examining the particle-scale behavior, uncertainties in melt fraction calculations due to particle effects like packing uncertainty can be estimated and propagated through to bed-scale models. This propagation of uncertainty from the particle scale to the bed scale is critical for making meaningful comparisons with experimental data.

Particle-scale models generally use domain sizes on the order of 100 μm and simulation times of no more than a few milliseconds, whereas bed-scale models employ length scales of O(100 mm) or more, and timescales of minutes or hours. Figure 1 shows the difference between the bed scale and the particle scale models. A number of groups have modeled the SLM process on the particle scale. Körner et al. [5] developed a two-dimensional lattice Boltzmann model of the melting of metal powder in SLM. Gürtler et al. [6] used a three-dimensional volume of fluid (VOF) method to model the melting process capable of simulating multiple laser passes. Khairallah and Anderson [7] developed a high-resolution SLM model considering a number of phenomena. However, to our knowledge, no groups have used the results from a particle-scale SLM simulation to calculate a temperature–melt fraction relationship for use in a bed-scale model, and thus infuse the bed-scale model with finite rate transport physics from the particle scale.

Powder bed melt behavior is modeled using a hybrid continuum-discrete approach. The initial positions of the powder particles before the laser is applied are determined using the discrete element method (DEM), which explicitly tracks the locations of particles in a domain under the influence of gravity and interparticle forces. The laser, heat transfer, and flow of the melted particles are then modeled by solving the governing differential equations on a continuum finite volume mesh. The length scale of the domain used is 200 μm × 200 μm × 100 μm and timescales on the order of milliseconds considered.

In order to calculate the temperature–melt fraction relationship, a powder bed is first generated using the discrete element method implemented in the open source software MFiX. In the DEM, particles are modeled as spheres, each with a position, radius, and velocity [8,9]. Particle packings are created by inserting a chosen number of particles in a domain and allowing them to interact with other particles and respond to gravitational forces. Particles interact with each other using a spring-dashpot model in which contact forces are generated based on the degree of overlap a particle has with its neighbors (described in detail by Garg et al. [8]). For the purposes of this work, the MFiX particle–particle interaction model is used simply as a means to generate a random packing of particles.

Once the particles settle, their positions and properties are used as an input for the melting model. The particle packing structure is placed on top of a finite volume mesh. For each mesh cell, the volume of overlap with each of the particles is calculated. The total volume of intersection with all the particles determines the volume faction, β, of solid material in that cell. Cells that are completely contained within a particle have solid volume fractions of 1.0. Cells not overlapping with particles at all have solid volume fractions of 0.0 and cells partially overlapped with particles have solid volume fractions between 0.0 and 1.0. In this way, the DEM representation of the particle bed is converted to a mesh representation that serves as an initial condition for the melting model. Figures 2 and 3 show a packing structure in both the DEM and mesh representations.

Here, γ is the volume fraction of the fluid, defined as 1.0−β, u_m is the mixture velocity, ρ_m is the mixture density, S is the source term due to melting, P is the pressure, τ is the fluid stress tensor, g is the gravitational acceleration, and F_s is the surface tension force. The $ε(1−γ)2/γ3um$ term is a momentum sink that drives the velocity of the mixture to zero in cells that are fully solid.

σ is the surface tension coefficient between the air and liquid metal. σ is calculated per cell based on the temperature, allowing for the Marangoni effect due to gradients in surface tension to be simulated.

Equations (3)–(5) are solved sequentially at each time-step using the OpenFOAM finite volume solvers [12] to evolve the mixture velocity and liquid volume fraction fields. The mixture properties required for the equations (i.e., ρ_m) are calculated by volume averaging using the liquid volume fraction field and the properties of the pure materials: ρ_m = αρ_l + (1.0 – α)ρ_g, where ρ_l and ρ_g are the densities of the pure liquid and gas, respectively.

A similar relationship is used for the enthalpy due to vaporization of the solid material.

Note that Eq. (8) is similar to Eq. (2) for the volume-averaged enthalpy of the powder. Here, however, g_metal(T) is a property of the bulk metal, not a volume-averaged effective property of a powder. Unlike for powders, the temperature–melt fraction relationship for bulk metals can be accurately approximated using a step function for pure elements or a lever rule between a solidus and liquidus temperature for alloys [3].

The energy equation is solved iteratively using the enthalpy method outlined by Swaminathan and Voller [3] in the OpenFOAM finite volume solver until the temperature and enthalpy fields are consistent. Finally, those cells containing metal whose temperatures cross the melting point of the metal are converted to solid or liquid, as appropriate, by explicitly updating the solid and fluid volume fraction fields (γ and β) and generating a source term for the VOF equation to update α, the liquid volume fraction field. For cells undergoing vaporization, sink terms are generated for the VOF and energy equations to remove the vaporized material and associated enthalpy. However, as vaporization is assumed to occur quickly in comparison to the melt flow dynamics, the location of vaporized material is not tracked and is assumed to immediately exit the domain.

$I(x,s)$ is the radiation intensity field, a_m is the mixture absorptivity, and σ is the Stefan–Boltzmann constant. The mixture absorptivity is calculated using the volume fraction fields: $am=(β+γα)ametal+γ(1.0−α)agas$⁠. This formulation of the RTE neglects radiation transport due to scattering. Gusarov and Smurov [13] found that scattering in SLM powder beds had two primary effects. The first is the reduction of laser power delivered to the bed due to reflection out of the domain. This effect is accounted for in the current model by including the emissivity of the material in the incident laser intensity in Eq. (12). Thus, the reflected energy is removed from the laser source before it enters the domain. The second effect is the broadening of the energy absorption region outside the radius of the laser due to scattering in the plane of the powder bed. This effect is not accounted for in the current work. However, Gusarov and Smurov [13] find that lateral scattering is most pronounced in highly reflective powders and so should be minimal for materials considered here.

G is the sum of the irradiation field over all of the solid angles.

ε is the emissivity of the metal, P is the laser power, ω is the characteristic radius, and r is the distance between the center of the laser and a given point on the boundary. The laser intensity is integrated over each cell on the boundary and applied as a fixed-value boundary condition in the RTE. As the laser beam is taken to be entering straight downward into the domain, the boundary condition is applied only to those solid angles for which $s·n$ is approximately −1. This causes the laser intensity to propagate down into the domain.

Uncertainty in the model is due to two main sources: input uncertainty and packing uncertainty. The input uncertainty is due to the uncertainties in the material and laser properties passed to the model as inputs. Packing uncertainty is due to the randomness inherent in the generation of the powder packing structures. Since the initial positions of the DEM particles in the domain are random, the resulting packing structure of the settled particles is also random.

In order to assess the impact of input uncertainties, the generalized polynomial chaos framework [15,16] is employed as a way of representing the stochastic relationship between the inputs and outputs of the model. The inputs are represented as probability density functions (PDFs) and the goal is to determine the probability density function of the output (melt track shape).

X is the quantity of interest, θ is a vector of model inputs, a_j are the coefficients of the expansion, Φ_j are the basis polynomials, and ζ is a vector of random variables representing the input probability distributions. Φ_j are a tensor product basis of Legendre polynomials for uniform random inputs and Hermite polynomials for Gaussian inputs [15].

A stochastic collocation technique [17,18] is used in which the model is solved deterministically at selected collocation points, sparsely distributed in the space of possible input parameters [19]. These collocation points are then used in interpolation schemes to reconstruct the coefficients of the polynomial expansion. In this paper, a Smolyak sparse grid is employed [19] with a second-order polynomial expansion for all output variables in terms of the input random variables.

Once the coefficients of the expansion are known, the resulting response surface may be sampled by drawing samples from the input distributions to predict a probability density function of the output. The standard deviation of the output is a measure of the overall uncertainty in the model prediction resulting from uncertainties in the inputs.

The inputs to the computational model are the laser properties (power, beam diameter, speed) and the material properties (emissivity, specific heat, thermal conductivity, density, surface tension, absorptivity, viscosity). A preliminary sensitivity analysis reveals that the system dynamics are dominated by the laser heat source. Therefore, variations in parameters related to the laser have the most impact on the output. These are laser power, beam diameter, speed, and emissivity (as emissivity scales the overall laser power). As laser speed and power can be accurately controlled, the largest input variations occur with beam diameter and emissivity. Thus, these are taken to be the uncertain inputs for all cases. Variations in other material properties are found to have negligible impact on the output of the model, and so are set to single, deterministic, values.

Packing uncertainties are estimated by generating ten different random powder bed configurations and using the average values of the input parameters to run melting simulations on them. The resulting standard deviation of the melt track shape is then a measure of the uncertainty in the model prediction due to the randomness of the powder bed. To estimate the overall uncertainty in the model prediction for cases involving powder, the standard deviation due to packing uncertainty is added to the standard deviation due to input uncertainty in quadrature. Assuming a Gaussian distribution, the resulting overall model uncertainty can then be computed.

In order to assess the ability of the melting model to accurately predict SLM melt pool geometries, simulation results are compared to experimental results for 316 L stainless steel and Ti-6Al-4V for both flat plates and powders.

For all simulations, a domain size of 200 μm × 200 μm × 100 μm is used consisting of 50 × 50 × 25 mesh cells, for a total of 62,500 cells. In order to assess mesh independence, a simulation is run for 27 μm average diameter 316 stainless steel particles for one packing structure under a 200 W, 54 μm beam diameter laser on a 50 × 50 × 25 mesh and a 100 × 100 × 50 mesh. Material properties used are given in Table 1. Calculated melt pool heights and widths differed by less than 1% between the two meshes and calculated maximum temperatures differed by less than 5%. As these uncertainties are smaller than other sources of uncertainty in the simulation, the 50 × 50 × 25 mesh is considered sufficient resolution. Also, as 54 μm is the smallest beam diameter used and 27 μm the smallest average particle diameter, this resolution is deemed sufficient for all cases.

Simulations are run for enough time to allow the laser to move in the x-direction from an initial x-position of 25 μm to a final x-position of 175 μm. Thus, simulation times are given by 150 μm/v_l, where v_l is the laser speed. Time stepping is controlled dynamically. At each time-step, the Courant number is calculated at the face of each mesh cell using the fluid velocity field and the current time-step size. The time-step size is then adjusted such that the maximum Courant number in the domain does not exceed 0.5. Time step sizes are limited to a maximum of 1 μs to adequately resolve the laser motion when no melt is present and thus the Courant number is zero everywhere. Once melt begins to form, time-step sizes are forced down into the nanosecond range to limit the Courant number.

Gusarov et al. [20] measured the melt pool depth and width on a 316 L stainless steel plate due to laser irradiation. They used an SLM machine to run a laser with different powers and speeds over a flat, stainless steel, plate with no powder. The resulting melt pool geometry was measured by observing the change in microstructure in the melted and resolidified regions.

In order to compare to Gusarov's results, simulations are conducted using a domain consisting of only blocks of air and solid stainless steel, with no super-imposed particles. Material properties for 316 stainless steel are given by Khairallah and Anderson [7] and Khairallah et al. [21] and summarized in Table 1. The emissivity of the plate is given as a uniform probability distribution as opposed to a fixed value, as this is an uncertain input. Surface impurities and finish have a large impact on the emissivity of metals and this range represents the set of possible values for stainless steel.

Table 2 shows the laser parameters used by Gusarov in the two cases simulated. In both cases, the laser beam diameter is given as a uniform probability distribution based on the uncertainties Gusarov measured on the laser used [20].

Using material emissivity and laser beam radius as the uncertain inputs, simulations are conducted to compute the probability distributions of the output melt pool width and depth. Once the PDFs are determined, a 90% confidence interval is calculated. This is done by integrating the PDFs from highest probability density to lowest probability density until a total probability of 90% is reached. The 90% confidence interval thus represents the region in which the true value is most likely to lie given the uncertainty in the inputs.

Figures 4–7 show the calculated PDFs of the melt pool width and depth for the two different laser power and speed combinations along with the 90% confidence interval and the value measured by Gusarov. As can be seen, the measured value falls within the 90% confidence interval for all the cases, indicating that the model is able to account for measurements taken given the uncertainties. However, uncertainties in the simulation are fairly large due mainly to a large uncertainty on the emissivity of the plate. A measurement of the emissivity of the plate taken before an experiment would reduce that uncertainty significantly and thus enable a better comparison.

Khairallah and Anderson [7] measured the height and width of a single melt track caused by a single laser scan line on a single layer of 316 L stainless steel particles on top of a 316 stainless steel plate. In order to compare to Khairallah's results, a domain is initialized consisting of a stainless steel block and a block of air. The air block is filled in with a random packing of stainless steel particles using the method described in Sec. 2.1. A laser scan is then simulated and the resulting melt track height and width calculated. The powder bed parameters used by Khairallah are a layer height of 40 μm and a Gaussian particle size distribution with average particle diameter of 27 μm and a standard deviation of 4.25 μm. The laser parameters are given in Table 3. The material properties of stainless steel are given in Table 1.

Thus, the 0.3–0.6 emissivity uncertainty range transforms to a 0.38–0.59 range, which is no longer uniformly distributed because of the nonlinearity of Eq. (14).

Using powder emissivity and laser diameter as uncertain inputs, simulations are conducted to compute the probability distributions of the output melt pool width and depth, and a mean and standard deviation is calculated for each. However, an additional packing uncertainty is present due to the random nature of the powder bed packing structure. This is estimated by running ten simulations at the mean value of emissivity and beam diameter with different powder bed structures. The standard deviations due to packing and inputs are summed in quadrature to produce an overall standard deviation. The calculated mean and standard deviations due to the different sources of uncertainty are shown in Table 4 along with Khairallah's measured values.

Assuming a Gaussian distribution for the result of the combined PDFs due to input and packing uncertainty, a 90% confidence interval is computed for both width and height and plotted along with Khairallah's measurements in Figs. 8 and 9. As can be seen, the measurements fall within the model 90% confidence interval for both width and height. Also, of note is the large uncertainty due to the powder packing. For both width and height, it is the dominant source of uncertainty. As the exact packing structure of a powder bed within an SLM machine cannot easily be determined, this indicates the presence of a large uncertainty that cannot be eliminated. However, certain factors, such as compaction, may impact the magnitude of the uncertainty.

As powder in an SLM bed may be compacted as a result of the powder spreading process, the effects of compaction on packing structure uncertainty are examined. To do this, powder bed structures are generated using DEM as discussed previously. Then, after a structure is generated, a force is applied to the top boundary, lowering it until it contacts the top of the powder bed, thus applying a pressure to it. Once the powder particles stop moving in response to the applied pressure, it is released. The compacted structure is considered stable if the particles do not move significantly after the pressure is released. As SLM powder beds are not constrained from the top during a build, the maximum applied pressure, which still generates a stable compacted structure, is found and the resulting structure is considered to be the most compacted that could realistically exist in an SLM machine. Ten of these structures are then generated, simulations conducted, and the variation in melt track width and height calculated. It is found that, for these compacted structures, the standard deviation of melt track width and height is reduced by almost half.

Thus, it can be concluded that the details of the powder packing can have a significant effect on the results and uncertainty estimate of a melting simulation. As Khairallah [7] only reports the results of a single experiment, it is impossible to determine the variation in the measured results due to packing and other uncertainties. A measurement of the powder bed emissivity directly prior to the experiment to reduce the uncertainty in emissivity and multiple experimental runs to assess the effect of powder packing would allow a better comparison between simulation and experimental results.

Gong et al. [24] measured the melt pool geometry and melt track geometry for both a Ti-6Al-4V titanium alloy plate and a single track, single layer of powder due to a laser. The experimental and simulation setup for both the plate and powder titanium cases are similar to those used previously for stainless steel plates and powders.

The material properties of Ti-6Al-4V are given in Ref. [4] and summarized in Table 5. The laser parameters used for both the flat plate and powder are given in Table 6. Like Khairallah, Gong gives no uncertainties on the laser beam diameter, so a 10% uncertainty similar to that of Gusarov is assumed.

Figures 10 and 11 show the calculated PDFs of the melt pool width and depth on a titanium flat plate along with the computed 90% confidence intervals and the measurements of Gong et al. [24]. Both measured values fall within the model's 90% confidence interval.

For the powder case, Gong uses titanium powder with a Gaussian size distribution with a mean diameter of 38 μm and standard deviation of 7.9 μm. The layer thickness is 30 μm and the laser parameters the same as those used for the flat plate. Again, Eq. (14) is used to convert the emissivity of the solid titanium to that of a powder.

Table 7 shows the calculated mean and standard deviation of the melt track height and width due to the uncertainty in the emissivity, beam diameter, and the packing structure. Assuming a Gaussian distribution for the combination of the input and packing PDFs, a 90% confidence interval is calculated and plotted in Figs. 12 and 13. As can be seen, both measurements from Gong fall within the 90% confidence interval. Similarly to the stainless steel cases, however, uncertainties in the simulation results are large due mainly to uncertainties in emissivity and, for the powder case, packing. A measurement of emissivity conducted prior to each experiment and multiple experimental runs would need to be conducted to enable to better comparison between model and experiment.

In order to demonstrate how a temperature–melt fraction curve can be calculated from the results of a powder melting simulation, a simulation is run with a 40 μm layer of powder. A 316 L stainless steel powder is assumed. The particle diameter distribution is assumed to be Gaussian, with a mean of 27 μm and a standard deviation of 4.25 μm. The powder is placed on top of a 316 L stainless steel substrate scanned with a 175 W laser at 2 m/s. Figure 14 shows the evolution of the material as it melts and resolidifies.

From the simulation results, a cubic domain element one layer thickness (40 μm) on a side located in the center of the domain is considered. The temperatures of all the cells making up the element are volume averaged together to create an average temperature for the element at each time-step. Similarly, the fraction of solid, unmelted material remaining (⁠$1−β/β0$⁠) is averaged over each cell in the element to create an average melt fraction for the element at each time-step. The resulting curve is shown in Fig. 15 along with the temperature–melt fraction curve of pure solid stainless steel.

For the bulk solid (diamond curve in Fig. 15), the melt fraction–temperature curve is a straight line between the solidus and liquidus temperatures listed in Table 1. However, the powder behaves very differently. Since the powder heats up unevenly due to the laser, small amounts of melt form at very low average temperatures, as only the very top of the material has reached the melt temperature while the material below is significantly cooler. The melting process also takes much longer to occur as the melt forming at the top shields the remaining material from the laser heat. Thus, it can be seen that approximating the temperature–melt fraction curve of a powder with a simple linear relationship between the solidus and liquidus temperatures of the bulk solid material, as done in Refs. [25] and [26], for example, would introduce significant errors in the melting behavior. For ease of use, we employ a linear fit to the average temperature-average melt fraction data from the powder melting calculation. This fit is shown on Fig. 15. Comparing the linear fit to the melting calculation data gives a maximum fitting error committed of 4%.

However, a temperature–melt fraction curve calculated in this way will not necessarily be a constant material property of the powder. In addition to the material properties of bulk stainless steel, the curve may depend on the laser power, scan speed, element orientation relative to the laser path, and powder packing structure. Thus, the effect of variations in all of these parameters will be investigated.

The effect of variations in the distance of the element from the laser path is investigated by shifting the center of the domain element considered by different distances away from the center of the laser scan line. This accounts for the fact that some areas of the powder bed will not fall directly under the center of the laser as it passes by, but will instead be irradiated by a weaker portion of the beam. The difference in calculated melt fraction curves is shown in Fig. 16 for distances of 5–25 μm. As the 4σ laser diameter is 54 μm, elements further away than 25 μm are not melted significantly. As can be seen, approximately the same curve is generated regardless of the element's distance from the laser path, even for elements that do not fully melt. Performing a single linear curve fit on the set of data generated from all of the shifted elements yields a maximum error of 7% when calculating the melt fraction.

The effect of packing structure variation is investigated by running simulations with the same laser and material parameters ten different times with different powder packing structures and then calculating the average temperature–average melt fraction relationship for the same element location in the center of each bed. The difference in calculated curves is shown in Fig. 17. Larger variations are caused by different packing structures than by different element positions. Using a single linear curve fit for all packing structures yields a maximum error of 19% when calculating the melt fraction.

The effect of laser power and speed variation is investigated by running simulations with various different laser powers and speeds on the same powder bed and calculating the temperature–melt fraction curves for elements directly under the laser path. Powers and speeds are chosen to cover a board range of processing parameters that might be used in a single build. The difference in calculated curves is shown in Figs. 18 and 19. Changing laser parameters causes a smaller variation in the melt fraction curves than changing the packing structure, and about the same level of variation as caused by distance away from the laser path. Using a single linear curve fit for all combinations of laser power and speed tested (125, 150, 175, and 200 W powers with 1.5, 1.75, 2, and 2.5 m/s speeds) gives a maximum error of 8% when calculating melt fraction. Combining this with the errors due to distance from the laser path and packing structure in quadrature gives a total estimated uncertainty in melt fraction prediction of 22% when using a single curve to calculate melt fraction regardless of laser parameters, element position, or packing structure.

Although the particle-scale model developed in this work is effective at predicting SLM melt track shapes, it is computationally impractical to use it to model millimeter or centimeter scale part builds. Generating a particle bed in MFiX and running a melting simulation in OpenFOAM requires hours or days of compute time for a 200 μm melt track. In contrast, a bed-scale model can be run for an equivalent melt track length in seconds or minutes. Thus, the goal of this work is to examine whether temperature–melt fraction relationships developed using the particle-scale model can be used to improve the predictions of a faster, reduced order bed-scale SLM model. To do this, two bed-scale calculations are done for 316 stainless steel powder with a 175 W, 2 m/s, 54 μm diameter laser. The results are compared to a particle-scale melting model, which provides “ground truth.” In the first bed-scale calculation, called the “bed-scale bulk material melting model,” bulk material melting is assumed, with melting occurring between the solidus and liquidus temperatures for stainless steel (1648–1673 K) through a linear relationship between the melt fraction and temperature, as done in Ref. [3]. In the second bed-scale calculation, referred to as the “bed-scale powder melting model,” the linear fit to the melt fraction–temperature curve obtained from the particle-scale melt model is used. Intercepts with melt fraction equal to 0 and 1 are found, and correspond to temperatures of 520 and 2740 K, respectively. The effective thermal conductivity of the powder used for the bed-scale simulations is 0.34 W/m K, calculated using the relation developed by Moser et al. [27]. The effective heat capacity used is given by $(1−ε)csteelρsteel+εcairρair$⁠, where c_steel and ρ_steel are given in Table 1, c_air is taken as 1005 J/Kg K, ρ_air is taken as 1 g/cm³, and ε is the porosity, taken to be 0.5. Material properties for the simulations aside from these are given by Table 1.

The bed-scale simulations are performed with a layer thickness of 40 μm on top of a substrate consisting of already consolidated 316 stainless steel using an adaptive mesh. In this way, the area of the powder bed near the laser, where small-scale heat transfer effects are important, is resolved more finely than the substrate, which is only undergoing conduction. The total domain size for both bed-scale simulations is 2.5 cm × 5 cm × 2.5 cm with a mesh size varying from 2.56 mm to 20 μm on a side. This is shown in Fig. 20. The finest mesh size is chosen such that the difference in maximum temperature under the laser calculated in the domain varies by less than 5% with continued refinement and yields a total of 168,453 cells for the domain. The time-step is chosen such that the laser moves at most one laser radius during a time-step, or 13.5 μs for a 54 μm laser diameter moving at 2 m/s. The bed-scale simulations are run for a single scan line in order to compare to the results of the particle-scale melting model.

These bed-scale simulations are compared with a particle-scale melting model. Particle-scale simulations are run with a 40 μm powder layer of 316 stainless steel deposited using the DEM. The particle diameter distribution is Gaussian, with a mean of 27 μm and a standard deviation of 4.25 μm, as used when calculating the temperature–melt fraction relation. The powder is placed on top of a 316 L stainless steel substrate in a 200 μm × 200 μm × 100 μm domain and scanned with a 175 W, 2 m/s, 54 μm diameter laser for 100 μs.

Finally, a 40 μm mesh is superimposed on the results of a particle-scale melting model and the average melt fraction in each cell calculated for comparison to the bed-scale models. This is done for ten different particle packing structure realizations and the results averaged across realizations. The resulting melt fraction predictions are shown in Fig. 21. The bed-scale results have been cropped to show the relevant portion of the domain. The cells marked 1, 2, and 3 in Fig. 21 all have the same centroid locations relative to the laser scan path and the same volume.

To examine the difference between the approaches to melt prediction quantitatively, the calculated melt fraction for the cells or groups of cells marked as 1, 2, and 3 in Fig. 21 are shown in Table 8. Cells in the center of the domain are chosen to avoid the effects of the domain boundaries and the beginning and end of the laser path. The error in Table 8 is defined as the difference between the melt fraction predicted by the particular bed-scale model and that predicted by the ensemble average of the ten different realizations used in the particle-scale melting model.

As can be seen, assuming that the melt properties of the powder are the same of that of the bulk leads to an underestimation of the melt fraction in the cells not directly under the laser path. This is because the bulk melt properties do not properly account for partial melt in the powders since melting can only occur when enough energy has been input to raise all the powder in the cell above the solidus temperature.

On the other hand, using the bed-scale powder melting model allows the uneven melting of the powder to be approximated. The partial melt occurring at low temperatures is representative of the fact that some material closer to the laser may be above the melt temperature even though the average temperature in a cell is not. Using the bed-scale powder melting model, all the bed-scale melt fraction predictions fall well within the expected 0.22 melt fraction uncertainty. Additionally, this uncertainty estimate incorporates powder effects such as packing structure variation and quantifies the uncertainty introduced by these effects. In contrast, bed-scale modeling based solely on bulk material properties provides no clear method by which this can be accomplished. Thus, the powder melting model provides an effective surrogate model at the bed-scale for predicting melt fractions, allowing substantially faster bed-scale models to accurately capture particle-scale effects.

A particle scale melting model of an SLM powder bed is developed and implemented using MFiX-DEM and OpenFOAM. MFiX-DEM is used to generate particle packing structures. These are then converted to a mesh representation and imported into the OpenFOAM finite volume framework in which the model governing equations for fluid flow, heat transfer, and radiation are solved. Results from the model are compared against experimental measurements for stainless steel and titanium plates and powders undergoing laser melting in SLM, and good agreement is obtained accounting for uncertainties in model inputs.

The model is then used to calculate effective, powder, temperature–melt fraction curves for stainless steel. The effects of variations in element position, packing structure, and laser parameters on these curves are examined and quantified. It is shown that using the melting properties calculated this way in a bed-scale model leads to better predictions of partial melting. This technique of using a detailed particle-scale melt model to calculate temperature–melt fraction curves can thus be used to improve the accuracy of melt predictions in SLM bed-scale models as well as provide uncertainty estimates due to powder effects not possible with other methods.

• National Science Foundation (Grant No. CNS-1239243).

a_j =

expansion coefficients

a_m =

mixture absorptivity

c =

heat capacity

c_l =

liquid heat capacity

c_s =

solid heat capacity

DEM =

discrete element method

f(x, y, z, t) =

heat source term

F_s =

surface tension force

g =

gravitational constant

g(T) =

metal temperature-melt fraction relationship

g_p(T) =

powder temperature-melt fraction relationship

G =

irradiation field

H =

enthalpy

H_m =

enthalpy of mixture

H_p =

enthalpy of powder

I =

radiation intensity

k_p =

thermal conductivity of powder

L =

latent heat of fusion

P =

laser power

P =

pressure

PDF =

probability density function

S =

mass source term

SLM =

selective laser melting

T =

temperature

T_ref =

reference temperature

Ti-6Al-4V =

titanium alloy

u_m =

mixture velocity

v_l =

laser speed

X(θ) =

quantity of interest

α =

liquid fluid volume fraction

β =

volume fraction of solid

β₀ =

initial volume fraction of solid

γ =

volume fraction of fluid

ε =

emissivity

ε =

momentum sink constant

ε =

porosity

ε_eff =

effective emissivity

ζ(θ) =

random variable

μ =

mean

ρ =

density

ρ_g =

gas density

ρ_l =

liquid density

ρ_m =

mixture density

σ =

Stefan–Boltzmann constant

σ =

standard deviation

σ =

surface tension coefficient

τ =

fluid stress tensor

ϕ_j =

basis polynomials

ω =

laser characteristic radius