1 Introduction

Additive manufacturing (AM) of metal components allows the generation of very complex parts even from high performance alloys such as titanium or nickel base superalloys.[1,2,3,4,5] Generally, powder bed fusion (PBF) based on the laser (L-PBF)[6] or the electron beam (E-PBF)[7] or direct energy deposition (DED)[8] also known as laser engineered net shaping (LENS) are used for high materials quality. The DED process is characterized by injecting metal powder into a melt pool generated by a coaxial laser. In contrast, the beam locally melts powder in a powder bed during PBF.

Generally, AM is involved with rapid and directional solidification. Typically, the solidification velocity decreases from L-PBF to E-PBF to DED. The component evolves during several hours in a line-by-line and layer-by-layer process defined by a variety of process parameters such as beam velocity, beam power, distance between lines, layer thickness, line length, etc. The resulting microstructure is governed by all these parameters since they determine the local solidification conditions, i.e., the solidification velocity \( v \) and the thermal gradient \( G \) at the solidification front.[9] Thus, the microstructure reflects the scanning strategy with the inherent periodicity and symmetry. We observe elongated grains following the solidification front as well as equiaxed microstructures.[10,11,12] The directional solidification conditions very often lead to strong texture formation, typically the <001> orientation is preferred.[13] Specific scanning strategies may even result in single crystals.[14]

As a consequence, the mechanical properties of additively build parts are dependent on the build direction.[15,16] In most cases, this is not desired and a variety of strategies has been developed to break up this kind of texture formation.[17,18] On the other hand, the possibility to manipulate locally the microstructure and thus the mechanical properties is also a unique opportunity inherent to AM. No other production technology is able to adjust locally the materials properties. Apart from the ability of producing near net shape parts with complex geometries, this is an important feature of AM methods, which has not been exploited so far and requires full understanding of microstructure evolution during the building process. Although microstructure evolution of single tracks is properly understood, the line-by-line and layer-by-layer process results in very complex patterns due to multiple remelting. Suitable simulation tools can only capture this complexity.

There are different physical and numerical models available to simulate microstructure evolution during solidification[19,20]: phase-field models,[21,22,23,24] cellular automaton models[25,26] or Monte Carlo[27] models. These models differ by their physical exactness and the necessary computational effort. The more physics is contained the higher is the computational effort.

Phase-field (PF) models consider the complex solidification pattern formation and segregation of an alloy based on a sound thermodynamic basis.[23] Due to the extreme high computational effort, PF models are usually restricted to alloys with only some components, typically two or three elements. In addition, only a small number of cells or dendrites is simulated. The lack of reliable material properties, e.g., the interface energy as a function of direction, complicates matters even more. Nevertheless, PF models not only give information about the grain structure but also the dendritic segregation patterns. This method gives the highest level of information but is restricted to very small sections of a part.

Cellular automation (CA) approaches have a coarser look on the cellular or dendritic microstructure. In this case, only the convex hull of a dendrite is considered.[28] There is no information of the internal segregation structure of a grain. The growth kinetics of the dendrite tips is governed by the local undercooling. Here, PF calculations may give information about the undercooling situation at the dendritic tips as well as the growth rate—undercooling dependence.[29]

Kinetic Monte Carlo models are based on curvature driven growth of phase boundaries.[30] Although large problems, i.e., dimensions in the order of a component, can be solved the correlation between the numerical model and the physics behind grain growth and grain selection is not quite clear.

All these methods show specific advantages and disadvantages.[31,32] Nevertheless, they have one common challenge – modeling of new grain nucleation.[28] Nucleation is a key problem and far from being solved. Nucleation models available in literature exhibit several free parameters. These parameters have substantial influence on the numerical result and are used to fit experimental and numerical results.[33]

The aim of this work is to review the theoretical approaches to simulate AM microstructures. A short overview of the most important phenomena observed in experiments gives us a sound physical basis for the different mechanisms governing grain structure evolution during AM. Eventually, the numerical models have to be validated against experimental findings. Basis for modeling the microstructure is a well-defined temperature field. Different approaches to model the temperature at the solidification front are described and their specific advantages and disadvantages are discussed. The main part of the review concerns the various approaches modeling grain structure evolution including a critical discussion on their predictive power to calculate microstructures in real components.

2 Experimental AM Microstructures

Basically, AM of metals is a welding process where a high power beam melts a powder layer or powder is injected into the melt pool produced by a laser beam. Thus, most of the phenomena observed during welding[34,35] are also governing the microstructure of AM components. Nevertheless, there is an important difference between welding processes and AM. During AM, several hundred of lines and several hundred of layers form the component over many hours. Thus, we observe transient temperature fields, superposition of temperature fields, in-situ heat treatment and geometric effects (Figure 1). Generally, the material is remolten and reheated many times.

Fig. 1
figure 1

AM temperature fields. (a) Superposition of the temperature field of neighboring lines. (b) Calculated evolution of the temperature for Ti-6Al-4V during E-PBF at different positions at and below the surface. Reproduced with permission.[36] Copyright 2017, Springer. (c) Geometric effects demonstrated using a ring structure

Full size image

