# Elasto-plastic residual stress analysis of selective sintered porous
materials based on 3D-multilayer thermo-structural phase-field simulations

Yangyiwei Yang  Somnath Bharech  Nick Finger  Xiandong Zhou<EMAIL_ADDRESS>Jörg Schröder  Bai-Xiang Xu<EMAIL_ADDRESS>Mechanics of Functional
Materials Division, Institute of Materials Science, Technische Universität
Darmstadt, Darmstadt 64287, Germany  Failure Mechanics and Engineering
Disaster Prevention Key Laboratory of Sichuan Province, College of
Architecture and Environment, Sichuan University, Chengdu 610207, China
Fakultät für Ingenieurwissenschaften, Abteilung Bauwissenschaften, Institut
für Mechanik, Universität Duisburg-Essen, Essen 45141, Germany

###### Abstract

Residual stress and plastic strain in additive manufacturing materials can
exhibit significant microscopic variation at the powder scale, profoundly
influencing the overall properties of printed components. This variation
depends on process parameters and stems from multiple factors, including
differences in powder bed morphology, non-uniform thermo-structural profiles,
and inter-layer fusion. In this research, we propose a powder-resolved
multilayer multi-physics simulation scheme tailored for porous materials
through the process of selective sintering. This approach seamlessly
integrates finite element method (FEM)-based non-isothermal phase-field
simulation with thermo-elasto-plastic simulation, incorporating temperature-
and phase-dependent material properties. The outcome of this investigation
includes a detailed depiction of the mesoscopic evolution of stress and
plastic strain within a transient thermal-microstructure, evaluated across a
spectrum of beam power and scan speed parameters. Simulation results further
reveal the underlying mechanisms. For instance, stress concentration primarily
occurs at the necking region of partially melted particles and the junctions
between different layers, resulting in the accumulation of plastic strain and
residual stress, ultimately leading to structural distortion in the materials.
Based on the simulation data, regression models were proposed to describe the
dependency of the residual stress and the plastic strain on densification or
beam energy input, along with the comparison to experimental results.

###### keywords:

additive manufacturing, powder bed fusion, heat transfer, residual stress,
microstructure

##  1 Introduction

In the recent decade, additive manufacturing (AM) emerged from a niche
technology to a widely applicable means of production. Among a variety of AM
techniques, powder bed fusion (PBF) stands out as a strong technique for
producing intricate geometries while maintaining good structural integrity and
superior material properties [1, 2, 3, 4, 5]. In PBF, layers of powder get
fused by a beam one after another to build up three-dimensional (3D)
geometries. Compared to conventional ones, this manufacturing process offers
several advantages, such as great design flexibility, reduction in waste,
enabling rapid tooling and decreasing production cycles.

Selective sintering (SS) is one prominent example of PBF, which utilizes the
tuned incident beam (mostly continuous laser scan or laser pulse) to bind the
free powders layerwisely. Due to its well-controlled powder bed temperature
compared to other PBF techniques like selective melting, where the significant
melting phenomenon can be captured, SS enables the partial melting of
particles and produces samples with relatively high porosity [6, 7, 8, 9, 10].
In this regard, SS has great potential in applications that require complex
geometries with high porosity. For instance, SS is applicable in manufacturing
porous components for medication applications, especially medical scaffolds
and artificial bones [11, 12, 13]. It is feasible in manufacuring functional
materials with a low processing temperature, such as ferromagnetic materials
[14, 15, 16, 17]. By precisely controlling the geometry and structure of the
printed material, it may also be possible to create the desired strain
gradients and electric polarization necessary to generate the flexoelectric
effect [18, 19, 20].

Residual stress has been a critical issue since it can affect the mechanical
properties, dimensional accuracy, corrosion resistance, crack growth
resistance, and performance of AM parts [21]. There are already some
researches investigating the residual stress distribution across the parts
developed using AM techniques, using experimental measurements and theoretical
estimations [22, 23, 24, 25]. Mercelis et al. explained the origins of
residual stress in parts produced by AM. They experimentally measured the
residual stress distribution in a selective laser sintering produced part and
compared it with the analytical and numerical solutions [22]. Pant et al.
employed a layer-by-layer finite element approach to predict the residual
stress and validated the model using measurements from neutron diffraction
[24]. Ibraheem et al. predicted the thermal profile and residual stress of SLS
processed H13 tool steel using a finite element model. High porosity was
observed in the range 25-40%, yet a homogenized powder bed was utilized to
simulate the thermal evolution and subsequently the residual stress [26].

The calculation of residual stress in literature is often performed in two
sequential steps: first the transient temperature is simulated in the entire
domain and then the thermal history is used to calculate the stress and strain
evolution in thermomechanical simulations. The accuracy of these calculations
critically rely on the quality of the thermal history. The temperature
gradient mechanism (TGM) model and the cool-down stage (CDS) model are two
commonly used models to explain the development of plastic strain and residual
stress formation mechanisms [22, 27, 28]. In TGM model, the deformation of the
material in the molten pool/fusion zone is restrained by surrounding materials
during both the heating and cooling stages due to a large temperature gradient
inside and outside the overheated region (where the on-site temperature is
beyond melting point). In the heating stage, the plastic strain is developed
in the surrounding material due to an expansion of the heated material. As the
heated material cools down, it starts to shrink but is again contracted by the
plastically deformed surrounding material. Thus a tensile residual stress is
formed. The CDS model, on the other hand, was proposed to elucidate the
residual stress due to the layerwise features of the AM process and the AM-
manufactured parts. This occurs because, in both the heating and cooling
stages, the deformation of the upper-layers is restrained by the fused lower-
layers and/or substrate. Meanwhile, melted material in the lower-layers will
undergo a remelting and re-solidification cycle. Both factors result in
tensile residual stress in the top layer due to the shrinkage of the material
during the cooling stage. These two models provide phenomenological aspects by
employing the idealized homogeneous layers for the understanding of the
formation mechanism of residual stresses. They have been widely adopted in
numerical analysis of residual stress at single layer scale models [29, 30,
31, 32] and part-scale multilayer models [33, 34, 35].

In these aforementioned works, homogeneous layers are used for analyses of the
thermal history and the mechanical response in an AM-manufactured part.
However, due to the complex morphology of the powder bed on the powder scale
and the resultant inhomogeneity of the local thermal properties, the thermal
profile is subjected to a high level of non-uniformity as well. This implies
heterogeneity of the thermal stress, the plastic deformation and the residual
stress on the powder scale. In fact, mechanical properties of AM-built parts,
such as fracture strength and hardening behavior, are significantly influenced
by the local defects and local stress. The diversity in the powder bed
structure and the packing density introduces thermal heterogeneity in the form
of mesoscopic high-gradient temperature profile and asynchronous on-site
thermal history. Such variability leads to varying degrees of thermal
expansion and, therefore, thermal stresses. Additionally, the presence of
stochastic inter-particle voids and lack-of-fusion pores contribute to
evolving disparities in thermo-mechanical properties [7, 36, 37]. Thermal
gradient across the powder particles is very sharp as per the modeling results
by Panwisawas et al. [38], which is experimentally observed in [39] and also
recaptured in our former numerical works [7, 36, 37]. Furthermore, the
layerwisely build-up process also results in local interaction between the
newly deposited and previously deposited layers with high surface roughness.
Thus, the residual stress formation on the powder scale should be revisited
for the SS process.

For this purpose, the powder bed morphology, the heat transfer and the pore
structure evolution during the SS process should be primarily considered. One
promising approach is the phase-field simulation. For instance, phase-field
simulations of the AM process can reveal the in-process microstructure
evolution, providing insights into key features such as porosity, surface
morphology, temperature profile, and geometry evolution [40, 41]. Based on a
non-isothermal phase-field multilayer simulation results of selective sintered
porous structures, thermo-elastic simulations and homogenization of the
elastic properties have been investigated in our previous work [36]. In
combination with the single layer phase-field thermo-microstructural
simulations, the first powder-resolved elasto-plastic simulations of the
residual stress in a SS-processed magnetic alloy has also been performed [37].
Results show that during the cooling stage, the partially melted and inter-
connected particles result in plastic deformation due to the shrinkage of the
fusion zone, but mostly at the necking region because of stress concentration
in porous microstructures. By combining the thermo-structural phase-field
methods and the thermomechanical calculations, this work demonstrates the
promising capability of the multi-physics simulation scheme for the local
residual stress analysis of AM-processed materials on the powder scale. In the
current article, we extend this multi-physics simulation scheme for elasto-
plastic multilayer SS process, providing systematic simulations of the local
plastic deformation and the residual stress. These allow us to reveal the
influential aspects related to the printing process and the powder bed.

##  2 Theoretical framework

The simulation workflow is arranged in three stages, as shown in the Fig. 1a.
The powders were first deposited into a prime simulation domain under a given
gravitational force based on the discrete element method (DEM), then non-
isothermal phase-field simulations were performed for the coupled thermo-
structural evolution during the SS process. Upon completion of a layer, the
resulting microstructure is voxelized and re-imported into the DEM program for
the deposition of the next powder layer until reaching the desired layer
number or height. Meanwhile, the subsequent thermo-elasto-plastic calculations
were performed to investigate the development of plastic strain and stress
under the quasi-static thermo-structures by mapping the nodal values of the
order parameters (OP, indicating the chronological-spatial distribution of the
phases) and temperature to a subdomain. It is worth noting that the thermo-
elasto-plastic calculation of next layer was restarted from the former one
with continuous history of the residual stress and accumulated plastic strain,
though the OPs and temperature were reiterated in the prime domain. As it is
usually done in the literature, we assume thereby that heat transfer is only
strongly coupled with microstructure evolution (driven by diffusion and
underlying grain growth) but is not impacted by mechanics.

###  2.1 Non-isothermal phase-field model for multilayer SS processes

The non-isothermal phase-field model employed in this work to simulate the
multilayer is based on our former works [7, 42, 36, 37]. For the sake of
completeness, in this section we summarize the essentials of the employed
model.

The model adopts both conserved and non-conserved OPs for representing the
microstructure evolution of polycrystalline material. The conserved OP $\rho$
indicates the substance, including the unfused and partially melted regions,
while the non-conserved OPs $\\{\eta_{i}\\},~{}i=1,2,\dots$ distinguish
particles with different crystallographic orientations. To guarantee that the
orientation fields ($\\{\eta_{i}\\}$) get valued only in the substance
($\rho=1$), a numerical constraint $\sum_{i}\eta_{i}+(1-\rho)=1$ is imposed to
the simulation domain. This constraint also applies to the substrate where
$\rho=1$ always holds.

The thermo-structural evolution is governed by

$\displaystyle\frac{\partial\rho}{\partial
t}=\nabla\cdot\mathbf{M}\cdot\nabla\left[\frac{\partial
f}{\partial\rho}-\underline{\kappa}_{\rho}(T)\nabla^{2}\rho\right],$ (1)
$\displaystyle\frac{\partial\eta_{i}}{\partial t}=-L\left[\frac{\partial
f}{\partial\eta_{i}}-\underline{\kappa}_{\eta}(T)\nabla^{2}\eta_{i}\right],$
(2) $\displaystyle c_{\mathrm{r}}\left[\frac{\partial T}{\partial
t}-\mathbf{v}\cdot\nabla T\right]=\nabla\cdot\mathbf{K}\cdot\nabla
T+q_{\mathbf{v}},$ (3)

where the local free energy density is formulated as

$\begin{split}f(T,\rho,\\{\eta_{i}\\})&=h_{\mathrm{ss}}(\rho)f_{\text{ht}}(T)+\underline{W}_{\mathrm{sf}}(T)\left[\rho^{2}(1-\rho)^{2}\right]+\\\
&\underline{W}_{\mathrm{gb}}(T)\left[\rho^{2}+6(1-\rho)\sum_{i}\eta_{i}^{2}-4(2-\rho)\sum_{i}\eta_{i}^{3}+3\left(\sum_{i}\eta_{i}^{2}\right)^{2}\right],\\\
f_{\mathrm{ht}}(T)&=c_{\mathrm{r}}\left[\left(T-T_{\mathrm{M}}\right)-T\ln\frac{T}{T_{\mathrm{M}}}\right]+f_{\mathrm{ref}}^{T_{\mathrm{M}}}-\frac{T-T_{\mathrm{M}}}{T_{\mathrm{M}}}\mathcal{L}_{\mathrm{M}}h_{\mathrm{M}}(T).\end{split}$
(4)

