# Crystal plasticity model of residual stress in additive manufacturing

Nicolò Grilli  Daijun Hu  Dewen Yushu  Fan Chen  Wentao Yan
<EMAIL_ADDRESS>Department of Mechanical Engineering, University of
Bristol, Queen’s Building, University Walk, Bristol, UK  Department of
Mechanical Engineering, National University of Singapore, Singapore 117575,
Singapore  Computational Mechanics and Materials Department, Idaho National
Laboratory, Idaho Falls, ID 83415, USA

###### Abstract

Selective laser melting is receiving increasing interest as an additive
manufacturing technique. Residual stresses induced by the large temperature
gradients and inhomogeneous cooling process can favour the generation of
cracks. In this work, a crystal plasticity finite element model is developed
to simulate the formation of residual stresses and to understand the
correlation between plastic deformation, grain orientation and residual
stresses in the additive manufacturing process. The temperature profile and
grain structure from thermal-fluid flow and grain growth simulations are
implemented into the crystal plasticity model. An element elimination and
reactivation method is proposed to model the melting and solidification and to
reinitialise state variables, such as the plastic deformation, in the
reactivated elements. The accuracy of this method is judged against previous
method based on the stiffness degradation of liquid regions by comparing the
plastic deformation as a function of time induced by thermal stresses. The
method is used to investigate residual stresses parallel and perpendicular to
the laser scan direction, and the correlation with the maximum Schmid factor
of the grains along those directions. The magnitude of the residual stress can
be predicted as a function of the depth, grain orientation and position with
respect to the molten pool.

###### keywords:

Additive manufacturing , 316 stainless steel , Residual stress , Finite
element method , Crystal plasticity , Grain structure

††journal:

##  1 Introduction

Additive manufacturing (AM) is becoming increasingly popular in industries.
Metal AM processes such as selective laser melting (SLM) [1], selective
electron beam melting (SEBM) [2] and selective laser sintering (SLS) [3] are
layered approaches fabricating the products with complex shapes [4]. One of
the bottlenecks for metal AM is the thermal stress caused by the high
temperature gradients during the high-frequency heating/cooling cycles. High
residual stress may result in defects [5, 6] such as the part distortion [7],
delamination [8] and poor fracture resistance [9, 10, 11]. The investigation
of temperature history and stress evolution is therefore necessary for
optimising the manufacturing parameters during the AM process.
Due to the large temperature gradients and high-frequency heating/cooling
cycles during the manufacturing process, it is difficult to monitor the
temperature and thermal stress in real time during experiments. For the
residual stress measurement, the experimental approaches, such as the neutron
diffraction [12, 13], X-ray diffraction [14], contour method [15] and semi-
destructive hole drilling method [16, 17], are costly and sometimes limited to
surface measurements. Numerical simulations are efficient to study the
mechanical behaviour during the manufacturing process. The finite element
method (FEM) has been widely used in simulating thermal stresses in AM [18].
One type of such models is based on the inherent strain method, where the
deposited regions are activated with predefined strains to calculate the final
deformation and residual stress of the part [19]. These models are
computationally efficient, but the accuracy may be compromised due to the
uncertainties of the inherent strains. Another type is thermo-mechanical
coupling, where the temperature profile is calculated by solving the governing
equations of heat transfer and then serves as the load for thermal stress
calculation [20, 21]. Thermo-mechanical models are mostly constructed based on
many geometrical simplifications and physical assumptions as the heat transfer
part does not incorporate the thermal-fluid flow effects [22]. Most of the FEM
simulations failed to resolve the track morphology and molten pool evolution
as all the tracks are assumed to be uniform. Therefore, those simulations are
only able to predict the overall part distortion at the macroscopic scales,
but cannot explain the thermal stress at the microscale.
An emerging type of the thermal stress models is the coupling between the
computational fluid dynamics (CFD) and FEM models [23, 24, 25]. These coupled
CFD-FEM models analyse the thermal stress and deformation using the high-
accuracy temperature profiles output from the CFD simulations, particularly
the powder-scale thermal-fluid flow simulations [26], instead of those from
simplified analytical models [27]. Cheon et al. [24] proposed a CFD-FEM
framework for welding to analyse the thermal stress along a single welding
line, but the model did not consider the residual stress reset due to melting.
Bailey et al. [23] implemented the temperature and surface profiles from
thermal-fluid flow simulations into an FEM residual stress model to simulate
the multi-track manufacturing cases of the laser directed energy deposition
(DED) process. The element activation and deactivation method in ABAQUS was
used to simulate the melting and solidification, which is troublesome to input
command lines to activate elements of the solidified materials in the input
file and frequently encounters problems on numerical convergence and
numerically-induced deformation of the newly-activated elements. In our
previous work [25], we simulated the multi-layer multi-track SLM process by
mapping the temperature profiles from the high-fidelity thermal-fluid flow
simulation into the FEM model, enabling the high resolution of the molten pool
evolution and the morphology (rough surfaces and inner voids) of the FEM
model. The quiet element method was used to simulate the melting and
solidification, where the melted materials are assigned with null material
properties.
There are two critical but common problems in the aforementioned thermo-
mechanical models and CFD-FEM models:

  1. 1.

How to represent the melting and solidification while avoiding numerically-
induced deformation and numerical divergence. In the previous models, two
methods are commonly used: first, the element activation method [28, 29, 30,
31, 32] can activate elements for solidified materials but cannot deactivate
elements again due to remelting, while the surface nodal temperature obtained
by the interpolation between active and inactive elements [33] can cause
significant inaccuracy and divergence problems. Second, the quiet element
method [32, 33, 34], which is essentially the same as the stiffness
degradation method commonly used in fracture mechanics [35, 36, 37, 38, 39],
can remove the effect of the liquid phase on the solidified region, however
upon solidification the liquid phase becomes solid with initial deformation
that is incompatible with the deformation of the surrounding solid region,
which induces un-physical deformation and may also cause numerical divergence.

  2. 2.

None of the previous models has considered the detailed grain structures and
the resultant temperature-dependent anisotropic mechanical properties at the
microscale. Using simplified homogenized constitutive laws, the previous
models cannot resolve the micro-scale residual stress [40], which can result
in unique mechanical properties of the AM parts significantly differing from
those by the conventional manufacturing [41].

In this paper, a crystal plasticity finite element model (CPFEM) [42, 43] is
developed, incorporating the temperature profile from the thermal-fluid flow
simulation [25] and grain structures from the grain growth simulation [44, 45]
and experimental data [40], as to be described in Section 2. Particularly, an
element elimination and reactivation method is proposed to model the melting
and solidification. More importantly, the algorithm can reinitialise the state
variables, as tested in a simple case in Section 3, such as the plastic
deformation of the solidifying elements, which is critical to avoid the pre-
existing plastic deformation just after solidification. Simulation cases of
SLM of 316 stainless steel are conducted as described in Section 4. The
proposed element elimination and reactivation method and the commonly used
stiffness degradation method are compared in terms of prediction accuracy,
convergence and computational efficiency. Based on the simulation results, the
formation of residual stresses in different grains and the correlation between
plastic deformation, depth and residual stress are discussed in Section 5.

##  2 Material model

###  2.1 Crystal plasticity framework

A crystal plasticity material model including thermal response is used for
this work. The deformation gradient $\bm{F}$ can be decomposed into elastic,
thermal and plastic parts [46, 47, 48, 49]:

$\bm{F}=\bm{F}_{e}\bm{F}_{th}\bm{F}_{p}\ ,$ (1)

where $\bm{F}_{e}$ is the elastic part including the stretching and rigid body
rotation, $\bm{F}_{th}$ is thermal part and $\bm{F}_{p}$ is the irreversible
plastic deformation, which evolves according to [50, 51]:

$\bm{L}_{p}=\dot{\bm{F}}_{p}\bm{F}_{p}^{-1}=\sum_{\alpha=1}^{n}\dot{\gamma}^{\alpha}\bm{m}^{\alpha}\otimes\bm{n}^{\alpha}\
,$ (2)

where $\dot{\gamma}^{\alpha}$ is the shear rate of slip system $\alpha$. The
vectors $\bm{m}^{\alpha}$ and $\bm{n}^{\alpha}$ are unit vectors which
describe the slip direction and normal of slip plane $\alpha$, respectively.
$n$ is the number of slip systems in active state.
The plastic strain rate on slip system $\alpha$ of face-centered cubic metals
can be calculated by a power-law equation [52, 53]:

$\dot{\gamma}^{\alpha}=\dot{\gamma}_{0}\left|\frac{\tau^{\alpha}}{\tau_{c}^{\alpha}(T)}\right|^{\frac{1}{m}}\rm{sgn}(\tau^{\alpha})\
,$ (3)

