A mesoscopic digital twin that bridges length and time scales for control of additively manufactured metal microstructures

The analytical model (i.e. Lipton–Glicksman–Kurz (LGK) model) developed by Lipton et al allows determination of a unique pair of $\left( {V,{ }{R_{{\text{tip}}}}} \right)$ for a given undercooling $\Delta T$, where $\Delta T \equiv {T_{{\text{liq}}}} - T$ with ${T_{{\text{liq}}}}$ and ${R_{{\text{tip}}}}$ being liquidus temperature and radius of curvature of the dendrite tip, respectively [32]. The total undercooling $\Delta T$ consists of three contributions: the thermal ($\Delta {T_{\text{T}}}$), solutal ($\Delta {T_{\text{C}}}$), and curvature ($\Delta {T_{\text{R}}}$) undercoolings. It is worth mentioning that the undercooling associated with attachment kinetics is neglected in the model. This assumption is thus valid for small undercooling. An expression for $\Delta T$ is given as follows:

$$\begin{equation}\Delta T = \Delta {T_{\text{T}}} + \Delta {T_{\text{C}}} + \Delta {T_{\text{R}}} = \frac{{{L_{\text{f}}}}}{{{c_{\text{p}}}}}{\text{I}}{{\text{v}}_{{\text{3D}}}}\left( {{\text{P}}{{\text{e}}_{\text{T}}}} \right) - {m_{\text{l}}}{C_0}\left[ {\frac{{\left( {1 - {k_0}} \right){\text{I}}{{\text{v}}_{{\text{3D}}}}\left( {{\text{P}}{{\text{e}}_{\text{C}}}} \right)}}{{1 - \left( {1 - {k_0}} \right){\text{I}}{{\text{v}}_{{\text{3D}}}}\left( {{\text{P}}{{\text{e}}_{\text{C}}}} \right)}}} \right] + \frac{{2{{{\Gamma }}_{{\text{sl}}}}}}{{{R_{{\text{tip}}}}}},\end{equation} \tag{ 1 }$$

where ${L_{\text{f}}}$ and ${c_{\text{p}}}$ are the latent heat and the specific heat, respectively, ${m_{\text{l}}}$ is the liquidus slope (assumed to be a constant), ${C_0}$ is the alloy composition, ${k_0}$ is the equilibrium partition coefficient (solute trapping effects are ignored), ${{{\Gamma }}_{{\text{sl}}}}$ is the Gibbs–Thomson coefficient (namely the ratio of the solid/liquid interface energy ${\gamma _{{\text{sl}}}}$ to the melting entropy $= {\gamma _{{\text{sl}}}}{T_{\text{f}}}/{\rho _{\text{s}}}{L_{\text{f}}}$), ${\text{P}}{{\text{e}}_{\text{T}}}$ is the thermal Peclet number: ${\text{P}}{{\text{e}}_{\text{T}}} = V{R_{{\text{tip}}}}/2{\alpha _{\text{l}}}$ with ${\alpha _{\text{l}}}$ the thermal diffusivity, and ${\text{P}}{{\text{e}}_{\text{C}}}$ is the solutal Peclet number: ${\text{P}}{{\text{e}}_{\text{C}}} = V{R_{{\text{tip}}}}/2{D_{\text{l}}}$ with ${D_{\text{l}}}$ the liquid diffusion coefficient. The Ivantsov function in 3D is defined as ${\text{I}}{{\text{v}}_{{\text{3D}}}}\left( {{\text{Pe}}} \right) \equiv {\text{Pe}}\exp \left( {{\text{Pe}}} \right){{\text{E}}_1}\left( {{\text{Pe}}} \right){ }\left( {{\text{Pe}} = {\text{P}}{{\text{e}}_{\text{T}}},{\text{ P}}{{\text{e}}_{\text{C}}}} \right)$, in which ${E_1}\left( {{\text{Pe}}} \right)$ is the exponential integral.

By assuming that ${R_{{\text{tip}}}}$ can be estimated from ${\lambda _{{\text{min}}}}$ derived in the linear stability analysis of a planar solidification front (namely ${R_{{\text{tip}}}} \approx {\lambda _{{\text{min}}}}$ with ${\lambda _{{\text{min}}}}$ being the minimum wavelength of instability), Lipton et al derived an expression for the tip radius