Here $f(T,\rho,\\{\eta_{i}\\})$ contains multiple local minima, representing
the thermodynamic stability of the states, i.e., substance, atmosphere/pore
and grains with distinct orientations, with the coefficients
$\underline{W}_{\mathrm{sf}}(T)$ and $\underline{W}_{\mathrm{gb}}(T)$ related
to the barrier heights among minima, varying with temperature [7, 37]. The
heat contribution $f_{\text{ht}}$ thereby manifests the stability of states by
shifting the minima according to the local temperature field. When the
material reached melting point $T_{\mathrm{M}}$, the extra contribution due to
the latent heat $\mathcal{L}_{\mathrm{M}}$ is mapped by the interpolation
function $h_{\text{ht}}$, which approaches unity when $T\rightarrow
T_{\text{M}}$ [7, 37]. It is worth noting that
$\underline{W}_{\mathrm{sf}}(T)$, $\underline{\kappa}_{\mathrm{sf}}(T)$,
$\underline{W}_{\mathrm{gb}}(T)$, and $\underline{\kappa}_{\mathrm{gb}}(T)$
inherit their temperature dependencies from the surface energy
$\gamma_{\mathrm{sf}}(T)$ and grain boundary energy $\gamma_{\mathrm{gb}}(T)$,
respectively (see Sec. 3.1). At the equilibrium, they can be related to
$\gamma_{\mathrm{sf}}(T)$, $\gamma_{\mathrm{gb}}(T)$, and a diffuse-interface
width of the grain boundary $l_{\mathrm{gb}}$, which have been explained in
our former work [7, 37].

The Cahn-Hilliard mobility employed in Eq. $($1$)$ equation adopts the
anisotropic form. Contributions considered in this work contain not only the
mass transfer through substance, atmosphere, surface and grain boundary, but
also the diffusion enhancement due to possible partial melting [42, 37, 43],
i.e.

$\begin{split}\mathbf{M}=\frac{1}{2(\underline{W}_{\mathrm{sf}}+\underline{W}_{\mathrm{gb}})}\left[h_{\mathrm{ss}}D_{\mathrm{ss}}\mathbf{I}+h_{\mathrm{at}}D_{\mathrm{at}}\mathbf{I}+h_{\mathrm{sf}}D_{\mathrm{sf}}\mathbf{T}_{\mathrm{sf}}+h_{\mathrm{gb}}D_{\mathrm{gb}}\mathbf{T}_{\mathrm{gb}}\right]+h_{\mathrm{M}}M_{\mathrm{M}}\mathbf{T}_{\mathrm{sf}},\end{split}$
(5)

where $D_{\mathrm{path}}$ is the effective diffusivity of the mass transport
via $\mathrm{path}=\mathrm{ss},~{}\mathrm{at},~{}\mathrm{sf},~{}\mathrm{gb}$,
and $h_{\text{ss}}$, $h_{\text{at}}$, $h_{\text{sf}}$ and $h_{\text{gb}}$ are
also interpolation functions to indicate substance (including solid and
liquid), atmosphere/pore, surface and grain boundary, respectively, which
obtain unity only in the corresponding region [42, 37, 43]. Mass transport
along surface and grain boundary is also regulated by corresponding projection
tensor $\mathbf{T}_{\mathrm{gb}}$ and $\mathbf{T}_{\mathrm{sf}}$. It is worth
noting that the partial melting contribution $M_{\text{ml}}$ is treated as an
enhanced surface diffusion when $T\rightarrow T_{\text{M}}$ due to the
assumption of a limited melting phenomenon around the surface of particles. In
this sense, formulations in Eqs. (1) and (5) disregard the contributions from
melt flow dynamics as well as from inter-coupling effects between mass and
heat transfer, i.e., Soret and Dufour effects [42]. The isotropic Allen-Cahn
mobility is explicitly formulated by the the temperature-dependent GB mobility
$G_{\text{gb}}$, GB energy $\Gamma_{\text{gb}}$ and gradient coefficient
$\underline{\kappa}_{\eta}$ as [44, 45]

$L=\frac{G_{\text{gb}}\Gamma_{\text{gb}}}{\underline{\kappa}_{\eta}}.$ (6)

The phase-dependent thermal conductivity tensor takes the continuity of the
thermal flux along both the normal and tangential directions of the surface
into account, and is formulated as [46, 47]

$\mathbf{K}=\left[h_{\mathrm{ss}}K_{\mathrm{ss}}+h_{\mathrm{at}}K_{\mathrm{at}}\right]\mathbf{N}_{\mathrm{sf}}+\left[\frac{K_{\mathrm{ss}}K_{\mathrm{at}}}{h_{\mathrm{ss}}K_{\mathrm{at}}+h_{\mathrm{at}}K_{\mathrm{ss}}}\right]\mathbf{T}_{\mathrm{sf}},$
(7)

where $K_{\mathrm{ss}}$ and $K_{\mathrm{at}}$ are the thermal conductivity of
the substance and pore/atmosphere. $\mathbf{N}$ is the normal tensor of the
surface [46, 37, 43]. It is worth noting that radiation contribution via
pore/atmosphere is considered in $K_{\mathrm{at}}$ as
$K_{\mathrm{at}}=K_{0}+4T^{3}\sigma_{\mathrm{B}}\ell_{\mathrm{rad}}/3$ with
$K_{0}$ the thermal conductivity of the Argon gas and $\sigma_{\mathrm{B}}$
the Stefan-Boltzmann constant. $\ell_{\mathrm{rad}}$ is the effective
radiation path between particles, which usually takes the average diameter of
the powders [48].

The beam-incited thermal effect is equivalently treated as a volumetric heat
source with its distribution along the depth direction formulated in a
radiation penetration fashion, as the powder bed is regarded as an effective
homogenized optical medium, i.e.

$q(\mathbf{r},t)=Pp_{xy}[\mathbf{r}_{O}(\mathbf{v},t)]\frac{\text{d}a}{\text{d}z},$

in which the in-plane Gaussian distribution $p_{xy}$ with a moving center
$\mathbf{r}_{O}(\mathbf{v},t)$. $P$ is the beam power and $\mathbf{v}$ is the
scan velocity with its magnitude $v=|\mathbf{v}|$ as the scan speed, which are
two major processing parameters in this work. The absorptivity profile
function along depth ${\text{d}a}/{\text{d}z}$ which is calculated based on
Refs. [7, 49]. It is obvious that when one ignores the effects of those latent
heat induced by microstructure evolution (i.e. the evolution of the
pore/substance as well as unique grains), Eq. (3) can be degenerated to the
conventional Fourier’s equation for heat conduction with an internal heat
source.

###  2.2 Elasto-plastic model for thermo-mechanical analysis

The transient temperature field and substance field $\rho$, resulting from the
non-isothermal phase-field simulations, are then imported into the quasi-
static elasto-plastic model as the thermal load as well as the phase indicator
for the interpolation of mechanical properties. The linear momentum equation
reads

$\nabla\cdot\bm{\sigma}=\mathbf{0},$ (8)

where $\bm{\sigma}$ is the $2^{\mathrm{nd}}$-order stress tensor. Taking the
Voigt-Taylor interpolation scheme (VTS), where the total strain is assumed to
be identical among pores/atmosphere and substance, i.e.
$\bm{\varepsilon}=\bm{\varepsilon}_{\mathrm{ss}}=\bm{\varepsilon}_{\mathrm{at}}$.
[50, 51, 52]. In this regard, the stress can be eventually formulated by the
linear constitutive equation

$\bm{\sigma}=\mathbf{C}(\rho,T):\left(\bm{\varepsilon}-\bm{\varepsilon}^{\mathrm{th}}-\bm{\varepsilon}^{\mathrm{pl}}\right)$
(9)

where the $4^{\mathrm{th}}$-order elastic tensor is interpolated from the
substance one $\mathbf{C}_{\mathrm{ss}}$ and the pores/atmosphere one
$\mathbf{C}_{\mathrm{at}}$, i.e.,

$\mathbf{C}(\rho,T)=h_{\mathrm{ss}}(\rho)\mathbf{C}_{\mathrm{ss}}(T)+h_{\mathrm{at}}(\rho)\mathbf{C}_{\mathrm{at}}.$
(10)

In this work, $\mathbf{C}_{\mathrm{ss}}$ is calculated from the temperature-
dependent Youngs’ modulus $E(T)$ and poison ratio $\nu(T)$, while
$\mathbf{C}_{\mathrm{at}}$ is assigned with a sufficiently small value.

The thermal eigenstrain $\bm{\varepsilon}^{\mathrm{th}}$ is calculated using
interpolated coefficient of thermal expansion, i.e.,

$\begin{split}\bm{\varepsilon}^{\mathrm{th}}&=\alpha(\rho,T)[T-T_{0}]\mathbf{I},\\\
\alpha(\rho,T)&=h_{\mathrm{ss}}(\rho)\alpha_{\mathrm{ss}}(T)+h_{\mathrm{at}}(\rho)\alpha_{\mathrm{at}}.\end{split}$
(11)

Here $\mathbf{I}$ is the $2^{\mathrm{nd}}$-order identity tensor. Meanwhile,
for plastic strain $\bm{\varepsilon}_{\mathrm{pl}}$ and isotropic hardening
model with the von Mises yield criterion is employed. The yield condition is
determined as

$f(\bm{\sigma},p_{\mathrm{e}})=\sigma_{\mathrm{e}}-[{\sigma}_{\mathrm{y}}(T)+\mathsf{H}(T)p_{\mathrm{e}}]\leq
0$ (12)

with

$\sigma_{\mathrm{e}}=\sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}},\quad
p_{\mathrm{e}}=\int\sqrt{\frac{2}{3}\mathrm{~{}d}\bm{\varepsilon}^{\mathrm{pl}}:\mathrm{d}\bm{\varepsilon}^{\mathrm{pl}}},$

where $\sigma_{\mathrm{e}}$ is the von Mises stress. $\mathbf{s}$ is the
deviatoric stress, and ${\sigma}_{\mathrm{{y}}}$ is the isotropic yield stress
when no plastic strain is present. The isotropic plastic modulus $\mathsf{H}$
can be calculated from the isotropic tangent modulus $E_{\mathrm{t}}$ and the
Young’s modulus $E$ as $\mathsf{H}=EE_{\mathrm{t}}/(E-E_{\mathrm{t}})$.
$p_{\mathrm{e}}$ is the effective (accumulate) plastic strain, which is
integrated implicitly by the plastic strain increment
$\bm{\varepsilon}^{\mathrm{pl}}$ employing radial return method and will be
elaborated in Sec. 3.2.

##  3 Simulation setup

###  3.1 Simulation domain and model parameters

In this work, SS processes with a constant beam spot diameter
$D_{\mathrm{L}}$, a beam power $P$ as well as a scan speed $v$ were simulated.
The processing window $P\in[15,~{}30]~{}$\mathrm{W}$$ and
$v\in[75,~{}150]~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ was set in order
to compare it with previous studies [7, 36].
$D_{\mathrm{L}}=D_{\mathrm{FWE2}}=200~{}$\mathrm{\SIUnitSymbolMicro m}$$
(i.e., full-width at $1/e^{2}$) was adopted as the nominal diameter of the
beam spot, within which around $86.5\%$ of the power is concentrated. The
full-width at half maximum intensity (FWHM) then equals to
$0.588D_{\mathrm{FWE2}}=117.6~{}$\mathrm{\SIUnitSymbolMicro m}$$,
characterizing $50\%$ power concentration within the spot. The powder bed is
preheated with a preheating temperature
$T_{\mathrm{0}}=0.4T_{\mathrm{M}}=680~{}$\mathrm{K}$$, which was embodied by
temperature initial condition (IC) and boundary condition (BC). Both powder
and substrate are SS316L. The powder size follows Gaussian distribution with a
mean diameter of $35~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, a standard
deviation of $10~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, and a cut-off
bandwidth of $15~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, which are
consistent with our former works [7, 36]. In-total four layers of SS processes
were performed for each pair of processing parameters (i.e., $P$ and $v$,
hereinafter as $P\text{-}v$ pair). After each layer’s beam scan (i.e., the
beam leaves the prime domain), the simulation continues, and the domain
undergoes a natural cooling stage, lasting twice as long as the scanning
period.

