A review of constitutive models used in macroscale finite element analysis of additive manufacturing and post-processing of additively manufactured components

ABSTRACT

The development of fast and reliable finite element (FE) simulation codes has increased the relevance of FE simulations in predicting the outcomes of additive manufacturing (AM) processes and in the post-processing of additively manufactured components. The accuracy of FE simulations relies heavily on the quality of inputs such as material data, physical properties, numerical settings, and simulation methods. We review models used for simulating AM processes and the post-processing of AM components based on macro-scale continuum FE simulations. Particular attention is given to phase transformation and microstructure-based models because they determine the properties of AM components along manufacturing process chains and their subsequent performance. However, it's noted that FE simulations of AM components usually only consider the post-processing of AM components rather than complete process chains that include both AM and post-processing. A process chain simulation framework is proposed to model the entire manufacturing chain of AM parts.

KEYWORDS:

• Additive manufacturing
• material models
• phase transformation
• process chain simulation

Abbreviations
FE	=	Finite element
AM	=	Additive manufacturing
ASTM	=	American society for testing and materials
PBF	=	Powder bed fusion
LPBF	=	Laser powder bed fusion (Selective laser melting)
ICME	=	Integrated computational material engineering
DED	=	Direct energy deposition
L-DED	=	Laser based DED
WA-DED	=	Wire Arc DED (wire arc additive manufacturing)
BJT	=	Binder jetting
GTAW	=	Gas tungsten arc welding
GMAW	=	Gas metal arc welding
HAZ	=	Heat-affected zone
JMAK	=	Johnson-Mehl-Avrami-Kolmogrov
JMA	=	Johnson-Mehl-Avrami
KM	=	Koistinen–Marburger
$\alpha$	=	Alpha phase
${\alpha }^{\mathrm{\prime }}$	=	Alpha prime phase
$\text{β}$	=	Beta phase
EDM	=	Electrical discharge machining
DRX	=	Dynamic recrystallization

1. Introduction

Additive manufacturing (AM) or 3D printing refers to a group of manufacturing processes that construct 3D components layer by layer from 2D prints. AM is a vital technology required to implement Industry 4.0 because it provides unparalleled design freedom and ease of process automation. Advances in AM have considerably enhanced process design by lowering costs and speeding up prototype generation. AM can produce parts with complex geometries, thereby replacing multiple components manufactured by traditional techniques, reducing lead time, and simplifying assembly. Manufacturing costs can be reduced in hot forging by printing near-net-shape preforms to ensure a single hit of components during forging, thereby lowering the number of dies in multi-stage forging. Furthermore, AM can replace metal casting in applications that require low porosity, short lead times, and complex geometries [1]. AM is expected to reach an estimated market value of $50 billion between 2029 and 2031 [2]. AM and conventional manufacturing processes can be combined to expand their applicability domain (hybrid-metal AM). Pragana et al. [3] conducted a review of hybrid-metal AM. They demonstrated that hybrid process routes increase flexibility, and the limitations of AM are mitigated by combining AM and traditional manufacturing operations. According to the ASTM F42 standard, the most prevalent methods for AM of metals are BJT, DED, and PBF [4].

Despite the enormous interest in AM technology, the reliability of AM components is generally low due to production flaws that can result in parts with inferior strength and fatigue performance compared to wrought components. During the solidification of the component, high-temperature gradients may prevail because of the re-melting of the deposited material and rapid cooling of some portions of the part. This leads to the formation of unusual microstructures and mechanical properties in AM components compared to components manufactured by traditional methods [4]. Fan et al. [5] highlighted the dependency of AM part microstructures on the cooling rate from molten temperature. Controlling the cooling rate is essential to achieving the desired microstructures, which directly influence the mechanical properties and performance of the final components. AM processes also introduce residual stresses and distortions in components [6]. High tensile residual stresses promote crack formation, especially in corrosive environments, and reduce fatigue resistance [7]. Multiple thermal cycles and the impracticality of directly measuring the temperature and microstructure evolution in AM experiments make experimental methods very expensive and time-consuming for AM process optimisation. For instance, Yadroitsev [8] stated that about 130 parameters affect LPBF, and about 20 are first-order parameters. This complexity makes experimental optimisation time-consuming and costly. The microstructures of AM components are usually altered in post-processing operations. As-built AM components usually have unique microstructures characterised by material discontinuities, highly textured columnar grains, complex phases, and compositional variation [9]. AM operations are usually part of process chains that include other manufacturing processes. Subsequent manufacturing operations are conducted on AM components to improve their mechanical properties, achieve grain refinement, close voids in components, reduce anisotropy, and minimise geometrical defects.

FE simulations are a low-cost method to understand the underlying physical processes in AM. FE simulations enable estimation of process outcomes before actual manufacturing, thus reducing the need for extensive trial designs and facilitating tailored experimental procedures. The outputs from FE simulations can be calibrated to mechanical, microstructural, and thermal properties obtained from experiments. To advance the adoption of FE simulation-based product development of AM components in the manufacturing industry, it is imperative to integrate microstructure evolution and phase change phenomena across entire manufacturing process chain simulations. According to existing literature, this integration entails performing FE simulations across three domains: process, material, and product [10]. FE modelling of process chains is performed by combining process models within an integrated computational materials engineering (ICME) framework. ICME integrates multi-scale (micro, meso, and macro) simulation results as inputs to subsequent simulations in a process chain, as shown in Figure 1. The goal of ICME is to create holistic systems that will accelerate material development, transform the engineering design optimisation process, and unify design and manufacturing [11]. ICME focuses on determining the properties of finished components while considering property alterations in the parts’ manufacturing history [12]. The processes in the simulation chain in Figure 1 can be simulated at all scales. Micro- and mesoscale simulations are time consuming at the part scale. To better model details on the microstructure, AM and surface hardening can be modelled at macro-, meso-, and micro-scales, while annealing is modelled using macro- and meso-scale simulations.

Figure 1. Illustration of ICME paths for simultaneous materials and process design of additively manufactured components.

1.1. Multiscale modelling of AM processes

Metal additive manufacturing processes are challenging to model because of the occurrence of multiple physical phenomena and the multiscale nature of these phenomena. In addition, coupling the interactions of mechanical properties and microstructure evolution to temperature change is very demanding due to the transient nature of AM [13]. To address these challenges, AM is modelled at micro, meso, and macro scales based on the domain of interest, as shown in Figure 2. The information from micro- and mesoscale models can be incorporated into the macroscale. Francois et al. [14] presented a review of AM simulation at different lengths and time scales and highlighted the challenges and opportunities in multiscale AM simulation. Similarly, Gatsos et al. [15] reviewed process-microstructure-property relationships in metal AM processes at different FE simulation scales and proposed creating a feedback-based approach to link the simulations and experiments.

Figure 2. Different scales of study in modelling metal additive manufacturing processes (top) and the fundamental physical phenomena occurring in and around the melt pool region at the scale of the powder (bottom) [16].

Microscale models.

Microscale simulations are essential for modelling physical and material properties. However, executing micro-scale simulations at the part scale is extremely time-consuming, which limits their practicality in industry settings. They can replicate properties such as surface tension anisotropies or kinematic coefficients in the rapid solidification regime [17]. At the solid–liquid interface, phase-field (PF) models are used to predict detailed morphologies and the dynamics of microstructure solidification [18]. PF models shed light on experimentally observed microstructures, produce quantitative predictions, and reveal the local grain orientation. The crystal structures of the phases are assigned Gibbs energy models and can be performed on CALPHAD and cellular automata. Monte Carlo models have also been employed for microstructure prediction [19]. Crystal plasticity FE models can link the results from PF and Monte Carlo micromodels to mesocale FE models [20].

Mesoscale models.

Mesoscale models are powerful at capturing the molten pool dynamics as a function of the powder, free surface flows, describing phase transformation, and determining the effects of processes and material parameters [21]. Mesoscale simulations can be calibrated to link the physical characteristics observed in micromodels to macroscale simulations [22]. The mesoscale model can be grouped as follows:

• Computational fluid dynamics (CFD) models: CFD models are good at predicting the thermal-fluid flow and geometry in AM; however, they are poor at capturing the structural response. As a result, coupling thermal-fluid results from CFD solvers to FEM to obtain the structural response has been explored. Employing the hybrid simulation method, Chen and Yan [23] proposed a high-fidelity model to simulate GMAW. Similarly, Mathews et al. [24] simulated L-DED of SS 316L and validated the residual stresses.

• Smoothed Particle Hydrodynamics (SPH) models: A recent review on numerical simulation of metal AM based on particle methods highlights the advantages of SPH methods, including the potential of particle methods to achieve realistic powder morphology, their strength at modelling AM of complex geometries, and their straight-forward implementation to solve multi-physics problems [16]. The research group also developed new algorithms based on a neighbour-search scheme [25] and dynamic refinement-coarsening [26] to speed up SPH simulations. However, SPH simulations suffer from numerical instability, poor surface representation, and boundary treatment.

• Lattice Boltzmann methods (LBM) models: LBM is the only method capable of modelling all the physical processes observed in the molten pool and is fully explicit, which offers massive parallelisation options [27]. LBM can model phenomena up to the mesoscale.

The review was conducted on models used in macro-scale (part-scale) FEM because they save time needed to model detailed phenomena in large industrial parts. Bayat et al. [28] stated that, given today’s computing capabilities and existing modelling approaches, it is currently impossible to simulate all relevant physical phenomena, including stress evolution and fluid flow, in macroscale simulations. Currently, it is only possible to simulate detailed phenomena for a few tracks at the mesoscale [26]. Additionally, manufacturing processes are optimised by part-scale simulations. Although the macroscale finite volume method (FVM) has proven to be faster at performing part scale AM simulations than FEM [29], FVM is poor at handling problems involving structural responses often present in post-AM processes.

1.2. Aim of the review

Currently, few studies have reported simulations of AM and the post-AM processing of components in a process chain. Complete process chain simulations are required for process optimisation, enhancing mechanical properties, and extending the service life of AM components. This can be achieved by phase change and microstructure-based simulation at each stage in process chains. Such an approach requires material models that are valid at each stage in a process chain. Thus, it is necessary to review material models used for AM and subsequent processes at a part-scale to establish frameworks for process chain simulation.

The aim of this review is fourfold: (1) To provide an overview of the models used in macro-scale FE simulations of AM processes. (2) To outline models used in FE simulation of the post-processing of AM components. (3) To review studies with process chain simulations. (4) To propose a framework to implement the processing-structure–property-performance of AM components and achieve AM process chain optimisation.

The review is structured as follows: The theoretical background and models used in FE simulations of AM processes are presented in Chapter 2. Chapter 3 presents models for simulating and evaluating post-AM components. Chapter 4 introduces studies with process chain simulations. In Chapter 5, we present a discussion and outlook on process chain simulation. Chapter 6 provides conclusions.

2. Macroscale finite element analysis of AM

2.1. Theoretical background of macroscale continuum AM simulation

Macroscale AM simulations can be modelled based on continuum mechanics by differential equations governing the conservation of mass, internal energy, momentum and current continuity [30]. (1) ${\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(\mathrm{\rho v}\right)=0}$(1) (2) $\rho \dot{w}=\sigma :D+\dot{q}-\mathrm{\nabla q}$(2) (3) $\rho \left({\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}+v\cdot \mathrm{\nabla v}}\right)=-{\mathrm{\nabla }}_{\mathit{p}}+\mu {\mathrm{\nabla }}^{2}v+{\displaystyle \frac{1}{3}\mu \mathrm{\nabla }\left(\mathrm{\nabla }\cdot v\right)+\mathrm{\rho f}}$(3) (4) ${\mathrm{\nabla }}^2\varphi +\sum =0$(4)

From the first law of thermodynamics, the general thermal energy balance in a closed system is determined by (5) $\dot{Q}=Q_{\mathrm{cd}}+Q_{\mathrm{cv}}+Q_r$(5) where $\dot{Q}$, $Q_{\mathrm{cd}}$, $Q_{\mathrm{cv}}$, and $Q_r$ represent the rate of heat generation, conduction, convection and radiation, respectively. For an isotropic volume, the heat conduction equation can be defined by Fourier’s law: (6) $\alpha \left({\displaystyle \frac{{\mathrm{\partial }}^{2}T}{\mathrm{\partial }{x}^2}+{\displaystyle \frac{{\mathrm{\partial }}^{2}T}{\mathrm{\partial }{y}^2}+{\displaystyle \frac{{\mathrm{\partial }}^{2}T}{\mathrm{\partial }{z}^2}}}}\right)=\rho C_p{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}}$(6) where $\alpha$ is the thermal conductivity, $\rho$ is the density, $C_p$ is the specific heat capacity, $T$ is the temperature, and $t$ is the time. The thermal boundary conditions for convection and radiation are calculated at the part surface by (7) $Q_{\mathrm{cv}}=h_{\mathrm{conc}}\left(T_s-T_0\right)$(7) (8) $Q_r=\mathrm{ϵ\sigma }\left[{\left(T_s-T_{\mathrm{abs}}\right)}^4-{\left(T_0-T_{\mathrm{abs}}\right)}^4\right]$(8) where $h_{\mathrm{conc}}$ refers to the heat transfer coefficient, $ϵ$ is the emissivity, $\sigma$ represents the Stefan–Boltzmann constant. $T_s$, $T_0$ and $T_{\mathrm{abs}}$ are the surface temperature, surrounding temperature, and absolute zero temperature respectively. The FE simulation programme works by solving the heat transfer problem: (9) $k{\mathrm{\nabla }}^{2}T-\rho C_p{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\dot{Q}=0}$(9) where $k{\mathrm{\nabla }}^2$ is thermal conduction, with an assumption of quasi-static equilibrium at the current time step. (10) $\mathrm{\nabla }\cdot \sigma =0$(10) where $\sigma$ refers to the stress. In Equation (9), convection is assumed to be zero for a solid body and heat loss by radiation is omitted. The constitutive law to solve the mechanical problem is defined by: (11) $\sigma =\mathrm{Cϵ}$(11) (12) $ϵ={ϵ}_e+{ϵ}_p+{ϵ}_t+{ϵ}_{\mathrm{tp}}$(12) where $C$ is a the forth order stiffness matrix tensor, $ϵ$, ${ϵ}_e$, ${ϵ}_p$, ${ϵ}_t$, and ${ϵ}_{\mathrm{tp}}$ represent the total strain, elastic strain, plastic strain, thermal strain, and strain due to phase transformation, respectively. The elastic strain is defined by ${ϵ}_e={\displaystyle \frac{\sigma }{E}}$. The plastic strain is calculated by ${ϵ}_p=\dot{\lambda }{\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\sigma }}$. Where $\dot{\lambda }$ is the plastic strain rate and $f$ is the yield function. The direction for plastic strain increment is given by ${\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\sigma }}$. The thermal strain is defined by ${ϵ}_t={\int }_{T_{\mathrm{ref}}}^T{\alpha }_T\left(T\right)\mathrm{dT}$. where $T_{\mathrm{ref}}$ is the temperature and ${\alpha }_T$ is the coefficient of thermal expansion. The phase transformation strain is determined based on the transformation law in section 2.6. The strains and displacements obey the following: (13) ${ϵ}_{\mathrm{ij}}={\displaystyle \frac{1}{2}\left[{\displaystyle \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{u}_j}+{\displaystyle \frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{u}_i}}}\right]}$(13) The continuum assumption erases the molecular discontinuities by averaging the microscopic quantities in small sampling volumes. The representative volume of the component is assumed to be large enough to erase the microscopic fluctuations but also small enough to point out the macroscopic variations in quantities of interest. All the macroscopic quantities (stresses, strains, temperature, etc.) are assumed to vary continuously from point to point within the bodies. The continuum assumption inherently limits the ability to simulate details of the process at lower scale lengths but provide the average values of quantities over the volume.