where $\tau^{\alpha}$ is the resolved shear stress, $\tau_{c}^{\alpha}$ is the
temperature dependent critical resolved shear stress (CRSS),
$\dot{\gamma}_{0}$ and $m$ are constants, which indicate the reference strain
rate and rate sensitivity of slip, respectively.
An exponential model calibrated by experiments can be used to calculate the
CRSS change with temperature decrease of each slip system [54, 55]:

$\tau_{c}^{\alpha}(T)=\left(k_{A}+k_{B}\cdot{\rm
exp}[-k_{C}(T-T_{0})]\right)\cdot\tau_{c}^{\alpha}(T_{0})\ ,$ (4)

where $\tau_{c}^{\alpha}(T)$ is the CRSS of slip system $\alpha$ at
temperature $T$, $T_{0}$ is the reference temperature, $k_{A}$, $k_{B}$ and
$k_{C}$ are constants. For slip system $\alpha$, its hardening behaviour is
influenced by other slip systems $\beta$ as [56]:

$\dot{\tau}_{c}^{\alpha}=\sum_{\alpha=1}^{n}h_{\alpha\beta}\left|\dot{\gamma}^{\beta}\right|\
,$ (5)

where $h_{\alpha\beta}$ is the hardening matrix given by [57, 58]:

$h_{\alpha\beta}=q_{\alpha\beta}\left[h_{0}\left(1-\frac{{\tau_{c}^{\alpha}}}{\tau_{\textrm{sat}}}\right)^{a}\right]\
,$ (6)

where $q_{\alpha\beta}$ is a hardening coefficient matrix for latent
hardening, $h_{0}$ is an initial hardening term, $\tau_{sat}$ is the
saturation slip resistance and $a$ is a constant. The parameters of the
crystal plasticity model are calibrated using the experimental data in [40].
For the thermo-mechanical response during the laser scan process, the Green-
Lagrange strain tensor can be expressed by the equation composed of thermal
and elastic deformation gradients as:

$\bm{E}_{e}=\frac{1}{2}\left(\bm{F}_{\textrm{th}}^{T}\bm{F}_{e}^{T}\bm{F}_{e}\bm{F}_{\textrm{th}}-\bm{I}\right)\
,$ (7)

where $\bm{I}$ is the identity matrix. To describe how the size of the
substance changes with the temperature, the linear thermal expansion
coefficient $\alpha_{l}$ of 316L SS is introduced [59]. The volumetric thermal
expansion coefficient can be expressed by $\alpha_{v}=3\alpha_{l}$ [60]. Since
the thermal expansion coefficient of 316L stainless steel varies slightly with
temperature [61, 62], and the temperature dependence is given by fitting the
experimental data into a linear equation [63]:

$\alpha_{v}=\alpha_{0}+\alpha_{1}(T-T_{0})\ ,$ (8)

where $\alpha_{0}$ denotes the thermal expansion coefficient at room
temperature and $\alpha_{1}$ is the linear coefficient of the temperature
dependence of $\alpha_{v}$. The infinitesimal volumetric change $dV$ from
${T_{0}}$ to ${T}$ can be described by [48]:

$\alpha_{v}\mathrm{d}T=\frac{\mathrm{d}V}{V}\ ,$ (9)

where $V$ is the volume of the region and $T_{0}$ is the reference temperature
(room temperature). Substituting equation (8) into equation (9), the thermal
eigenstrain tensor $\bm{\alpha}$ can be expressed by [64]:

$\bm{\alpha}=\frac{1}{2}\left[{\rm
exp}\left(\frac{1}{3}\alpha_{1}(T-T_{0})^{2}+\frac{2}{3}\alpha_{0}(T-T_{0})\right)-1\right]\bm{I}\
.$ (10)

The 2nd Piola-Kirchhoff stress can be calculated by:

$\bm{S}=\mathbb{C}(\bm{E_{e}}-\bm{\alpha})\ ,$ (11)

where $\mathbb{C}$ is the rank four elasticity tensor of 316L stainless steel
[65, 66]. The elasticity tensor of the solid is temperature dependent [54],
which can be described by a linear function for the components:

$\mathbb{C}_{ij}^{\rm{solid}}(T)=\mathbb{C}_{ij}^{\rm{solid}}(T_{0})+\frac{\mathrm{d}\mathbb{C}_{ij}^{\rm{solid}}}{\mathrm{d}T}(T-T_{0})\
,$ (12)

where $\mathbb{C}_{ij}^{\rm{solid}}(T_{0})$ is the elasticity tensor component
at room temperature, and the derivative term indicates the change of the
components caused by temperature change. Phase change can be modelled by
changing the stiffness tensor $\mathbb{C}$, as to be explained in section 2.2.
Since positive strain represents tension and negative strain represents
compression, the same applies to the stress tensor components. In the
following, compression will appear as a negative stress component while
tension will appear as a positive stress component.
The temperature $T$ is obtained from thermal-fluid flow simulations [25]. It
should be mentioned that the temperature mapping from thermal-fluid to
mechanical simulation is one-way, and the heat transfer analysis is merely
conducted in the high-fidelity thermal-fluid flow simulation, from which the
thermal information are extracted for the crystal plasticity calculation.
Thus, deactivation of the melted parts will not affect the temperature
distribution in the representative volume.

###  2.2 Residual stiffness method

During the additive manufacturing process, the phase transitions cause
significant change in the mechanical properties, which should be taken into
consideration. Since in our simulation the temperature profiles are obtained
by mapping and interpolating the thermal-fluid flow simulation results to the
finite element model in time and space [67, 25], there are liquid phase and
gas phase existing in the model, which represent the melted material and gas
above molten pool, respectively. Therefore, a residual stiffness method based
on the temperature field is applied here to prevent sudden changes in
mechanical properties caused by the temperature dependence of stiffness and
avoid the convergence problem of simulation. The specific method is to set a
temperature range between the melting point $T_{m}$ of 316L stainless steel
and the gas temperature $T_{g}$ set in the thermal-fluid flow simulation. In
this temperature range, the stiffness tensor of the solid material is
$\mathbb{C}_{ij}^{\rm{solid}}(T)$, which conforms to equation (12).
The temperature from the thermal-fluid flow simulation is recorded every $45$
$\mu$s and these time points are called the “CFD output” or “temperature
output”. The CPFEM simulations have smaller time steps, therefore an
interpolation is necessary. For instance, in the time interval
$\left[t_{1},t_{2}\right]$ the temperature is interpolated linearly from
$T_{1}$ to $T_{2}$, as shown in Figure 1.
A phase identification algorithm has been implemented: if the temperature is
higher than the melting temperature $T_{m}$ or lower than the gas temperature
$T_{g}$, it means that the material is in liquid or gas phase. This algorithm
transforms the stiffness matrix $\mathbb{C}_{ij}$ linearly with time between
those of phase $1$ at $t_{1}$ and phase $2$ at $t_{2}$ in Figure 1. The
formulas to find the stiffness tensor are listed in Table.1. A residual
stiffness coefficient $q_{r}$ is introduced to describe the residual stiffness
tensor of the material:

$\mathbb{C}_{ij}^{\rm{residual}}=q_{r}\cdot\mathbb{C}_{ij}^{\rm{solid}}(T_{0})\
.$ (13)

The stiffness tensor $\mathbb{C}_{ij}$ in Table 1 is used in equation (11) and
the interpolation avoids discontinuities in the stiffness tensor that would
create convergence problems.

Figure 1:  Schematic of the phase identification algorithm to implement the residual stiffness method for liquid and gas phases. Table 1: Time evolution of the stiffness tensor depending on the phase state at time $t_{1}$ and $t_{2}$. Phase 1 | Phase 2 | Stiffness
---|---|---
Solid | Solid | $\mathbb{C}_{ij}=\mathbb{C}_{ij}^{\mathrm{solid}}\left(T\right)$
Solid | Gas/Liquid | $\mathbb{C}_{ij}=\frac{(t_{2}-t)}{(t_{2}-t_{1})}\mathbb{C}_{ij}^{\mathrm{solid}}\left(T_{1}\right)+\frac{(t-t_{1})}{(t_{2}-t_{1})}\mathbb{C}_{ij}^{\mathrm{residual}}$
Gas/Liquid | Solid | $\mathbb{C}_{ij}=\frac{(t_{2}-t)}{(t_{2}-t_{1})}\mathbb{C}_{ij}^{\mathrm{residual}}+\frac{(t-t_{1})}{(t_{2}-t_{1})}\mathbb{C}_{ij}^{\mathrm{solid}}\left(T_{2}\right)$
Gas/Liquid | Gas/Liquid | $\mathbb{C}_{ij}=\mathbb{C}_{ij}^{\mathrm{residual}}$