$$\begin{equation}{R_{{\text{tip}}}} = 4{\pi ^2}{{{\Gamma }}_{{\text{sl}}}}{\left[ {\frac{{{\text{P}}{{\text{e}}_{\text{T}}}{L_{\text{f}}}}}{{{c_{\text{p}}}}} - \frac{{{\text{P}}{{\text{e}}_{\text{C}}}{m_{\text{l}}}{C_0}\left( {1 - {k_0}} \right)}}{{\left\{ {1 - \left( {1 - {k_0}} \right){\text{I}}{{\text{v}}_{{\text{3D}}}}\left( {{\text{P}}{{\text{e}}_{\text{C}}}} \right)} \right\}}}} \right]^{ - 1}}.\end{equation} \tag{ 2 }$$

For a given $\Delta T$, equations (1) and (2) can be numerically solved using e.g. the Newton–Raphson method to find a unique pair of $\left( {V,{ }{R_{{\text{tip}}}}} \right)$.

The LGK model outlined above loses validity at large undercoolings and thus may be inaccurate under some AM conditions. PF simulations provide an alternative that imposes fewer approximations and restrictions on the solidification process, permitting more accurate predictions of $V$ versus $\Delta T$ for CA simulations at greater computational cost than the LGK model. Figure 2 shows a comparison of PF and LGK results for Ti–50Nb (in at. %) [33]. The agreement is very good for $\Delta T\, \lesssim \,\Delta {T_{\text{f}}}\, = \,{T_{{\text{liq}}}}\, - \,{T_{{\text{sol}}}}$, where ${T_{{\text{sol}}}}$ is the solidus temperature, but the LGK model significantly overpredicts $V$ relative to the PF results at larger undercoolings. Such differences may propagate into CA-predicted microstructures, and in such cases PF predictions will typically be more justifiable on physical grounds.

Zoom In Zoom Out Reset image size

The PF simulations used to produce the data in figure 2 were implemented within the adaptive mesh phase evolution (AMPE) code. Details of the implemented model can be found in prior publications [34, 35]. In general, this model permits simulating solidification microstructure under wide ranges of cooling rates and thermal gradients, which are able to cover casting and AM conditions. It is important to note here that the model is parameterized to account for phase-, temperature-, and composition-dependent Gibbs free energies and atomic mobilities based on information from the CALPHAD framework. This PF–CALPHAD integration is key for quantitatively simulating thermally or chemically driven-microstructure evolution. The AMPE code implementing this integration has been successfully demonstrated for rapid solidification of Cu–Ni [36] and Ti–Nb [37] alloys.

The grain topology, size, and crystallographic texture of AM metal are determined by local melting and solidification behavior during a laser scan. Accordingly, it is critically important for the CA model to be informed by accurate prediction of local temperature evolution during the laser scan [24–27]. The AM solidification proceeds via two key mechanisms: (a) epitaxial growth from partially melted pre-existing grains in the substrate; and (b) nucleation and growth of new grains in the solidifying melt pool. Both mechanisms are well described within the CA modeling formalism. Details of the associated algorithms to account for fundamental processes leading to solidification microstructures can be found in supplementary material S1 (available online at stacks.iop.org/JPMATER/4/034012/mmedia).

The CA method requires two major inputs: (a) the time-dependent spatial thermal profile; and (b) dendrite tip velocity as a function of undercooling. Within our DT framework, the thermal profile evolution is simulated by the high-fidelity hydrodynamics model for melt pool dynamics implemented within the ALE3D code. In particular, our model accounts for effects of laser beam shapes (e.g. circular Gaussian (CG) beam, transverse elliptical (TE) beam, longitudinal elliptical beam, etc [26, 29, 38]) for L-PBF AM. The simulated profile on the coarse finite element (FE) mesh is projected onto the finer CA mesh employing an interpolation scheme based on the 'scatteredInterpolant' function in MATLAB. The dendrite tip velocity can be computed from either the analytical model (LGK model) or quantitative PF simulations as detailed in section 2.1. We present two examples based on these different models for calculating the dendrite tip velocity.

