In this work, a discrete element model (DEM) is developed and implemented in the open source flow solver MFiX to simulate the effective thermal conductivity of powder beds for selective laser sintering (SLS) applications, considering scenarios common in SLS such as thin beds, high temperatures, and degrees of powder consolidation. Random particle packing structures of spherical particles are generated and heat transfer between the particles is calculated. A particle–particle contact conduction model, a particle–fluid–particle conduction model, and a view factor radiation model using ray-tracing for calculation of view factors and assuming optically thick particles are used. A nonlinear solver is used to solve for the particle temperatures that drive the net heat transfer to zero for a steady state solution. The effective thermal conductivity is then calculated from the steady state temperature distribution. Results are compared against previously published experimental measurements for powder beds and good agreement is obtained. Results are developed for the impacts of very high temperatures, finite bed depth, consolidation, Young's modulus, emissivity, gas conductivity, and polydispersity on effective thermal conductivity. Emphasis is placed on uncertainty quantification in the predicted thermal conductivity resulting from uncertain inputs. This allows SLS practitioners to control the inputs to which the thermal response of the process is most sensitive.

Introduction

The free-form fabrication techniques reduce the costs of creating prototypes or small batch parts by creating parts directly from computer-aided design models. SLS is a promising free-form fabrication process as it works with a wide variety of materials and does not require the creation of support structures to hold the part during build. SLS produces a solid object by selectively fusing successive layers of 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, the piston holding the part is lowered, 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]. SLS processing parameters (laser power and speed, scan pattern, preheat temperature, etc.) have a strong influence on the quality of the resulting part. However, it is difficult to determine the optimal processing parameters for a given material and geometry. Thus, experimentation is often required when using new materials or geometries to determine the parameters needed to produce an acceptable part. Testing of each part is often needed if the part is required to meet certain quality standards.

Accurate computational models of the SLS process have the potential for providing solutions to the problems of parameter determination and quality assurance. 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 on the same scale as the part being built without incurring prohibitive computational expense. Continuum models describe the heat transfer in the powder bed using the heat equation
(1)

While density can be easily measured and is often provided by the powder manufacturer, specific heat does not change when the solid material is powderized, and laser properties are either known or controllable, bulk properties such as the effective thermal conductivity cannot be so easily inferred. While the thermal conductivity of the solid may be available, it is expected to decrease considerably when the material is in powdered form due to the introduction of contact resistance between the powder particles.

Numerous works have investigated, experimentally, analytically, and computationally the effective thermal conductivity of powdered materials [2]. Masamune and Smith [3], Cheng and Vachon [4], Gusarov et al. [5], and Slavin et al. [6,7] all developed analytical models of conductivity of powder beds. Xue and Barlow [8,9] and Yuan et al. [10] both measured the effective thermal conductivities in Nylon-12 powder beds and developed empirical relations from this for SLS applications. Sih [11] measured the effective thermal conductivities, as well as other bulk powder properties, of a variety of metals for use in SLS. Vargas and McCarthy [12] used a particle dynamics simulation to predict the effective conductivity of granular media. Zhang et al. [13], Tsory et al. [14], Widenfeld et al. [15], and Feng et al. [16] all used DEM to calculate the thermal conductivity of powder beds. However, none of these works consider powder bed scenarios commonly found when performing SLS, specifically finite bed depths, when a thin layer of powder is present on top of a solid, already sintered, surface, and temperatures close to the melting temperature where radiation effects become increasingly important. This work uses DEM to examine the impact of these, as well as other, effects on the effective thermal conductivity of a powder bed.

Modeling Approach

In order to calculate effective thermal conductivity, a powder bed is represented using the DEM of the open source software MFiX as a series of spherical particles, each with a position, radius, emissivity, and solid conductivity [17,18]. As the particles are made up of solid material, solid properties can be used for the particle conductivity. Each particle is modeled as a lumped capacitance control volume with a single, uniform temperature. Temperature gradients within particles are not resolved. Heat can be exchanged between particles via particle–particle conduction (Q˙r), particle–fluid–particle conduction (Q˙pfp), and radiation (Q˙rad). The net heat source in a particle is then given by the sum of all the heat being exchanged with all other particles across all three mechanisms as illustrated in Fig. 1. Heat transfer between the particles and the background gas is neglected as the gas conductivity is an order or magnitude less than the calculated effective thermal conductivities. We have performed scoping calculations which establish that this approximation has negligible impact on our results.

Fig. 1

Fixed temperatures are set at two opposite walls of the domain, as shown in Fig. 2. All particles are given an initial temperature and net heat sources are calculated for each particle. Powell's method [19], with Jacobians calculated by finite differencing, is then used to determine a temperature for each particle such that Q˙r+Q˙pfp+Q˙rad=0, thus putting the particle bed in steady state. Once the steady-state particle temperatures are known, the heat fluxes from the walls are calculated and from these the effective thermal conductivity is determined using Fourier's law.

Fig. 2

Powder Bed Generation.

Particle packings are generated by inserting a chosen number of particles per time at a random location on a boundary with an initial velocity and allowing them to interact with other particles in the domain and respond to gravitational forces [18]. 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. [17]). For the purposes of this work, the MFiX particle–particle interaction model was used only as a means to generate a random packing of particles. Particles are injected into the domain from the top boundary and allowed to fall under the influence of gravity and interact with other particles already in the domain. Once the particles settle, their positions and properties are used as an input for the heat transfer model.