Figure 2 gives an overview of important experimental observations. During AM, the previous layer is partially remolten and solidification is very often dominated by epitaxial growth combined with directional solidification in building direction (Figure 2(a), left). Grain selection often leads to columnar grains associated with fiber texture (for cubic crystals the <001>-direction is oriented in building direction).[37] This columnar grain structure tends to show secondary alignment and texture formation within the building plane (x-y-plane).[38,39] The latter is dependent on the scanning strategy and reflects the symmetry of the scanning pattern. Suitable scanning strategies support secondary grain selection in such a way that even single crystalline structures evolve.[14,40,41,42]

Fig. 2
figure 2

Adapted under the terms of the CC BY 3.0 license.[46] Copyright 2013, Elsevier

Experimental observations with respect to microstructure evolution: (a) columnar to equiaxed transition by increasing the energy input (E-PBF, IN718), (b) new grain formation. Left: new grains at the bottom of the melt pool. Right: new grains everywhere due to the addition of TiB2 particles (L-PBF, Al-12Si). Reproduced with permission.[47] Copyright 2019, Elsevier. (c) Influence of the scanning strategy. Left: hatch direction changes in each layer by 90 deg. Right: hatch direction changes only after 10 layers (E-PBF, IN718). Reproduced under the terms of the CC BY 4.0 license.[48] Copyright 2014, EDP Sciences. (d) Geometric effects (E-PBF, Ti-6Al-4V).

Full size image

Very often, new grain formation at the bottom of the melt pool is observed (Figure 2(b), left). Several experimentalists reported but not investigated this phenomenon in detail until today. One reason of grain nucleation is that a segregated microstructure is remolten. During remelting it may happen that some precipitates or particles do not fully dissolve and serve as heterogeneous nuclei (Figure 2(b)). In addition, the rough interface between cells or dendrites and melt pool in combination with a more or less homogenous melt (convection) may lead to a strong undercooling situation, especially within the interdendritic regions. The latter phenomenon is almost unexplored and needs detailed investigations in the future.

Properly understood is the well-known columnar-equiaxed transition (CET) during solidification.[43,44] During CET new nuclei form in front of the solidification front and the columnar grain structure transfers to an equiaxed one. This phenomenon is governed by the temperature gradient \( G \) and the solidification velocity \( v \).[45] The CET is initiated by strong undercooling ahead of the solidification front. Figure 2(a) shows this transition for IN718 (E-PBF) provoked by increasing the total energy input. Nevertheless, standard processing conditions for L-PBF or E-PBF typically do not lead to CETs.

In AM, the scanning strategy has strong influence on the microstructure (Figure 2(c)). In this case, the beam parameters and the energy input are identical. The only difference between these two experiments is that the scanning direction changed by 90 deg for each layer for the equiaxed sample and every 10 layers for the columnar one. This example demonstrates that the change of the main solidification direction in each layer has strong influence on the appearance of new grains. Either more new grain nuclei are formed and/or the present nuclei have better growth conditions due to the change of the solidification direction.

Finally, surface effects have strong influence on the microstructure (Figure 2(d)). For small geometries, e.g., thin walls or lattice structures, surface effects may be even dominant. New grains origin from partially molten powder particles at the surface.[46] In addition, the temperature field near the surface differs from that within the bulk. Thus, the properties of very fine structures are quite different from that of the bulk material.

3 Physical Models and Numerical Methods

During AM, rapid and directional solidification determines the resultant microstructure of the metallic deposit.[46] Speed, power and size of the heat source govern melt pool geometry and solidification kinetics.[49] Microstructures observed[50] are the result of a complex interplay between solidification conditions, thermal environment, thermal cycling, phase transformations and their kinetics. Macroscale components (hundreds of mm in dimension) develop from a submillimeter melt pool formed from powder particles with about 101-102 µm diameter. Hundreds of layers between 20 µm and 500 µm are necessary to realize components whereas the microstructure is on the micron or submicron scale.[51] The time scale for microstructure evolution is 10−3 s or less whereas the total building time takes several hours. These time and length scales across many orders of magnitude constitute an enormous challenge for modeling both, the thermal field as well as the resulting microstructure.

3.1 Temperature Fields

Basis for microstructure evolution is the transient temperature field (Figure 1).[52] The melt pool depth is strongly dependent on the scan length. There is a superposition of the temperature fields of neighboring lines (Figure 1(a)). The melt pool depth normally increases from line to line until a quasi-stationary value is reached. Shorter scan lengths lead to larger melt pools since the beam return time decreases. At turning points, the energy input depends on the AM machine. Typically, the melt pool depth at turning points is larger compared to the bulk. Finally, boundary effects have strong influence on the temperature field since the surroundings (powder bed in PBF, atmosphere in DED) are nearly perfect thermal isolators.

Although the basic equation for the temperature is a parabolic scalar equation, which is relatively easy to solve, many simplifications are still necessary to make the whole problem feasible for microstructure evolution calculations. Generally, individual powder particles are not considered. The powder bed is treated in a homogenized way as second material with adapted properties.[53] Hydrodynamics is in most cases not considered. Exceptions can be found in Rai et al.[54] albeit in two dimensions. Evaporation effects are usually ignored. Even with these strong simplifications, the calculation of predictive temperature fields is still a challenge:

  1. 1.

    The absorbed beam energy and energy loss by evaporation is not exactly known.

  2. 2.

    Many material parameters such as the specific heat or the thermal diffusivity are unknown at high temperatures. These parameters are simply estimated or extrapolated.

  3. 3.

    The temperature field is not stationary but transient. It is influenced by initial effects, turning points, length of scan vectors, etc. (Figure 1(c)).