For the first example, we used a power law function fit to the LGK model results, which were informed by the materials parameters for 316L stainless steel [39]. A CA grid based on a regular lattice of cubic cells is defined and superimposed onto the irregular FE mesh. In this work, we take examples for two beam shapes, CG and TE. The sizes of the CA grids are $700\,\Delta x \times 250\,\Delta y \times 134\,\Delta z$ for the CG beam and $700\,\Delta x \times 300\,\Delta y \times 134\,\Delta z$ for the TE beam with a grid size of $\Delta x = \Delta y = \Delta z = 2{ }\mu {\text{m}}$, which correspond to total volumes of $1.4 \times 0.5 \times 0.268{\text{ mm}}$ and $1.4 \times 0.6 \times 0.268{\text{ mm}}$, respectively. The cell size of $2\,\mu {\text{m}}$ is small enough to resolve the branching of dendrite arms and grain competition [30]. The initial substrate grain structure is fed to the CA model using DREAM.3D [40] informed from an experimental characterization of the microstructure. Grain size distribution and crystallographic texture information were experimentally collected by electron backscatter diffraction and used as an input for DREAM.3D to generate a statistically equivalent synthetic grain structure in the substrate. The laser power and scanning speed were chosen to be $350$ W and $750$ mm s⁻¹, respectively. For the nucleation model, we employed the maximum nucleation density ${N_0} = {10^{15}}\;{{\text{m}}^{ - 3}}$ and mean (or critical) undercooling $\Delta {T_{\text{N}}} = 10\,\;{\text{K}}$. Figure 9 (left, top) shows an example of simulated AM microstructures for 316L stainless steel melted using the TE beam. These digital microstructures for the two beam shapes served as input microstructures for the mechanical response modeling discussed in section 2.4.

Our second example utilizes the simulated liquid–solid interface dynamics in Ti–50Nb computed by PF simulations (see figure 2) for parameterizing the dendrite tip velocity for a given undercooling. As a demonstration, we employed two scanning speeds: 500 and 200 mm s⁻¹ for the given laser power 60 W and the CG beam shape. Single-track CA simulations for Ti–50Nb in $700\,\Delta x \times 250\,\Delta y \times 134\,\Delta z$ grids (with the cell size of $2\,\mu {\text{m}}$) are shown in figure 3. For both processing conditions, grains tend to grow epitaxially, resulting in columnar grains. However, our simulations indicate that grain texture is affected by the laser scanning speed.

Zoom In Zoom Out Reset image size

To capture more detailed microstructural features during solidification, we also combined thermal modeling (ALE3D melt pool model) and PF modeling (by the AMPE code). Here we demonstrate two such examples for Ti alloy solidification. As the first example, a larger-scale multiphysics/thermally-informed PF modeling scheme was applied to simulating Ti–50Nb AM solidification microstructure formation as shown in figure 4 [33]. The same model parameters for the Ti–Nb system as introduced in section 2.1 were used for these simulations. Specifically, PF simulations of directional solidification were used to construct comprehensive maps of steady state solidification mode, cell/dendrite spacing, and microsegregation in $G - V$ space (figure 4(a)). ALE3D simulations of single-track laser scans on bare Ti–50Nb plates were also conducted to predict thermal gradients along the melt pool solidification front (figure 4(b)). By assuming that the corresponding experiments involve similar thermal gradients (∼$1 \times {10^7} - 4 \times {10^7}$ K m⁻¹), PF predictions of steady-state directional growth characteristics could be compared with single-track experiments, and the simulations were found to provide fair-to-good predictive accuracy in most respects. Figure 4(a) also shows the corresponding solidification mode boundaries predicted by the Kurz–Giovanola–Trivedi (KGT) theory of directional dendrite/cell growth [41]. Agreement between analytical KGT predictions, PF simulations, and experimental results is good.

Zoom In Zoom Out Reset image size

For more direct microstructure-level simulations of laser scans, the thermal profiles computed with ALE3D were also used as inputs for 2D PF simulations of full tracks with time-dependent thermal/laser driving (figure 4(c)). These simulations describe microstructure formation across >100 μm track lengths with sub-10 nm resolution (the scale of the liquid–solid interface). Single crystal and polycrystal materials were studied to isolate and characterize effects (e.g. from crystalline anisotropy and competitive grain growth). Coupled PF and multiphysics approaches of this type could soon become reliable predictive tools for tailoring microstructure in the context of both AM process design and AM alloy design strategies [33].