The equations above are implemented into the MOOSE finite element framework
[68], where a Newton-Raphson approach to find the Cauchy stress tensor
$\bm{\sigma}$ is applied. To obtain the 2nd Piola-Kirchhoff stress in equation
(11), a function $\bm{\psi}=\bm{S}-\mathbb{C}(\bm{E_{e}}-\bm{\alpha})$ is
defined, the absolute value of which must be minimised. An iterative return
mapping algorithm, in which the stress $\bm{S}$ is the variable is
implemented. The stress increment during the return mapping algorithm can be
written as [69, 70]:

$\bm{S}_{i+1}=\bm{S}_{i}-\left(\frac{\partial\bm{\psi}}{\partial\bm{S}}\right)^{-1}\cdot\bm{\psi}\
.$ (14)

where $i$ and $i+1$ indicate two subsequent iterations. The Cauchy stress
$\bm{\sigma}$ is obtained by [71]:

$\bm{\sigma}=\frac{1}{\rm{det}(\bm{F_{e}})}\cdot\bm{F_{e}}\cdot\bm{S}\cdot\bm{F_{e}^{T}}\
.$ (15)

The Cauchy stress is used by the finite element solver to find the
displacement field $\mathbf{u}$ that leads to stress equilibrium.
The parameters used in the simulations and the elastic constants of 316L
stainless steel in Voigt notation are listed in Table 2.

Table 2: Material and model parameters used in the simulations [54, 40, 63, 72, 73] Reference plastic strain rate ($\dot{\gamma}_{0}$) |  $1\times 10^{-7}$ s-1
---|---
Plastic strain exponent constant ($m$) | 0.1
Hardening matrix ($h_{\alpha\beta}$) | 3839 MPa
Hardening exponent ($a$) | 2.5
Reference temperature ($T_{0}$) | 303 K
Melting temperature ($T_{m}$) | 1648.15 K
Gas temperature ($T_{g}$) | 298.1 K
Critical resolved shear stress at $T_{0}$ ($\tau_{c}^{\alpha}$) | 213 MPa
Saturation slip resistance ($\tau_{sat}$) | 302 MPa
Volumetric thermal expansion coefficient at $T_{0}$ ($\alpha_{0}$) | 44.73$\times 10^{-6}$ K-1
Derivative of volumetric thermal expansion coefficient ($\alpha_{1}$) | 0.01011$\times 10^{-6}$ K-2
Temperature dependence of CRSS ($k_{A}$) | 0.53
Temperature dependence of CRSS ($k_{B}$) | 0.47
Constant in temperature dependence eq.(4) of CRSS ($k_{C}$) | 0.008
Elastic constants at $T=T_{0}$ ($\mathbb{C}_{11}$) | 204.6 GPa
Elastic constants at $T=T_{0}$ ($\mathbb{C}_{12}$) | 137.7 GPa
Elastic constants at $T=T_{0}$ ($\mathbb{C}_{44}$) | 126.2 GPa
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{11}/{\mathrm{d}T}$) | -90.33 MPa K-1
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{12}/{\mathrm{d}T}$) | -45.10 MPa K-1
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{44}/{\mathrm{d}T}$) | -51.78 MPa K-1
Residual stiffness coefficient ($q_{r}$) | 0.01 and 0.001
Laser scan speed | 700 mm/s
Laser power | 150 W

###  2.3 Implementation of the element elimination and reactivation method

To model the melting and solidification of material in the SLM process, an
element elimination and reactivation approach is employed. This approach is
based on a moving subdomain paradigm that is originally developed [74], which
is based on MOOSE’s subdomain restricted system [68]. Specifically, the
physical domain is divided into an active subdomain ($\Omega_{a}$) and an
inactive subdomain ($\Omega_{i}$). The thermo-mechanical finite element
computation is defined and carried out in $\Omega_{a}$, while $\Omega_{i}$
only carries the geometric information. In other words, the elements in the
inactive domain, or inactive elements, do not contribute to the equation
system that calculates stress equilibrium.
The active and inactive subdomains are determined only by the prescribed
temperature field. At a typical FEM time step, any active element
$\mathcal{T}\in\Omega_{a}$ with an average temperature that is above the
melting temperature (i.e., $T_{\text{avg}}>T_{m}$) or below the gas
temperature (i.e., $T_{\text{avg}}<T_{g}$) is moved from the active subdomain
to the inactive one, i.e.:

$\Omega_{a}\rightarrow\Omega_{a}\backslash\mathcal{T},\quad\Omega_{i}\rightarrow\Omega_{i}\cup\mathcal{T}.$
(16)

This process mimics the solid material being eliminated from the physical
domain, which happens when the material is melted around the laser beam, and
is therefore referred to as element elimination. Similarly, when the average
temperature of one inactive element $\mathcal{T}\in\Omega_{i}$ falls to the
solid temperature range (i.e., $T_{g}\leq T_{\text{avg}}\leq T_{m}$), it is
moved to the active subdomain, i.e.:

$\Omega_{a}\rightarrow\Omega_{a}\cup\mathcal{T},\quad\Omega_{i}\rightarrow\Omega_{i}\backslash\mathcal{T}.$
(17)

This mimics the process of material being added back to the solid physical
domain, which happens when the material solidifies behind the laser beam, thus
is referred to as element reactivation.

After elements being transferred between the two subdomains, several tasks
need to be accomplished before the subsequent computation. First, the boundary
information is updated based on the coverage of the new subdomains. Here, the
subdomain boundaries are updated for both $\Omega_{a}$ and $\Omega_{i}$. The
boundary conditions are projected to the updated boundaries correspondingly.
Second, the initial conditions of the displacement variables and the material
properties that define the crystal plasticity behaviour, such as the plastic
deformation (see Section 2.1), are projected to the quadrature points of the
reactivated elements. Projection of the initial condition indicates that the
displacement, strain, and stress calculations are reset to zero. This allows
for the restart of the plastic deformation calculation during the simulation,
which cannot be implemented easily in the stiffness degradation method [75].
Indeed, if the elements are not eliminated and reactivated, they can still
carry a certain amount of load. In the stiffness degradation method, a sudden
change of the plastic deformation due to melting and solidification causes
convergence problems because an elastic deformation jump, and consequent
stress discontinuity, occur.
Finally, the element elimination and reactivation method resolves the issue
that is brought by large deformation from the reactivated elements if they are
set to carry the displacement history upon reactivation.

##  3 Test case: element elimination during tension

The element elimination and reactivation method is tested on a simple example
geometry, as shown in Figure 2. The aim of this test is to demonstrate that
inactive elements do not contribute to the stress equilibrium calculations,
i.e. they do not exert stress on the active elements, and to prove that
boundary conditions are implemented correctly on the continuously changing
free surface.
A cubic representative volume is constituted of 10 $\times$ 10 $\times$ 10
elements and has a size of 1 $\mu$m $\times$ 1 $\mu$m $\times$ 1 $\mu$m.
Boundary condition $\mathbf{u}_{x}=0$ is applied on the surface $x=0$,
$\mathbf{u}_{y}=0$ on $y=0$ and $y=1$ $\mu$m, $\mathbf{u}_{z}=0$ on $z=0$ and
$z=1$ $\mu$m. Displacement boundary condition is imposed on the surface $x=1$
$\mu$m. $\mathbf{u}_{x}$ is increased up to $0.001$ $\mu$m on that surface
before the element elimination algorithm starts.
The element elimination algorithm is based on a prescribed temperature field,
which is chosen in such a way that elements from the surface $z=1.0$ $\mu$m
are progressively removed, down to the surface $z=0.5$ $\mu$m. The
displacement $\mathbf{u}_{x}$ on the surface $x=1$ $\mu$m is kept constant
during element elimination, as shown in Figure 2(a)-(c), therefore
$\bm{\varepsilon}_{xx}$ remains constant at the value of $0.001$.
Plasticity and thermal eigenstrain are not activated in this simulation case,
and the elastic constants are set as $\mathbb{C}_{11}=1.5$ GPa and
$\mathbb{C}_{12}=0.75$ GPa. $\bm{\varepsilon}_{yy}$ and the off diagonal
strain components are zero. When no element elimination has taken place, the
boundary condition $\mathbf{u}_{z}=0$ on the surface $z=1$ $\mu$m holds,
therefore $\bm{\varepsilon}_{zz}=0$ and $\mathbf{u}_{z}=0$ everywhere in the
representative volume, as shown in Figure 2(d). This lateral constraint is
reflected by the large positive value of the stress component
$\bm{\sigma}_{zz}$ in Figure 2(g). Therefore, a tensile stress is present on
that surface.

(a) (a)

(b) (b)

(c) (c)

(d) (d)

(e) (e)

(f) (f)

(g) (g)

(h) (h)

(i) (i)