Whereas the first and second aspects require suitable experimental input, e.g., to adapt the energy absorption coefficient, the last aspect is a question of computational effort. High deflection speeds combined with small melt pool dimensions constitute an enormous challenge on the computational resources when modeling the resulting thermal fields. This leads to the fact that quite different approaches, i.e., analytical or numerical solutions, are applied to generate temperature fields as input for microstructure calculations.

Many researchers address this challenge by developing analytical or quasi-analytical solutions for the transient heat conduction problem. These approaches are based on the fundamental solution of the heat conduction equation, i.e., the solution for a spatial and temporal point-line source.[55,56,57,58] Rosenthal calculated the quasi-stationary thermal field surrounding a moving point source.[55] Eagar and Tsai replaced Rosenthal’s point source with a Gaussian heat distribution.[57] Yajun combined a point- and line source approach for modeling welding melt pools.[58] Analytical solutions are reaching their limit at phase transitions or for real geometries. Promoppatum et al. compared the results of Rosenthal’s solution and a numerical finite element (FE) model.[59] While the authors find good agreement between both models for moderate energy input values, Rosenthal’s solution overestimates the melt pool dimensions in comparison to the numerical results when applied to line energies larger than 0.4 J/mm. One reason for that is that Rosenthal’s equation is not able to treat phase transformations and latent heat effects.

Solving the heat diffusion equation numerically involves the discretization of the underlying governing equation. Besides finite element (FE), the most commonly used approaches are Finite-Difference (FD) and Finite-Volume (FV) techniques.[60] An efficient combination of thermal analysis and fluid dynamics for modeling of PBF provides the Lattice-Boltzmann (LB) method.[61] These methods incorporate the energy input of the laser- or electron beam by an additional source term, in most cases a Gaussian distribution analogous to the work of Eagar and Tsai.[57] Goldak implemented a double ellipsoidal power density distribution for analyzing both shallow and deep penetration welds.[62] A survey of the use of the FE method for welding is given by Lindgren.[63] The small length and time scales encountered in fusion welding enforce a very fine spatial and temporal resolution. Maintaining these extraordinary resolutions for an additively build part is completely unfeasible even for the most capable super computers available today. Thus, King et al. suggest the separation of different scales for generating insights in the AM process.[64]

Summarizing, most approaches in literature to calculate temperature fields as input for microstructure evolution calculations use stationary solutions of the temperature field, i.e., the temperature field is constant in a coordinate system moving with the same velocity as the beam moves. These fields result either from analytical solutions of a moving heat source or numerical approaches. Analytical solutions are cheap with respect to numerical effort but show specific deficits to model the correct melt pool geometry. Numerical approaches for the stationary temperature field result in a more realistic melt pool geometry. Based on these stationary solutions for the temperature field, different process strategies can be simulated resulting in bulk microstructures. Although analytical solutions allow the superposition of temperature fields of neighboring lines, this approach fails for real part geometries such as depicted in Figure 1(c). That is, the full temperature field, which is the result of a specific process strategy and the geometry, is necessary to predict the microstructure of real parts. Köpf et al. already demonstrated the consideration of the full temperature field for simple geometries as input for microstructure simulation.[65]

3.2 Microstructure Modeling

Each model is only able to capture part of the full physics of microstructure evolution. Before we discuss the different microstructural modeling approaches in detail, we review the fundamental processes at the solidification front, which eventually determine the microstructure. This is essential to understand the different model assumptions, their advantages and disadvantages, and eventually their predictive power.

Figure 3 shows schematically the situation at the solidification front. Solidification of alloys is characterized by planar, cellular or dendritic solidification.[45] Grain growth velocity in z-direction is given by the isotherm velocity \( v_{T} \). Undercooling is necessary for each grain to reach \( v_{T} \) in solidification direction (here z-direction). The growth velocity in cubic crystal systems is highest in <001>-direction. Unfavorably oriented grains with angle \( \phi \) between <001> orientation and z-direction have to grow faster with velocity \( v_{\phi } = v_{T} /{ \cos }\left( \phi \right) \) in order to reach the velocity \( v_{T} \) in z-direction. As a consequence, undercooling at unfavorably oriented grains has to be higher than at favorably oriented ones. Thus, unfavorably oriented grains fall behind favorably oriented ones. This mechanism is essential for grain selection since favorably oriented grains are able to overgrow unfavorably oriented ones. As a result, grains oriented with their preferred growth orientation aligned with the thermal gradient will outgrow grains with other orientations and the microstructure will evolve a fiber texture.

Fig. 3
figure 3

Schematic of the situation at the solidification front (in analogy to Gandin and Rappaz[66])

Full size image