For the non-isothermal phase-field simulations, the powders are firstly
deposited into a $250\times 500\times 400$
$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$ prime simulation domain under a
given gravitational force based on the discrete element method (DEM). The
prime domain contains the free space of
$150~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ height for positioning powders
and a substrate of $250~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ thickness.
Together with the governing equations shown in Eqs. $($1$)$-$($3$)$, the
following BCs are also employed for the process simulations

$\displaystyle\left.\nabla\rho\right|_{\Gamma}\cdot\hat{\mathbf{n}}=0,$ (13)
$\displaystyle\left.\mathbf{K}\cdot\nabla
T\right|_{\Gamma_{\mathrm{S}}\cup\Gamma_{\mathrm{T}}}\cdot\hat{\mathbf{n}}=h_{\mathrm{at}}\left[\underline{h}\left(\left.T\right|_{\Gamma_{\mathrm{S}}\cup\Gamma_{\mathrm{T}}}-T_{0}\right)+\varepsilon\sigma_{\mathrm{B}}\left(\left.T\right|_{\Gamma_{\mathrm{S}}\cup\Gamma_{\mathrm{T}}}^{4}-T_{0}^{4}\right)\right],$
(14) $\displaystyle\left.T\right|_{\Gamma_{\mathrm{B}}}=T_{0}$ (15)

with the convectivity $\underline{h}$, Stefan-Boltzmann constant
$\sigma_{\mathrm{B}}$, the hemispherical emissivity $\varepsilon$, and the
pre-heating (environmental) temperature $T_{0}$. $\Gamma_{\mathrm{T}}$ and
$\Gamma_{\mathrm{B}}$ are correspondingly the top and bottom boundaries of the
simulation domain, and $\Gamma_{\mathrm{S}}$ is the set of all surrounded
boundaries.
$\Gamma_{\mathrm{B}}\cup\Gamma_{\mathrm{S}}\cup\Gamma_{\mathrm{T}}=\Gamma$, as
the schematic shown in the inset of Fig. 1b. $\hat{\mathbf{n}}$ is the normal
vector of the boundary. Eq. $($13$)$ physically shows the close condition for
the mass transfer, restricting no net mass exchange with the environment. Eq.
$($14$)$ shows the heat convection and radiation, which are only allowed via
the pore/atmosphere at the boundary (masked by $h_{\mathrm{at}}$). Eq.
$($15$)$ emulates a semi-infinite heat reservoir under the bottom of the
substrate with constant temperature $T_{0}$, consistently draining heat from
the simulation domain and reducing its temperature back to $T_{0}$.

As summarized in Sec. 2.1, this simulation requests following
parameters/properties: the thermodynamic parameters
$\mathop{}\\!\underline{W}_{\mathrm{sf}}$,
$\mathop{}\\!\underline{W}_{\mathrm{gb}}$,
$\mathop{}\\!\underline{\kappa}_{\rho}$, and
$\mathop{}\\!\underline{\kappa}_{\eta}$; the kinetic properties (diffusivities
and GB mobility) $D_{\mathrm{path}}$ (path = ss, at, sf, gb) and
$G_{\mathrm{gb}}$; and the thermal properties $K_{\mathrm{ss}}$ and
$K_{\mathrm{at}}$. Among them, $\mathop{}\\!\underline{W}_{\mathrm{sf}}$,
$\mathop{}\\!\underline{W}_{\mathrm{gb}}$,
$\mathop{}\\!\underline{\kappa}_{\rho}$, and
$\mathop{}\\!\underline{\kappa}_{\eta}$ are parameterized by the temperature-
dependent surface as well as GB energies, and a nominal diffuse interface
width $l_{\mathrm{gb}}$ as

$\begin{split}&\gamma_{\mathrm{sf}}(T)=\frac{\sqrt{2}}{6}\tau_{\mathrm{sf}}(T)\sqrt{({W}_{\mathrm{sf}}+7W_{\mathrm{gb}})(\kappa_{\rho}+\kappa_{\eta})},\\\
&\gamma_{\mathrm{gb}}(T)=\frac{2\sqrt{3}}{3}\tau_{\mathrm{gb}}(T)\sqrt{{W}_{\mathrm{gb}}\kappa_{\eta}},\\\
&l_{\mathrm{gb}}\approx\frac{2\sqrt{3}}{3}\sqrt{\frac{\kappa_{\eta}}{{W}_{\mathrm{gb}}}}\end{split}$
(16)

with normalized tendencies $\tau_{\mathrm{sf}}(T)$ and $\tau_{\mathrm{gb}}(T)$
that reach unity at $T_{\mathrm{M}}$.
$\mathop{}\\!\underline{W}_{\mathrm{sf}}={W}_{\mathrm{sf}}\tau_{\mathrm{sf}}(T)$,
$\mathop{}\\!\underline{W}_{\mathrm{gb}}={W}_{\mathrm{gb}}\tau_{\mathrm{gb}}(T)$,
$\mathop{}\\!\underline{\kappa}_{\rho}=\kappa_{\rho}\tau_{\mathrm{sf}}(T)$,
and
$\mathop{}\\!\underline{\kappa}_{\eta}=\kappa_{\eta}\tau_{\mathrm{gb}}(T)$.
Constants ${W}_{\mathrm{sf}}$, $W_{\mathrm{gb}}$, $\kappa_{\rho}$,
$\kappa_{\eta}$ also satisfy a constraint
$({W}_{\mathrm{sf}}+{W}_{\mathrm{gb}})/\kappa_{\rho}=6W_{\mathrm{gb}}/\kappa_{\eta}$
derived from the coherent diffuse-interface profile at equilibrium
(Supplementary Note 1 of Ref. [7]). These parameter/properties are
collectively shown in Table 1.

For the thermo-elasto-plastic calculations, a $250\times 250\times 250$
$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$ subdomain was selected from the
center of the prime domain to eliminate the boundary effects (Fig. 1a). The
above momentum balance is subjected to the following rigid support BC as

$\left.\mathbf{u}\right|_{\Gamma_{\mathrm{S}}\cup\Gamma_{\mathrm{B}}}\cdot\hat{\mathbf{n}}=0,$
(17)

restricting the displacement $\mathbf{u}$ along the normal direction of all
boundaries except the top ($\Gamma_{\mathrm{T}}$), which is traction free. As
summarized in Sec. 2.2, temperature-dependent $E$, $\nu$, $\alpha$,
$E_{\mathrm{t}}$, and $\sigma_{\mathrm{y}}$ are as listed in Table 2, where
piecewise linear interpolation was employed to implement their temperature
dependence.

###  3.2 Numerical implementation and parallel computing

The theoretical models are numerically implemented via the finite element
method within the program NIsoS, developed by authors based on MOOSE framework
(Idaho National Laboratory, ID, USA) [53]. 8-node hexahedron Lagrangian
elements are chosen to mesh the geometry. A transient solver with
preconditioned Jacobian-Free Newton-Krylov method (PJFNK) and backward Euler
algorithm was employed for both problems.

For non-isothermal phase-field simulations, the Cahn–Hilliard equation in Eq.
$($1$)$ was solved in a split way, The constraint of the order parameters was
fulfilled using the penalty method. To reduce computation costs, h-adaptive
meshing and time-stepping schemes are used. The additive Schwarz method (ASM)
preconditioner with incomplete LU-decomposition sub-preconditioner was also
employed for parallel computation of the vast linear system, seeking the
balance between memory consumption per core and computation speed [54]. Due to
the usage of adaptive meshes, the computational costs vary from case to case.
The peak DOF number is on order of 10,000,000 for both the nonlinear system
and the auxiliary system. The peak computational consumption is on the order
of 10,000 core-hour. More details about the FEM implementation are shown in
the Supplementary Note 4 of Ref. [7].

For thermo-elasto-plastic simulations, a static mesh was utilized to avoid the
hanging nodes generated from h-adaptive meshing scheme. In that sense, the
transient fields $T$ and $\rho$ of each calculation step were uni-
directionally mapped from the non-isothermal phase-field results (with
h-adaptive meshes) into the static meshes, assuming a weak coupling between
thermo-structural and mechanical problems in this work. This is achieved by
the MOOSE-embedded SolutionUserObject class and associated functions. The
parallel algebraic multigrid preconditioner BoomerAMG was utilized, where the
Eisenstat-Walker (EW) method was employed for determining linear system
convergence. It is worth noting that a vibrating residual of non-linear
iterations would show without the employment of EW method for this work. The
DOF number of each simulation is on the order of 1,000,000 for the nonlinear
system and 10,000,000 for the auxiliary system. The computational consumption
is on the order of 1,000 CPU core-hour.

A modified radial return method [55, 56] was employed to calculate the
integral of the plasticity as well as to determine the yield condition during
the process with the temperature-dependent elasto-plastic properties at any
time $t$ with the time increment $\Delta t$. This method employs the trial
stress $\bm{\sigma}^{\star}$ calculated assuming an elastic new strain
increment $\Delta\bm{\varepsilon}^{\star}$

$\bm{\sigma}^{\star}=\mathbf{C}(T_{t}):\left[\bm{\varepsilon}^{\mathrm{el}}_{t}+\Delta\bm{\varepsilon}^{\star}\right],$
(18)

where the elasticity tensor is obtained under the current temperature field
$T_{t}$. Once the trial stress state is outside of the yield condition (Eq.
$($12$)$), i.e., the plastic flow exists, the stress is then projected onto
the closet point of the expanded yield surface with the normal direction
determined as
$\hat{\mathbf{n}}_{\mathrm{y}}={3\mathbf{s}^{\star}}/{2\sigma_{\mathrm{e}}^{\star}}$
with the von Mises $\sigma_{\mathrm{e}}^{\star}$ and deviatoric trial stress
$\mathbf{s}^{\star}$. Meanwhile, assuming isotropic linear hardening under
$T_{t}$ of every timestep, the amount of the effective plastic strain
increment $\Delta p$ for returning the stress state back to the yield surface
is calculated in an iterative fashion adopting Newton’s method

$\displaystyle\mathrm{d}\Delta
p=\frac{\sigma_{\mathrm{e}}^{\star}-3G(T_{t})\Delta
p_{t}-\left[\mathsf{H}(T_{t})p_{t}+\sigma_{\mathrm{y}}(T_{t})\right]}{3G(T_{t})+v(T_{t})},$
(19) $\displaystyle\Delta p_{(t+\Delta t)}=\Delta p_{t}+\mathrm{d}\Delta p,$
(20) $\displaystyle p_{(t+\Delta t)}=p_{t}+\Delta p_{(t+\Delta t)}$ (21)

where $p_{t}$ and $\Delta p_{t}$ are updated as $p_{(t+\Delta t)}$ and $\Delta
p_{(t+\Delta t)}$ at the end of the timestep ($t+\Delta t$). $G(T_{t})$ and
$\mathsf{H}(T_{t})$ are shear and isotropic plastic modulus calculated under.
Here $G(T_{t})=E(T_{t})/[2+2\nu(T_{t})]$ with $E(T_{t})$ and $\nu(T_{t})$
correspondingly the temperature-dependent Young’s modulus and Poisson ratio.
Knowing that the plastic strain increment
$\Delta\bm{\varepsilon}^{\mathrm{pl}}=\Delta p\hat{\mathbf{n}}_{\mathrm{y}}$
following the normality hypothesis of plasticity [55], the updated stress and
plastic strain at the end of timestep are thereby obtained by

