A Novel Approach to Predict the Process-Induced Mechanical Behavior of Additively Manufactured Materials

Abstract

The grain structure and texture of additively manufactured materials depend strongly on the local temperature gradients during the solidification of the material. These grain structures and textures influence the mechanical behavior, ranging from isotropy to transversal and orthotropic symmetry. In the present contribution, a cellular automaton is used to model the grain growth during selective electron beam melting. The resulting grain structures and textures serve as input for a mesoscopic mechanical model. The mechanical behavior on the mesoscale is modeled by means of gradient-enhanced crystal plasticity, applying the finite element method. Computational homogenization is applied to determine the resulting macroscopic elastic and plastic properties of the additively manufactured metals. A general orthotropic yield criterion is identified by means of the initial yield loci computed with mesoscopic simulations of representative volume elements. The numerical results are partly validated with experimental data.

Introduction

Additive manufacturing (AM) of metals like selective laser melting or selective electron beam melting (SEBM) is widely used to manufacture prototypes and becomes more and more attractive for low volume production and spare parts (Ref 1). Freedom of design and no need of part-specific tools make AM very interesting for these applications. Furthermore, the mechanical properties of the material can be adjusted by tailoring its mesostructure using appropriate process strategies (Ref 2,3,4,5). Applying this approach to design parts with customized mechanical properties necessitates a modeling and simulation method which is able to predict the macroscopic mechanical behavior depending on the mesostructure and, thus, on the process strategy. A purely experimental approach to determine the relationship between the process and macroscopic properties would be very complex and expensive. To reduce the experimental effort, in this contribution a meso-mechanically motivated approach is presented to identify the macroscopic material behavior. A crystal growth simulation developed in previous works (Ref 6, 7) and a crystal plasticity-based mechanical simulation implemented in an in-house code using the finite element library deal.II (Ref 8) are linked to compute the process-induced mechanical behavior of additively manufactured materials.

Three major approaches have been established for microstructure prediction during solidification in recent years, namely phase field methods (PF), cellular automaton (CA) as well as Monte Carlo (MC)-based methods (Ref 9). While the numerical costs of PF-models restrict their application to the microscale (Ref 10) and MC- methods lack reliable physical coherencies, CA-models have prevailed for the prediction of the resulting grain structure. Ao et al. (Ref 11) combined the thermal field resulting from a moving laser beam with a two-dimensional cellular-automaton (CA) model to calculate the microstructure and thermal characteristics of the SLM process using AlSi10Mg. With this model, they investigated the effect of scanning speed, scanning spacing, and pre-heating temperature on the resulting microstructure. Zinovieva et al. (Ref 12) applied a finite-difference (FD) heat solver in combination with a three-dimensional CA model to analyze the grain structure and texture of SS316L resulting from SLM. They found a predominating Goss-texture as well as a weaker cube texture developing in their simulations. Li and Tan (Ref 13) investigated the influence of grain nucleation in direct laser deposition processes by combining a mesoscale CA model with a macroscale thermal model. An overview on modeling of microstructure–mechanical property relationships in AM of metals is given in Ref 14. In the context of additive manufacturing, crystal plasticity is applied, e.g., in Ref 15 where a local crystal plasticity formulation and numerical homogenization is used to model mechanical behavior of Inconel and results are compared with uniaxial tensile tests at different temperatures and in different directions. Similar applications applying local crystal plasticity formulations are given in Ref 16 or Ref 17 where the tool DAMASK is applied.

Here, a cellular automaton-based method is used to predict the resulting mesostructure of a part built additively from Inconel 718 by SEBM. The simulated crystal structure and texture serve as a basis to generate representative volume elements (RVEs) modeling the columnar grained mesostructure. The finite element method is applied to solve the balance of linear momentum and in addition a micro-force balance resulting from the gradient-enhancement of the applied single crystal plasticity model. The transition from meso to macroscale is achieved by computational homogenization. The gradient-enhanced crystal plasticity model was introduced in Ref 18 and is based on Ref 19,20,21. In general, local and non-local formulations of crystal plasticity can be distinguished. Local formulations as proposed in Ref 22, 23 consider only quantities at the present material point for the local constitutive behavior. Thereby size effects can only be considered with approaches like core–shell like strategies (Ref 24) or grain size dependent locally adjusted material parameters (Ref 25). In contrast, non-local formulations (Ref 19, 26,27,28,29) introduce a length-scale that enables to consider size effects directly but increase numerical cost. Furthermore, the grain boundaries can be modeled based on misorientations of the lattice structure (Ref 20, 30, 31). An approach to reduce the numerical cost is using an equivalent strain gradient model (Ref 32) instead of the usual slip gradient model. This approach reduces the number of additional scalar fields from 12 for a face centered cubic material to only one.

The macroscopic behavior of a material can be derived from simulating RVEs at the mesoscale and identifying a suitable macroscopic model and its parameters by fitting the macroscopic stress–strain responses. In Ref 33, a theory is developed to determine the yield locus of poly-crystalline materials from an analytical point of view based on Schmid’s law. Examples in two-dimensional stress space with $\langle 100\rangle$ textured materials, similar to columnar grained structures in additive manufacturing, are given in Ref 34. As a result, the fundamental form of an fcc polycrystals yield locus in the deviatoric plane is a four-sided-shape. According to Ref 35, the yield locus of a single surface yield criterion lies in between a hexagon inscribed into and another one circumscribing the von-Mises yield circle in the deviatoric plane. The common anisotropic yield criteria scale various axes in the stress space differently. Thus, it is not possible to match the theoretically derived four-sided shape of the yield locus of an fcc material by means of an anisotropic single surface yield criterion. Another approach is choosing a multi-surface-based anisotropic yield criterion. A detailed discussion on the modeling of fcc metals with a multi-surface-based anisotropic $J_2$ yield criterion is given in Ref 36. In Ref 37, 38, the same principle is used to model the plastic deformation of paperboard. A different approach to define an anisotropic yield criterion is the concept of Isotropic Plasticity Equivalent material as for example given in Ref 39. The idea is to separate the anisotropy from the yield function and to use a fourth-order projection tensor to transform the stresses from an anisotropic stress space into an isotropic one. Since the homogenized results from RVE simulations reproduce qualitatively the analytical solution from Ref 33, a multi-surface-based anisotropic yield criterion is favored here.