A further, very important mechanism is new grain formation within the undercooled zone in front of the solidification front. The size of the undercooled zone depends on the temperature gradient \( G \) and the solidification velocity \( v = v_{T} \). As long as the necessary undercooling for nucleation ΔTn, is higher than the undercooling in front of the solidification front, new grain formation is suppressed. Nucleation of new grains becomes more and more important with decreasing \( G \) and increasing \( v \). Eventually, epitaxial grain growth transforms to an equiaxed microstructure. This is the well-known CET.

Nucleation has a pronounced effect on the resultant microstructure. Generally, phenomenological approaches model nucleation according to experimental observations. Usually, the probabilistic phenomenological approach based on the work of Gandin and Rappaz is used[66]

$$ N = \frac{{N_{0} }}{{\Delta T_{\sigma } \sqrt {2\pi } }}\mathop \smallint \limits_{0}^{\Delta T} e^{{ - \left( {\frac{{\Delta T^{\prime} - \Delta T_{\text{c}} }}{{\sqrt 2 \Delta T_{\sigma } }}} \right) }} {\text{d}}(\Delta T^{\prime}), $$
(1)

where \( N_{0} \), ΔTc and ΔT are nucleation parameters which have to be determined experimentally. Physically, the three parameters have different meanings: N0 is a maximum nuclei density, ΔTc is the critical undercooling necessary for nucleation (it corresponds to ΔTn in Figure 3), ΔT describes a temperature interval where nucleation takes place.

The following sections discuss the microstructure modeling approaches, namely phase-field, cellular automaton and Monte Carlo Potts models. Figure 4 shows schematically the basic idea behind each of these models.

Fig. 4
figure 4

Schematic of the different modeling approaches: Left: Phase-field method, middle: Cellular automaton, right: Monte Carlo Potts model

Full size image

3.2.1 Phase-Field Model

The phase-field (PF) method is an important approach to describe phase transformation induced microstructure evolution.[67,68] PF methods were originally applied to model solidification[69] but the area of application has spread to, e.g., solid-state phase transformation, recrystallization, grain coarsening or heat treatment. The fundamental principle is a smoothly varying function called ‘phase-field’ describing the interface between two mobile phases. One important advantage of the PF method to so-called sharp interface models[70,71,72,73] is the avoidance of explicitly tracking the interface between different phases reducing computational costs.[74,75] The temporal evolution of the PF rests upon a sound thermodynamic and kinetic basis and represents the phase evolution, i.e., microstructure evolution. Thus, PF modeling of alloys is always linked to thermodynamic modeling (CALPHAD).

PF models make use of a scalar field, e.g., \( \phi \in \left[ { - 1;1} \right] \), to distinguish between two phases, e.g., the solid (\( \phi = 1 \)) and the liquid (\( \phi = - 1 \)) phase.[76] Typically, a phenomenological free energy functional \( F \) is applied to describe the two-phase system [76,77]

$$ F\left( {\phi ,c,T} \right) = \mathop \smallint \limits_{V}^{{}} \left[ {f\left( {\phi , c,T} \right) + \frac{{\smallint_{\phi }^{2} }}{2}\left| {\nabla \phi } \right|^{2} + g\left( \phi \right)} \right]{\text{d}}V, $$
(2)

where c is a concentration field, T is the temperature field, V a volume, f is the thermodynamic potential and \( \varepsilon_{\phi }^{2} /2 \cdot \left| {\nabla \phi } \right|^{2} \) corresponds to the interfacial energy with parameter \( \varepsilon_{\phi } \) and g is a potential function.

The thermodynamic potential f is an interpolation function between the free energy densities of the various phases. The construction of this interpolation function is a crucial step for setting up a PF model. On the one hand, the bulk free energy densities have to be consistent with values obtained from thermodynamic databases.[76] On the other hand, there is a freedom in the choice of the interpolation function because it only has an influence on the interface region because of its dependence on \( \phi \).[76] The role of the potential function g, usually the double well potential, is to stabilize the two phases \( \phi \; = \; \pm \;1 \).[78] The combination of \( f \) and\( g \) is the so-called free energy density function.[77]

During any process, the free energy functional decreases monotonously. The equations of motion describing the time evolution of the system are gained by calculus of variations of the functional \( F \).[67] The first equation is the Allen–Cahn equation

$$ \frac{\partial \phi }{\partial t}\; = \; - M_{\phi } \frac{\delta F}{\delta \phi }, $$
(3)

which describes the time evolution of the phase field. The parameter \( M_{\phi } \) is a mobility related to the interface kinetic coefficien[77] which is usually temperature dependent.[76] The second equation is the so-called Cahn–Hilliard equation

$$ \frac{\partial c}{\partial t} = \nabla \cdot \left( {M_{\text{c}} \nabla \frac{\delta F}{\delta c}} \right), $$
(4)

where Mc is the mobility of solute atoms. The Cahn–Hillard equation describes the process of phase separation. There also exist multi-component and multi-phase PF models applying several phase-fields \( \phi_{i} \in \left[ {0;1} \right] \), which sum up to unity.[68,77] These models allow more than one solid phase, e.g., to model precipitations in the interdendritic regions.[79]