$\begin{split}\bm{\sigma}_{(t+\Delta t)}&=\left\\{\begin{aligned}
&\bm{\sigma}^{\star}-2G(T_{t})\Delta p_{(t+\Delta
t)}\hat{\mathbf{n}}_{\mathrm{y}},\quad&T<T_{\mathrm{A}},\\\
&\mathbf{0},\quad&T\geq T_{\mathrm{A}},\end{aligned}\right.\\\
\bm{\varepsilon}^{\mathrm{pl}}_{(t+\Delta t)}&=\left\\{\begin{aligned}
&\bm{\varepsilon}^{\mathrm{pl}}_{t}+\Delta p_{(t+\Delta
t)}\hat{\mathbf{n}}_{\mathrm{y}},\quad&T<T_{\mathrm{A}},\\\
&\mathbf{0},\quad&T\geq T_{\mathrm{A}},\end{aligned}\right.\end{split}$ (22)

in which the vanishing of the stress as well as the plastic strain beyond a
stress-free temperature $T_{\mathrm{A}}$ (in this work
$T_{\mathrm{A}}=T_{\mathrm{M}}$) is also considered.

Eqs. $($18$)$-$($22$)$ are sequentially executed on every timestep and update
the quantities under $T_{t}$. It is worth noting that the linear-interpolated
$\mathsf{H}$ is implemented for both $p$\- and $T$-dependence

$\mathsf{H}(p_{t},T_{t})=(1-\tau_{t})f(\hat{T}_{i},p_{t})+\tau_{t}f(\hat{T}_{i+1},p_{t})$
(23)

with

$\displaystyle\tau_{t}=\frac{T_{t}-\hat{T}_{i}}{\hat{T}_{i+1}-\hat{T}_{i}},$

where $f(\hat{T}_{i},p_{t})$ is a piecewise function with the grids
$\hat{T}_{i}$ ($i=1,2,...$) and $T_{i}$ is inside the section bound by
$\hat{T}_{i}$ and $\hat{T}_{i+1}$, i.e.,
$T_{t}\in(\hat{T}_{i},\hat{T}_{i+1}]$. Considering reduction of the non-
linearity, the $p$-independent $\mathsf{H}(T)$ was practically employed in the
simulations, as formulated in Eq. $($19$)$.

##  4 Results

###  4.1 Thermo-structural evolution during multilayer SS

Fig. 2a1-a4 present the transient thermo-structural profiles for a single scan
of $P=20~{}$\mathrm{W}$$ and
$v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ where the beam spot is
consistently positioned along scan direction (SD) across various building
layers. In the overheated region ($T\geq T_{\mathrm{M}}$), particles can
undergo complete or partial melting. This prompts molten materials to flow
from convex to concave points, making the fusion of particles possible and,
therefore, forging the fusion zone. In regions where the temperature remains
below the melting point ($T<T_{\mathrm{M}}$), however, the temperature is
still high enough to cause diffusion as the diffusivity grows exponentially
w.r.t. temperature. This can be evidenced by the necking phenomenon among
neighboring particles. Temperature profiles demonstrate a strong dependence on
the local morphology. Notably, the concentrated temperature isolines can be
observed around the neck region among particles, indicating an increased
temperature gradient (up to around 50-100
$\mathrm{K}~{}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-1}$, comparing to 1
$\mathrm{K}~{}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-1}$ in the densified
region). It is worth noting that such thermal inhomogeneity induced by
stochastic transient morphology can be hardly resolved by the simulation works
employing homogenized powder bed [28, 57, 25].

The heat dissipation in the powder bed also changes considerably as upper-
layers are continuously built and processed. This can be noticed by a varying
shape and significance of the overheated region across layers. As a
comparison, the overheated region at 4 layer (Fig. 2a1) is greater than at
initial layer (Fig. 2a4). This can simply be reasoned by increasing porosity
in fused lower-layers that block the heat dissipation. At the initial layer,
this dissipation is affected the least as there is no lower-layer but
substrate. We further probed the temperature history at surface points located
at the center of the scan path, as shown in Fig. 2c1, where the recurring
peaks in every single-layer SS stage (shaded section) present the thermal
cycle. The first peaks of distinctive cycles (probed on different surface
points) nearly match. In an identical cycle (e.g., C1), the peaks gradually
decrease as the SS processes advance. Once the spot leaves the probing point,
the temperature drop is rapid at first (mostly due to the heat conduction
driven by the high-gradient temperature profile) and then becomes moderate
(due to heat convection and radiation) as the cooling stage begins, when the
heat convection and radiation are effective. This steep temperature drop
gradually disappears at points from fused lower-layers as the upper-layers are
continuously built, comparing the thermal cycles at point C1. Meanwhile, the
temperature drop at distinct points gradually coincides as the cooling stage
proceeds, comparing the thermal cycles at points C1-C4.

The multilayer thermo-structural evolution is also significantly affected by
processing parameters, as the profiles of varying beam powers and scan speeds
presented in Fig. 2b1-b4. Increasing the beam power and/or decreasing the
scanning speed is observed to intensify heat accumulation at the beam spot,
resulting in an enlargement of the overheated region. When $v$ is held
constant at $100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$, decreasing $P$
from $30~{}$\mathrm{W}$$ (Fig. 2b1) to $15~{}$\mathrm{W}$$ (Fig. 2b2) results
in a more pronounced overheated region. On the other hand, maintaining a
constant $P=20~{}$\mathrm{W}$$ and increasing $v$ from
$75~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ (Fig. 2b3) to
$150~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ (Fig. 2b4) leads to a
reduction in the overheated region. It is also evident that increasing $P$
and/or decreasing $v$ produces less porosity in fused lower-layers, which may
further change the thermal conditions and improve the heat dissipation. Fig.
2c2 presents the probed thermal cycle of point C1 under various $P$ with
$v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ maintained. One can tell
that the temperature of the case $P=30~{}\mathrm{W}$ quickly drops from a
higher peak to the value coincided with one of the case $P=20~{}\mathrm{W}$ at
the end of scan duration, implying an improved heat dissipation. It should
also be noted that surface point C1, located at the initial layer, experienced
three peaks of temperature that are beyond $T_{\mathrm{M}}$ due to the
extended depth of the fusion zone at higher beam power. In Sec. 5.1 we will
continue the discussion regarding the relation between porosity and processing
parameters.

###  4.2 Stress and plastic strain evolution during multilayer SS process

To analyze the overall stress and plastic strain evolution during SLS, the
average quantities (incl. von Mises stress $\sigma_{\mathrm{e}}$, effective
plastic strain $p_{\mathrm{e}}$, and temperature $T$) in the powder bed are
defined as

$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}=\frac{\int_{\Omega^{\prime}}\rho\sigma_{\mathrm{e}}\mathrm{d}\Omega}{\int_{\Omega^{\prime}}\rho\mathrm{d}\Omega},\quad\bar{p}_{\mathrm{e}}^{\mathrm{P}}=\frac{\int_{\Omega^{\prime}}^{\mathrm{P}}\rho
p_{\mathrm{e}}\mathrm{d}\Omega}{\int_{\Omega^{\prime}}\rho\mathrm{d}\Omega},\quad\overline{T}^{\mathrm{P}}=\frac{\int_{\Omega^{\prime}}^{\mathrm{P}}\rho
T\mathrm{d}\Omega}{\int_{\Omega^{\prime}}\rho\mathrm{d}\Omega},$ (24)

where $\Omega^{\prime}$ is the volume of the simulation domain without the
substrate meshes, and $\rho$ is the OP indicating the substance with $\rho=1$.
They are hereinafter termed as PB-averaged quantities. Fig. 3a presents that
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ develops as the SS process advances
on distinctive layers. When the beam spot enters the domain, accompanied by
the appearance of the overheated region,
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ begins to decrease while
$\overline{T}^{\mathrm{P}}$ continues to rise to its peak. This is attributed
to losing structural integrity caused by full/partial melting within the
overheated region, resulting in zero stress as shown in Fig. 3b1. The areas
surrounding the overheated region also exhibit relatively low stress due to a
significant reduction in stiffness at elevated temperatures. Stress
accumulates faster around concave features, such as surface depressions and
particle sintering necks, compared with the stress around convex features as
well as unfused powders away from the fusion zone. This is due to localized
temperature gradients, as explained in Sec. 4.1. As cooling progresses, stress
decreases in convex features and unfused powders, yet remains concentrated
around concave features (see Figs. 3b2-3b3). The stress at the end of each
cooling stage (Figs. 3b2-3b6) serves as the residual stress of corresponding
processed layers. After the deposition of new layers,
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ has an instant drop due to the
addition of stress-free substances (powders), then a new cycle begins. It is
worth noting that the peak values of the stress cycle (which are also residual
stresses of each processed layer) are almost identical, implying that each
single-layer SS generates nearly the same amount of residual stress by average
in the powder bed. For upper-layers, $\sigma_{\mathrm{e}}$ experiences an
extra reduction in the cooling stage, which may be attributed to the more
delicate variations in the on-site stress components and will be discussed in
Fig. 4.

The accumulated plastic strain $p_{\mathrm{e}}$, on the other hand, presents
an overall increasing tendency vs time during the single-layer SS and cooling
stage, as shown in Fig. 3c. This distinguishes the plastic strain history from
the stress cycle, which suffers a reduction during SS due to loss of
structural integrity. The continuous accumulation of $p_{\mathrm{e}}$ leads to
a locally concentrated profile in the vicinity of the melt zone (see Fig.
3d2-3d3). The high temperature gradient at the front and bottom of the
overheated region during single layer SS also results in the localized
$p_{\mathrm{e}}$ surrounding the fusion zone (see Fig. 3d2-3d6), where the
existing high thermal stress locally activates the plastic deformation. A
similar explanation applies to the concentrated $p_{\mathrm{e}}$ around pores
and concave features such as sintering necks near the fusion zone. In
contrast, unfused powders and the substrate away from the melt zone show
minimal $p_{\mathrm{e}}$, since the thermal stress at these locations is
insufficient to induce material plastification due to lower local
temperatures.

We probed the history of the stress components at six points, where
L${}_{1}0$-L${}_{4}0$ locates at the center of 1-4 layers, L${}_{1}1$ locates
at the boundary of the fusion zone, and L${}_{1}2$ locates outside the fusion
zone in the 1 layer, as shown in the inset of Fig. 4. The results mainly show
the difference of the stress field inside and outside the fusion zone. Since
the temperature cools down from above the melting temperature to pre-heating
temperature, the significant increase of the elasticity leads to a significant
increase of $\sigma_{\mathrm{e}}$ for temperature above $0.65T_{M}$. As a
comparison, the increment of the $\sigma_{\mathrm{e}}$ becomes slower for the
point $L_{2}0$ and the point $L_{2}0$. It can also be seen that
$\sigma_{\mathrm{e}}$ for the mentioned three points during SS of 2- and 3
layers are very similar when they are all outside the fusion zone. Similar
results can be observed for the development of normal stress components
$\sigma_{ii}$ ($i=x,~{}y,~{}z$). The main differences in normal stress for the
three points are mostly for the SS process of the 1 layer. $\sigma_{ii}$ are
tensile in the fusion zone while compressive outside the fusion zone. This is
due to the thermal stress in the fusion zone being tensile while it is
compressive outside the fusion zone. Curves for the 2 and 3 layers are similar
because the three points are all outside the fusion zone. Shear stress
components $\sigma_{ij}$ ($i,j=x,~{}y,~{}z$; $i\neq j$) remain small inside
and outside the fusion zone. The variation of the shear stress is more visible
at the boundary of the fusion zone, as shown in Fig. 4a2. This is due to a
large gradient of the thermal strain and material mechanical properties across
the boundary of the fusion zone. With a larger variation of the shear stress
components, variation of $\sigma_{\mathrm{e}}$ at the highest temperature
during SS process of each layer is also more visible. The saturation of the
stress field can be observed after the 2 layer. It means that the upper layer
has limited influence on the plastic deformation of previous layers which
experience lower temperatures and lower temperature gradients outside the
fusion zone.