In Sect. 2, a cellular automaton-based method to model crystal growth during selective beam melting processes is described. Section 3 summarizes the applied gradient enhanced single-crystal plasticity model. In Sect. 4, the numerical homogenization is summarized, and the transition of results to the macroscale and the anisotropic yield criterion are given. Numerical results are given in Sect. 5. First, the results from grain growth simulations are presented and compared with experimental data. Subsequently, a representative volume element is derived and the macroscopic mechanical behavior is investigated. Finally, results using a Voronoi tessellation-based generated RVE are compared with the cellular automaton-based one. In Sect. 6, the results are summarized and discussed.

Grain Growth Model

The crystal growth model utilized in this work for predicting the microstructure of a part built additively from IN718 powder is based on a cellular automaton (CA) model originally developed by Gandin and Rappaz (Ref 40). With both the genuine model as well as its application to additive manufacturing being well documented in the literature (Ref 6, 7, 41), the following section concentrates on briefly summarizing the key aspects. In Sect. 5, both the experimental as well as the numerical parameters used within this work are presented followed by the resulting microstructure predicted by the software.

CA are mathematical procedures describing the discrete evolution (both in space and time) of complex systems by applying local transformation rules to the cells of a lattice (Ref 42). Gandin and Rappaz (Ref 43) utilized the method for modeling the dendritic growth of face-centered-cubic (fcc) metals during solidification. Based on experimental observations of fcc-crystals forming an octahedral shape when growing in the absence of a thermal gradient, the model approximates the envelope of a dendritic grain by an octahedral object. The vertices of the octahedron are representing the tips of the dendrites primary arms growing along the preferred growth directions (${<}100{>}$ in the case of fcc-crystals), mimicking its orientation in three-dimensional space.

Each cell in the CA lattice being part of a growing grain is assigned an octahedron, rotated according to the orientation of the dendrites. Solidification is modeled by step-wise growing of the octahedra with a velocity v depending on the undercooling $\Delta T$ prevailing at the cell it is assigned to. Consequently, the size of an octahedron $L^t$ at time t is calculated by

$$\begin{aligned} L^t = L^{t_c} + \frac{1}{3} \int _{t_c}^{t}v( {\Delta }T({\tau })) \, \mathrm {d}\tau \end{aligned}$$

(1)

with $\hbox {L}^{{\mathrm{t}}_{\mathrm{c}}}$ representing the size of the octahedron at the time $t_c$. The undercooling $\Delta T$ at the specific cells is calculated either analytically (Ref 6) or by numerical means (Ref 7). Based on the model of Kurz, Giovanola and Trivedi (KGT) (Ref 44), the relationship between v and $\Delta T$ is approximated by a quadratic expression (Ref 45)

$$\begin{aligned} v = A \, {\Delta }T^2 \end{aligned}$$

(2)

As soon as a growing octahedron comprises the center of an adjacent cell not yet part of a grain itself, this cell is captured. A new octahedron is generated, truncated according to specific rules (explained in detail by Gandin et al. (Ref 43)) and starts to grow on its own with a velocity determined by the undercooling prevailing at the center of the newly incorporated cell. This procedure is iterated until all fluid cells are captured and solidification stops.

Mesoscopic Mechanical Model

The mechanical behavior of the grain structure which is simulated by means of the cellular automaton is modeled with the gradient-enhanced crystal plasticity model introduced in Ref 18. In the following section, the governing equations of gradient-enhanced single crystal plasticity are shortly given. A detailed discussion can be found in Ref 18 and references therein.

Kinematics

The strain ${\varvec{\varepsilon }}$ is defined as the symmetric gradient of the displacement ${\varvec{u}}$ at a material point ${\varvec{x}}$ in the body $\mathcal {B} \in \mathcal {R}^3$ and additively decomposed into elastic and plastic parts

$$\begin{aligned} {\varvec{\varepsilon }} = \left[ \nabla {{\varvec{u}}} \right] ^ {sym} = {\varvec{\varepsilon }}^e + {\varvec{\varepsilon }}^p. \end{aligned}$$

(3)

The plastic strain is composed of slips $\gamma ^{\alpha }$ in predefined slip systems $\alpha$ in the so-called closest packed planes with normals ${\varvec{m}}^{\alpha }$ and slip directions ${\varvec{s}}^{\alpha }$. By summing over all slip systems, the plastic strain reads

$$\begin{aligned} {\varvec{\varepsilon }}^p = \sum _\alpha \gamma ^{\alpha } \left[ {\varvec{s}}^{\alpha } \otimes {\varvec{m}}^{\alpha } \right] ^{sym} = \sum _\alpha \gamma ^{\alpha } {\varvec{p}}^{\alpha }. \end{aligned}$$

(4)

The gradient enhancement of the model is motivated by geometrically necessary dislocations (GNDs). According to Ref 46, GNDs are the non-integrability of plastic distortions and thus the entries of the so-called Nye tensor. Due to the definition of plastic strain and the curl operation, only gradients tangential to the slip planes contribute to the dislocation densities

$$\begin{aligned} \rho _\perp ^\alpha = {\varvec{s}}^{\alpha } \cdot \nabla {\gamma }^\alpha , \quad \rho _\odot ^\alpha = {\varvec{l}}^{\alpha } \cdot \nabla {\gamma }^\alpha . \end{aligned}$$

(5)