The PF method describes the evolution of the phase and the concentration field by using a coupled temperature model. The simplest coupling is to apply defined temperature conditions, like a constant undercooling, cooling rate or temperature gradient. The assumption of defined temperature conditions is only valid for the simulation of few dendrites on short time scales. Additionally, a quasi-stationary evolution of the solidification front in the melt pool is required. The defined values are often gained by analytical solutions or macroscopic simulation results of temperature fields and are provided within a small subset as boundary conditions.[80,81,82] A second case is the direct coupling to a numerical temperature field which is computed parallel to the PF.[83] The option to couple the concentration fields back to the temperature solver is usually not used. Sometimes the coupling is extended by melt pool motion, i.e., the melt step is additionally solved (usually on the powder scale) and the melt pool dynamics are coupled to the PF model.[83,84]

PF models require many different material as well as model parameters. The applied material parameters are often unknown, such as exact temperature dependent mobility coefficients. However, a reasonable approximation is due to the relation to real physics often possible. In addition, atomistic simulations are suitable to gain material parameters like kinetic coefficients.[29,85,86] It is important to notice, that some model parameters used in the free energy function can arbitrarily be chosen because they are not related to physical properties. One example is the interface width, which can have an influence on the numerical result. Exemplarily, Fallah et al. show, that the interface width has to be small enough to achieve a correspondence to experimental results.[80] On the one hand, these additional model parameters allow an exact calibration of the PF model to experimental results. On the other hand, a careful validation of the model after calibration with more experimental results is necessary to ensure correct numerical predictions.

PF models are usually applied on the dendritic scale (Figure 5(a)), i.e., the spatial resolution is of the order of nanometers to resolve dendritic structures and concentration fields. The maximum domain size is therefore limited to the order of micrometers due to computational restrictions. In the recent years, numerical approaches like mesh refinement and massively parallelization improved the computational efficiency of PF methods.[85] However, the necessary fine resolution near the dendritic tips makes an application on the part scale impossible.

Fig. 5
figure 5

Adapted with permission.[98] Copyright 2018, Laboratory for Freeform Fabrication and University of Texas at Austin

PF simulation results. (a) Solute concentration field of rapid solidification on the dendritic scale (resolution of 30 nm) of Ti-6Al-4V with a cooling rate of 50,000 K/s. Reproduced with permission.[83] Copyright 2019, Elsevier. (b) Phase field information on the grain scale (resolution 100 nm) for E-PBF of Ti-6Al-4V with 420 W beam power and 0.8 m/s scan speed.

Full size image

The dendritic scale is suitable to predict the primary dendrite arm spacing for different materials under different process conditions.[80,81,82,83,87,88] The numerical predictions are in good agreement with experimental results as well as with phenomenological equations (e.g., Hunts model[89]). Due to the fine spatial resolution, the concentration in the dendritic core as well as segregations in the interdendritic regions are resolved. This allows the relation between process parameters and the resulting changes in the dendrite morphology as well as in the concentrations.[81,90]

A further application is to identify growth rates of the dendritic tips.[29,82,91,92] The identified relations are very useful for mesoscopic interface tracking approaches (CA-Model in Sect Cellular Automata). Once multi-phase PF methods are applied, it is also possible to predict additional phases during solidification like precipitates in the interdendritic regions.[79]

The fine spatial resolution of the PF approach on the dendritic scale prevents the use on the part scale. In order to bridge these scales, there are two main approaches. One approach is the dendritic needle network model, where a hierarchical network of needle-like branches interact through the long-range diffusion field.[85,93,94] The computational effort decreases due to the simplified modeling of dendrites by needles. This approach has already been used for directional solidification,[94] but a validation for rapid solidification processes is missing. A different approach is the PF method on the grain scale, where the free energy functional is modified, that whole grains instead of single dendrites are modeled (Figure 5(b)).[84,95,96,97,98] Each PF represents an individual grain. The number of grain orientations is limited by the number of different PF variables.[97] As a consequence to the high number of PFs and the coarser resolution, the evolution of the concentration field is no longer tracked. The approach mimics different grain orientations by different mobility constants for different grains.[84] But, these constants are not dependent on the temperature gradient, i.e., independent from the favorable growth direction, a certain grain is always preferred due to a higher growth velocity. Although this approach is already applied on AM,[84,98] a broad experimental validation is missing.

PF models on the grain scale increase the order of the spatial resolution to micrometers. Consequently, they can operate on the same scale than CA approaches and may reach similar results. However, the PF method needs more model parameters, which implies a more careful validation and consequently raises issues regarding numerical predictability. In addition, regarding computation time, CA models are more efficient.[99]

Summarizing, on the dendritic scale, the PF method is important to bridge between atomistic and mesoscopic models and to predict the microstructure of few dendrites. Regarding the grain scale, we recommend the application of CA instead of PF approaches because of less model parameters, a broader experimental validation and a better numerical efficiency.

3.2.2 Cellular Automata

Gandin et al. employed successfully a CA model coupled with an FE heat model for the prediction of the microstructure in complex castings (CAFE model).[100] The challenge to use this model for AM is melting and solidification of hundreds of layers with hundreds of lines with extremely varying temperature fields. On the one hand, this is a question of computational power. On the other hand, and this is different from castings, a specific fraction of material is remolten during AM.

