I. INTRODUCTION
Direct metal laser sintering (DMLS) is an additive manufacturing technique that uses a laser as the power source to sinter powdered material. In this process, laser automatically aims at points in space defined by a three-dimensional (3D) model, binding the material together through sintering process to create a solid structure. 1 An important type of metal used in DMLS technique is nickel based alloy, e.g., Inconel 625 and 718. 2 Inconel alloys are oxidation and corrosion resistant materials which are well suited for service in extreme environments subject to pressure and heat. 2 Therefore, it is important to understand the sintering behavior and the resultant mechanical properties of nickel based metals.
There are several modeling studies about sintering of nickel particles. Xu et al. developed a molecular dynamics simulation method to investigate the sintering of nickel nanoparticles in nickel and yttria-stabilized zirconia (Ni/YSZ) anode of solid oxide fuel cell. Reference Xu, Sakanoi, Higuchi, Ozawa, Sato, Hashida and Kubo3 The model reveals that the YSZ nanoparticle framework suppresses the sintering of nickel nanoparticles by disrupting the growth of the neck between two nickel nanoparticles. Reference Xu, Sakanoi, Higuchi, Ozawa, Sato, Hashida and Kubo3 Yang studied the neck growth in the sintering of different sized gold (100) nanoparticles under different heating rates. Reference Yang4 The results show that a stable neck width increases as the size of the other nanoparticle increases. Reference Yang4 Cheng et al. investigated the structural change, grain growth, and mass transfer mechanisms of copper nanoparticles sintering by molecular dynamics simulation. Reference Cheng and Ngan5,Reference Cheng and Ngan6 Henz et al. conducted molecular dynamics simulations to gain understanding of alloying reaction of Ni and Al nanoparticles during self-propagating high-temperature synthesis. Reference Henz, Hawa and Zachariah7 The effect of primary particle size on sintering time and sintering temperature for separate nanoparticles was investigated. Zhang et al. conducted sintering simulation of iron particles in DMLS process. Reference Zhang, Wu, El-Mounayri, Brand and Zhang8 The diffusion coefficient and mechanical strength of sintered iron parts were studied.
Most of the previous efforts were focused on the configuration changes during sintering process. The mechanical aspects of the sintered particles are rarely reported. This work is aimed to provide insights to answer a fundamental question in direct DMLS, “What is the effect of heating rate due to sintering laser beam power on the mechanical strength of a 3D printed part?” Although it is impossible to answer the above question thoroughly in a single study, the models developed in this work may provide some guidance of laser sintering processing parameters to design printed parts with desired mechanical properties. Three cases of simulations will be conducted to study the sintering phenomena and the resultant mechanical properties, including an assembly of 41 Ni nanoparticles, bulk Ni, and a Ni–Ni nanoparticle pair. The bulk Ni and the Ni–Ni nanoparticle pair are used for comparison propose. For the Ni–Ni nanoparticle pair, the activation energy for the sintering process is also investigated.
II. SIMULATION DETAILS
Molecular dynamics software package large-scale atomic/molecular massively parallel simulator (LAMMPS) Reference Plimpton9 was used for all simulations. For Ni particles, the embedded-atom method (EAM) Reference Bonny, Pasianot and Malerba10,Reference Daw and Baskes11 was used for describing the calculation of atomic interaction. In EAM, the total energy E i of an atom i is given by
$${E_i} = {F_\alpha }\left( {\mathop \sum \limits_{j \ne i} {\rho _\beta }\left( {{r_{ij}}} \right)} \right) + {1 \over 2}\mathop \sum \limits_{j \ne i} {\emptyset _{\alpha \beta }}\left( {{r_{ij}}} \right)\quad ,$$
where F is the embedding energy that represents the energy required to place atom i of type α into the electron cloud. ρβ is the contribution to the electron charge density from atom j of type β at the location of atom i, ϕαβ is a pair-wise potential interaction, and α and β are the element types (both were Ni in this work) of atoms i and j, r ij is the distance between atoms i and j. The multi-body nature of the EAM potential is a result of the embedding energy term. Both summations in the formula are over all neighbors j of atom i within the cutoff distance (4.6 Ǻ was used in this work).
Three groups of atomistic models were developed: (i) case 1 is a Ni–Ni nanoparticle pair; (ii) case 2 is an assembly of 41 Ni nanoparticles; and (iii) case 3 is a Ni face centered cubic (FCC) crystal bulk. Sintering and tensile test simulations were conducted on case 1 and case 2, and tensile test on case 3. For the 41 Ni nanoparticles, each particle has a diameter of 2.8 nm. The simulation block has 3 layers of particles with total 52,670 atoms. Using this model, the mechanical response of the sintered material was studied. For comparison, a bulk nickel FCC crystal and a two-particle pair model were also studied. Using the Ni–Ni nanoparticle pair model, the diffusion process and activation energy of this two-particle model during solid sintering process can be calculated.
A. Case 1: diffusion and sintering study of nickel–nickel nanoparticle pair
Sintering kinetics has been studied by many investigators from both experimental and theoretical points of view. Coble divided the sintering process into three successive elementary stages. Reference Coble12 The initial stage corresponds to neck formation and growth, the intermediate stage to the growth of particle vacancies, and the final stage to diffusion and disappearance of the vacancies. In particular, the initial stage of sintering is when the inter-particle contact area increases from zero to ≅0.2 of the cross-sectional area of the particle. This initial stage is frequently referred as neck growth. The intermediate stage begins after grain growth and pore shape change causes transition to a pore and grain boundary matrix consisting of equilibrium dihedral angles formed on the solid–vapor (pore) surface at the locus of intersections with the solid–solid (grain boundary) surfaces with singly curved surfaces between intersections. The final stage begins when the pore phase is eventually pinched off. Reference Coble12
1. Atomistic model set up
For the two-particle pair model, a spherical nickel nanoparticle was cut out from a bulk FCC nickel crystal. Another nickel nanoparticle was duplicated and displaced to a contact distance with respect to the first one. The Ni–Ni pair model was built by adding these two particles together. Each particle has a diameter of 2.8 nm with 1265 atoms. An initial particle center to particle center distance was set to 2.8 nm. Figure 1 shows the nanoparticle pair atomistic model.
Initial configuration of the Ni–Ni nanoparticle pair model.
2. Solid state sintering of Ni–Ni nanoparticle pair
A 3D simulation box was created with periodic boundary conditions in each direction. All atoms in the atomistic model were randomized at 300 K, with an initial velocity following the Maxwell–Boltzmann distribution. To get the optimized lattice parameter and relaxed structure, the particles were equilibrated at room temperature (300 K) with a time step of 1 × 10−15 s for the equation of motion integration. During this process, the initial neck growth occurred due to the strong local potential gradients. Reference Pan, Ko and Grigoropoulos13 After the model equilibrium, an external heat was applied to the model. The temperature was changed from 300 K to 834 K. For this part of calculation, the system was constrained at a fixed heating rate. The temperature is linearly increased from 300 K to 834 K in 30 ps. The dimension changes after the simulations were extracted and analyzed.
3. Diffusion of Ni–Ni nanoparticle pair
To calculate the diffusion phenomena during sintering, the mean-squared displacement (MSD) needs to be calculated. 14 The atoms in a nanoparticle were divided into two groups: in (particle core) and out (particle shell). MSD allows defining a computation that calculates the MSD of the group of atoms, including all effects due to atoms passing through periodic boundaries. A vector of four quantities is calculated by this compute. The first 3 elements of the calculated vector are the squared d x , d y , d z displacements, summed and averaged over atoms in the group. The 4th element is the total squared displacement, i.e., (d x × d x + d y × d y + d z × d z ), summed and averaged over atoms in the group. The slope of the MSD versus time is proportional to the diffusion coefficient of the diffusing atoms. 14
The diffusion coefficient is assumed to follow the Boltzmann–Arrhenius dependency Reference German and Munir15,Reference Elliot and Munir16 :
$$D = {D_0}\exp \left( { - {{{E_{\rm{a}}}} \over {RT}}} \right)\quad ,$$
where D 0 is the pre-exponential diffusion coefficient and E a is the activation energy, R is the gas constant. A plot of the natural logarithm of D/D 0 at various temperatures versus 1/T should give a straight line whose slope is −E a/R and can be used to derive the activation energy E a.
B. Case 2 and case 3: tensile test simulations of sintered nickel nanoparticles and bulk nickel
1. Atomistic model set up
To investigate the mechanical properties of the structure after sintering, a larger model with 41 particles was created for MD calculations. Figure 2 shows the structure of this model from top and side views. 41 spherical particles with a diameter of 2.8 nm were cut out of a bulk FCC nickel block. Three layers of particles were constructed, with the top and bottom layers of 4 × 4 particles, and the middle layer of 3 × 3 particles. Additionally, the two-particle models used in Sec. II. A (Fig. 1) and a bulk nickel FCC crystal with the same simulation box dimensions were calculated, to compare the tensile response results with the larger 41-particle simulation model.
Atomistic model of tensile simulation of the 41-particle Ni sintered structure. (a) Top view; (b) side view.
2. Modeling of sintering and tensile test of sintered structure
Three processes were modeled in this simulation: (i) sintering simulations: which are NPT processes. NPT refers to constant number of atoms (N), constant pressure (P), and temperature or temperature changing rate (T). In NPT calculations implemented LAMMPS, temperature can be fixed or temperature changing rate can be fixed. In this work, a constant temperature increase rate was used through using LAMMPS command line “fix all npt temp Tinitial Tfinal 1 iso 0 0 1 drag 0.1”, where T initial is the initial temperature, and T final is the final temperature. In terms of temperature uniformity, the 41 Ni nanoparticle system has a cross section area of approximately 12 × 12 nm2, which is much smaller than a typical laser beam spot size of 70 μm. Reference Qian, Taimisto, Lehti, Piili, Nyrhilä, Salminen and Shen17 Therefore the heating can be reasonably assumed to be uniform. In this study, the temperature of the system increased from room temperature (300 K) to the sintering temperature of nickel particles at constant pressure. According to traditional laser sintering, the heating rate is in the range of 0.89 × 1011 to 2.20 × 1011 K/s with a laser intensity of 1010 W/m2. Reference Pan, Ko and Grigoropoulos13 In the sintering simulations, several different heating rates were applied to study the effect of heating rate on stress–strain response. To this end, different final temperatures with the same simulation time of 6000 ps were used. The time step for motion equation of integration was set to 0.01 ps. The total sintering simulation times were 600,000 steps.
(ii) Equilibrium at room temperature: since the tensile test simulation was conducted at room temperature (300 K), whereas the final temperature of the sintering simulation was finished at a higher temperature. It is necessary to stabilize the sintered structure at room temperature before the tensile calculations. This part of simulation also used NPT for the entire system, with a constant temperature of 300 K, without external pressures. The system reached a stable state after 30 ps.
(iii) Tensile test simulations: as the sintered structure reached a stable configuration at room temperature, an elongation rate of 10−9/s was applied to the model in the elongation direction. The Ni nanoparticles were pulled in the tensile direction to a strain till 0.2. The strain was increased for a specified number of times in a loop, and the stress was calculated at each point. In LAMMPS, the stress tensor (pressure) of the entire system of atoms is calculated by using command “compute pressure”. 14 A symmetric pressure tensor, stored as a 6-element vector, was also calculated by this compute command. 14 The 6 components of the vector are ordered xx, yy, zz, xy, xz, and yz. In this study, only the xx component is used to derive the stress. The stresses in yy and zz direction were minimal. The NPT condition was again applied to the whole system with a constant temperature of 300 K and 0 external pressures. During the calculations, the initial model length in elongation direction was recorded, and the lengths at each 100 steps were calculated to get the strain that was applied to the model at each time. In the meantime, the pressure of the system along the elongation direction was calculated to extract the stress value of the model using pressure unit conversion from bar to GPa. At the end, the stress–strain curves of the system were plotted, and the results from different heating rate simulations were compared.
III. RESULTS AND DISCUSSION
A. Temperature and configuration change of Ni–Ni nanoparticle pair sintering
In the simulations, the temperature was increased from 300 K to 834 K of Ni–Ni nanoparticle pair in 30,000 time steps with each step of 0.001 ps. The total simulation time was 30 ps. The temperature profile of this simulation is shown in Fig. 3. The initial temperature was set to 300 K. Small oscillations can be seen at the beginning (before 3 ps) and end of the simulation (after 20 ps). Overall, the temperature is increased linearly from 300 K to 834 K.
Temperature history of the Ni–Ni nanoparticle pair sintering simulation.
Figure 4 shows the configuration changes of the Ni–Ni nanoparticle pair sintering model during the 30 ps period. As shown in Fig. 4(a), the initial structure of these two nanoparticles was barely contacted each other. As shown in Fig. 4(b), when temperature was increased, the atoms near the interface of the two particles migrated together to minimize the surface energy of the system, leading to a neck growth. Because the movement of atoms occurred near the neck, the atoms on the surface of the particle lost their crystalline structures and formed a pre-melt layer. Reference Moitra, Kim, Kim, Park, German and Horstemeyer18 As shown in Fig. 4(c), when the temperature is kept increasing, the neck size became larger, and a thicker pre-melt layer is formed. Also as shown in Fig. 4(c), some Ni atoms near the boundary diffused to the other particle at 30 ps.
Configuration changes of Ni–Ni pair during sintering. (a) Initial configuration, (b) 15 ps, (c) 30 ps.
B. Diffusion of Ni–Ni nanoparticle pair
In Fig. 5, the diffusivities (D) at five temperatures, 700, 800, 900, 1000, and 1100 K were calculated using the MSD. The diffusivity of each temperature was obtained by fitting the stabilized region of the MSD (220–320 ps in our simulations) to avoid the uncertainty caused by initial particle interaction. As expected, the diffusivities increased with increasing temperature for both in (particle core) and out (particle shell) regions. The out region had a higher diffusivity than the in region due to combined surface and grain boundary diffusions, while the in region had less diffusion due to mainly volume or lattice diffusion mechanisms which were effective inside of the particle.
Diffusivities of In (particle core) and Out (particle shell) regions at the temperature range of 700 K–1100 K. The linear fitted lines and corresponding regression equations are also shown.
Several possible mass transfer mechanisms have been proposed for sintering: viscus flow, evaporation condensation, plastic dislocation deformation, and diffusion. Reference German19 For nano-size particles sintering, only diffusion mechanism is typically considered. Reference Pan, Ko and Grigoropoulos13 To understand the diffusion mechanisms in this work, we compared the activation energy of particle surface (surface diffusion) and particle core (volume diffusion) using Eq. (2). The calculated activation energies of Ni diffusion is 6.10 kJ/mol in the particle core, and 6.24 kJ/mol on the particle surface, which were lower than the experimentally reported value of 7.89 kJ/mol. Reference Elliot and Munir16 Although the surface activation energy was higher than that in the core, the surface diffusion was still higher than the core (Fig. 5), due to a larger Do, Do is the pre-exponential diffusion coefficient [Eq. (2)].
Surface diffusion and volume diffusion have similar activation energy, suggesting that mass transfer in our model can't be simply explained by only one single dominating diffusion mechanism. As suggested by Cheng et al., Reference Cheng and Ngan20,Reference Cheng and Ngan21 sintering of nano-sized gold particles is reasonably explained by a dislocation-free mechanism involving shear and rigid-body rotation processes caused by elastic waves. Movements of atoms are mostly driven by crystal structure transformation rather than any conventional diffusion mentioned above. Reference Cheng and Ngan20,Reference Cheng and Ngan21 In our work, the vector plot of atom displacements is shown in Fig. 6(a) to illuminate atom movements during sintering. By analyzing movement of atoms, neither surface diffusion nor grain boundary diffusion pattern can be identified. Atoms move in a relatively small range to change the internal crystal structure. This result is consistent with the analysis in Cheng's work. Reference Cheng and Ngan20
(a) Displacement vector plot during the sintering simulation, and the sintered structures with following sintering temperatures: (b) 1000 K, (c) 1200 K, and (d) 1300 K.
At the end of sintering, different dislocation-free twinning structures were observed [Figs. 6(b)–6(d)]. The finalized sintered structures at temperature range of 700–1200 K were quenched to 0 K for structural analysis. The color shown is defined by localized common neighbor analysis, where green color represents FCC and red represents HCP. The structures of sintering at lower temperatures (700–1100 K) resulted in almost single crystal FCC structure. For sintering at temperature of 1200 K, the sintered structure ended up with a 5-fold twinning structure. Sintering at 1300 K generated two parallel twin boundaries as shown in Fig. 6(d). All of the structures are in agreement with the simulated result by Cheng. Reference Cheng and Ngan20 Thus, it verifies that the mass transfer in our two-particle model is due to the dislocation-free elastic wave.
C. Mechanical response of sintered nanoparticles
For all simulations in the 41 particles model, the initial temperature and simulation time were fixed, namely 300 K and 6000 ps, respectively. Temperatures in the simulations were carefully controlled using the laser heating rates from 0.89 × 1011 K/s to 1.48 × 1011 K/s corresponding to a laser intensity of 1010 W/m2. Reference Pan, Ko and Grigoropoulos13 Different heating rates were achieved by adjusting the final temperature at the end of simulation. Table I shows the final temperatures according to the heating rates of the laser beam. The temperature evolutions of the five cases are plotted in Fig. 7. The temperature briefly oscillated at the beginning of the simulation and then converged quickly. Temperatures increased linearly during the 6000 ps period, and reached the designed final temperature.
Five heating rates applied to the 41 Ni nanoparticle model with their final heating temperatures.
Temperature evaluations of the 41 Ni nanoparticles sintered at different heating rates and final heating temperatures.
Figure 8 shows an example of the 41-Ni nanoparticles configurations at different stages of sintering at 834 K. The model has 3-layer nanoparticles. Each particle has the same size and a perfect crystal structure. After the sintering, a neck growth between each particle was observed, and atoms filled up the space near the initial contact regions.
Sintering of the 41 Ni nanoparticles (a) initial configuration; (b) after sintering at 834 K.
For the MD tensile test simulation, a tensile strain was applied in the horizontal direction of the model. The pressures of the simulation systems were extracted, and the stress data were plotted against strain from 0 to 0.2. Both the two-particle model and 41-particle model under tensile loading were simulated. Dislocation and common neighbor analysis were conducted to analyze the structural change during the tensile process as shown in Fig. 9 using displacement vectors. Color curves in the figure show dislocation lines based on localized structure. Reference Alexander and Karsten22 The two-particle model's plastic deformation began at strain of 0.10, as shown in Fig. 9(a) of initial occurrence of dislocations. In Fig. 9(b), the number of dislocation segments increased rapidly during strain of 0.10–0.11, with a strong atomic motions along dislocation lines, indicating large plastic deformation. With further application of strain, the dislocations were annihilated by combining dislocation lines as shown in Fig. 9(c).
Structural change of two-particle model under tensile load at various strains: (a) 0.10, (b) 0.11, and (c) 0.20. Colors show the dislocation type: blue is 1/2 〈110〉 perfect, green is 1/6 〈112〉 Shockley, purple is 1/6 〈110〉 stair-rod, and red is other.
Similar process of structural change of the 41-particle model was observed in Fig. 10. For this model, although there were some dislocations at zero strain after sintering, the number of dislocation segments did not increase until the plastic deformation. After this point, the number of dislocation was kept almost constant.
Configurational change during simulated tensile test of the 41-nanoparticle model sintered at 834 K. (a) Initial configuration after sintering; (b) strain 0.11; (c) strain 0.2. Colors represent crystal structures: red is HCP, and green is FCC.
To get a deeper understanding on the effect of the generation and evolution of dislocation during the tensile load, we plotted the dislocation length simultaneously with the tensile stress–strain curve in Fig. 11. From the stress–strain curve, one can easily observe that linear elastic deformation occurred at the strain less than ∼0.06. In this range, the total length of dislocation was small (less than 750 Å) and almost did not change. At the strain range of 0.06–0.10, the stress–strain reached to plastic yield strength, while the corresponding dislocation length sharply increased to 2700 Å, illustrating that the plastic deformation was caused by dislocation multiplication. After strain greater than 0.10, the stress decreased, while the dislocation length remained almost constant, suggesting the deformation in this stage was not controlled by dislocation, but primarily dominated by increasing damaging voids.
Total dislocation length and stress with respect to tensile strain of sintered 41 Ni particles at 834 K.
Figure 12 shows the stress–strain curves of the all cases with different heating rates, including the two-particle model, the 41-partile model, and FCC bulk Ni crystal for comparison. At lower strain range, all the models had linear increments of stress, until the yield stress was reached at the strain range of 0.05–0.1. Comparing the results, the sintered 41-particle structures had lower mechanical strengths than the bulk crystal, due to their porous structure and non-periodic nature of the nickel atoms. The Ni–Ni two-particle pair showed the lowest strength among all models, since the small neck size between two nickel particles. By comparing the sintered structures at different heating rates, the model with a higher heating rate showed a higher mechanical strength. The reason is that during the same sintering time, higher heating rates promote densification to reduce porosity, and resultant density is more close to the bulk structure. Oscillations of the stress–strain curves after yield were observed in Fig. 12. Due to relatively small volume size of the simulation system, the calculated stresses after yield are greatly affected by the increased porosity and the cut-off distance used in the EAM potentials.
Stress–strain curve of sintered nickel nanoparticles with various heating rates.
The bulk Ni simulation in Fig. 9 shows a calculated Young's modulus of 127 GPa, which is lower than the reported 200 GPa in polycrystalline FCC Ni. 23 The underestimation is likely due to the accuracy of EAM potential used in this work.
IV. CONCLUSION
The sintering process and the resultant mechanical strength of nickel nanoparticles are studied using molecular dynamics simulations. The conclusions are summarized below.
-
(1) Nickel nanoparticle pair diffusion coefficient as well as diffusion activation energy during sintering are studied. The diffusion of particle surface is higher than the particle core. The calculated activation energy of nickel particle diffusion is 6.10 kJ/mol in the particle core, and 6.24 kJ/mol on the particle surface, which are reasonably in agreement with experimental data of 7.89 kJ/mol. The sintered model shows a 5-fold twinning structure at 1200 K, and two parallel twin boundaries at 1300 K.
-
(2) Atomistic models of tensile tests of sintered material are developed. Uniaxial tensile responses of sintered nickel material with different heating rates are simulated and compared against the bulk nickel model. Deformation mechanism of sintered structure under tensile load is illustrated from the correlation between stress and dislocation evolution. The results show that sintered nickel nanoparticles have lower mechanical strengths than the bulk nickel crystal because of their porous structures. Higher heating rate leads to a higher mechanical strength because of accelerated sintering rates.
ACKNOWLEDGMENTS
J.Z. acknowledges Walmart Foundation's support of research project “Optimal Plastic Injection Molding Tooling Design and Production through Advanced Additive Manufacturing”. We also thank the critical comments by the anonymous reviewer.