Figure 2:  Test case for the element elimination and reactivation method.
(a),(b),(c) Firstly, a displacement $\mathbf{u}_{x}$ is applied and elements
are eliminated. (d) The lateral displacement starts at zero, (e)-(f) then
lateral contraction appears due to the newly created free surface
perpendicular to the $z$ axis. For the same reason, (g) the lateral stress
$\bm{\sigma}_{zz}$ starts at 750 MPa, (h)-(i) then $\bm{\sigma}_{zz}$ becomes
zero.

After the element elimination algorithm starts, the top surface perpendicular
to the $z$ axis becomes unconstrained. Therefore, lateral contraction takes
place because of the constant displacement along the $x$ axis. The lateral
contraction can be calculated as follows:
$\bm{\varepsilon}_{zz}=-\left(\mathbb{C}_{12}/\mathbb{C}_{11}\right)\bm{\varepsilon}_{xx}=-0.0005$.
For example, in Figure 2(e), 8 elements are still present along the $z$ axis,
therefore the lateral displacement on the top surface is given by
$\mathbf{u}_{z}=\bm{\varepsilon}_{zz}\cdot 0.8$ $\mu$m $=-0.0004$ $\mu$m. The
same holds in Figure 2(f), in which the lateral displacement on the top
surface is $\mathbf{u}_{z}=-0.00025$ $\mu$m because 5 elements are still
present along the $z$ axis.
As shown in Figures 2 (h)-(i), there is no lateral stress $\bm{\sigma}_{zz}$
after the element elimination algorithm starts. This proves that the method
developed creates new free surfaces in the representative volume and the
eliminated elements do not apply stress on the active elements.

##  4 Simulation of selective laser melting

###  4.1 Simulation setup

Crystal plasticity simulations of selective laser melting are carried out
using the representative volume in Figure 3. A single scan track is modelled.
The dimensions along the X, Y and Z axes are $200$ $\mu$m, $160$ $\mu$m and
$108$ $\mu$m respectively. The FEM model represents a small region inside a
larger sample. However, it is representative of the texture because of the
large number of grains included.
The temperature field is obtained from thermal-fluid flow simulations [76, 77,
25]. The governing equations are the conservation of mass, momentum and
energy. Most of the physical factors are incorporated, including the ray-
tracing heat source model for laser, thermal conduction, surface radiation and
convection, latent heat, evaporation, recoil pressure, viscosity, buoyancy
force, surface tension and Marangoni effect. The details of the model and
material parameters are reported in our previous works [76, 77].
The grain structure is obtained from a phase field model for grain growth
[44], where heterogeneous nucleation [78] and the initial grain structures in
powder particles and substrate are incorporated, and the same temperature
profile from the thermal-fluid flow simulation is used. The phenomena,
including grain nucleation and growth, competitive growth, epitaxial growth
from powder particles and substrate, and grain coarsening in heat affected
zones (HAZs), are comprehensively considered and experimentally validated. The
details of the model and parameters are reported in our previous work [44].
A total of 4522 grains are present with different orientations. The voxels
obtained from the phase field model for grain growth are mapped directly into
the FEM mesh, therefore the shape of the grains is reproduced accurately. The
Euler angles $\varphi_{1}$, $\Phi$ and $\varphi_{2}$ are assigned based on
electron backscatter diffraction (EBSD) experiments provided by Dr. Yin Zhang
and published elsewhere [40] to reproduce the realistic texture of AM 316
stainless steel. This is important because the grain orientation determines
both the macroscopic and microscopic mechanical properties of the sample. On
the other hand, the 3D shape of the grains cannot be fully extracted from EBSD
measurements, therefore simulation results from the grain growth model are
adopted. Figure 3 shows the value of $\varphi_{1}$. The average grain size is
approximately $16$ $\mu$m. However, grains in the centre of the laser track
tend to be elongated along the $z$ axis, as shown in Figure 3.
The temperature field obtained from thermal fluid-flow simulations is
projected on the mesh of the FEM simulations. Melting in the crystal
plasticity simulations is described either by the stiffness degradation model
described in section 2.2 or by the element elimination and reactivation method
described in section 2.3. The temperature fields at different time steps are
shown in Figure 4 (a)-(c), in which the eliminated elements are hidden from
view. A structured mesh made of first order hexahedral elements is used.
Simulations are carried out with the element size being $4$ $\mu$m or $8$
$\mu$m, respectively. The numbers of elements are 54000 and 7000 respectively.

Figure 3:  Representative volume for the selective laser melting simulations
and grain structure implemented from the phase field simulation for the
crystal plasticity model.

The boundary conditions used in the simulations represent the constraints that
are applied by the substrate and by the surrounding material. The displacement
components are zero on the surfaces $x=0$, $x=200$ $\mu$m, $y=0$ and $y=160$
$\mu$m. The vertical displacement $\mathbf{u}_{z}$ is zero on the bottom
surface, $z=0$, to avoid translation along the $z$ axis. The displacement is
free on the upper surface ($z=108$ $\mu$m), where the laser beam is present.

(a) (a)

(b) (b)

(c) (c)

(d) (d)

(e) (e)

(f) (f)

Figure 4:  Snapshots of the temperature field using the element elimination
and reactivation method at (a) $t=99$ $\mu$s, (b) $t=226$ $\mu$s, (c) $t=451$
$\mu$s; (d)-(f) corresponding stress component $\bm{\sigma}_{xx}$ using the
stiffness degradation method.

###  4.2 Simulation convergence

The convergence speed of the simulations using the stiffness degradation
method depends mainly on the residual stiffness parameter $q_{r}$. A value
smaller than 0.001 significantly increases the simulation time, but there is
very small difference in the plastic deformation between the simulations with
$q_{r}=0.01$ and $q_{r}=0.001$, as will be shown in the following. Simulations
with both $q_{r}=0.01$ and $q_{r}=0.001$ are carried out.
If the element elimination method is used, the time step must be decreased
immediately before and after the temperature steps at which element
elimination takes place, as described in section 2. In this way, good
convergence can be obtained. Specifically, the time step is reduced to $0.1$
$\mu$s in an interval of $2$ $\mu$s around each temperature step, i.e. every
$45$ $\mu$s. For instance, considering the first temperature step, the time
step is reduced to $0.1$ $\mu$s at $t=44$ $\mu$s and it is increased again to
$1$ $\mu$s after $t=46$ $\mu$s. This time step reduction applies only during
the laser scan, while time step is kept constant at $1$ $\mu$s during the
subsequent cooling stage, which is the second part of the simulation, carried
out after the laser beam has exited the geometry. By contrast, in the
stiffness degradation method, the value of the time step is always $1$ $\mu$s.
Since the eliminated elements do not contribute to the residual or Jacobian
calculation of the solver, the computation of each time step is slightly
faster if the element elimination method is used. However, a reduction of the
time step size before and after the element elimination steps is necessary in
that case. Therefore, the simulations carried out with the element elimination
method are not necessarily faster, unless the number of eliminated elements is
a significant fraction of the total number of elements. This is shown by the
computational costs reported in table 3.

(a) (a)

(b) (b)

(c) (c)

(d) (d)

Figure 5:  Vertical displacement $\mathbf{u}_{z}$ at $t=771$ $\mu$s: (a) just before element reactivation, (c)-(d) just after element reactivation in case the reinitialisation $\mathbf{u}_{z}=0$ or $\mathbf{u}_{z}=3$ $\mu$m is used respectively. (b) Element distortion due to the reinitialisation $\mathbf{u}_{z}=0$. Table 3: Computational costs of the simulations. Mesh size | CPUs | Element elimination | Stiffness degradation
---|---|---|---
4 $\mu$m | 16 | 99 hours | 79 hours
8 $\mu$m | 12 | 9.5 hours | 9.5 hours

As stated before, the element elimination steps require a reduction of the
time step, however, they do not cause further convergence problems. The time
steps at which element reactivation takes place leads to the following
problem. Because of thermal expansion, there is a significant displacement
along the $z$ axis on the top surface, as shown in Figure 5(a). Normally,
reactivated elements have zero displacement because the displacement variable
is reinitialised. In the MOOSE framework, the displacement in the
reinitialised nodes is assigned in such a way as to get an average zero
displacement at the eight integration points, as shown in Figure 5(b).
Therefore, when an element is reactivated with zero displacement, the
excessive distortion leads to divergence, as shown in Figure 5(c).
Several attempts have been carried out in which the initial conditions on the
displacement after element reactivation are changed as a function of the
coordinates. Surprisingly, a simple solution is to reinitialise
$\mathbf{u}_{z}=3$ $\mu$m when elements are reactivated, while
$\mathbf{u}_{x}$ and $\mathbf{u}_{y}$ are reinitialised to zero everywhere.
This displacement value is closer to the average displacement of the molten
pool region on the top surface in the simulation using the residual stiffness
method, while it also leads to good convergence. Other reinitialisation values
for $\mathbf{u}_{z}$, such as 1 $\mu$m, 2 $\mu$m, 5 $\mu$m and $6$ $\mu$m have
been used, but they have lead to excessive element distortion just after
reinitialisation. Only values between 3 $\mu$m and 4 $\mu$m remove the element
distortion problem, as shown in Figure 5(d). This implementation is made
possible by the MOOSE framework because the initial conditions can be set as a
function of time. At $t=0$, the displacement is initialised to zero, while for
$t>0$, the components $\mathbf{u}_{x}$, $\mathbf{u}_{y}$ are reinitialised to
zero and $\mathbf{u}_{z}$ is reinitialised to $3$ $\mu$m everywhere in the
geometry. This strategy can be applied to any additive manufacturing
simulation and can be extended by setting the initial conditions as a more
complex function of space and time.