Particle–Particle Conduction.

The particle–particle conduction model accounts for the heat transfer due to particles being in direction contact. The MFiX–DEM particle–particle conduction model uses a modified Batchelor and O'Brien [20] method in which the heat transfer between two particles in contact is proportional to the radius of contact, as shown in Fig. 3 [18].

Fig. 3
Particle–particle conduction
Fig. 3
Particle–particle conduction
Close modal
(2)
(3)

Gusarov et al. [5] and Vargas and McCarthy [12] used a similar model in their efforts. A similar equation applies for a particle in contact with a wall, which is modeled as a flat plane. In that case, li,j is defined as the distance between the particle center and the wall and the contact radius becomes Rc=Ri2li,j2.

When using the MFiX–DEM spring–dashpot model to generate bed packing structures, the size of the time step that the explicit solver can take is governed by the value of the particle spring constant used in the model [17]. Thus, if spring constants were set to values reflective of the actual material Young's modulus, the necessary time step size would be prohibitively small. To counteract this, particles are generally modeled as being softer than they actually are in order to reduce the computational cost. Reducing spring coefficients to 105 N/m is found to have little impact on the resulting packing structures as compared to much stiffer coefficients and results in greatly reduced computational times.

However, when using the heat transfer models, this leads to an overestimation of the overlap between particles and thus an overestimation of particle–particle conduction. Therefore, a correction is applied to Rc by equating the normal contact force calculated by the spring–dashpot model (Fn=knδs) with the contact force calculated using Hertzian contact theory (Fn=(4/3)EeRe1/2δh3/2), where δs is the particle overlap predicted by the spring–dashpot model. The equation is solved for δh, the overlap predicted by Hertzian theory, and which is used to calculate a new radius of contact to use in the heat transfer model [13]
(4)
(5)
(6)

Thus, while softer particles are used for bed generation, true particle material properties are used by the heat transfer model.

In this analysis, the contacts between particles are assumed to be perfectly smooth. Thus, while contact resistance is modeled using the contact radius, the effect of surface roughness, which can reduce the effective contact radius further, is not included. By its nature, the DEM approximates particles as perfectly spherical and is not expected to yield good predictions for nonspherical particles or particles with highly irregular surfaces. As is discussed in Secs. 3.1 and 3.2, the smooth, spherical particle model outlined here yields good results for metal powders, but ceramic powders require calibration.

Particle–Fluid–Particle Conduction.

The particle–fluid–particle model accounts for the heat transfer between the particles through the gas layer around the particles which occurs when two particles are close by, as shown in Fig. 4. This mechanism differs from the conduction of heat from a particle or wall to a background gas and then back to a particle or wall some distance later. The particle–fluid–particle mechanism accounts for heat transfer directly between particles close to each other across a small interstitial gas path. The MFiX–DEM particle–fluid–particle conduction model uses a modified Rong and Horio [21] method

Fig. 4
Particle–fluid–particle conduction
Fig. 4
Particle–fluid–particle conduction
Close modal
(7)
(8)

lcond is the conduction distance between the two particles at a given r value and Rf is the fluid radius, the radius of the circle formed by the overlap of the two particle's fluid lenses. A fluid lens is a layer of fluid surrounding a particle in which particle–fluid–particle heat transfer can occur. By default, MFiX–DEM uses a value of 0.2R for this parameter [18]. Equation (7) contains a singularity at the point of contact between two particles as the conduction distance goes to zero. This is remedied in MFiX–DEM by imposing a minimum conduction distance, which is set to 1 μm by default [18].

Radiation.

In this paper, the heat transfer models of MFiX–DEM are extended by the addition of a view factor radiation model. The view factor model is only valid in the geometric optic limit, so (πDp/λ)1 must be satisfied for it to be used, where Dp is the average particle diameter and λ is the peak wavelength emitted by the particles at the temperature of interest. Additionally, the model assumes all particles to be completely opaque, so it cannot be applied to transparent or partially transparent materials.

The radiation heat transfer between two particles is given by the equation for reradiating surfaces
(9)

A similar equation applies for heat transfer between a particle and a wall due to radiation.

The view factors between the particles and other particles and particles and walls are determined using a Monte Carlo method. For every particle, a number of rays are fired from the particle from a random location on its surface in a random direction and the first particle or wall they intersect with is recorded. The view factors between a particle and any other particle or wall are then given by the ratio of the number of rays fired from the first particle that intersect the second particle or wall to the total number of rays fired from that particle.

The Marsaglia rejection method [22] is used to pick a random point on the surface of a particle from which to fire a ray. This method allows a point on a sphere to be randomly selected without the use of trigonometric functions. Two uniform random numbers, x1 and x2, between −1 and 1, are selected and the sum of their squares computed. If the sum of squares is greater than 1, the set is discarded and two new numbers are selected. If it is less then 1, x, y, and z coordinates of a random point on the sphere can be calculated from the pair using algebraic manipulations [22]
(10)
(11)
(12)