CA were introduced by von Neumann.[101] The physical system of interest is divided into cells, where each cell is allowed to interact iteratively with its neighbors according to some specific rules.[102] The strength of a CA model in comparison to Potts Monte Carlo- or Ginzburg-Landau type PF models lies in the combination of computational simplicity and scalability with physical stringency.[103]

The fundamental idea of the model lies in reproducing the envelope of a growing dendritic grain by the superposition of many small geometric objects (squares in 2D, octahedrons in 3D), each controlled individually by a cell of the CA (Figure 6). A serious issue of modeling solidification microstructures with CA lies in their tendency to reflect the underlying grid network.[28] Gandin and Rappaz solved this problem by developing a sophisticated CA crystal growth model for cubic metals both in two- as well as three dimensions.[25] The model avoids the lattice dependency by tracking the dendritic growth front at the level of the typical secondary arm spacing (decentered square model).

Fig. 6
figure 6

Schematic of the basic CA approach to simulate grain growth. (a) 2D-dendritic growth at T = const. with a square as envelope, (b) 2D-dendritic growth within an temperature gradient G and envelopes, (c) representation of the envelope by squares defined for each cell forming the envelope, (d) 3D-representation of a dendrite (Adapted under the terms of the CC BY 3.0 license.[104] Copyright 2012, IOP Publishing Ltd) by an octahedron and representation of a 3D-grain

Full size image

Based on observations by Ovsienko et al.[105] of cyclohexanol crystals exhibiting square contours when growing in a uniform thermal field, the two-dimensional model approximates the envelope of a growing dendrite by a square. The vertices of the square represent the tips of the primary dendrite arms growing along the dendrites main growth directions (<10> for cubic metals, Figure 6(a)). The orientation of the dendrite is considered by the rotation of primary growth directions about an angle\( \varTheta \) against the lattice coordinate system.

For grain nucleation in a distinct cell, a square is positioned with its center matching the center of the nucleation cell. Subsequently, the square grows with velocity \( v\left( {\Delta T} \right) \) depending on the local undercooling \( \Delta T \)of the cell. With the temperature assumed constant throughout the cell, growth of the square is realized by moving all vertices along the respective diagonals of the square with the same velocity \( v \). Once a square reaches the center of an adjacent liquid cell, this cell is incorporated into the growing dendrite. The newly initiated square matches the grain number and -orientation of the parental grain. It is truncated and repositioned to ensure the envelope of the grain remains intact according to distinct rules explained in detail by Gandin and Rappaz.[25] Due to this repositioning, the algorithm is denoted decentered squares (DS).

Different temperatures of the cells in the presence of a thermal gradient \( G \) result in varying growth velocities of the corresponding squares. This growth results in a distortion of the grain growing faster in the direction of increasing undercooling (Figure 6(b)). Growth of the dendrite is accelerated in the direction opposing the direction of \( G \) while maintaining its preferred growth orientation. Figure 6(c) illustrates the envelope of the dendrite calculated by this algorithm. The active cells (red) at the border of the dendrite form the envelope with their corresponding squares (blue).

The three-dimensional analogy of the DS algorithm is denoted decentered octahedra (DO). The algorithm is quite similar to the two-dimensional case. Octahedra are growing depending to the undercooling of the associated (cubic) cells. When the envelope of an octahedron comprises the center of an adjacent cell, the cell is incorporated in the growing grain. Figure 6(d) shows an octahedron enveloping a three-dimensional dendrite and the result of the DO algorithm for a three-dimensional dendrite growing within an undercooled melt.

As mentioned above, the growth velocity \( v\left( {\Delta T} \right) \) of the square or octahedron associated with a distinct cell is a function of the cells undercooling. The growth velocity correlates the CA model with the thermodynamics and kinetics of the liquid-solid phase transformation. For this reason, the choice of \( v\left( {\Delta T} \right) \) is crucial. In their original work, Gandin and Rappaz[28] applied a model form Kurz, Giovanola and Trivedi[106] (KGT) for dendritic growth rates in the range of the limit of absolute stability for calculation of the dendritic growth kinetics. Foundation of this analysis represents the theory of morphological instability. This theory constitutes a dynamical approach on calculating the stability of perturbed surfaces in contrast to the static approach used by the theory of constitutional supercooling.[107] Consideration of the complex relationships addressed in the KGT model for prediction of crystal growth using a CA at the mesoscale is computationally expansive. Gandin et al.[26,108] use a simple polynomial function to speed up the calculations:

$$ v\left( {\Delta T} \right)\; = \;a_{2} \cdot \left( {\Delta T} \right)^{2} + a_{3} \cdot \left( {\Delta T} \right)^{3} , $$
(5)

where \( a_{2} \) and \( a_{3} \) are parameters to be defined. Other authors, e.g., Guillemot et al.,[109] use the form

$$ v\left( {\Delta T} \right) = A \cdot \left( {\Delta T} \right)^{n} , $$
(6)

where the parameters \( A \) and \( n \) have to be adapted to the investigated material. The variety of implemented approximations shows the arbitrariness of the solutions. As presented, a reasonable approach for modeling the growth kinetics at the microscale is the PF method.