For points in the layer-wise direction ($L_{2}^{0}$ to $L_{4}^{0}$ at the
boundary of each layer), the results in Fig. 4b show a direct comparison
between the development of the stress field of a point in the microstructure
during SS process of different layers. It is worth noting that the development
of the $\sigma_{\mathrm{e}}$ of each point is very similar at the beginning.
The saturation of the $\sigma_{\mathrm{e}}$ for each point can be observed in
the layer-wise direction when the points are outside the fusion zone. However,
variation of the $\sigma_{\mathrm{e}}$ is much more visible when the fusion
zone is moving close to the point. This indicates that the point is located
around the boundary of the fusion zone, which encounters a large variation of
the mechanical properties of the material w.r.t. the high-gradient temperature
temperature, and thus undergoes a large variation of the shear stress.

###  4.3 Powder-resolved mesoscale residual stress formation mechanism

For SS-processed porous microstructures, the residual stress is generally
related to stress concentration under thermal loading in the microstructure.
In Fig. 5a3 and 5a5, the particles undergo a low degree of fusion (due to low
level of overheating) due to relatively low beam power or high scan speed. The
stress concentration mostly occurs at the necking region, leading to
accumulated plastic strain and residual stress. In Fig. 5a1 and 5a4, with a
higher degree of fusion due to a relatively high beam power or low scan speed,
accumulated plastic strain and residual stress tends to concentrate on both
necking region and inter-layer region (as indicated by the black lines in Fig.
5c and 5d). Note that the inter-layer boundary has a high surface roughness
which also causes stress concentrations. In Fig. 5a2, with the highest beam
power within the processing window, residual stress can be observed in the
whole fusion zone, which makes it difficult to identify the residual stress
concentration. Besides, accumulated plastic strain and residual stress can be
observed on the boundary of the substrate and the 1 layer, especially at high
beam power. This is because a high beam power generates a much larger fusion
zone, which even penetrates into the substrate. When the substrate with
continuum material is melted, it generates a large plastic deformation and
thus the residual stress, which can be significantly larger than those in the
porous microstructure, as shown in Fig. 5a2. These results show that the
residual stress is directly related to the degree of fusion and thus related
to SS processing parameters.

To gain a better understanding of the distribution of residual stress, we
analyzed the stress state of a few representative positions in the
microstructure, as shown in Fig. 5c and 5d. The former surfaces of fused
layers are marked by black lines, and the bottom of the fused strut, a.k.a.,
the fusion zone boundary (FZB) of the 1 layer, is indicated by white lines.
Fig. 5c1 and 5d1 show that the residual stress and the accumulated plastic
strain are solely concentrated in the necking region and inter-layer region
for the case with a low degree of fusion. In Fig. 5c2, with the high degree of
fusion (the smallest porosity acheived in the selected processing window), the
residual stress and the accumulated plastic strain still tend to be solely
concentrated in the inter-layer region, yet less distinguishable compared to
the former case. By comparing Lamé’s stress ellipsoids at six points on the
boundaries of different layers, it is worth noting that the stress states at
the boundary of the two layers are very similar, e.g. the stress states of P2
and P3 are similar to that of P5 and P6, no matter how high the degree of
fusion is. It means that the residual stress is still formulated due to stress
concentration at the inter-layer region. On the other hand, the stress states
of P1 and P4 at the top layer are quite different. This is because the
residual stress of the top layer is directly determined by the thermal loading
which depends on the surface morphology and morphology-induced thermal
inhomogeneity.

Based on the above observations, we propose a powder-resolved mesoscale
residual stress formation mechanism, which is summarized in the schematic
illustrated in Fig. 6. The residual stress formulated in porous
microstructures manufactured by a multilayer SS process contains two primary
contributions:

  1. 1.

Residual stress is directly caused by the stress concentration at the necking
region of partially melted particles under thermal loading. The partially
melted particles are inter-connected via necking regions, which is the weak
link of the microstructure. Therefore, both the thermal expansion in the
overheated region and the thermal contraction during the cooling stage cause
severe plastic deformation at the necking region, which is one primary source
of the accumulated plastic strain and residual stress.

  2. 2.

Residual stress due to interaction between the upper- and lower-/fused layers
in the layerwisely build-up process. In the cooling stage, the shrinkage of
the upper-layer after overheating results in tensile stress on itself and
compressive stress on the lower-/fused layer, which causes plastic deformation
of the porous structure, especially at the inter-layer region with stress
concentration due to a high surface roughness, which is the other primary
source of the accumulated plastic strain and residual stress, as indicated by
the white lines in Fig. 6a.

The proposed mechanism is based on the detailed simulation results using the
powder-resolved thermo-mechanical model, which evidently demonstrates the
formation and distribution of residual stress in the porous structure. The
accumulated plastic strain results in structural distortion of the fused
strut, which is schematically illustrated in Fig. 6a and demonstrated by the
contour plot of deformation component $u_{z}$ in Fig. 6b. Compared to the
aforementioned TGM and CDS models proposed by Mercelis and Kruth [22]
correspondingly, which ignored the difference between the part manufactured by
SM and SS, the partial melting as the main fusion mechanism plays an important
role in the residual stress in the microstructure. There are some
similarities, the concept of the TGM model also applies to the proposed
mechanism because the temperature gradient indeed directly leads to residual
stress. However, in the SS process, since the power bed is fully relaxed
before fusion, the particle provides very weak constrain on the fusion zone
and thus with limited plastic deformation in the particle. The plastic
deformation is majorly accumulated during the cooling stage due to thermal
contraction and stress concentration at the necking region. For the CDS model,
the proposed model also considered the shrinkage of the top layer resulting in
a tensile residual stress in the top layer and a compressive residual stress
in the previous layers. However, the top layer in the simulated SS process has
very complex boundaries with previous layers. These boundaries are the source
of the stress concentration which leads to the dominant residual stress at the
inter-layer regions.

Based on the TGM model and the CDS model, Mercelis and Kruth [22] also
introduced a theoretical model to predict the relationship between residual
stress and the part height. Similarly, we can use the proposed powder-resolved
model and mesoscale residual stress formation mechanism to predict the
dependence of the residual stress on the microstructure porosity, which is
also a characteristic geometric parameter of a porous microstructure. Since
the porous structure induces complex inhomogeneity of both material properties
and thermal and mechanical loading, and we also need to consider the porosity
is directly related to SS process parameters, it is more straightforward to
propose phenomenological models based on our simulation results to predict the
relationship between the residual stress and SS process parameters, as will be
discussed in the next section.

##  5 Discussion

###  5.1 Phenomenological relation for porosity control

Controlling the porosity of processed sample through the tuning of the
processing parameters plays a central role in tailoring the end-up performance
of the porous materials in AM, as many properties are determined by (or
related to) porosity, including but not limited to mechanical strength,
permeability, acoustic/optical absorptivity and various effective
conductivities. We start the discussion with phenomenologically relating the
porosity to the processing parameters (in this work $P$ and $v$). In this
work, a nominal porosity is defined as

$\varphi=\frac{\int_{\Omega^{\prime\prime}}\rho\sigma_{\mathrm{e}}\mathrm{d}\Omega}{\int_{\Omega^{\prime\prime}}\rho\mathrm{d}\Omega}$
(25)

with the substance OP $\rho$. $\Omega^{\prime\prime}$ represents the volume of
a post-processed simulation domain with surface and unforced powders
sufficiently removed (termed as “virtual polishing”). Fig. 7a presents the
local porosity that evaluated segment-wise along $z$\- (BD) and $y$-direction
(perpendicular to SD) to help identify the representative domain with a width
$W_{\mathrm{R}}$ and a height $H_{\mathrm{R}}$ for porosity calculation, while
the complete length (250 $\mathrm{\SIUnitSymbolMicro}\mathrm{m}$) along
$x$-direction (SD) is included due to exhibited relatively minor variations in
local porosity [36]. We then conduct the virtual polishing on simulated
multilayer SS-processed microstructures, and proceed calculating their
porosity using Eq. $($25$)$. The resulting $P\text{-}v$ map of porosity is in
shown Fig. 7c. Selected microstructures illustrated in Figs. 7b2-7b4 for the
cases varying $P$ while maintaining $v$; Figs. 7b5-7b7 for the cases varying
$v$ while maintaining $P$. Fig. 7b1 is the microstructure under a reference
processing parameter ($P=20~{}$\mathrm{W}$$,
$v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$), while Figs. 7b8 and 7b9
are respectively the one with maximum and minimum porosity. It demonstrates
that improved densification can be achieved by either increasing $P$ or
decreasing $v$, where smoothed surface morphology implies an enhanced partial
melting. Porosity drops from 31% to 18% for the increase of $P$ from 15 to 30
$\mathrm{W}$ ($v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$) and from 29%
to 21% for the decrease of $v$ from 150 to 75
$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$ ($P=20~{}$\mathrm{W}$$). The
presented tendencies imply a possible allometric relation $\varphi\propto
P^{-m}v^{n}$ with positive indices $m$ and $n$. On the other hand, combining
beam power and scan speed as one characteristic quantity, $P/v$ is widely
adapted to define a specific energy density, notably the volumetric energy
input as

$U=\frac{P}{HWv},$ (26)

where $H$ is the powder layer thickness (here
$H=60~{}$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$) and $W$ is the scan track
width (here $W$ takes $D_{\mathrm{FWHM}}$). Evidently, $\varphi$ exhibits an
overall decrease tendency with rising $U$, meaning lower porosity can be
achieved as higher specific energy input.

We then examined the simulation results with a proposed phenomenological
relation following Refs. [58, 7], read as

$\ln[1-\varrho(\varphi)]=-K_{\varrho}U$ (27)

with a densification factor defined as

$\quad\varrho(\varphi)\equiv\frac{\varphi_{0}-\varphi}{\varphi_{0}-\varphi_{\text{min}}}.$

$\varrho$ indicates the ratio between a reduced porosity w.r.t. the maximum
porosity reduction in the chosen processing window, as $\varphi_{0}$ and
$\varphi_{\text{min}}$ represent respectively the initial and minimum achieved
porosity. $\varphi_{\text{min}}$ usually varies between 3% to 30% for metallic
materials [58]. In pursuit of uniformity, we have chosen $\varphi=3\%$ in this
study, aligning with our previous single-layer SS simulation [7].
$\varphi_{0}$ is, however, difficult to be determined for multilayer SS
simulation since there is sorely one or two particles in thickness for every
new layer while the old layers have been fused. In this sense, we evaluate
$\varphi_{0}$ by three routes:

  1. 1.

Assuming that the porosity of the powder bed region, which is away from the
beam spot (i.e., without significant densification), is the initial one of the
powder bed, $\varphi_{0}$ can be read as the converged value from segment-wise
porosity evaluation along $y$-direction (Fig. 7a2), which is
$\varphi_{0}=27.8\%\sim 33.6\%$.

  2. 2.

As the particle size distribution in this work is assumed to be Gaussian-type,
$\varphi_{0}$ is then estimated by statistic random-close-packing model of
spherical particle as
$\varphi_{0}=0.366-0.0257\left({\varsigma_{d}}/{\mu_{d}}\right)$ with
$\mu_{d}$ and $\varsigma_{d}$ the mean and standard deviation of the particle
size [59, 57]. Taking $\mu_{d}=20~{}$\mathrm{\SIUnitSymbolMicro m}$$ and
$\varsigma_{d}=5~{}$\mathrm{\SIUnitSymbolMicro m}$$, it can be calculated as
$\varphi_{0}=36.0\%$.

  3. 3.

We also piled multiple densification-free powder stacks using DEM method with
the same domain volume fraction of the particles (48.0%) as the overall one of
the multilayer SS simulation, in which each layer deposits powders of 12%
domain volume fraction. The measured nominal porosity is thereby regarded as
the initial one, as $\varphi_{0}=40.0\%\sim 40.5\%$.