We also considered Ti64 (Ti–6Al–4V in wt.%) as the second example, utilizing a similar integrated modeling approach. Rapid solidification microstructures of Ti64 were simulated by the AMPE code informed by a Ti64 CALPHAD database for temperature- and composition-dependent thermodynamic free energies and atomic diffusion mobilities. Like the case for Ti–Nb, the thermal boundary condition imposing the thermal gradient was informed by ALE3D melt pool dynamics simulations. More specifically, Ti64 melt pool dynamics simulations were performed for the following laser parameters: power 200 W, scan speed 1.3 m s⁻¹, beam size (D4σ) 57 µm. The predicted temperature profiles are shown in figure 5, from which time-dependent boundary conditions were extracted for AMPE simulations.

Zoom In Zoom Out Reset image size

To quantify the underlying thermodynamic forces driving phase stability and transformation within the AMPE code, we developed a comprehensive CALPHAD thermodynamic database including all phases for the Ti–Al–V system based on assessment of the relevant binary databases available in the literature: Al–Ti [42]; Al–V [43]; and Ti–V [44]. The developed database well reproduces the literature-reported property diagram (figure 6(a)) and isopleth section (figure 6(b)), demonstrating the accuracy of the phase-dependent thermodynamic inputs for simulating solidification (liquid to BCC) microstructures in the AMPE code. In addition, a database for composition-, temperature-, and phase-dependent atomic diffusion mobility was developed for the Ti–Al–V system and implemented within the AMPE framework. In particular, the BCC atomic mobilities were taken from [45–49].

Zoom In Zoom Out Reset image size

Due to the lack of information, a constant diffusivity (e.g. 10⁻⁹ m² s⁻¹) is usually assumed in the liquid phase for all the elements within the conventional approach. However, as pointed out by Liu et al [50], the determination of diffusion coefficients of elements in liquid metals is of vital importance when investigating metallurgical processes involving molten metals. Thus, Liu et al [50] developed a predictive diffusion coefficient equation, aimed at describing the diffusivity of solute atoms in liquid metals in the case where experimental data are scarce. We adopted this approach to compute the self- and impurity diffusion coefficients in the liquid phase as reported in figure 6(c).

The integrated AMPE code allows us to run quantitative PF simulations of rapid solidification of Ti64 under relevant AM processing conditions. Figure 7 shows an example of the simulated solidification front of three separate grains, exhibiting planar (no interface instabilities, see figure 7(a)) and partition-less (no segregation of elements, see figure 7(b)) solid–liquid interfaces. In addition, as shown in figure 7(c), the interface migration rate increases with time, from the initial bottom ($v$ = 0 m s⁻¹) towards the top of the melt pool ($v$ = 1.4 m s⁻¹ ≈ V_beam). The characterized maximum rate slightly exceeds the laser scan speed. No grain selection, coarsening, or merging has been observed in the simulations. No nucleation of equiaxed grains is expected because of the very small constitutional undercooling in front of the solid/liquid interface and the high thermal gradient in the liquid phase. Overall, our simulation indirectly indicates that the underlying substrate microstructure should directly impact the primary β grain size as a result of epitaxial growth from pre-exiting substrate grains. In fact, our simulation results are in agreement with the widely accepted solidification theory for Ti64. Indeed, the partition coefficients of Al and V in Ti are very close to unity, meaning that almost no solute is rejected from the solid phase during solidification. In addition, the freezing range (i.e. the temperature difference between the liquidus and solidus curves) of Ti64 is extremely small (less than 1 K in our model). These factors altogether lead to a very small constitutional undercooling in front of the solid/liquid interface, thus preventing the development of interfacial instabilities such as cellular or dendritic structures. Moreover, the high thermal gradient in the liquid phase (G_L ∼ 2 × 10⁷ K m⁻¹ at the beginning of the simulation) prevents nucleation of new grains in the liquid phase. Therefore, these thermal conditions favor the development of large continuous epitaxial β grains exhibiting a planar interface without significant partitioning, as observed in our simulations.

Zoom In Zoom Out Reset image size

