Comparison of AM microstructure simulation methods.
In this chapter, a predictive multiscale model based on a cellular automaton (CA)-finite element (FE) method has been developed to simulate thermal history and microstructure evolution during metal solidification for direct metal deposition (DMD) process. The macroscopic FE calculation that is validated by the thermocouple experiment is developed to simulate the transient temperature field and cooling rate of single layer and multiple layers. In order to integrate the different scales, a CA-FE coupled model is developed to combine with thermal history and simulate grain growth. In the mesoscopic CA model, heterogeneous nucleation sites, grain growth orientation and rate, epitaxial growth, remelting of preexisting grains, metal addition, grain competitive growth and columnar to equiaxed phenomena are simulated. The CA model is able to show the entrapment of neighboring cells and the relationship between undercooling and the grain growth rate. The model predicts the grain size and morphological evolution during the solidification phase of the deposition process. The developed “decentered polygon” growth algorithm is appropriate for the nonuniform temperature field. Finally, the single- and multiple-layer DMD experiments are conducted to validate the characteristics of grain features in the simulation.
- finite element
- cellular automata
- grain morphology
- direct metal deposition
- thermal modeling
Compared with the conventional subtractive manufacturing technologies, additive manufacturing (AM) has unique advantages including low heat input, small heat-affected zone, solid-free-form fabrication, near-net-shape and so on. Direct metal deposition (DMD), a rapid developing AM technique, is able to manufacture a fully dense metal part without intermediate steps, which is especially appropriate for the heterogeneous components manufacturing. During the deposition process, solidification thermodynamics determined by a series of process parameters affect microstructure evolution, which directly affects materials’ mechanical properties. The temperature field history and the cooling rate are the key factors to control the solidification microstructure after DMD process . Several approaches, including stochastic and deterministic, have been taken to model solidification microstructure evolution. Anderson and Srolovitz et al. [2, 3] developed a Monte Carlo (MC) stochastic method to simulate the grain growth, grain size distribution, curvature and growth rate as well as their interrelationships. Saito and Enomoto  incorporated the anisotropy of the grain boundary energy, the pinning effect of precipitates on growth kinetics into the MC simulation. Another idea of modeling is the deterministic approach. Chen  investigated a phase field (PF) method to model and to predict the mesoscale morphological and microstructure evolution in materials. C.E. Krill III, Böttger B, and Moelans N et al. [6, 7, 8] developed PF to simulate 2D grain growth, 3D grain growth and equiaxed solidification. However, a phase field model usually carries a very high computational cost because of a requirement for a particularly fine computational grid.
In order to reduce the computational cost, Rappaz and Gandin  put forward a 2D cellular automaton approach to model the grain structure formation in the solidification process. The model includes the mechanisms of heterogeneous nucleation and grain growth during the casting process. Nucleation occurring at the solid/liquid interface and the liquid bulk is treated by using different nucleation sites preference. The crystallographic orientation and locations of the grains are randomly selected among a certain number of orientation classes and millions of CA cells, respectively. However, the model was only applied to uniform temperature field. In order to develop the nonuniform temperature prediction, Gandin and Rappaz  proposed a 2D cellular automaton (CA) technique for the simulation of grain formation during solidification. The nonuniform temperature situation was fully coupled to finite element (FE) heat flow calculation with enthalpy. This progress made it possible to combine the temperature field history with the microstructure evolution. The coupled CA-FE model is applied to A1-7 wt% Si alloy. A 3D CA-FE model was analyzed for the prediction of dendritic grain structures formed during solidification . The potentiality of the CA-FE model is demonstrated through the predictions of typical grain structures formed during the investment casting and continuous casting processes. Based on the features of several developing approaches, Choudhury et al.  compared a CA model with a PF model for dendritic solidification of an Al-4 wt%Cu alloy, 2D and 3D at different undercooling conditions. In 2D case, there is a very good agreement of the simulated tip properties. At high undercooling, the CA model becomes more favorable, as its reproduction of the theoretical behavior is improved. Since the CA model can simulate at coarse scales during a relatively short time, its output can be employed as the input for a PF simulation in order to resolve finer details of microstructure formation within grains. This can be utilized to build a hybrid model to integrate CA high efficiency and PF accuracy. Dore  investigated quantitative prediction of microsegregation during solidification of the ternary alloy system, which is applied to solidification of Al–Mg–Si. Jarvis et al.  firstly compared 1D, 2D and 3D cellular automaton finite difference (CA-FD) simulations of nonequilibrium solidification in Al–3.95Cu–0.8Mg ternary alloy. It has been demonstrated that there is a good agreement between all the CA-FD models in terms of primary -Al phase. However, final dendrite arm spacings of 2D and 3D are slightly overestimated.
High cooling rate and nonequilibrium are typical characteristics of the DMD technique comparing conventional casting process and simulation. Grujicic et al.  proposed a modified CA-based method to investigate the evolution of solidification grain microstructure during the LENS rapid fabrication process. This research established the relationship between the process parameters (e.g., laser power, laser velocity) and solidification microstructure in binary metallic alloy. The finite difference analysis was also coupled with the modified CA to calculate the temperature field as the input of microstructure prediction. Kelly and Kampe [16, 17] developed the thermal history in DMD of Ti6Al4V and microstructural characterization. Nie et al.  developed a multiscale model to simulate microstructure evolution during laser additive manufacturing solidification. The study presented the relationship between the solidification conditions and the resultant microstructure, especially Laves phase particles in Ni-based superalloy. Rodgers et al.  proposed a 2D mesoscale model to simulate grain structure near a moving heat source with kinetic Monte Carlo simulator during electric beam melting (EBM) process. The method is capable of simulating both singlepass and multipass welds grain morphology. It also investigates the influence of initial substrate grain size on HAZ and FZ grain shape and size. Rai et al.  coupled a lattice Boltzmann (LB) and cellular automaton (CA) to simulate the microstructure evolution during electron beam melting. Initial grain selection at the base plate, grain boundary perturbation, grain nucleation due to unmolten powder particles in the bulk, grain penetration can all be simulated. The influence of process parameters on the final grain structure and texture evolution is analyzed. Keller et al.  investigated aspects of microstructure and microsegregation during rapid solidification in a laser powder bed fusion additive manufacturing process. Finite element analysis is employed to simulate the laser melt pool and temperature field. Microsegregation between dendritic arms is calculated by using the Scheil-Gulliver solidification model and DICTRA software. Phase field is developed to produce microstructures with primary cellular/dendritic arm spacing. However, there are few investigations on microstructure evolution prediction based on substrate and fusion zone during DMD process. Compared to other powder bed additive manufacturing process, there is different thermal cycle and the cooling rates for DMD process, which results in different microstructure. This part-level simulation on microstructure is critical because it provides the foundation for the prediction and control of mechanical properties.
CA simulation is appropriate for mesoscale modeling of grain structure because it does not consider much details inside a specific grain such as secondary dendritic arm spacing (SDAS) and microsegregation. Since it belongs to mesoscale model and does not cost as much computational resources as other microscale models, such as phase field and molecular dynamic, these characteristics of CA model make it appropriate for simulating the part-level grain structure instead of a very small region including less than a hundred grains. Thus, it can be used to predict and control the mechanical properties of the whole part based on the part-level grain structure under different parameters. Molecular dynamic (MD) simulation of microstructural evolution during additive manufacturing  is focused on a microscalability, which is between nanometer and micrometer. MD simulation can provide a method to investigate the crystallization process within the HAZ and clarify its crystallization mechanism because it is difficult to observe directly the crystallization process in the HAZ during the cyclic heating and cooling process. Even though MD can investigate microstructure evolution on a molecular level, it will cost too much resource to simulate the whole part structure and properties. PF model [21, 23], a microscopic one, can be used to simulate the solute concentration and phase transformation by solving the potential equation. The coefficients in the evolution equation of phase-field variables are related to the material parameters so that it can quantitatively simulate grain growth within a finer scale compared to CA and MC. Lattice Boltzmann (LB) method is adopted to numerically simulate the solute transport within the melt pool domain because it is appropriate for the complex geometry shape and is built from the temporal and spatial discretized grid, avoiding solving macroscopic N-S equations. The computational domain is discretized into regular lattices with the same cell size as the CA model. The governing equation and boundary conditions for transport process are described in detail in Refs. [24, 25, 26, 27]. The comparison of AM microstructure simulation methods is shown in Table 1. The current CAFÉ model can be used to consider multiple components if it considers the solute concentration and there is no chemical reaction and intermetallic phase formation during the solidification process. However, it is not capable of determining the mechanical properties directly if it is not incorporated with other models such as Hall-Petch model.
|Empirical microstructure models|
In this study, a predictive multiscale model based on a cellular automaton (CA)-finite element (FE) method has been developed to simulate thermal history and microstructure evolution during metal solidification for a laser-based additive manufacturing process shown in Figure 1. ABAQUS was used to calculate the temperature field of the whole part, which offers the macroscopic FE nodes’ temperature. In order to integrate the different scales, a coupled model is developed to combine with thermal history and simulate nucleation site, grain growth orientation and rate, epitaxial growth of new grains, remelting of preexisting grains, metal addition and grain competitive growth. Interpolation was utilized to obtain the finer nodes’ temperature based on the FE nodes result. The temperature field was validated by the type K thermocouples. The CA model, which was able to show the entrapment of neighboring cells and the relationship between undercooling and the grain growth rate, was built to simulate the microstructure information such as the grain size and columnar grain orientation. The developed “decentered polygon” algorithm is more appropriate for grain structure development in the highly nonuniform temperature field. This simulation will lead to new knowledge that simulates the grain structure development of single- and multiple-layer deposition during DMD process. The microstructure simulation results were validated by the experiment. The model parameters for the simulations were based on Ti-6Al-4V material (Figure 1).
2. Mathematical model
2.1. Ti6Al4V transient temperature field during the deposition process
In the direct metal deposition (DMD) process, the temperature history of the whole domain directly influences the deposition microstructure, which is critical to mechanical properties . In order to obtain the microstructure information during the solidification process, the temperature field must be known at each time step. The transient temperature field throughout the domain was obtained by solving the 3D heat conduction Eq. (1), in the substrate, along with the appropriate initial and boundary conditions .
where T is the temperature, is the density, is the specific heat, is the heat conductivity and Q is the internal heat generation following certain energy distribution per unit volume.
The initial conditions applied to solve Eq. (1) were:
where is the ambient temperature. In this study, was set as room temperature, 298 K. The boundary conditions, including thermal convection and radiation, are described by Newton’s law of cooling and the Stefan-Boltzmann law, respectively. The laser heating source term, in Eq. (1), was also considered in the boundary conditions as a surface heat source. The boundary conditions then could be expressed as 
where k, T, and Q bear their previous definitions, n is the normal vector of the surface, is the heat convection coefficient, ε(Τ) is the emissivity, σ is the Stefan-Boltzmann constant, which is W/, Γ represents the surfaces of the work piece and Λ denotes the surface area irradiated by the Gaussian laser beam.
In order to simulate the thermal history during the direct metal deposition more efficiently and reduce the computational cost, some assumptions were taken into account. In the experiment, a Gaussian distributed laser beam was utilized to melt the substrate vertically with a nonuniform power density . Thus, the transverse intensity variation is described as Eq. (4):
where is the laser absorption coefficient, P is the power of the continuous laser and is the distance from the beam axis where the optical intensity drops to (≈ 13.5%) of the value on the beam axis. was set as 0.4 based on numerical experiments in the LAMP laboratory and is 1 mm in this simulation. The motion of laser beam was simulated by adjusting the position of beam center R with programming a user subroutine “DFLUX” in ABAQUS. The formula of R is as follows:
where and are the spatial coordinates of the Gaussian laser beam center, and and are the laser moving velocities.
The Marangoni effect caused by the thermocapillary phenomena can directly influence the temperature field in the whole domain, so it is considered to obtain more accurate thermal history during DMD . The artificial thermal conductivity was put forward to address the Marangoni effect in the finite element method 
where is the modified thermal conductivity and is the liquidus temperature.
In the FEA model, the powder addition was simulated by activating elements in many small steps . The width of the deposit area is assumed to be the same as the Gaussian laser beam. The thickness of each layer is calculated by transverse speed, powder feed rate and powder absorption efficiency. The deposit geometry, boundary condition and heat flux were updated after each step.
2.2. Ti6Al4V morphology prediction after solidification
Heterogeneous nucleation occurs nearly instantaneously at a characteristic undercooling. The locations and crystallographic orientation of the new nuclei are randomly chosen at the surface or in the liquid. As explained by Oldfield , the continuous nucleation distribution, , which characterizes the relationship between undercooling and the grain density, is described by a Gaussian distribution both at the mold wall and in the bulk liquid. The parameters of these two distributions, including maximum nucleation density , the mean undercooling and the standard deviation of the grain density distribution , can be obtained from experiments and grain size measurements. The grain density, , is given by Eq. (7):
where is the maximum nucleation density of nucleation grains, which is obtained by the integral of the nucleation distribution (from zero undercooling to infinite undercooling). and are the mean undercooling and standard deviation of the grain density distribution, respectively. Here, all temperatures are in Kelvin.
Undercooling is the most important factor in the columnar and dendrite growth rate and grain size. The total undercooling of the dendritic tip consists of three parts such as solute undercooling, thermal undercooling and curvature undercooling. For most metallic alloys, the kinetic undercooling for atom attachment is small, so it is neglected . The total undercooling can be calculated as follows:
where m is the liquidus slope, is the Gibbs-Thomson coefficient, is the solute concentration in the liquid far from the solid-liquid interface, and are the thermal and solutal Pclet numbers, respectively, k is the solute partition coefficient at the solid-liquid interface, equals , is the unit thermal undercooling (and R is the radius of the dendritic tip.
For the laser deposition process, the rapid solidification condition corresponds to a high Peclet number at which the dendritic tip radius is given by Eq. (9)
where is the marginal stability constant, approximately equals , and and are the effective temperature gradient and concentration gradient, respectively.
2.3. Coupling macroscopic FE and mesoscopic CA models
The temperature field result can be used to calculate enthalpy increment, which is necessary to calculate enthalpy at each time step. A linearized implicit FE enthalpy formulation of the heat flow equation can be given 
where is the mass matrix, is the conductivity matrix, is the boundary condition vector and and are the temperature and enthalpy vectors at each node of the FE mesh, respectively. The Newton method and Euler implicit iteration are included in (10). This set of equations can be solved using the Gauss elimination method for .
Thus, the next time-step enthalpy can be obtained by the relationship of . The new temperature field can be obtained from the coupling model using (11). is the latent heat of fusion per unit volume. represents the fraction of solid. can be calculated as in .
In the FE macroscopic model, the temperature field was calculated on a relatively coarse mesh, but the solidification microstructure had to be developed on a finer regular CA mesh with a cell size in the order of the secondary dendrite arm spacing (SDAS). Figure 2 indicates the interpolate relationship between coarse FE nodes and fine CA cells. The known temperature and the volumetric enthalpy variation were interpolated into the CA network by the linear interpolation in Eqs. (12) and (13). is the interpolation coefficient. Every CA cell temperature in the calculation domain can be obtained with this interpolation.
The finer temperature, , and enthalpy variations, , in regular CA cells were used in Eq. (13) to yield the temperature in the next microtime step. After a few microtime steps, the temperature field in the CA network could be substituted into the coarser nodes of the macroscopic model. The interpolated temperature field is employed as the model input. Heterogeneous nucleation, grain growth orientation and grain growth are solved in the CA-FE model in terms of nucleation location distribution, random crystallographic orientation and CA cells capture. Figure 3 indicates the flow chart of coupling FE-CA model. The details of the CA growth algorithm are shown in Figure 4.
Figure 5 illustrates the conventional and modified cell capture algorithm. For the conventional method, the vertices of the square envelope move along the diagonal, and the growth of the square envelope is determined by the center cell temperature, not local temperature, at each time step, which results in the same growth rate for the four vertices. The grain orientation will be along with the axis of computational domain after a few time steps, thus, losing its original orientation information. The modified “decentered polygon” algorithm is implemented to control the grain growth within the melt pool and at the sold/liquid interface. Compared to the traditional “decentered square” algorithm of cell capturing, the modified “decentered polygon” algorithm does not need to create square for each cell when it begins to grow. Only the decentered polygon of a starting nucleated cell is tracked during the grain growth process, which reduces the computational cost. Besides, the modified algorithm can prevent grain orientations from realigning with x axis after a few growing steps because each cell will stop growing when Von Neumann and Moore neighbors are both solid. The controlling point growth rate is determined by the local cell temperature. Therefore, the region with higher thermal gradient will solidify faster along the steepest thermal gradient.
3. Results and discussion
3.1. Single-layer temperature and grain structure
The deposition temperature field and grain morphology were simulated first only in one layer. Figure 6 shows thermal history of the whole block during the DMD process. Figure 6(a) indicates the temperature field of the whole block when laser beam is passing along the x direction at time = 1.0 s, while Figure 6(b) shows the temperature field when substrate cools down with laser off at time = 29.0 s. The total physical time of single-layer laser deposition is 2 s, while the cooling time is 28 s in the simulation. For each step, the step time is 0.1 s when the laser is shot on the surface of the deposited material. After 30 s of cooling down, the temperature distribution is more uniform. Figure 7 indicates the thermal history of two nodes, which locate at the center point in the deposit and 1 mm away from the deposit. The result shows that the highest temperature in the deposit is approximately 2884 K, which occurred at the center of the Gaussian beam. The center node at 1 mm away from the deposit arrives at peak temperature of 1126 K that cannot melt the Ti6Al4V substrate. Based on every node’s thermal history, the undercooling (discrepancy between liquidus temperature and current temperature) that is critical to resulting in grain nucleation and growth rate can be determined.
In order that the input of microstructure model is reliable, the temperature field is validated with four type-K thermocouples. The locations are shown in Figure 8. One is located at the starting of laser path, which distance to the laser is approximately 3–3.5 mm. Another three points are located by one side of laser path, which distance to the laser center is approximately 2 mm. Arduino device is used to sample the temperature data. A laser deposition experiment is conducted with the power of 750 W, scanning speed of 600 mm/min and 2 g/min for single-layer deposit. The difference between the experiment and the FEM modeling is less than 10°C shown in Figure 9. In Figure 9(a), the delay between the simulated temperature and the thermocouple itself is more visible than Figure 9(b)–(d) because the distance of first thermocouple point is further than other three ones. In the real experiment, the substrate is fixed by the metal fixture, which resulting in the more heat conduction than the FEM model. Because of argon gas, forced convection occurred in the real experiment. This also causes lower cooling rate in the temperature simulation. Because the difference between experiment and simulation is smaller than 10%, the current FEA modeling is still considered as a reasonable simulation of temperature field, which can provide the reliable thermal input for the CA model.
A laser deposition experiment is conducted with the power of 700 W, scanning speed of 600 mm/min and 2 g/min for single-layer deposit. For this case, the cross section shown in the figure is the computational domain. The cell size for this simulation is 6 × 6 μm. X and Y axes represent the number of cells. The simulation result from conventional method is shown in Figure 10. It can be observed that even though different grains own diverse orientation at the very beginning, the crystallographic orientation preference tends to be along with the axis after several time steps. Here, different colors represent various grain orientations. Finally, the equiaxed grains dominate the fusion zone. The original grain orientations are not kept during the solidification process. It does not agree well with the single-layer experimental result shown in Figure 12.
The developed CA grain growth method is implemented under the same condition. According to the developed CAFÉ simulation, the single layer simulation result is shown in Figure 11. The grain keeps its original crystallization orientation when grain growth is modeled. The columnar grain can be identified from the solid/liquid interface. When grains continue to grow toward melt pool center, some grains overgrow other grains such that there are fewer grains further away from the solid/liquid interface.
Three samples of single-layer deposits are prepared with EDM cutting, grinding, polishing and etching. The optical microscope is shown in Figure 12. The comparison between simulation and experimental results is shown in Figure 13. An average of 20 measurements per sample is performed to determine the average grain size. It compares the experimental average grain size with the predicted one. The shown data suggest that a 15% error between experimental measurements and predictions is present. This can be considered as a reasonable prediction of grain morphology and size.
3.2. Multilayer temperature and grain structure
Figure 14 depicts the temperature field of the substrate and deposited material, including the 25-layer deposition materials added on the substrate when the laser is moved forward and backward. The laser deposition of multiple-layer Ti-6Al-4V was conducted with the power of 750 W, scanning speed of 200 mm/min and powder delivery of 2 g/min. The elemental size is nonuniform along the three directions because it is not necessary to apply fine elements to where the location is far from the molten pool. Figure 15 shows that the thermal history and peak temperature of different layers are not identical. The higher layer performs higher thermal history because the higher layer accumulates more heat than the lower one, and it is closer to heat source.
Figure 16 shows Ti-6Al-4V deposition grain microstructure. The cross-sectional dimension of deposit region is 1.8 1.9 mm, which is close to 2 2 mm assumption in the simulation. In Figure 16, it can be observed that at the bottom deposition, crystallographic orientation is not only limited to the vertical direction. It can also be observed that columnar grains dominate in the laser deposition area. Figure 16(a) and (b) indicates the whole deposition region at different magnification and the locations of top and bottom region, while Figure 16(c) and (d) shows the grain size and shape with higher magnification. Under the same condition, the experiment is conducted, and the optical microscope images are taken. Figure 16(e) shows multiple layers of the Ti-6Al-4V grain morphology under the laser deposition process. Irregular grain shape and size can be obtained. When more layers were deposited, prior columnar grains began to dominate, while equiaxed grains began disappearing. As the solidification process continues, competitive growth among different grains occurs. Therefore, the size of columnar grain increases, and the number of grains goes down. The orientations of the columnar grains were almost perpendicular to the laser motion’s direction because the grains grew along the steepest thermal gradient direction. This phenomenon verifies the columnar grain orientation in the simulation result. The domain size in the CA model was After measurement of grain size, it can be found in Figure 17 that in the simulation, the grain size ranges from 113 to 346 μm. For the experiment, the grain size ranges from 156 to 599 μm. The grain size at the bottom and top is larger than the simulation. This may be because it does not consider the cyclic heating and cooling process’ effect on the solidified grain evolution. Usually, cyclic heating will coarsen the grain and make the grain become larger. This effect will be solved in the future research task.
Figure 18 presents the simulated grain structure from Rai et al.  during powder bed additive manufacturing. It can be seen that some grains overgrow others at the top layers, and most surviving grains have negative misorientations indicating grain orientation is aligned well with the beam scanning direction. The detailed local grain boundary misorientation is determined by local thermal gradient and the neighboring grains’ orientation. The rate of overgrowth process also has an effect on the grain boundary angle. Compared to multiple layer results in this investigation, it shows the similarity of grain overgrowth mechanism and misorientation distribution. The grain size between two results is not similar because the thermal gradient and cooling rate are different between powder bed-based additive manufacturing and DMD process.
The transient temperature field of single-layer and multiple-layer deposition of Ti-6Al-4V was simulated with finite element method. The simulation result was validated by thermocouple experiment. The FE model provides the temperature at a relatively coarse scale (200 ), and interpolation algorithm was used to scale the temperature field to match that of the CA model. The FE-CA model predicts grain morphology evolution as the deposition cools down. Hence, the instantaneous nucleation law, grain growth and crystallographic orientation were modeled in this study. It has been found that the developed “decentered polygon” growth method is more appropriate for the highly nonuniform temperature field, and the simulation result is closer to the real experimental measurement compared to the conventional growth method. For multilayer deposit, columnar grains dominated in the 25-layer deposition in the simulation. The grain size becomes larger when the position is closer to the top area of the deposition, which matches well with the optical microscopic result. The grain size of single and multiple layers between simulation and experiment is similar. It demonstrates that this FE-CA simulation can reasonably predict thermal history and grain morphology during this case of direct metal deposition.
This work was funded through NASA’s Fundamental Aeronautics Program, Fixed Wing Project, under NRA NNX11AI73A.