${\varvec{l}}^{\alpha }$ is the second tangential basis vector in a slip system obtained by ${\varvec{l}}^{\alpha } = {\varvec{m}}^{\alpha } \times {\varvec{s}}^{\alpha }$. Geometrically, $\rho _\perp ^\alpha$ represent edge dislocation densities, the off diagonal entries of Nye’s tensor, and $\rho _\odot ^\alpha$ stand for screw dislocation densities, the diagonal entries of Nye’s tensor. For a more detailed discussion, the reader is referred to Ref 47 for example.

Balance Equations

Based on the principle of virtual power, equations for the displacement field and the crystallographic slips are derived. In the following, the balance equations are given directly, for a short recapitulation see Ref 18, whereas a detailed derivation is given in Ref 19. The energetically conjugated pairs are elastic strains and stresses ${\varvec{\varepsilon }}^e$—${\varvec{\sigma }}$, as well as crystallographic slips and corresponding gradient extended Schmid stresses $\gamma ^{\alpha }$—$\pi ^\alpha$. The energetically conjugated quantity to the tangential part of the microscopic slip gradients $\nabla {\gamma }^\alpha$ is a microstress ${\varvec{\xi }}^{\alpha }$. Both balance equations are coupled by the so-called Schmid stresses, shear stresses that induce crystallographic slips. The local Schmid stresses are obtained by a projection of the elastic stresses to the slip systems

$$\begin{aligned} \tau ^\alpha = {\varvec{\sigma }} : {\varvec{p}}^{\alpha }. \end{aligned}$$

(6)

The standard balance of linear momentum and the non-standard microscopic balance equation for each slip system $\alpha$ read

$$\begin{aligned}&\mathrm {div}\,{{\varvec{\sigma }}} + {\varvec{b}} = {\varvec{0}}, \end{aligned}$$

(7)

$$\begin{aligned}&\mathrm {div}\,{{\varvec{\xi }}^{\alpha }} + \tau ^\alpha - \pi ^\alpha = 0. \end{aligned}$$

(8)

Note that the elastic stresses ${\varvec{\sigma }}$ and the microstresses ${\varvec{\xi }}^{\alpha }$ are of energetic nature, whereas the gradient extended Schmid stress $\pi ^\alpha$ is of dissipative nature. The micro-force balance (8) results in a spatial regularization of the crystallographic slips and as a consequence allows for the definition of grain boundary conditions.

Constitutive Equations: Bulk

The balance equations are completed by constitutive equations that couple the kinematic quantities with their conjugated counterparts. A free energy $\Psi \left( {\varvec{\varepsilon }}^e,{\varvec{\rho }} \right)$ is defined depending on the elastic strains and all edge and screw dislocation densities $\rho _\perp ^\alpha$ and $\rho _\odot ^\alpha$, respectively, collected in the vector ${\varvec{\rho }}$ contributing to a defect energy $\Psi ^d \left( {\varvec{\rho }} \right)$

$$\begin{aligned} \Psi \left( {\varvec{\varepsilon }}^e,{\varvec{\rho }} \right) = \frac{1}{2}{\varvec{\varepsilon }}^e : {\varvec{\mathcal {C}}}^e : {\varvec{\varepsilon }}^e + \underbrace{\frac{1}{2}\pi _0 l_e ^2 \sum \limits _\alpha \left[ [\rho _\perp ^\alpha ]^2+[\rho _\odot ^\alpha ]^2\right] }_{\Psi ^d \left( {\varvec{\rho }} \right) }. \end{aligned}$$

(9)

${\varvec{\mathcal {C}}}^e$ denotes the fourth-order elasticity tensor with the symmetry class of the actual unit cells, cubic symmetry. The critical local Schmid stress is given by $\pi _0$, and $l_e$ is an energetic length, the only parameter that contributes to the model due the gradient enhancement. The elastic stress is obtained by the derivative of the free energy with respect to the elastic strain and reads

$$\begin{aligned} {\varvec{\sigma }} = \frac{\partial \Psi }{\partial {\varvec{\varepsilon }}^e}={\varvec{\mathcal {C}}}^e : {\varvec{\varepsilon }}^e. \end{aligned}$$

(10)

The derivative of the free energy with respect to the slip gradients leads to the micro-stresses

$$\begin{aligned} {\varvec{\xi }}^{\alpha } = \frac{\partial \Psi ^d }{\partial \nabla {\gamma }^\alpha } = \pi _0 l_e ^2 \left[ - {\varvec{s}}^{\alpha } \rho _\perp ^\alpha + {\varvec{l}}^{\alpha } \rho _\odot ^\alpha \right] . \end{aligned}$$

(11)

In order to overcome the problem of searching for active slip systems, the slips are regularized in time, applying a creep law. The often used powerlaw is numerically unstable close to the quasistatic limit. To deal with this issue, a hyperbolic tangent as proposed in Ref 21 is applied and the gradient extended Schmid stress reads

$$\begin{aligned} \pi ^\alpha = \left[ \pi _0 + g^\alpha \right] \mathrm {tanh}\left( \frac{\dot{\gamma }^\alpha }{\dot{\gamma }_0} \right) . \end{aligned}$$

(12)

To improve the stability of the numerical algorithm, the regularization parameter $\dot{\gamma }_0$ is adaptively adjusted and a point-wise damping of the solution scheme is applied. Self- and latent hardening are taken into account by a simple linear hardening relation with the hardening slope $h_0$ and a parameter $q > 0$ describing the ratio of self and latent hardening. The hardening stress $g^\alpha$ evolves as

$$\begin{aligned} \dot{g}^\alpha = \sum _\beta \underbrace{h_0 \left[ q + \left[ 1 - q \right] \delta ^{\alpha \beta } \right] }_{h^{\alpha \beta }} \dot{\gamma }^\beta \, \quad \text {with} \quad g^\alpha \left( t_0\right) = 0 . \end{aligned}$$

(13)

Grain Boundaries