###  4.3 Comparison between the computational methods

The two computational methods presented in sections 2.2 and 2.3 are compared.
The plastic deformation is a very important quantity because it determines the
residual stress after the laser scan. The diagonal components of
$\bm{F}_{p}-\bm{I}$ are averaged over the central part of the representative
volume, constituted of the highlighted elements in Figure 6(d). This procedure
is carried out for simulations using the element elimination method and the
stiffness degradation method.
Figures 6(a)-(c) show the comparison of simulation results by the two methods.
During laser scan, the hot region undergoes thermal expansion and applies
compression on the neighbouring regions. Therefore, plastic compression along
the $x$ and $y$ axes takes place in the central region. Since plastic
deformation is isochoric, the elements expand along the $z$ axis, as shown in
Figure 6(c).
The plastic deformation is always larger using the element elimination and
reactivation method. If the stiffness degradation method is used, the residual
plastic deformation component along the $y$ axis, after laser scan and
cooling, is underestimated by about 20%. The difference is particularly
important at time $400$-$500$ $\mu$s, when the laser just passes the
representative volume but the melting pool is still deep. In this time
interval, the stiffness degradation model underestimates the plastic
deformation by a factor 2.

(a) (a)

(b) (b)

(c) (c)

(d) (d)

Figure 6:  (a)-(c) Plastic strain components in the centre of the
representative volume, averaged over the highlighted elements in (d).

Overall, in the stiffness degradation method, the presence of low stiffness
elements prevents the elements near the melting pool to accommodate part of
the plastic deformation. Simulations with both $q_{r}=0.01$ and $q_{r}=0.001$
are carried out and the plastic deformation along the $y$ axis is shown in
Figure 7 when the laser beam is approximately in the centre of the
representative volume. The strong similarity between the plastic deformation
in the two simulations confirms that the range of values chosen for the
parameter $q_{r}$ does not affect the simulation results. The other components
of the plastic deformation gradient are also very similar, therefore a value
$q_{r}=0.001$ is suitable for this kind of simulations.

(a) (a)

(b) (b)

Figure 7:  Component $\left(\bm{F}_{p}\right)_{yy}-1$ of the plastic
deformation gradient at time $t=179$ $\mu$s using the residual stiffness
method with (a) $q_{r}=0.01$ and (b) $q_{r}=0.001$.

###  4.4 Grain orientation and plastic deformation

The element elimination and reactivation method will be used in the following
simulation results to understand the correlation between plastic deformation,
residual stress and grain orientation. Firstly, the correlation between
residual plastic deformation (after the laser scan and cooling) and grain
orientation is analysed. The scatter plots in Figure 8(a)-(c) are obtained as
follows: the plastic deformation components are averaged over each element in
the representative volume and one element is selected every 16 $\mu$m and
represented as a point in the scatter plots. The distribution of these
elements is uniform in the representative volume and this algorithm allows to
select a representative set of elements that do not fill completely the
scatter plot. Given the Euler angles in each selected element, the maximum
Schmid factor among the slip systems is calculated for load along a particular
direction ($x$, $y$ or $z$ axes respectively) [79]. This is motivated by the
fact that, despite of the complex load undergone by different regions, the
stress component along one direction should in principle induce plastic
deformation along that particular direction. The colour of the points in
Figure 8(a)-(c) corresponds to the position of the element along the $z$ axis.

(a) (a)

(b) (b)

(c) (c)

Figure 8:  Correlation between maximum Schmid factor among the slip systems
and plastic deformation after laser scan.

(a) (a)

(b) (b)

(c) (c)

Figure 9:  Correlation between maximum Schmid factor among the slip systems
and the range of plastic deformation after laser scan.

As shown in Figure 8, the maximum plastic deformations occur on some points
that are at an intermediate depth (about 40 $\mu m$ to 80 $\mu m$). This is
related to the solidification and thermal stress around the molten pool, as
also revealed in Figure 11. Moreover, it is clear that a higher number of
grains with larger maximum Schmid factors along a specific direction ($x$, $y$
or $z$) tends to have larger plastic deformation along that direction.
However, there are grains with large maximum Schmid factors that show little
plastic deformation. Many of this kind of points are dark blue in Figure
8(a)-(c), which means they are near the substrate. Because the substrate is
fixed during simulation, the plastic deformation near the substrate is
limited. Nevertheless, there are still some points at intermediate depth that
do not have much plastic deformation, this is quite surprising and reflects
the complexity of the thermal load in different regions. This behaviour will
be clarified in the following, when the plastic deformation during laser scan
is shown near the melting pool. There are also some grains with larger Schmid
factor that show little plastic deformation and are located near the top
stress-free surface (dark red points), as shown in Figure 8(a)-(c).
To investigate the directionality of the plastic deformation, the range of the
plastic deformation components is investigated, which shows a greater
correlation with the maximum Schmid factor [80], as shown in Figure 9(a)-(c).
The range is calculated as the difference between maximum and minimum values
of the plastic strain components for a given maximum Schmid factor. This
correlation does not depend strongly on the depth $z$. The majority of the
grains undergo plastic compression along the $x$ and $y$ axes, while plastic
tension is observed along the $z$ axis. This is consistent with the average
plastic strain in Figure 6.

###  4.5 Plastic deformation and residual stress

The residual stress, after laser scan and cooling, is due to the heterogeneity
of the plastic strain in different grains. In order to maintain strain
compatibility, the elastic deformation must compensate the difference in
plastic deformation, therefore residual stresses are generated, as modelled by
equations (1) and (11). The correlation between residual stress and plastic
deformation can be visualised using the scatter plots in Figure 10. Each point
represents an element in the representative volume, selected every 16 $\mu$m,
with its corresponding residual stress and plastic strain, averaged over the
integration points. The colour of each point represents the $z$ coordinate of
the corresponding element.
It is certainly true that there is a correlation between residual stress and
plastic strain. Compressive plastic strain along the $x$ and $y$ axes
corresponds to tensile residual stress components along those directions,
while tensile plastic strain along the $z$ axis corresponds to compressive
residual stress along $z$. However, few points show the opposite behaviour:
for instance, a tensile residual stress $\bm{\sigma}_{yy}$ is associated with
a small tensile plastic strain along the $y$ axis in some regions. Overall,
the tensile residual stress components $\bm{\sigma}_{xx}$ and
$\bm{\sigma}_{yy}$ are similar and large, therefore these pre-tensions can
induce fracture along both the $x$ and $y$ directions.

(a) (a)

(b) (b)

(c) (c)

Figure 10:  Correlation between the components of the residual stress and
components of the plastic deformation gradient.

The plastic deformation depends also on the depth. In Figure 10, it is evident
that the regions at the top of the representative volume, represented by dark
red points, are not the ones accommodating the largest plastic strain. The
points with the largest plastic deformation are the ones with depth in the
interval 30 $\mu$m $<z<80$ $\mu$m. These are the points that are closer to the
bottom of the melting pool during laser scan. However, the points at the top
surface can have residual stresses that are similar to points that undergo
larger plastic deformations. This shows that the influence of the plastic
deformation in the substrate is very important; it will be discussed in more
details in section 5. It is important to note that some elements show values
of the residual stress components that are larger than the critical resolved
shear stress to induce plastic deformation. This is due to large volumetric
stress in those elements. Moreover, some dark red points in Figure 10(c) have
residual stress $\sigma_{zz}$ that is different from zero. These points are
picked just below the surface, in the coordinate interval 90 $\mu$m $<z<100$
$\mu$m. Indeed the stress component $\sigma_{zz}$ on the free top surface is
very close to zero, as shown in Figure 12(c).

(a) (a)

(b) (b)

(c) (c)

Figure 11:  Components of the plastic deformation gradient $\bm{F}_{p}$ at the
bottom of the melting pool at time $t=179$ $\mu$s.

(a) (a)

(b) (b)

(c) (c)

Figure 12:  Residual stress components after the cooling stage of the
simulation.