Once a point on a particle is chosen, a coordinate shift is then applied and centered about that point to align the z-direction with the outward facing normal. Then, the direction of the fired ray is selected by choosing a random point on the unit sphere centered at the point chosen on the particle. The Marsaglia rejection method [22] is used again with an additional rejection criterion applied that the chosen point must have a non-negative z-component to ensure that a point on the outward facing hemisphere is chosen. Finally, the vector determined by connecting the chosen point on the surface of the particle with the chosen point on the outward facing unit hemisphere is shifted back to the original coordinate system [23].

Uncertainty Quantification.

Uncertainty in the model prediction is due to three different sources: input uncertainty, bed generation uncertainty, and consolidation uncertainty. The input uncertainty is due to the uncertainties in the material properties passed to the model as inputs. These are gas conductivity, solid particle conductivity, Young's modulus, Poisson's ratio, and solid density. Bed generation uncertainty is due to the inherent randomness in the process used to generate the powder beds. Since the initial positions of the particles when they enter the domain are random, the final resulting powder bed structure once the particles settle is also random. Consolidation uncertainty occurs because the degree of consolidation in a given real powder bed is unknown. Consolidation occurs when particles in a powder bed are rearranged to give a more compact structure and commonly occurs as the result of a pressure being applied to the bed and then released or the bed being shaken. Consolidation results in a decrease in the porosity of a powder bed and thus changes the effective thermal conductivity.

In order to assess the impact of input uncertainties, the generalized polynomial chaos framework [24,25] is employed as a way of representing the stochastic relationship between the inputs and outputs of the model. The model inputs are represented as probability density functions and the goal is to determine the probability density function of the output (effective thermal conductivity). The output is expressed as a polynomial expansion in the input random variables and a stochastic collocation technique [26,27] is used in which the model is solved deterministically at selected collocation points, sparsely distributed in the multidimensional random parameter space [28]. These collocation points are then used in interpolation schemes to reconstruct coefficients of a polynomial expansion. In this paper, a Smolyak sparse grid is employed [28] with a second-order polynomial expansion. 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.

Bed generation uncertainties are estimated by generating multiple different random powder bed configurations and averaging the calculated effective thermal conductivity across all of them. The resulting standard deviation is then a measure of the uncertainty in the model prediction due to the randomness of the powder bed.

Consolidation uncertainties are estimated by lowering the top boundary of a powder bed to the point where it overlaps with the top layer of particles by a small amount (on the order of the particle radius). The resulting forces then cause the powder bed to rearrange itself into a more compact structure. Finally, the top boundary is raised again to relieve any stresses that could not be relieved by particle rearrangement. The effective thermal conductivity is then calculated for the unconsolidated powder bed as well as for four different consolidated beds produced by lowering the top boundary by different amounts. Results are averaged across all five beds and a standard deviation calculated. This standard deviation is an estimate of the uncertainty due to the unknown degree of bed consolidation, or, equivalently, the unknown porosity.

Results and Discussion

Model Validation for Metal Powders.

In order to assess the impact of finite bed depth, polydispersity, and high temperatures on the effective thermal conductivity of powder beds, the model is first validated against data available in the literature for metal particles. Cheng and Vachon [4] obtained the experimental data for both steel of an unknown alloy and lead particles at room temperature, Bala et al. [29] measured the effective thermal conductivity of copper particles at room temperature, and Widenfeld et al. [15] presented the data for steel particles of an unknown alloy at approximately 30 °C. Monodisperse particle beds are generated for each material in a cubic domain 10D on each side. This is determined to be large enough that further increasing the size of the domain has negligible effect on the results. Table 1 shows the data from these three sources compared against the predictions of the current model. kmodel is the effective thermal conductivity predicted by the model, σmodel is the standard deviation of the model prediction, and kexp is the reported experimental effective thermal conductivity.

Table 1

Comparison of model predictions with experimental data

MaterialDp (mm)kmodel (W/m K)σmodel (W/m K)kexp (W/m K)
Copper0.250.6270.0510.652 [29]
Copper0.150.5760.0610.546 [29]
Lead1.60.4570.0220.418 [4]
Steel1.00.3330.0220.34 [15]
Steel3.20.3970.0220.4–0.6 [4]
MaterialDp (mm)kmodel (W/m K)σmodel (W/m K)kexp (W/m K)
Copper0.250.6270.0510.652 [29]
Copper0.150.5760.0610.546 [29]
Lead1.60.4570.0220.418 [4]
Steel1.00.3330.0220.34 [15]
Steel3.20.3970.0220.4–0.6 [4]

Table 2 shows the material properties used for the three different metals and Table 3 shows how many particles are in the beds generated for each of the five cases along with the corresponding average porosity across the different levels of consolidation. The standard deviation of the porosity is included to show the degree of uncertainty in the porosity as a result of consolidation. In each case, particles are added to completely fill the domain. Young's modulus, particle conductivity, and gas conductivity are represented as uniform probability density functions for the purposes of calculating input uncertainty. The uniform probability density functions are chosen to capture the high degree of uncertainty in the steel alloys being used and the exact temperature experiments are run at.

Table 2

Model input parameters

Materialks (W/m K)kg (W/m K)E (GPa)νρ (g/cm3)
Steel12.11–45.00.025–0.027190–2060.287.7
Copper393–4090.025–0.027110–1280.348.96
Lead34.24–36.360.025–0.02713.8–160.4411.36
Materialks (W/m K)kg (W/m K)E (GPa)νρ (g/cm3)
Steel12.11–45.00.025–0.027190–2060.287.7
Copper393–4090.025–0.027110–1280.348.96
Lead34.24–36.360.025–0.02713.8–160.4411.36
Table 3