2.2. An overview on macro-scale FE simulation

The basics of FE modelling of AM processes draw upon the works of Rosenthal and Goldak et al. on modelling of welding operations that include the energy source, filler material, and build structure [31,32]. Owing to the similarities between multi-pass welding and AM, AM is simulated like multi-pass welding, and data from welding simulations is often used in AM simulations. The temperature is highest in the molten region where the material is deposited. Heat transfers from the newly deposited material to the previous layers mainly by conduction and causes thermal cycles in the previously deposited layers. High temperatures cause changes in the mechanical and metallurgical properties of the molten area and heat-affected zone (HAZ). Figure 3 presents a flowchart summarising AM simulation process. The algorithm works by solving the thermal and mechanical problems with dedicated solvers. Apart from the heat source model, the rest of the algorithm is required for post AM simulations.

Figure 3. Flowchart of a thermo-metallurgical-mechanical coupling model for the SLM process [33].

In the rest of chapter 2, thermal, mechanical, phase transformation, and microstructural models used in the macroscale FE simulation of AM are presented. The effects of model parameters on FE simulation, advanced modelling techniques, and the interactions between parameters of different phases are highlighted.

2.3. Thermal features

2.3.1. Heat source models

The Goldak (volumetric heat source) model is often used to describe the heat source in AM and welding processes for macro-scale continuum FE simulations [32]. The Goldak is usually used for AM modelling due to its capacity to estimate the molten pool shape. The Goldak model describes the front and rear parts of the ellipsoid formed by the moving heat source. (14) $\dot{Q}\left(x,y,z,t\right)={\displaystyle \frac{6\sqrt{3}{f}_{f,r}Q}{\mathrm{ab}{c}_{f,r}\pi \sqrt{\pi }}{e}^{-3{\left(\frac{x}{w}\right)}^2}{e}^{-3{\left(\frac{y}{d}\right)}^2}{e}^{-3{\left(\frac{z+v\left(\tau -t\right)}{c_{f,r}}\right)}^2}}$(14) where $Q$ is the total induced power, $f_f$ and $f_r$ are the portions of the heat generated in the frontal and rear parts of the double ellipsoid ($f_1+f_2=2$). $a,b,C_f$ and $C_r$ are geometric parameters that define the size and shape of the ellipses. The weld pool energy is determined by $\dot{Q}=\mathrm{\mu VI}$ where $\mu$ is the efficiency of the heat source, $V$ and $I$ are the voltage and current of the heat source. Recent reviews present the state of the art on heat source models for macroscale continuum FE simulations [34,35]. Such heat source models do not reveal process details; however, they are the basis for calculating the temperature, residual stresses, and fusion zone shape and size. Continuum FE simulations of multiple layers and hatches are ‘ … critical to understanding the structure and properties of the entire build … ’ in AM [20].

2.3.2. Coefficient of thermal expansion

Heat input in AM processes causes a temperature rise that leads to thermal expansion. AM parts contract upon cooling. The temperature drop causes residual stresses to develop and distort components. The coefficient of thermal expansion defines the increase in length of a specimen per unit temperature increase. The coefficient of thermal expansion increases with temperature [36,37]. Blumm and Henderson [38] measured the volumetric coefficient of thermal expansion of alloy 718 in the temperature range $0<T<{1600}^{\circ }\mathrm{C}$ in a dilatometer test. The volumetric thermal coefficient was transformed into a linear coefficient. At high temperatures, the thermal expansion coefficient of the material is considered zero because the material is assumed to be in a liquid state and no straining occurs. Zhang et al. [39] modified the values of the thermal expansion coefficient from the ASTM standards to obtain the distortion measured in the welding experiments of two AF1410 steel plates. Perić et al. [40] demonstrated that the coefficient of thermal expansion has a significant influence on part deflection in AM simulation. In the same light, Park et al. [41] showed that distortion increased with an increase in the coefficient of thermal expansion using experimental and numerical analyses in the WA-DED of a group of 8 metals.

Some studies consider the influence of multiple phases on the coefficient of thermal expansion. Zhang et al. [6] implemented multiple thermal expansion coefficients in modelling the transformation plasticity of multi-phase AF1410 steel. Babu et al. [42] calculated the average coefficient of thermal expansion of the individual phases using the rule of mixture. Caron et al. [43] revealed that the residual stresses in FE simulations depend on the portion of formed martensite.

2.3.3. Thermal conductivity

Thermal conduction is the process of heat transfer from the hotter to the colder ends of a material and is quantified by thermal conductivity. The thermal conductivity of bulk materials is significantly higher than that of powder materials due to the close arrangement of particles. Zhu and Chao [44] proved that thermal conductivity significantly affects the transient temperature distribution. The thermal conductivity gradually increases with temperature up to the solid-to-liquid transition temperature and is followed by a sudden increase at temperatures where the material is completely liquid [36,45]. The high values at the melting point are explained by the effect of convection stirring in the molten pool [46]. Draxler et al. [47] determined the conductivity of the mushy zone by extrapolating the thermal conductivity of the metal and liquid phases using a linear mixture rule. (15) $k_{\mathrm{mush}}=k_s{f}_s+k_l\left(1-f_s\right)$(15) where $f_s$ is the volume fraction of the solid, $k_s$ and $k_l$ are extrapolated thermal conductivities of the solid and liquid phases. It was assumed that the temperature of the mushy zone lies between the temperature of the solid and liquid phases. Zhang et al. [6] modelled the transformation plasticity of multi-phase AF1410 steel using the multiple thermal conductivities of different phases. In some studies, the thermal conductivity reduces with an increase in temperature until the solid-to-liquid transformation, before a sudden rise when the material is fully liquid [36,48,49]. Ding et al. [50] described a heat conduction equation that considers the direction of thermal conductivity. Perić et al. [40] demonstrated that thermal conductivity influences the heat-affected zone (HAZ) in welding simulations. High thermal conductivities lead to large HAZs, and vice versa.

2.3.4. Thermal radiation and convection

The freshly deposited material loses heat to the environment through radiation and convection during cooling of the AM component. The thermal problem becomes non-linear because of the high temperatures, so it's critical to have a proper estimate of the thermal radiation and convection at the component surface [51]. Thermal radiation and convection increase with temperature [37]. Tran et al. [52] performed a parametric study to determine an accurate convection coefficient (natural convection) in the laser cladding of Ti-6Al-4 V. Furthermore, they adjusted the critical phase transformation temperature to obtain the fusion zone and HAZ measured by optical microscopy.

2.4. Mechanical features

2.4.1. Elastic modulus