$\varphi_{0}$ from route (i) directly reflects the porosity of on-site unfused
powders, but the influences from the potential necked particles due to thermal
processes cannot be well eliminated. $\varphi_{0}$ from route (ii) is based on
the statistics of random-closed-packed particles, which is rather ideal
compared with the practical powder bed created by powder spreading. Since
route (iii) creates particle stack via simulating the deposition in powder
spreading, the calculated $\varphi_{0}$ may be still inflated as the
morphological variability of the deposition surface is missing. Considering
all of these factors, we have selected $\varphi_{0}=36.0\%$, which stands at
the midpoint among all the assessed values.

In Fig. 7d, linear regression between $-\ln(1-\varrho)$, calculated from
simulated $\varphi$, and $U$ is presented. For comparison, regressions on data
by single-layer SS simulation and experiment are also illustrated [58, 7]. The
regressed line together with the 95% confidence interval
($\mathrm{CI}_{95\%}$) of the multilayer SS simulations is right in between
those of the single-layer SS simulation and experiments. As the coefficient
$K_{\varrho}$ is related to the material and the size distribution of the
powders, the multilayer result $K_{\varrho}=0.016\pm
0.001~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$ demonstrates an improved
coherence with the experimental one
$K_{\varrho}=0.013~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$ compared
the single-layer one
$K_{\varrho}=0.019~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$, which can
suffer from inflated porosity mostly due to insufficient volume microstructure
for porosity calculation.

###  5.2 Phenomenological relations for residual stress and plastic strain
control

Taking the PB-averaged residual von Mises stress
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and effective plastic strain
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ (defined in Eq. $($24$)$) at the end of
the simulations, Fig. 8 presents the distributions of
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ w.r.t. the $P$ and $v$ in the chosen
processing window. Generally, the rise in
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ corresponds to an increase in the
specific energy input $U$ (i.e., increase $P$ with maintained $v$ or increase
$v$ with maintained $P$), resulting in further concentrated residual stress
around around concave features (surface depressions and particle sintering
necks) and across layers, as already presented in Fig. 5a and 5b. To
understand the relationship between $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$
and $U$, we interpreted $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ as the stored
mechanical energy density after the multilayer SS processes, which can be
regarded as the residue of the energy density imported via the beam-induced
thermal effect in the powder bed after all types of in-process dissipation. In
this regard, the nonlinear regression analysis of an energy conversion law
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}=\sigma^{\mathrm{P}}_{\infty}\\{1-\exp[{{-\frac{1}{U^{\mathrm{P}}_{\sigma}}(U-U^{\mathrm{P}}_{\mathrm{th}})}}]\\}$
was conducted on the simulation data. The parameters
$U^{\mathrm{P}}_{\mathrm{th}}$ and $\sigma^{\mathrm{P}}_{\infty}$ adopt the
physical meanings as the volumetric energy input at stress-free state and the
saturated stress at infinity energy input, respectively, and
$\sigma^{\mathrm{P}}_{\infty}/U^{\mathrm{P}}_{\sigma}$ is the increasing rate
(slope) of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ vs. $U$ at stress-free
state, as shown in the inset of Fig. 8c. The result in Fig. 8c presents a high
correlation ($R^{2}=98.41\%$) between simulated
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ vs. $U$ with narrow confidence
interval, demonstrating the applicability of proposed energy conversion law in
predicting the residual stress in the powder bed.

Unlike the residual stress, the effective plastic strain $p_{\mathrm{e}}$ is a
measure of cumulative plastic deformation at any given moment during the
process and lacks a clear physical picture of its relationship with volumetric
energy input $U$. Therefore, the nonlinear regression analysis of an
allometric scaling law
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}=C^{\mathrm{P}}_{p}(U)^{I^{\mathrm{P}}_{p}}$
with parameters $C^{\mathrm{P}}_{p}$ and $I^{\mathrm{P}}_{p}$ was conducted to
phenomenologically relate $\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ to $U$, as shown
in Fig. 8d. The analysis gives a relatively low correlation coefficient
$R^{2}=87.4\%$ of the $\bar{p}_{\mathrm{e}}^{\mathrm{P}}(U)$ relation compared
with one of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}(U)$, with an expanded
confidence interval at high $U$ range. Regressed $I^{\mathrm{P}}_{p}=1.17\pm
0.13$ indicates an almost linear scalability of accumulated plastic strain in
the processed powder bed on energy input via beam scan. It is worth noting
that the proposed scaling law can be challenged as $U$ appears to not being
able to uniquely identify $\bar{p}_{\mathrm{e}}^{\mathrm{P}}$, in other words,
$U$ of certain value can correspond to multiple $p_{\mathrm{e}}$ value, as
also depicted in Fig. 8d. Nonetheless, this scaling law can find it feasible
in estimating strength of in-process plastification of the microstructure at
given specific energy input.

Since both $\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ and
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ consider a complete powder bed with
unfused particles, we further defined two average quantities that only take
the residual stress and effective plastic strain inside the fused strut into
account, denoted respectively $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$. A historical indicator $\xi$ was added in
the system to emulate the phenomenological fusion of the strut during
multilayer SS process following our former work [57, 37], which is initialized
as zero and irreversibly transitions as one once $T\geq T_{\mathrm{M}}$.
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ are then calculated as

$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}=\frac{\int_{\Omega}\xi\sigma_{\mathrm{e}}\mathrm{d}\Omega}{\int_{\Omega}\xi\mathrm{d}\Omega},\quad\bar{p}_{\mathrm{e}}^{\mathrm{S}}=\frac{\int_{\Omega}\xi
p_{\mathrm{e}}\mathrm{d}\Omega}{\int_{\Omega}\xi\mathrm{d}\Omega}.$ (28)

Fig. 9 presents the distribution of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$
and $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ w.r.t. the $P$ and $v$ in the chosen
processing window. Contrasting with those shown in Fig. 8, a notable
distinction in the maps presented in Fig. 9c and 9d is the emergence of
regions where the selected $P$ and $v$ fail to form continuous fused strut.
Meanwhile, the profiles of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ receive significant influences from the
size of the strut. Comparing Fig. 9a1 with 9a2-9a3 and 9a4-9a5, the increase
of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ follows the direction of
increasing $U$ as well, which simultaneously leads to the enlargement of the
fused strut. As depicted in Figs. 5c and 9a, concentrated residual stress is
primarily found within the upper-layers of the strut. When the strut size is
enlarged, it encompasses a larger volume with elevated maximum
$\sigma_{\mathrm{e}}$ and high-$\sigma_{\mathrm{e}}$ region, resulting in an
increased $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ with the rise in $U$.
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ exhibits similar tendency w.r.t. $P$ and
$v$ as $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$, with an increase in $U$
resulting in an increase in $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$. Nonetheless,
regions with highly accumulated $p_{\mathrm{e}}$ locates at the bottom of each
layer’s fusion zone, as illustrated in Figs. 5d and 9b. At high $U$, such
accumulation intensifies, especially at the bottom of the strut (which is also
the bottom of the initial layer’s fusion zone). Simultaneously, enlarged depth
of a fusion zone indicates an extended remelting in the fused lower-layers,
which removes some accumulated $p_{\mathrm{e}}$ by the process or former layer
within the strut, as shown in Fig. 6a. It also concentrates the
high-$p_{\mathrm{e}}$ region further to the bottom of the strut, notably the
profile in Fig. 9b4. Eventually, rise of $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$
w.r.t. $U$ is relatively less “rapid” then one of
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$, evident in slightly sparser contours
in Fig. 9d compared to 9c.

Nonlinear regression analyses were also conducted on
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ employing the proposed energy conversion
law for residual stress and scaling law for effective plastic strain. Results
are correspondingly presented in Fig. 9e and 9f. Notably, comparing to those
of $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}(U)$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}(U)$, the correlation of
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}(U)$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}(U)$ present decline, respectively, with
enlarged confidence interval at both low and high $U$ ranges. This attributes
to the removal of the influences from the unfused particles. Moreover, for
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}(U)$, regressed threshold energy input
presents an negative value as $U_{\mathrm{th}}^{\mathrm{S}}=-11.11\pm
17.43~{}$\mathrm{J}~{}\mathrm{m}\mathrm{m}^{-3}$$, indicating a required
energy output to achieve stress-free state. This is explainable as the
$\sigma_{\mathrm{e}}$ already exists in a just-formed strut. In other word,
primary SS process (i.e., the SS without subsequent thermal post-process to
release residual stress) cannot achieve samples with stress-free strut.
Regressed saturated stress $\sigma^{\mathrm{S}}_{\infty}=260.05\pm
24.5~{}$\mathrm{M}\mathrm{P}\mathrm{a}$$ also presents decline comparing to
the regressed one on PB-averaged data. For
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}(U)$, a reduced regressed index of the
scaling law $\bar{p}_{\mathrm{e}}^{\mathrm{S}}(U)$, i.e.,
$I^{\mathrm{R}}_{p}=0.30\pm 0.06$ in Fig. 9f, demonstrates a sublinear
scalability of $\bar{p}_{\mathrm{e}}^{\mathrm{S}}(U)$ compared with one of
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}(U)$, which is almost linear. It also
reflects a reduced growth rate of $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ at high
$U$, and the underlying intensified remelting, which removes some accumulated
$p_{\mathrm{e}}$ within the strut while further concentrates $p_{\mathrm{e}}$
around strut bottom, and the comparably faster increase in strut size shall be
the reason. Nonetheless, information conveyed by
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}(U)$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}(U)$ is more relevant for practical
application, as unfused particles shall be removed during post-process of a
practical SS. Further examination and validation of the proposed laws for
residual stress and plastic strain w.r.t. specific energy input are expected
in the future numerical and experimental studies.

##  6 Conclusion

In this work, we proposed a powder-resolved multilayer simulation scheme for
producing porous materials using selective sintering, combing FEM-based non-
isothermal phase-field simulation and themo-elasto-plastic simulation with
temperature-dependent material properties. This work has presented the
mesoscopic evolution of stress and plastic strain on a transient thermal-
microstructure under various beam power ($P$) and scan speed ($v$). Process-
property relationships between porosity, residual stress and effective plastic
strain and the volumetric energy input ($U\propto P/v$) have also been
demonstrated and discussed. The following conclusions can be drawn from the
present work:

  1. 1.

We proposed in this work a new powder-resolved mesoscopic residual stress
formation mechanism for porous materials manufactured by the SS process, which
collectively lead to the structural distortion appeared in fused strut. It was
demonstrated with simulation results that the stress concentration at the
necking region of the partially melted particles and inter-layer region
between different layers provide dominant accumulated plastic strain and
residual stress in the porous material.

  2. 2.

Based on the proposed residual stress formation mechanism, we examined the
proposed phenomenological relation between the porosity (densification) and
the volumetric energy input $U$. Regression analysis on the resulting porosity
from multilayer SS simulations suggested an improved coherence with the
experimental data, as the regressed densification coefficient
$K_{\varrho}=0.016\pm 0.001~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$ in
this work comparing to the experimental
$K_{\varrho}=0.013~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$ and one
from our former single-layer simulation
$K_{\varrho}=0.019~{}$\mathrm{m}\mathrm{m}^{3}~{}\mathrm{J}^{-1}$$.

  3. 3.

Two types of average quantities, namely PB-averaged
($\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}$) and strut averaged ones
($\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$), were defined to characterize the
residual stress and plastic strain within the powder bed and fused strut,
correspondingly. The relationships between these quantities and volumetric
energy input ($U$) are unveiled by conducting nonlinear regression analyses.
The average residual stress ($\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$) relates to $U$ by the energy
conversion law, while the average effective plastic strain
($\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ and $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$)
by the allometric scaling law.

  4. 4.

Attributing to the removal of influences from unfused particles, the
correlation of the relations $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ present drops compared with one of
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}(U)$ and
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}(U)$, respectively. Saturation behavior is
observed on both $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}(U)$ and
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}(U)$, while the linear scalability in
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}(U)$ degenerates into sublinear one in
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}(U)$, demonstrating a reduced growth
rate of $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ at high $U$.