The two standard boundary conditions at the grain boundaries are either of Neumann type, namely micro-free ${\varvec{\xi }}^{\alpha } \cdot {\varvec{n}} = 0$, or micro-hard $\gamma ^{\alpha } = 0$, a Dirichlet type boundary condition. In Ref 18, it is shown that micro-hard grain boundaries lead to overestimated macroscopic stiffness. A third option is to define a free energy between two grains to derive a constitutive equation for the micro-tractions between two grains, depending on their plastic slips. The equations based on the work in Ref 20 are shortly summarized. The free energy is based on the Burgers tensor on a grain boundary and is defined by

$$\begin{aligned} \Psi ^{gb} \left( {\varvec{\gamma }} \right)=\, & {} {\frac{1}{2} \lambda ^{gb} \biggl | \left[ {\varvec{\varepsilon }}^p _B - {\varvec{\varepsilon }}^p _A \right] \times {\varvec{n}}^{gb}_{A} \biggr |^2} \nonumber \\=\, & {} \frac{1}{2} \lambda ^{gb} \biggl | \underbrace{\left[ \left( {\varvec{\varepsilon }}^p _B \right) _{ij} - \left( {\varvec{\varepsilon }}^p _{A} \right) _{ij} \right] e_{jkl} \left( {\varvec{n}}^{gb}_{A} \right) _{k} \,{\varvec{e}}_i \otimes {\varvec{e}}_l}_{{\varvec{\mathcal {G}}}} \biggr |^2 \end{aligned}$$

(14)

with all slips $\gamma ^{\alpha }$ arranged in a vector ${\varvec{\gamma }}$, the third-order Levi–Civita tensor represented by $e_{ijk}$ and $\lambda ^{gb}$ a grain boundary stiffness. The subscripts A and B represent the grains at the grain boundary. By deriving the energy with respect to the slips, a respective constitutive relation for the energetically conjugated micro-traction $\chi ^\alpha _A$ on the grain boundary is obtained that reads

$$\begin{aligned} \chi ^\alpha _A={\varvec{\xi }}^{\alpha } _A \cdot {\varvec{n}}^{gb}_{A} = \frac{\partial \Psi ^{gb} \left( {\varvec{\gamma }} \right) }{\partial \gamma ^{\alpha } _A} = - \lambda ^{gb} \sum _\beta \left[ \gamma ^{\beta } _A C^{\alpha \beta }_{AA} + \gamma ^{\beta } _B C^{\alpha \beta }_{AB} \right] . \end{aligned}$$

(15)

Mind that $\lambda ^{gb}$ is the only parameter in the formulation. The energy will be used as a penalty formulation that ensures compatible plastic deformations at the grain boundaries. The geometric quantities $C^{\alpha \beta }_{XY}$ depend only on geometrical quantities as well as on the grain orientations and read

$$\begin{aligned} C^{\alpha \beta }_{XY} = \left[ {\varvec{s}}^{\alpha } _X \cdot {\varvec{s}}^{\beta } _Y \right] \left[ {\varvec{m}}^{\alpha } _X \times {\varvec{n}}^{gb}_{A} \right] \cdot \left[ {\varvec{m}}^{\beta } _Y \times {\varvec{n}}^{gb}_{B} \right] . \end{aligned}$$

(16)

Identification of Macroscopic Mechanical Behavior

Computational Homogenization

Computational homogenization is used to transfer the results from the mesoscale—the scale of particular grains—to the scale of parts—the macroscale. The following subsection presents the basics of computational homogenization for a detailed discussion see, e.g., Ref 48. The first-order-homogenization ansatz relates meso- and macroscopic quantities with

$$\begin{aligned} {\varvec{u}} = \bar{{\varvec{\varepsilon }}} \cdot {\varvec{x}} + \tilde{{\varvec{u}}};\quad {\varvec{\varepsilon }} = \bar{{\varvec{\varepsilon }}} +\tilde{{\varvec{\varepsilon }}}, \end{aligned}$$

(17)

with $\left( \bar{\bullet } \right)$ denoting macroscopic quantities and $\left( \tilde{\bullet } \right)$ denoting mesoscopic fluctuations. The macroscopic stresses $\bar{{\varvec{\sigma }}}$, strains $\bar{{\varvec{\varepsilon }}}$ and dissipation $\bar{\mathcal {D}}$ are obtained by averaging their microscopic counterparts

$$\begin{aligned} \bar{{\varvec{\varepsilon }}} = \displaystyle \frac{1}{\mathcal {V}}\int \limits _\mathcal {B} {\varvec{\varepsilon }} \,\mathrm {d}V;\quad \bar{{\varvec{\sigma }}} = \displaystyle \frac{1}{\mathcal {V}}\int \limits _\mathcal {B} {\varvec{\sigma }} \,\mathrm {d}V;\quad \bar{\mathcal {D}} = \displaystyle \frac{1}{\mathcal {V}}\int \limits _\mathcal {B} \mathcal {D} \,\mathrm {d}V \end{aligned}$$

(18)

$\mathcal {V}$ denotes the volume of the RVE. According to Ref 49, the meso- and macroscopic stress power have to be equivalent. This results in an extended Hill–Mandel condition giving a requirement for the fluctuations of the displacements and additionally for the crystallographic slips

$$\begin{aligned} \frac{1}{\mathcal {V}}\int \limits _{\partial \mathcal {B}} {\varvec{t}}\cdot \dot{\tilde{{\varvec{u}}}}\,\mathrm {d}A\ + \ \frac{1}{\mathcal {V}} \sum _\alpha \int \limits _{\partial \mathcal {B}} \chi ^\alpha \dot{\gamma }^\alpha \mathrm {d}A = 0. \end{aligned}$$

(19)

(19) is satisfied by a proper choice of boundary conditions such as periodic–antiperiodic displacement fluctuation - traction and slip - microstress boundary conditions