Due to the unique thermal history and non-trivial elemental redistribution during the L-PBF AM process, non-conventional solid-state phase transformation microstructures may develop in metallic systems. For instance, during L-PBF AM of Ti64, the local temperature in the printed layer often cyclically varies across the $\beta$ transus (i.e. the $\beta$-to-$\alpha$ transformation temperature) due to the thermal impact of successive layer melting. This strongly non-linear dynamic temperature profile leads to complex microstructures incorporating hierarchical and multiscale features ranging from nm to μm. Specifically, the primary $\beta$ solidification microstructure (e.g. grain size, morphology, and crystallographic texture) together with the $\alpha \left( {\alpha {^{^{\prime}}}} \right)$ second-phase precipitate characteristics (e.g. volume fraction, size, shape, orientation, coherency state, and spatial distribution) dictate the mechanical properties (e.g. strength, ductility, fatigue, and fracture toughness) of the alloy. Therefore, it is important to isolate the quantitative relationship between the temperature history and key multiscale microstructural features in order to relate processing parameters and resultant AM microstructures [14]. To this end, within our DT framework, we have formulated PF models for non-isothermal solid-state phase transformations in metallic alloys involving multiple solid-state allotropic phases for the AM processes. For Ti64 for example, the framework is based on a well-trained PF formalism for the precipitation of $\alpha$ phase, coupled with reliable CALPHAD database that provides temperature- and composition-sensitive thermodynamic and diffusion mobility parameters. The cooling and heating rates are informed by the location-dependent time-temperature profile from macro-scale process models for the given processing conditions. To account for physically relevant phase nucleation events, our PF framework incorporates the explicit nucleation algorithm (ENA) based on a probabilistic Poisson seeding process [51]. For the ENA, the nucleation rate for the given temperature and local composition is calculated by the classic nucleation theory, for which thermodynamic driving force, activation energy barrier for nucleation, and kinetics factors (e.g. Zeldovich factor and the rate at which atoms or molecules are attached to the critical nucleus) are informed by both well-assessed thermodynamic and kinetics databases. Therefore, concurrent nucleation, growth, coarsening, and dissolution of different precipitate phase during rapid cyclic heat treatment, with full-field microstructural information, can be treated cohesively in a single framework within the PF context. More details can be found in our previous work [27]. Figure 8 includes an example wherein the size scale, number density, and spatial distribution of $\alpha$ precipitates in a bicrystal $\beta$ matrix vary significantly with the cooling rate from $\beta$ transus temperature [27].

Zoom In Zoom Out Reset image size

We employed microelasticity theory [52] to compute the effective elastic moduli [27, 53] of digital representations of sampled AM microstructures in $246\,\Delta x \times 96\,\Delta y \times 24\,\Delta z$ grids (with the cell size of 2 μm). The mechanical equilibrium equation was numerically solved for AM microstructures under applied strains. We used the Fourier-spectral iterative-perturbation method [54, 55] to efficiently solve the equation with the position-dependent elastic moduli [27].