Particle count and porosity

MaterialNo. of particlesAvg. porosityStd. dev.
Steel (1 mm)10390.410.01
Steel (3.2 mm)10900.4180.009
Copper (0.25 mm)10950.4150.009
Copper (0.15 mm)10910.4160.008
Lead (1.6 mm)10930.4160.009
MaterialNo. of particlesAvg. porosityStd. dev.
Steel (1 mm)10390.410.01
Steel (3.2 mm)10900.4180.009
Copper (0.25 mm)10950.4150.009
Copper (0.15 mm)10910.4160.008
Lead (1.6 mm)10930.4160.009

The other inputs (density and Poisson's ratio) are either well known or determined to have little impact on the calculated effective thermal conductivity and thus are assumed to be fixed values.

Table 4 shows the input and consolidation uncertainties calculated for each material. Input uncertainties represent the standard deviation in the output due to all uncertain inputs. Bed generation uncertainty is calculated for 1 mm steel particles averaged across eight powder beds, and the resulting standard deviation is 0.002 W/m K. As this is an order of magnitude less than the largest source of uncertainty, bed generation uncertainty is neglected for the remaining cases. The input and consolidation uncertainties are added in quadrature to yield the values shown in Table 1 for the total uncertainty. As can be seen from Table 1, the model shows good agreement with the experimental results as all but one case is within one standard deviation of the predicted value and all cases are within two.

Table 4

Input and consolidation uncertainties

MaterialDp (mm)σinput (W/m K)σcons (W/m K)
Copper0.250.0090.05
Copper0.150.0080.06
Lead1.60.0080.02
Steel1.00.010.02
Steel3.20.020.01
MaterialDp (mm)σinput (W/m K)σcons (W/m K)
Copper0.250.0090.05
Copper0.150.0080.06
Lead1.60.0080.02
Steel1.00.010.02
Steel3.20.020.01

Figure 5 shows the variation of the particle bed temperature between the two fixed temperature walls for 1 mm steel particles. Since particles are modeled using a lumped capacity method, no temperature gradient is maintained in the particles. Thus, the temperature variation in the bed is calculated using a binning method. The domain is divided into ten bins based on location. All particles whose centers lie within a given bin have their temperatures averaged together to give an average temperature for the bin. These average temperatures are then plotted against bin location to give a temperature profile. As can be seen, this bin-averaged temperature profile is approximately linear.

Fig. 5
Particle bed temperature profile
Fig. 5
Particle bed temperature profile
Close modal

Model Validation for Ceramic Powders.

Model predictions are also validated for ceramic particles against the data of Slavin et al. [7] for packed alumina particles in helium. We find that if the procedure described for metallic particles is used for ceramic particles, predicted effective thermal conductivity deviates from the measurements by 30%. This is likely due to either the actual particles not being near-spherical, as is assumed in the model, the importance of surface roughness, or the assumed Hertzian contact mechanics not being a good approximation of the actual contact behavior in ceramic beds. Therefore, calibration parameters are introduced. The first parameter allows the contact area between two particles to be adjusted by scaling Eq. (4). The new contact radius is then given by the below equation:
(13)

The second parameter is the minimum conduction distance, lmin, which impacts Eqs. (7) and (8). It removes the singularity that occurs when particle–fluid–particle conduction is calculated. For metals, the default MFiX value of 1 μm works well, but for ceramics it is used as a calibration parameter.

Taken together, Ac and lmin allow the relative importance of particle–particle and particle–fluid–particle conduction to be adjusted. These parameters are calibrated using the data of Slavin et al. [7] for temperatures up to 550 K by means of a Levenberg–Marquardt least squares algorithm [30]. The calculated parameter values are Ac=5.981 and lmin=2.419×105m. Since Ac is larger than 1 and lmin is larger than 1 μm, this means that particle–particle contact conduction plays a larger role in the heat transfer than is predicted by the default MFiX–DEM model parameters.

The calibrated model is then validated against the alumina data for the remaining temperatures measured by Slavin et al. [7], 550–750 K. At the first temperature, 337 K, model input uncertainty is estimated using the same method as for metals, with the material properties of alumina given in Table 5. The calculated standard deviation due to input uncertainty is 0.02 W/m K. Averaging across different random powder bed structures produced a bed generation uncertainty of 0.02 W/m K and averaging across consolidated and unconsolidated beds yielded an average porosity of 0.403 with a standard deviation of 0.008 and a consolidation uncertainty of 0.5 W/m K. Therefore, since consolidation uncertainties are an order of magnitude greater than other uncertainties, only consolidation uncertainty is considered at all other temperatures. The predicted conductivities and the experimental results are shown in Fig. 6. As can be seen, the measurements fall within the model uncertainties for all data points. Thus, the model can be used to predict the thermal conductivity of ceramic powder beds as well as metal ones, but calibration is necessary.

Fig. 6
Comparison of model predictions with experimental results of Slavin et al. [7]
Fig. 6
Comparison of model predictions with experimental results of Slavin et al. [7]
Close modal
Table 5

Alumina material properties

PropertyValue
ks (W/m K)22.7–26.0
E (GPa)344.8–409
kg (W/m K)0.160–0.166
v0.21
ρ (g/cm3)3.95
PropertyValue
ks (W/m K)22.7–26.0
E (GPa)344.8–409
kg (W/m K)0.160–0.166
v0.21
ρ (g/cm3)3.95

Model Predictions.

The model is then used to asses the effect of several parameters on the thermal conductivity of metal powder beds. Monodisperse, 1 mm steel particles with a bed height of 10 particle diameters, a bed temperature of 300 K ((σT3Dp/kg)=0.06), and properties given in Table 2 are used as an example material in all the following investigations unless otherwise indicated. Uncertain properties are taken to be exactly at the mean value and uncertainties in the outputs are thus due entirely to consolidation uncertainty.

Effect of Bed Temperature.

First, the effect of bed temperature is examined. In SLS builds using metal powders, temperatures greater than 1000 K are common and radiation heat transfer is expected to play a large role in the effective conductivity. In order to assess this, the bed temperature is increased while all other inputs except the gas conductivity are held constant. As gas conductivity varies with temperature and most SLS builds are done in air, the gas conductivity is set to be the conductivity of air at the operating temperature. The relation kg(W/mK)=6.566×1012T33.386×108T2+9.426×105T+7.505×104, obtained from curve fitting a cubic polynomial to the thermal conductivity of air for temperatures ranging from 175 to 1900 K, is used. The effective thermal conductivity is then calculated for each case. At each temperature, the model is also run with the radiation module switched off for comparison. The results are shown in Fig. 7. Radiation begins to play a significant role in the heat transfer for ratios of σT3Dp/kg greater than 0.5. For the steel particles used in the model, this corresponds to a temperature of about 1000 K. For particles on the order 0.1 mm, this will be around 1300 K in air and for 0.01 mm particles this will be almost 3000 K. Thus, radiation heat transfer will need to be considered in the effective thermal conductivity when performing SLS with metals or other similarly high melting point materials with particle sizes on the order 0.1 mm or above. For smaller particles, the material will likely melt before radiation heat transfer becomes significant. Additionally, even at higher temperatures, the conductivity of the powder is one to two orders of magnitude smaller than that of the solid, indicating that conduction is limited by the contact between the particles and may thus be significantly increased because of the parallel path provided by radiation at higher temperatures.

Fig. 7
Variation of effective thermal conductivity with temperature
Fig. 7
Variation of effective thermal conductivity with temperature
Close modal

Effect of Gas Conductivity.

The model is then similarly used to examine the effect of interstitial gas conductivity. Although some SLS builds may be performed in a vacuum or near vacuum, most commonly they are conducted in air or nitrogen at atmospheric pressure. Thus, the gas conductivity contributes to the overall effective thermal conductivity through the particle–fluid–particle (Q˙pfp) conduction pathway. Particle–gas interactions are neglected. All other model inputs are held constant and the gas thermal conductivity is increased across a broad range of values and the effective thermal conductivity calculated. The thermal conductivity of the beds with the gas conductivity switched off is also calculated for comparison. The results are shown in Fig. 8. As can be seen, the conductivity of the gas begins to have a significant effect around a kg/ks ratio of 0.001 and becomes dominant as the ratio goes to 0.1. As most metals have a thermal conductivity on the order of 10–100 W/m K and air and nitrogen have conductivity on the order to 0.01 W/m K, gas conductivity will have to be considered when determining the effective thermal conductivity of SLS powder beds for some cases.

Fig. 8
Variation of effective thermal conductivity with gas conductivity
Fig. 8
Variation of effective thermal conductivity with gas conductivity
Close modal

Effect of Young's Modulus.

Next, the effect of the material's Young's modulus is considered by varying the Young's modulus while holding all other model inputs constant. The results are shown in Fig. 9. The effective thermal conductivity is relatively insensitive to changes in Young's modulus, even when varied by several orders of magnitude. As most metals have a Young's modulus on the order of 100 GPa, variations in Young's modulus between materials would generally not have to be taken into account when determining the effective thermal conductivity of an SLS powder bed.

Fig. 9
Variation of effective thermal conductivity with Young's modulus
Fig. 9
Variation of effective thermal conductivity with Young's modulus
Close modal

Effect of Emissivity.

The effect of emissivity is also considered by varying the emissivity between 0.1 and 1 while holding all other inputs constant. As emissivity only plays a role in the radiation calculation, the model is run at a temperature of 1000 K, which corresponds to an σT3Dp/kg of about 0.7, so that the impact of emissivity variations can be seen. The results are shown in Fig. 10. As can be seen, the effective thermal conductivity is also not particularly sensitive to emissivity, even at temperatures where radiation heat transfer is significant. This can be explained by examining Eq. (9). Since view factors between particles that can see each other are generally in the order of 0.1, varying emissivity between 0.1 and 1.0 causes the radiative flux to vary at most by a factor of 3. As radiative flux only accounts for about 30% of the total flux at this temperature, the overall sensitivity to emissivity is low.

Fig. 10
Variation of effective thermal conductivity with emissivity
Fig. 10
Variation of effective thermal conductivity with emissivity
Close modal

Effect of Bed Height.

The effect of bed height is also examined as the SLS machines use very thin layers during a part build and thus it is possible for a thin single powder layer to be spread on top of an already formed solid surface, resulting in a thin powder bed. In the model, this is simulated by discarding all particles which extend above a certain height in each powder bed geometry and then by calculating the effective thermal conductivity while holding all other parameters constant. The results are shown in Fig. 11. For very thin beds, on the order of a particle diameter, the effective thermal conductivity can be almost cut in half. However, this effect diminishes quickly and by bed heights of 6–8 particle diameters the conductivity is within the uncertainty of the infinitely deep bed value. Depending on the ratio of SLS layer thickness to powder particle diameter, bed height may play a role in the thermal conductivity of the powder for some applications.

Fig. 11
Variation of effective thermal conductivity with bed height
Fig. 11
Variation of effective thermal conductivity with bed height
Close modal

Effect of Polydispersity.

Finally, the effect of powder polydispersity is examined. In all previous simulations, the powder particles were assumed to have a uniform diameter. For real powder particles, however, this will generally not be the case and there will be some variation of particle diameter within the powder bed. In order to simulate the effect of powder polydispersity, packings are generated with particles of different sizes. When a particle is added to the domain, it is added at a random location with a random diameter. Although the model can be used to simulate other distributions, Gaussian particle size distributions are assumed as only one additional input parameter (standard deviation) is introduced. The Gaussian distribution is split into nine bins, each with a representative diameter and a fraction of the total number of particles determined by the Gaussian. Particles entering the domain are then assigned one of the nine possible diameters based on the probability of each. Table 6 shows the average porosity and standard deviation due to consolidation on the polydisperse powder beds.

Table 6

Porosity of polydisperse beds

σAvg. porosityStd. dev.
0.05Dp0.4230.007
0.125Dp0.4190.007
0.25Dp0.4210.007
0.375Dp0.4150.007
0.5Dp0.4280.007
σAvg. porosityStd. dev.
0.05Dp0.4230.007
0.125Dp0.4190.007
0.25Dp0.4210.007
0.375Dp0.4150.007
0.5Dp0.4280.007

Once the polydisperse beds are generated, the effective thermal conductivity of each is calculated using the same method as for monodisperse beds. The effect of polydispersity on the calculated conductivity can be seen in Fig. 12, with the x-axis being the standard deviation of the particle diameter distribution normalized by the average particle diameter. There is about a 30% increase in the effective conductivity across the investigated range of standard deviations. This is likely due to the fact that the smaller particles in the distribution are able to fill in the gaps in the packing structure, leading to more contact between particles. Thus, polydispersity can impact the conductivity of a powder bed, even for small deviations, although the impacts are not as large as those caused by gas conductivity or temperature. The uncertainties due to consolidation also increase with increasing standard deviation as the variation in particle size adds extra degrees-of-freedom to the powder bed layout.

Fig. 12
Variation of effective thermal conductivity with particle size standard deviation
Fig. 12
Variation of effective thermal conductivity with particle size standard deviation
Close modal

Discussion and Empirical Relation.

Of the parameters investigated, effective thermal conductivity is most dependent on temperature and gas conductivity. At low temperatures and low gas conductivities, conduction is primarily through the particle–particle path and thus is contact controlled. However, as temperature increases or gas conductivity increases (either as a result of temperature increase or due to the use of a higher conducting gas), parallel paths are provided for heat transfer through radiation and particle–fluid–particle conduction. Temperature and gas conductivity, therefore, have the largest impact of the effective conductivity as they control the availability of these parallel pathways. Thus, an empirical correlation is developed to estimate the effective conductivity based on these quantities. A second-order polynomial expansion is developed and the coefficients determined using collocation points calculated from a Smolyak sparse grid [28]. Temperature and gas conductivity are nondimensionalized as θ=(σT3Dp/ks) and κ=(kg/ks). The correlation is then given below:
(14)

The correlation is developed for a θ range of 0.00035–0.023 (representing 500 K–2000 K for 1 mm steel particles) and a κ range of 0.00125–0.006. Material properties for monodisperse 1 mm steel particles were used for all other inputs. However, since variations in emissivity between 0.1 and 1.0 only lead to approximately an 8% change in conductivity even at high temperatures and variations in Young's modulus across several orders of magnitude only produce a 2% change in conductivity, the correlation should be applicable to a wide range of metals. Polydispersity up to 0.5Dp can lead to a change in conductivity of 30% and bed heights less than 6Dp lead to 100% changes. Therefore, for highly polydisperse powders and very thin beds this correlation cannot be used and the full model must be employed directly. Degree of bed consolidation is found to be the dominate uncertainty in calculating the effective thermal conductivity, leading to a model uncertainty of up to 11%.

Conclusions

A particle model of an SLS powder bed is developed and implemented using MFiX–DEM. The existing MFiX–DEM heat transfer modules are enhanced by the addition of a ray-tracing technique for particle–particle and particle–wall radiation. The model is used to calculate thermal conductivity of metal powder beds. Quantification of uncertainties in model inputs as well as bed layout and consolidation is included. The DEM model necessarily assumes beds made up of smooth, spherical particles. For metal powder beds, the assumption is good and calculated results are in good agreement with experimental data without the need for calibration. For ceramics, this is not the case, and calibration is necessary, but good agreement with experimental results is achieved after calibration of contact conduction and particle–fluid–particle models. The variation of the effective thermal conductivity with temperature, Young's modulus, gas conductivity, emissivity, finite bed depth, and polydispersity is calculated. As temperature and gas conductivity are found to have the most impact, a correlation is proposed to calculate the effective thermal conductivity from these two quantities for metal beds in air. These results have the potential to aid in the quantification of uncertainties in continuum SLS models by providing the methods to calculate effective conductivity in a variety of scenarios important in SLS where other methods of estimating conductivity cannot be applied. However, care should be taken when applying the model to beds known to consist of highly nonspherical or very rough particles as results may not be accurate and calibration with experimental measurements may be required.

Acknowledgment

We would like to acknowledge Oak Ridge National Laboratory for their support and the contributions of Aaron Morris. The support by the National Science Foundation under Grant No. CNS-1239243 is gratefully acknowledged.

Nomenclature

Ac =

contact area correction parameter

Ai =

particle surface area

cp =

specific heat

CAD =

computer-aided design

Dp =

particle diameter

DEM =

discrete element model

E =

Young's modulus

f(x,y,z,t) =

laser heat source

Fij =

view factor from particle i to j

g =

gravitational constant 9.81 m/s2

k =

thermal conductivity

kg =

gas thermal conductivity

ks =

solid thermal conductivity

kN =

spring constant for DEM spring–dashpot model

lcond =

conduction distance for particle–fluid–particle conduction

li,j =

distance between the centers of particles i and j

lmin =

minimum conduction distance for particle–fluid–particle conduction

Q˙pfp =

particle–fluid–particle heat flow

Q˙r =

particle–particle heat flow

Q˙rad =

radiation heat flow

Rc =

particle–particle contact radius

Rc,e =

correct particle–particle contact radius

Re =

effective radius

Rf =

fluid radius for particle–fluid–particle conduction

Ri =

radius of particle i

SLS =

selective laser sintering

T =

temperature

ε =

emissivity

σ =

standard deviation

σ =

Stefan–Boltzmann constant

ν =

Poisson's ratio

ρ =

density

References

1.
Nelson
,
J. C.
,
Xue
,
S.
,
Barlow
,
J. W.
,
Beaman
,
J. J.
,
Marcus
,
H. L.
, and
Bourell
,
D. L.
,
1993
, “
Model of the Selective Laser Sintering of Bisphenol-A Polycarbonate
,”
Ind. Eng. Chem. Res.
,
32
(
10
), pp.
2305
2317
.
2.
van Antwerpen
,
W.
,
du Toit
,
C. G.
, and
Rousseau
,
P. G.
,
2010
, “
A Review of Correlations to Model the Packing Structure and Effective Thermal Conductivity in Packed Beds of Mono-Sized Spherical Particles
,”
Nucl. Eng. Des.
,
240
(
7
), pp.
1803
1818
.
3.
Masamune
,
S.
, and
Smith
,
J. M.
,
1963
, “
Thermal Conductivity of Beds of Spherical Particles
,”
Ind. Eng. Chem. Fundam.
,
2
(
2
), pp.
136
143
.
4.
Cheng
,
S. C.
, and
Vachon
,
R. I.
,
1969
, “
Thermal Conductivity of Packed Beds and Powder Beds
,”
Int. J. Heat Mass Transfer
,
12
(
9
), pp.
1201
1206
.
5.
Gusarov
,
A. V.
,
Laoui
,
T.
,
Froyen
,
L.
, and
Titov
,
V. I.
,
2003
, “
Contact Thermal Conductivity of a Powder Bed in Selective Laser Sintering
,”
Int. J. Heat Mass Transfer
,
46
(
6
), pp.
1103
1109
.
6.
Slavin
,
A. J.
,
Arcas
,
V.
,
Greenhalgh
,
C. A.
,
Irvine
,
E. R.
, and
Marshall
,
D. B.
,
2002
, “
Theoretical Model for the Thermal Conductivity of a Packed Bed of Solid Spheroids in the Presence of a Static Gas, With No Adjustable Parameters Except at Low Pressure and Temperature
,”
Int. J. Heat Mass Transfer
,
45
(
20
), pp.
4151
4161
.
7.
Slavin
,
A. J.
,
Londry
,
F. A.
, and
Harrison
,
J.
,
2000
, “
A New Model for the Effective Thermal Conductivity of Packed Beds of Solid Spheroids: Alumina in Helium Between 100 and 500 C
,”
Int. J. Heat Mass Transfer
,
43
(
12
), pp.
2059
2073
.
8.
Xue
,
S.
, and
Barlow
,
J. W.
,
1990
, “
Thermal Properties of Powders
,”
Solid Freeform Fabrication Proceedings
, pp.
179
185
.
9.
Xue
,
S.
, and
Barlow
,
J. W.
,
1991
, “
Models for the Prediction of the Thermal Conductivities of Powders
,”
Solid Freeform Fabrication Proceedings
, pp.
62
69
.
10.
Yuan
,
M.
,
Diller
,
T. T.
,
Bourell
,
D.
, and
Beaman
,
J.
,
2013
, “
Thermal Conductivity of Polyamide 12 Powder for Use in Laser Sintering
,”
Rapid Prototyping J.
,
19
(
6
), pp.
437
445
.
11.
Sih
,
S. S.
,
1996
, “
The Thermal and Optical Properties of Powders in Selective Laser Sintering
,” Ph.D. thesis, Ann Arbor, MI.
12.
Vargas
,
W. L.
, and
McCarthy
,
J. J.
,
2002
, “
Conductivity of Granular Media With Stagnant Interstitial Fluids Via Thermal Particle Dynamics Simulation
,”
Int. J. Heat Mass Transfer
,
45
(
24
), pp.
4847
4856
.
13.
Zhang
,
H.
,
Zhou
,
Q.
,
Xing
,
H.
, and
Muhlhaus
,
H.
,
2011
, “
A DEM Study on the Effective Thermal Conductivity of Granular Assemblies
,”
Powder Technol.
,
205
(
1–3
), pp.
172
183
.
14.
Tsory
,
T.
,
Ben-Jacob
,
N.
,
Brosh
,
T.
, and
Levy
,
A.
,
2013
, “
Thermal DEMCFD Modeling and Simulation of Heat Transfer Through Packed Bed
,”
Powder Technol.
,
244
, pp.
52
60
.
15.
Widenfeld
,
G.
,
Weiss
,
Y.
, and
Kalman
,
H.
,
2003
, “
The Effect of Compression and Preconsolidation on the Effective Thermal Conductivity of Particulate Beds
,”
Powder Technol.
,
133
(
1–3
), pp.
15
22
.
16.
Feng
,
Y.
,
Han
,
K.
,
Li
,
C.
, and
Owen
,
D.
,
2008
, “
Discrete Thermal Element Modelling of Heat Conduction in Particle Systems: Basic Formulations
,”
J. Comput. Phys.
,
227
(
10
), pp.
5072
5089
.
17.
Garg
,
R.
,
Galvin
,
J.
,
Li
,
T.
, and
Pannala
,
S.
,
2012
, “
Documentation of Open-Source MFIX-DEM Software for Gas-Solids Flows
,” https://mfix.netl.doe.gov/documentation/qmomk_doc_2012-1.pdf
18.
Musser
,
J. M. H.
,
2011
, “
Modeling of Heat Transfer and Reactive Chemistry for Particles in Gas-Solid Flow Utilizing Continuum-Discrete Methodology (CDM)
,” Ph.D. thesis, West Virginia University, Morgantown, WV.
19.
Powell
,
M. J. D.
,
1964
, “
An Efficient Method for Finding the Minimum of a Function of Several Variables Without Calculating Derivatives
,”
Comput. J.
,
7
(
2
), pp.
155
162
.
20.
Batchelor
,
G. K.
, and
O'Brien
,
R. W.
,
1977
, “
Thermal or Electrical Conduction Through a Granular Material
,”
Proc. R. Soc. London, Ser. A
,
355
(
1682
), pp.
313
333
.
21.
Rong
,
D.
, and
Horio
,
M.
,
1999
, “
DEM Simulation of Char Combustion in a Fluidized Bed
,”
Second International Conference on CFD in the Minerals and Process Industries
, Melbourne, Australia, Dec. 6–8, pp.
65
70
.
22.
Marsaglia
,
G.
,
1972
, “
Choosing a Point From the Surface of a Sphere
,”
Ann. Math. Stat.
,
43
(
2
), pp.
645
646
.
23.
Zhou
,
J.
,
Zhang
,
Y.
, and
Chen
,
J. K.
,
2009
, “
Numerical Simulation of Laser Irradiation to a Randomly Packed Bimodal Powder Bed
,”
Int. J. Heat Mass Transfer
,
52
(
1314
), pp.
3137
3146
.
24.
Xiu
,
D.
, and
Karniadakis
,
G.
,
2002
, “
The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations
,”
SIAM J. Sci. Comput.
,
24
(
2
), pp.
619
644
.
25.
Xiu
,
D.
, and
Karniadakis
,
G. E.
,
2003
, “
A New Stochastic Approach to Transient Heat Conduction Modeling With Uncertainty
,”
Int. J. Heat Mass Transfer
,
46
(
24
), pp.
4681
4693
.
26.
Xiu
,
D.
, and
Hesthaven
,
J.
,
2005
, “
High-Order Collocation Methods for Differential Equations With Random Inputs
,”
SIAM J. Sci. Comput.
,
27
(
3
), pp.
1118
1139
.
27.
Ganapathysubramanian
,
B.
, and
Zabaras
,
N.
,
2007
, “
Sparse Grid Collocation Schemes for Stochastic Natural Convection Problems
,”
J. Comput. Phys.
,
225
(
1
), pp.
652
685
.
28.
Smolyak
,
S. A.
,
1963
, “
Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions
,”
Sov. Math., Dokl.
,
4
, pp.
240
243
.
29.
Bala
,
K.
,
Pradhan
,
P. R.
,
Saxena
,
N. S.
, and
Saksena
,
M. P.
,
1989
, “
Effective Thermal Conductivity of Copper Powders
,”
J. Phys. D: Appl. Phys.
,
22
(
8
), pp.
1068
1072
.
30.
Press
,
W. H.
,
Flannery
,
B. P.
,
Teukolsky
,
S. A.
, and
Vetterling
,
W. T.
,
1992
,
Numerical Recipes in C: The Art of Scientific Computing
, 2nd ed., Cambridge University Press, New York.