Despite the feasibility of the multilayer simulation scheme in recapitulating
mesoscopic formation of porosity, residual stress and plastic strain under
given processing parameters; and the proposed mechanism in explaining the
structural distortion of SS-produced samples, several points should be further
examined and discussed in future works:

  1. 1.

The present work omits the consideration of chronological-spatial distribution
of the thermal-elasto-plastic properties among polycrystals, as the properties
such as elasticity and crystal plasticity vary from grain to grain with
distinct orientations.

  2. 2.

The present findings are examined at relatively low specific energy input $uv$
and, correspondingly, low generated residual stress and accumulated plastic
strain. Pore formation is also limited to lack-of-fusion mechanism. It is
anticipated to conduct further simulations with relatively high energy input,
where the mechanisms like keyholing co-exist with high residual stress and
plastic strain. Extension of the proposed mechanism into the high-$U$ range
together with the following examination and validation are also expected.

## Acknowledgements

Authors acknowledge the financial support of German Science Foundation (DFG)
in the framework of the Collaborative Research Centre Transregio 270 (CRC-TRR
270, project number 405553726, sub-projects A06 and A07), the Research
Training Groups 2561 (GRK 2561, project number 413956820, sub-project A4), the
Priority Program 2122 (SPP 2122, project number 493889809). X. Zhou
acknowledges the support from the National Natural Science Foundation of China
(project number 12302231), Sichuan Science and Technology Program (project
number 2023NSFSC0910), and China Postdoctoral Science Foundation (project
number 2023M732433). The authors also greatly appreciate the access to the
Lichtenberg II High-Performance Computer (HPC) and the technique supports from
the HHLR, Technische Universität Darmstadt. The computating time on the HPC is
granted by the NHR4CES Resource Allocation Board under the project
“special00007”. Y. Yang also greatly thanks Dr. Binbin Lin for helping with
the setup of layerwise powder deposition.

## Data Availability

The authors declare that the data supporting the findings of this study are
available within the paper. Source codes of MOOSE-based application NIsoS and
related utilities are cured in the online repository
bitbucket.org/mfm_tuda/nisos.git. The simulation results, statistics and
metadata are cured in the online dataset (DOI: xxxx/zenodo.xxxxxx).

Table 1: Material properties of the bulk SS316L used in the non-isothermal
phase-field simulations. Here $\mathfrak{R}$ represents the ideal gas
constant.

Properties | Expressions ($T$ in K) | Units | References |
---|---|---|---|---
$T_{\text{M}}$ | $\sim 1700$ | K |  |
$\gamma_{\mathrm{sf}}$ | $10.315-5.00\mathop{}\\!\times 10^{-3}T$ | J/m2 | [60] |
$\gamma_{\mathrm{gb}}$ | $13.018-7.50\times\mathop{}\\!\times 10^{-3}T$ | J/m2 | [60] |
$l_{\mathrm{gb}}$ | $2\mathop{}\\!\times 10^{-6}$ | m |  |
$D_{\text{sf}}$ | $0.40\text{exp}\left(-2.200\times 10^{5}/\mathfrak{R}T\right)$ | m2/s | [61] |
$D_{\text{gb}}$ | $2.40\times 10^{-3}\text{exp}\left(-1.770\times 10^{5}/\mathfrak{R}T\right)$ | m2/s | [61] |
$D_{\text{ss}}$ | $2.17\times 10^{-5}\text{exp}\left(-2.717\times 10^{5}/\mathfrak{R}T\right)$ | m2/s | [62] |
$G_{\text{gb}}$ |  $3.26\times 10^{-3}\text{exp}\left(-1.690\times 10^{5}/\mathfrak{R}T\right)$ * | m4/(J s) |  |
$M_{\text{M}}$ |  $\sim 3.45\times 10^{-13}$ ** | m5/(J s) |  |
$K_{\text{ss}}$ | $10.292+0.014T$ | J/(s m K) | [63] |
$K_{\text{at}}$ | $\sim 0.06$ | J/(s m K) | [64] |
$c_{\text{ss}}$ | $3.61\times 10^{6}+1272T$ | J/(m3 K) | [63] |
$c_{\text{at}}$ | 717.6 | J/(m3 K) | [65] |
$\mathcal{L}_{\text{M}}$ | $2.4\times 10^{9}$ | J/m3 | [63] |

  * *

Activation energy is obtained from [66] while the prefix factor is estimated
as unity at $T_{\text{M}}$ after normalization.

  * **

Estimated as
$100D_{\text{sf}}/2(\mathop{}\\!\underline{W}_{\mathrm{sf}}+\mathop{}\\!\underline{W}_{\mathrm{gb}})$.

Table 2: Temperature-dependent mechanical properties of the bulk SS316L used in the thermo-elasto-plastic simulations [67, 68] $T$ (K) | 298 | 873 | 1073 | 1473 | 1623 | $\geq 1773$
---|---|---|---|---|---|---
$\nu$ | 0.33 | 0.35 | 0.36 | 0.38 | 0.39 | 0.40
$E$ (MPa) | $2.00\mathop{}\\!\times 10^{5}$ | $1.35\mathop{}\\!\times 10^{5}$ | $7.75\mathop{}\\!\times 10^{4}$ | $1.21\mathop{}\\!\times 10^{4}$ | $6.14\mathop{}\\!\times 10^{3}$ | $200$
$\sigma_{\mathrm{y}}$ (MPa) | $3.45\mathop{}\\!\times 10^{2}$ | $2.12\mathop{}\\!\times 10^{2}$ | $1.99\mathop{}\\!\times 10^{2}$ | $100$ | $50$ | 5
$E_{\mathrm{t}}$ (MPa) | $5.89\mathop{}\\!\times 10^{3}$ | $1.70\mathop{}\\!\times 10^{3}$ | $1.40\mathop{}\\!\times 10^{3}$ | $100$ | $10$ | $1$
$\mathsf{H}$ (MPa) | $6.07\mathop{}\\!\times 10^{3}$ | $1.72\mathop{}\\!\times 10^{3}$ | $1.43\mathop{}\\!\times 10^{3}$ | $101$ | 10 | 1
$\alpha$ (1/K) | $1.20\mathop{}\\!\times 10^{-5}$ | $1.30\mathop{}\\!\times 10^{-5}$ | $1.32\mathop{}\\!\times 10^{-5}$ | $1.36\mathop{}\\!\times 10^{-5}$ | $1.38\mathop{}\\!\times 10^{-5}$ | $1.40\mathop{}\\!\times 10^{-5}$

Figure 1: (a) The workflow of multilayer selective sintering (SS) simulation,
incl. the layerwise powder deposition, the non-isothermal phase-field
simulation and subsequent thermo-elasto-plastic calculation. (b) Prime
simulation domain with four deposited-processed powder layers. Inset: domain
boundaries and their denotations.  Figure 2: (a1)-(a4) Transient thermo-
structural profiles under a beam power $20~{}$\mathrm{W}$$ and a scan speed
$100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$ with the beam spot
consistently positioned across various layers. The beam spot size
characterized by $D_{\mathrm{L}}$ and $D_{\mathrm{FWHM}}$ is indicated.
(b1)-(b4) Transient thermo-structural profiles of the 4 layer under various
beam power and scan speed with the beam spot consistently positioned. The
overheated regions ($T\geq T_{\mathrm{M}}$) are illustrated by a continuous
color map. (c1)-(c2) Temperature history of the selected surface points across
layers at various beam power with scan speed maintained at 100
$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$. Single-layer SS (shaded sections)
and cooling phase are also denoted. Inset: the location of the points.  Figure
3: (a) Calculated von Mises stress ($\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$)
and temperature ($\overline{T}^{\mathrm{P}}$) in powder-bed average vs. time
with the profile of $\sigma_{\mathrm{e}}$ at the denoted states shown in
(b1)–(b6). (c) Calculated effective plastic strain in powder-bed average
($\bar{p}_{\mathrm{e}}^{\mathrm{P}}$) vs. time with the profile of
$p_{\mathrm{e}}$ at the denoted states shown in (d1)–(d6). Figure 4: Probed
history of (a1)-(a3) stress components at surface points L${}_{1}0$,
L${}_{1}2$ and L${}_{1}3$, and (b) von Mises stress $\sigma_{\mathrm{e}}$ at
surface points L${}_{1}0$, L${}_{2}0$, L${}_{3}0$ and L${}_{4}0$ across
layers. Inset: Location of the selected points. Fushion zone boundary (FZB) of
the initial layer is also denoted. Figure 5:  Simulated profiles of (a1)-(a5)
von Mises stress $\sigma_{\mathrm{e}}$ and (b1)-(b5) effective plastic strain
$p_{\mathrm{e}}$ in the SS-processed four-layer powder bed with varying beam
power $P$ and scan speed $v$. Sectional profiles of (c1)-(c2) residual von
Mises stress $\sigma_{\mathrm{e}}$ and (d1)-(d2) accumulated plastic strain
$p_{\mathrm{e}}$ among SS-processed four-layer powder bed with (c1)-(d1)
$P=20~{}$\mathrm{W}$$ and $v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$
and (c2)-(d2) $P=30~{}$\mathrm{W}$$ and
$v=75~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$, respectively. Former
surfaces of each layer are denoted by discontinuous black lines, and fused
strut bottom is denoted with a solid white line. The Lamé’s stress ellipsoids,
representing the stress state at selected points P1-P6, are also illustrated
with the principle stresses denoted and colored (marine: tension; red:
compression). Figure 6: (a) Schematic of the mechanism of powder-resolved
mesoscopic residual stress formation, and (b) distorted strut captured in the
simulation, colormapped with the displacement along $z$-direction. These
schematics are based on the simulated microstructure with
$P=30~{}$\mathrm{W}$$ and $v=100~{}$\mathrm{m}\mathrm{m}~{}\mathrm{s}^{-1}$$.
Figure 7: (a) Segment-wise porosity of the sample produced under various
process parameters along (a1) building direction (BD) and (a2) scan direction
(SD), where the range of the substrate as well as the beam spot size
($D_{\mathrm{FWHM}}$) are denoted. Representative height ($H_{\mathrm{R}}$)
and width ($W_{\mathrm{R}}$) for porosity calculation are also selected.
(b1)-(b9) microstructures of SS-processed four-layer powder bed with varying
beam power and scan speed, which are marked as points in the porosity
$P\text{-}v$ map (c). The dotted lines represent the volumetric energy input
$U$ isolines and the dash-dotted line represents the median porosity isoline
($24.1\%$). (d) Phenomenological relation between densification factor
$\varrho$ (calculated from porosity) and $U$. Figure 8: $P\text{-}v$ maps of
(a) PB-averaged residual stress $\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and
(b) PB-averaged plastic strain $\bar{p}_{\mathrm{e}}^{\mathrm{P}}$. The dotted
lines maps represent the volumetric energy input $U$ isolines. $P\text{-}v$
pairs for simulated profiles in Fig. 5a and 5b are also denoted
correspondingly. Nonlinear regressions of (c)
$\bar{\sigma}_{\mathrm{e}}^{\mathrm{P}}$ and (d)
$\bar{p}_{\mathrm{e}}^{\mathrm{P}}$ on $U$ with the regression parameters
indicated correspondingly. Figure 9: Simulated profiles of (a1)-(a5) residual
von Mises stress $\sigma_{\mathrm{e}}$ and (b1)-(b5) effective plastic strain
$p_{\mathrm{e}}$ in the fused strut with varying beam power $P$ and scan speed
$v$, which are marked as points in the $P\text{-}v$ maps of (c) strut-averaged
residual stress $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and (d) strut-
averaged plastic strain $\bar{p}_{\mathrm{e}}^{\mathrm{S}}$, respectively. The
dotted lines maps represent the volumetric energy input $U$ isolines.
Structural distortion of each fused strut are also illustrated. Nonlinear
regressions of (e) $\bar{\sigma}_{\mathrm{e}}^{\mathrm{S}}$ and (f)
$\bar{p}_{\mathrm{e}}^{\mathrm{S}}$ on $U$ with the regression parameters
indicated correspondingly.