In order to understand the plastic deformation in the representative volume
induced by the laser scan, the components of $\bm{F}_{p}$ are shown in Figure
11(a)-(c) when the laser beam is approximately in the centre of the
representative volume. The thermal expansion ahead of the laser beam leads to
plastic compression along the $x$ axis. This is visible on the two sides and
at the bottom of the melting pool in Figure 11(a), where blue areas are
present. By contrast, the plastic deformation along the $y$ axis is
compressive at the bottom of the melting pool but it is tensile on the two
sides, as shown by the red regions in Figure 11(b). Since plastic deformation
is isochoric, the region at the bottom of the melting pool expands along the
$z$ axis, as shown in Figure 11(c). Because of the subsequent solidification,
the melting pool becomes more shallow and the plastic deformation
concentration shown in Figure 11(a)-(c) spreads in the central region of the
representative volume. This is the reason why the regions with the largest
plastic deformation are concentrated at depth 30 $\mu$m $<z<80$ $\mu$m, as
shown in Figure 10. At the top free surface, the plastic deformation is
slightly smaller than in the centre of the representative volume, as shown in
Figure 8(a)-(c). This is the reason why not all grains with larger Schmid
factor at the top surface can accommodate large plastic deformations, as
stated in section 4.4.
The diagonal components of the residual Cauchy stress tensor are shown in
Figures 12(a)-(c). $\bm{\sigma}_{xx}$ and $\bm{\sigma}_{yy}$ are heterogeneous
at the top surface and in the centre of the representative volume, as shown in
Figures 12(a) and (b). This shows that they are strongly affected by the grain
orientation and that the substrate is applying an important constraint along
$x$ and $y$ on the top surface. In fact, even if the top surface is free,
large values of $\bm{\sigma}_{xx}$ and $\bm{\sigma}_{yy}$ are present. The
compressive stress $\bm{\sigma}_{zz}$ is concentrated in the substrate but it
is not very large at the top surface because of the free boundary condition.

##  5 Discussion

The simulations suggest the following mechanism for the build up of the
residual stresses, as shown in Figure 13. The thermal expansion around the
molten pool induces plastic compression along the horizontal directions, $x$
and $y$, and plastic tension along the out-of-plane $z$ direction. This is
particularly obvious at the bottom of the melting pool, as shown in Figure 11.
After the laser scan is over, tensile residual stresses build up along $x$ and
$y$ directions to compensate the plastic compression. This effect is more
relevant in the centre of the representative volume, however constraints are
present between the regions close to the top surface and the centre. This
induces residual stresses along the $x$ and $y$ directions also close to the
top surface, as shown in Figure 12.

Figure 13:  Mechanism for the build up of the residual stresses.

The different grain orientations cause the heterogeneity of the residual
stress. Grains with slip systems favourably oriented for plastic slip during
loads along the $x$, $y$ and $z$ directions are more likely to develop plastic
deformation along the corresponding directions. The simulations show that the
maximum Schmid factor along the three axes is a suitable quantity to predict
the corresponding plastic deformation and consequent residual stress, as shown
in Figures 8 and 10.
The simulations show that, in the stiffness degradation method, the presence
of regions with a small residual stiffness imposes a constraint on the regions
near the melting pool, preventing more plastic deformation to take place. This
is the origin of the underestimation of the plastic deformation and of the
residual stress using the stiffness degradation method. By contrast, the
element elimination and reactivation method predicts higher plastic
deformation in selective laser melting simulations. This difference between
the two methods must be taken into account in future studies.
Regions with larger tensile residual stresses along specific directions are
more likely to undergo fracture perpendicular to those directions. Therefore,
the present simulations suggest that fracture perpendicular to the $x$ and $y$
axes is more probable because of the magnitude of the corresponding residual
stresses $\bm{\sigma}_{xx}$ and $\bm{\sigma}_{yy}$. These two directions have
similar residual stresses, as shown in Figure 10. This is consistent with
laser-induced microcracking [81], in which fracture is observed at the top
surface perpendicular to the laser scan direction.

##  6 Conclusions

A modelling approach to predict residual stresses at the grain scale in
additive manufactured 316 stainless steel is developed. Thermal-fluid flow and
grain growth simulations are used to provide the temperature field and the
grain structure respectively, which are implemented in crystal plasticity
finite element simulations. A method for element elimination and reactivation
is developed to simulate melting and solidification during the laser scan.
This method allows to reinitialise state variables, such as the plastic
deformation, when elements are reactivated and it represents a step forward
compared with previous studies.
Simulations show that residual stresses are strongly correlated with the
plastic deformation, which in turn depends on the grain orientation. The
directions with largest residual stresses are on the laser scan plane, both
parallel and perpendicular to the laser scan direction. This is consistent
with the microcracks observed experimentally. Even though this study is
focused on 316 stainless steel and selective laser melting, the computational
method developed can be applied to other metals and other additive
manufacturing processes like directed energy deposition and electron beam
additive manufacturing. More detailed simulations of residual stresses will be
conducted in the future with experimental validation to systematically
investigate the relationships between the manufacturing parameters, molten
pool flow, grain structures, and residual stresses.

## Acknowledgement

Thanks to Dr. Yin Zhang from Georgia Institute of Technology for providing
EBSD data used to obtain the Euler angles for grain orientations. Thanks to
Prof. Yinmin Morris Wang for useful discussion on his experiments on AM 316L
stainless steel. This research is supported by the Ministry of Education,
Singapore, under its Academic Research Fund Tier 2 (MOE-T2EP50120-0012).