The computed effective elastic moduli provide useful information for analyzing overall elastic responses of AM 316L stainless steel microstructures. In fact, the computed effective moduli of the microstructures are given as full $C_{ij}^{{\text{eff}}}$ tensors, symmetry of which slightly deviates from cubic (i.e. orthotropic symmetry is expressed) even though the input reference modulus for a single crystal incorporates cubic symmetry. Therefore, it is not obvious to define a physical quantity that comprehensively characterizes the properties of the computed effective elastic modulus tensor. It should be noted that for the sake of simple comparison between different AM microstructures in this study, we project the computed effective modulus tensor into the hypothetical cubic symmetry by defining the characteristic modulus components ${\bar C_{11}}$ (= $(C_{11}^{{\text{eff}}} + C_{22}^{{\text{eff}}} + C_{33}^{{\text{eff}}})/3$), ${\bar C_{12}}$ (= $(C_{12}^{{\text{eff}}} + C_{13}^{{\text{eff}}} + C_{23}^{{\text{eff}}})/3$), and ${\bar C_{44}}$ (= $(C_{44}^{{\text{eff}}} + C_{55}^{{\text{eff}}} + C_{66}^{{\text{eff}}})/3$). These indicate the average values of similar types of components in the modulus tensor. In addition, to quantify the deviation from cubic symmetry of the computed effective modulus tensor, we introduce a new metric, the so-called 'orthotropicity' factor ${O_{\text{h}}}$ = $\sqrt {{{\left[ {\Delta \left( {{C_{11}}} \right)/{{\bar C}_{11}}} \right]}^2} + {{\left[ {\Delta \left( {{C_{12}}} \right)/{{\bar C}_{12}}} \right]}^2} + {{\left[ {\Delta \left( {{C_{44}}} \right)/{{\bar C}_{44}}} \right]}^2}}$, where $\Delta ({C_{ij}})$ represents a standard deviation of tensor components that are used for computing ${\bar C_{ij}}$ (e.g. ${O_{\text{h}}}$ = 0: perfect cubic symmetry; ${O_{\text{h}}}$> 0: orthotropic symmetry). Table 1 includes computed characteristic effective modulus components (${\bar C_{11}}$, ${\bar C_{12}}$, ${\bar C_{44}}$) and the orthotropicity factor (${O_{\text{h}}}$) for the substrate and AM microstructures for the two selected beam shapes in comparison with the input reference values for a single crystal [56]. All of the fundamental elastic modulus parameters can be derived from ${\bar C_{11}}$, ${\bar C_{12}}$, and ${\bar C_{44}}$, including the Zener anisotropy factor (${A_Z} = 2{\bar C_{44}}{ }/\left( {{{\bar C}_{11}} - {{\bar C}_{12}}} \right)$), Young's modulus (E), bulk modulus (K), shear modulus (μ), and Poisson's ratio (ν). We caution that these parameters may not fully represent all the features of the computed $C_{ij}^{{\text{eff}}}$ due to our intended definition of ${\bar C_{11}}$, ${\bar C_{12}}$, and ${\bar C_{44}}$. However, we claim that these may well capture difference in elastic responses of complex AM microstructures with different features.

Table 1. Computed effective elastic properties of 316L stainless steel microstructures produced by single-track scans.

Property	Before scan (substrate)	After scan (CG beam)	After scan (TE beam)	Reference value
${\bar C_{11}}$ (GPa)	281.47	282.31	277.24	247.50
${\bar C_{12}}$ (GPa)	123.32	122.85	122.21	140.30
${\bar C_{44}}$ (GPa)	97.00	96.23	97.12	123.74
Az	1.227	1.207 (LAZ: 1.185)	1.253 (LAZ: 1.295)	2.309
Oh	0.039	0.020	0.050	0.00
E (GPa)	230.31	229.88	228.69	243.02
K (GPa)	176.03	176.01	173.89	176.03
μ (GPa)	89.83	89.63	89.27	95.68
ν	0.282	0.282	0.281	0.270

As expected, the effective parameters deviate from the single crystal counterpart. However, no noticeable differences in the parameters were observed in the initial polycrystalline substrate or after scanning with both CG beam and TE beam, except for the Zener anisotropy and the orthotropicity factors. We note that the effective properties include both contributions from LAZ and non-LAZ regions in the simulated AM microstructures. In particular, for the Zener factor, we extracted the effective A_Z only for the LAZ ($A_{\text{Z}}^{{\text{LAZ}}}$) according to the formula: $A_{\text{Z}}^{{\text{LAZ}}} = \left[ {A_{\text{Z}}^{{\text{tot}}} - A_{\text{Z}}^{{\text{sub}}}\left( {1 - {f_{{\text{LAZ}}}}} \right)} \right]/{f_{{\text{LAZ}}}}$, where $A_{\text{Z}}^{{\text{tot}}}$ is the effective anisotropy factor for the total volume of the extracted microstructure sample, $A_{\text{Z}}^{{\text{sub}}}$ is the anisotropy factor for the substrate before the laser scan, and ${f_{{\text{LAZ}}}}$ is the volume fraction of LAZ, assuming that the anisotropy factor for the total volume is determined by a volume-averaged value. These calculations indicate that the microstructure for the TE beam exhibits stronger elastic anisotropy than the CG beam counterpart. In addition, the computed ${O_{\text{h}}}$ values indicate that the TE beam makes the microstructure have stronger orthotropic symmetry than the substrate before the laser scan, whereas the CG beam leads to more cubic-like symmetry. These arise from the microstructural differences for the two beam cases in terms of grain morphology and texture as analyzed above.

