# Metallurgically-driven thermomechanical analysis of multiple side-to-side laser melting on a 316L substrate

Zixuan Li, Michel Bellet, Charles-André Gandin, Manas Vijay Upadhyay, Yancheng Zhang

► To cite this version:

Zixuan Li, Michel Bellet, Charles-André Gandin, Manas Vijay Upadhyay, Yancheng Zhang. Metallurgically-driven thermomechanical analysis of multiple side-to-side laser melting on a 316L substrate. Additive Manufacturing, 2025, 112, pp.104991. 10.1016/j.addma.2025.104991. hal-05326559

HAL Id: hal-05326559

https://hal.science/hal-05326559v1

Submitted on 4 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# Metallurgically-driven thermomechanical analysis of multiple side-to-side laser melting on a 316L substrate

Zixuan Li<sup>a</sup>, Michel Bellet<sup>aID,*</sup>, Charles-André Gandin<sup>aID</sup>, Manas Vijay Upadhyay<sup>b</sup>, Yancheng Zhang<sup>aID</sup>

## ARTICLE INFO

**Keywords:**

Grain structure

Crystal plasticity

Finite elements

Residual stress

Additive manufacturing

## ABSTRACT

In additive manufacturing, the solidification grain structure has a significant influence on the properties of as-built material. In this context, the solidification grain structure and internal stress evolution during laser scanning of polycrystalline 316L stainless steel are simulated. A strongly coupled crystal viscoplasticity model is developed and integrated with a cellular automaton–finite element (CAFE) approach to accurately capture grain structure and stress evolution, where the CAFE model is validated based on a literature experiment. The crystal viscoplasticity model is calibrated using stress–strain curves of annealed 316L from experiments considering small thermo-elasto-viscoplastic (TEVP) deformations. The resolution algorithm dynamically couples heat transfer, melting and solidification simulations while concurrently computing stress and strain evolution within the grain structure. Four scanning strategies are simulated using the coupled CAFE–crystal viscoplasticity approach, capturing stress evolution during grain growth. This enables the simultaneous thermo-viscoplastic modeling in the mushy zone and TEVP modeling in the solid, providing insights into stress evolution and grain orientation over a large domain. The melting-solidification process involves variations in compression and tension, leading to stress concentration within neighboring grains with significant orientation differences, extending along elongated grains. A framework for multiscale process-structure-mechanical investigation is established based on microscale stress evolution in additive manufacturing.

## 1. Introduction

Laser Powder Bed Fusion (L-PBF) consists of high-frequency cycles of melting and solidification [1], which result in extremely high-temperature gradients (10–20 K µm⁻¹) [2] and cooling rates (10²–10⁶ K s⁻¹) [3] at the periphery of the melt pool. Such cycling, with track-by-track repetitions, gives rise to heterogeneous stresses within a non-equilibrium microstructure, leading to defects, particularly cracking that includes solidification cracking [4] in the melt zone and cold cracking [5] in the solid state. The microstructure features play a critical role in the initiation and propagation of solidification cracking [6–8]. Grain size and misorientation influence crack propagation, with high-angle grain boundaries (HAGB) promoting cracks [9]. In addition, the crystalline materials exhibit anisotropic properties due to the growth of columnar dendrites-based microstructure. To address the relationship between stresses and grain growth in L-PBF, it is crucial to gain a comprehensive understanding of their origin at the grain level and their potential impact on the anisotropic mechanical behavior [10].

Identifying critical process parameters to minimize defect formation in L-PBF remains challenging due to the complexity of multi-physics phenomena occurring over short timescales. Experimental methods, often trial-and-error, are costly and limited by the inability of *in situ* techniques to simultaneously capture microstructure and stress evolution. In addition, grain structure analysis via scanning electron microscopy (SEM) or electron backscatter diffraction (EBSD) [11] is mostly limited to two-dimensional measurements and cannot operate simultaneously during the solidification process. Therefore, numerical modeling has emerged as a cost-effective alternative, enabling simulations at the track level to predict stress distribution and grain structure evolution.

The metallurgical simulation of solidification grain structures aims to enhance our understanding of its relation with defect formation. Given the complexity of dendritic morphology, models often avoid detailed descriptions of the solid–liquid interface, focusing instead on the grain envelope [12]. This is found in the Monte Carlo method [13–16] and in the cellular automaton (CA) method used for simulating

<footnote><sup>a</sup> MINES Paris, PSL University, Centre de Mise en Forme des Matériaux (CEMEF), CNRS UMR 7635, 1 Rue Claude Daunesse, Sophia Antipolis, 06904, France</footnote>

<footnote><sup>b</sup> Ecole Polytechnique, Institut Polytechnique de Paris, Laboratoire de Mécanique des Solides (LMS), CNRS UMR 7649, Route de Saclay, Palaiseau Cedex, 91128, France</footnote>

<footnote>* Corresponding author.
E-mail address: michel.bellet@mines-paristech.fr (M. Bellet).</footnote>
additively manufactured alloys [17–29]. The latter is believed to be more physics-based as it makes use of growth kinetics linked to the local temperature [30–32] and has been shown to replicate texture formation during directional solidification [33,34]. One could also notice tentative modeling using the phase field method applied to additive manufacturing, yet limited to small domain size [35–37]. In particular, Chadwick et al. [37] performed a one-to-one comparison of their phase-field model predictions with laser scanning experiments performed using a unique coupling between a continuous-wave laser and a scanning electron microscope (CWLASER-SEM) [38]. Some experimental results from this work are also used to validate the model proposed in this work.

For thermo-mechanical simulations, the finite element method (FEM) has been widely used to simulate the construction process at the part scale. For mesoscale models, appropriate material constitutive law for high temperature thermo-mechanics should be used. Beghini et al. [39] integrated a thermo-elasto-viscoplastic (TEVP) model to assess material stress and microstructure evolution in mesoscale during directed energy deposition processes. Darabi et al. [40] coupled the thermo-elasto-plastic model with a phase-field analysis of fracture considering shrinkage, thermal strain, and the rapid cooling stage in solidification crack evolution.

To address the anisotropic material behavior and the grain morphology, the crystal plasticity (CP) finite element method has emerged as a suitable approach, as demonstrated in recent studies. Based on thermal-fluid flow and grain growth simulations, Hu et al. [41,42] exploited the temperature distribution and the grain structure in a CP model using an element elimination technique and a reactivation method to model melting and solidification, as well as to reinitialize state variables. Lindroos et al. [43] introduced a 2D thermo-mechanical CP model aimed at scrutinizing the thermal and microstructural evolution of the melt pool area, alongside residual stresses, within the context of single-track scanning. This analysis takes into account phase transitions occurring in the solid state. Mohanan et al. [44] combined simulation and laser scanning experiments to analyze intergranular stresses and plastic strains in stainless steel by investigating the impact of microstructure refinement and plastic heterogeneity through a TEVP CP model. They also performed laser scanning experiments using the CWLASER-SEM device [38] and validated their model via a statistical comparison of the predicted Nye's tensor with the experimentally measured one from electron backscatter diffraction data. To predict residual stress and strain evolution at the microstructure resolved polycrystalline length scale, Kuna et al. [45] developed a framework to directly model the effect of processing parameters used during LPBF on the thermomechanical response of the material during solidification. Model validation was performed on various features of intragrain misorientation compared to electron backscatter diffraction data of a 316L thin wall geometry. In addition, correlations between processing parameters, temperature field, grain texture and the resulting as-built mechanical response were determined and discussed.

In light of the above, the present work uniquely focuses on laser scanning of a substrate under process conditions without a powder bed, an essential but often overlooked step toward the L-PBF process. Unlike previous studies, this work couples the CAFE framework with crystal plasticity, enabling the simultaneous prediction of stress evolution during grain growth in the mushy zone and in the solid state. The paper is structured as follows: Section 2 describes the model, focusing on the CAFE framework, on the coupling strategy with the CP method, and the rules to be applied upon transition between solid and liquid. The resolution of global and local non-linear systems is also be discussed. Model calibration of the parameters, experimental validation using the results from [37], and the basic setting of calculation are presented in Section 3. Section 4 discusses the application of the model to multiple passes to reveal grain growth and the stress evolution at the grain structure level as a function of the process parameters. Section 5 summarizes the results and discusses the potential for future insights into crack formation due to temperature gradients and cooling processes. The alloy considered in this study is the 316L stainless steel.

## 2. Numerical model description

The present model is based on the FE solution for heat transfer and the CA method for the evolution of grain structure. Mechanical simulation is carried out simultaneously using the CP method, incorporating grain structure and thermal data through coupling with the CAFE framework.

### 2.1. Thermal modeling

The temperature field is calculated by the finite element method, through the resolution of the following equation for energy conservation, incorporating phase changes:

$$ \frac{\partial \langle \rho h \rangle}{\partial t} + \nabla \cdot (\langle \rho h \rangle \langle v \rangle) - \nabla \cdot (\langle \lambda \rangle \nabla T) = 0 \quad (1) $$

where $\langle \cdot \rangle$ means the spatial averaging of quantities over all phases, $T$ [K] is the temperature, $t$ [s] is the time, $\rho$ [kg m⁻³] is the density, $\lambda$ [W m⁻¹ K⁻¹] is the thermal conductivity, $v$ [m s⁻¹] is the velocity and $h$ [J kg⁻¹] is the specific enthalpy including solid and liquid therefore including the latent heat of fusion. The net input power includes the input power and the heat loss at the system boundaries. An appropriate heat source model is required to simulate the energy applied by the laser: a Gaussian source term $\dot{Q}_{laser}$ [W m⁻²] is considered at the top surface of the domain $\Omega$:

$$ \dot{Q}_{laser}(\mathbf{x}) = \eta \frac{2P_{laser}}{\pi r_{laser}^2} \exp \left( -2 \frac{(\mathbf{x} - \mathbf{x}_{laser})^2}{r_{laser}^2} \right) \quad (2) $$

where $\eta$ is the laser absorptivity of the surface, $P_{laser}$ [W] is the power of the laser beam, $\mathbf{x}$, $\mathbf{x}_{laser}$ and $r_{laser}$ [m] represent the position of selected point on the substrate top surface, the position of the laser beam axis and the nominal laser radius, respectively.

The heat loss model includes a convection term at the lateral and bottom boundaries of the substrate with parameters specified in Section 3. Preliminary laser scanning studies indicate that radiation heat loss has a negligible impact on temperature variation [44]. Therefore, the loss term $\dot{Q}_{loss}$ [W m⁻²] writes:

$$ \dot{Q}_{loss} = h_c(T_{ext} - T) \quad (3) $$

where $h_c$ [W m⁻² K⁻¹] is the convection coefficient and $T_{ext}$ [K] is the external temperature.

Fig. 1 describes the partial equilibrium (PE) transformation path [46], for the alloy 316L, computed with the TCFE9 thermodynamic database [47], with activation of the peritectic transformation including liquid (l), ferrite (a), and austenite (γ). Compared to the classical Gulliver-Scheil (GS) model, the chemical potential of the interstitial elements (e.g., C) are equal in all phases, therefore assuming high diffusivity. The other assumptions of the GS are maintained, i.e. no diffusion of substitutional elements in the solid phases and complete mixing of the liquid. The detailed algorithm to compute the PE solidification path with peritectic reaction is given in Ref. [46]. Fig. 1 shows the evolution with temperature of the volume fraction of liquid, $g^l$, and the total volume fraction of the solid mixture, $g^s$, with $g^s + g^l = 1$. The solid mixture is a combination of α-ferrite and γ-austenite, $g^s = g^a + g^\gamma$, with primary solidification of α-ferrite from the liquidus temperature, $T_{liquidus}$, and formation of γ-austenite below the peritectic temperature, $T_{peritectic}$. The average volumetric enthalpy of the solid mixture writes $g^s(\rho h)^s = g^a(\rho h)^a + g^\gamma(\rho h)^\gamma$, hence the expression of the global average volumetric enthalpy:

$$ \langle \rho h \rangle = g^l \langle \rho h \rangle^l + g^a \langle \rho h \rangle^a + g^\gamma \langle \rho h \rangle^\gamma . \quad (4) $$

This approach was validated by Xue et al. [48] for welding conditions. The use of the PE approximation with activation of the peritectic reaction is believed to be superior to the more classical lever rule or GS approximations as it predicts the predominance of austenite after
<figure><img src="image_3.png" /><figcaption>Fig. 1. Temperature evolution of (a) the volumetric enthalpy of the liquid and solid phases involved during solidification and (b) their volume fraction assuming partial equilibrium [46] for the alloy 316L, based on thermodynamic data provided by the TCFE9 thermodynamic database [47]. The solid phase is made of a mixture of α ferrite and γ austenite depending on the temperature [48].</figcaption></figure>

completion of solidification (below the solidus temperature, 1574 K), in accordance with experimental observations [49]. Ferrite only forms in a narrow temperature interval between the liquidus of the alloy, 1725 K, and the peritectic temperature, 1699 K. In the following solid mechanical analysis, the focus of the study is on austenite with a FCC structure.

## 2.2. Cellular automaton method

To model the grain structure, a CA grid is established. It consists of a structured array of cubic cells uniformly distributed over a CA mesh. Topological links are defined between grid and mesh. Both the CA mesh and the CA grid are illustrated in Figs. 2(c) and 2(d), respectively. Given the comprehensive documentation of the CA methodology in existing literature, including its use for welding and L-PBF [28,29,48, 50], the present discussion serves only to reiterate its basic principles. The CA grid tracks mesoscopic grain evolution. This is related to the temperature field computed on a FE mesh (Fig. 2(a)). Similarly, the grain structure is used by a crystal plasticity model solved a CP mesh Fig. 2(b). The CA mesh serves as an auxiliary mesh to distribute the information from the CA grid to the FE/CP meshes, and vice-versa, as schematized in Fig. 2. It avoids creating two topological links with the CA grid: one with the FE mesh plus one with the CP mesh.

Each CA cell is assigned a state index indicating its current state, categorized as “liquid”, “mushy” or “solid”. In the context of melting, i.e., above the liquidus temperature, the cell index is “liquid”. A cell is “mushy” when its temperature falls below the liquidus temperature of the alloy and is surrounded by at least one neighboring liquid cell. At this point, it begins to grow, endowed with a small octahedral structure. In the context of solidification, a cell could remain liquid below the liquidus temperature and may be captured by a growing “mushy” cell based on a dedicated algorithm [51].

The CA method assumes that primary solidification takes place with a dendritic microstructure along the main (100) crystallographic directions. Mushy cells are initialized with a small octahedron centered within the cell. For clarity, in the simplified 2D schematic diagram of Fig. 2(e), it is represented as a quadrilateral with (10) half-diagonals perpendicular to each other. The diagonals of the octahedron are aligned with the (100) crystallographic directions of the grain. Dendrite tip growth kinetics is used to update the half-diagonal [30–32]. The dendrite tip velocity $v_{(100)}$ [mm s⁻¹] is approximated by a power law [48]:

$$v_{(100)} = A \cdot (T_{liquidus} - T_{cell,(100)})^{2.7} \quad (5)$$