$$\begin{aligned} \begin{array}{lrclrcl} &{}\tilde{{\varvec{u}}} ^+ &{}=&{}\tilde{{\varvec{u}}} ^-\,;&{}\quad {\varvec{t}}^+&{}=&{}-{\varvec{t}}^-\\ \text {and} \quad &{}\left[ \gamma ^{\alpha } \right] ^+&{} = &{}\left[ \gamma ^{\alpha } \right] ^-\,;\quad &{} \left[ \chi ^\alpha \right] ^+&{}=&{} - \left[ \chi ^\alpha \right] ^- \quad \text {on} \quad \partial \mathcal {B}=\partial \mathcal {B}^+\cup \partial \mathcal {B}^- . \end{array} \end{aligned}$$

(20)

Computation of Macroscopic Properties

The macroscopic material properties are determined from the homogenization procedure using the homogenized stress $\bar{{\varvec{\sigma }}}$, the macroscopic strains $\bar{{\varvec{\varepsilon }}}$ and the macroscopic dissipation $\bar{\mathcal {D}}$. The macroscopic stiffness tensor ${\bar{ \mathcal {C} }}^e$ is obtained by means of the simulation of 6 unit loadcases with fixed slip fields $\gamma ^{\alpha } = 0 \quad \forall \quad {\varvec{x}} \in \mathcal {B}$ applying the macroscopic elasticity relation

$$\begin{aligned} \bar{{\varvec{\sigma }}} ={\bar{ \mathcal {C} }}^e : \bar{{\varvec{\varepsilon }}}^e. \end{aligned}$$

(21)

In order to identify the initial point of yielding, a dissipation-based criterion similar to the one given in Ref 50 is defined. The macroscopic dissipation $\mathcal {D}^\star$ for a simulated uniaxial tensile test with a plastic strain $\bar{\varvec{\varepsilon }} ^p_y$ of 0.2%

$$\begin{aligned} \bar{{\varvec{\varepsilon }}}^p = \bar{{\varvec{\varepsilon }}} - \left[ {\bar{ \mathcal {C} }}^e \right] ^{-1} : \bar{{\varvec{\sigma }}} = \bar{\varvec{\varepsilon }} ^p_y \end{aligned}$$

(22)

is identified and defined as criterion for the initial point of yielding. For other loadpaths, the macroscopic stress at initial yielding is evaluated at time $t^\star$, at which the identical macroscopic dissipation $\mathcal {D}^\star$ is obtained, by

$$\begin{aligned} \bar{\mathcal {D}} = \int \limits _0^{t^\star } \int \limits _{\mathcal {B}} \sum _\alpha \dot{\gamma }^\alpha \pi ^\alpha \, \mathrm {d}V \, \mathrm {d}t {\mathop {\,=\,}\limits ^{!}} \mathcal {D}^\star . \end{aligned}$$

(23)

Macroscopic Yield Function

Due to its ability to reproduce the characteristic biaxial behavior of fcc polycrystalline structures, a multi-surface-based anisotropic yield criterion is used to model the homogenized macroscopic yielding of the RVE. The particular yield function is briefly introduced here, for a detailed discussion see Ref 36. Other applications are given in Ref 37, 38. According to Ref 36, three yield surfaces as a basis are accurate enough to model the orthotropic plastic behavior of polycrystalline aggregates in engineering applications. The criterion is based on three surfaces each only considering the $J_2$-type invariants of the stress deviator. The first invariant is neglected because volumetric stresses do not cause yielding in metals what is in accordance with Schmid’s law given in (4) and (6). The consideration of the third invariant would lead to different behavior in tension and compression which is not observed for metals and in accordance with the model on the mesoscale where the hyperbolic tangent in (12) does not distinguish between positive and negative values. A yield function can be written like