We note that full $C_{ij}^{{\text{eff}}}$ tensors of the effective moduli can serve as input parameters for the larger part-scale thermomechanical analysis during metal AM [57]. As we discussed, the effective property parameters of polycrystalline AM or substrate microstructures are noticeably different from the parameters for a single crystal. In particular, by employing the established AM processing–microstructure relationship (sections 2.2 and 2.3), the local elastic property in the part-scale model can be continuously updated by the microstructure-aware effective moduli over the course of the part-scale AM processing simulations, which may enable in silico and in situ monitoring of site-specific, thermally induced residual stress during the process.

We analyzed the spatial distributions of the local internal stress within the same sampled AM microstructures under applied loads. In particular, we computed microstructure-level von Mises stress (${\sigma _{{\text{vM}}}}$) profiles based on microelasticity theory [58, 59]. Figure 10 insets show the spatial distributions of ${\sigma _{{\text{vM}}}}$ in the microstructures after the laser scan by TE beam (figure 10(a)) and CG beam (figure 10(b)) as examples. It is worth noting that this micromechanical modeling naturally accounts for microscopic defects such as pores at the microstructure level, which cannot be captured in larger-scale simulations (e.g. effective medium approach in which material microstructure is homogenized). For instance, the simulated AM microstructure for the TE beam scan contains pores (e.g. the sampled microstructure in figure 9 incorporates 0.48% pores in the volume fraction) and its mechanical response reflects the presence of those defects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

More quantitatively, we compared the statistics of ${\sigma _{{\text{vM}}}}$ before and after the laser scan. To focus on the impacts of laser melting on the local stress profile, we extracted the von Mises stress (${\sigma _{{\text{vM}}}}$) statistics only from the LAZ. For comparison, we also extracted the statistics from the same region in the substrate. By directly comparing the stress profiles from the same region (before and after the laser scan), we can explicitly isolate the impact of the AM process. We considered applied strains ($\epsilon = 0.0023$) along all three orthogonal directions (see figure 9 for the definition of x-, y-, and z-directions). As clearly shown in the statistics (figure 10(a)), the two microstructures (i.e. before and after the scan) exhibit different stress profiles. For all three applied strain directions, the TE beam tends to move the ${\sigma _{{\text{vM}}}}$ distribution toward the lower stress regime. We speculate that this is caused by the presence of relatively smaller equiaxed grains [26]. For comparison, we also conducted the same analysis for the CG beam case. Interestingly, the CG beam results in an increased population of higher ${\sigma _{{\text{vM}}}}$ in the microstructure (see figure 10(b)).

To utilize the stress profiles for diagnosing the mechanical failure tendency, we define so-called 'stress hot spots' with a simple criterion: ${\sigma _{{\text{vM}}}} > {\sigma _{\text{Y}}}$, where ${\sigma _Y}$ is the yield strength (=500 MPa [60] was used for this study). In fact, the stress hot spots do not necessarily point to the locations where mechanical failure occurs. However, these may indicate the probable sites with a propensity of mechanical failure due to the formation of structural defects. Figure 11(a) shows an example of stress hot spots for microstructures produced by the TE beam, which are under applied strains ($\epsilon = 0.0023$) along the three orthogonal directions. We observed different hot spot densities for the three different loading directions, indicating mechanical response anisotropy. We further quantified the stress hot spots. In particular, we monitored volume fractions of stress hot spots with increasing applied loads. To assess the impact of the laser scan, we first compared the hot spot evolution in the substrate microstructure and the TE beam-induced microstructure for applied strains along the three orthogonal directions (see the results in figures 11(b) and (c), respectively). Note that we only tracked the microstructure regions defined by the LAZ in both cases. Due to the innate anisotropy in system size and shape, as well as grain texture therein, the substrate microstructure exhibits directionally anisotropic hot spot evolution (see figure 11(b)). Interestingly, after the TE beam scan (see figure 11(c)), the evolution behavior along the lateral x- and y-directions (i.e. laser scanning and transverse directions, respectively) is somehow regulated and becomes isotropic, deviating from the z-directional (i.e. build direction) behavior.

Zoom In Zoom Out Reset image size

For more comprehensive and comparative analysis of the different beam shaping effects, we present the hot spot evolutions in the microstructures before and after the CG and TE beam scans for the individual applied strain directions. These results are shown in figures 11(d)–(f). Our results indicate that the CG beam generally modifies the microstructure in a way that the hot spot evolution becomes more sensitive (i.e. steeper slope) to the applied strain (for all three directions) as compared to the TE beam. This tendency is more apparent for strain applied along the build direction (z-direction) as shown in figure 11(f). In fact, the TE beam decreases the sensitivity of the microstructure to strain applied along the z-direction than the substrate, indicating an improved mechanical response in comparison to the CG beam.