The high temperatures in AM cause a reduction in the force required to achieve deformation. As a result, the elastic modulus decreases with increasing temperature. The elastic modulus is defined at specific temperatures, ranging from room temperature to molten temperature. Some studies express the elastic modulus as a function of temperature [42,47]. Based on experiments conducted by Outinen and Mäkeläinen [53], Perić et al. [40] developed formulas to describe the evolution of the elastic modulus with temperature. Two equations were used to describe the elastic modulus above and below 600°C. Similarly, Draxler et al. [47] proposed an isotropic ‘soft’ material model to determine the elastic modulus when the metal turns to a liquid state. (16) $E_l=3K_l\left(1-2v\right)$(16) where $v$ is the Poison’s ratio and $K_l$ refers to the bulk modulus. The elastic modulus of the mushy zone was approximated using a linear mixture rule for temperatures above the coherent temperature ($T_c$ = 1278°C). The mushy zone was assumed to have the same elastic modulus as the liquid state at temperatures above $T_c$ since the mush cannot transmit tensile loads above the coherent temperature. The elastic modulus of the mush zone was written as (17) $E_{\mathrm{mush}}=\{\begin{array}{c}{E}_s{f}_s+\left({\displaystyle \frac{E_l-E_s\left(T_c\right)f_s\left(T_c\right)}{1-f_s\left(T_c\right)}}\right)\left(1-f_s\right),T\le T_c\\ E_l,T>T_c\end{array}$(17) where $E_s$ is the elastic modulus of the solid phase extrapolated to the mushy zone and $f_s$ is the volume fraction of the solid. Zhang et al. [6] included the effect of transformation plasticity in the welding of AF410 steel to improve the prediction of the residual stresses and distortion. The linear mixture rule was used to model the elastic moduli of the phases. (18) $X_{\mathrm{mix}}=\sum _1^4{X}_k{Z}_k\sum _1^4{Z}_k=1,k=1,2,3,4$(18) $k$  = 1,2,3,4 are phase fractions of the initial state, austenite, bainite and martensite.

2.4.2. Yield stress

Similar to the elastic modulus, the yield stress reduces with an increase in temperature as lower forces are required to deform materials at high temperatures [54,55]. Park et al. [41] showed that distortion moderately increased with an increase in elastic modulus and yield stress. Chiumenti et al. [36] coupled a temperature-dependent J2-yield surface and temperature evolution in multi-pass welding of nickel super-alloy 718. The yield surface contracted with an increase in temperature. Zhang et al. [56] showed that the dependence of the yield strength for metallic materials on temperature and strain is described by: (19) ${\sigma }_y\left(T,\dot{ϵ}\right)=\sigma \left(T_0,{\dot{ϵ}}_0\right){\left[{\displaystyle \frac{T_m-T}{T_m-T_0}\times {\displaystyle \frac{E_T}{E_{T_0}}\times {\displaystyle \frac{1+v_{T_0}}{1+v_T}}}}\right]}^{\frac{1}{2}}{\left({\displaystyle \frac{T_m-T}{T_m-T_0}}\right)}^n\left[1+{\left({\displaystyle \frac{\dot{ϵ}}{D}}\right)}^{1/p}\right]$(19) where $\sigma \left(T_0,{\dot{ϵ}}_0\right)$ refers to the reference yield strength, at room temperature $T_0$ and strain rate ${\dot{ϵ}}_0$. The melt, current, and reference temperatures are represented by $T_m$, $T$, and $T_0$, respectively. $E_{T_0}$ and $E_T$ are the elastic modulus at temperatures $T$ and $T_0$. $v_{T_0}$ and $v_T$ are the Poisson’s ratio at temperatures $T$ and $T_0$. The current strain rate is defined by $\dot{ϵ}$. $n$, $D$, and $p$ are material parameters obtained from fitting experimental data. Ali and Han [57] implemented the law of linear weight to obtain the average yield strength of individual phases in WA-DED of EH36 steel. Zhang et al. [6] included the yield strength of different phase fractions.

2.4.3. Flow stress

In AM simulations, the flow stress is usually defined till the yield point to make the analysis easier. This is because little plastic deformation is thought to happen in AM, and part distortion is thought to be mostly caused by elastic springback [37,58]. The plastic part of the flow stress and yield functions is included in AM FE simulations to improve the accuracy of predicted residual stresses and distortion [59]. Manurung et al. [60] used mechanical parameters determined by JMATpro software to predict distortion in the WA-DED of SS3I6L with good accuracy. Deng et al. [61] performed FE simulations to calculate the residual stresses of a multi-pass butt-welded joint using several material laws (perfect plastic, isotropic strain hardening, kinematic strain hardening, and mixed isotropic-kinematic strain hardening). The residual stresses were influenced by the yield strength and the strain-hardening model. Temperature-dependent flow stresses were implemented by Malmelöv et al. [62] to introduce the temperature dependence of mechanical properties. Chen et al. [63] used the Prager-Ziegler kinematic hardening model and the Von-Mises yield criterion in modelling the yield behaviour of different phases in the powder L-DED of a low-temperature transformation pre-alloyed powder.

2.4.4. Strain rate sensitivity

The strain rate sensitivity is known to be profoundly influenced by temperature [64]. Gkatzogiannis [65] considered the strain rate dependency in welding simulations. Although not as high as the strain rates observed in crash analyses$\left(\dot{ϵ}\gg 100\right)$, the strain rate deviates from the static case $\left(\dot{ϵ}\to 0\right)$. Strain rates approaching $\dot{ϵ}>0.1$ were observed at the interface of the heat-affected zone (HAZ) and the fusion zone in the welding of structural steel S355. The strain rate dependency of the yield stress in the structural analysis was determined using the model proposed by Perzyna [66]. (20) $\dot{\widehat{{ϵ}_{\mathrm{pl}}}}=\gamma {\left({\displaystyle \frac{\sigma }{{\sigma }_0}-1}\right)}^{1/\mathit{m}}$(20) where $\dot{\widehat{{ϵ}_{\mathrm{pl}}}}$ is the equivalent plastic strain rate, ${\sigma }_0$ refers to the static yield stress, $\sigma$ defines the yield stress at the investigated strain rate, and $\gamma$ and $m$ are coefficients with no physical meaning. A deviation of 150 MPa was observed in the residual stresses of the FE simulation with the strain rate dependency compared to the FE simulation with a static model. The results imply that the strain rate dependency should be included when high accuracy is required; otherwise, it can be neglected since the deviation had a negligible influence on the overall residual stresses of the profile. Draxler et al. [47] predicted the crack susceptibility in the welding of Varestraint tests of alloy 718. The simulation was performed using data from different strain rates and temperatures. A rate-independent plasticity model was used between$20<T<{700}^{\circ }\mathrm{C}$, the Chen’s rate-dependent flow stress model [67] at the temperature range${700}^{\circ }\mathrm{C}<T<{1050}^{\circ }\mathrm{C}$, a simple rate-dependent material model based on the work of Antonsson and Fredriksson [68] was used between ${1050}^{\circ }\mathrm{C}<T\le T_c,$ and a flow stress model for soft isotropic solids was used when $T>T_c$.

2.4.5. Poisson’s ratio

The Poisson’s ratio relates the transverse strain to the axial strain. The Poisson’s ratio increases with temperature [37,54]. This is explained by a change in the flow stress of most materials from elastoplastic behaviour at low temperatures to visco-plastic behaviour at high temperatures. Babu et al. [42] approximated the variation of the Poisson’s ratio ($\mu$) and the temperature as: (21) $\mu =0.34+6.34e^{-5}T$(21) where $T$ refers to the temperature in Celsius. A constant Poisson’s ratio is used in some studies because $\mu$ has little influence on plastic deformation compared to the Young’s modulus [47,49].

The overall effects of mechanical properties are summarised in [44]. Apart for the yield strength, the mechanical properties at room temperature reasonably predict the transient temperature fields, residual stress, and distortion. The variation of mechanical and thermal features with temperature is well understood and replicable in the low-temperature range. However, measurements become very difficult close to the melting temperatures of materials due to contamination, reliability, and difficulties in fixing specimens.

Some raw data available for the temperature-dependent properties is provided in Table 1 to serve as a reference for finite element modelling. The raw data is mainly provided by government agencies and manufacturers. Temperature-dependent properties of several nickel-base, cobalt–nickel and iron-nickel super alloys are provided on the website of Special Metals Corporation and the nickel institute [69].

Table 1. List of publications with raw data on temperature dependent physical properties.

References	Material	EM	YS	FS	SR	PR	TE	TC	SH	CL	MP
[70]	Alloy 182	x		x			x	x	x	x	x
[70]	Low alloy steel	x		x			x	x	x	x	x
[70]	Stainless steel	x		x			x	x	x	x	x
[71]	Alloy 82/182	x		x			x	x	x
[71]	SA-105 CS	x		x			x	x	x
[71]	Stainless steel	x		x			x	x	x
[72]	Ti–6Al–4V			x	x
[73]	Ti–6Al–4V	x	x				x
[74]	316 Stainless Steel	x	x			x	x	x	x
[75]	316 Stainless Steel	x	x			x	x	x	x
[76]	Ti–6Al–4V	x	x				x	x	x
[77]	316 Stainless Steel	x	x			x	x	x	x
[78]	AA5052							x	x		x
[78]	AA6061							x	x
[79]	AA5083-O		x					x	x		x
[80]	AA 6061-T6	x	x				x	x	x
Key
Elastic modulus = EM	Thermal expansion = TE
Yield stress = YS	Thermal conductivity = TC
Flow stress = FS	Specific heat = SH
Strain rate sensitivity = SR	Creep law = CL
Poisson’s ratio = PR	Miscellaneous properties = MP (density, latent heat, solidus, liquidus, anneal temperate)

2.5. Phase transformation

Layer-by-layer material deposition in AM causes repeated melting and re-melting of the material, as indicated in Figure 4. The numerous solid-state and liquid–solid phase transformations modify the microstructure and properties of AM components. The thermal history of the part is different at each point and depends on the AM process, deposition strategy, and process parameters. This leads to the creation of multi-phase components with a heterogeneous distribution of the phases, microstructure, and residual stresses [81]. Phase transformations are predicted by density-based models that couple macroscale heat transfer to create the process-microstructure simulation rather than predict the actual grain morphology evolution [82]. The volumetric phase fractions are derived from phase transformation kinematic models that use volumetric field data (time, temperature, and strain rate).

Figure 4. Diagram of physical state change between raw/liquid/solid states (RLS) [83].

The choice of phase transformation models depends on the AM material and process conditions. In the following sections, we present a summary of materials commonly used in AM and semi-empirical phase transformation models.

2.5.1. Commonly used materials in AM

FE simulations of metal AM has mainly concentrated on steels, titanium, nickel, and aluminium alloys. AM also offers the possibility of developing new alloy systems by mixing material powders or rods and altering process conditions to achieve desired properties. It is widely accepted that the next critical step in the advancement of metal AM will be the development of process-specific alloys [84]. Figure 5 shows the different metal alloys investigated in AM processes in the decade between 2013–2023 (source: clarivate web of science).

• Steels

Figure 5. A pie chart of publications of the recent 10 years (2013–2023) for different AM alloy systems.

Steels are among the most widely used materials in AM due to their low cost, good corrosion resistance, high ductility, and reasonable strength. As the material solidifies at the liquidus temperature, the martensite start-temperature is determined by the carbon content of the steel, while the ferrite and bainite transformation kinematics are mainly influenced by the alloying elements (manganese, nickel, and chromium). The alloying elements control the driving force for austenite decomposition. During rapid weld thermal cycles, alloying elements and precipitates fail to completely solutionize and cause inhomogeneity in the austenite. Hence, continuous cooling transformation (CCT) diagrams derived from short-time holds in the austenitic temperature range are more suitable for predicting the HAZ_CCT behaviour in AM than the conventional approach of holding the specimen for over an hour above the austenite temperature [85]. Due to high solidification rates, columnar cellular or dendritic microstructures have been observed in AM steel components, compared to equiaxed grain structures in the wrought or annealed state [86].

• Titanium

Most research on FE simulations of AM of titanium has been performed on the Ti-6Al-4 V alloy because of its beneficial properties, such as lightweight, high strength, good corrosion resistance, and excellent mechanical performance at high temperatures. At room temperature, Ti-6Al-4 V comprises a hexagonal closed-packed α phase ($X_{\alpha }$) and a 5-10% body-centred cubic β phase ($X_{\text{β}}$). At temperatures approaching the transus temperature (882°C), the α phase converts to the body-centred cubic β phase. The phase transformation that occurs as the material temperature decreases below the transus temperature depends on the cooling rate. For example, the $x_{\text{β}}$ phase of Ti-6Al-4 V changes to $x_{\alpha }$ phase for diffusionless transformation and $x_{\alpha +\text{β}}$ phase for diffusion transformation, as shown in Figure 6. Pegues et al. [87] present some commonly observed microstructures in the AM of Ti-6Al-4 V and their influence on mechanical properties.

• Aluminum

Figure 6. (a) Phase transformations of titanium alloys. (b) Overview of phase change model for Ti-6Al-4 V.[88]

The commercial usage of AM components made from aluminium is limited compared to parts made from steel and titanium alloys. Nevertheless, the cost of AM aluminium components is lower compared to the other alloys. Aluminium is hardened mainly by the formation of precipitates (solid-state transformation) that are formed in different alloy systems. The main alloy systems in most AM studies consist of the Al–Si (Mg) system, the Al-Cu system, the Al-Zn-Mg system, and the Al-Mg-Sc-Zr system [89]. The dissolution of the alloying elements occurs by holding the material at elevated temperatures. The formation of phases after cooling depends on the ageing temperature. At first, small clusters of solute atoms are formed, which evolve into larger clusters. The clusters induce elastic straining in the surrounding matrix, which strengthens the alloy. Kannan et al. [90] reviewed the AM of AA2024 and AA6061, two heat-treatable alloys with high strength and a good strength-to-weight ratio.

• Inconel

Inconel is a nickel-based superalloy extensively used in aerospace, energy, and other industries because of its excellent properties (tensile strength, creep strength, weldability, and corrosion resistance at elevated temperatures). Repeated heating in the AM produces components with substantially different microstructures compared to parts made by casting and forging. Inconel is hardened by precipitation. Lu et al. [91] studied the evolution of precipitates in multi-layer and multi-pass WA-DED using thermo-mechanical-metallurgical phase transformation kinematics. The microstructure changes include nucleation, growth, and collision of grains. Gao et al. [92] investigated the effects of different solution heat treatments on the microstructure and high-temperature mechanical properties of Inconel 718. The distribution and morphology of the δ phase in the γ matrix determined the mechanical properties.

2.5.2. Semi-empirical phase transformation equations

In this section, mainly semi-empirical models that predict the decomposition of the austenite into different phases upon cooling are presented, since the heating process to austenitic temperatures typically has negligible influence on the final phase distribution. The semi-empirical models are derived from CCT diagrams. The models used in the heating process are presented when provided in the literature. Above the austenite temperature, the parent (austenite) phase is assumed to constitute 100% of the material. The phase fracture of the products depends on the cooling rate. Commonly used phase transformation models in macroscale simulations are presented in Table 2.

Table 2. Phase transformation kinematic models.

Models	References	Models
Koinstinen-Marburg (KM)	[93]	$f_{{\alpha }^{\mathrm{\prime }}}\left(t,T\right)=\left[1-\exp \left(-k_{\text{β}\to {\alpha }^{\mathrm{\prime }}}\left(M_s-T\right)\right)\right]f_{\text{β}}\left(t,T\right)$ where $f_{{\alpha }^{\mathrm{\prime }}}$ is the fraction of ${\alpha }^{\mathrm{\prime }}$, $f_{\text{β}}$ is the fraction of $\text{β}$, $k_{\text{β}\to {\alpha }^{\mathrm{\prime }}}$ is the material dependent KM parameter, $M_s$ is the martensite start temperature and$T$ is the temperature.
JMA	[94]	$f\left(t\right)=1-\exp \left(-k{t}^n\right)$ where $f\left(t\right)$, $k$ and $n$ are the fraction of phase after time t, the reaction rate constant and the Avrami exponent, respectively.
JMAK	[93]	$f_{\alpha }\left(t,T\right)=\left[1-\exp \left(-k_{\text{β}\to \alpha }\left(T\right)t^{n_{\text{β}\to \alpha }}\right)\right]f_{\alpha }^{\mathrm{eq}}\left(T\right)$ where is $f_{\alpha }$ the fraction of the $\alpha$ phase at time t, $f_{\alpha }^{\mathrm{eq}}$ is the equilibrium fraction of the $\alpha$ phase, and $k_{\text{β}\to \alpha }$ and $n_{\text{β}\to \alpha }$ are the temperature dependent and independent JMAK parameters that define the kinetics of the $\text{β}$ to $\alpha$ transformation.
Waeckel	[95]	$\dot{\xi }\left(t\right)=f\left(T,T,\dot{\xi },M_s;d\right)=f\left(T,T,\dot{\hat{\xi }};d\right){\displaystyle \frac{〈T-M_S〉}{T-M_S}}$ where ξ are constituents that transform with a diffusional mechanism, and $X$ is the positive part of $X$. $T-M_S$ states the martensitic transformation start temperature (Ms).
Leblond	[96]	${\displaystyle \frac{d{Z}_{\gamma }}{\mathrm{dt}}={\displaystyle \frac{Z_{\mathrm{eq}}\left(T\right)-Z_{\gamma }}{\tau \left(T\right)}}}$ $Z_{\mathrm{eq}}\left(T\right)$ is the equilibrium proportion, is reached at constant temperature, τ indicates the characteristic time necessary for the phase to reach an equilibrium state and T is the temperature.

Most studies use the KM model for diffusionless phase transformation and the JMAK (JMA) for diffusion phase transformation. Diffusion phase transformation refers to phase transformations that occur by atomic diffusion, i.e. an interchange of the positions of different atom species in the crystal lattice. Diffusional transformations usually happens when the cooling rate from the austenitic temperature is slow (lower than ∼400$\mathrm{K}/\mathrm{s}$); otherwise, non-diffusional transformations occur. A few studies use the Leblond model to simulate diffusional phase transformation. In most studies, austenite decomposition is modelled by a combination of diffusional and non-diffusional transformations to replicate experimental findings. The linear mixture rule [42], the nonlinear mixture rule [57], the rule of mixtures [97], the lever rule, or Scheil’s additive rule are commonly used to define the mixture rule of phase transformations. Table 3 presents phase transformation models and materials in the AM literature.

Table 3. A summary of phase transformation models in AM simulation.

References	Material	KM	JMAK	JMA	Leblond
[65]	steel S355, VATM1 and VATM2				X
[97]	DP 980		X
[57]	EH36 steel	X			X
[98]	ER90S-B91 steel	X
[6]	AF1410 steel	X
[99]	Ti-6Al-4V	X		X
[42]	Ti-6Al-4V			X
[100]	Ti-6Al-4V			X
[101]	Ti-6Al-4V			X
[5]	Ti-6Al-4V			X
[39]	AF1410 steel	X			X
[102]	Ti-6Al-4V	X	X
[103]	Ti-6Al-4V			X
[83]	Ti-6Al-4V+steel 5140	X		X
[104]	Ti-6Al-4V	X		X
[88]	Ti-6Al-4V	X	X
[105]	AISI 420	X
[94]	Ti-6Al-4V	X		X
[93]	Ti-6Al-4V	X	X
[63]	LTT alloy	X			X
[106]	TA6V	X
[43]	S355J2	X		X
[107]	Ti-6Al-4V	X		X
[82]	Ti-6Al-4V	X	X
[108]	6MND5 or SA533	X			X
[109]	Ti-6Al-4V	X	X
[110]	Ti-6Al-4V	X	X
[111]	Ti-6Al-4V	X	X
[112]	Ti-6Al-4V	X	X
[91]	Inconel		X
[113]	Ti-6Al-4V	X	X
[114]	Ti-6Al-4V	X	X
[115]	Ti-6Al-4V		X
[116]	Ti-6Al-4V	X	X

2.5.2.1. Diffusional phase transformations

• Johnson-Mehl-Avrami-Kolmogrov (JMAK) model

Diffusion transformation in AM is widely modelled using the austenite decomposition procedure proposed by Johnson and Mehl [117], Avrami [118], and Kolmogorov (1937). The model and its derivatives is usually abbreviated as JMA, JMAK, or KJMA and describes the time evolution of transformed phase fractions of materials (Ti-6Al-4 V, steels, and Inconel). Fan et al. [5] modelled the phase transformation in the laser AM of Ti-6Al-4 V. The heating phase was performed using the JMA model. Superheating and undercooling were controlled by the k and n parameters of the JMA model. Depending on the rate of cooling, the β-phase is transformed into the α and/or αʹ-phase(s). Vastola et al. [103] compared the evolution of the microstructure of Ti-6Al-4 V manufactured by electron beam-DED and LPBF. At fast cooling rates, the dissolution of the β-phase into the αʹ-phase was described as in [119]. (22) ${\chi }_r\left(t\right)=1-e^{-\text{β}{t}^k}$(22) The distribution of the αʹ-phase was more abrupt in LPBF compared to electron beam-DED, which is explained by the lower build temperature and higher thermal gradients in LPBF, which limit the areas for martensite formation.

Murgau et al. [99] coupled the phase fraction and microstructure evolution (α lath size) in FE simulation to solve the temperature field, phase change, and microstructural evolution problems in GTAW of Ti-6Al-4 V. It was envisaged that the microstructure from the FE simulations could predict the future material properties of AM components. Irwin [101] implemented Kelly's microstructure model and Charles’ lath width model to study the influence of process parameters (dwell time) on the microstructure in AM of Ti-6Al-4 V. The growth and nucleation of the α-phase were described by the JMA model. The α lath width is defined by the work of Charles, following the Arrhenius equation. (23) ${\omega }_{\mathrm{eq}}=k_w\exp \left(-T_{\mathrm{act}}/T_i\right)$(23) where k_w is the prefactor, −T_act is the activation temperature and T_i is the actual temperature.

Babu et al. [42] presented a microstructural thermo-mechanical model to predict residual stresses and distortion in the AM of Ti-6Al-4 V. The phase evolution was defined by the JMA model. An advanced isotropic plasticity model (with long- and short-range stress components, a hardening process, and the evolution of immobile dislocation density) was used to define the flow stress. Beyond the fourth welding pass, there was no change in the material’s microstructure. This is corroborated by results from the literature.

• Leblond

The Leblond diffusion-controlled phase transformation model [96] has been implemented in the FE simulation of AM. Chen et al. [63] investigated the microstructure evolution and resulting residual stresses in the AM of low-temperature transformation (LTT) pre-alloyed powders. The Leblond model was used to describe the ferrite-austenite transformation during heating and cooling. The KM phase transformation model expressed the inverse displacive-controlled transformation from austenite to martensite. Gkatzogiannis [65] examined the effects of strain rate in the FE investigation of welding and included the effect of Leblond phase transformation. Ali and Han [57] investigated the influence of five scanning strategies on the residual stresses and warpage in WA-DED with a single-layer deposition. Two material models (with Leblond phase change and without phase change) were simulated to determine the influence of phase change on the predictions.

2.5.2.2. Diffusionless transformations

• Koistinen–Marburger (KM)

The KM model mainly describes rapid cooling (diffusionless transformation) in the literature and is valid for many ferrous alloys [120]. Jimenez et al. [98] compared the residual stresses (XRD) in WA-DED of B91 steel to data from thermo-mechanical simulations. The residual stress prediction improved in the phase transformation model by over 230%, demonstrating that phase transformation should be included in the simulation of WA-DED in the process development stage. Zhang et al. [6] modelled the phase transformation of AF1410 steel based on continuous cooling transformation (CCT) curves obtained by thermomechanical tests. The residual stresses from the FE simulations were validated by the advanced contour method and X-ray diffraction (XRD). In the same light, Zhang et al. [39] showed that the microstructure of AF1410 steel in multi-pass welding is mainly composed of martensite. Phase transformation caused deformation and residual stresses during cooling. Depradeux et al. [106] performed thermo-mechanical-metallurgical FE analyses of the rapid additive forging of TA6 V.

2.5.2.3. Combining phase transformation equations

Most AM studies with phase transformations include diffusional and diffusionless phase transformations. Zhang et al. [83] presented a combined experimental and numerical methodology to estimate the process-structure–property-performance relationships of AM components. The KM model was used to simulate rapid cooling, while the JMA model was used to predict slow cooling and (re)heating operations. The methodology was verified on Ti-6Al-4 V and steel 5140. Salsi et al. [102] combined the JMAK and KM phase change laws to predict the phases in the LPBF of Ti-6Al-4 V. The JMAK model was used to model the slow cooling of the β phase to alpha Widmanstätten $\left({\alpha }_w\right)$, while the K-M model was used to predict fast cooling (formation of alpha martensite, ${\alpha }_m$). The additivity rule is used to describe phase evolution. Reheating was performed to account for incomplete and multiple-phase transformations in AM. Baykasoğlu et al. [82] simulated the microstructure of Ti-6Al-4 V produced by powder L-DED. The metallurgical phase transformation from the beta-transus temperature included the effects of Widmanstätten colony alpha, basketweave alpha, grain boundary at slow cooling rates, and martensite alpha at high cooling rates. The microstructure evolution model was validated by comparing the volumetric phase fractions of the model to the literature. Murgau et al. [109] developed analytical models that predict solid-state transformations of Ti-6Al-4 V subjected to arbitrary temperature changes. The $\text{β}$ phase dissolves into three $\alpha$ phases (widmannstätten α_wid, grain boundary α_gb, and martensite) after cooling. The diffusionless transformation (${\alpha }^{\mathrm{\prime }}$) was modelled in the KM model. The JMAK model was modified to include incomplete transformation, multiple-phase transformation, and multiple heating and cooling cycles.

2.5.2.4. Phase transformation in multiple heating and cooling

The phase composition of a material may alter if the material is heated to sufficiently high temperatures by successive AM layers. The previous heating cycles have an effect on the phase transformation during reheating. In light of this, Xie et al. [104] implemented a general physics-based framework to perform thermal and mechanical analysis of AM. Rapid cooling was modelled with the KM model. The JMA model was used for slow cooling and re-heating. Similarly, Crespo [94] modified the diffusional transformations in the AM of Ti-6Al-4 V with repeated heating and cooling cycles. If the deposited material is composed mainly of$\alpha +\text{β}$, reheating leads to diffusion-controlled $\alpha \to \text{β}$ transformation. Alternatively, if the microstructure is martensitic (${\alpha }^{\mathrm{\prime }}$), at temperatures above 400°C ${\alpha }^{\mathrm{\prime }}$ converts to $\alpha +\text{β}$. The second diffusion transformation of the fraction of the α phase $f_{\alpha }\left(t_0\right)$ is described by the JMA equation: (24) $f_{\alpha }\left(t_1\right)=\left\{1-\exp \left[-k_0{\left(t_1-t_0\right)}^{n_0}\right]\right\}.f_{\alpha }^{\mathrm{eq}}\left(T_0\right)$(24) where $k_0$ is the reaction rate constant and $n_0$ is the Avrami exponent at the temperature $T_0$ for the time interval, $\left[t_0,t_1\right[$. In the next reheat interval, the transformation temperature is assumed as$T_1$. If we consider the phase fraction of α in the previous step$f_{\alpha }\left(t_1\right)$, the α phase in the next time interval $\left[t_1,t_2\right[$ is obtained by (25) $f_{\alpha }\left(t_2\right)=\left\{1-\exp \left[-k_1{\left(t_1^f+t_2-t_1\right)}^{n_1}\right]\right\}.f_{\alpha }^{\mathrm{eq}}\left(T_1\right)$(25)

The generalised form of the equation of the α phase fraction at an arbitrary time step $\left[t_s,t_{s+1}\right[$ at temperature $T_s$ is written as (26) $f_{\alpha }\left(t_{s+1}\right)=\left\{1-\exp \left[-k_s{\left(t_s^f+t_{s+1}-t_s\right)}^{n_s}\right]\right\}\cdot f_{\alpha }^{\mathrm{eq}}\left(T_s\right)$(26) If martensitic decomposition is incomplete, the microstructure consists of ${\alpha }^{\mathrm{\prime }}+\alpha +\text{β}$. During cooling, the $\text{β}$ phase converts to $\alpha$ phase if the cooling rate is lower than ∼400$\mathrm{K}/\mathrm{s}$. At cooling rates higher than 410°C/s, the $\text{β}$ phase is retained or converted to ${\alpha }^{\mathrm{\prime }}$ phase depending on the volume fraction of$\text{β}$. The martensitic transformation was modelled using the KM transformation. Confer [94] for the full definition of terms in Equations (24)–(26). Klusemann and Bambach [110] used the modified JMAK in [109], the KM model, and the parabolic growth rate to predict the martensite phase in laser metal deposition of Ti-6AL-4 V. Figure 7 shows the phase evolution with temperature in the AM of Ti-6Al-4 V.

Figure 7. (a) Temperature variations vs time. (b) Corresponding simulated α_gb, α_wid, and β phase fractions vs time. (c) Corresponding simulated α_gb, α_wid, and β phase fractions vs. temperature. [99].

The studies presented so far are mainly based on CCT and TTT. However, AM processes cause isothermal, non-isothermal, and kinematic transformations due to multiple heating and cooling. Tchuindjang et al. [100] proposed a new concept to simplify phase change modelling that considers the time-temperature transformation (TTT), continuous cooling transformation (CCT), and continuous heating transformation (CHT). They introduced a time-phase transformation block (TTB) to account for TTT, CCT, and CHT since the temperature and process parameters in AM are constantly changing and affect the microstructure. The different phase transformation mechanisms were activated depending on the temperature value and the heating or cooling rate. The microstructure can be predicted based on a good understanding of the phase changes predicted by phenomenological models like JMAK, KM, or neural network-based models.

2.5.3. Hardness prediction from phase transformation

The hardness of the material is thought to be the same as the weighted average hardness of its different phases. This is based on mixture laws that explain how different fractions of a single phase affect hardness. Costa et al. [105] identified martensite and austenite as the main phases in AISI 420 steel. They used the following equation to derive the hardness of the AM component: (27) $H=f_{\mathit{\gamma }}\cdot H_{\gamma }+\left(1-f_{\mathit{\gamma }}\right)\cdot H_M$(27) where $H_{\gamma }$ and $H_M$ refer to the hardness of austenite and martensite respectively, and $f_{\gamma }$ is the volume fraction of austenite. The softening of the material caused by non-isothermal tempering of martensite was approximated by Réti et al. [121] using the following equation: (28) $H_M=H_M\left(T=0\right)-A\cdot {\left[{\int }_0^t{e}^{-Q/R\cdot T\left(t\right)}\mathrm{dt}\right]}^m$(28) where $H_M\left(T=0\right)$ defines the hardness of martensite, $R$ is the universal gas constant, $Q$ refers the activation energy, and $A$ and $m$ are fitting constants. Using the same formulation, Costa et al. [122] coupled the hardness and volume fraction of martensite to the temperature evolution in laser powder deposition of a thin wall of AFNOR NF 50CV2 steel. Similarly, Crespo and Vilar [107] investigated the phase transformation of Ti-6Al-4 V using the JMA equation. Anisothermal conditions and the additivity rule were used to obtain the phase concentrations. A rule of mixtures is used to calculate hardness.

2.6. Coupling of physical metallurgical effects in macroscale AM simulations

High temperatures in AM processes lead to changes in the microstructure of materials, affecting the mechanical properties of AM components [88]. Achieving a homogeneous microstructure in entire AM parts facilitates process design; however, certain components may require varied microstructures in various portions. Thermo-mechanical FE simulations can be conducted to match the microstructure and temperature evolution of AM components without explicitly coupling them. For instance, Tran et al. [52] correlated the microstructure and temperature evolution at 3 points in laser cladding involving re-melting stages. Rapid heating rates caused a shift of the transformation points, including${\text{β}}_{\mathrm{transus}}$, and solidus, while re-melting caused pseudo-isothermal annealing. While such models improve the predictability of FE simulations, coupling microstructure to FE simulations is essential for process chain simulations and performance evaluation.

2.6.1. Physically based flow stress

Physically-based flow models can model the properties of individual phases and microstructure evolution, thereby providing a better description of phenomena in AM simulations. In the last decade, microstructure-coupled FE simulations have been employed to predict AM component microstructure, minimise statistical variability, and simplify design [123]. Physically based flow stress models of materials are determined using a material volume that represents the slip systems, mobile dislocations, grain size, etc. The models assume that the motions of dislocations are the dominant source of plastic strain. The physically-based flow stress models can reveal intrinsic physical mechanisms and represent multifaceted microstructure evolution [124,125]. The models can be found in [42,88,126,127]. The flow stress consists of the long-and short-range stress components. (29) ${\sigma }_y={\sigma }_G+{\sigma }_{\mathrm{HP}}+{\sigma }^{\ast }$(29) where, ${\sigma }_G$ refers to the interactions with the dislocation substructure (forest dislocations), ${\sigma }_{\mathrm{HP}}$ is the Hall-Petch effect, and ${\sigma }^{\ast }$ is the friction stress needed to move dislocations through the lattice and to go pass short-range obstacles (solutes or precipitates).

2.6.1.1. Long-range stress component

The long-range stress component${\sigma }_G$ is the physical basis for the strain hardening of metals and is defined by [128] as (30) ${\sigma }_G=\mathrm{m\alpha Gb}\sqrt{{\rho }_i}$(30) where $m$ is the Taylor orientation factor, $\alpha$ refers to the proportionality factor, $G$ represents the shear modulus, $b$ is the Burger’s vector and ${\rho }_i$ represents the immobile dislocation density. The evolution of the dislocation density, ${\rho }_i$ in Equation (30) is computed to describe the hardening and restoration processes. (31) ${\dot{\rho }}_i={\dot{\rho }}_i^{(+)}+{\dot{\rho }}_i^{(-)}$(31)

During hardening, mobile dislocations are assumed to move a distance $\mathrm{\Lambda }$ (mean free path) on average before they are annihilated or immobilised. The mobile dislocation densities and their velocities are assumed to be proportional to the effective plastic strain rate $\left({\dot{\overline{ϵ}}}^p\right)$ [129]. (32) ${\dot{\overline{ϵ}}}^p={\displaystyle \frac{{\rho }_m\overline{v}b}{M}}$(32) where ${\rho }_m$ refers to the density of mobile dislocations, $\overline{v}$ describes the average velocity of dislocation densities, and $M$ is the Taylor lattice factor. The immobile dislocation density ${\dot{\rho }}_i^{(+)}$is also assumed to follow the same relationship as the mobile dislocation densities [130]. Thus, ${\dot{\rho }}_i^{(+)}$ is defined as (33) ${\dot{\rho }}_i^{(+)}={\displaystyle \frac{m}{b}{\displaystyle \frac{1}{\mathrm{\Lambda }}{\dot{\overline{ϵ}}}^p}}$(33)

The mean free path is determined by the grain size (g) and dislocation sub cell or sub grain diameter (s) as follows: (34) ${\displaystyle \frac{1}{\mathrm{\Lambda }}=\left({\displaystyle \frac{1}{g}+{\displaystyle \frac{1}{s}+\mathrm{others}}}\right)}$(34) others here represent obstacles like precipitates, martensite lathes, interstitial elements, etc. The formation and evolution of the sub cells are based on a relationship proposed by [131]. (35) $s=K_c{\displaystyle \frac{1}{\sqrt{{\rho }_i}}}$(35)

• The restoration processes

Babu et al. [42] considered the mechanisms in restoration processes as a combination of dislocation glide, dislocation climb, and globularization. (36) ${\dot{\rho }}_i^{(-)}={\dot{\rho }}_i^{\left(\mathrm{glide}\right)}+{\dot{\rho }}_i^{\left(\mathrm{climb}\right)}+{\dot{\rho }}_i^{\left(\mathrm{globularization}\right)}$(36) Dislocation glide is modelled based on the Bergström formulation. (37) ${\dot{\rho }}_i^{\left(\mathrm{glide}\right)}=\mathrm{\Omega }{\rho }_i{\dot{\overline{ϵ}}}^p$(37) where $\mathrm{\Omega }$ refers to a temperature dependent function. Militzer et al. [132] proposed a model for dislocation climb based on the works of [133] and [134]. (38) ${\dot{\rho }}_i^{\left(\mathrm{climb}\right)}=2C_{\gamma }{D}_{\mathrm{app}}{\displaystyle \frac{G{b}^3}{\mathrm{kT}}\left({\rho }_i^2-{\rho }_{\mathrm{eq}}^2\right)}$(38) where $C_{\gamma }$ is a material coefficient, $D_{\mathrm{app}}$ is the apparent diffusivity that includes the weighted factions of the phases, and ${\rho }_{\mathrm{eq}}^2$defines the equilibrium value of the dislocation density. Sargent et al. [135] observed static and dynamic coarsening of grains in Ti-6Al-4 V at elevated temperatures (850°C). Coarsening results in the grains becoming more globular, hence the term globularization. Babu and Lindgren [126] modelled the evolution of the dislocation density during globularization. The effect of grain growth on the reduction of flow stress is considered above the critical value of dislocation density ${\rho }_{\mathrm{cr}}$ $\begin{array}{rl}& \mathrm{if}{\rho }_i\ge {\rho }_{\mathrm{cr}}\\ & {\dot{\rho }}_i^{\left(\mathrm{globularization}\right)}=\psi {\dot{X}}_g\left({\rho }_i-{\rho }_{\mathrm{eq}}\right);\mathrm{until}{\rho }_i\le {\rho }_{\mathrm{eq}}\\ & \mathrm{else}\\ & {\dot{\rho }}_i^{\left(\mathrm{globularization}\right)}=0\end{array}$where$\psi$ is a calibration constant and ${\dot{X}}_g$ is the rate of globularization. Based on the observation that during deformation and shortly afterwards, globularization is due to dislocation sub-structuring resulting in boundary split and edge spheroidization in Ti-6Al-4 V by Semiatin et al. [136], Thomas and Semiatin [137] modelled a two-stage process of static and dynamic recrystallization. They proposed Equation (39) to model globularization since there are similarities between recrystallization and globularization. (39) $X_g=X_d+\left(1-X_d\right)X_s$(39) where $X_g$, $X_d$ and $X_s$ are the volume fractions of the total, dynamic and static components of the globularization respectively. If the grain growth and static recrystallization are assumed to be driven by the same driving force, the static recrystallization can be determined by (40) ${\dot{X}}_s=M{\displaystyle \frac{\dot{g}}{g}}$(40) where $M$ is a material parameter. The rate of globularization is modelled based on [137] (41) ${\dot{X}}_d={\displaystyle \frac{B{k}_p{\dot{\overline{ϵ}}}_p}{{\overline{ϵ}}_p^{1-k}{e}^{B{\overline{ϵ}}_p^{k_p}}}}$(41) where $B$ and $k_p$ are material parameters.

2.6.1.2. Hall-Petch

The Hall-Petch effect leads to the activation of slip systems in neighbouring grains and is derived from (42) ${\sigma }_{\mathrm{HP}}={\displaystyle \frac{k_{\mathrm{HP}}}{\sqrt{g}}}$(42) where $g$ refers the grain size and $k_{\mathrm{HP}}$ is a coefficient that depends on the critical resolved shear stress of the material and shear modulus. The Hall-Petch effect has an indirect temperature dependence due to changes in elastic properties and is defined as (43) ${\sigma }_{\mathrm{HP}}={\displaystyle \frac{G\left(T\right)}{G_{\mathrm{RT}}}{\displaystyle \frac{k_{\mathrm{HP}}}{\sqrt{g}}}}$(43) where $G\left(T\right)$ defines the shear modulus at the current temperature and $G_{\mathrm{RT}}$ is the shear modulus at room temperature. However, the temperature dependency of ${\sigma }_{\mathrm{HP}}$ is not valid for very small grain sizes.

2.6.1.3. Short-range stress component

The short-range stress component ${\sigma }^{\ast }$is the material’s resistance to plastic deformation due to short-range interactions formulated by Mecking and Kocks [130] (44) ${\sigma }^{\ast }={\tau }_{0}G{\left[1-{\left[{\displaystyle \frac{\mathrm{kT}}{\mathrm{\Delta }{f}_{0}G{b}^3}\ln \left({\displaystyle \frac{{\dot{ϵ}}^{\mathrm{ref}}}{{\dot{\overline{ϵ}}}^p}}\right)}\right]}^{1/q}\right]}^{1/p}$(44) where $\mathrm{\Delta }{f}_{0}G{b}^3$ describes the activation energy needed to overcome lattice resistance, ${\dot{ϵ}}^{\mathrm{ref}}$ is the reference strain rate, ${\tau }_{0}G$ refers to the shear strength in the absence of thermal energy, ${\dot{\overline{ϵ}}}^p$ represents the dislocation velocity, $T$ is the temperature in Kelvin, $k$ is the Boltzmann constant and $p,q,f_0$ are calibration parameters. Furthermore, phenomena like precipitation hardening, vacancy annihilation and a fine α lath width increase ${\sigma }^{\ast }$.

2.6.1.4. Precipitation hardening

Lindgren et al. [88] included the contribution of precipitate hardening in the AM of Inconel. The stress required to move a dislocation past an obstacle is defined based on [138] by: (45) ${\sigma }_p={\displaystyle \frac{\overline{F}M}{b{l}_p}}$(45) where $\overline{F}$ refers the mean obstacle strength and $l_p$ is the average particle spacing in the slip plane.

2.6.1.5. Vacancy evolution and annihilation

Materials generate excess vacancies when subjected to deformation or temperature changes. Assuming that only long-range stress causes vacancy formation and considering the contribution of temperature change, Militzer et al. [132] proposed the following model for excess vacancy concentration with generation and annihilation components: (46) ${\dot{c}}_v^{\mathrm{ex}}=\left[\chi {\displaystyle \frac{\mathrm{m\alpha G}{b}^2\sqrt{{\rho }_i}}{Q_{\mathrm{vf}}}+\zeta {\displaystyle \frac{c_j}{4b^2}}}\right]{\displaystyle \frac{{\mathrm{\Omega }}_0}{b}\dot{\overline{ϵ}}{D}_{\mathrm{vm}}\left[{\displaystyle \frac{1}{s^2}+{\displaystyle \frac{1}{g^2}}}\right]\left(c_v-c_v^{\mathrm{eq}}\right)+c_v^{\mathrm{eq}}\left({\displaystyle \frac{Q_{\mathrm{vf}}}{k{T}^2}}\right)\dot{T}}$(46) where ${\mathrm{\Omega }}_0$ is the atomic volume, $\chi$ is the fraction of mechanical energy spent on vacancy generation, and $\zeta$ is the neutralisation effect by vancancy emitting and absorbing jogs.

2.6.1.6. α lath width prediction

The α lath width is a key microstructural parameter that determines the mechanical properties of Ti-6Al-4 V. Irwin [101] modelled the α lath width by (47) $w_i={\displaystyle \frac{1}{F_{i,\alpha }}\left[w_{i-1}{F}_{i-1,\alpha }+w_{\mathrm{eq}}\left(F_{i,\alpha }-F_{i-1,\alpha }\right)\right]}$(47) where $F_{i,\alpha }$ refers to the α phase fraction, and $w_{\mathrm{eq}}$ describes the equilibrium α lath width, determined by, (48) $w_{\mathrm{eq}}=k_w\exp \left(-T_{\mathrm{act}}/T_i\right)$(48) where $k_w$ describes the Arrhenius pre-exponential factor and $T_{\mathrm{act}}$is the activation energy. Furthermore, Yang et al. [113] modelled the dislocation density at the interior and boundary of the ${\alpha }^{\mathrm{\prime }}$ lath using 2 separate equations to consider the pile-up effect.${\alpha }^{\mathrm{\prime }}$ lath interior ($\mathrm{li}$) was defined by (49) ${\rho }_{\mathrm{li}}={\displaystyle \frac{12E_{{\alpha }^{\mathrm{\prime }}}{s}_{{\alpha }^{\mathrm{\prime }}}}{\left(1+2v^2\right){\mu }_{{\alpha }^{\mathrm{\prime }}}{b}_{{\alpha }^{\mathrm{\prime }}}}{\displaystyle \frac{{ϵ}_{\mathrm{trans}}^2}{w_{{\alpha }^{\mathrm{\prime }}}^2}}}$(49) where $E_{{\alpha }^{\mathrm{\prime }}}$ refers to the Youngs modulus of ${\alpha }^{\mathrm{\prime }}$,$s_{{\alpha }^{\mathrm{\prime }}}$ defines the thickness of $\mathrm{lb}$, $v$ is the Poisson’s ratio, ${\mu }_{{\alpha }^{\mathrm{\prime }}}$ represents the shear modulus of ${\alpha }^{\mathrm{\prime }}$, $b_{{\alpha }^{\mathrm{\prime }}}$is the burgers vector of ${\alpha }^{\mathrm{\prime }}$, ${ϵ}_{\mathrm{trans}}$ describes the transformation strain from β.${\alpha }^{\mathrm{\prime }}$ lath boundary ($\mathrm{lb}$) was defined by (50) ${\rho }_{\mathrm{lb}}={\displaystyle \frac{{\mathrm{\Lambda }}_l{P}_l}{2A_L}}$(50) where ${\mathrm{\Lambda }}_l$ is the number of $\mathrm{lb}$ dislocations per metre of boundary, $P_l$ represents the lath perimeter and $A_L$ is the $\mathrm{lb}$ section area. Other definitions used in the model are found in [113].

2.6.2. AM studies with physically based flow stress

Microstructure-coupled FE simulations have been performed principally on Ti-6AL-4 V due to the highly desired characteristics of AM components made from Ti-6AL-4 V. Lindgren et al. [88] implemented a multi-class approach to modelling the nucleation, growth, and coarsening of precipitates. The method predicts the distribution and evolution of particle size in spatial space based on the dislocation density. The athermal long-range contribution of the flow stress was modelled while considering the interactions of dislocations and dislocation substructures (strain hardening) and the interactions between mobile and immobile dislocation densities (the Hall-Petch effect).

The size dependency of AM components increases the necessity of determining the microstructure by FE simulation before production. Promoppatum et al. [112] investigated the effect of feature size on the cooling rate, phase composition, and grain size of laser-deposited Ti-6AL-4 V. The Hall-Petch model determined the grain size. Overall, there was a slight change in the yield strength and phase composition due to part size. Mirkoohi et al. [139] compared the residual stresses predicted by a physical model and the Johnson-Cook empirical flow stress model in the direct metal deposition of Ti-6Al-4 V. The physically-based model outperformed the Johnson-Cook model because it included material attributes like grain size and dislocation density in the flow stress. To reduce the simulation time, Lundbäck and Lindgren [140] applied the lumping of welds with a physically based material model. Overall, the computed temperatures and part displacements were close. However, microstructural investigations were not performed to validate the approach.

The FE determination of the mechanical properties of AM components made from materials other than Ti-6AL-4 V is in the initial stages and not widely reported in the literature. Eskandari Sabzi and Rivera-Díaz-del-Castillo [141] proposed a methodology to predict the yield strength of three different alloys (Ti-6Al-4 V, IN718, and 316L) produced by LPBF. Hardening and softening were revealed by four principal mechanisms: Kock-Mecking’s contribution ${\displaystyle \frac{\mathrm{d\rho }}{\mathrm{dϵ}}}$, Taylor strengthening ${\displaystyle \frac{d{\sigma }_Y}{\mathrm{d\rho }}}$, thermal straining ${\displaystyle \frac{\mathrm{dϵ}}{\mathrm{dT}}}$, and heat dissipation cooling ${\displaystyle \frac{\mathrm{dT}}{d{H}_n}}$. Kock-Mecking’s, Taylor strengthening, and thermal straining are the mechanisms that control superior mechanical properties. Kock-Mecking’s contribution had the most significant effect on yield strength.

It remains challenging to predict the mechanical properties of AM components because of the influence of multiple process parameters on the microstructure. A summary of the literature on physically based flow stress models used in AM process simulation is provided in Table 4.

Table 4. Summary of physically based models in AM simulation.

References	Materials	Globularization	Precicipitation	Vancancy	Dislocaation	Rule of mixture	Grain size	Morphology	Stress	Homogenisation
[62]	Alloy 625			X	X
[97]	DP 980					X	X
[88]	Ti-6Al-4 V, IN718		X		X		X		X	X
[99]	Ti-6Al-4V			X	X				X
[42]	Ti-6Al-4V	X	X	X	X	X	X		X
[140]	Ti-6Al-4V			X	X				X
[142]	Ti-6Al-4V			X	X		X		X
[126]	Ti-6Al-4V	X	X	X	X		X		X
[111]	Ti-6Al-4V	X			X	X	X		X
[113]	Ti-6Al-4V				X			X
[139]	Ti-6Al-4V				X		X		X
[115]	Ti-6Al-4V				X	X	X		X
[141]	Ti-6Al-4 V, IN718, 316L				X

3. Modelling the post-processing of AM components

Models used in the post-processing of AM are divided into the following categories: flow stress models for forming operations, performance analysis models, and models used in machining operations. FE simulations of the processes require material models that account for thermal, mechanical, and microstructural phenomena. Coret and Combescure [108] described the coupling mechanisms in material models for traditional manufacturing processes in Figure 8.

Figure 8. Coupling mechanisms [108].

3.1. Forming AM components

Post-forming of AM components is usually executed using empirical flow stress models [4,143]. Empirical flow stress models assume that AM components are homogenous with no metallurgical alterations (phases or microstructures). However, AM parts undergo significant thermal fluctuations in production, resulting in enormous microstructure variation within a single part and across production batches.

Hot forging is commonly applied to large AM components to reduce processing and solidification defects, attain anisotropic properties, and modify the microstructure to improve the mechanical properties. Empirical flow stress models based on experimental data have been developed to simulate hot forging [144–146]. Niu et al. [147] compared the predictive power of empirical (Johnson-Cook, Fields-Backofen, and strain-compensated Arrhenius) models in predicting flow stresses during the hot forming of powder L-DED Ti6-Al-4 V alloy. Motallebi et al. [148] reviewed the hot deformation behaviour and empirical flow stress models used in the deformation of AM components. In this section, only semi-empirical and physically-based flow stress models are reviewed to propose process chain simulation frameworks that include the effects of microstructure and phase change.

• Semi-empirical and physically based models

Few studies have used semi-empirical or physically-based flow stress models in the FE simulation of AM components. Bambach and Sizova [149] presented a semi-empirical model of the flow stresses in the hot working of Inconel 718 specimens produced by LPBF, powder L-DED, and a conventional wrought alloy that considers DRX. The flow stress due to strain hardening or softening is defined by (51) ${\sigma }_0={\sigma }_p{\left[{\displaystyle \frac{ϵ}{{ϵ}_p}\cdot \exp \left(1-{\displaystyle \frac{ϵ}{{ϵ}_p}}\right)}\right]}^C$(51)

They also defined the peak strain, maximum low stress, critical strain for initiation of DRX, and phase change DRX kinematics. The rule of mixture was used to define the interaction of the phases:

(52) ${\sigma }_Y=\{\begin{array}{c}{\sigma }_0\mathrm{if}ϵ<{ϵ}_{\mathrm{cr}}\\ \sum _{i=0}^{n-1}\left(X_i-X_{i+1}\right){\sigma }_i+X_n{\sigma }_n\mathrm{if}ϵ\ge {ϵ}_{\mathrm{cr}}\end{array}$(52) In another work, Bambach et al. [150] correlated the microstructure and flow stress of Ti-6AL-4 V processed by powder L-DED, LPBF, and the wrought material. The model considers the microstructure (α+β phases) and predicts the flow stresses for the individual processes. The flow stress of the β-phase is derived from the following equation: (53) ${\sigma }_{\text{β}}={\sigma }_{1\text{β}}-{\sigma }_{2\text{β}}$(53) ${\sigma }_{1\text{β}}$ and ${\sigma }_{2\text{β}}$ are defined as: (54) ${\sigma }_{1\text{β}}=\{\begin{array}{c}{\sigma }_p{\left[{\displaystyle \frac{ϵ}{{ϵ}_p}\exp \left(1-{\displaystyle \frac{ϵ}{{ϵ}_p}}\right)}\right]}^c\mathrm{for}ϵ<{ϵ}_p\\ {\sigma }_p\mathrm{for}ϵ\ge {ϵ}_p\end{array}$(54) (55) ${\sigma }_{2\text{β}}=\{\begin{array}{c}0\mathrm{for}ϵ<{ϵ}_p/\text{β}\\ \left({\sigma }_p-{\sigma }_s\right)\left[1-\exp \left(-k{\left({\displaystyle \frac{ϵ-{\displaystyle \frac{{ϵ}_p}{\text{β}}}}{{ϵ}_{\mathrm{ss}}-{\displaystyle \frac{{ϵ}_p}{\text{β}}}}}\right)}^n\right)\right]\mathrm{for}ϵ>{ϵ}_p/\text{β}\end{array}$(55) where ${\sigma }_p$ refers to the peak stress, ${ϵ}_p$ is the peak strain, ${\sigma }_s$ represents the steady-state stress, ${ϵ}_{\mathrm{ss}}$ is the steady-strain strain, and c, $k,n$, $\text{β}$ are constants. The flow stress of the α-phase was modelled based on the physical flow stress model implemented in [151]. The flow stress ${\sigma }_{\alpha }$ is defined by: (56) ${\sigma }_{\alpha }=M\left({\tau }^{\ast }+{\tau }_{\mu }\right)$(56) where $M$ represents the average Taylor factor, ${\tau }_{\mu }$ refers to the athermal long-range component of the flow stress which includes the effects of grain boundaries and stress field of dislocation forests, and ${\tau }^{\ast }$ is the thermally activated stress. The microstructure model was performed by defining the dislocation density, grain size, and globularization. The Hall-Petch coefficient $K_{\mathrm{HP}}$ was defined and the dislocation density was calculated using the Kocks-Mecking model, which considers dynamic recovery and work hardening. The complete definition of terms is found in [150].

Other investigations have been done on forming post-AM components using physically-based flow stress models. Galindo-Fernández et al. [152] developed a physically based model to describe the flow stress and yield strength of Ti-6Al-4 V with seven different microstructures (3 for cast and wrought components and 4 AM components). The model considered relevant microstructural features by linking the strengthening mechanisms, grain size, the former grain boundary morphology, and phase volume fractions. The yield strength was mainly determined by the relative fractions and morphologies of the α, ${\alpha }^{\mathrm{\prime }}$ and β phases. Hayes et al. [153] developed a physically based model to predict the yield strength of Ti-6Al-4 V across 3 heat treatments (stress relief annealing in the α+β phase field, hot isostatic pressing in the α+β phase field, and β-annealing) to an accuracy of 5%. A summary of studies using these models is presented in Table 5.

Table 5. Summary of physical models used in modelling forming operations of post-AM components.

References	Material	Globularization	Dislocaation	Rule of mixture	Grain size	Morphology	Flow stress
[149]	Inconel			X			X
[150]	Ti-6Al-4V	X	X	X	X		X
[152]	Ti-6Al-4V		X	X	X		X
[153]	Ti-6Al-4V		X	X	X		X
[115]	Ti-6Al-4V			X	X		X
[154]	Ti-6Al-4V		X	X	X	X	X

3.2. Machining of AM components

Subtractive manufacturing operations (machining) are usually performed after AM to achieve part accuracy and precision. However, machining may negatively affect the residual stress state on the part surface and reduce fatigue life [155]. Hence, FE simulations on the machining of AM components are required to develop strategies to minimise unfavourable residual stresses and geometrical deviations induced in AM and to study complete process chains that involve machining.

Salonitis et al. [156] predicted the residual stresses and geometry of a structural steel manufactured by the laser cladding process and subsequently processed by high-speed machining (HSP) in FE simulations. The thermal and mechanical loads imposed by the tool were neglected. Only elastic properties were assigned to the AM component. Machining was performed by deleting elements in the target area. The FE simulation results were not validated with experiments. Koenis et al. [157] performed EDM after wire L-DED to achieve part geometrical accuracy. The mesh of the final component was included in the original mesh and used throughout the process chain simulation to facilitate the transfer of previous simulation results by mapping the original mesh, thus limiting errors. Machining was performed as proposed by Salonitis et al. [156]. The process chain simulation results matched the temperature and deformation profiles measured in the experiments. Furthermore, Afazov et al. [158] developed a method for the precision manufacture of nuzzles with high geometrical accuracy and surface tolerance. The additive and subtractive processes were performed using element birth-and-death techniques. Additive-subtractive hybrid manufacturing was performed on thin-walled 316 steel parts [77,159,160]. Experiments and simulations revealed that machining induces compressive stresses, which reduce the residual stresses from the AM process. Luo et al. [160] present good agreements between temperature and deviation for AM experiments and simulations. The residual stresses post-milling mostly followed the same trend in the experiment and simulation; however, the forces in the milling operation were different in the experiment and FE simulation. The flow chart implemented to simulate the machining operation in He et al. [159] is presented in Figure 9.

Figure 9. Flowchart of additive manufacturing and machining FE simulations [159].

3.3. Performance analysis

Repeated heating in AM may alter the microstructure, cause phase transformation, and modify the crystallographic texture [161]. In some instances, moderate to superior mechanical properties have been observed in AM components compared to wrought materials. For instance, Spierings et al. [162] revealed that the fatigue life of SS316L produced by LPBF was 25% lower than the conventional material at low-stress values, and the fatigue life is comparable at high-stress. Bronkhorst et al. [163] showed that as-built AM 316L steel has 1.7 times the flow stress of wrought and heated-treated samples. The flow stress difference was attributed to the complex grain morphology and high initial dislocation density (2.5 times) in the as-built material compared to the wrought material. Riemer et al. [164] found that the fatigue performance of as-built 316L steel produced by LPBF compared to conventionally processed material.

Porosity, incomplete fusion, and other defects may negatively affect the performance of AM components. To utilise AM components in load-bearing applications, it is essential to understand their mechanical response under static and dynamic loading before they can be adopted in engineering material selection and processing. In addition, the effects of deposition strategies and other process parameters on the part’s performance need to be identified. Experiments and FE simulations have been performed to describe the deformation and failure of AM components.

3.3.1. Forming limits of AM components

The forming limits of AM components have been investigated using in-plane tests, uniaxial tensile tests, notched tension, equibiaxial tension, and shear tests. The lower bound of test results is often considered the forming limit of AM components due to the wide scatter of mechanical properties [165]. The mechanical properties of the components are usually defined by flow curves derived from empirical laws and yield criteria. The forming limits are defined by various ductile failure formulations. Stress-based forming limit diagrams have been investigated due to the dependency of failure on the stress state. A few studies are reviewed in the following paragraphs to indicate commonly implemented failure models used in forming limit predictions in FE simulations of AM components (see Table 6).

Table 6. Flow laws, yield criteria and failure models used in forming limit prediction of AM components.

References	Process	Material	Flow law + hardening type	Yield criterion	Failure model
[1]	LPBF	AlSi10Mg	Stochastic Voce hardening laws.	3D extension of the YLD2000–2D	Probabilistic Hosford-Coulomb fracture model.
[4]	LPBF	Stainless steel 316L	modified Johnson-Cook type plasticity model	Hill’48 non-associated flow rule	Hosford-Coulomb fracture initiation model.
[168]	LPBF	Ti-6Al-4V	associated flow rule+ swift hardening law	Hill’48 anisotropic yield
[170]	LPBF	Ti-6Al-4V	associated flow rule+ swift hardening law	Hill’48 anisotropic yield	Constant equivalent strain, Johnson-Cook fracture criterion, two-branch empirical fit fracture criterion, maximum shear stress fracture criterion, modified Mohr-Coulomb fracture criterion, Hosford-Coulomb fracture criterion.
[169]	LPBF	Austenitic stainless steel 316L	swift hardening law	Hill’48 anisotropic yield	Anisotropic Hosford-Coulomb fracture criterion, Anisotropic Mohr-Coulomb fracture criterion,
[166]	LPBF	Ti6Al4 V, 17-4PH and AlSi10Mg	$\sigma =K{ϵ}^n+\mathrm{Aϵ}$		Maximum Equivalent Strain (MES), Rice and Tracey (RT), Mohr and Coulomb (MC), Coppola Cortese and Folgarait (CCF).
[167]	BJT LPBF	Stainless steel 316L	Linear (Swift) and exponential (Voce)	Hill’s 48 non-associated flow role	Hosford-Coulomb ductile model.
[171]	LPBF	Ti6Al4V	Ramberg-Osgood	Von Mises	Rice and Tracey.
[172]	LPBF	Inconel 625	$\sigma =\mathrm{nA}{\left({ϵ}_0+\stackrel{.p}{\overline{ϵ}}\right)}^{n-1}$	Hill’s 48 model

Roth et al. [1] compared the microstructure and mechanical properties of AM (LPBF) and cast AlSi10Mg. FE simulations were performed with a stochastic Voce hardening law, a non-quadratic plasticity model with in-plane isotropy, and the probabilistic Hosford-Coulomb fracture initiation model. The ductility was lower for LPBF-printed parts. The proposed stochastic models can be used to assess the performance gap when replacing castings with LPBF parts. Along similar lines, Nalli et al. [166] investigated the ductile fracture of wrought and LPBF-produced AlSi10Mg, Ti-Al-4 V, and 17-4PH in different stress triaxialities. The maximum formability was calibrated for four failure criteria. The wrought parts had higher formability than the AM components, and formability increased for the machined AM parts.

Formability is strain rate-dependent for some materials. Li et al. [4] compared the rate-dependent plasticity of stainless steel 316L produced by the LPBF process and wrought material, see Figure 10. The Hill yield parameter and a combination of the Swift and Voce models were used to describe the hardening law in FE simulations to match the tension tests. Failure was defined by the Johnson and Cook failure model. Similarly, Tancogne-Dejean et al. [167] determined the rate-dependent mechanical properties and the fracture behaviour of binder-jetted stainless steel 316L. The modified Johnson-Cook flow stress model, Hill’s 48 non-associated flow role, and the Hosford-Coulomb ductile failure model were used to describe the stresses and strains in the tension experiments. Binder-jetted parts had equivalent plastic strains that were consistently lower than wrought (cold-rolled) material and LPBF parts. The low formability was explained by the porosity of parts, which was revealed by the fractured surfaces.

Figure 10. Comparison of the properties of additively manufactured (AM, solid lines) with cold-rolled (dashed lines) stainless steel 316L: (a) Engineering stress-strain curves as obtained from quasi-static uniaxial tension experiments at nominal strain rate of 0.001/s at 0°, 45° and 90° with respect to the building or rolling direction. (b) Hosford-Coulomb fracture envelopes for plane stress conditions in the equivalent plastic strain versus stress triaxiality plane and (c) fracture surface as a function of triaxiality and Lode angle parameter [4].

The build orientation and process parameters may cause anisotropy in AM components that influence flow stress and failure. Wilson-Heid et al. [168] proposed a plasticity model for Ti-6Al-4 V LPBF that predicts the anisotropy and stress state-dependent mechanical behaviour. The Swift law was used to describe the evolution of yield strength before necking. The Hill 1948 anisotropic yield function was used to model the multiaxial plastic behaviour. In another study, Wilson-Heid et al. [169] compared the accuracy of the anisotropic Hosford-Coulomb and anisotropic Mohr-Coulomb fracture criteria to predict the failure of LPBF austenitic stainless steel 316L in two orientations. Both models predict the shape of the forming limit diagram. However, the anisotropic Hosford-Coulomb fracture criterion is more performant because it is more sensitive to stress triaxiality and the Lode parameter, and it depends on all three principal stresses. Furthermore, Wilson-Heid and Beese [170] evaluated the performance of six failure models in predicting stress-state-dependent failure of LPBF Ti-6Al-4 V in two build orientations. The modified Mohr-Coulomb (MMC) failure criterion was the best model for predicting the forming limits; however, the MMC requires many tests for calibration. The maximum shear test can predict the failure strains using one test.

Some studies were conducted to investigate the effect of heat treatment on the formability of AM components. Yang et al. [171] characterised the forming limits of annealed and solution-aged Ti-Al-4 V produced by LPBF in different build directions. The annealed material had higher strength and better ductility than the solution-aged material. Similarly, Qin et al. [172] investigated the forming limits of LPBF Inconel 625 in the as-built and stress-relieved states. They showed that stress relief annealing decreased the fracture strain by 45% in the direction perpendicular to the build and increased the fracture strain by 36% in the build direction. The formation of the 10% δ -phase during stress relief annealing explains the difference.

3.3.2. Fatigue and creep limits of AM components

The mechanical properties like yield stress, ultimate tensile strength, and maximum elongation of AM components may be comparable to those of wrought materials [173]. However, process-induced defects usually reduce the fatigue and creep performance of AM parts. This arises from a heterogeneous microstructure distribution in different directions, randomly dispersed defects, geometry/size effects, etc., observed in AM components [174]. Furthermore, anisotropy and post-AM heat treatment may influence the fatigue behaviour of AM components [175]. Hence, estimating the reliability and lifetime of AM components becomes necessary for part safety and certification, especially in components subjected to oscillating stress and high temperatures like in-service aircraft components.

• Fatigue

Fatigue failure is caused by damage accumulation or subcritical crack growth that leads to part failure at stress states below the material’s tensile strength. Although the monotonic stress-stress properties of AM components may be superior to those of the rough part, the mechanical properties under cyclic loading need to be examined. Fotovvati et al. [176] reviewed the influence of the microstructure, residual stresses, and surface roughness on the fatigue performance of AM parts. The studies revealed that the fatigue crack growth properties of LPBF and wrought Ti-6Al-4 V are the same in the steady-state crack growth phase for heat-treated and hot isostatic pressed parts. The fracture toughness value can be improved by heat treatment. Similarly, Molaei et al. [177] presented a review of the fatigue performance of AM Ti-6Al-4 V considering the effects of build orientation, process parameters, fabrication machine, heat treatment, and surface condition. Torries et al. [178] presented models to predict fatigue using microstructure and defect-sensitive models. The review contains mostly experimental works, although a sub-section is dedicated to FE-based fatigue prediction. In this section, we present recent models on FE-based fatigue prediction of AM components, emphasising physically based models and the transferability of results from test coupons to AM components.

Lindström et al. [179] proposed a model to simulate the mid-life behaviour of a nickel-based superalloy subjected to thermomechanical fatigue and rate-independent low-cycle fatigue as in gas turbines operating in high-temperature environments. Mooney et al. [180] demonstrated that the anisotropic Hill model was more accurate in describing the cyclic loading than the isotropic plasticity model. Kourousis et al. [181] compared the performances of four advanced rate-independent plasticity kinematic models (multicomponent Armstrong-Frederick, multicomponent Armstrong-Frederick with a threshold term, Ohno-Wang, and multipoint Armstrong-Frederick with a multiplier) in modelling the cyclic elastoplastic behaviour of LPBF Ti-6Al-4 V. All four models demonstrate satisfactory performance. Lindström et al. [182] estimated the fatigue life of LPBF nickel-based superalloys subjected to isothermal cyclic loading at different temperatures using the Ohno-Wang model developed to describe cyclic elastoplastic midlife fatigue behaviour. A cycle-jumping procedure is used to simulate mid-life behaviour.

The geometric dependency of fatigue has been investigated. Zhang et al. [183] demonstrated using FE simulations of WA-DED parts that the residual stresses near the crack tip in compact tension specimens affect the crack growth rate. High stress intensity was recorded at the crack tip. Huynh et al. [184] showed that microtruss structures have lower fatigue performance compared to dog bone specimens produced by direct metal laser deposition of IN718 under the same process conditions.

Li et al. [185] introduced probabilistic dependency to the fatigue limit prediction of AM components to model the dispersion observed in the failure of AM components. Trends in the median fatigue performance of test coupons were transferred to AM components for fatigue prediction. Hypothetically increasing crack initiation sites in the fatigue model of the material by 200% led to a decrease in the median fatigue performance of the test coupon and component by 10% and 5%, respectively. The material changed from a homogenous stress state to a heterogeneous stress state for the coupon and component. This indicates that the results of test coupons can be implemented in designs after considering the size effect. Arbeiter et al. [186] proposed a proof-of-concept to enable rapid, easy, and reliable estimation of the lifetime of AM components from test coupons. FE simulations and fatigue tests were performed to estimate the lifetime of AM coupons under different fatigue conditions. The fatigue life is estimated based on the stress intensity factor $K$, defined as: (57) $K=\sigma \cdot Y\cdot \sqrt{\pi \cdot a}$(57) where $\sigma$ refers to the applied nominal stress, $Y$ defines a geometry dependent function, and $a$ is the crack length. The Paris’ law was used to describe crack advancement ($\mathrm{da}$) for an increment of cycles ($\mathrm{dN}$) as a function of the stress intensity factor ($K_I$, ‘I’ for tension opening mode). (58) ${\displaystyle \frac{\mathrm{da}}{\mathrm{dN}}=A\cdot (K_I{)}^m}$(58)

The estimated number of cycles before failure is given by Paris’ law. (59) $N_{\mathrm{prop}}={\int }_{a_{\mathrm{ini}}}^{a_f}{\displaystyle \frac{\mathrm{da}}{A\cdot {\left[K_I\left({\displaystyle \frac{a}{W}}\right)\right]}^m}}$(59) where $A$, $m$ and $K_I$ are dependent on the geometry $\left(a/w\right)$. $a_{\mathrm{ini}}$ and $a_f$ define the initial and final crack lengths measured postmortem on the fractured surface. The approach provides a means to upscale fatigue results from test coupons to AM components. Although the methodology was applied to the lifetime prediction of an AM-processed wrench made from poly (lactic acid), it can be adapted for AM of metallic components. A schematic representation of the procedure is presented in Figure 11.

Figure 11. (a) Fracture mechanics-based lifetime estimation of an arbitrary component based on material parameters. (b) Finite element analysis [186].

Physical models have been implemented to predict fatigue. Torries et al. [187] proposed a microstructure-sensitive fatigue model for powder L-DED Ti-6Al-4 V based on pore size, second-phase particle size, grain size, porosity, and other factors. The microstructure-sensitive fatigue model was calibrated for the average size of the inclusions and fitted for as-built and heat-treated specimens. Rutherford et al. [188] quantified the relationship between the microstructure and mechanical properties (tensile and strain-control fatigue) of AA6061 processed by additive friction stir-deposition. The component’s fatigue life was expressed as a sum of three parts: (1) The crack incubation stage (inc). (2) Microstructurally small crack growth (MSC) and physically small crack (PSC) growth. (3) Long crack (LC) growth. (60) $N_{\mathrm{total}}=N_{\mathrm{inc}}+N_{\mathrm{MSC}/\mathrm{PSC}}+N_{\mathrm{LC}}$(60)

The incubation stage ($N_{\mathrm{inc}}$) is defined by the modified Coffin-Manson expression as: (61) $C_{\mathrm{inc}}{N}_{\mathrm{inc}}^{\alpha }={\displaystyle \frac{\mathrm{\Delta }{\gamma }_{\max }^{p^{\ast }}}{2}=\psi \overline{\mathrm{{\rm Y}}}[{ϵ}_a-{ϵ}_{\mathrm{th}}{]}^q}$(61)

The microstructurally-physically small crack stage (62) ${\left({\displaystyle \frac{\mathrm{da}}{\mathrm{dN}}}\right)}_{\mathrm{MSC}}=\chi \left(\mathrm{\Delta CTD}-\mathrm{\Delta CT}{D}_{\mathrm{th}}\right),a_i=0.625D$(62)

The long crack (LC) growth was modelled based on linear elastic fracture mechanics. Using a physically based lifetime estimation, Blinn et al. [175] investigated the effect of heat treatment on the fatigue performance of LPBF AISI 316L stainless steel. They estimated fatigue caused by surface, subsurface, and internal defects.

• Creep

Models have been developed to describe the influence of heat treatment and high temperatures on the creep performance of AM components. Pröbstle et al. [189] revealed that LPBF Inconel 718 shows superior creep strength in all heat treatment conditions compared to conventional cast and wrought materials. Higher volume fractions of the strengthening phases and process-induced sub-grains lead to superior creep of the LPBF material. Kim et al. [190] controlled the microstructure of the LPBF Ti-Al6-4 V alloy to improve high-temperature creep properties. The $\alpha /\text{β}$ interface acts as an obstacle in dislocation motion during creep, thereby increasing the creep resistance. The heat-treated specimen had better creep resistance than the as-fabricated specimens due to the development of the Widmanstätten structure. The creep time was estimated by the Garofali-Dorn model, and the Sherby-Dorn equation calculated the steady-state creep rate. (63) ${\dot{ϵ}}_{\mathrm{SS}}=K{\sigma }^n\exp \left({\displaystyle \frac{-Q}{\mathrm{RT}}}\right)$(63) where $K$ represents the Dorn constant, $\sigma$ refers to the applied stress, $n$ is the stress exponent,$Q$ describes the activation energy for creep,$T$ is the ideal gas constant and $R$ is the gas constant. A similar formulation was employed by Uzan et al. [191] in the LPBF of AlSi10Mg. Cardon et al. [192] added more parameters to the model and included temperature dependence in steady-state creep rate prediction. He et al. [193] examined the creep performance of LPBF Ti6Al4 V at various temperatures and build directions. The creep property was worse at 0° compared to 45^° and 90^°. The LPBF samples had a longer creep rupture time than the forged samples but lower ductility. Furthermore, they modified the Wen-Tu creep model to describe the anisotropic behaviour: (64) ${\dot{ϵ}}^c={\displaystyle \frac{3}{2}A{\sigma }_e^{n-1}\exp \left[{\displaystyle \frac{2\left(n+1\right)}{\pi \sqrt{1+3/n}}{\omega }^{3/2}}\right]t^{m}M:s}$(64) (65) ${\sigma }_e={\left({\displaystyle \frac{3}{2}s:M:s}\right)}^{{\displaystyle \frac{1}{2}}}$(65) (66) $\dot{\omega }={\dot{ϵ}}^c/{ϵ}_f$(66) (67) ${\displaystyle \frac{{ϵ}_f^{\ast }}{{ϵ}_f}=\exp \left[{\displaystyle \frac{2}{3}\left({\displaystyle \frac{n-0.5}{n+0.5}}\right)}\right]/\exp \left[2\left({\displaystyle \frac{n-0.5}{n+0.5}}\right){\displaystyle \frac{{\sigma }_m}{{\sigma }_e}}\right]}$(67)

where $A$, $n$, $m$ are materials constants, $M$ refers to the anisotropic tensor, $t$ represents the creep time, ${\sigma }_e$ is the equivalent stress, $s$ is the deviatoric stress tensor, $\dot{\omega }$ is the damage rate, ${\sigma }_m$ is the hydrostatic stress, ${\dot{ϵ}}^c$ is the creep strain rate, ${ϵ}_f^{\ast }$ is the multiaxial failure strain, and ${ϵ}_f$ is the uniaxial creep ductility.

Physically-based constitutive models have been developed for creep prediction. Paoletti et al. [194] described the creep of LPBF AlSi10Mg and considered the coarsening and ripening phenomena of Si-particles. The time for the minimum creep rate of the as-built specimen was estimated by (68) $d^3=d_0^3+K_g{D}_{\mathrm{pL}0}{t}_m\exp \left({\displaystyle \frac{Q_{\mathrm{Ls}}}{\mathrm{RT}}}\right)\left[1+10\exp \left({\displaystyle \frac{Q_{\mathrm{Ls}}-Q_{\mathrm{ps}}}{\mathrm{RT}}}\right)B{\dot{ϵ}}_m\right]$(68) The model considers Ostwald ripening phenomena that describe Si particle evolution during annealing and the accelerating effects of strain on the progressive softening of the alloy. The material model was described by these sets of equations: (69) $\sigma ={\sigma }_0+{\sigma }_p={\sigma }_0+\mathrm{\alpha mGb}\sqrt{\rho }$(69) (70) ${\dot{ϵ}}_m={\displaystyle \frac{2M_{\mathrm{cg}}{\tau }_l\mathrm{bL}}{m}{\left({\displaystyle \frac{{\sigma }_{\rho }}{\mathrm{\alpha mGb}}}\right)}^4}$(70) (71) $M_{\mathrm{cg}}\cong {\displaystyle \frac{D_{0L}b}{\mathrm{kT}}\exp \left({\displaystyle \frac{{\sigma }_{\rho }{b}^3}{\mathrm{kT}}}\right)\exp \left\{-{\displaystyle \frac{Q_L}{\mathrm{RT}}\left[1-{\left({\displaystyle \frac{{\sigma }_{\rho }}{R_{\max }}}\right)}^2\right]}\right\}\exp \left({\displaystyle \frac{U_{\mathrm{ss}}}{\mathrm{RT}}}\right)}$(71) (72) $R_{\max }=1.5\left(R_{\mathrm{UTS}}^a+{\sigma }_{\mathrm{Or}}\right){\displaystyle \frac{G_T}{G_{\mathrm{RT}}}}$(72)

Furthermore, Paoletti et al. [195] applied the model to predict the effect of low-temperature annealing (225 °C/2 h) on the creep response of heat-treated LPBF AlSi10Mg and compared the mechanical response of the as-built material and heat-treated specimens. Kim et al. [196] proposed another set of physically based creep prediction models to describe the compressive creep properties of an LPBF-built equiatomic CoCrFeMnNi high-entropy alloy. They employed the Sherby-Dorn equation to describe the strain rate in steady-state creep (${\dot{ϵ}}_{\mathrm{ss}}$). They estimated the activation volume by Equation (73) to have a better insight into the creep mechanism: (73) $V^{\ast }=-\mathrm{MkT}\left({\displaystyle \frac{d\left(\ln {\dot{ϵ}}_{\mathrm{ss}}\right)}{\mathrm{d\sigma }}}\right)$(73)

where $M$ represents the Taylor factor and $k$ is the Boltzmann’s constant. They further stated that viscous glide is responsible for creep deformation in the material at low strain rates and dislocation climb at high strain rates. The viscous glide was modelled by (74) ${\dot{ϵ}}_{\mathrm{SS}}={\displaystyle \frac{\pi \left(1-v\right)}{6}{\left({\displaystyle \frac{\mathrm{kT}}{e^{2}c}}\right)}^2\left({\displaystyle \frac{\tilde{D}\mathrm{Gb}}{\mathrm{kT}}}\right){\left({\displaystyle \frac{\sigma }{G}}\right)}^3}$(74) where $v$ refers the Poisson’s ratio, $e$ represents the solute-solvent size difference, $c$ is the solute concentration, $\tilde{D}$ is the chemical diffusivity of the solute atoms, $b$ is the Burgers vector and $G$ is the shear modulus.

4. Studies on AM process chain simulation

FE simulation of AM is still maturing, and far less research has addressed the process chains involving AM and post-processing of AM components [157]. Process chain simulations provide parametric-associative modelling and knowledge integration that can assist engineers in designing components and assemblies, thus supplying information flow within manufacturing process chains [197]. Current research primarily focuses on low-fidelity and fast simulation approaches that ignore AM process details such as molten pull fluid flow, phase change, etc. The same mesh is usually used throughout process chain simulations to limit errors when transferring results between processes [198]. Afazov et al. [199] developed a Finite Element Data Exchange System (FEDES) for manufacturing process chain simulation. FEDES can map FEA data between six commercial FE packages and different mesh types generated by 2D and 3D linear or quadratic element types. Studies that include AM process chain simulations and welding include the following:

Fisk and Lundbäck [127] predicted the residual stresses of Inconel (Alloy 718) after welding, followed by global and local heat treatment using physically based flow-stress data defined by the dislocation density. Deepu and Phanikumar [97] performed process chain simulations to predict the temperature and microstructure evolution in gas arc metal welding of DP980 and performed post-welding tensile tests. Ahn et al. [93] demonstrated that implementing phase transformation models had a marginal influence on the residual stress prediction and distortion in fibre laser welding of Ti-6Al-4 V. They performed post-welding heat treatment to reduce the tensile residual stresses observed at the HAZ.

Studies on process chain simulations of AM components include: Zhang et al. [115] modelled the LPBF of Ti-6Al-4 V and performed tensile tests on the AM specimens. The yield strength of the component was modelled based on a Hall-Petch relationship that is mainly dependent on the martensitic ${\alpha }^{\mathrm{\prime }}$ and ${\alpha }_w$ phases. After the LPBF simulation, the part with its local mechanical properties was imported into a mechanical solver to conduct tensile tests. The yield and ultimate stresses were converted to true stress and strains by fitting the Ramberg-Osgood power-law equation. Koenis et al. [157] developed a virtual framework to perform process chain simulations of AM components to determine the influence of all the steps in the manufacturing chain on the final product (see Figure 12). The final part dimension was obtained by EDM because EDM has a minimal effect on the surface residual stresses. The virtual framework thus provides a fast method to determine the influence of the various operations in a process chain.

Figure 12. Evolution of the process induced residual stress field for different process chains obtained by the process chain simulations. [157].

The inherent strain method that directly applies strains in the FE model was used in AM process chain simulation instead of the traditional method of predicting process temperature and using a thermal expansion model to compute thermal strains [200]. O’Brien et al. [201] performed complete manufacturing process chain simulations of an impeller made from Inconel 718 using the inherent strain method. The process chains involved LPBF, heat treatment (HT), shot peening (SP), or laser shock peening (LSP). The inherent strain method is used in simulating the LSP. The Johnson-Cook flow stress model was used in the LSP and SP simulations. Similarly, Yaghi et al. [202] adopted the inherent strain method to simulate LPBF. Distortion and stresses in AM and post-AM machining were predicted by FE simulations and verified with experiments.

The feasibility of utilising process chain simulation in industrial settings was analysed to reduce the experimental effort in AM product design. Yaghi et al. [203] assessed the maturity of modelling and simulation methods used in AM (LPBF) based on the Modelling and Simulation Readiness Level (MSRL) scale derived from the Technology Readiness Level (TRL) of NASA. Seven relevant categories and maturity factors were considered and analysed on a scale ranging from level 0 (no capability exists) to level 5 (operational), as shown in Figure 13. The processes in the manufacturing process chain were: LPBF, machining, heat treatment, laser shock peening, or shot peening. The overall MSRL score of the process chain simulation was 3 (rounded from 2.7), which indicates that the maturity is between ‘applied R&D’ and ‘instrumental.’ This means that LPBF process chain modelling and simulation are not common in industry and are mostly performed by academia and research and technology organisations (RTOs). In addition, ‘LPBF process chain modelling can support low-risk development with a high reliance on simulation experts.’ Thus, FE simulations can be used for qualitative assessment of process chains, defining design structure, constraining design space, ranking alternative designs and solutions, and replacing some experimental trials.

Figure 13. MSRL of the reported LPBF process chain [203]

Afazov et al. [22] presented an overview of the modelling techniques used in PBF process chains involving AM, heat treatment, machining operations, and surface engineering. They highlighted the need to validate process chain simulations with experimental results and develop accurate models that predict the mechanical response of macro-FE simulations to reduce errors. A summary of studies with process chain simulations involving AM and welding are presented in Table 7.

Table 7. Summary of process chain simulations involving AM and welding.

References	Process chain					Material models
	AM process	Heat treatment	Forming	Machining	Fatigue/Creep
[157]	wire L-DED			Wire-EDM	Fatigue	Temperature dependent thermal and mechanical features. Norton-Bailey creep
[159]	LPBF			Milling
[160]	LPBF			Milling
[77]	LPDB			Milling
[203]	LPBF	Heat treatment		laser shock peening/ shot peening
[115]	LPBF		Tensile test			Physically based flow stress model with Hall-Petch Phase JMAK model. Ramberg-Osgood power-law equation
[93]	Laser welding	Heat treatment				Temperature dependent thermal and mechanical features KM and JMAK phase change model.
[97]	GMAW		Tensile test
[156]	Laser cladding			High speed machining		Temperature dependent thermal and mechanical features.
[22]	LPBF Electron beam PBF	Heat treatment		Machining Surface engineering		Mechanical features
[158]	LPBF			Micro-milling
[201]	LPBF	Heat treatment		laser shock peening/ shot peening		The inherent strain method. The Johnson-Cook flow stress.
[202]	LPBF			Machining		Mechanical features
[127]	GTAW	Heat treatment				Physically based flow stress model

5. Discussion and outlook

The macro-scale FE simulations presented in this review provide information on the global stress state, temperature distribution, component geometry, phase transformations, etc. However, including the process details is required to improve the prediction accuracy of the simulations. In this chapter, the challenges to macroscale simulations are out lined. This is followed by future development directions and an outlook on performing process chain simulation and component performance analysis.

5.1. Challenges to the implementation of macro scale process chain simulation frameworks

Numerical simulations of traditional manufacturing processes have been investigated for several decades and can be considered mature. In this section, mainly the challenges to macroscale AM simulation are discussed and how simulation results can be linked to other simulations to form process chains. The challenges are summarised as phenomenological and methological in Figure 14.

Figure 14. Challenges to the implementation of macroscale AM process chain simulations

Phenomenological challenges

• The phenomena in AM simulation occur at multiple scales, as highlighted in section 1.1. Some phenomena occurring at lower length and time scales cannot be replicated at the part scale, which reduces AM simulation accuracy. For example, a deeper understanding of powder properties, a description of properties at the molten state, and the interactions of the solid–liquid-gas in AM are required to improve simulation accuracy.

• The fact that some events in AM are transient and occur for very short periods of time makes it challenging to replicate them in part-scale simulations.

• AM processes have inherent variation, even for the same set of parameters. This stochastic nature increases the effort needed to model the processes.

Methodological challenges

• Complete material properties are often rare to find in databases due to cost factors. In addition, material data in the simulations is often derived from the literature, which may be different from the investigated material. Furthermore, simulations are usually performed using the properties of the bulk material in powder-based processes, which induce errors due to variations in bulk and powder material properties.

• Highly qualified personnel are required to implement process chain simulation in industry, which limits their application potential. However, the high interest and progress in AM process modelling promise better integration in future. The integration of AM in process chains can be facilitated by developing AM modules in existing CAD/CAM and FE simulation packages.

• Powder-based AM processes are modelled as a continuum solid with isotropic properties. Such oversimplifications, the multi-physics nature of AM, and a limited understanding of the intricate evens in AM processes reduce accuracy. In addition, the effects of gravity on the displacement and distortion of the molten metal are often omitted.

• The actual scanning or material deposition path is often replaced by complete deposition layers to reduce computational time in AM simulations, which may lead to errors in phase prediction.

• Heat loss from the free surfaces to the surrounding area caused by thermal radiation and conduction at the bottom surface of the substrate is often neglected, which reduces the accuracy of simulations. In addition, macro-scale models usually ignore the details of the heat focus, mass transfer, the role of powders in deposition, porosity, material interactions, and fluid flow in the fusion zone during AM.

• Finding and coupling interactions between the multi-physics problems in AM is challenging.

• Software incomparability may result from using different software packages for modelling the processes in the chain. In addition, transferring simulation results to other process simulations may cause problems due to the need to perform remeshing. This can be eliminated by using the same software package or by developing platforms for the exchange of data between softwares.

5.2. Future development directions in AM process chain simulations

Linking data from micro, meso, and macroscale FE simulations in an AM ICME process chain framework remains a major challenge since many process parameters may not be within the material calibration domain and the parameters may have different time and length scales [204]. In addition, a huge amount of experimental and computational effort is required to develop ICME models. Since simulating at lower scales is currently very time-consuming, it seems more feasible to perform mainly macroscale simulations. The studies with process chain simulations presented in section 4 mainly employ empirical models and omit phase transformation and microstructure. In this section, we propose some strategies to enhance the results of macroscale simulations and achieve the goal of process chain simulation.

5.2.1. Implementing physically based material models

Physically based models offer the possibility of describing process-structure–property-performance relationships for AM components at the macroscale. Such models can describe variations in the microstructural properties with changes in temperature and stress. Some researchers have addressed this by developing multiscale, physically based models for AM components and also studying the influence of post-AM processes [205–208]. However, these models are often complex and have several variables. It is required to identify important parameters and reduce variables [209].

5.2.2. Including mesoscale features in macroscale FE simulations

Some features that are typically observed in mesoscale simulations can be simplified and investigated with macroscale FE simulations. For example, in powder melting in LPBF, a statistical distribution of the powder is often used to describe the area of interest rather than the actual distribution of the powder in the process [22]. Pores that can be observed in mesoscale simulations can be detected by X-ray computed tomography (X-ray CT) and can be included in the mesh of macroscale simulations [210] [211].

5.2.3. Integrating artificial intelligence and big data

Many parameters and interactions between parameters influence the properties and performance of AM components in a process chain. Artificial intelligence models can identify and predict the properties of such complex relationships. Machine learning models should be trained based on the inputs, process parameters, and outputs of FE simulations and experiments to identify the relevant features required to predict component properties and performance. Predictions from a process should be robust before being used as inputs for subsequent predictions. The data from the previous processing stage is indicated in red in Figure 15. The FE simulations and experimental results can also be integrated into digital manufacturing tools and cloud platforms to create digital production chains.

Figure 15. The principal inputs and outputs from AM process chain simulations

5.3. Outlook on process chain simulation and performance analysis

A framework for process chain simulation that considers the influence of AM on the component microstructure is proposed in Figure 16. In the framework, FE simulations of individual processes (Pn) in the process chain are performed and calibrated to the results of experiments. FE simulations of subsequent processes use the outputs of the previous FE simulation as inputs until the final component (${\mathrm{P}}_{\mathrm{End}}$) is obtained. A material modelling approach that meets the requirements is described in [111]. The thermomechanical modelling of the flow curve is based on the phase fractions, which calculates the yield stress using the mixing rule. The special feature is that all the usual phase transformation equations and the phase decomposition kinetics arising under reheating are implemented in an extensible manner. The kinetics as well as the yield stresses themselves should be described depending on superordinate defined parameters such as dislocation densities or grain sizes. To establish the complete process chain simulation, the feedback loop should be executed and boundary conditions adapted until the FE model reasonably reproduces the final properties of the AM component obtained in experiments (microstructure, fatigue, forming limits, etc.). The outputs from AM simulation are required for associative modelling and knowledge integration in component and process design. The process chain simulation sequence permits process control, optimisation, inverse simulation of experiments and a reliable tool for predicting manufacturing quality and can serve in digital twin frameworks. An additional facet of future simulation will include the virtual comparison of various production chains via the use of multifunctional benchmark tests [111, 212].

Figure 16. Process-chain simulation with AM indicating the material properties, process parameters and evaluation.

Advances in phase transformation modelling and microstructure prediction in macroscale AM simulations are promising. As such, AM simulations should be included in subsequent process simulations to achieve process chain simulation.

6. Conclusions

Process chain simulation is required to facilitate the design of AM components. However, most research studies separate AM process simulation from the simulation of post-AM components. This review aimed to present the state of the art in macro-scale AM process chain simulation and explore means for AM process chain simulation. The following conclusions were drawn from the review:

• Research should be dedicated to developing physics-based and phase change models that describe the microstructure evolution of materials at different process conditions rather than using empirical models that have no physical meaning. Physics-based models permit the transfer of microstructure properties through process chain simulations. The development of physical models is also required to simulate process chains without AM.

• Given the large amount of available simulation and experimental data for some materials, like Ti-6Al-4 V, Inconel 718 and their industrial relevance, AM process chain simulations should be conducted and validated to benchmark standards for process simulation of such materials.

• Physically based models often contain many parameters, which, in combination with the multiple forming steps, would create a huge set of material input parameters. By employing parameter identification and reduction techniques, the number of input parameters can be reduced while maintaining an acceptable level of accuracy.

• Pores and other defects observed in experiments and mesoscale simulations can be calibrated in the performance evaluation of macroscale FE simulations to reduce the experimental and numerical effort.

• AM is characterised by randomness. Stochastic parameters should be included in AM modelling to describe the randomness.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.