where $A = 1.8 \cdot 10^{-6}$ [mm s⁻¹ K⁻².⁷] is the kinetic coefficient and $T_{cell,(100)}$ [K] is the temperature at the dendrite tip, i.e. at the (100) vertices of the octahedron associated to the CA cell. Upon growth of the octahedron over the center of an adjacent liquid cell, the latter undergoes a transition to a “mushy” state, inheriting the same crystallographic orientation as the capturing cell. If there are no liquid cells in the neighborhood of a mushy cell, its growth ceases, leading to a transition of the cell to a “solid” state. The methodology to initialize the octahedron of the newly captured cell and more details of the mechanism are described elsewhere [51].

A grain is represented by a cluster of CA cells that share a unique crystallographic orientation defined by a triplet of Euler angles $(\phi_1, \phi, \phi_2)$ in the Bunge convention [52]. The coupling scheme is used in the CAFE method as shown in Fig. 2. The temperature field computed at nodes n of the FE mesh, $T_{n,FE}$, is projected to the nodes of the CA mesh, $T_{n,CA}$, and then this nodal field is projected to the cells v by linear interpolation, $T_{v,cell}$. For the triplet of Euler angles, the value of the CA mesh element, $(\phi_1, \phi, \phi_2)_{v,CA}$, is chosen as the set of Euler angles of the nearest CA cell to the centroid of the element, $(\phi_1, \phi, \phi_2)_{v,cell}$. Then the triplet values are interpolated from the center of the CA element to the center of the CP element, $(\phi_1, \phi, \phi_2)_{v,CP}$.

## 2.3. Mechanical modeling

The objective of this Section is to construct a mechanical model, capable of capturing the thermo-mechanical behavior after solidification induced by laser scanning.

According to Dantzig and Rappaz [53], the mechanical behavior of mushy metal shows significant differences above and below a critical point called the coherency point, which is defined by a critical liquid fraction $g'_{coh}$. Below this threshold (where $g' < g'_{coh}$), the mushy metal can be considered homogeneous but is inherently weakened by the presence of isolated liquid pockets and localized interdendritic films. This point can be defined as the liquid fraction below which the semi-solid medium can transmit stresses, due to the appearance of a continuity of the solid skeleton. Specifically, this means that elastic effects begin to be perceptible in the mushy zone at high solid fraction. With reference to the experimental studies of Vicente-Hernandez [54] using needle indentation tests of mushy steels, the selected value for $g'_{coh}$ varies between 0.3 and 0.4 in the different grades studied. In the present study, the selected value for $g'_{coh}$ is 0.3 as shown in Fig. 1(b), corresponding to a critical temperature $T_{coh} = 1690$ K. To account for the different mechanical response above and below this critical liquid fraction, two constitutive equations are applied:
<figure><img src="image_4.png" /><figcaption>Fig. 2. Schematics of the coupling strategies developed for the Cellular Automaton (CA) - Finite Element (FE) - Crystal Plasticity (CP) model.</figcaption></figure>

- A solid-like thermo-elasto-viscoplastic (TEVP) constitutive model when the liquid fraction is lower than the critical value $g'_{coh}$.

- A fluid-like thermo-viscoplastic (TVP) model when $g'$ is higher than $g'_{coh}$.

This separate treatment approach is derived from previous studies to analyze the solidification process [55,56] depending on the local state of the metal.

The TVP approach is described first in Section 2.3.1. Then, the constitutive law of the TEVP model, which is enriched with crystal plasticity, is described in Section 2.3.2. Finally, the non-linear global and local resolutions are explained in Section 2.3.3.

### 2.3.1. Thermo-viscoplastic model

As for the TVP model in the mushy state, the metallic alloy exhibits a characteristic two-phase composition, comprising both liquid and solid phases. In the present approach, the mushy zone is approximated as a homogenized “single phase” continuum with a generalized non-Newtonian fluid behavior:

$$ \dot{\epsilon} = \dot{\epsilon}^{vp} + \dot{\epsilon}^{th} \qquad (6) $$

$$ \dot{\epsilon}^{vp} = \frac{1}{2K_{vp}} (\sqrt{3\bar{\epsilon}})^{1-m_{vp}} \mathbf{s} \qquad (7) $$

$$ \dot{\epsilon}^{th} = -\frac{1}{3\rho} \frac{\partial \rho}{\partial T} \frac{dT}{dt} \mathbf{I} = \alpha \bar{T} \mathbf{I} \qquad (8) $$

where $\dot{\epsilon}$ is the strain rate tensor, symmetric part of the velocity gradient: $\dot{\epsilon} = 1/2(\nabla v + (\nabla v)^T)$ decomposed into a viscoplastic part and a thermal part. The former is proportional to the stress deviator $\mathbf{s} = \sigma + \rho\mathbf{I}$ (with pressure $p = -\frac{1}{3}\text{tr}(\sigma)$), through a power law where $K_{vp}$ is the viscoplastic consistency, and $m_{vp}$ is the strain rate sensitivity. The later part, $\dot{\epsilon}^{th}$, includes dilatation and solidification shrinkage where $\alpha$ denotes the linear thermal expansion coefficient and $\mathbf{I}$ is the identity tensor.

In the present study, our primary focus is on the solid region below the coherency point. The viscoplastic mechanical behavior occurring above this point falls outside the primary scope of our investigation. Therefore, the Newtonian behavior is adopted by taking $m_{vp} = 1$ and $K_{vp}$ is the dynamic viscosity of the homogenized semi-solid continuum, with a higher value than the liquid viscosity for reasons of numerical stability.

### 2.3.2. Thermo-elasto-viscoplastic model

As for the TEVP model, crystal elasto-viscoplasticity (CEVP) is adopted. So, the solution is conducted in the crystal reference frame $(x, y, z)|_C$ aligned with the $\langle 100 \rangle$ directions of the 316L cubic crystal lattice, through the rotation matrix defined by the three Euler angles $(\phi_1, \phi, \phi_2)$ of Bunge's convention (see Box I):

Tensors and vectors in the different reference frames write $\mathbf{d} = \mathbf{M}\mathbf{d}_C$ and $\mathbf{A} = \mathbf{M}\mathbf{A}_C\mathbf{M}^T$ where $\mathbf{d}$ and $\mathbf{d}_C$ denote the column-vectors respectively in the global reference frame $(x, y, z)$ and in the crystal reference frame $(x, y, z)|_C$, and similarly for tensors $\mathbf{A}$. The expression of strain rate and stress tensor used in CEVP hereafter are in the crystal reference, so the subscript $C$ is not applied for the sake of readability.

Owing to the small strains characterizing L-PBF, the additive decomposition of the strain rate tensor into elastic, viscoplastic, and thermal contributions is considered:

$$ \dot{\epsilon} = \dot{\epsilon}^{el} + \dot{\epsilon}^{vp} + \dot{\epsilon}^{th} \qquad (10) $$

$$ \dot{\epsilon}^{el} = \mathbb{C}^{-1} \dot{\sigma} \qquad (11) $$

where $\mathbb{C}$ is the rank four anisotropic elastic stiffness tensor with $\mathbb{S} = \mathbb{C}^{-1}$ being the elastic compliance tensor. The deviatoric part of the preceding expressions write:

$$ \dot{\epsilon} = \dot{\epsilon}^{el} + \dot{\epsilon}^{vp} \qquad (12) $$

$$ \dot{\epsilon}^{el} = \text{dev}(\mathbb{C}^{-1}\dot{\sigma}) \qquad (13) $$

where $\dot{\epsilon}$ denotes the deviatoric part of the strain rate tensor: $\dot{\epsilon} = \dot{\epsilon} - (1/3)\text{tr}(\dot{\epsilon})\mathbf{I}$, and where $\text{dev}(\mathbf{A})$ denotes the deviatoric part of a tensor $\mathbf{A}$. Note that $\dot{\epsilon}^{vp} \equiv \dot{\epsilon}^{th}$.

Denoting $\tilde{\mathbf{A}}$ as the Voigt's representation of tensors and in the context of cubic anisotropy, we have:

$$ \tilde{\mathbf{C}} = \begin{bmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \end{bmatrix} \qquad (14) $$

Besides, we have:

$$ \text{dev}(\tilde{\mathbf{A}}) = \tilde{\mathbb{D}}\tilde{\mathbf{A}} \quad \text{with} \quad \tilde{\mathbb{D}} = \frac{1}{3} \begin{bmatrix} 2 & -1 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ -1 & -1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3 \end{bmatrix} \qquad (15) $$
$$ \mathbf{M} = \begin{bmatrix} \cos \varphi_1 \cos \varphi_2 - \sin \varphi_1 \cos \phi \sin \varphi_2 & -\cos \varphi_1 \sin \varphi_2 - \sin \varphi_1 \cos \phi \cos \varphi_2 & \sin \varphi_1 \sin \phi \\ \sin \varphi_1 \cos \varphi_2 + \cos \varphi_1 \cos \phi \sin \varphi_2 & -\sin \varphi_1 \sin \varphi_2 + \cos \varphi_1 \cos \phi \cos \varphi_2 & -\cos \varphi_1 \sin \phi \\ \sin \varphi_2 \sin \phi & \cos \varphi_2 \sin \phi & \cos \phi \end{bmatrix} \quad (9) $$

Box I.

Note that $\tilde{\mathbf{D}}$ is a singular matrix, not invertible. However, it is easy to show that $\tilde{\mathbf{D}}\tilde{\mathbf{C}}^{-1} = \tilde{\mathbf{C}}^{-1}\tilde{\mathbf{D}}$. Therefore, Eq. (13) writes:

$$ \tilde{\mathbf{e}}^l = \tilde{\mathbf{D}}\tilde{\mathbf{C}}^{-1}\tilde{\mathbf{\sigma}} = \tilde{\mathbf{C}}^{-1}\tilde{\mathbf{D}}\tilde{\mathbf{\sigma}} = \tilde{\mathbf{C}}^{-1}\tilde{\mathbf{s}} \quad (16) $$

The viscoplastic part of the strain rate tensor is expected by use of crystal viscoplasticity, in the context of FCC lattice. The relationship between $\dot{\epsilon}^{vp}$ and the shear rates of the different slip systems is classically expressed by [57]:

$$ \dot{\epsilon}^{vp} = \sum_{\alpha=1}^{N} \dot{\gamma}^{\alpha} \mathbf{P}^{\alpha} \quad (17) $$

where $N$ is the number of slip systems and $\dot{\gamma}^a$ is the shear rate of slip system $a$. $\mathbf{P}^a$ is the symmetrized Schmid tensor, expressed as:

$$ \mathbf{P}^a = \frac{1}{2} (\mathbf{m}^a \otimes \mathbf{n}^a + \mathbf{n}^a \otimes \mathbf{m}^a) \quad (18) $$

where vectors $\mathbf{m}^a$ and $\mathbf{n}^a$ are unit vectors which describe the slip direction within the slip plane and the normal to the slip plane, respectively.

The plastic shear on a slip system $\alpha$ is expressed by a Norton-type power-law equation with threshold initially proposed by Cailletaud [58], also used in [59,60]:

$$ \dot{\gamma}^{\alpha} = \dot{\gamma}_0 \left[ \frac{|\tau^{\alpha}| - \tau_c^{\alpha}(T)}{K(T)} \right]_+^{\frac{1}{m(T)}} \operatorname{sgn}(\tau^{\alpha}) = \dot{\gamma}_0 \left[ \frac{\|T_{RSS}^{\alpha}\| - \tau_c^{\alpha}(T)}{K(T)} \right]_+^{\frac{1}{m(T)}} \frac{T_{RSS}^{\alpha} \cdot \mathbf{m}^{\alpha}}{\|T_{RSS}^{\alpha}\|} \quad (19) $$

where $[x]_+$ is an alternative Macaulay notation to take the positive part of the scalar $x$, $|x|$ is the absolute value of the scalar $x$ and $\|x\|$ is the Euclidean norm of the vector $x$. $\tau^a$ is the resolved shear stress and $T_{RSS}$ denotes the resolved shear surface force:

$$ T_{RSS}^{\alpha} = \tau^{\alpha} \mathbf{m}^{\alpha} = (\mathbf{P}^{\alpha} : \boldsymbol{\sigma}) \mathbf{m}^{\alpha} = (\mathbf{P}^{\alpha} : \mathbf{s}) \mathbf{m}^{\alpha} \quad (20) $$

$\tau_c^a$ is the critical resolved shear stress for the slip system $\alpha$. $K$ is the viscous parameter for slip. $\dot{\gamma}_0$ and $m$ are material constants, which indicate the reference strain rate and rate sensitivity of slip, respectively.

Here, an incremental hardening model [57,61,62] is introduced:

$$ \tau_c^a(T) = \tau_0(T) + \int_0^t \sum_{\beta=1}^N h_{a\beta} |\dot{\gamma}^\beta| dt \quad (21) $$

$$ h_{a\beta} = q_{a\beta} \left[ h_{final} \left( 1 + (h_{ratio} - 1) e^{-h_{exp}} |\dot{\gamma}^\beta| \right) \right] \quad (22) $$

where $q_{a\beta}$ is a matrix of coefficients for latent hardening and $\dot{\gamma}^a = \int_0^t \dot{\gamma}^a dt$ the accumulated shear for slip system $\alpha$. $h_{final}$, $h_{ratio}$ and $h_{exp}$ are coefficients that determine the strength of the exponential part.

All the parameters relevant to the simulations are detailed in Table 1.

### 2.3.3. Numerical implementation

The momentum conservation should be satisfied everywhere in the domain $\Omega$ and at every instant as shown in the first equation of Eq. (23). The mechanical resolution is in the velocity/pressure framework, where the pressure $p$ is retained as a primary variable. The tensor $s$ is associated with the velocity field $v$ and only involves the deviatoric component of the constitutive equation. Therefore, an additional equation should be considered that includes the spherical component of the

constitutive equation, as shown in the second equation of Eq. (23).

$$ \left\{ \begin{array}{l} \nabla \cdot \mathbf{s} - \nabla p + \rho g = 0, \\ \nabla \cdot \mathbf{v} = -H(g_{coh}^l - g^l) \frac{\dot{p}}{\chi} - \frac{1}{\rho} \frac{\partial \rho}{\partial T} \frac{dT}{dt}, \end{array} \right. \quad (23) $$

where $\chi = (C_{11} + 2C_{12})/3$ is the elastic compressibility modulus. Note that it can be shown that the second relationship in Eq. (23) is mass conservation, so that Eq. (23) concurrently encompasses momentum and mass conservation.

The standard Heaviside function $H$ is introduced as an indicator of the utilization of the CEVP constitutive model.

$$ H(g_{coh}^l - g^l) = \begin{cases} 1 & \text{if } g^l < g_{coh}^l \\ 0 & \text{if } g^l \ge g_{coh}^l \end{cases} \quad (24) $$

To implement Eq. (23), the liquid fraction $g^l(T)$ and the density $\rho(T)$ are determined from the thermodynamic data base by pre-computed temperature-dependent relations (see Fig. 1(b) for the former; $\rho(T)$ in Eq. (34)). These material properties are then accessed via look-up tables that store their one-to-one correspondence with temperature.

The weak form of Eq. (23) is solved using the finite element method with a mixed velocity/pressure formulation, using $P1+/P1$ tetrahedral elements. These elements, initially proposed by Arnold et al. [63] and Fortin and Fortin [64], have an enriched velocity discretization ($P1+$, the enrichment consisting of a “bubble”-type function) and a linear ($P1$) pressure interpolation, so that they comply with the Brezzi-Babuska stability conditions. Because the additional degrees of freedom are internal to each element, they can be condensed after the assembly process and prior to resolution to keep only 4 unknowns per node: 3 velocity components and 1 pressure value. This formulation was previously developed for conventional elasto-viscoplastic behavior laws and is fully documented in Refs. [65,66].

In this context, the stress deviator $s$ must be expressed as a function of the total strain rate tensor $\dot{\epsilon}$, which is linked to the velocity field $v$. In addition, the Newton-Raphson algorithm is used to solve the discretized weak form of Eq. (23). This requires calculating the tangent stiffness, involving the tangent modulus $\partial s/\partial \dot{\epsilon}$ in each finite element in CEVP state (that is for which, at element center, $g^l < g_{coh}^l$). These two issues (establishing the expression of $s(\dot{\epsilon})$ and $\partial s/\partial \dot{\epsilon}$) are examined below for such elements. For elements following a Newtonian law (that is for which, at element center, $g^l \ge g_{coh}^l$) see [65,66].

In each finite element governed by CEVP behavior, the three Euler angles defining crystal orientation are determined using data from the CA method. For substrate elements, these angles are initially set based on existing texture. Upon solidification, crystal orientation is assigned according to CA predictions for epitaxial growth, whereas elements that melt or remelt have their Euler angles reset to zero as soon as $g^l$ exceeds $g_{coh}^l$. However these Euler angles are ignored since elements such as $g^l > g_{coh}^l$ are treated with the viscoplastic TVP constitutive model. Thus, the Euler angles are consistently updated at each simulation time step. Note, however, that the lattice rotation resulting from the combined shears in the different slip systems is not taken into account: consistent with the small-deformation framework, the CEVP resolution does not induce any change in the Euler angles.

Since calculating the global residual vector requires $s$ in the global frame, while CEVP is solved in the crystal reference frame $(x, y, z)|_C$, the sequence is as follows:
* The local strain rate tensor $\dot{\epsilon}$ from the velocity field $v$ and the previous stress deviator $s^i$ are transformed by M-transformation (Eq. (9)) to $\dot{\epsilon}_{|C}$ and $s^i_{|C}$ in the crystal reference frame.

* The CEVP law is applied in the crystal reference frame $(x, y, z)|_C$ to compute the new stress deviator estimate $s|_C$ and the local tangent modulus $(\partial s/\partial \dot{\epsilon})|_C$.

* $s|_C$ and $(\partial s/\partial \dot{\epsilon})|_C$ are transformed back to the global frame for use in the global residual vector and stiffness matrix.

In the next Sections, the focus is made on the second step of this sequence.

### 2.3.3.1 Local Newton-Raphson algorithm

The objective of the local resolution is to determine, in the crystal reference frame of each element, the deviatoric stress tensor $s$ which satisfies Eqs. (12),(13) and (17), in which $\dot{\epsilon}$ represents the strain rate tensor corresponding to the current estimate of the velocity field $v$. From these three equations, the stress tensor $s$ should satisfy the following equation, where Voigt's notation is used to get a vector form:

$$ \tilde{\mathbf{R}}(\tilde{s}) = \tilde{\mathbf{e}} - \tilde{\mathbf{S}} \frac{(\tilde{s} - \tilde{s}^{\text{r}})}{\Delta t} - \sum_{\alpha=1}^{N} \dot{\gamma}^{\alpha}(\tilde{s}) \tilde{\mathbf{P}}^{\alpha} = 0 \quad (25) $$

where $\mathbf{R}$ is the residual. It is worth noting that in most CP numerical algorithms reported in the literature, the plastic shear rates in each slip system $\dot{\gamma}^a$ (or their corresponding increments $\Delta\dot{\gamma}^a$) are treated as primary unknowns [67–69]. Note that in the present work, the number of primary unknowns is minimized, since they are limited to the components of the deviatoric stress tensor $s$. Considering the zero-trace property of $s$, there are only 5 independent primary unknowns.

This vector equation can be solved by a Newton-Raphson iterative algorithm. At each iteration $v$, the updated estimation of $s$ is calculated by:

$$ \left(\frac{\partial \tilde{\mathbf{R}}}{\partial \tilde{s}}\right) \partial \tilde{s}^{(v)} = -\tilde{\mathbf{R}}(\tilde{s}^{(v-1)}) \quad \text{and} \quad \tilde{s}^{(v)} = \tilde{s}^{(v-1)} + \delta \tilde{s} \quad (26) $$

The local tangent matrix $\partial \tilde{\mathbf{R}}/\partial \tilde{s}$ is straightforwardly expressed from Eq. (25). By components, we have:

$$ \frac{\partial \tilde{R}_i}{\partial \tilde{s}_j} = -\frac{\tilde{S}_{ij}}{\Delta t} - \sum_{a=1}^{N} \frac{\partial \dot{\gamma}^a}{\partial \tilde{s}_j} \tilde{P}_i^a \quad (27) $$

where $\tilde{R}_i$, $\tilde{s}_j$, $\tilde{S}_{ij}$ and $\tilde{P}_i^a$ are the Voigt's components of $\mathbf{R}$, $\mathbf{s}$, $\mathcal{S}$ and $\mathbf{P}^a$. The convergence criterion is defined as $\max_j (\delta \tilde{s}_j^{(v)}/\tilde{s}_j^{(v)}) < 10^{-2}$. As shown in Eq. (27), the completion of $\partial \mathbf{R}/\partial s$ requires the calculation of the compliance tensor $\partial \dot{\gamma}^a/\partial s$, detailed hereafter.

### 2.3.3.2. Calculation of the CEVP compliance tensor $\partial \dot{\gamma}^a/\partial \tilde{s}$

In the definition of $\dot{\gamma}^a$ (Eq. (19)), the Macauley bracket is approached by a derivable expression, inspired by the work of Zhang and Oskay [70], but expressed in a different way on a neighborhood $[-c, +c]$ around zero value:

$$ [x]_+ = \begin{cases} 0 & \text{if } x < -c \\ \frac{x}{2} + \frac{x^2}{4c} + \frac{c}{4} & \text{if } |x| \le c \\ x & \text{if } x > c \end{cases} \quad (28) $$

Note that $c$ has the dimension of a stress, and that it should be sufficiently lower than $\tau_c^a$. This is why $c$ is written as $c = \beta \tau_c^a$, with $\beta$ a small positive constant. In addition, the values of the different $\tau_c^a$ are taken at the beginning of the time increment, so that during the local Newton-Raphson iterations, their derivative value is constantly

null. After calculation (see Appendix A), we have:

$$ \frac{\partial \dot{\gamma}^a}{\partial \tilde{s}_j} = \begin{cases} 0, & \text{if } \|T_{RSS}^a\| - \tau_c^a < -\beta \tau_c^a \\ \frac{\dot{\gamma}_0}{2\beta \tau_c^a m K^{1/2m}} [\|T_{RSS}^a\| - \tau_c^a]_+^{\frac{1}{m}-1} \times (\|T_{RSS}^a\| - \tau_c^a(1-\beta)) \tilde{P}_j^a, & \text{if } -\beta \tau_c^a \le \|T_{RSS}^a\| - \tau_c^a \le \beta \tau_c^a, \\ \frac{\dot{\gamma}_0}{m K^{1/m}} [\|T_{RSS}^a\| - \tau_c^a]_+^{\frac{1}{m}-1} \tilde{P}_j^a, & \text{if } \|T_{RSS}^a\| - \tau_c^a > \beta \tau_c^a, \end{cases} \quad (29) $$

### 2.3.3.3. Expression of the local tangent modulus $\partial s/\partial \dot{\epsilon}$

The solution of Eq. (25) also provides, as a by-product, the tangent modulus $\partial s/\partial \dot{\epsilon}$ used in the expression of the global mechanical tangent stiffness matrix.

First, let us calculate $\partial \dot{\epsilon}/\partial s$. From Eq. (12), we get:

$$ \frac{\partial \tilde{\mathbf{e}}_i}{\partial \tilde{s}_j} = \frac{\partial \tilde{\mathbf{e}}_i^l}{\partial \tilde{s}_j} + \frac{\partial \tilde{\mathbf{e}}_i^{lp}}{\partial \tilde{s}_j} \quad (30) $$

Based on Eqs. (11), (17) and (29), this leads (in the case $\|T_{RSS}^a\| - \tau_c^a > \beta \tau_c^a$) to:

$$ \frac{\partial \tilde{\mathbf{e}}_i}{\partial \tilde{s}_j} = \frac{1}{\Delta t} \tilde{S}_{ij} + \sum_{\alpha=1}^{N} \frac{\dot{\gamma}_0}{m K^{1/m}} [\|T_{RSS}^a\| - \tau_c^a]_+^{\frac{1}{m}-1} \tilde{P}_i^a \tilde{P}_j^a \quad (31) $$

Therefore:

$$ \tilde{d\mathbf{s}} = \left[ \frac{1}{\Delta t} \tilde{\mathbf{S}} + \tilde{\mathbf{G}} \right]^{-1} \tilde{d\dot{\mathbf{e}}} = \left[ \frac{1}{\Delta t} \tilde{\mathbf{S}} + \tilde{\mathbf{G}} \right]^{-1} \tilde{D} \tilde{d\dot{\mathbf{e}}} \quad (32) $$

Finally, the tangent modulus $\partial s/\partial \dot{\epsilon}$ can be written as:

$$ \frac{\partial \tilde{s}}{\partial \tilde{\dot{\epsilon}}} = \left[ \frac{1}{\Delta t} \tilde{\mathbf{S}} + \tilde{\mathbf{G}} \right]^{-1} \tilde{D} \quad (33) $$

In the local resolution, $s$ and the incremented variables $(\dot{\gamma}^a$ and $\tau_c^a)$ are conserved and carried forward to the next time step.

### 2.3.3.4. Numerical treatment of state variables in phase transition

During the phase transition from liquid to solid (numerically triggered when the liquid volume fraction falls below the critical threshold $g_{coh}^l$), the critical resolved shear stresses $\tau_c^a$ are reinitialized. As shown in Eq. (21), each $\tau_c^a$ comprises two components: a temperature-dependent basal term $\tau_0(T)$ and a contribution from the accumulated slip in crystal slip systems. Here, the shear slip accumulations are reset to zero for the newly solidified grains. Thus, the reinitialized $\tau_c^a$ depend solely on temperature (i.e., $\tau_0(T)$). Consistently, the transient plastic strain rate in liquid region is not accumulated as permanent deformation. This treatment guarantees thermo-mechanics consistency across moving phase boundaries, while avoiding numerical instability.

## 3. Calibration, first validation and simulation setup

### 3.1. Material parameters

The thermal parameters of 316L stainless steel are deduced from [47, 73]:

$$ \begin{aligned} \lambda &= 0.0143T + 10.076 \, [\text{W K}^{-1} \text{m}^{-1}] \\ \alpha &= 5.655 \times 10^{-9}T + 1.556 \times 10^{-5} \, [\text{K}^{-1}] \\ \rho &= 7928 - 0.6087T + 5.989 \times 10^{-4}T^2 \, [\text{kg m}^{-3}] \end{aligned} \quad (34) $$

Regarding mechanical calculations, the data are given in Table 1. The components of the elastic stiffness tensor are assumed to vary linearly with temperature [71], and this linear relationship is extrapolated to higher temperatures until $g^l = g_{coh}^l$:

$$ C_{ij}(T) = C_{ij}(T_0) + \frac{dC_{ij}}{dT}(T - T_0) \quad (35) $$
<table><caption>**Table 1**

Mechanical parameters for 316L stainless steel.</caption><tr><td>Elastic constant C<sub>11</sub>(T<sub>0</sub>)</td><td>[71]</td><td>204.6 ×10<sup>3</sup> MPa</td></tr><tr><td>Elastic constant C<sub>12</sub>(T<sub>0</sub>)</td><td>[71]</td><td>137.7 ×10<sup>3</sup> MPa</td></tr><tr><td>Elastic constant C<sub>44</sub>(T<sub>0</sub>)</td><td>[71]</td><td>126.2 ×10<sup>3</sup> MPa</td></tr><tr><td>Elastic constant derivative dC<sub>11</sub>/dT</td><td>[71]</td><td>-90.33 MPa K<sup>-1</sup></td></tr><tr><td>Elastic constant derivative dC<sub>12</sub>/dT</td><td>[71]</td><td>-45.10 MPa K<sup>-1</sup></td></tr><tr><td>Elastic constant derivative dC<sub>44</sub>/dT</td><td>[71]</td><td>-51.78 MPa K<sup>-1</sup></td></tr><tr><td>Reference temperature T<sub>0</sub></td><td></td><td>273.15 K</td></tr><tr><td>Boltzmann constant k<sub>B</sub></td><td></td><td>1.381 × 10<sup>-23</sup> J K<sup>-1</sup></td></tr><tr><td>Burgers vector magnitude b<sub>B</sub></td><td></td><td>2.48 × 10<sup>-10</sup> m</td></tr><tr><td>Reference plastic strain rate γ'<sub>0</sub></td><td></td><td>1.0 × 10<sup>-7</sup> s<sup>-1</sup></td></tr><tr><td>Critical resolved shear stress τ<sub>0</sub>(T<sub>0</sub>) (Eq. (37))</td><td>[71]</td><td>68.5 MPa</td></tr><tr><td>Viscous parameter for slip K(T<sub>0</sub>) (Eq. (37))</td><td></td><td>17.8 MPa</td></tr><tr><td>Hardening parameter h<sub>final</sub>(T<sub>0</sub>) (Eq. (37))</td><td></td><td>1.45 × 10<sup>5</sup></td></tr><tr><td>Hardening parameter h<sub>exp</sub></td><td></td><td>1200</td></tr><tr><td>Hardening parameter h<sub>ratio</sub></td><td></td><td>35</td></tr><tr><td>Hardening coefficient matrix q<sub>aa</sub></td><td>[72]</td><td>1.0</td></tr><tr><td>Hardening coefficient matrix q<sub>aβ</sub></td><td>[72]</td><td>1.4</td></tr></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Shift functions $g_{\tau_0}(T)$ [78], $g_h(T)$ and $g_K(T)$ associated with the temperature variations of the viscoplastic characteristics.</figcaption></figure>

For the viscoplastic part, a formulation referring to [44,74,75], is used to express the temperature-dependent sensitivity:

$$m = \frac{9k_B T}{b_B^3 \sqrt{2C_{44} \times (C_{11} - C_{12})}} \quad (36)$$

where $k_B$ is the Boltzmann constant and $b_B$ is the Burgers vector magnitude.

In addition, the yield parameter $\tau_0$, the hardening parameter $h_{final}$ and the viscous parameter $K$ are temperature dependent by multiplying their value at reference temperature $T_0$ by shift functions $g(T)$:

$$\begin{align}
\tau_0(T) &= \tau_0(T_0)g_{\tau_0}(T) \nonumber \\
h_{final}(T) &= h_{final}(T_0)g_h(T) \tag{37} \\
K(T) &= K(T_0)g_K(T) \nonumber
\end{align}$$

Experimental data from [71,73,76] were used to determine these parameters. Depradeux [73] performed uniaxial tensile tests on solution annealed 316L stainless steel, which exhibits relatively low yield strength [77], over a wide temperature range (293 K to 1273 K) at two strain rates ($2.5 \times 10^{-4} \text{ s}^{-1}$ and $2.5 \times 10^{-3} \text{ s}^{-1}$). At room temperature, Daymond and Bouchard [71] provided $\tau_0^a$ values, while Thiercelin et al. [78] later proposed shift functions based on Depradeux's data (Fig. 3).

In this study, $g_h(T)$, $g_K(T)$, $K(T_0)$ and $h_{final}(T_0)$ were calibrated. A mesh with randomly distributed Euler angles ($\phi_1, \phi, \phi_2$) was generated, with $\phi_1$ and $\phi_2$ varying randomly in $[0, \pi]$ and $\cos(\phi)$ in $[-1, 1]$. Simulated uniaxial tensile tests replicated the experimental conditions of Depradeux [73]. The resulting parameters $K(T_0)$ and $h_{final}(T_0)$ are listed in Table 1 and $g_h(T)$ and $g_K(T)$ are in Fig. 3. Figs. 4(a) and (b) show strong

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Comparison between the stress-strain computed responses and the experimental data across temperature and for two specified strain rates, namely $2.5 \times 10^{-4} \text{ s}^{-1}$ and $2.5 \times 10^{-3} \text{ s}^{-1}$. The experimental data [73] are denoted by dot markers, while the computed responses are illustrated by solid lines. Computations are based on Table 1, Eq. (37) and Fig. 3.</figcaption></figure>

agreement between calculated and experimental stress-strain responses across the temperature range and strain rates.

Depradeux's experiment [73] focused only on strains below 2%. To validate the calibrated parameters beyond this range, simulations were compared to uniaxial tensile test results of solution annealed 316L [76] covering strains up to 25% at room temperature (Fig. 5). The simulated results closely match the experimental data, confirming that the calibrated parameters accurately represent the properties, even outside the initial calibration range.

## 3.2. Thermal and CA model validation

The experimental results presented in Chadwick et al. [37] on single-track laser scanning of 316L stainless steel substrates are used to evaluate the capabilities of the present models for heat transfer and grain growth. The data has been made available upon request (see acknowledgments). The experiments had been performed using a unique coupling between continuous-wave laser and a scanning electron microscope (CW laser-SEM) [38]. For more details on the design of this device and laser specifications, please see the indicated reference.

The substrate used in the experiments was hot-rolled and annealed 316L stainless steel. The grain size was determined to be 19.3 µm by measuring the average equivalent sphere diameter of the grains. With this grain size, the corresponding grain density was measured to be $2.65 \times 10^5 \text{ mm}^{-3}$ and the initial texture of the substrate can be directly extracted from EBSD measurements. This information is used to generate
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Stress–strain response of an axial tensile test: comparison between experiment [76] and simulation using the previously identified parameters.</figcaption></figure>

the initial equiaxed grain structure [48]. For comparative analysis, the laser spot size is fixed at 60 µm. Fig. 6 presents the experimental results together with the corresponding simulated grain structures. Scans 1 and 2 were performed along the +x direction at a constant scanning speed of 500 mm s⁻¹ with laser power 66 W and 24 W, respectively. A fixed laser power of 24 W was maintained for Scans 2 and 3 while employing scanning speeds 500 mm s⁻¹ and 150 mm s⁻¹, respectively. In this simulation validation, all mesh settings are identical to those used in Section 3.3: see Table 3.

First, a comparison is made regarding the melt pool dimensions reported in Fig. 6. The melt pool dimensions in all simulations are deduced from the liquidus temperature iso-surface after the simulations reach steady-state conditions. The absorption coefficients of the material during simulation are deduced from Trapp et al. [79], with values 0.3, 0.28, and 0.27 for Scans 1, 2 and 3, respectively. In Scan 1, due to the highest laser power, both the width, *w*, and the depth, *d*, of the melt pool are the largest, whereas Scan 2 exhibits the smallest melt pool dimensions. As the scanning speed decreases, the melt pool size in Scan 3 increases compared to that in Scan 2. The simulation results reported in Fig. 6 demonstrate good agreement with measurements.

The simulated and experimental microstructures also show good agreement in overall morphology. In Scan 1, the largest melt pool size leads to the most extensive grain remelting and resolidification, while Scan 2 shows the least amount of grain remelting. Notably in Scan 2 and 3, as the melt pool becomes more circular in shape and the ratio *w*/*d* decreases, the direction of thermal gradient at the solidification front gradually changes from the ±*y* direction (at the maximum width of the melt pool) to the +x scanning direction. This occasionally results in the preferential survival of grains aligned with the scanning direction, which are ultimately but briefly retained along the centerline. It is also worth mentioning the very weak variation of texture associated to the present single-track laser experiments due to the low ratio between the melt pool dimensions and the initial size of the grain structure.

### 3.3. Calculation design

Following this first validation, the objective is to apply the simulation to successive melting and epitaxial grain growth using a laser heat source to scan a substrate, as well as studying the role of the scanning strategy. In this Section, the meshes, the boundary conditions and the scan paths are presented and discussed.

Simulation parameters are listed in Tables 2 and 3. The calculation is conducted on a substrate of dimensions 1.4 mm × 0.7 mm × 0.3 mm shown in Fig. 7. Laser melting passes are performed along the x-axis with different hatch spacings, as detailed in Table 2. Cases 1 and 2

<caption>**Table 2**</caption>

Definition of the laser scanning strategy for the different cases. All scans paths are parallel to the x-axis. For each pass, the y coordinate of scan line is given in the first column, while the scanning direction (SD)±x is given in a second column. Scans are between position x = 0.3 mm and x = 1.1 mm. All other parameters are maintained the same for all cases (see text).

<table><caption>**Table 3**</caption><thead><tr><th rowspan="2"></th><th colspan="2">Case 1</th><th colspan="2">Case 2</th><th colspan="2">Case 3</th><th colspan="2">Case 4</th></tr><tr><th>y [mm]</th><th>SD</th><th>y [mm]</th><th>SD</th><th>y [mm]</th><th>SD</th><th>y [mm]</th><th>SD</th></tr></thead><tbody><tr><td>Pass 1</td><td>0.23</td><td>+x</td><td>0.23</td><td>+x</td><td>0.23</td><td>+x</td><td>0.23</td><td>+x</td></tr><tr><td>Pass 2</td><td>0.28</td><td>+x</td><td>0.28</td><td>-x</td><td>0.30</td><td>+x</td><td>0.32</td><td>+x</td></tr><tr><td>Pass 3</td><td>0.33</td><td>+x</td><td>0.33</td><td>+x</td><td>0.37</td><td>+x</td><td>0.41</td><td>+x</td></tr><tr><td>Pass 4</td><td>0.38</td><td>+x</td><td>0.38</td><td>-x</td><td>0.44</td><td>+x</td><td>0.50</td><td>+x</td></tr><tr><td>Pass 5</td><td>0.43</td><td>+x</td><td>0.43</td><td>+x</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Process and numerical parameters for the multiple linear passes simulations.

<table><tr><td>Calculation domain Ω<sub>cal</sub></td><td>1.4 × 0.7 × 0.3 mm<sup>3</sup></td></tr><tr><td>CA cell size l<sub>cell</sub> (Fig. 2(e))</td><td>1.5 µm</td></tr><tr><td>CA element size l<sub>CA</sub> (Fig. 2(c))</td><td>5 µm</td></tr><tr><td>FE element size l<sub>FE</sub> (Fig. 2(a))</td><td>5 µm</td></tr><tr><td>CP element size at mesh refined zone l<sub>CP,refined</sub> (Fig. 2(b))</td><td>5 µm</td></tr><tr><td>CP element size at mesh coarse zone l<sub>CP,coarse</sub> (Fig. 2(b))</td><td>20 µm</td></tr><tr><td>Time step Δt</td><td>10<sup>-5</sup> s</td></tr><tr><td>Initial grain density N<sub>0</sub></td><td>9.77 × 10<sup>5</sup> mm<sup>-3</sup></td></tr><tr><td>Laser velocity v<sub>laser</sub></td><td>400 mm s<sup>-1</sup></td></tr><tr><td>Beam radius r<sub>laser</sub></td><td>0.04 mm</td></tr><tr><td>Laser power P<sub>laser</sub></td><td>100 W</td></tr><tr><td>Convective heat transfer coefficient h<sub>c</sub></td><td>100 W m<sup>-2</sup> K<sup>-1</sup></td></tr><tr><td>External temperature T<sub>ext</sub></td><td>500 K</td></tr><tr><td>Initial substrate temperature T<sub>init</sub></td><td>500 K</td></tr><tr><td>Absorption coefficient η</td><td>0.65</td></tr><tr><td>Solidus temperature T<sub>solidus</sub></td><td>1574 K</td></tr><tr><td>Liquidus temperature T<sub>liquidus</sub></td><td>1725 K</td></tr></table>

have a hatch spacing $h_p = 0.05$ mm. However, Case 1 scans solely in the +x direction, while Case 2 alternates between +x and -x. Cases 3 and 4 scan only in the +x direction, with hatch spacings increased to 0.07 mm and 0.09 mm, respectively. Each laser pass takes 4 ms: 1 ms of static preheating, 2 ms of movement from x = 0.3 mm to x = 1.1 mm, and 1 ms of cooling. The laser parameters are listed in Table 3: a Gaussian heat source with radius $r_{laser} = 0.04$ mm (Eq. (2)) moves at velocity $v_{laser} = 400$ mm s⁻¹ with power $P_{laser} = 100$ W and an absorption coefficient $\eta = 0.65$ [79].

The initial grain structure features randomly distributed Euler angles with a grain density of $9.77 \times 10^5$ mm⁻³. A typical Mackenzie distribution of the grain disorientation is found [80]. This was verified by comparing the actual distributions for the 4 cases (Table 2) with the exact distribution, as reported by the white distribution in Fig. B.18 of Appendix B. For the CP mesh, the calculation focuses on the melting and solidification zone, where a refined mesh is necessary to accurately capture grain boundaries. In this region (0.18 mm < x < 1.22 mm, 0.14 mm < y < 0.56 mm, -0.08 mm < z < 0 mm), the element size is 5 µm, smaller than the average grain size to capture boundary disorientation, as shown in Fig. 7. Outside this zone, the element size is coarser (20 µm) to reduce computation time. This coarser mesh in regions with random crystallographic orientations supports isotropic behavior. The total CP mesh comprises 1.7 million elements.

To verify the mesh independence, we conducted tests on three different CP mesh resolutions for the first pass, with mesh size of 3.5, 5, and 7 µm in the refined region. As shown in Fig. 8 the stress distributions for the 5 µm and 3.5 µm simulations are very similar. Therefore, for the sake of both computational accuracy and efficiency, a mesh size of 5 µm is adopted for the CP mesh in the following calculations.

As far as heat transfer along the surface is concerned, laser heating takes place along the top surface while convection applies to all other surfaces, as specified in Eq. (3). For mechanical boundary conditions,
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Comparative analysis of the top surface microstructure as obtained from the (left columns) CAFE simulation and (right column) experimental observation [37] for three different scan configurations. The melt pool size from simulation and experiment is given for each scan configuration.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Initial grain structure and laser scanning strategy (arrows) for Case 1 (see Table 2), together with close-up views with CP mesh. The black dashed box on the z = 0 mm top surface indicates the mesh refined region.</figcaption></figure>

the velocity field is fixed to 0 at the bottom surface, z = 0.3 mm, while stress-free conditions are applied to the other surfaces.

The calculations are performed on two AMD Rome 7502 processors, each comprising 32 cores. With the given number of elements and calculation platform, the total computation time is approximately 2 days per case to solve the heat flow, grain structure, and mechanical equilibrium.

# 4. Results

## 4.1. Grain growth

In this study it is assumed that there is no nucleation within the melt pool [37], as demonstrated by the experimental validation in Section 3.2 for laser scanning of a substrate. The grain structure is solely due to the partial remelting of the previously solidified substrate and

its resolidification. As the heat source passes through, the melt zone cools and solidifies. Epitaxial growth begins at the grain/pool interface, resulting in grain growth competition. As the growth front is forced to mainly develop along the temperature gradient, grain competition occurs in this direction, resulting in elongated grains.

The grain structures obtained by CAFE simulations for 4 cases studied are shown in Fig. 9 by means of Inverse Polar Figure (IPF) plots on CA cells. In Fig. 9(a), the solidus and liquidus temperature contours reveal the direction of the temperature gradient at the top surface (z = 0 mm). The semi-width of the melt pool is $r_{z=0,pool}$ = 0.08 mm, which is twice the radius of the laser. Grain structures for Case 1 and Case 3 (hatch spacing of 0.05 mm and 0.07 mm, respectively) are shown in Figs. 9(b) and (d). These structures consist of highly elongated grains extending across multiple passes, a result of epitaxial growth facilitated by remelting in successive passes. When the hatch distance is less than half the melt pool width, as in Cases 1 and 3, a significant portion of the previous pass is remelted, allowing the grains to grow uninterrupted along the y-axis. However, in Case 4, when the hatch distance increases by more than half the melt pool width ($h_p$ = 0.09 mm > $r_{z=0,pool}$), the grains located in the upper part of the previous pass, e.g. $y > y_{pass1}$) are not fully remelted. As a result, a discontinuity in grain morphology is observed as the selected (100) growth directions are not the same as in the previous pass. Thus, the grain structure appears “segmented” in Fig. 9(e). In Case 2, as the scanning direction alternates between the +x and -x directions, zigzag grains are observed in Fig. 9(c). This is the result of successive changes of the orientation of the temperature gradient between passes.

On the top surface of the 3D simulated domains, the melt zones are highlighted in black. The volumetric pole figures consider all grains encountered in the remelted and solidified volumes of (Cases 3 and 4) the 4 passes and (Cases 1 and 2) the 5 passes. As the grains preferentially develop along the (100) primary direction closest to the temperature gradient, the grains elongate predominantly in a direction nearly perpendicular to the heat source motion, stopping midway along the heat source trajectory, explaining the (100) fiber texture of Fig. 9(a) near the y-direction. However, the laser velocity has a significant impact on the grain morphology on the top surface [81]. Lower laser velocities shorten the drop-shaped melt pool, causing the temperature
<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of the computed $\sigma_{yy}$ profile along the line defined by $y = 0.23$ mm and $z = 0$ mm at the end of Pass 1 in Case 1 using different mesh resolutions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Grain structure for multiple linear passes using parallel scanning strategies along the x axis has been developed for different cases: (a) The temperature distribution, with the red and black contours representing the liquidus and solidus temperatures, respectively, and a melt pool width of 0.16 mm. Grain structures for (b) Case 1, (c) Case 2, (d) Case 3, and (e) Case 4. The pole figures correspond to the entire three-dimensional melt zone, with the volume texture index, $J$, explicitly annotated.</figcaption></figure>

gradient to deviate from the $y$ direction, which explains the observed grain shapes that are not well aligned with the $y$-direction. The (100) texture becomes weaker as the hatch spacing increases or the scanning direction changes, as shown by the $J$ texture index values in Fig. 9. Case 1 exhibits the strongest texture index of 1.2017 among the four cases, although it still represents a relatively weak texture. The texture indices for Cases 3 and 4 become progressively weaker with increasing hatch spacing. Notably, Case 2 demonstrates that the back-and-forth

scanning method leads to a more pronounced reduction of the texture index compared to Case 3.

## 4.2. Stress evolution

Figs. 10 and 11 present simulation results at three successive times during Pass 4 for Case 1: (a) the temperature field and melt pool contour, (b) the grain structure, and (c) and (d) stress distributions.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Evolution of (a) the temperature field, (b) the grain structure and the stress components (c) $\sigma_{yy}$ and (d) $\sigma_{xx}$ when the laser is at the beginning, midway and at the end of Pass 4 in Case 1. The black solid line represents the critical temperature at which the mechanical behavior changes, $T_{coh}$.</figcaption></figure>

Presentations are given along the top surface (Fig. 10) and along the longitudinal section containing the scan line (Fig. 11).

The temperature increases ahead of the melt pool due to the steep temperature gradient created by rapid heating of the laser beam and the relatively slow heat transfer by conduction. The thermal expansion of the heated solid is constrained by the colder, stiffer underlying material, creating compressive stress. As the laser moves along its path, rapid cooling induces shrinkage, which is constrained by the substrate, resulting in horizontal tensile stresses. These well-known general phenomena can be seen in Figs. 10(c-d) where the distribution of stress components $\sigma_{yy}$ and $\sigma_{xx}$ along the top surface is shown. The laser remelts grains possibly formed in the previous passes, modifying the tensile stresses caused by cooling-induced shrinkage in previously remelted zones. As the crystalline structure disappears during melting, the accumulated plastic shears $\gamma^a$ are reset to zero. Therefore, as described by Eq. (21), the critical resolved shear stresses $\tau_c^a(T, \gamma^a)$ reach minimum values immediately after solidification. This leads to a reduction in stress values immediately behind the melt pool.

As shown in Fig. 10, an asymmetry of the stress distribution is observed in the wake of the laser trajectory with respect to the scan line. In the left-side (toward the +y direction with respect to the scan line), the substrate is remelted for the first time and higher stresses are observed. In the right-side (toward the -y direction with respect to the scan line), the remelted zone of the previous pass (Pass 3) is remelted again. This asymmetry is observed for both $\sigma_{yy}$ (Fig. 10(c)) and $\sigma_{xx}$ (Fig. 10(d)). Note that one of the main differences between the left and right sides of the scan line is that a columnar structure forms on the right side by successive epitaxial growth from previous scans, whereas the initial equiaxed structure of the substrate is remelted on the left side. Regarding $\sigma_{xx}$ (Fig. 10(d)), a band pattern is visible, with bands parallel to the scan line, on its right side, in the neighborhood of the newly remelted zone, expressing the legacy of the previous remelting/solidification events during the previous scans. In regions

where marked columnar solidification takes place, such as in zoom ①, elongated grains can be observed in Fig. 10(b), which correlate strongly with the $\sigma_{yy}$ distribution.

The longitudinal section (Figs. 11(a-b)) is also noteworthy. Although at 13 ms the current grains appear to be uniformly distributed across the section, a high stress distribution with a quasi-periodic pattern is revealed in Fig. 11 at zoom ①. In fact, the structure in the upper region of Fig. 11(b) is inherited from melting and solidification during Pass 3 and is already columnar before melting in Pass 4, as seen in Fig. 10(b) at the same time 13 ms. After laser scanning, at time 15 ms, the $\sigma_{zz}$ component exhibits distinct discontinuities along the x-direction in Fig. 11(c) at zoom ②, strongly correlating with the columnar structure now visible in Fig. 11(b). Tensile and compressive $\sigma_{zz}$ concentration regions extend along elongated grains with different orientations. During solidification, tensile $\sigma_{zz}$ stresses develop below the melt pool and extend to the bottom of the substrate, transitioning to compressive regions further down. In contrast to the distribution of $\sigma_{zz}$, the distribution of the axial stress $\sigma_{xx}$ (Fig. 11(d)) exhibits minimal correlation with the underlying grains. Instead, $\sigma_{xx}$ forms a compressive region extending deeper below the melt pool, with tensile stresses appearing in its wake.

To investigate the relationship between stress and thermal history, a virtual sensor is considered, at the coordinates (0.7, 0.23, 0.0) mm, corresponding to the mid-point of Pass 1 in Case 1, to record the evolution of stress and temperature. In Fig. 12, a solid black line and a red dashed line show the evolution of the axial stress and temperature at the sensor, respectively. Each scan pass is outlined by a solid blue line to indicate its duration.

During Pass 1, as the laser beam approaches the sensor location, the axial stress becomes increasingly compressive (even before local heating), up to a peak just before melting. Over the temperature $T_{coh}$, the material behavior becomes Newtonian, resulting in a very low stress level. Subsequently, as the temperature decreases, the axial stress
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Evolution of (a) the temperature field, (b) the grain structure and the stress components (c) $\sigma_{zz}$ and (d) $\sigma_{xx}$ when the laser is at the beginning, midway and at the end of the Pass 4 in Case 1. The black solid line represents the critical temperature at which the mechanical behavior changes, $T_{coh}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Evolution of the stress component $\sigma_{xx}$ (black solid line) and the temperature (red dashed line) at point (0.7, 0.23, 0.0) mm for Case 1.</figcaption></figure>

becomes tensile, forming a peak before decreasing as the laser beam moves away. During Pass 2, the temperature and stress evolution show similar trends, including a shorter remelting period. During Passes 3 to 5, the increasing distance from the scan line mitigates the compression effect and prevents remelting. However, thermal diffusion continues to influence the stress evolution, which remains tensile.

### 4.3. Quantitative analysis of stress and grain structure

In this Section, the grain boundary character is quantified using the disorientation angle ($\Delta\theta$), which represents the minimum orientation difference between adjacent grains within the standard stereographic triangle. For cubic systems, $\Delta\theta$ ranges up to 62.8°. The analysis considers crystal symmetry operations to identify the fundamental zone representation of disorientations. A detailed mathematical formulation of the quaternion-based calculation appears in Appendix B.

#### 4.3.1. Accumulation of stress in cross section

To understand how the stress accumulates, we plot in Fig. 13 the $\sigma_{zz}$ component along the line $z = -0.036$ mm at the cross section yz plane ($x = 0.95$ mm) at the end of Passes 1 to 5 in Case 1 corresponding to times 4 ms, 8 ms, 12 ms, 16 ms and 20 ms, respectively. The black solid line represents the profile of $\sigma_{zz}$ stress. The remelted zone for each pass is visualized by vertical blue dotted lines. Along the black dashed line, we have selected Points 1 to 4, which are locations of accumulated stress at the end of each pass.

The most significant stress changes occur within the remelted zone. At $t = 4$ ms, Points 1 to 4 exhibit low and similar stress levels as they are far from the remelted zone. In the following Pass 2 ($t = 8$ ms) and Pass 3 ($t = 12$ ms), the stress previously built-up are affected by substantial changes within the remelted zone. After Pass 4 ($t = 16$ ms), the stress profile shows a highest value at Point 2. A progressive increase takes place as the laser approaches the location while it is followed by annealing in the wake of the laser (e.g. stress at Point 2 first increases
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Profiles of the stress component $\sigma_{zz}$ along the line $z = -0.036$ mm in the yz-plane defined by $x = 0.95$ mm at times 4 ms, 8 ms, 12 ms, 16 ms and 20 ms corresponding to the end of Pass 1 to 5. The profile obtained with the effective Case 1 grain structure is in plain black. The profile obtained for two single crystal orientations SX1 (0°, 0°, 0°) and SX2 (45°, 54.74°, 45°), are in plain red and plain blue, respectively. Labels 1 to 4, also reported in Fig. 14(a), permit correlation with the grain structure. The dashed line in the insert locates the position of the profiles in (plain black) the rectangular simulation domain. The vertical dotted blue lines in the graphs represent the limits of the remelted zones for each pass along the profile, with full geometry of the remelted zones in the transverse section sketched at the top of graphs.</figcaption></figure>

up to 500 MPa after melting and solidification at time 16 ms, and then decreases to around 350 MPa as Point 2 is only reheated during the last pass.) It is worth noting that the stress level is the same as in the work of [45], and also corresponds to the macro residual stress measurements on AM 316L parts for different geometries [82–84].

To investigate the influence of grain orientation on the stress distribution, the simulation is repeated with a uniform mechanical behavior corresponding to a fixed orientation of a single crystal (SX). Two orientations are considered. The first SX orientation (labeled SX1 in the figure) has Euler angles of (0°, 0°, 0°), indicating that the (100) crystal directions are all aligned with the *x*, *y*, and *z* axes. The second SX orientation (labeled SX2) has Euler angles of (45°, 54.74°, 45°), i.e. a (111) direction aligned with the *z* direction.

Among the three simulation cases, SX1 with (100) along the *z*-axis shows lower stress, while SX2 with (111) along the *z*-axis exhibits a higher magnitude of $\sigma_{zz}$. This is consistent with the 316L single crystal tensile tests by Karman et al. [85], where (111) gives a stronger stress response than (100). In addition, the stress concentrations are found to accumulate along the edges of the melt pool in each scanning track, indicating that the magnitude of the stress is influenced not only by the crystallographic orientation of the grains but also by the thermal history induced by laser scanning. The only difference between Case 1 and the single crystal cases is the polycrystalline grain structure. The significant variations between the stress profiles in Fig. 13 demonstrate that the orientation differences between the grains play an important role in the stress distribution.

### 4.3.2. Stress distribution and grain orientation

After the material has cooled down to temperature $T_{ext}$, stress heterogeneity remains between the different grains. Figs. 14 (a) to (d), show the grain structure and $\sigma_{zz}$ distribution for the 4 different cases at cross section $x = 0.95$ mm in the zoomed region $0.15 < y < 0.58$ mm and $-0.08 < z < 0$ mm. In all 4 cases, the positions of the selected points, labeled 1 to 6, are fixed but the corresponding grain number varies between cases. The grains and points are labeled in Figs. 14(a)–(d). For Case 4, additional Points 5 and 6 have been added due to the larger hatch spacing.

To quantitatively analyze the relationship between grain orientation and stress response, the stereographic triangle with respect to different directions is usually used to represent the one-point relationship between grain orientation and stress response. However, the texture formed in this work is weak as indicated by the texture index. The reason for this is that we only perform laser melting on the substrate. Therefore, the stress distribution within the stereographic triangle does not follow the strong trend observed in the literature [41,45]. Accordingly, in the following Sections, we will primarily investigate the influence of point disorientation on the stress difference [41].

We project the above points onto a stereographic triangle with respect to the *y* direction. The stress $\sigma_{zz}$ at each sample point as well as the disorientation angle between the corresponding grains are also labeled in Fig. 15. As shown in Fig. 15, Point 1 exhibits negative stress while Point 2 shows positive stress in all cases and the difference remains large in Cases 1 to 3 (as the colors are deep red and deep
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Distribution of (top) the grain structure and (bottom) the stress component $\sigma_{zz}$ at cross section yz plane ($x = 0.95$ mm) for (a) Case 1, (b) Case 2, (c) Case 3 and (d) Case 4. The region $0.15 \text{ mm} < y < 0.58 \text{ mm}$ and $-0.08 \text{ mm} < z < 0 \text{ mm}$ is chosen to observe the grain morphology, with points of interest labeled 1 to 6 within the grains labeled a to h along a unique dashed line defined by $z = -0.036 \text{ mm}$ in yz plane ($x = 0.95 \text{ mm}$). The positions of both the dashed line and (dotted rectangle) the region of the grain structures with respect to (plain rectangle) the simulation domain are schematized below the images. Points 1 to 4 in Fig. 13 are selected to quantitatively analyze the relationship between grain disorientation and stress distribution.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Concurrent representation of (colors) the stress component $\sigma_{zz}$ within a grain and (angle values) the disorientation for (connectors) different pairs of points using an IPF triangle with respect to the y-direction for (a) Case 1, (b) Case 2, (c) Case 3 and (d) Case 4. The selected points are the same as those shown in Fig. 14. Points and grains are labeled by a number and a letter, respectively.</figcaption></figure>

blue). The disorientation angle between Points 1 and 2 is 44.2° in Case 1–2 and 52.3° in Case 3, which belong to HAGB (above 15°). The different hatch spacing makes grain e appear at the position of Point 2 in Case 3. Although the corresponding grains have some differences, their distances on the stereographic triangle are very large, and their hue values are essentially opposite (i.e., they are located on opposite sides of the triangular white center). For Point 3–4, the disorientation angle between them is 54.7° in Case 1–2 and 57.8° in Case 3 belonging to HAGB. In Case 4, because of different hatch spacing, the stress contrasts between Points 1 and 2 and between Points 3 and 4 are not as obvious as in Cases 1 to 3. However, for Points 5 and 6, which have a large stress difference, they still follow this pattern, i.e. they have a large disorientation angle between them, and their hue values on the triangle are significantly different.

### 4.3.3. Statistical study

To statistically study the impact of disorientation angle on the stress distribution, we proceed with the following steps:

1. We set up a regular Cartesian grid of 215 × 40 points (8600 points) in the part of the transverse mid-section indicated by dotted lines at the bottom right of Fig. 13 ($0.15 \text{ mm} < y < 0.58 \text{ mm}$ and $-0.08 \text{ mm} < z < 0 \text{ mm}$) with a spacing of 2 µm in both the y and z directions. This spacing is smaller than the initial grain size of 10 µm.

2. From this grid, we select all pairs of points such that the two points are respectively located in neighboring grains (grains sharing a grain boundary), excluding pairs that fall within the same grain, with the additional requirement that their coordinate differences in both the y and z directions do not exceed 40 µm, a distance that approximately corresponds to the diameter of the largest grains observed in the remelted region.

3. We then calculate the disorientation angle $\Delta\theta$ and the stress difference $\Delta\sigma_{zz}$ for each pair of points. The range of disorientation angle is divided into four intervals ($0° - 15°$, $15° - 30°$, $30° - 45°$, $45° - 62.8°$). To focus the analysis on large stress differences, only
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Box plots of stress difference ($\Delta\sigma_{zz} > 400$ MPa) as a function of the grain boundary disorientation ($\Delta\theta > 0$) for point pairs, grouped in (0° – 15°, 15° – 30°, 30° – 45°, 45° – 62.8°). As indicated in the right bottom part of the figure, candidate points for point pairs are located in the part of the transverse mid-section indicated by dotted lines at the bottom right of Fig. 13 (0.15 mm < y < 0.58 mm and -0.08 mm < z < 0 mm). Selected point pairs are such that the two points respectively belong to neighboring grains and the distance between them does not exceed 40 µm. Box plots are shown for (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4. The blue boxes indicate the interquartile range (IQR, corresponding to the 25th–75th percentile), with whiskers extending to 1.5xIQR. Within the blue box, the horizontal black line shows the median value of the stress difference and the white square indicates the mean value.</figcaption></figure>

point pairs such that $\Delta\sigma_{zz} > 400$ MPa are retained, this value being arbitrarily considered as large.

Cases 1 and 2 involve five laser passes, while Cases 3 and 4 involve four. To neutralize this difference, the above operations are performed after the completion of the fourth scan track, i.e., at 16 ms, for all cases. Table 4 shows the number of point pairs, $N_{pp}$, with condition ($\Delta\theta > 0$) as well as $N_{pp}$ with the combined condition ($\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa). Under condition ($\Delta\theta > 0$), there are between 917807 and 947648 events for the 4 cases. The small variation reflects only a slight variation in the total length of the grain boundaries. However, under the combined conditions ($\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa), the $N_{pp}$ value for cases 1 to 4 is 10662, 7923, 3098, 1386, respectively. This reflects a significant decrease in the occurrence of high stress

differences at neighboring grains when moving from Case 1 to Case 4.

Fig. 16 shows the box plots of stress difference ($\Delta\sigma_{zz} > 400$ MPa) as a function of the disorientation ($\Delta\theta > 0$) for point pairs for the four cases. For all of them, almost no point pairs are located in the low disorientation angle region ($\Delta\theta \le 15°$). Most point pairs with $\Delta\sigma_{zz} > 400$ MPa are found on the high disorientation angles ($\Delta\theta > 15°$).

In Cases 3 and 4 (Figs. 16c and d), compared to Case 1 (Fig. 16(a)), the same unidirectional scanning strategy is used, while gradually increasing the hatch spacing. The maximum stress difference at high disorientation angle is found to be reduced. This suggests that adjusting the laser scanning strategy by increasing the hatch spacing can mitigate the stress variations within grain structures. This finding is consistent
<caption>**Table 4**</caption>

Statistical results for the 4 scanning strategies defined in Table 2 where $N_{pp}(\Delta\theta > 0)$ represents the number of point pairs with condition $\Delta\theta > 0$. $N_{pp}(\Delta\theta > 0$, $\Delta\sigma_{zz} > 400$ MPa) the number of point pairs of neighboring grains with condition $\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa. NMI is the Normalized Mutual Information index [87,88].

<table><thead><tr><th></th><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 4</th></tr></thead><tbody><tr><td>$N_{pp}(\Delta\theta > 0)$</td><td>947,648</td><td>917,807</td><td>941,420</td><td>938,217</td></tr><tr><td>$N_{pp}(\Delta\theta > 0$, $\Delta\sigma_{zz} > 400$ MPa)</td><td>10,662</td><td>7923</td><td>3096</td><td>1386</td></tr><tr><td>NMI</td><td>0.418</td><td>0.400</td><td>0.358</td><td>0.453</td></tr></tbody></table>

with the results reported by Balbaa and Elbestawi [86], even if their analysis was more focused on residual stresses.

We employ the Normalized Mutual Information (NMI) to quantify the correlation between variables $\Delta\theta$ and $\Delta\sigma_{zz}$ for point pairs with ($\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa). As proposed by Strehl and Ghosh [87,88], NMI is a normalized metric for measuring the amount of shared information between two variables. Its value ranges from 0 to 1, with values close to 1 indicating a higher degree of shared information between the two variables. Conversely, a value close to 0 means that there is no correlation. Table 4 reports the NMI analysis results for Case 1 (0.418), Case 3 (0.358), and Case 4 (0.453). These results indicate that $\Delta\theta$ and $\Delta\sigma_{zz}$ exhibit a certain degree of correlation, demonstrating the existence of a statistical relationship.

Finally, compared to Case 1, the back-and-forth scanning strategy used in Case 2 (Fig. 16(b) is found to slightly reduce the maximum stress difference. However, a high $N_{pp}$ value is still observed for stress differences greater than 400 MPa in comparison with Cases 3 and 4. Again, the NMI analysis for case 2 (0.400) reveals correlation between disorientation and stress difference.

As explained in Section 3.3, the substrate is characterized by a weak texture (white distribution in Fig. B.18). Despite the selection of grains shown in Fig. 9, texture remains weak. This is verified by comparing the distribution of $\Delta\theta$ of adjacent grains with the Mackenzie distribution [80] for all cases, as reported in Fig. B.18. But these distributions also mean that, for each of the disorientation classes, the number of point pairs differ: there is more probability to find a point pair with high disorientation (e.g. 45°) in comparison with a point pair with low disorientation (e.g. 5°). To eliminate the influence of the distribution of point pairs disorientation, we normalize the distributions of point pairs with $\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa by the total number of point pairs within each 5° interval, from 0° to 62.8°, reported in Fig. B.18, for each case. Note that this is very similar to using the Mackenzie distribution for normalization. In the absence of correlation between disorientation and stress difference, a uniform distribution is now expected. As shown in Fig. 17, the distribution is far from being uniform. In all 4 cases, the proportion of point pairs with $\Delta\sigma_{zz} > 400$ MPa is minimal when $\Delta\theta \le 25°$. When $\Delta\theta$ exceeds 25°, the proportion of point pairs increases, with maxima within the 30° to 35° classes and then for $\Delta\theta > 55°$. In Fig. 17, we also calculate the proportion of all point pairs with $\Delta\theta > 0$ and $\Delta\sigma_{zz} > 400$ MPa relative to the total number of point pairs $\Delta\theta > 0$, i.e. using values in Table 4. It is clearly observed that the proportion for Case 4 (0.15%) is lower compared to Case 3 (0.33%). These results express the effect of hatch spacing on the magnitude of stress difference at grain boundaries previously discussed. A larger remelted region between passes (Case 1) tends to increase the magnitude of stress and further amplifies the effect of disorientation between neighboring regions on stress differences. Case 2 (0.86%) is only slightly lower than Case 1 (1.13%), again showing the minor role of the back-and-forth scanning strategy in the present investigations.

## 5. Conclusion and perspectives

This study presents a computational framework to dynamically couple thermal resolution, grain growth, and mechanical response for

<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Proportion of point pairs with stress variation larger than 400 MPa to the total number of point pairs for Cases 1 to 4.</figcaption></figure>

multi-track laser scanning of polycrystalline 316L, where a crystal visco-plasticity model is integrated with the evolved grain structure predicted by the CAFE model to simulate the mechanical response during process. The developed model can predict the evolution of stress and strain at the grain scale during laser scanning processes. The simulation results elucidate the relationship between stress evolution and the orientation of grains, as well as the influence of different scanning strategies on grain structure and residual stresses.

The computational framework is based on a thermo-elasto-viscoplastic model, in small deformation. The crystal viscoplasticity model was calibrated using stress-strain curves obtained from two distinct experiments conducted on heat-treated 316L specimens.

Using this computational framework, four different laser scanning cases were studied to investigate the relationship between grain structure and residual stress under various conditions. The simulation results reveal the mechanisms of stress evolution and how they relate to the laser scanning process. Around the melt pool, laser scanning induces complex patterns of compressive and tensile stresses. The results also highlight the relationship between stress distribution and grain orientation. Notably, stress concentrations observed in cross-sections occur preferentially where adjacent grains with significantly different orientations come into contact during competitive grain growth. This trend is more pronounced when the hatch spacing is increased, although the maximum stress concentration is then reduced.

The framework provides robust support for microscale mechanical analysis in the field of additive manufacturing and offers a valuable computational tool for multi-scale process-structure-property investigations. Furthermore, the directionality of the residual stresses associated with the laser tracks may provide a basis for studying the mechanisms of micro-crack formation in additive manufacturing.

## CRediT authorship contribution statement

Zixuan Li: Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. Michel Bellet: Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Investigation, Conceptualization. Charles-André Gandin: Writing – review & editing, Supervision, Project administration, Methodology, Conceptualization. Manas Vijay Upadhyay: Writing – review & editing, Supervision, Project administration, Conceptualization. Yancheng Zhang: Writing – review & editing, Validation, Supervision, Software, Methodology, Investigation, Conceptualization.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are grateful to Dr. Juan Guillermo Santos Macías (currently at IMDEA Materiales, Spain) for providing the experimental results from [37] that were used to validate the model predictions in the present work. The authors thank the Scientific Interest Group “High Energy in Additive Manufacturing” (GIS HEAD) for its support through the grant GISTEAM3. MVU is grateful to the European Research Council (ERC) for their support through the European Union’s Horizon 2020 – EXCELLENT SCIENCE – research and innovation program (grant agreement number 946959) and to LMS, United Kingdom for providing financial support to build the CW laser-SEM system. The authors thank the French National Research Agency (ANR) for the financial support of the fundamental research project GRAMME within the frame of the JCJC program (ANR-21-CE08-0031).

Appendix A. Obtention of the expression of the CEVP compliance modulus $\partial\dot{\gamma}^\alpha/\partial\tilde{s}$

For a given slip system $\alpha$ (hereafter omitted for brevity), we start from Eq. (19), repeated here:

$$ \dot{\gamma} = \dot{\gamma}_0 \left[ \frac{\|\mathbf{T}_{RSS}\| - \tau_c(T)}{K(T)} \right]_+^{\frac{1}{m(T)}} \frac{T_{RSS}}{\|\mathbf{T}_{RSS}\|} \cdot \mathbf{m} \quad (A.1) $$

Deriving with respect to component $\tilde{s}_j$, we have:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \dot{\gamma}_0 \left[ \frac{\left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m}} \frac{\partial}{\partial \tilde{s}_j} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+ + \frac{\mathbf{T}_{RSS} \cdot \mathbf{m}}{\|\mathbf{T}_{RSS}\|}}{K^{1/m}} + \frac{\left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m}} \frac{\partial}{\partial \tilde{s}_j} \left( \frac{\mathbf{T}_{RSS} \cdot \mathbf{m}}{\|\mathbf{T}_{RSS}\|} \right)}{K^{1/m}} \right] \quad (A.2) $$

First, let us show the last derivative is null:

$$ \frac{\partial}{\partial \tilde{s}_j} \left( \frac{\mathbf{T}_{RSS} \cdot \mathbf{m}}{\|\mathbf{T}_{RSS}\|} \right) = - \|\mathbf{T}_{RSS}\|^{-2} \frac{\partial \|\mathbf{T}_{RSS}\|}{\partial \tilde{s}_j} (\mathbf{T}_{RSS} \cdot \mathbf{m}) + \frac{1}{\|\mathbf{T}_{RSS}\|} \frac{\partial}{\partial \tilde{s}_j} (\mathbf{T}_{RSS} \cdot \mathbf{m}) \quad (A.3) $$

Using Eq. (20), we have:

$$ - \|\mathbf{T}_{RSS}\|^{-2} \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{P} : \mathbf{s}) \tilde{P}_j (\mathbf{T}_{RSS} \cdot \mathbf{m}) + \frac{1}{\|\mathbf{T}_{RSS}\|} \tilde{P}_j = 0 \quad (A.4) $$

This leads to:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{m K^{1/m}} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m} - 1} \frac{\partial}{\partial \tilde{s}_j} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+ + \frac{\mathbf{T}_{RSS} \cdot \mathbf{m}}{\|\mathbf{T}_{RSS}\|} \quad (A.5) $$

Let us calculate now $\frac{\partial}{\partial \tilde{s}_j} [\|\mathbf{T}_{RSS}\| - \tau_c]_+$. There are three cases, according to the approximation of the Macauley bracket (Eq. (28)).

• Case 1: $-\beta\tau_c \le \|\mathbf{T}_{RSS}\| - \tau_c < +\beta\tau_c$

$$ \frac{\partial}{\partial \tilde{s}_j} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+ = \frac{\partial}{\partial \tilde{s}_j} \left( \frac{\|\mathbf{T}_{RSS}\| - \tau_c}{2} + \frac{(\|\mathbf{T}_{RSS}\| - \tau_c)^2}{4\beta\tau_c} + \frac{\beta\tau_c}{4} \right) \quad (A.6) $$

After some calculation, we get:

$$ \frac{\partial}{\partial \tilde{s}_j} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+ = \frac{1}{2\beta\tau_c} \frac{\partial \|\mathbf{T}_{RSS}\|}{\partial \tilde{s}_j} (\|\mathbf{T}_{RSS}\| - \tau_c(1 - \beta)) \quad (A.7) $$

Using again Eq. (20), we have

$$ \frac{\partial}{\partial \tilde{s}_j} [\|\mathbf{T}_{RSS}\| - \tau_c]_+ = \frac{1}{2\beta\tau_c} (\|\mathbf{T}_{RSS}\| - \tau_c(1 - \beta)) \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{P} : \mathbf{s}) \tilde{P}_j \quad (A.8) $$

Injecting this result in Eq. (A.5), we get:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{2\beta\tau_c m K^{1/m}} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m} - 1} (\|\mathbf{T}_{RSS}\| - \tau_c(1 - \beta)) \times \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{P} : \mathbf{s}) \tilde{P}_j \quad (A.9) $$

$$ (\mathbf{P} : \mathbf{s}) \tilde{P}_j \frac{1}{\|\mathbf{T}_{RSS}\|^2} (\mathbf{T}_{RSS} \cdot \mathbf{m}) $$

After simplification, we get finally:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{2\beta\tau_c m K^{1/m}} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m} - 1} (\|\mathbf{T}_{RSS}\| - \tau_c(1 - \beta)) \tilde{P}_j \quad (A.10) $$

• Case 2: $\|\mathbf{T}_{RSS}\| - \tau_c < -\beta\tau_c$, that is $\|\mathbf{T}_{RSS}\| < (1 - \beta)\tau_c$

In this case, $[\|\mathbf{T}_{RSS}\| - \tau_c]_+ = 0$ and $\frac{\partial}{\partial \tilde{s}_j} [\|\mathbf{T}_{RSS}\| - \tau_c]_+ = 0$

So, we have:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = 0 \quad (A.11) $$

• Case 3: $\|\mathbf{T}_{RSS}\| - \tau_c > \beta\tau_c$, that is $\|\mathbf{T}_{RSS}\| > (1 + \beta)\tau_c$

In this case, $[\|\mathbf{T}_{RSS}\| - \tau_c]_+ = \|\mathbf{T}_{RSS}\| - \tau_c$

So, we have:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{m K^{1/m}} \left[ \|\mathbf{T}_{RSS}\| - \tau_c \right]_+^{\frac{1}{m} - 1} \frac{\partial}{\partial \tilde{s}_j} (\|\mathbf{T}_{RSS}\| - \tau_c) \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{T}_{RSS} \cdot \mathbf{m}) \quad (A.12) $$

Using again Eq. (20), we get:

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{m K^{1/m}} (\|\mathbf{T}_{RSS}\| - \tau_c)^{\frac{1}{m} - 1} \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{P} : \mathbf{s}) \tilde{P}_j \frac{1}{\|\mathbf{T}_{RSS}\|} (\mathbf{T}_{RSS} \cdot \mathbf{m}) \quad (A.13) $$

Finally,

$$ \frac{\partial \dot{\gamma}}{\partial \tilde{s}_j} = \frac{\dot{\gamma}_0}{m K^{1/m}} (\|\mathbf{T}_{RSS}\| - \tau_c^\alpha)^{\frac{1}{m} - 1} \tilde{P}_j \quad (A.14) $$

Comparing Eqs. (A.9), (A.11), and (A.14), the continuity of the derivatives at both ends of the interval $[-\beta\tau_c, +\beta\tau_c]$ is easily observed. This is clearly a consequence of the continuity of the derivative of the approximation function defined in Eq. (28).

Appendix B. Crystal symmetry in quaternion and Mackenzie distribution for disorientation angle

B.1. Symmetry operators

The orientation difference between two grains is calculated using quaternion algebra. A quaternion $q = \{q_0; q\}$ represents crystal orientation, where $q_0 = \cos(\omega/2)$ and $q = \sin(\omega/2)$ $n$, with $\omega$ being the rotation angle and $n = (n_1, n_2, n_3)$ the rotation axis related to Euler angles $(\phi_1, \phi, \phi_2)$:

$$ q = \{ \cos(\phi/2) \cos(\frac{\phi_1+\phi_2}{2}); \sin(\phi/2) \cos(\frac{\phi_1-\phi_2}{2}), \sin(\phi/2) \cos(\frac{\phi_1-\phi_2}{2}) \} \quad (B.1) $$

The quaternion difference $\Delta q$ between two grains A and B is:

$$ \Delta q(q^A, q^B) = (q^A)^{-1} \cdot q^B = q_0^A q_0^B + q^A \cdot q^B; q_0^A q^B - q_0^B q^A - q^A \times q^B \quad (B.2) $$

where $q^{-1} = \{q_0; -q\}$. Considering crystal symmetry operators $S$, the disorientation $D = \Delta q(S^A \cdot q^A, (S^B \cdot q^B))$ gives the minimum rotation angle in the fundamental zone, with $\Delta\theta = \arccos(2D_0^2 - 1)$.

For cubic system, there are 24 crystal symmetry in quaternion (see Box II):
$$
\begin{align*}
S_1 &= (1; 0, 0, 0) & S_2 &= (0; 1, 0, 0) & S_3 &= (0; 0, 1, 0) & S_4 &= (0; 0, 0, 1) \\
S_5 &= (0.5; 0.5, 0.5, 0.5) & S_6 &= (0.5; -0.5, -0.5, -0.5) & S_7 &= (0.5; 0.5, -0.5, 0.5) & S_8 &= (0.5; -0.5, 0.5, -0.5) \\
S_9 &= (0.5; -0.5, 0.5, 0.5) & S_{10} &= (0.5; 0.5, -0.5, -0.5) & S_{11} &= (0.5; 0.5, 0.5, -0.5) & S_{12} &= (0.5; -0.5, -0.5, 0.5) \\
S_{13} &= (1/\sqrt{2}; 1/\sqrt{2}, 0, 0) & S_{14} &= (1/\sqrt{2}; 0, 1/\sqrt{2}, 0) & S_{15} &= (1/\sqrt{2}; 0, 0, 1/\sqrt{2}) & S_{16} &= (1/\sqrt{2}; -1/\sqrt{2}, 0, 0) \\
S_{17} &= (1/\sqrt{2}; 0, -1/\sqrt{2}, 0) & S_{18} &= (1/\sqrt{2}; 0, 0, -1/\sqrt{2}) & S_{19} &= (0; 1/\sqrt{2}, 1/\sqrt{2}, 0) & S_{20} &= (0; -1/\sqrt{2}, 1/\sqrt{2}, 0) \\
S_{21} &= (0; 0, 1/\sqrt{2}, 1/\sqrt{2}) & S_{22} &= (0; 0, -1/\sqrt{2}, 1/\sqrt{2}) & S_{23} &= (0; 1/\sqrt{2}, 0, 1/\sqrt{2}) & S_{24} &= (0; -1/\sqrt{2}, 0, 1/\sqrt{2}) \tag{B.3}
\end{align*}
$$

Box II.

<figure><img src="image_21.png" /><figcaption>Fig. B.18. Distribution of disorientation angle $\Delta\theta$ for adjacent grains in the selected regions of Fig. 14 (dotted rectangle) and comparison with the Mackenzie distribution.</figcaption></figure>

B.2. Mackenzie distribution

For random textures, the distribution disorientation angles for neighboring grains follows the Mackenzie function [80]. We conducted a study on the disorientation angle distribution of all adjacent grains in the selected region of Fig. 14 (dotted rectangle), as shown in the histogram in Fig. B.18. The same analysis is performed for the initial grain structure reported in Fig. 7. We can observe that for all cases, their $\Delta\theta$ distributions are essentially similar to those of a arbitrary texture. This implies that there is essentially no texture formation on this selected region. The probability density of the $\Delta\theta$ gradually increases and reaches its maximum around 45°.

Data availability

Data will be made available on request.

References

[1] Y. Zhang, Q. Chen, G. Guillemot, C.-A. Gandin, M. Bellet, Numerical modelling of fluid and solid thermomechanics in additive manufacturing by powder-bed fusion: Continuum and level set formulation applied to track- and part-scale simulations, C. R. Méc. 346 (11) (2018) 1055-1071, http://dx.doi.org/10.1016/j.cme.2018.08.008.

[2] Y. Li, D. Gu, Parametric analysis of thermal behavior during selective laser melting additive manufacturing of aluminum alloy powder, Mater. Des. 63 (2014) 856-867, http://dx.doi.org/10.1016/j.matdes.2014.07.006.

[3] Z. Zhang, Y. Wang, P. Ge, T. Wu, A review on modelling and simulation of laser additive manufacturing: Heat transfer, microstructure evolutions and mechanical properties, Coatings 12 (9) (2022) http://dx.doi.org/10.3390/coatings12091277.

[4] Q. Han, R. Mertens, M.L. Montero-Sistiaga, S. Yang, R. Setchi, K. Vanmeensel, B. Van Hooreweder, S.L. Evans, H. Fan, Laser powder bed fusion of Hastelloy X: Effects of hot isostatic pressing and the hot cracking mechanism, Mater. Sci. Eng.: A 732 (2018) 228-239, http://dx.doi.org/10.1016/j.msea.2018.07.008.

[5] T. Kannengiesser, T. Boellinghaus, Cold cracking tests—an overview of present technologies and applications, Weld, the World 57 (2013) 3-37, http://dx.doi.org/10.1007/s40194-012-0001-7.

[6] M. Rappaz, J.-M. Drezet, M. Gremaud, A new hot-tearing criterion, Met. Mater. Trans. A 30 (2) (1999) 449-455, http://dx.doi.org/10.1007/s11661-999-0334-z.

[7] P. Kontis, E. Chauvet, Z. Peng, J. He, A.K. da Silva, D. Raabe, C. Tassin, J.-J. Blandin, S. Abed, R. Dendievel, B. Gault, G. Martin, Atomic-scale grain boundary engineering to overcome hot-cracking in additively-manufactured superalloys, Acta Mater. 177 (2019) 209-221, http://dx.doi.org/10.1016/j.actamat.2019.07.041.

[8] E. Chauvet, P. Kontis, E.A. Jägle, B. Gault, D. Raabe, C. Tassin, J.-J. Blandin, R. Dendievel, B. Vayre, S. Abed, Hot cracking mechanism affecting a non-weldable Ni-based superalloy produced by selective electron Beam Melting, Acta Mater. 142 (2018) 82-94, http://dx.doi.org/10.1016/j.scriptamat.2018.03.041.

[9] N. Wang, S. Mokadem, M. Rappaz, W. Kurz, Solidification cracking of superalloy single- and bi-crystals, Acta Mater. 52 (11) (2004) 3173-3182, http://dx.doi.org/10.1016/j.actamat.2004.03.047.

[10] W. Yan, S. Lin, O.L. Kafka, Y. Lian, C. Yu, Z. Liu, J. Yan, S. Wolff, H. Wu, E. N dip-Agbor, M. Mozaffar, K. Ehmann, J. Cao, G.J. Wagner, W.K. Liu, D ata-driven multi-scale multi-physics models to derive process-structure-property relationships for additive manufacturing, Comput. Mech. 61 (5) (2018) 521-541, http://dx.doi.org/10.1007/s00466-018-1539-z.

[11] Q. Wang, S. Zhang, C. Zhang, C. Wu, J. Wang, J. Chen, Z. Sun, Microstructure evolution and EBSD analysis of a graded steel fabricated by laser additive manufacturing, Vacuum 141 (2017) 68-81, http://dx.doi.org/10.1016/j.vacuum.2017.03.021.

[12] C. Körner, M. Markl, J. Koepf, Modeling and simulation of microstructure evolution for additive manufacturing of metals: a critical review, Met. Mater. Trans. 51A (2020) 4970-4983, http://dx.doi.org/10.1007/s11661-020-05946-3.

[13] T.M. Rodgers, J.D. Madison, V. Tikare, M. Maguire, Predicting mesoscale microstructural evolution in electron beam welding, JOM 68 (2016) 1419-1426, http://dx.doi.org/10.1007/s11837-016-1863-8.

[14] T.M. Rodgers, J.D. Madison, V. Tikare, Simulation of metal additive manufacturing microstructures using kinetic Monte Carlo, Comput. Mater. Sci. 135 (2017) 78-89, http://dx.doi.org/10.1016/j.commatsci.2017.03.053.

[15] T.M. Rodgers, D. Moser, F. Abdeljawad, O. Underwood Jackson, J. Carroll, B. Jared, D. Bolintineanu, J. Mitchell, J.D. Madison, Simulation of powder metal additive manufacturing microstructures with coupled finite difference-Monte Carlo method, Addit. Manuf. 41 (2021) 101953, http://dx.doi.org/10.1016/j.addma.2021.101953.

[16] S. Sunny, H. Yu, R. Mathews, A. Malik, W. Li, Improved grain structure prediction in metal additive manufacturing using a dynamic kinetic Monte Carlo framework, Addit. Manuf. 37 (2021) 101649, http://dx.doi.org/10.1016/j.addma.2020.101649.

[17] J.A. Koepf, M.R. Gotterbarm, M. Markl, C. Körner, 3D multi-layer grain structure simulation of powder bed fusion additive manufacturing, Acta Mater. 152 (2018) 119-126, http://dx.doi.org/10.1016/j.actamat.2018.04.030.

[18] A. Rai, M. Markl, C. Körner, A coupled Cellular Automaton-Lattice Boltzmann model for grain structure simulation during additive manufacturing, Comput. Mater. Sci. 124 (2016) 37-48, http://dx.doi.org/10.1016/j.commatsci.2016.07.005.

[19] A. Zinoviev, O. Zinovieva, V. Ploshikhin, V. Romanova, R. Balokhonov, Evolution of grain structure during laser additive manufacturing. Simulation by a cellular automata method, Mater. Des. 106 (2016) 321-329, http://dx.doi.org/10.1016/j.matdes.2016.05.125.

[20] O. Zinovieva, A. Zinoviev, V. Ploshikhin, Three-dimensional modeling of the microstructure evolution during metal additive amnufacturing, Comput. Mater. Sci. 141 (2018) 207-220, http://dx.doi.org/10.1016/j.commatsci.2017.09.018.

[21] J. Yang, H. Yu, H. Yang, F. Li, Z. Wang, X. Zeng, Prediction of microstructure in selective laser melted Ti-6Al-4V alloy by cellular automaton, J. Alloys Compd. 748 (2018) 281-290, http://dx.doi.org/10.1016/j.matdes.2016.06.117.
[22] J.A. Koepf, D. Soldner, M. Ramsperger, J. Mergheim, M. Markl, C. Körner, Numerical microstructure prediction by a coupled finite element cellular automaton model for selective electron beam melting, Comput. Mater. Sci. 162 (2019) 148-155, http://dx.doi.org/10.1016/j.commatsci.2019.03.004.

[23] Y. Lian, Z. Gan, C. Yu, D. Kats, W.K. Liu, G.J. Wagner, A cellular automaton finite volume method for microstructure evolution during additive manufacturing, Mater. Des. 169 (2019) 107672, http://dx.doi.org/10.1016/j.matdes.2019.107672.

[24] Y. Zhang, J. Zhang, Modeling of solidification microstructure evolution in laser powder bed fusion fabricated 316L stainless steel using combined computational fluid dynamics and cellular automaton, Addit. Manuf. 28 (2019) 750-765, http://dx.doi.org/10.1016/j.addma.2019.06.024.

[25] M. Mohebbi, V. Ploshikhin, Implementation of nucleation in cellular automaton simulation of microstructural evolution during additive manufacturing of Al alloys, Addit. Manuf. 36 (2020) 101726, http://dx.doi.org/10.1016/j.addma.2020.101726.

[26] O. Zinovieva, A. Zinoviev, V. Romanova, R. Balokhonov, Three-dimensional analysis of grain structure and texture of additively manufactured 316L austenitic stainless steel, Addit. Manuf. 36 (2020) 101521, http://dx.doi.org/10.1016/j.addma.2020.101521.

[27] K. Teferra, D.J. Rowenhorst, Optimizing the cellular automata finite element model for additive manufacturing to simulate large microstructures, Acta Mater. 213 (2021) 116930, http://dx.doi.org/10.1016/j.actamat.2021.116930.

[28] T. Camus, D. Maisonnette, O. Baulin, O. Senninger, G. Guillemot, C.-A. Gandin, Three-dimensional modeling of solidification grain structures generated by laser powder bed fusion, Materialia 30 (2023) 101804, http://dx.doi.org/10.1016/j.mtla.2023.101804.

[29] Y. Zhang, G. Guillemot, T. Camus, O. Senninger, M. Bellet, C.-A. Gandin, Part-scale thermomechanical and grain structure modeling for additive manufacturing: status and perspectives, Metals 14 (2024) 1173, http://dx.doi.org/10.3390/met14101173.

[30] W. Kurz, B. Giovanolla, R. Trivedi, Theory of microstructural development during rapid solidification, Acta Metall. 34 (1986) 823-830, http://dx.doi.org/10.1016/0001-6160(86)90056-8.

[31] O. Hunziker, Theory of plane front and dendritic growth in multicomponent alloys, Acta Mater. 49 (20) (2001) 4191-4203, http://dx.doi.org/10.1016/S1359-6454(01)00313-5.

[32] G. Guillemot, O. Senninger, C.A. Hareland, P.W. Voorhees, C.-A. Gandin, Thermodynamic coupling in the computation of dendrite growth kinetics for multicomponent alloys, Calphad: Comput. Coupling Phase Diagrams Thermochem. 77 (2022) 102429, http://dx.doi.org/10.1016/j.calphad.2022.102429.

[33] A. Pineau, G. Guillemot, D. Tourret, A. Karma, C.-A. Gandin, Growth competition between columnar dendritic grains - cellular automaton versus phase field modeling, Acta Mater. 155 (2018) 286-301, http://dx.doi.org/10.1016/j.actamat.2018.05.032.

[34] E. Dorari, K. Ji, G. Guillemot, C.-A. Gandin, A. Karma, Growth competition between columnar dendritic grains - the role of microstructural length scales, Acta Mater. 223 (2022) 117395, http://dx.doi.org/10.1016/j.actamat.2021.117395.

[35] S. Sahoo, K. Chou, Phase-field simulation of microstructure evolution of Ti-6Al-4V in electron beam additive manufacturing process, Addit. Manuf. 9 (2016) 14-24, http://dx.doi.org/10.1016/j.addma.2015.12.005.

[36] S.M. Elahi, R. Tavakoli, I. Romero, D. Tourret, Grain growth competition during melt pool solidification - comparing phase-field and cellular automaton models, Comput. Mater. Sci. 216 (2023) 111882, http://dx.doi.org/10.1016/j.commatsci.2022.111882.

[37] A.F. Chadwick, J.G. Santos Macías, A. Samaei, G.J. Wagner, M.V. Upadhyay, P.W. Voorhees, On microstructure development during laser melting and resolidification: An experimentally validated simulation study, Acta Mater. 282 (2025) 120482, http://dx.doi.org/10.1016/j.actamat.2024.120482.

[38] J.G. Santos Macías, K. Chen, A. Tanguy, N. Isac, M. Vallet, L. Cornet, V. Michel, M.V. Upadhyay, High-vacuum laser treatments enhance strength, ductility and fatigue limit of additively manufactured stainless steel, Mater. Des. 254 (2025) 114064, http://dx.doi.org/10.1016/j.matdes.2025.114064.

[39] L.L. Beghini, M. Stender, D. Moser, B.L. Trembacki, M.G. Veilleux, K.R. Ford, A coupled fluid-mechanical workflow to simulate the directed energy deposition additive manufacturing process, Comput. Mech. 67 (4) (2021) 1041-1057, http://dx.doi.org/10.1007/s00466-020-01960-9.

[40] R. Darabi, E. Azinpour, A. Reis, J.C. de Sa, Elasto-plastic phase-field analysis of thermal induced-cracking and its application towards metal additive manufacturing, Eur. J. Mech. A Solids 107 (2024) 105369, http://dx.doi.org/10.1016/j.euromechsol.2024.105369.

[41] D. Hu, N. Grilli, L. Wang, M. Yang, W. Yan, Microstructure residual stresses in additively manufactured stainless steel: Computational simulation, J. Mech. Phys. Solids 161 (2022) 104822, http://dx.doi.org/10.1016/j.jmps.2022.104822.

[42] N. Grilli, D. Hu, D. Yushu, F. Chen, W. Yan, Crystal plasticity model of residual stress in additive manufacturing using the element elimination and reactivation method, Comput. Mech. 69 (3) (2022) 825-845, http://dx.doi.org/10.1007/s00466-021-02116-z.

[43] M. Lindroos, T. Pinomaa, A. Antikainen, J. Lagerbom, J. Reijonen, T. Lindroos, T. Andersson, A. Laukkanen, Micromechanical modeling approach to single track deformation, phase transformation and residual stress evolution during selective laser melting using crystal plasticity, Addit. Manuf. 38 (2021) 101819, http://dx.doi.org/10.1016/j.addma.2020.101819.

[44] N. Mohanan, J.G.S. Macías, J. Bleyer, T. Helfer, M.V. Upadhyay, Intergranular stress and plastic strain formation during laser scanning of additively manufactured stainless steel: An experimentally-driven thermomechanical simulation study, Materialia 34 (2024) 102082, http://dx.doi.org/10.1016/j.mtla.2024.102082.

[45] L. Kuna, A.J. Birnbaum, K. Teferra, A framework for predicting residual mechanical response fields for laser powder bed fusion at the polycrystalline length scale, Addit. Manuf. 78 (2023) 103830, http://dx.doi.org/10.1016/j.addma.2023.103830.

[46] T. Koshikawa, C.-A. Gandin, M. Bellet, H. Yamamura, M. Bobadilla, Computation of phase transformation paths in steels by a combination of the partial- and para-equilibrium thermodynamic approximations, ISIJ Int. 54 (6) (2014) 1274-1282, http://dx.doi.org/10.2355/isijinternational.54.1274.

[47] Thermo-Calc Software, TCFE Steels/Fe-alloys Database version 9.3, 2020, https://thermocalc.com/products/databases/steel-and-fe-alloys/.

[48] C. Xue, N. Blanc, F. Soulié, C. Bordreuil, F. Deschaux-Beaume, G. Guillemot, M. Bellet, C.-A. Gandin, Structure and texture simulations in fusion welding processes - comparison with experimental data, Materialia 21 (2022) 101305, http://dx.doi.org/10.1016/j.mtla.2021.101305.

[49] X. Li, F. Gao, J. Jiao, G. Cao, Y. Wang, Z. Liu, Influences of cooling rates on delta ferrite of nuclear power 316H austenitic stainless steel, Mater. Charact. 174 (2021) 111029, http://dx.doi.org/10.1016/j.matchar.2021.111029.

[50] S. Chen, G. Guillemot, C.-A. Gandin, Three-dimensional cellular automaton-finite element modeling of solidification grain structures for arc-welding processes, Acta Mater. 115 (2016) 448-467, http://dx.doi.org/10.1016/j.actamat.2016.05.011.

[51] C.-A. Gandin, M. Rappaz, A 3D cellular automaton algorithm for the prediction of dendritic grain growth, Acta Mater. 45 (5) (1997) 2187-2195, http://dx.doi.org/10.1016/S1359-6454(96)00303-5.

[52] H.-J. Bunge, Texture Analysis in Materials Science: Mathematical Methods, Elsevier, 2013, http://dx.doi.org/10.1016/C2013-0-11769-2.

[53] J.A. Dantzig, M. Rappaz, Solidification, Revised & Expanded Edition, EPFL Press, 2016.

[54] P. Vicente-Hernandez, Vers une approche thermomécanique du procédé de fonderie: caractérisation rhéologique et microstructurale d'alliages métalliques à l'état semi-solide (Ph.D. thesis), École Nationale Supérieure des Mines de Paris, 1994.

[55] M. Bellet, V.D. Fachinotti, ALE method for solidification modelling, Comput. Methods Appl. Mech. Engrg. 193 (39) (2004) 4355-4381, http://dx.doi.org/10.1016/j.cma.2003.11.016.

[56] M. Bellet, O. Jaouen, I. Poitault, An ALE-FEM approach to the thermomechanics of solidification processes with application to the prediction of pipe shrinkage, Internat. J. Numer. Methods Heat Fluid Flow 15 (2) (2005) 120-142, http://dx.doi.org/10.1108/09615530510578410.

[57] D. Peirce, R. Asaro, A. Needleman, An analysis of nonuniform and localized deformation in ductile single crystals, Acta Metall. 30 (6) (1982) 1087-1119, http://dx.doi.org/10.1016/0001-6160(82)90005-0.

[58] G. Cailletaud, J.-L. Chaboche, S. Forest, L. Rémy, On the design of single crystal turbine blades, Met. Res. Technol. 100 (2) (2003) 165-172, http://dx.doi.org/10.1051/metal:2003137.

[59] A. Musienko, Plasticité cristalline en présence de grandes déformations et d'endommagement (Ph.D. thesis), École Nationale Supérieure des Mines de Paris, 2005.

[60] M. Lindroos, T. Pinomaa, K. Ammar, A. Laukkanen, N. Provatas, S. Forest, Dislocation density in cellular rapid solidification using phase field modeling and crystal plasticity, Int. J. Plast. 148 (2022) 103139, http://dx.doi.org/10.1016/j.ijplas.2021.103139.

[61] B. Clausen, T. Lorentzen, T. Leffers, Self-consistent modelling of the plastic deformation of f.c.c. polycrystals and its implications for diffraction measurements of internal stresses, Acta Mater. 46 (9) (1998) 3087-3098, http://dx.doi.org/10.1016/S1359-6454(98)00014-7.

[62] R. Hill, Generalized constitutive relations for incremental deformation of metal crystals by multislipj, J. Mech. Phys. Solids 14 (2) (1966) 95-102, http://dx.doi.org/10.1016/0022-5096(66)90040-8.

[63] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (4) (1984) 337-344, http://dx.doi.org/10.1007/BF02576171.

[64] M. Fortin, A. Fortin, Experiments with several elements for viscous incompressible flows, Internat. J. Numer. Methods Fluids 5 (10) (1985) 911-928, http://dx.doi.org/10.1002/fld.1650051005.

[65] S. Zhang, G. Guillemot, C.-A. Gandin, M. Bellet, A partitioned two-step solution algorithm for concurrent fluid flow and stress-strain numerical simulation in solidification processes, Comput. Methods Appl. Mech. Engrg. 356 (2019) 294-324, http://dx.doi.org/10.1016/j.cma.2019.07.006.

[66] O. Jaouen, Modélisation tridimensionnelle par éléments finis pour l'analyse thermo-mécanique du refroidissement des pièces coulées, No. 1998ENMP0826 (Ph.D. thesis), École Nationale Supérieure des Mines de Paris, 1998.
[67] S.N. Kuchnicki, A.M. Cuitiño, R.A. Radovitzky, Efficient and robust constitutive integrators for single-crystal plasticity modeling, Int. J. Plast. 22 (10) (2006) 1988–2011, http://dx.doi.org/10.1016/j.ijplas.2006.02.008.

[68] K. Terada, I. Watanabe, Computational aspects of tangent moduli tensors in rate-independent crystal elastoplasticity, Comput. Mech. 40 (3) (2007) 497–511, http://dx.doi.org/10.1007/s00466-006-0123-0.

[69] J.-M. Scherer, V. Phalke, J. Besson, S. Forest, J. Hure, B. Tanguy, Lagrange multiplier based vs micromorphic gradient-enhanced rate-(in)dependent crystal plasticity modelling and simulation, Comput. Methods Appl. Mech. Engrg. 372 (2020) 113426, http://dx.doi.org/10.1016/j.cma.2020.113426.

[70] X. Zhang, C. Oskay, Polycrystal plasticity modeling of nickel-based superalloy IN 617 subjected to cyclic loading at high temperature, Modelling Simul. Mater. Sci. Eng. 24 (5) (2016) 055009, http://dx.doi.org/10.1088/0965-0393/24/5/055009.

[71] M.R. Daymond, P.J. Bouchard, Elastoplastic deformation of 316 stainless steel under tensile loading at elevated temperatures, Met. Mater. Trans. A 37 (6) (2006) 1863–1873, http://dx.doi.org/10.1007/s11661-006-0129-4.

[72] W. Chen, T. Voisin, Y. Zhang, J.-B. Forien, C.M. Spadaccini, D.L. McDowell, T. Zhu, Y.M. Wang, Microscale residual stresses in additively manufactured stainless steel, Nat. Commun. 10 (1) (2019) 4338, http://dx.doi.org/10.1038/s41467-019-12265-8.

[73] L. Depradeux, Simulation numérique du soudage-acier 316L: validation sur cas tests de complexité croissante (Ph.D. thesis), INSA, Lyon, 2004.

[74] U.F. Kocks, Laws for work-hardening and low-temperature creep, J. Eng. Mater. Technol. 98 (1) (1976) 76–85, http://dx.doi.org/10.1115/1.3443340.

[75] E.D. Cyr, M. Mohammadi, R.K. Mishra, K. Inal, A three dimensional (3D) thermo-elasto-viscoplastic constitutive model for FCC polycrystals, Int. J. Plast. 70 (2015) 166–190, http://dx.doi.org/10.1016/j.ijplas.2015.04.001.

[76] Y. Zhu, G. Kang, Q. Kan, O.T. Bruhns, Y. Liu, Thermo-mechanically coupled cyclic elasto-viscoplastic constitutive model of metals: Theory and application, Int. J. Plast. 79 (2016) 111–152, http://dx.doi.org/10.1016/j.ijplas.2015.12.005.

[77] J. Lei, Y. Ge, T. Liu, Z. Wei, Effects of heat treatment on the microstructure and mechanical properties of selective laser melting 316L stainless steel, Shock. Vib. 2021 (2021) 6547213, http://dx.doi.org/10.1155/2021/6547213.

[78] L. Thiercelin, F. Praud, F. Meraghni, E. Fleury, Thermally-activated hardening recovery in viscoplastic materials with kinematic hardening at high temperatures, Mech. Mater. 180 (2023) 104636, http://dx.doi.org/10.1016/j.mechmat.2023.104636.

[79] J. Trapp, A.M. Rubenchik, G. Guss, M.J. Matthews, In situ absorptivity measurements of metallic powders during laser powder-bed fusion additive manufacturing, Appl. Mater. Today 9 (2017) 341–349, http://dx.doi.org/10.1016/j.apmt.2017.08.006.

[80] J.K. Mackenzie, Second paper on statistics associated with the random disorientation of cubes, Biometrika 45 (1–2) (1958) 229–240, http://dx.doi.org/10.1093/biomet/45.1.2.229.

[81] S. Chen, G. Guillemot, C.-A. Gandin, Three-dimensional cellular automaton-finite element modeling of solidification grain structures for arc-welding processes, Acta Mater. 115 (2016) 448–467, http://dx.doi.org/10.1016/j.actamat.2016.05.011.

[82] A.S. Wu, D.W. Brown, M. Kumar, G.F. Gallegos, W.E. King, An experimental investigation into additive manufacturing-induced residual stresses in 316L stainless steel, Met. Mater. Trans. A 45 (13) (2014) 6260–6270, http://dx.doi.org/10.1007/s11661-014-2549-x.

[83] E.K. Karasz, J. Taylor, D.M. Autenrieth, P. Reu, K. Johnson, M. Melia, P. Noell, Measuring the residual stress and stress corrosion cracking susceptibility of additively manufactured 316L by ASTM G36-94, Corrosion 78 (1) (2021) 3–12, http://dx.doi.org/10.5006/3894.

[84] R.C. Laurence, D. Canelo-Yubero, E. Maawad, G.A. Faria, P. Staron, N. Schell, R.S. Ramadhan, S. Cabeza, A. Paecklar, T. Pirling, M. Sánchez Poncela, J.M. Martinez, M.F. Slim, T. Buslaps, P.J. Withers, R. Matthew, Determination of residual stress in additively manufactured 316L stainless steel benchmark parts through synchrotron X-ray diffraction and neutron diffraction, Strain 61 (2) (2025) e70005, http://dx.doi.org/10.1111/str.70005.

[85] I. Karaman, H. Sehitoglu, H. Maier, Y. Chumlyakov, Competing mechanisms and modeling of deformation in austenitic stainless steel single crystals with and without nitrogen, Acta Mater. 49 (19) (2001) 3919–3933, http://dx.doi.org/10.1016/S1359-6454(01)00296-8.

[86] M. Balbaa, M. Elbestawi, Multi-scale modeling of residual stresses evolution in laser powder bed fusion of inconel 625, J. Manuf. Mater. Process. 6 (1) (2022) 2, http://dx.doi.org/10.3390/jmmp6010002.

[87] J. Vang, Data science topics, 2019, URL https://datascience.oneoffcoder.com/migaussian.html.

[88] A. Strehl, J. Ghosh, Cluster ensembles—a knowledge reuse framework for combining multiple partitions, J. Mach. Learn. Res. 3 (Dec) (2002) 583–617.