We argue here that the evolution of hot spot density with increasing strain may indicate the mechanical response of a specific site upon mechanical loading induced by rapid local thermal cycles (i.e. thermally induced residual stress) during AM. As mentioned for the static hot spot analysis above, having the local stress exceed the yielding strength (i.e. hot spots) may not directly lead to mechanical failure. In fact, upon rapid thermal cycling, the local microstructure may experience relatively large deformation without introducing fully resolved plasticity (due to time and mechanical constraints). However, the hot spots can be considered to be regions more vulnerable to mechanical attack. Hence, we hypothesize that hot spot evolution could indicate the microstructure-aware propensity for mechanical failure under the thermal cycling conditions during the AM processes. Therefore, minimizing the hot spot density and reducing its sensitivity to increasing load would be preferred. Note that even though this modeling neglects detailed behavior beyond yielding, this computational diagnosis is simpler and more efficient for considering complex AM microstructures compared to other sophisticated mechanical models that account for detailed plastic deformation modes. We also claim that this method provides information about early-stage deformation signatures affected by microstructural features, which enables early detection of site-specific mechanical failure. For better fidelity, accurately relating this model to mechanical failure modes in a large variety of sampled AM microstructures is additionally required.

To predict the mechanical response behavior beyond the yielding point, we also ran crystal plasticity simulations. This is far more sophisticated in its consideration of deformation modes and provides more direct information about deformation behavior beyond yielding, while more advanced numerical algorithms and extended computational resources are required. For capturing the elasto-viscoplastic response of a polycrystal, crystal-mechanics-based constitutive models are used in a standard FE formulation. The constitutive model formulation is similar to ones described in [61, 62], and the use of such models has been valuable for predicting mechanical response and for comparing these predictions to detailed experimental measurements [63–66]. Simulation throughput is facilitated by high-performance computing and the use of parallel solvers for the velocity field at each time step using the ALE3D code [67]. Figure 12 compares the computed elasto-plastic behavior of the AM 316L stainless steel microstructures for single tracks of CG and TE beams. We employed exactly the same microstructures used for the microelasticity modeling above.

Zoom In Zoom Out Reset image size

Two interesting features in simulated plastic behavior are observed in figure 12 in terms of comparison between the two beam shape cases. Note that we divided the simulated elasto-plastic behavior into three regimes for comparison. First, as shown in figure 12(a) (regime I), the microstructure with the TE beam is slightly more compliant (i.e. the slope is shallower) in the elastic regime (I), which is consistent with our microelasticity modeling of effective elastic moduli (see table 1). On the other hand, the CG beam case exhibits a slightly higher yield strength than the TE beam counterpart according to the predicted stress–strain curves in figure 12(a). Second, for the CG beam case, the dislocation density initially evolves more slowly with applied strain (see the early stage in regime II), whereas it becomes more rapid at later stages (from the later stage in regime II) than the TE case. Eventually, the microstructure produced by the CG beam allows more deformation modes to be activated in regime III. We speculate that the microstructure for the CG beam incorporates more heterogeneity in the deformation field that may accommodate more dislocation activity. Combining the simulated hot spot evolution (figure 11) and dislocation density evolution (figure 12), we claim that the microstructure for the TE beam single-track has a lower propensity for mechanical failure under identical applied loading conditions. However, to certify the accuracy and relevance of the predictions, more rigorous comparison with experimental characterization and testing are required. Moreover, more foundational investigation to connect this elasto-plastic deformation behavior to actual failure modes at the microstructure level is required.

Additionally, it is interesting to note that, in both cases, the overall microstructural setups and the corresponding applied loading conditions are configured in a way that the dislocation density does not make much contribution to the flow strength, resulting in fairly small amount of hardening for small amounts of strain in this study. Although we demonstrated the difference in elasto-plastic responses between CG and TE beam cases in figure 12, more quantitative understanding requires the microscopic analysis of the operating deformation modes and their activities/sensitivities to the applied loading conditions and grain/cell boundaries in the microstructure, which are currently underway.