$$\begin{aligned} f \left( J_2^{'(1)}, J_2^{'(2)}, J_2^{'(3)} \right) = 0. \end{aligned}$$

(24)

The relations between the stresses in the three different yield systems are defined similarly as in Ref 36 and result in the following definition of the three $J_2$-type invariants

$$\begin{aligned} \begin{bmatrix} J_2^{'(1)}\\ J_2^{'(2)}\\ J_2^{'(3)} \end{bmatrix} = \frac{1}{6} \underbrace{ \begin{bmatrix} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}\\ a_{11}^{-2}&{} a_{22}^{-2}&{} a_{33}^{-2}&{} a_{12}^{-2}&{} a_{13}^{-2}&{} a_{23}^{-2}\\ b_{11}^{-2}&{} b_{22}^{-2}&{} b_{33}^{-2}&{} b_{12}^{-2}&{} b_{13}^{-2}&{} b_{23}^{-2}\\ \end{bmatrix}}_{\mathbf {M}} \begin{bmatrix} \left[ \sigma _{11} - \sigma _{22} \right] ^2\\ \left[ \sigma _{22} - \sigma _{33} \right] ^2\\ \left[ \sigma _{33} - \sigma _{11} \right] ^2\\ 6\sigma _{12}^2\\ 6\sigma _{13}^2\\ 6\sigma _{23}^2\\ \end{bmatrix}. \end{aligned}$$

(25)

The yield criterion is then formulated as a quartic polynomial of the $J_2$-type invariants

$$\begin{aligned} f \left( {\varvec{\sigma }} \right) = \left( J_2^{'(1)}\right) ^4 + \left( J_2^{'(2)}\right) ^4 + \left( J_2^{'(3)}\right) ^4 - 1. \end{aligned}$$

(26)

This is different to Ref 36 where a quadratic polynomial was used. However, Eq 26 is formulated in accordance with Ref 51 which states that a yield function with exponents of the stresses of 8–10 gives the best results for FCC materials. The total number of parameters to identify is 13, $a_{ij}$, $b_{ij}$ and c. It is assumed here that the mechanical behavior of the material is invariant for rotations of $90^\circ$ about the building direction. Thus, the yield function (26) needs to be invariant, e.g., for substituting $\sigma _{11}$ and $\sigma _{22}$. Further sets of interchangeable stress components are given in Table 1.

Table 1 Symmetry assumptions for the material behavior. Interchangeable pairs and the resulting constraints are given

By means of comparison of coefficients, one obtains a reduced set of 7 parameters and the matrix M which defines the three $J_2$-type invariants reads

$$\begin{aligned} \mathbf {M}= \begin{bmatrix} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}&{} c^{-2}\\ a_{11}^{-2}&{} a_{22}^{-2}&{} a_{33}^{-2}&{} a_{12}^{-2}&{} a_{13}^{-2}&{} a_{23}^{-2}\\ a_{11}^{-2}&{} a_{33}^{-2}&{} a_{22}^{-2}&{} a_{12}^{-2}&{} a_{23}^{-2}&{} a_{13}^{-2}\\ \end{bmatrix}. \end{aligned}$$

(27)

Numerical Results

In this section, the homogenized mechanical behavior of selective electron beam molten Inconel 718 is investigated. First, results from the grain growth simulation are presented and compared with experimental measurements. Based on the grain growth results, representative volume elements are constructed and material parameters are inversely identified using the homogenized mechanical behavior of the RVE and experimental data from macroscopic uniaxial tensile tests. Following this, the macroscopic elastic tangent and the initial yield locus are identified from uni- and biaxial loading conditions on the RVE, and the parameters of the macroscopic yield criterion are identified. Finally, the results for CA simulated and Voronoi-tessellation-based generated RVEs are compared.

Grain Growth for Inconel 718

In a previous work (Ref 6), the capabilities of the crystal growth model outlined in Sect. 2 were demonstrated by calculation of the resulting microstructure of a part built additively via SEBM from IN718 powder. A cross-snake hatching strategy was applied with the beam propagation direction rotating by $90^{\circ }$ every new layer. The predicted microstructure was in excellent accordance with experimentally derived EBSD measurements. For deriving the mechanical properties, a similar setup was used in this work, albeit with a higher spatial resolution. The numerical domain was sufficiently dimensioned to encase the RVE used for deriving the mechanical properties. The experimental setup reproduced in the simulation as well the numerical parameters used for modeling are listed in Table 2. To avoid distorting boundary effects, only the core of the resulting microstructure consisting of $160\times 160\times 160$ cells was employed as RVE for the subsequent mechanical characterization.

Table 2 Parameters used for the grain growth model

The resulting grain structure is depicted in Fig. 1. To reduce the computational load, the microstructure is calculated within the re-molten zone on top of the specimen solely. This approach enables the simulation domain to remain in a constant size mimicking the melt pool depth regardless of the final build height. After completion of the simulation, the resulting microstructures of the individual layers are combined, forming a consistent data set.

Fig. 1

(a) Snapshots of developing microstructure during the melting of the first four layers, (b) resulting three-dimensional microstructure after the build, (c) two-dimensional cut through the simulated domain, (d) experimentally derived EBSD-measurement (adapted from Ref 6) under the term of CC BY-NC-ND license

Figure 1(a) shows the simulation domain during the melting of the first four layers. Initially, the developing persistent melt pool traverses the domain from back to front. Due to the rotation of the beam direction by $90^{\circ }$ every layer, the melt pool traverses from right to left in the second layer, from front to back during the third and finally from left to right in the fourth layer. Longitudinal cuts through the assembled microstructure of the first three millimeters are shown in Fig. 1(b) in 3D as well as 2D (c). Figure 1(d) shows an EBSD measurement taken from the same section using a FEI Helios 600i dual-beam SEM with a step size of $5\,\upmu {\hbox {m}}$. The speckled clusters of new grains visible in the EBSD-measurements originate from impurities within the powder or interferences during the powder deposition. Both phenomena are not considered in the model up to now.

A detailed analysis of both the average grain size as well as the prevailing grain orientation is shown in Fig. 2. Starting from the initial grain size of $30\,\upmu {\hbox {m}}$ the grains are broadening linear with the build height for the first 5 mm (red dotted line). After reaching an average grain diameter of $65\,\upmu {\hbox {m}}$, the coarsening ceases, resulting in a constant grain size from 5 mm upward.

Fig. 2

Grain size progression derived by the CA simulation accompanied by pole figures representing the distinct texture at build heights indicated by the dotted lines

The pole figures plotted in Fig. 2 clearly visualize the developing texture in the specimen. From the initially randomly distributed grains of the start configuration soon grains emerge with their (100)—orientation best aligned in build direction, resulting in a fiber texture. From these, the grains with their secondary dendrite arms aligned with the direction of the persistent melt pool prevail. From this time on, the grain size remains constant for the remainder of the build.

Representative Volume Elements

The considered RVEs represent the columnar grained Inconel 718 depicted in Fig. 1. An inner cutout of the domain at the top surface when grain size and texture are converged is chosen. To obtain a periodic geometry from the grain growth simulations, small modifications are made at the boundaries of the RVE: grains are forwarded periodically and grain boundaries are smoothed. Care is taken not to generate too large non-representative grains. In Fig. 3, the result of a grain growth simulation (a), its simplification (b) and the resulting RVE for the mesoscopic mechanical simulations (c) are shown. To represent the columnar grained structure the cross section is extruded in building direction mimicking infinitely long grains approximating the high aspect ratio from experimental results. For an easier interpretation of the results, the axes of the RVE are aligned with the secondary anisotropy of the grown crystals, depicted in the polfigures in Fig. 2. In Ref 18, it could be shown that the exact shape of the grains has only minor influence on the homogenized mechanical behavior of columnar grained material. In contrast, it is very sensitive to changes in the texture. As a result, the geometry of the RVE is kept constant, but the grain orientations will be varied to investigate the statistical behavior. Different numerically grown crystal structures are used to obtain sets of unit cell orientations for the generation of different RVEs. Periodic boundary conditions are applied on opposite sides for the fluctuations of the displacement and the slips as well.

Fig. 3

(a) Cross section of columnar grains from CA simulations, (b) CA result and RVE grain boundaries, (c) CA-based (simulated) RVE including mesh. Grains are colored to distinguish between different grains

Inverse Parameter Identification

The parameters identified in Ref 18 by means of inverse parameter identification comparing the homogenized results from RVE simulations with experimental data from uniaxial tensile tests are used. The experimental results of the uniaxial tensile tests are given in Ref 52. Specimens with three different orientations in the building chamber—parallel, perpendicular and under $45^\circ$ to BD—were produced and tested. The orientations correspond to the softest, an intermediate and the stiffest directions of the material. The parameters are given in Table 3. Thereby, $\nu$, $q$, $\lambda ^{gb}$, $l_e$ and $\dot{\gamma }_0$ were fixed to reasonable values and $E$, $\mu$ $\pi _0$ and $h_0$ were inversely identified by means of the experimental data. In Fig. 4, the numerical results are compared with the experimental tensile tests under loading perpendicular (left) and parallel (right) to the building direction.

Table 3 Parameters of the model

Fig. 4

Stress strain curves due to uniaxial tension parallel and perpendicular to the BD with the generated 12 crystal RVE in comparison with experimental data from Ref 52

Macroscopic Material Behavior

In contrast to experimental investigations of the material, the numerical model is not restricted due to the complexity of experimental set-ups or, e.g., buckling in compression, but all possible load cases can be simulated. This advantage is exploited in the following to identify the macroscopic elastic and plastic properties of the additively manufactured material by means of applying various load cases to the RVEs.

Elastic Material Behavior

First, the elastic tangent is evaluated as described in Sect. 4.2 using 100 different RVEs with constant geometry but varying sets of crystal orientations from the corresponding polfigures. Figure 5 shows the Young’s moduli of the RVEs in all directions as the radius of the figure and additionally coded in color. It can clearly be seen that the Young’s modulus of the columnar grained material is mainly transversal isotropic but has additionally a slight secondary anisotropy. This is due to the cross snake scan strategy with rotating beam directions every 10 layers and thus the growth of grains and selection of orientations with a second preferred direction, compare the polfigure in Fig. 2. As a result, the BD corresponds to the softest and the scan vector directions to the second softest orientations. The stiffest direction is under $45^\circ$ against BD and scan vector axis.

Fig. 5

Three-dimensional representation of the orientation dependent Young’s modulus

This behavior can be utilized to maximize the stiffness of a part. Since the preferred direction is the building direction of the process, changing the orientation of parts in the building chamber changes the stiffest direction in parts. In Ref 53, a topology optimization of additively manufactured parts was performed, accounting for this material anisotropy as an additional degree of freedom. An increase of up to 30% in stiffness could be achieved depending on the loadcase.

Macroscopic Yield Locus

As expected and visualized in Fig. 4, additionally to the Young’s modulus (slope in the elastic region) investigated in the previous section, the plastic behavior like the yield stress also depends on the loading direction and type. A sample uniaxially loaded parallel to the BD is less stiff but has a higher yield stress than loaded perpendicular to BD. Applying the ansatz described in Sect. 4.2, the yield locus for the generated RVE in the $\sigma _{11} / \sigma _{22}$ biaxial stress sub-space is evaluated for ten realizations of the RVE. In Fig. 6, the results for all ten realizations (left) as well as their mean values (red line) and ± standard deviation (grey) are given (right). It can be seen that ten realizations of the RVE with 12 crystals are sufficient.

Fig. 6

Points of yielding in biaxial loading in the $\sigma _{11} / \sigma _{22}$ stress sub-space for 10 different realizations of simulated RVEs

The initial yield loci for all biaxial combinations of principal stresses are depicted in Fig. 7. All combinations of biaxial stresses are calculated and compared with the symmetry assumptions given in Sect. 4.3, proofing that the assumptions are correct. Additionally the yield function (26) could be fitted to the results and shows a good agreement with the RVEs behavior. The characteristic rhomboidal shape is obtained. The parameters of the determined yield function are given in Table 4.

Table 4 Parameter of the determined yield function

Fig. 7

Yield locus for biaxial principal stress states determined for a CA-based RVE (red curves). The initial yield locus of the macroscopic anisotropic yield criterion is given for comparison (blue curves). The green dashed line in the $\bar{\sigma }_{22}-\bar{\sigma }_{33}$ diagram is a copy of the red curve in the $\bar{\sigma }_{11}-\bar{\sigma }_{33}$ diagram and is inserted to verify the symmetry assumptions from Table 1. Indeed both lines are exactly spot on (Color figure online)

As discussed, it is rather complicated to test the other biaxial combinations of stresses in experiments. Thus, from an experimental point of view the yield function (26) matches the material’s behavior perfectly. But, the numerical model easily allows for all kinds of multi-axial loadings. The results and the corresponding yield loci are given in Fig. 8. Like for the principal biaxial stress combinations, the symmetry considerations are given additionally proving the symmetry considerations to be correct. When comparing the macroscopic yield function with the RVE results, the function does not fit as well when shear stresses are involved.

Fig. 8

Yield locus under all combinations of biaxial stress states in CA-based RVE (red curves). The initial yield locus of the macroscopic anisotropic yield criterion is given for comparison (blue curves). The green dashed lines are again copies of the red curves, corresponding to the symmetry assumptions in Table 1 and are inserted to verify these assumptions (Color figure online)

Comparison with Voronoi-Based Generated RVEs

For comparison, a RVE is generated by means of a centroidal Voronoi tessellation, compare (Ref 54). In Fig. 9, the RVE from grain growth simulation (b) is compared against the generated (Voronoi Tessellation based) RVE (c). The orientations of the unit cells are depicted in Fig. 9(d), for comparison the polfigure for the simulated structure is given in (a). In the generated RVE, the unit cells are randomly tilted up to $5^\circ$ against building direction and are randomly distributed. Thus, the generated texture shows no secondary anisotropy but a strong anisotropy in BD. Periodicity is enforced by construction allowing for periodic boundary conditions. To represent the columnar grained structure, the cross section is extruded in building direction like for the simulated RVE. The geometry is again kept constant but the grain orientations will be varied. Periodic boundary conditions are applied on opposite sides for the fluctuations of the displacement and the slips as well.

Fig. 9

(a) and (b) CA-based (simulated) polfigure and RVE including mesh, (c) and (d) Voronoi Tessellation-based (generated) RVE and polfigure. Grains are colored to distinguish between different grains

Using the same ansatz as in Sect. 5.4, the elastic behavior and initial yield locus are identified for the generated RVE. In Fig. 10, the Young’s modulus of generated (left) and simulated (right) RVE’s are given. The radius of the figure is scaled to the Young’s modulus in the corresponding direction and additionally coded in color. It can clearly be seen that the generated RVE with idealized texture in Fig. 10(b) is transversal isotropic symmetric with the BD as preferred direction. In comparison, the Young’s modulus of the simulated RVEs with simulated texture, depicted in Fig. 10(a), has, due to the scan strategy, in addition a secondary anisotropy. Both RVEs have a strong primary anisotropy with the lowest stiffness of $\approx$ 125 GPa in BD and the stiffest under $45^\circ$ with $\approx$228 GPa. Like for the simulated RVE, the yield loci is calculated and averaged over 10 RVEs. In Fig. 11, the mean values and standard deviations of the yield locus under biaxial loading conditions in the plane perpendicular to BD are given. The maximum of the standard deviation is similar for both RVEs, indicating that averaging over 10 RVEs with 12 grains is sufficient.

Fig. 10

Three-dimensional representation of the orientation dependent Young’s modulus

Fig. 11

Points of yielding (mean ± standard deviation) in biaxial loading in the 11–22 stress sub-space for 10 different RVEs (left generated, right simulated)

The yield loci of all biaxial combinations of stresses are calculated, and the results are depicted in Fig. 12. For comparison, the results from the cellular automaton-based RVE are given in addition. Both RVEs consist of columnar grains with a similar grain size, but differ in grain shape and texture. As a result, the homogenized behavior of both RVEs is similar in general but different in some details. When uniaxially loaded parallel to BD the narrower distribution of the unit cell orientations around BD in the Voronoi-based RVEs results in a slightly higher yield strength, as can be seen in the $\sigma _{11}/\sigma _{22}$ diagram. Another difference is when both RVEs are loaded with $\sigma _{11} = - \sigma _{22}$, the transversely isotropic symmetry of the Voronoi generated RVEs results in a higher yield strength. Since $\sigma _{11} = - \sigma _{22}$ is equal to a uniaxial shear stress in the 12 plane, this can also be seen in all diagrams containing $\sigma _{12}$. With pure shear stresses comprising the BD $\sigma _{13}/\sigma _{23}$, the influence of the secondary anisotropy can be seen directly. As a result, the circular shape in the case of transversal isotropy is slightly shifted towards a square when secondary anisotropy is involved as for the CA simulated RVEs. However, the overall plastic behavior of the CA simulated and Voronoi-based generated RVEs is very similar. Only the artificial texture which does not include the secondary anisotropy leads to small differences in the yield stresses. This confirms the observation from Ref 18 that the texture is very significant for the macroscopic behavior but the exact shapes of the grains are not.

Fig. 12

Yield locus for all combinations of biaxial stress states in a Voronoi-based RVE with idealized texture in comparison with the simulated RVE and according texture. 9 of the 15 diagrams are simulated, the remaining 6 are constructed using the previously verified symmetry assumptions from Table 1

Conclusion

A combination of a grain growth simulation and a multi-scale mechanical model was introduced to identify the process-dependent macroscopic mechanical behavior of additively manufactured IN 718.

Both, the grain growth simulation and the multi-scale mechanical model required innovative improvements to be applicable for the presented approach. The former was improved by the application of sophisticated runtime optimization algorithms and execution on a high-performance cluster computer using Message-Passing-Interface (MPI)- parallelization. These improvements enabled the numerical prediction of the resulting microstructure in unprecedented level of detail and in reasonable computation time. The application of the gradient-enhanced crystal plasticity model, validated with experimental data, on a size scale relevant for AM processes is novel and required several new strategies to ensure convergence.

The grain growth simulation was based on simulated temperature gradients during solidification and implemented in an in-house cellular automaton code. The mechanical model built on the simulated grain structures and combined gradient-enhanced crystal plasticity and numerical homogenization. This multi-scale approach was implemented in an in-house Finite Element code. Both models were validated against experimental data, and it has been shown that they are suitable for predictive simulations.

Different grain structures were simulated with the cellular automaton for the SEBM process of IN 718. The resulting textures and grain structures were used as an input to define Representative Volume Elements for the meso-mechanical simulations. Based on different realizations of the RVE, the macroscopic elastic and plastic properties were successfully identified. The macroscopic yield behavior was described by an anisotropic yield criterion whose 7 independent parameter could be identified by means of the RVE simulations. The macroscopic yield criterion was in perfect agreement with the RVE results for biaxial loading, but when shear stresses were involved, initial yielding was less well captured. This identified macroscopic mechanical behavior of IN718 depends (via the grain growth simulation and the numerical homogenization) on the process strategy during the SEBM process. Therefore, the combined numerical approach can predict the relationship between process and mechanical properties. Additionally, the results of Voronoi-based generated RVEs were evaluated in comparison with the CA simulated RVEs. It was observed that mainly the texture influences the macroscopic behavior, the grain shapes were not of major importance. Therefore, the combination of the grain growth simulation and the mechanical simulation can be simplified in a way that only the polfigures are used as an input in the RVE simulations and the grains within the RVEs are generated based on Voronoi tessellation.

So far, we focused on bulk material neglecting the influence of defects like pores and surface roughness typical for the AM process, as discussed in Ref 55, which represent other sources of size dependence. Including such geometric imperfections, surface roughness due to partially molten particles and other surface defects in the current approach, is the focus of our current research and will be discussed in a forthcoming publication.