Macroscopic heat conduction calculations usually provide the local undercooling ΔT necessary for the calculation of the dendrite growth velocity. For simple scan paths and geometries, analytical solutions for the heat conduction equation were successfully employed for microstructure simulations by Köpf et al. for E-PBF.[110] Zinovieva et al. applied Goldak’s heat source model with CA for calculating the three-dimensional microstructure evolution for L-PBF.[111] While analytical methods are computationally efficient, their use is limited to model bulk materials microstructures. Modeling complex beam paths of real components requires numerical solutions. A combination of FE heat analysis and CA crystal growth modeling is shown by Koepf et al.[65] In this work, the authors resolved the individual scan paths and the resulting thermal field for an additively build cylinder (160 layers, 30 lines in each layer) and validated the result with optical micrographs (Figure 7). Li and Tan combined a mesoscale CA with a macroscale thermal model for to investigate the influence of nucleation in direct laser deposition processes.[33] By adopting the nucleation model originally introduced by Rappaz et al.,[28] the authors found the parameters used in the nucleation model are significantly affecting the resulting grain structure, albeit without experimental validation.

Fig. 7
figure 7

Reproduced under the terms of the CC BY 4.0 license.[65] Copyright 2012, Elsevier

3D section view of the numerically predicted microstructure (left) and comparison of the simulation result with a longitudinal section cut of the superalloy CMSX-4 (right).

Full size image

The repeated beam scanning during processing of many layers remelts previously consolidated material. Modeling of this phenomenon in CA necessitates the dynamical adjustment of all octahedra constituting the partly remolten grains. There is a high uncertainty how to describe this effect within the CA model. Köpf et al introduced a factor for dynamical size adaption of the grains, but kept it at unity for the time being.[110]

3.2.3 Potts Model

The Potts model shows the highest level of abstraction with respect to the physical processes. Dendrites, segregation and undercooling during solidification are not considered. Figure 4, right, shows a schematic of the approach for modeling grain structure evolution during solidification. All cells within the melt pool have a different state (“spin”). During solidification, the structure starts coarsening due to surface energy effects. Coarsening only talks place within a heat affected zone (HAZ). At lower temperatures, i.e., outside of the HAZ, coarsening stops and the grain structure is stable.

The Potts model is a generalization of the Ising model, which describes the interaction of spins on a lattice.[112,113] The strength of Potts model is that it can be generalized to model different problems in physics, e.g., coarsening in grain structures or in foams. Potts models to simulate grain evolution during AM are derived from the Potts Kinetic Monte Carlo (KMC) method.[114] The KMC method aims to simulate the time evolution of processes based on given transition rates among states. The problem is defined on a lattice and the driving force for grain boundary movement and coarsening is the curvature of the interface between different grains. Each lattice site is assigned a state, the so-called “spin”. Neighboring lattice sites with the same spin form a grain. The total energy of the system is the sum of the interaction energy between neighboring lattice sites with different spin[34]

$$ E\; = \;\frac{1}{2}\mathop \sum \limits_{i = 0}^{N} \mathop \sum \limits_{j = 1}^{L} \left( {1 - \delta \left( {q_{i} , \;q_{j} } \right)} \right), , $$
(7)

where \( N \) is the number of lattice sites, \( L \) is the number of neighbors at each site, \( q_{i} \) is the spin at lattice site \( i \), \( \delta \) is the Kronecker delta function which is 1 when \( q_{i} = q_{j} \) and 0 otherwise.

Curvature driven changes of the spin of each lattice site lowers the total energy of the system. A random process changes the spin of a lattice site to that of a neighboring site. As a result, the total energy changes. Whether this new spin is accepted is decided with the help of the Metropolis algorithm where a random number between 0 and 1 is compared with an acceptance probability P, given by[34]