## References

  * Smith et al. [2016] J. Smith, W. Xiong, W. Yan, S. Lin, P. Cheng, O. L. Kafka, G. J. Wagner, J. Cao, W. K. Liu,  Linking process, structure, property, and performance for metal-based additive manufacturing: computational approaches with experimental support,  Computational Mechanics 57 (2016) 583–610.
  * Schwerdtfeger and Körner [2014] J. Schwerdtfeger, C. Körner,  Selective electron beam melting of Ti–48Al–2Nb–2Cr: Microstructure and aluminium loss,  Intermetallics 49 (2014) 29–35.
  * Kruth et al. [2005] J.-P. Kruth, P. Mercelis, J. Van Vaerenbergh, L. Froyen, M. Rombouts,  Binding mechanisms in selective laser sintering and selective laser melting,  Rapid prototyping journal (2005).
  * Paul et al. [2014] R. Paul, S. Anand, F. Gerner,  Effect of thermal deformation on part errors in metal powder based additive manufacturing processes,  Journal of manufacturing science and Engineering 136 (2014).
  * Li et al. [2018] C. Li, Z. Liu, X. Fang, Y. Guo,  Residual stress in metal additive manufacturing,  Procedia Cirp 71 (2018) 348–353.
  * Taheri et al. [2017] H. Taheri, M. R. B. M. Shoaib, L. W. Koester, T. A. Bigelow, P. C. Collins, L. J. Bond,  Powder-based additive manufacturing-a review of types of defects, generation mechanisms, detection, property evaluation and metrology,  International Journal of Additive and Subtractive Materials Manufacturing 1 (2017) 172–209.
  * Xie et al. [2020] C. Xie, M. Xiang, J. Du, W. Yan, L. Shen, J. Chen, X. Shu, Q. Fang,  Asymmetric yield effect evolving with internal variables during continuous large deformations and its fem validation,  International Journal of Plasticity 130 (2020) 102723.
  * Mukherjee et al. [2017] T. Mukherjee, W. Zhang, T. DebRoy,  An improved prediction of residual stresses and distortion in additive manufacturing,  Computational Materials Science 126 (2017) 360–372.
  * Mori [2006] K. Mori,  Finite element simulation of powder forming and sintering,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 6737–6749. Computational Metal Forming.
  * Grilli et al. [2021a] N. Grilli, E. Tarleton, A. C. Cocks,  Coupling a discrete twin model with cohesive elements to understand twin-induced fracture,  International Journal of Fracture 227 (2021a) 173 – 192.
  * Grilli et al. [2021b] N. Grilli, E. Tarleton, A. C. Cocks,  Neper2CAE and PyCiGen: Scripts to generate polycrystals and interface elements in abaqus,  SoftwareX 13 (2021b) 100651.
  * Grilli et al. [2020] N. Grilli, A. C. Cocks, E. Tarleton,  Crystal plasticity finite element modelling of coarse-grained $\alpha$-uranium,  Computational Materials Science 171 (2020) 109276.
  * Grilli et al. [2019] N. Grilli, A. Cocks, E. Tarleton,  Crystal plasticity finite element simulations of cast $\alpha$-uranium,  in: E. Onate, D. Owen, D. Peric, M. Chiumenti (Eds.), Computational plasticity XV: fundamentals and applications. 15th International Conference on Computational Plasticity - Fundamentals and Applications (COMPLAS), Barcelona, Spain, Sep 03-05, 2019.
  * Hocine et al. [2020] S. Hocine, H. Van Swygenhoven, S. Van Petegem, C. S. T. Chang, T. Maimaitiyili, G. Tinti, D. Ferreira Sanchez, D. Grolimund, N. Casati,  Operando x-ray diffraction during laser 3d printing,  Materials Today 34 (2020) 30–40.
  * Vrancken et al. [2014] B. Vrancken, V. Cain, R. Knutsen, J. Van Humbeeck,  Residual stress via the contour method in compact tension specimens produced via selective laser melting,  Scripta Materialia 87 (2014) 29–32.
  * Mathar et al. [1934] J. Mathar, et al.,  Determination of initial stresses by measuring the deformation around drilled holes,  Trans. ASME 56 (1934) 249–254.
  * Martínez-García et al. [2019] V. Martínez-García, G. Pedrini, P. Weidmann, A. Killinger, R. Gadow, W. Osten, S. Schmauder,  Non-contact residual stress analysis method with displacement measurements in the nanometric range by laser made material removal and slm based beam conditioning on ceramic coatings,  Surface and Coatings Technology 371 (2019) 14–19.
  * Ganeriwala et al. [2019] R. Ganeriwala, M. Strantza, W. King, B. Clausen, T. Q. Phan, L. E. Levine, D. W. Brown, N. Hodge,  Evaluation of a thermomechanical model for prediction of residual stress during laser powder bed fusion of Ti-6Al-4V,  Additive Manufacturing 27 (2019) 489–502.
  * Liang et al. [2018] X. Liang, L. Cheng, Q. Chen, Q. Yang, A. C. To,  A modified method for estimating inherent strains from detailed process simulation for fast residual distortion prediction of single-walled structures fabricated by directed energy deposition,  Additive Manufacturing 23 (2018) 471–486.
  * Schoinochoritis et al. [2017] B. Schoinochoritis, D. Chantzis, K. Salonitis,  Simulation of metallic powder bed additive manufacturing processes with the finite element method: A critical review,  Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 231 (2017) 96–117.
  * Qian et al. [2014] G. Qian, V. González-Albuixech, M. Niffenegger,  Probabilistic assessment of a reactor pressure vessel subjected to pressurized thermal shocks by using crack distributions,  Nuclear Engineering and Design 270 (2014) 312–324.
  * Yan et al. [2020] W. Yan, Y. Lu, K. Jones, Z. Yang, J. Fox, P. Witherell, G. Wagner, W. K. Liu,  Data-driven characterization of thermal models for powder-bed-fusion additive manufacturing,  Additive Manufacturing 36 (2020) 101503.
  * Bailey et al. [2017] N. S. Bailey, C. Katinas, Y. C. Shin,  Laser direct deposition of aisi h13 tool steel powder with numerical modeling of solid phase transformation, hardness, and residual stresses,  Journal of Materials Processing Technology 247 (2017) 223–233.
  * Cheon et al. [2016] J. Cheon, D. V. Kiran, S.-J. Na,  Thermal metallurgical analysis of GMA welded AH36 steel using CFD–FEM framework,  Materials & Design 91 (2016) 230–241.
  * Chen and Yan [2020] F. Chen, W. Yan,  High-fidelity modelling of thermal stress for additive manufacturing by linking thermal-fluid and mechanical models,  Materials & Design 196 (2020) 109185.
  * Yan et al. [2018] W. Yan, Y. Qian, W. Ge, S. Lin, W. K. Liu, F. Lin, G. J. Wagner,  Meso-scale modeling of multiple-layer fabrication process in selective electron beam melting: Inter-layer/track voids formation,  Materials & Design 141 (2018) 210–219.
  * Wei et al. [2021] H. Wei, T. Mukherjee, W. Zhang, J. Zuback, G. Knapp, A. De, T. DebRoy,  Mechanistic models for additive manufacturing of metallic components,  Progress in Materials Science 116 (2021) 100703.
  * Stender et al. [2018] M. E. Stender, L. L. Beghini, J. D. Sugar, M. G. Veilleux, S. R. Subia, T. R. Smith, C. W. San Marchi, A. A. Brown, D. J. Dagel,  A thermal-mechanical finite element workflow for directed energy deposition additive manufacturing process modeling,  Additive Manufacturing 21 (2018) 556–566.
  * Montevecchi et al. [2016] F. Montevecchi, G. Venturini, A. Scippa, G. Campatelli,  Finite element modelling of wire-arc-additive-manufacturing process,  Procedia Cirp 55 (2016) 109–114.
  * Yang et al. [2016] Q. Yang, P. Zhang, L. Cheng, Z. Min, M. Chyu, A. C. To,  Finite element modeling and validation of thermomechanical behavior of Ti-6Al-4V in directed energy deposition additive manufacturing,  Additive Manufacturing 12 (2016) 169–177.
  * Lindgren et al. [1999] L.-E. Lindgren, H. Runnemalm, M. O. Näsström,  Simulation of multipass welding of a thick plate,  International journal for numerical methods in engineering 44 (1999) 1301–1316.
  * Lindgren and Hedblom [2001] L.-E. Lindgren, E. Hedblom,  Modelling of addition of filler material in large deformation analysis of multipass welding,  Communications in numerical methods in engineering 17 (2001) 647–657.
  * Michaleris [2014] P. Michaleris,  Modeling metal deposition in heat transfer analyses of additive manufacturing processes,  Finite Elements in Analysis and Design 86 (2014) 51–60.
  * Ales [2018] T. K. Ales, An integrated model for the probabilistic prediction of yield strength in electron-beam additively manufactured Ti-6Al-4V, Ph.D. thesis, Iowa State University, 2018\.
  * Miehe et al. [2010] C. Miehe, F. Welschinger, M. Hofacker,  Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations,  International Journal for Numerical Methods in Engineering 83 (2010) 1273–1311.
  * Borden et al. [2018] M. Borden, T. Hughes, C. Landis, A. Anvari, I. Lee, Phase-Field Formulation for Ductile Fracture, In: Oñate E., Peric D., de Souza Neto E., Chiumenti M. (eds) Advances in Computational Plasticity. Computational Methods in Applied Sciences, vol 46. Springer, Cham., 2018.
  * Grilli et al. [2018] N. Grilli, C. A. Duarte, M. Koslowski,  Dynamic fracture and hot-spot modeling in energetic composites,  Journal of Applied Physics 123 (2018) 065101.
  * Duarte et al. [2018] C. A. Duarte, N. Grilli, M. Koslowski,  Effect of initial damage variability on hot-spot nucleation in energetic materials,  Journal of Applied Physics 124 (2018) 025104.
  * Grilli and Koslowski [2019] N. Grilli, M. Koslowski,  The effect of crystal anisotropy and plastic response on the dynamic fracture of energetic materials,  Journal of Applied Physics 126 (2019) 155101.
  * Chen et al. [2019] W. Chen, T. Voisin, Y. Zhang, J.-B. Florien, C. M. Spadaccini, D. L. McDowell, T. Zhu, Y. M. Wang,  Microscale residual stresses in additively manufactured stainless steel,  Nature communications 10 (2019) 1–12.
  * Wang et al. [2020] G. Wang, H. Ouyang, C. Fan, Q. Guo, Z. Li, W. Yan, Z. Li,  The origin of high-density dislocations in additively manufactured metals,  Materials Research Letters 8 (2020) 283–290.
  * Grilli et al. [2015] N. Grilli, K. G. Janssens, H. Van Swygenhoven,  Crystal plasticity finite element modelling of low cycle fatigue in fcc metals,  Journal of the Mechanics and Physics of Solids 84 (2015) 424–435.
  * Grilli et al. [2020] N. Grilli, P. Earp, A. C. Cocks, J. Marrow, E. Tarleton,  Characterisation of slip and twin activity using digital image correlation and crystal plasticity finite element simulation: Application to orthorhombic $\alpha$-uranium,  Journal of the Mechanics and Physics of Solids 135 (2020) 103800.
  * Yang et al. [2021] M. Yang, L. Wang, W. Yan,  Phase-field modeling of grain evolutions in additive manufacturing from nucleation, growth, to coarsening,  npj Computational Materials 7 (2021) 56.
  * Yan et al. [2018] W. Yan, Y. Lian, C. Yu, O. L. Kafka, Z. Liu, W. K. Liu, G. J. Wagner,  An integrated process–structure–property modeling framework for additive manufacturing,  Computer Methods in Applied Mechanics and Engineering 339 (2018) 184–204.
  * Asaro and Rice [1977] R. J. Asaro, J. Rice,  Strain localization in ductile single crystals,  Journal of the Mechanics and Physics of Solids 25 (1977) 309–338.
  * Kalidindi [1998] S. R. Kalidindi,  Incorporation of deformation twinning in crystal plasticity models,  Journal of the Mechanics and Physics of Solids 46 (1998) 267–290.
  * Vujošević and Lubarda [2002] L. Vujošević, V. Lubarda,  Finite-strain thermoelasticity based on multiplicative decomposition of deformation gradient,  Theoretical and applied mechanics (2002) 379–399.
  * Grilli et al. [2020] N. Grilli, A. C. Cocks, E. Tarleton,  A phase field model for the growth and characteristic thickness of deformation-induced twins,  Journal of the Mechanics and Physics of Solids 143 (2020) 104061.
  * Roters et al. [2010] F. Roters, P. Eisenlohr, L. Hantcherli, D. D. Tjahjanto, T. R. Bieler, D. Raabe,  Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications,  Acta Materialia 58 (2010) 1152–1211.
  * Roters et al. [2019] F. Roters, M. Diehl, P. Shanthraj, P. Eisenlohr, C. Reuber, S. Wong, T. Maiti, A. Ebrahimi, T. Hochrainer, H.-O. Fabritius, S. Nikolov, M. Friák, N. Fujita, N. Grilli, K. Janssens, N. Jia, P. Kok, D. Ma, F. Meier, E. Werner, M. Stricker, D. Weygand, D. Raabe,  DAMASK – The Düsseldorf Advanced Material Simulation Kit for modeling multi-physics crystal plasticity, thermal, and damage phenomena from the single crystal up to the component scale,  Computational Materials Science 158 (2019) 420–478.
  * Rice [1971] J. R. Rice,  Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity,  Journal of the Mechanics and Physics of Solids 19 (1971) 433–455.
  * Peirce et al. [1983] D. Peirce, R. J. Asaro, A. Needleman,  Material rate dependence and localized deformation in crystalline solids,  Acta metallurgica 31 (1983) 1951–1976.
  * Daymond and Bouchard [2006] M. Daymond, P. Bouchard,  Elastoplastic deformation of 316 stainless steel under tensile loading at elevated temperatures,  Metallurgical and Materials transactions A 37 (2006) 1863–1873.
  * Grilli et al. [2020] N. Grilli, E. Tarleton, P. D. Edmondson, M. N. Gussev, A. C. F. Cocks,  In situ measurement and modelling of the growth and length scale of twins in $\alpha$-uranium,  Phys. Rev. Materials 4 (2020) 043605.
  * Kalidindi and Anand [1992] S. R. Kalidindi, L. Anand,  An approximate procedure for predicting the evolution of crystallographic texture in bulk deformation processing of fcc metals,  International journal of mechanical sciences 34 (1992) 309–329.
  * Diehl et al. [2017] M. Diehl, M. Wicke, P. Shanthraj, F. Roters, A. Brueckner-Foit, D. Raabe,  Coupled crystal plasticity–phase field fracture simulation study on damage evolution around a void: Pore shape versus crystallographic orientation,  JOM 69 (2017) 872–878.
  * Irastorza-Landa et al. [2016] A. Irastorza-Landa, H. Van Swygenhoven, S. Van Petegem, N. Grilli, A. Bollhalder, S. Brandstetter, D. Grolimund,  Following dislocation patterning during fatigue,  Acta Materialia 112 (2016) 184–193.
  * Callen [1998] H. B. Callen, Thermodynamics and an introduction to thermostatistics, 1998.
  * Dubrovinsky [2002] L. Dubrovinsky,  Thermal Expansion and Equation of State,  Encyclopedia of Materials: Science and Technology (2002) 1–4.
  * Yadroitsev and Yadroitsava [2015] I. Yadroitsev, I. Yadroitsava,  Evaluation of residual stress in stainless steel 316l and ti6al4v samples produced by selective laser melting,  Virtual and Physical Prototyping 10 (2015) 67–76.
  * Grimvall [1999] G. Grimvall, Thermophysical properties of materials, Elsevier, 1999.
  * Jiang et al. [2012] W. Jiang, Y. Zhang, W. Woo,  Using heat sink technology to decrease residual stress in 316l stainless steel welding joint: Finite element simulation,  International journal of pressure vessels and piping 92 (2012) 56–62.
  * Grilli and Koslowski [2018] N. Grilli, M. Koslowski,  The effect of crystal orientation on shock loading of single crystal energetic materials,  Computational Materials Science 155 (2018) 235–245.
  * Reed and Horiuchi [1982] R. Reed, T. Horiuchi, Austenitic steels at low temperatures, Plenum press, New York, 1982\.
  * Grilli et al. [2018] N. Grilli, K. Janssens, J. Nellessen, S. Sandlöbes, D. Raabe,  Multiple slip dislocation patterning in a dislocation-based crystal plasticity finite element method,  International Journal of Plasticity 100 (2018) 104–121.
  * Yan et al. [2018] W. Yan, S. Lin, O. L. Kafka, C. Yu, Z. Liu, Y. Lian, S. Wolff, J. Cao, G. J. Wagner, W. K. Liu,  Modeling process-structure-property relationships for additive manufacturing,  Frontiers of Mechanical Engineering 13 (2018) 482–492.
  * Permann et al. [2020] C. J. Permann, D. R. Gaston, D. Andrš, R. W. Carlsen, F. Kong, A. D. Lindsay, J. M. Miller, J. W. Peterson, A. E. Slaughter, R. H. Stogner, R. C. Martineau,  MOOSE: Enabling massively parallel multiphysics simulation,  SoftwareX 11 (2020) 100430.
  * Adhikary et al. [2016] D. P. Adhikary, C. Jayasundara, R. K. Podgorney, A. H. Wilkins,  A robust return-map algorithm for general multisurface plasticity,  International Journal for Numerical Methods in Engineering 109 (2016) 218–234.
  * Chockalingam et al. [2013] K. Chockalingam, M. Tonks, J. Hales, D. Gaston, P. Millett, L. Zhang,  Crystal plasticity with jacobian-free newton–krylov,  Computational Mechanics 51 (2013) 617–627.
  * Lee [1969] E. H. Lee,  Elastic-plastic deformation at finite strains (1969).
  * Irastorza-Landa et al. [2017] A. Irastorza-Landa, N. Grilli, H. Van Swygenhoven,  Effect of pre-existing immobile dislocations on the evolution of geometrically necessary dislocations during fatigue,  Modelling and Simulation in Materials Science and Engineering 25 (2017) 055010\.
  * Clausen et al. [1998] B. Clausen, T. Lorentzen, T. Leffers,  Self-consistent modelling of the plastic deformation of F.C.C. polycrstals and its implications for diffraction measurements of internal stresses,  Acta Materialia 46 (1998) 3087–3098.
  * Yushu et al. [2021] D. Yushu, W. Jiang, C. Tan, C. Sun, D. Schwen, B. Spencer,  A thermal-mechanical model for directed energy deposition process with element activation,  (submitted manuscript) (2021).
  * Duarte et al. [2019] C. A. Duarte, M. Koslowski, N. Grilli,  Modeling the effect of plasticity and damage in $\beta$-hmx single crystals under shock loading,  21st Biennial Conference of the APS Topical Group on Shock Compression of Condensed Matter 64 (2019).
  * Yan et al. [2017] W. Yan, W. Ge, Y. Qian, S. Lin, B. Zhou, W. K. Liu, F. Lin, G. J. Wagner,  Multi-physics modeling of single/multiple-track defect mechanisms in electron beam selective melting,  Acta Materialia 134 (2017) 324–333.
  * Wang et al. [2020] L. Wang, Y. Zhang, W. Yan,  Evaporation model for keyhole dynamics during additive manufacturing of metal,  Physical Review Applied 14 (2020) 064039.
  * Simmons et al. [2000] J. Simmons, C. Shen, Y. Wang,  Phase field modeling of simultaneous nucleation and growth by explicitly incorporating nucleation events,  Scripta materialia 43 (2000) 935–942.
  * Grilli [2016] N. Grilli, Physics-based constitutive modelling for crystal plasticity finite element computation of cyclic plasticity in fatigue, Ph.D. thesis, École Polytechnique Fédérale de Lausanne, 2016.
  * Irastorza-Landa et al. [2017] A. Irastorza-Landa, N. Grilli, H. Van Swygenhoven,  Laue micro-diffraction and crystal plasticity finite element simulations to reveal a vein structure in fatigued Cu,  Journal of the Mechanics and Physics of Solids 104 (2017) 157–171.
  * Vrancken et al. [2020] B. Vrancken, R. K. Ganeriwala, M. J. Matthews,  Analysis of laser-induced microcracking in tungsten under additive manufacturing conditions: Experiment and simulation,  Acta Materialia 194 (2020) 464–472.