## References

  * Gu et al. [2021] D. Gu, X. Shi, R. Poprawe, D. L. Bourell, R. Setchi, J. Zhu, Science 372 (2021) eabg1487.
  * Ladani and Sadeghilaridjani [2021] L. Ladani, M. Sadeghilaridjani, Metals 11 (2021) 1391\. doi:10.3390/met11091391.
  * Zhang et al. [2017] Y. Zhang, L. Wu, X. Guo, S. Kane, Y. Deng, Y.-G. Jung, J.-H. Lee, J. Zhang, Journal of Materials Engineering and Performance 27 (2017) 1–13. doi:10.1007/s11665-017-2747-y.
  * Panwisawas et al. [2020] C. Panwisawas, Y. T. Tang, R. C. Reed, Nature communications 11 (2020) 1–4.
  * Körner et al. [2020] C. Körner, M. Markl, J. A. Koepf, Metallurgical and Materials Transactions A 51 (2020) 4970–4983.
  * Kim et al. [2020] W. R. Kim, G. B. Bang, O. Kwon, K.-H. Jung, H.-K. Park, G.-H. Kim, H.-T. Jeong, H. G. Kim, Materials & Design 195 (2020) 109035.
  * Yang et al. [2019] Y. Yang, O. Ragnvaldsen, Y. Bai, M. Yi, B.-X. Xu, npj Comput. Mater. 5 (2019) 1–12.
  * Gu et al. [2012] D. Gu, W. Meiners, K. Wissenbach, R. Poprawe, International materials reviews 57 (2012) 133–164.
  * Li et al. [2018] J. Q. Li, T. H. Fan, T. Taniguchi, B. Zhang, International Journal of Heat and Mass Transfer 117 (2018) 412–424. doi:10.1016/j.ijheatmasstransfer.2017.10.001.
  * Li et al. [2010] R. Li, J. Liu, Y. Shi, M. Du, Z. Xie, Journal of Materials Engineering and Performance 19 (2010) 666–671.
  * Shuai et al. [2019] C. Shuai, G. Liu, Y. Yang, F. Qi, S. Peng, W. Yang, Z. Liu, Ceramics International 45 (2019) 20234–20242.
  * Duan et al. [2010] B. Duan, M. Wang, W. Y. Zhou, W. L. Cheung, Z. Y. Li, W. W. Lu, Acta biomaterialia 6 (2010) 4495–4505.
  * Eshraghi and Das [2010] S. Eshraghi, S. Das, Acta biomaterialia 6 (2010) 2467–2476.
  * Mapley et al. [2019] M. Mapley, J. P. Pauls, G. Tansley, A. Busch, S. D. Gregory, Scripta Materialia 172 (2019) 154–158.
  * Wendhausen et al. [2017] P. Wendhausen, C. Ahrens, A. Baldissera, P. Pavez, J. Mascheroni, in: 2017 IEEE International Magnetics Conference (INTERMAG), IEEE, pp. 1–1.
  * Sing et al. [2017] S. L. Sing, W. Y. Yeong, F. E. Wiria, B. Y. Tay, Z. Zhao, L. Zhao, Z. Tian, S. Yang, Rapid Prototyping Journal 23 (2017) 611–623.
  * Jhong and Lee [2016] K. J. Jhong, W. H. Lee, in: Industrial Technology (ICIT), 2016 IEEE International Conference on, IEEE, pp. 1115–1118.
  * Zhang et al. [2022] B. Zhang, C. Guo, X. Cao, X. Yuan, X. Li, H. Huang, S. Dong, Z. L. Wang, K. Ren, Nano Energy 104 (2022) 107897. doi:10.1016/j.nanoen.2022.107897.
  * Zhang et al. [2021] M. Zhang, D. Yan, J. Wang, L.-H. Shao, Journal of the Mechanics and Physics of Solids 151 (2021) 104396. doi:10.1016/j.jmps.2021.104396.
  * Yan et al. [2023] D. Yan, J. Wang, J. Xiang, Y. Xing, L.-H. Shao, Science Advances 9 (2023). doi:10.1126/sciadv.adc8845.
  * Chen et al. [2022] S.-g. Chen, H.-j. Gao, Q. Wu, Z.-h. Gao, X. Zhou, et al., Journal of Materials Research and Technology 17 (2022) 2950–2974.
  * Mercelis and Kruth [2006] P. Mercelis, J.-P. Kruth, Rapid prototyping journal 12 (2006) 254–265.
  * Roberts et al. [2010] I. Roberts, C. Wang, K. Kibble, M. Stanford, D. Mynors, in: Proceedings of the 36th International MATADOR Conference, Springer, pp. 535–538.
  * Pant et al. [2021] P. Pant, S. Sjöström, K. Simonsson, J. Moverare, S. Proper, S. Hosseini, V. Luzin, R. Peng, Metals 11 (2021) 861\.
  * Mukherjee et al. [2017] T. Mukherjee, W. Zhang, T. DebRoy, Computational Materials Science 126 (2017) 360–372.
  * Ibraheem et al. [2002] A. K. Ibraheem, B. Derby, P. J. Withers, MRS Online Proceedings Library (OPL) 758 (2002) LL1–8.
  * Simson et al. [2017] T. Simson, A. Emmel, A. Dwars, J. Böhm, Addit. Manuf. 17 (2017) 183–189.
  * Takezawa et al. [2022] A. Takezawa, H. Guo, R. Kobayashi, Q. Chen, A. C. To, Addit. Manuf. 60 (2022) 103194.
  * Yi et al. [2019] M. Yi, B.-X. Xu, O. Gutfleisch, Computational Mechanics 64 (2019) 917–935.
  * Hussein et al. [2013] A. Hussein, L. Hao, C. Yan, R. Everson, Materials & Design (1980-2015) 52 (2013) 638–647.
  * Parry et al. [2016] L. Parry, I. Ashcroft, R. D. Wildman, Additive Manufacturing 12 (2016) 1–15.
  * Vastola et al. [2016] G. Vastola, G. Zhang, Q. Pei, Y.-W. Zhang, Additive Manufacturing 12 (2016) 231–239.
  * Chiumenti et al. [2017] M. Chiumenti, E. Neiva, E. Salsi, M. Cervera, S. Badia, J. Moya, Z. Chen, C. Lee, C. Davies, Additive Manufacturing 18 (2017) 171–185.
  * Ganeriwala et al. [2019] R. Ganeriwala, M. Strantza, W. King, B. Clausen, T. Q. Phan, L. E. Levine, D. W. Brown, N. Hodge, Additive Manufacturing 27 (2019) 489–502.
  * Carraturo et al. [2020] M. Carraturo, J. Jomo, S. Kollmannsberger, A. Reali, F. Auricchio, E. Rank, Additive Manufacturing 36 (2020) 101498.
  * Zhou et al. [2021] X. Zhou, Y. Yang, S. Bharech, B. Lin, J. Schröder, B.-X. Xu, GAMM-Mitteilungen 44 (2021) e202100017.
  * Yang et al. [2023] Y. Yang, T. D. Oyedeji, X. Zhou, K. Albe, B.-X. Xu, npj Computational Materials 9 (2023) 103.
  * Panwisawas et al. [2015] C. Panwisawas, C. Qiu, Y. Sovani, J. Brooks, M. Attallah, H. Basoalto, Scripta Materialia 105 (2015) 14–17.
  * Arısoy et al. [2019] Y. M. Arısoy, L. E. Criales, T. Özel, Optics & Laser Technology 109 (2019) 278–292.
  * Zhang et al. [2022] Z. Zhang, Y. Wang, P. Ge, T. Wu, Coatings 12 (2022) 1277. doi:10.3390/coatings12091277.
  * Yang et al. [2019] Y. Yang, O. Ragnvaldsen, Y. Bai, M. Yi, B.-X. Xu, npj Computational Materials 5 (2019). doi:10.1038/s41524-019-0219-7.
  * Yang et al. [2020] Y. Yang, T. D. Oyedeji, P. Kühn, B.-X. Xu, Scripta Materialia 186 (2020) 152–157.
  * Oyedeji et al. [2023] T. D. Oyedeji, Y. Yang, H. Egger, B.-X. Xu, Physical Review E 108 (2023) 025301.
  * Moelans et al. [2008] N. Moelans, B. Blanpain, P. Wollants, Physical Review B 78 (2008) 024113.
  * Turnbull [1951] D. Turnbull, JOM 3 (1951) 661–665.
  * Yang et al. [2022] Y. Yang, M. Fathidoost, T. D. Oyedeji, P. Bondi, X. Zhou, H. Egger, B.-X. Xu, Scr. Mater. 212 (2022) 114537.
  * Nicoli et al. [2011] M. Nicoli, M. Plapp, H. Henry, Phys. Rev. E 84 (2011) 046707.
  * Denlinger et al. [2016] E. R. Denlinger, V. Jagdale, G. Srinivasan, T. El-Wardany, P. Michaleris, Addit. Manuf. 11 (2016) 7–15.
  * Gusarov et al. [2009] A. V. Gusarov, I. Yadroitsev, P. Bertrand, I. Smurov, Journal of Heat Transfer 131 (2009) 072101.
  * Voigt [1889] W. Voigt, Annalen der Physik 274 (1889) 573–587.
  * Schneider et al. [2015] D. Schneider, O. Tschukin, A. Choudhury, M. Selzer, T. Böhlke, B. Nestler, Computational Mechanics 55 (2015) 887–901.
  * Durga et al. [2015] A. Durga, P. Wollants, N. Moelans, Computational Materials Science 99 (2015) 81–95.
  * Permann et al. [2020] C. J. Permann, D. R. Gaston, D. Andrvs, R. W. Carlsen, F. Kong, A. D. Lindsay, J. M. Miller, J. W. Peterson, A. E. Slaughter, R. H. Stogner, R. C. Martineau, SoftwareX 11 (2020) 100430.
  * Balay et al. [2019] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al. (2019).
  * Dunne and Petrinic [2005] F. Dunne, N. Petrinic, Introduction to Computational Plasticity, Oxford University Press on Demand, 2005.
  * Simo and Hughes [2006] J. C. Simo, T. J. Hughes, Computational inelasticity, volume 7, Springer Science & Business Media, 2006\.
  * Yang et al. [2022] Y. Yang, A. Großmann, P. Kühn, J. Mölleney, L. Kropholler, C. Mittelstedt, B.-X. Xu, Journal of Materials Processing Technology 299 (2022) 117316\.
  * Simchi [2006] A. Simchi, Materials Science and Engineering: A 428 (2006) 148–158.
  * Desmond and Weeks [2014] K. W. Desmond, E. R. Weeks, Physical Review E 90 (2014) 022204.
  * Price et al. [1964] A. Price, H. Holl, A. Greenough, Acta metallurgica 12 (1964) 49–58.
  * Blakely and Mykura [1963] J. Blakely, H. Mykura, Acta Metallurgica 11 (1963) 399–404.
  * Mead and Birchenall [1956] H. Mead, C. Birchenall, JOM 8 (1956) 1336–1339.
  * Liu et al. [2012] F. Liu, Q. Zhang, W. Zhou, J. Zhao, J. Chen, Journal of Materials Processing Technology 212 (2012) 2058–2065.
  * Hoshino et al. [1986] T. Hoshino, K. Mito, A. Nagashima, M. Miyata, International Journal of Thermophysics 7 (1986) 647–662.
  * Chase Jr [1998] M. W. Chase Jr, J. Phys. Chem. Ref. Data, Monograph 9 (1998).
  * Di Schino et al. [2002] A. Di Schino, J. Kenny, G. Abbruzzese, Journal of materials science 37 (2002) 5291–5298. doi:10.1023/A:1021068806598.
  * Mills [2002] K. C. Mills, Recommended values of thermophysical properties for selected commercial alloys, Woodhead Publishing, 2002.
  * Khairallah and Anderson [2014] S. A. Khairallah, A. Anderson, Journal of Materials Processing Technology 214 (2014) 2627–2636.