$$ P = \left\{ {\begin{array}{*{20}c} {e^{{ - \frac{\Delta E}{{k_{\text{B}} T_{\text{s}} }}}} ,} & {\Delta E > 0} \\ {1,} & {E \le 0} \\ \end{array} } \right. $$
(8)

where \( \Delta E \) the energy change of the system, \( k_{B} \) is Boltzmann’s constant, TS is the simulation temperature.

If the energy of the system is reduced, the new state is always accepted. If the energy is increased, the new state is accepted with the probability P. During one time step, all lattice sites are updated. It is important to note that the simulation temperature Ts does not represent the real temperature of the system. The energy \( k_{\text{B}} T_{\text{s}} \) defines the thermal fluctuations present in the simulation.

Rogers et al. use the KMC method to simulate microstructure evolution during AM.[115] The heat source is not directly simulated but imposed as a molten zone surrounded by a heat affected hone (HAZ). The molten zone is simulated by assigning random spins to lattice sites within the melt pool. Local grain boundary mobility is only simulated within the HAZ. The grain boundary mobility M within the HAZ surrounding the melt pool is given by

$$ M\left( T \right) = M_{0} e^{{ - \frac{Q}{RT}}} , $$
(9)

where Q is the activation energy for grain boundary movement, R is the gas constant, M0 is a constant. With this, the probability of executing a spin change becomes[34]

$$ P = \left\{ {\begin{array}{*{20}c} {M\left( T \right) e^{{ - \frac{\Delta E}{{k_{\text{B}} T_{\text{s}} }}}} ,} & {\Delta E > 0} \\ {M\left( T \right), } & {\Delta E \le 0} \\ \end{array} } \right., $$
(10)

Small grains just at the boundary between the melt pool and the HAZ appear. These grains are quickly consumed by larger grains or grow to from larger grains for themselves.

Using KMC model, grain structure formation takes place within the HAZ and is governed by surface curvature effects inducing grain coarsening. If we compare this mechanism to generate grain structures with the real physical grain growth and selection process (Figure 3), it becomes clear, that this approach does not cover the fundamental physical mechanism. Solidification is a highly non-equilibrium process governed by the kinetics of diffusion and direction dependent growth kinetics. The complex dendritic structure demonstrates that interfacial energies do not govern the grain selection process during solidification. Consequently, the current KMC approach is not able to model grain selection and texture formation. There is some similarity between experimental grain structures and simulated structures. This similarity results from the direction of solidification. Li and Soshi compare aspect ratio, grain orientation and size with experimental results.[116] Nevertheless, the Potts Model does not represent the mechanism of grain coarsening by competitive grain growth. The same is true for texture formation.

In the current form, the KMC method is not able to predict correctly AM grain structures. Nevertheless, the KMC method is an interesting approach to model in-situ heat treatment effects within the solidified microstructure, i.e., grain coarsening in the HAZ.

3.2.4 Influence of the Nucleation Parameters

Modeling of nucleation of new grains is crucial for microstructure evolution. It is important to note, that Eq. [1] is not derived from fundamental physics but is a phenomenological description. Thus, the three parameters have to be adapted for each alloy under investigation. This is a practicable approach as long as one parameter set is used for one specific alloy for all experimental situations. Actual state-of-the-art is that the influence of different nucleation parameters are numerically investigated by different groups.[33,117] Figure 8 demonstrates the influence of the nuclei density \( N_{0} \) and critical undercooling Tc on the resulting grain structure for a CA model for DED.

Fig. 8
figure 8

Reproduced with permission.[117] Copyright 2019, IOP Publishing Ltd

Cross-sectional views of the SD-BD plane for each microstructural domain obtained by varying N0 and ΔTc. (a) Domain A: N0 = 1013 m−3, ΔTc = 10 K; (b) domain B: N0 = 1014 m−3, ΔTc= 0K; (c) domain C: N0 = 1015 m−3, ΔTN = 5 K; (d) domain D: N0 = 1015 m−3, ΔTN = 0 K. The black dashed line indicates the location of the fusion line from the first layer.

Full size image

This example demonstrates a fundamental problem. Nearly all kind of results can be generated by adapting the nucleation parameters. A sound correlation between experimental observations (Figure 2) and nucleation parameters is still missing and one of the important tasks for the future.

4 Summary

The prediction of the microstructure of AM components is essential for a broad industrial use of this technology in the future.

The starting point for the calculation of the microstructure is the temperature field. Bulk microstructures can be predicted based on analytical or quasi-analytical models for the temperature. Nevertheless, to capture all effects involved with real parts, such as surface effects, turning points, etc., simplifications of the temperature field are unrewarding.

Phase field (PF) models have a sound physical basis including thermodynamics and kinetics of microstructure evolution. Complex solidification microstructures such as cells or dendrites including segregation effects develop from this approach on the dendritic scale. Thus, PF models represent the gold standard among all modeling approaches for the microstructure. Nevertheless, PF models are extremely CPU-intensive. Thus, only a small number of dendrites is accessible. Even 1000 dendrites is still small for rapid solidification conditions where dendrite arm spacing are of the order of microns or even less. Thus, microstructure calculations on the scale of real parts will not be accessible in the near future. In addition, one should not forget that PF models are based on a variety of uncertain data. This concerns thermodynamic as well as kinetic data. Different databases lead to different results. Thus, reliable predictions by PF models are not invariably the case. Nevertheless, PF models bring important insights into the situation at the solidification front, e.g., the correlation between undercooling and dendrite tip velocity. In addition, segregation effects at rapid solidification become more transparent. This is very important for the appearance of certain phases[88] or the evolution of hot cracks.[118]

The cellular automaton (CA) approach describes microstructure evolution on the scale of the grains. The convex envelope of cells or dendrites is modeled whereas cellular or dendritic structures and segregations are not captured. Thermodynamics comes into the model by the correlation of the dendrite tip velocity and the local undercooling. CA models have the potential to predict grain structure evolution on part dimension. They also have a high power to predict texture formation during AM and its dependence on the process strategy.

The Kinetic Monte Carlo model in its current form is not able to predict grain structure evolution since the underlying mechanism is not representative for the real one. Grain coarsening is treated within the KMC model not by competitive grain growth but is driven by the curvature of the grain boundaries. Texture is not predicted and even the appearance of the grain structure is not reproduced. Thus, the KMC model is not recommended for AM grain structure evolution. Nevertheless, KMC models might be very useful to model grain coarsening effects during AM induced by in-situ heat treatment.

One essential problem for microstructure simulation, which is independent from the model, is new grain formation by nucleation. In principle, by adapting the nucleation parameters nearly all kind of grain microstructures evolve. This is a fundamental problem because it weakens the predictive force of all models. Thus, one task for the future is to determine the nucleation parameters for a variety of alloys based on experimental data. This would help to assess the predictive force of microstructure evolution models.