# Thermo-mechanical simulation of L-PBF process at part-scale by coupling grain structure calculation and crystal viscoplasticity

Trung-Chien Vo, Zixuan Li, Michel Bellet, Charles-André Gandin, Yancheng Zhang

► To cite this version:

Trung-Chien Vo, Zixuan Li, Michel Bellet, Charles-André Gandin, Yancheng Zhang. Thermo-mechanical simulation of L-PBF process at part-scale by coupling grain structure calculation and crystal viscoplasticity. Computer Methods in Applied Mechanics and Engineering, 2026, 448, pp.118429. 10.1016/j.cma.2025.118429. hal-05342162

HAL Id: hal-05342162

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

Submitted on 1 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>

# Thermo-mechanical simulation of L-PBF process at part-scale by coupling grain structure calculation and crystal viscoplasticity

Trung-Chien Vo ID, Zixuan Li ID, Michel Bellet ID, Charles-André Gandin ID, Yancheng Zhang ID*

MINES Paris PSL University Centre de Mise en Forme des Matériaux (CEMEF), CNRS UMR 7635, CS 10297, 1 rue Claude Daunesse, Sophia Antipolis, 06904, France

## ARTICLE INFO

**Keywords:**

Laser powder bed fusion (L-PBF)

Cellular automaton (CA)

Grain structure

Crystal plasticity

Part-scale thermo-mechanical modeling

## ABSTRACT

An integrated numerical framework is proposed to perform a thermomechanical analysis of Laser Powder Bed Fusion (L-PBF) at the part scale, covering the entire construction process. The novelty lies in the fact that the thermomechanical analysis is driven by the evolving grain structure, also predicted at the part scale. The grain structure is first generated using a hybrid - Cellular Automaton (CA) method, which enables the generation of the grain structure on a CA grid at part-scale while incorporating the detailed scanning trajectories. For the thermo-mechanical simulation, a layer-by-layer thermal analysis includes the non-exposed powder, whereas the mechanical analysis only considers the substrate and the part under construction. A crystal plasticity model is employed for the nickel-based superalloy Inconel 718 (IN718), utilizing the CA grain structure referenced by Euler angles. A reduced grain approach is proposed, wherein Euler angles from the CA grid are projected onto the mechanical mesh. The mechanical mesh size is chosen appropriately to balance computational efficiency and accuracy. The calibration of the parameters of the crystal plasticity laws is conducted from room temperature up to 1100 °C using tests performed on homogenized IN718 material. Tensile tests on the representative volume elements (RVEs) of grain structures with different laser scanning trajectories are performed to study the mechanical response for each specified texture. Finally, thermo-mechanical simulations are applied to the L-PBF construction of a propeller with single-crystal and polycrystalline grain structures. At each time increment, the temperature field is transferred to the mechanical mesh, and grains are activated according to the part under construction. The influence of local grain texture on anisotropic stress distribution and distortion is analyzed and discussed.

## 1. Introduction

Additive manufacturing (AM) is the state of the art in generating complex-mechanical components, and Laser Powder Bed Fusion (L-PBF) is one of the promising processes in the AM process. L-PBF uses a fast-moving heat source to melt the material powder and form an object through a selective laser path. Customized products, therefore, are fabricated without requiring the traditional manufacturing tools or further manufacturing processes. In addition, L-PBF techniques can be effectively applied in scientific and medical fields because of their ability to produce prototypes quickly and cost-effectively. Therefore, this manufacturing approach has great potential and is being applied in many industries such as aviation, military, medical, and automation [1] (Table 1).

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: yancheng.zhang@minesparis.psl.eu (Y. Zhang).</footnote>
<table><caption>**Table 1**
Nomenclature.</caption><thead><tr><th colspan="4">Variables for leve-set method and thermal model</th></tr></thead><tbody><tr><td>t + Δt</td><td>Current time</td><td>Δt</td><td>Time increment</td></tr><tr><td>t'<sub>scan</sub></td><td>Scanning time</td><td>t'<sub>heat</sub></td><td>Heating time</td></tr><tr><td>t'<sub>cool</sub></td><td>Cooling time</td><td>t'<sub>dwell</sub></td><td>Dwell time</td></tr><tr><td>(x, y, z)</td><td>Global frame</td><td>(X, Y, Z)</td><td>Reference frame</td></tr><tr><td>V<sup>l</sup></td><td>Scanned volume</td><td>S<sup>l</sup></td><td>Area of a layer</td></tr><tr><td>Δz<sup>l</sup></td><td>Thickness of one layer</td><td>φ<sub>L</sub></td><td>Diameter of the laser beam</td></tr><tr><td>(ψ, φ<sub>m</sub>, φ, φ<sub>sl</sub>)</td><td>level-set functions</td><td>v<sub>L</sub></td><td>Laser scanning velocity</td></tr><tr><td>α<sub>h</sub></td><td>Artificial constant</td><td>h<sub>x</sub></td><td>Hatch spacing</td></tr><tr><td>P<sub>L</sub></td><td>Power of the laser beam</td><td>A</td><td>Absorption coefficient</td></tr><tr><td>Q̇<sub>L</sub></td><td>Volumetric heat source</td><td>ρ</td><td>Density</td></tr><tr><td>T</td><td>Temperature</td><td>T<sub>ext</sub></td><td>Air temperature</td></tr><tr><td>h</td><td>Enthalpy</td><td>n</td><td>Unit outward normal vector</td></tr><tr><td>κ</td><td>Thermal conductivity</td><td>h<sub>conv</sub></td><td>Heat exchange coefficient</td></tr><tr><td>{•}</td><td>Averaged operator</td><td>Ω</td><td>Computational domain</td></tr><tr><td>χ</td><td>Global properties</td><td>H</td><td>Heaviside function</td></tr><tr><td>D<sub>m</sub></td><td>Material domain</td><td>D<sub>g</sub></td><td>Gas domain</td></tr></tbody><thead><tr><th colspan="4">Variables for elastic-viscoplastic model</th></tr></thead><tbody><tr><td>ε̇</td><td>Total strain rate</td><td>ε̇<sup>el</sup></td><td>Elastic strain rate</td></tr><tr><td>ε̇<sup>vp</sup></td><td>Viscoplastic strain rate</td><td>ε̇<sup>th</sup></td><td>Thermal strain rate</td></tr><tr><td>T<sub>c</sub></td><td>Critical temperature</td><td>T<sub>v</sub></td><td>solidus temperature</td></tr><tr><td>T<sub>f</sub></td><td>Liquidus temperature</td><td>S</td><td>Compliant tensor</td></tr><tr><td>P</td><td>Symmetrized Schmid tensor</td><td>I</td><td>Identity tensor</td></tr><tr><td>α<sub>T</sub></td><td>One-dimensional</td><td>α</td><td>Slip system</td></tr><tr><td>γ̇<sup>a</sup></td><td>Thermal expansion coefficient</td><td>(φ<sub>1</sub>, φ, φ<sub>2</sub>)</td><td>Euler angles</td></tr><tr><td>m</td><td>Shear strain rate</td><td>γ̇<sub>0</sub></td><td>Reference shear strain rate</td></tr><tr><td>τ<sup>a</sup><sub>c,s</sub></td><td>Rate sensitivity</td><td>τ<sup>a</sup></td><td>Resolved shear stress</td></tr><tr><td>H<sub>αβ</sub></td><td>Thermally activable slip resistance</td><td>τ<sup>a</sup><sub>c,at</sub></td><td>Athermally activable slip resistance</td></tr><tr><td>h<sub>f</sub></td><td>The hardening modulus</td><td>q<sup>l</sup><sub>at</sub></td><td>The latent-hardening coefficients</td></tr><tr><td>h<sub>f</sub><sub>asympl</sub></td><td>Instantaneous hardening modulus</td><td>γ<sup>a</sup></td><td>Accumulated shear strain</td></tr><tr><td>ε̇</td><td>Asymptotic hardening modulus</td><td>h<sub>ratio</sub></td><td>Ratio between the initial and the final hardening</td></tr><tr><td>h<sub>exp</sub></td><td>Equivalent viscoplastic strain rate</td><td>s</td><td>Deviatoric stress</td></tr><tr><td>σ̇</td><td>Strength of the exponential part</td><td>K</td><td>Dynamic viscosity</td></tr><tr><td>Ω</td><td>Equivalent stress is calculated as</td><td>V<sup>c</sup></td><td>Volume of an element</td></tr><tr><td>E</td><td>Young's modulus</td><td>G</td><td>Shear modulus</td></tr><tr><td>E<sup>pre</sup></td><td>Predicted Young's modulus</td><td>E<sup>iso</sup></td><td>Isotropic Young's modulus</td></tr><tr><td>E<sub>ref</sub></td><td>Reference Young's modulus</td><td>E<sup>iso</sup><sub>ref</sub></td><td>Reference isotropic Young's modulus</td></tr><tr><td>(e<sub>h2</sub>, e<sub>h1</sub>, e<sub>h4</sub>)</td><td>Error between coarse meshes and reference mesh</td><td>(h<sub>1</sub>, h<sub>2</sub>, h<sub>3</sub>, h<sub>4</sub>)</td><td>Mesh types</td></tr></tbody><thead><tr><th colspan="4">Variables for printing process</th></tr></thead><tbody><tr><td>L<sub>x</sub>, L<sub>xy</sub></td><td>Scanning strategies</td><td>P<sub>laser</sub></td><td>Laser power</td></tr><tr><th colspan="4">Abbreviation</th></tr><tr><td>AM</td><td>Additive manufacturing</td><td>L-PBF</td><td>Laser powder bed fusion</td></tr><tr><td>CA</td><td>Cellular automaton</td><td>CP</td><td>Crystal plasticity</td></tr><tr><td>CAFE</td><td>Cellular automaton finite element</td><td>CPFE</td><td>Crystal plasticity finite element</td></tr><tr><td>DED</td><td>Directed energy deposition</td><td>FCC</td><td>Face-centered cubic</td></tr><tr><td>VPSC</td><td>Viscoplastic self-consistent</td><td>RVEs</td><td>Representative volume elements</td></tr><tr><td>EBSD</td><td>Electron backscatter diffraction</td><td>EVP</td><td>elastic-viscoplastic</td></tr></tbody></table>

To generate the microstructure during the solidification process, Gandin and Rappaz [2] coupled the cellular automaton (CA) method with finite element method (CAFE) to simulate grain growth in mushy zones of melt pools. Camus et al. [3] proposed a hybrid methodology to predict grain structure evolution with thermo-hydraulic simulations for L-PBF process at part-scale [4]. This approach predicts the melt pool shape and the local temperature field of the moving laser spot, coupling this temperature field with the CA model to construct full grain structure through scanning tracks. Recently, Kats et al. [5] have proposed a physics-informed machine learning algorithm to predict the grain structure in the directed energy deposition (DED) process, in which a neural network model is developed to identify the correlation between the local thermal features and their corresponding grain structure features, such as grain size.

Several macroscopic thermal simulations for constructing large parts in L-PBF process have been developed. Chiumenti et al. [6] introduced the layer-by-layer and multi-layered approaches to reduce the computational cost compared to the scanning track strategy. Later, Zhang et al. [7] proposed a series of methods to reduce the computational cost while maintaining the accuracy of thermal evolution predictions. They first introduced a strategy with the decomposition of a scanning layer by fractions, where the local temperature distribution and evolution can be simulated correctly. Another strategy proposed by Zhang et al. [8] is the super-layer
method with multi-shot method. This method demonstrated that the temperature evolution and distribution were appropriate and accurate compared to layer-based simulation.

To comprehensively investigate the mechanical behavior of polycrystalline alloy structures, the crystal plasticity finite element (CPFE) method is utilized. This method effectively accounts for crystallographic lattice orientations and grain morphology [9,10]. The accuracy of CPFE models in predicted mechanical behavior largely depends on the constitutive model. Early pioneering work by Peirce et al. [11,12] introduced the application of CPFE with a power-law model for face-centered cubic (FCC) metallic crystals. Their models numerically analyzed the deformation behavior of ductile single-crystals, explicitly accounting for lattice rotations. Subsequently, Cailletaud [13] proposed a power-law model by incorporating back stress effects for kinematic hardening. This enhanced formulation was later adopted by Anahid et al. [14] within the CPFE framework to develop a crystal plasticity-based crack nucleation criterion, enabling the prediction of the number of cycles leading to the crack nucleation in dwell fatigue experiments. To explicitly incorporate temperature effects to account for temperature-sensitive phenomena, Kocks et al. [15] introduced an Arrhenius-type power-law model based on activation energy. Leveraging this model, Cheng and Ghosh [16] addressed local thermal and athermal barriers to simulate heterogeneous deformation and crack nucleation driven by dislocation. The power law models capture deformation mechanisms driven by dislocation slip and are suitable for strain rates below $10^4 s^{-1}$.

Conversely, various homogenization methods have been integrated with CPFE to accurately predict the macroscopic mechanical response of the polycrystalline aggregates. One of the earliest approaches is the grain interaction model proposed by Taylor [17], which assumes a uniform macroscopic strain field identical to the local strain field within each grain. However, this method results in inconsistencies in the microscopic stress state at grain boundaries due to variations in grain orientations, violating stress equilibrium and leading to inaccuracies in the global mechanical behavior. To address these limitations, Molinari et al. [18] introduced the one-site viscoplastic self-consistent (VPSC) method to homogenize the mechanical response of polycrystalline material under large deformations. This approach builds on the concept proposed by Eshelby [19], embedding individual grains within a uniform matrix treated as a homogeneous medium. The VPSC method considers grains with elliptical shapes to model the interaction between grains and the macroscopic domain in elastic-viscoplastic regime. Later, Lebensohn and Tomé [20] developed a numerical formulation to simulate plastic deformation and texture evolution in polycrystalline alloys. This theory accommodates finite strain increments and infinitesimal elastic strains.

To characterize the mechanical response of crystals as a function of temperature, researchers employ crystallographic RVEs. For instance, Hu et al. [21] investigated the anisotropic mechanical properties of an aluminum alloy at elevated temperatures, incorporating the effects of damage. They utilized a crystal plasticity model coupled with thermal and damage modules to determine thermal material parameters. These parameters were calibrated through simulations on polycrystalline RVEs and validated using strain rate jump tests. Similarly, Li et al. [22] developed a thermo-mechanical framework based on an exponential crystal plasticity model. This approach accounted for thermal strains, temperature-dependent crystal behavior, and heat generation caused by dissipation during plastic slip. In contrast, Ghorbanpour et al. [23] studied the mechanical behavior of IN718 from room temperature to 550°C, which accounts for solid solution, precipitate shearing, grain size and shape contributions to initial slip resistance, its evolution with the dislocation density-based hardening law, and non-Schmid effects for the activation stress.

In the context of applying crystal plasticity to study mechanical behavior in AM, thermo-mechanical processes are typically examined at the meso-scale for single track simulations. Grilli et al. [24] developed a CPFE model coupled with element elimination and reactivation methods to simulate residual stress during the formation of grain structures in the L-PBF process for a single scan track. Their work focused on investigating the residual stress distribution along the laser scan direction. Subsequently, Hu et al. [25] validated this CPFE model by analyzing and comparing the lattice strain results with experimental data [26]. In a related work, Lindroos et al. [27] proposed a thermo-mechanical approach coupled with crystal plasticity to study residual stress states during phase transformations from austenite to martensite of H13 tool steel in single track.

This paper introduces a purely numerical study within a comprehensive framework for thermo-microstructure-mechanical simulations at the part-scale. The novelty of the work lies in the development of a macroscopic thermo-mechanical framework integrated with crystal plasticity during the simulation process, enabling the simulation of mechanical behavior with complex grain structures generated by the CAFE model. To address computational challenges associated with large-scale simulations, a reduced grain approach is proposed to simplify the grain representation in large-scale grain structures. In section 2, a specimen with a high number of randomly oriented grains is used to calibrate the parameters of a crystal plasticity model for IN718 material, based on a phenomenological elastic-viscoplastic model and material parameters from [28]. The impact of the reduced grain approach on simulation outcomes is also evaluated. In section 3, RVEs generated from L-PBF simulations with various scanning strategies [3] are employed to study the mechanical response. In section 4, the thermo-mechanical simulation is extended to a full grain structure by the proposed reduced grain method. A level-set and layer-by-layer methodology is incorporated within the thermal model to precisely control the macroscopic energy source deposition with a complex mesh. The stress response and total deformation of the part are analyzed on a simple mesh, highlighting the influence of the complex grain structure on the mechanical response, thus demonstrating the efficiency of the proposed thermo-mechanical framework. Conclusions are presented in section 5. All of the above simulations are implemented using an in-house code called Cimlib [3,4,7,8,29,30].

## 2. Thermo-mechanics based on grain structure at part-scale

The objective of the present study is to develop a thermo-mechanical simulation strategy tailored to capture the mechanical response of a part under construction, taking into account the formation of its grain structure. Considering the state of the art presented previously, it could be envisaged to couple CAFE grain structure calculation [2], which is based on the heat transfer resolution at
Track-scale thermo-hydraulic analysis for macro-scale grain growth

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Methodology for a weak coupling of the thermo-mechanical simulation with grain structure includes: track-scale thermo-hydraulic analysis for macro-scale grain growth using (a) CAFE model with hybrid methodology with (i) melt pool extracted at steady state, (ii) melt pool applied along the laser scanning path, (iii) grain structure; and thermo-mechanical simulation considering micro-structure at macro-scale consists of (b) reduced grain approach with (i) part-scale grain structure, (ii) mechanical mesh, (iii) mapped grain structure; (c) thermal simulation with (i) thermal mesh, (ii) thermal response; and (d) mechanical simulation with (i) grain information activated based on level-set ψ, (ii) mechanical response.</figcaption></figure>

part-scale [6–8], with the mechanical resolution at part-scale. However, this is not possible because the temperature fields generated by the part-scale thermal solvers in the vicinity of the construction front are not precise enough and generally lead to underestimated maximum temperatures in this region, due to the use of (multi-layer) flash heating methods [6–8]. In fact, this was the motivation for the development of the hybrid approach by Camus et al. [3,4] to predict the grain structure at the part-scale. Therefore, in the present study, the following approach is chosen. In a preliminary step, the grain structure is calculated at the part-scale using the hybrid approach. Once this is achieved, a thermo-mechanical simulation of the construction is conducted at the part-scale [31] in which, for each deposition of a new layer, the grain structure previously calculated by the hybrid approach is mapped onto the newly activated region around the construction front.

The present thermo-mechanical approach differs from [31] in that it uses a bi-mesh framework, consisting of separate thermal and mechanical meshes. Fig. 1 outlines the methodology for the strategy of thermo-mechanical simulation during the construction of a part, which includes the following steps. Details will be given in Sections 2.1 and 2.2.

* **Step 0 - Grain structure generation.** In a preliminary step to the thermo-mechanical analysis, the grain structure is generated at part-scale by using CAFE modeling with the hybrid approach developed by Camus et al. [3] (Fig. 1.a). This approach integrates track-scale thermo-hydraulic simulations to account for hydrodynamic effects, and allows extracting a consolidated track shape and temperature field at a steady-state (Fig. 1.a.i). These parameters are then applied along the effective laser scan path of the fabrication (Fig. 1.a.ii). The grain structure forms through the competition between nucleation and epitaxial growth of grains, governed by the temperature field. This process is modeled by coupling a CA grid, represented as a lattice of cubic cells, to simulate grain evolution. Using this methodology, the grain structure of either a RVE (Fig. 1.a.iii) or the whole part (Fig. 1.b.i) can be predicted.
* **Step 1 - Grain structure mapping.** The part-scale grain structure obtained by CAFE model with hybrid methodology (Fig. 1.b.i) is subsequently mapped onto the mechanical mesh shown in Fig. 1.b.ii, using the reduced grain approach. In this process of data transport, the grain information (represented by Euler angles), which is assigned to each finite element of the mechanical mesh, is the one of the CA cell containing the element barycenter (Fig. 1.b.iii). During part construction, the grain information on the mechanical mesh is activated according to the upward displacement of the construction front, which is defined by a level-set function (Fig. 1.d.i). Details on grain structure mapping will be given in Section 2.2.6.

* **Step 2 - Thermal analysis.** It is conducted using a thermal mesh that consists of the part under construction, the non-exposed powder, and the gas subdomain represented through the level-set method (Fig. 1.c.i). A non-linear thermal solver operated in this level-set framework allows computing the temperature distribution on the whole mesh. This temperature field is then transferred from the thermal mesh to the mechanical mesh for subsequent mechanical analysis. Details of the thermal analysis will be explained in Section 2.1.

* **Step 3 - Mechanical analysis.** Finally, the mechanical simulation (Fig. 1.d) is conducted on a mesh that encompasses the part under construction and the gas subdomain, but excluding the non-exposed powder (Fig. 1.b.ii). The mechanical solver integrates the temperature and grain structure data from the previous steps to simulate, in a level-set approach similar to the one in Step 2, the mechanical response of the part. Detailed information on mechanical analysis will be provided in Section 2.2.

For each new layer deposition, associated with the evolution of the level-set function, Steps 1 to 3 are successively operated, until the end of the construction, followed by the final cooling.

## 2.1. Heat transfer modeling at part-scale

In the thermal model, the computational domain $\Omega$ consists of a single unstructured mesh composed of tetrahedral elements. This domain $\Omega$ represents both material domain $D_m$ and gas domain $D_g$. Within the material domain, the dense part region comprises, on one hand, the part under construction and the substrate, and on the other hand the non-exposed powder. Thanks to the level-set framework, the interface between the material and gas domains can be efficiently captured, even for the complex geometry. The inclusion of non-exposed powder also eliminates the need to define a boundary condition along the lateral surface of the part. Additionally, the mesh is refined along the interfaces between these domains through level-set based mesh adaptation. This results in a smooth transition of the thermophysical properties from one domain to another, using the Heaviside function of the level-set [32]. Furthermore, thermal interface resistance [33,34] can be applied between two different domains, although this is not activated in the present study. A layer-by-layer approach is used to apply an appropriate heat source to each newly deposited layer. The thermal analysis provides the distribution of temperature $T$, temperature rate $\dot{T}$, and the level-set function which defines the construction progress, hereunder denoted $\psi$. These results are then transferred to the mechanical model at each time step. The level-set function $\psi$ will also define the interface between the grain structure in the part under construction and the gas zone within the mechanical mesh, ensuring accurate representation of both geometry evolution and material properties.

### 2.1.1. Level-set method and dynamic adaptive control

The interface between material and gas domains in the thermal mesh is described by the above mentioned level-set function $\psi$ (Fig. 1.c.i). The level-set is defined such that $\psi = 0$ at the interface, $\psi > 0$ in the gas domain, and $\psi < 0$ in the material domain. When a new layer of thickness $\Delta z^l$ is deposited, the function $\psi$ is updated by:

$$ \psi^{l+\Delta t} = \psi^l - \Delta z^l \quad (1) $$

where $t + \Delta t$ is the current time, and $\Delta t$ is the time increment. In the material zone, the constructed part and the powder are separated by an intermediate level-set function $\varphi_m$. The level-set $\varphi_m$ is defined such that $\varphi_m < 0$ indicates the powder zone, and $\varphi_m > 0$ indicates the part zone. To track the geometry of the constructed part (Fig. 1.c.i) separately from the gas domain and powder zone, a second level-set $\varphi$ is defined as a combination of $\psi$ and $\varphi_m$. Consequently, $\varphi > 0$ in the built part, and $\varphi < 0$ in the powder and gas domains, with:

$$ \varphi = \min(\varphi_m, -\psi) \quad (2) $$

To determine the material zone of the current layer being constructed, the level-set function for the construction front $\varphi_l$, is defined as:

$$ \varphi_l^{l+\Delta t} = \min(\varphi^{l+\Delta t}, \psi^{l+\Delta t} + \Delta z^l, -\psi^{l+\Delta t}) \quad (3) $$

where $\varphi_l$ is positive within the considered layer, and $\varphi_l$ is negative in the domain outside that layer. In addition, to efficiently integrate multiple physical domains within a single computational domain, the averaging operator $\{\cdot\}$ is introduced. This operator defines the averaging of global properties $\chi$ between material and gas domains by the smooth Heaviside function $H(\psi)$ [32] with:

$$ \{\chi\} = H(\chi)^{D_g} + (1 - H)(\chi)^{D_m} \quad (4) $$

where $H(\psi) > 0$ in the gas domain $D_g$ of the thermal mesh. Mesh refinement is applied at the interface between powder bed and gas, and between powder and part, and in the vicinity of the construction front to maintain the precision in the areas of interest. The size of anisotropic mesh used around interfaces is driven by a metric tensor $M$ [35], which is defined by:
$$ M = \begin{bmatrix} \frac{1}{h_{xj}^2} & 0 & 0 \\ 0 & \frac{1}{h_{yj}^2} & 0 \\ 0 & 0 & \frac{1}{h_{zl}^2} \end{bmatrix} \qquad (5) $$

where $h_{xj}$, $h_{yj}$, and $h_{zl}$ represent the mesh sizes along the three dimensions of a local reference frame that is tangent to the interfaces between domains in the thermal mesh.

### 2.1.2. Layer-by-layer method

The layer-by-layer method with flash heating, previously presented in [7,8], is the strategy adopted for part construction in this study. During the L-PBF process, the addition of a new layer is decomposed into three steps: heating, cooling, and dwell times. The dwell time $t'_{dwell}$ is a fixed parameter, depending on the printing machine, and corresponding to the time required to spread a new powder layer. The cooling time $t'_{cool}$ is calculated as follows:

$$ t'_{cool} = t'_{scan} - t'_{heat} \qquad (6) $$

where $t'_{scan}$ is the scanning time and $t'_{heat}$ is the heating time. $t'_{scan}$ is calculated from the area of the layer $S'$, the hatch spacing $h_s$, and the speed of the laser beam $v_L$ as:

$$ t'_{scan} = \frac{S'}{h_s v_L} \qquad (7) $$

The heat source to be assigned to each layer of thickness $\Delta z'$ is determined as follows. First, the heating time $t'_{heat}$ for the whole layer is expressed as:

$$ t'_{heat} = \alpha_h \frac{\phi_L}{v_L} \qquad (8) $$

where $\phi_L$ is the diameter of the laser beam and $v_L$ is the laser scanning velocity. Eq. (8) expresses that each printed layer is flash-heated as a whole, during a short time that is proportional to the exposure time $\frac{\phi_L}{v_L}$. The role of the multiplier $\alpha_h$ is to increase this heating time in order to avoid diffusive thermal shocks without altering the solution.

The value heat source $\dot{Q}'_L$ the results from the energy balance of each layer:

$$ AP_L t'_{scan} = \dot{Q}'_L V^l t'_{heat} \qquad (9) $$

where $P_L$ is the power of the laser beam, $A$ is the absorption coefficient, and $V^l = \Delta z' S^l$ is the scanned volume receiving the energy. Hence, $\dot{Q}'_L$ is also supplied as:

$$ \dot{Q}'_L = \frac{AP_L}{\alpha_h h_s \phi_L \Delta z'} \qquad (10) $$

Note that an increase in $t'_{heat}$, accordingly to $\alpha_h$, leads to a decrease in the heat source $\dot{Q}'_L$.

### 2.1.3. Thermal resolution

The thermal resolution is based on the equation for energy conservation, as follows:

$$ \frac{\partial \{\rho h\}}{\partial t} = \nabla \cdot (\{\kappa\} \nabla T) + \dot{Q}'_L \qquad (11) $$

where $T$ is the temperature, $\rho$ is the density, $h$ is the specific enthalpy, $\kappa$ is the thermal conductivity, and $\dot{Q}'_L$ is the volume heat source which is applied during the heating time $t_{heat}$ of each layer. The convective heat exchange is considered at the bottom surface of the substrate by:

$$ -\{\kappa\} \nabla T \cdot \mathbf{n} = h_c (T - T_{ext}) \qquad (12) $$

where $h_c$ is the heat exchange coefficient, $T_{ext}$ is the air temperature, and $\mathbf{n}$ is the unit outward normal vector.

To account for the effect of argon gas flow within the chamber, the temperature after the cooling step is reinitialized to room temperature $T_0$ using Heaviside function $H$, as described in [7], where:

$$ T^{l+\Delta t} = T_0 H(\psi') + T^l(1 - H(\psi')) \qquad (13) $$

The adiabatic boundary condition is applied to the top surface of the entire domain. The simulation domain being large enough, the adiabatic boundary condition is also applied to the lateral surfaces. No specific condition is applied along the powder bed/gas, powder bed/part, and the part/gas interfaces. As explained in more detail in [7,8], the Galerkin weak form of Eq. (11) is spatially discretized with linear P1 tetrahedral elements, with an implicit time scheme. This results in a set of non-linear equations which is solved for the nodal temperature values using an iterative Newton-Raphson method. At each Newton-Raphson iteration, corrections to the nodal temperature estimates are computed by an iterative solver. A detailed description can be found in [40]. The solution scheme also includes the Cond_Split technique dedicated to the treatment of diffusional thermal shocks [41]. Due to mesh refinement within each layer during thermal simulation, the element size along the energy deposition zone is relatively small to avoid thermal shocks. The properties of IN718 material, powder bed, and argon gas for thermal simulation are presented in Table 2.
<caption>**Table 2**</caption>

Properties of IN718 alloy (dense state and powder bed), and of argon gas for thermal simulation.

<table><thead><tr><th>Characteristics</th><th>Variables</th><th>Units</th><th>Values</th><th>Ref</th></tr></thead><tbody><tr><td>Material solidus temperature</td><td>T<sub>s</sub></td><td>[°C]</td><td>1170</td><td>[36]</td></tr><tr><td>Material liquidus temperature</td><td>T<sub>l</sub></td><td>[°C]</td><td>1336</td><td>[36]</td></tr><tr><td>Room temperature</td><td>T<sub>0</sub></td><td>[°C]</td><td>20</td><td></td></tr><tr><td>Solid density of material</td><td>ρ<sub>s</sub></td><td>[kg.m<sup>-3</sup>]</td><td>8190</td><td>[36]</td></tr><tr><td>Liquid density of material</td><td>ρ<sub>l</sub></td><td>[kg.m<sup>-3</sup>]</td><td>7300</td><td></td></tr><tr><td>Porosity of powder bed</td><td>φ<sub>p</sub></td><td>[-]</td><td>0.5</td><td></td></tr><tr><td>Density of powder bed</td><td>ρ<sub>p</sub></td><td>[kg.m<sup>-3</sup>]</td><td>ρ<sub>s</sub>(1 - φ<sub>p</sub>)</td><td></td></tr><tr><td>Density of gas</td><td>ρ<sub>g</sub></td><td>[kg.m<sup>-3</sup>]</td><td>1.3</td><td>[37]</td></tr><tr><td>Heat capacity of solid material and powder bed</td><td>C<sub>p,s</sub></td><td>[J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>512</td><td>[37]</td></tr><tr><td>Heat capacity of liquid material</td><td>C<sub>p,l</sub></td><td>[J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>720</td><td>[37]</td></tr><tr><td>Heat capacity of gas</td><td>C<sub>p,g</sub></td><td>[J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>520</td><td>[37]</td></tr><tr><td>Absorption coefficient</td><td>A</td><td>[-]</td><td>0.3</td><td>[38]</td></tr><tr><td>Heat exchange coefficient</td><td>h<sub>c</sub></td><td>[W.m<sup>-2</sup>.K<sup>-1</sup>]</td><td>150</td><td></td></tr><tr><td>Thermal conductivity of material</td><td>κ<sub>m</sub></td><td>[W.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>(Appendix A)</td><td>[36]</td></tr><tr><td>- powder bed</td><td>κ<sub>p</sub></td><td>[W.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>(Appendix A)</td><td>[39]</td></tr><tr><td>- gas</td><td>κ<sub>g</sub></td><td>[W.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>0.025</td><td>[37]</td></tr></tbody></table>

## 2.2. Mechanical modeling at part-scale

Neglecting gravity and inertia effects, the mechanical equilibrium equation is expressed as follows:

$$ \nabla \cdot \{\sigma\} = 0 \qquad (14) $$

where {$\sigma$} denotes the stress tensor, expressed in the global reference frame, and mixed in the level-set sense according to Eq. (4), while based on the element. In the material the stress results from the crystal elastic-viscoplastic behavior law, as explained in the next sub-section, while in the gas domain, the stress results from an arbitrarily adapted Newtonian behavior law, not presented here. The weak form of Eq. (14) is considered in a mixed pressure/velocity framework. It is discretized with stable bubble-type linear tetrahedral elements P1+/P1, with an implicit time scheme, resulting in a set of non-linear equations that are solved for nodal velocities and pressures by a Newton-Raphson algorithm. At each Newton-Raphson iteration, corrections to the nodal estimates are computed by a preconditioned iterative solver. Detailed information on FE discretization and solution procedures can be found in [42].

### 2.2.1. Elastic-viscoplastic model

Several constitutive laws of crystal plasticity (CP) were investigated to study grain structures [11,12,43,44]. Among these, the widely used power law model proposed by Hutchinson [44] remains a foundational model. Under the assumption of small strains, the constitutive equations of the elastic-viscoplastic (EVP) model are the followings. They are expressed in each finite element in a local reference frame aligned with the <100> direction, according to the grain structure mapping defined in Step 1 of Section 2.

$$ \dot{\epsilon} = \dot{\epsilon}^{el} + \dot{\epsilon}^{vp} + \dot{\epsilon}^{th} \quad \text{where} \quad T \le T_c \qquad (15) $$

where $\dot{\epsilon}$, $\dot{\epsilon}^{el}$, $\dot{\epsilon}^{vp}$, $\dot{\epsilon}^{th}$ are the total, elastic, viscoplastic, and thermal strain rates, respectively, and $T_c$ is a critical temperature for the consideration of EVP law. The expressions of the different components are:

$$ \begin{aligned} \dot{\epsilon}^{el} &= \mathbf{S}^{el} \dot{\sigma} \\ \dot{\epsilon}^{vp} &= \sum_{\alpha=1}^{n} \dot{\gamma}^{a} \mathbf{P}^{\alpha} \\ \dot{\epsilon}^{th} &= \alpha_{T} \mathbf{T} \mathbf{I} \end{aligned} \qquad (16) $$

where $\mathbf{S}^{el}$ is the 4th order elastic compliance tensor, which characterizes the anisotropic elastic properties of IN718 at room temperature (20 °C), as reported by Kumara et al. [45]. $\dot{\sigma}$ is the stress rate tensor, $n$ is the total number of slip systems, $\mathbf{I}$ is the identity tensor, $\alpha_T$ is the one-dimensional thermal expansion coefficient, and $\mathbf{P}^{\alpha}$ is the symmetrized Schmid tensor of the slip system $\alpha$. $\mathbf{P}^{\alpha}$ is expressed by:

$$ \mathbf{P}^{\alpha} = \frac{1}{2}(\mathbf{S}^{\alpha} + (\mathbf{S}^{\alpha})^{T}) \qquad (17) $$

where $\mathbf{S}^{\alpha}$ is the Schmid tensor, generated by the normal $n^{\alpha}$ of the slip plane and the slip direction $m^{\alpha}$ as:

$$ \mathbf{S}^{\alpha} = \mathbf{m}^{\alpha} \otimes \mathbf{n}^{\alpha} \qquad (18) $$
The power-law model [44] is applied for material calibration, using the following formula:

$$ \dot{\gamma}^{\alpha} = \dot{\gamma}_0 \left( \frac{|\tau^{\alpha}|}{\tau_c^{\alpha}} \right)^{\frac{1}{m}} \mathrm{sign}(\tau^{\alpha}) \qquad (19) $$

For a slip system $\alpha$, $\dot{\gamma}^\alpha$ represents the shear rate, $\dot{\gamma}_0$ denotes the reference shear rate, $m$ is the strain rate sensitivity, $\tau^\alpha$ indicates the resolved shear stress obtained by $\tau^\alpha = \sigma : P^\alpha$, and $\tau_c^\alpha$ is the critical resolved shear stress. The hardening rule encompasses the thermally activable slip resistance part $\tau_{c,*}^\alpha$ and the athermally activable slip resistance $\tau_{c,a}^\alpha$. This formulation can be expressed as:

$$ \tau_c^\alpha = \tau_{c,*}^\alpha(T) + \tau_{c,a}^\alpha \qquad (20) $$

where the expression of $\tau_{c,a}^\alpha$ is taken from [46,47]:

$$ \tau_{c,a}^{\alpha} = \int_{0}^{t} \sum_{\beta=1}^{n} H^{\alpha\beta} |\dot{\gamma}^{\beta}| dt \qquad (21) $$

where the coefficients $H^{\alpha\beta}$ represent the hardening modulii, which follow:

$$ H^{\alpha\beta} = h_{\gamma}(q^{\mathrm{lat}} + (1 - q^{\mathrm{lat}})\delta^{\alpha\beta}) \qquad (22) $$

where $q^{\mathrm{lat}}$ stands for the hardening coefficients that determine the degree of latent-hardening. In which, $\alpha = \beta$ provides self-hardening, and $q^{\mathrm{lat}} = 1.4$ is for the latent-hardening when $\alpha \neq \beta$, $\delta^{\alpha\beta}$ is Kronecker's delta, and $h_\gamma$ is the instantaneous hardening modulus. $h_\gamma$ depends on the accumulated shear strain $\gamma^\alpha$ following:

$$ h_{\gamma} = h_{\mathrm{asympl}}(1 + (h_{\mathrm{ratio}} - 1)e^{-h_{\mathrm{exp}}|\gamma^{\alpha}|}) \qquad (23) $$

where $h_{\mathrm{asympl}}$ is the asymptotic hardening modulus, $h_{\mathrm{ratio}}$ is the ratio between the initial and the asymptotic hardening, and $h_{\mathrm{exp}}$ is the strength of the exponential part.

### 2.2.2. Viscoplastic law

When the temperature is greater than the critical temperature $T_c$ but remains below the liquidus temperature $T_l$, the material behavior of the mushy zone is assumed to be governed by the Norton-Hoff model. The decomposition of the strain rate tensor is now:

$$ \dot{\mathbf{e}} = \dot{\mathbf{e}}^{\mathrm{vp}} + \dot{\mathbf{e}}^{\mathrm{th}} \quad \text{where} \quad T_c < T < T_l \qquad (24) $$

where the viscoplastic constitutive law for this model is as follows:

$$ \dot{\mathbf{e}}^{\mathrm{vp}} = \frac{1}{2K} \left( \sqrt{3} \dot{\mathbf{e}} \right)^{1-m_{\mathrm{vp}}} \mathbf{s} \qquad (25) $$

where $\mathbf{s} = \sigma - \frac{1}{3} \mathrm{tr}(\sigma)\mathbf{I}$ represents the deviatoric stress, $\dot{\mathbf{e}}$ is the generalized viscoplastic strain rate defined as $\dot{\mathbf{e}} = \sqrt{\frac{2}{3}} \dot{\mathbf{e}}^{\mathrm{vp}} : \dot{\mathbf{e}}^{\mathrm{vp}}$. The parameters $K$ and $m_{\mathrm{vp}}$ refer to the viscoplastic consistency and viscoplastic strain rate sensitivity, respectively, both of which are illustrated in Fig. A.2. Consequently, the von Mises stress $\bar{\sigma} = \sqrt{\frac{2}{3}} \mathbf{s} : \mathbf{s}$ is expressed as:

$$ \bar{\sigma} = K \left( \sqrt{3} \right)^{m_{vp}+1} \dot{\mathbf{e}}^{m_{vp}} \qquad (26) $$

When $T > T_l$, the Newtonian behavior is applied, with $m_{vp} = 1$ in the previous equation, leading to:

$$ \dot{\mathbf{e}}^{\mathrm{vp}} = \frac{1}{2K} \mathbf{s} \qquad (27) $$

$K$ appearing then as the dynamic viscosity. In the mechanical resolution, elements having their mean temperature greater than $T_c$ have a viscoplastic or Newtonian behavior. They are included like other elements in the global assembly and resolution process. Details can be found in [48,49]. It should be noted that the gas domain is also resolved by assuming the Newtonian behavior, with an arbitrary $K$ value of 0.1 Pa.s. Based on modules presented, discretization analysis for thermo-mechanical model is presented in Appendices B and C.

### 2.2.3. Calibration of a temperature-dependent crystal plasticity model for IN718

The single-crystal properties appearing in the power-law (Eq. (19)) and in the hardening law (Eq. (20)) are calibrated for the IN718 alloy. For the calibration process, a reference set of stress-strain curves is generated across a range of temperatures. These curves (labeled IS) are numerically produced using the phenomenological elastic-viscoplastic model described in [31], with model parameters previously identified in [28] based on experimental data. In contrast, another set of homogenized stress-strain curves (labeled PC) is obtained by simulating uniaxial tensile tests on a polycrystalline cube using the previously presented crystal elastic-viscoplasticity model. The single-crystal mechanical properties of the IN718 alloy are calibrated by comparing the PC curves with the IS reference curves. The polycrystalline cube with dimensions 0.1 × 0.1 × 0.1 mm³ is subjected to a temperature range from 20 °C to 1100 °C, with the bottom surface fixed and a strain rate $\dot{\epsilon}_{\mathrm{imp}} = 10^{-3} \mathrm{s}^{-1}$ applied along the z direction. The polycrystalline cube has 3213 grains distributed randomly to ensure isotropy. The volume-averaged strain of the polycrystal cube is calculated through time
<table><caption>**Table 3**
Elastic constants of single-crystal IN718.</caption><thead><tr><th>S<sub>11</sub> [GPa<sup>-1</sup>]</th><th>S<sub>12</sub> [GPa<sup>-1</sup>]</th><th>S<sub>44</sub> [GPa<sup>-1</sup>]</th><th>Ref</th></tr></thead><tbody><tr><td>0.007</td><td>-0.003</td><td>0.010</td><td>[51]</td></tr><tr><td>0.008</td><td>-0.004</td><td>0.009</td><td>[50,52]</td></tr><tr><td>0.008</td><td>-0.003</td><td>0.009</td><td>[53]</td></tr><tr><td>Average</td><td>0.008</td><td>-0.003</td><td>0.009</td><td>[45]</td></tr></tbody></table>

<table><caption>**Table 4**
The modified single-crystal-elastic constants.</caption><thead><tr><th></th><th>S<sub>44</sub> [GPa<sup>-1</sup>]</th><th>G [GPa]</th></tr></thead><tbody><tr><td>Calculated [45]</td><td>0.009</td><td>110.8</td></tr><tr><td>Modified</td><td>0.010</td><td>104.0</td></tr></tbody></table>

integration of the volume averaged strain rate. The formulation for the volume-averaged stress and strain components in the tensile direction are expressed as follows:

$$ 
\begin{cases}
\langle \sigma_{zz} \rangle = \frac{1}{\sum_e V_e} \sum_e \sigma_{zz_e} V_e \\
\langle \epsilon_{zz} \rangle = \int_0^t \left( \frac{1}{\sum_e V_e} \sum_e \dot{\epsilon}_{zz_e} V_e \right) dt
\end{cases}
\qquad (28) $$

where $\Omega$ is the computational domain, and $V^e$ is the volume of a finite element $e$.

### 2.2.4. Calibration for elastic constants of the IN718 single-crystal

The elastic parameters of IN718 single-crystal were calibrated by Kumara et al. [45] and later utilized by Chen et al. [50] in crystal plasticity model for AM. These parameters are summarized in Table 3. The Young's modulus $E$, the Poisson's ratio $\nu$, and the shear modulus $G$ are derived from the anisotropic compliance matrix aligned with $\langle 100 \rangle$ crystallographic direction by the formula:

$$ S_{\langle 001 \rangle}^{el} = \begin{bmatrix} \frac{1}{E_{\bar{\nu}}} & \frac{-\bar{\nu}}{E_{\bar{\nu}}} & \frac{-\bar{\nu}}{E_{\bar{\nu}}} & 0 & 0 & 0 \\ \frac{-\bar{\nu}}{E_{\bar{\nu}}} & \frac{1}{E_{\bar{\nu}}} & \frac{-\bar{\nu}}{E_{\bar{\nu}}} & 0 & 0 & 0 \\ \frac{-\bar{\nu}}{E_{\bar{\nu}}} & \frac{-\bar{\nu}}{E_{\bar{\nu}}} & \frac{1}{E_{\bar{\nu}}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G} \end{bmatrix} \qquad (29) $$

$E$, $\nu$, and $G$ are then used in the polycrystalline tensile test to compare with the tensile test using material properties in [28]. As a result, the shear modulus $G$, which is equal to $1/S_{44}$, is modified by approximately 7% to fit the polycrystalline tensile tests with the isotropic tensile test. The computed parameters, therefore, are shown in Table 4.

Data of Young's modulus for the isotropic IN718 in [28] is provided in a range from room temperature to solidus one. Based on these data, Young's modulus values below solidus temperature are adopted by a scaling operation:

$$ E^{sca}(T) = E^{iso}(T) \frac{E_{ref}}{E_{ref}^{iso}} \qquad (30) $$

where $E^{sca}(T)$, $E^{iso}(T)$ are respectively the scaled, isotropic Young's modulus at a certain temperature $T$, and $E_{ref}^{iso}$, $E_{ref}$ are respectively the isotropic, reference Young's modulus at 20°C. A similar equation is used to define the shear modulus $G(T)$ at temperature $T$.

### 2.2.5. Calibration for IN718 parameters of power-law model

Using IN718 data in [28], the parameters of the elastic-viscoplastic model are calibrated to accurately describe the material behavior. The reference plastic shear rate $\dot{\gamma}_0$ is set to $10^{-4} \text{ s}^{-1}$, the exponential parameter $h_{exp}$ is assigned a value of 2000, and the ratio $h_{ratio}$ is fixed at 35 across all temperatures. The values of $h_{exp}$ and $h_{ratio}$ are adjusted by trial and error to identify the mechanical response along crystallographic directions $\langle 100 \rangle$, $\langle 110 \rangle$, and $\langle 111 \rangle$. While $\dot{\gamma}_0$, $h_{exp}$, and $h_{ratio}$ remain constant, the rate sensitivity $m$, the thermally activable slip resistance $\tau_{c,*}^\alpha$, and the asymptotic hardening coefficient $h_{asymp}$ are calibrated based on the experimental data, as illustrated in Fig. 2.

According to the kinetic model proposed by Kocks [54], the rate sensitivity coefficient $m$ is temperature-dependent, with $m$ equal to 0.04 at room temperature and 0.45 at 1100°C. For numerical convergence, $m$ is adjusted: it is set to a small value (0.05) at room temperature and increases linearly from 500°C to 1100°C to account for the significant variation in yield stress across this temperature range, ensuring consistency with the data from [28]. Regarding the decreasing functions $\tau_{c,*}^\alpha(T)$ and $h_{asymp}(T)$, $\tau_{c,*}^\alpha$ is calibrated to quantify the stress magnitude as a function of strain, while $h_{asymp}$ is tuned to model the hardening behavior. The IN718 material properties used in the crystal viscoplastic model for the mechanical simulation are summarized in Table 5.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Calibrated values versus temperature:(a) m, (b) $\tau_{c,s}^a$, and (c) $h_{asymp}$.</figcaption></figure>

<caption>**Table 5**</caption>

Properties of IN718 alloy used in the crystal elastic-viscoplastic model and in the viscoplastic model for the mechanical simulation.

<table><thead><tr><th>Characteristics</th><th>Variables</th><th>Units</th><th>Values</th></tr></thead><tbody><tr><td>Critical temperature</td><td>T<sub>c</sub></td><td>[°C]</td><td>1170</td></tr><tr><td>Viscoplastic law over T<sub>c</sub></td><td></td><td></td><td></td></tr><tr><td>Consistency</td><td>K</td><td>[kPa.s]</td><td>Fig. A.2.a</td></tr><tr><td>Strain rate sensitivity</td><td>m<sub>cp</sub></td><td>[-]</td><td>Fig. A.2.b</td></tr><tr><td>Crystal elastic-viscoplastic law below T<sub>c</sub></td><td></td><td></td><td></td></tr><tr><td>Elasticity parameters</td><td>S<sub>11</sub>, S<sub>12</sub>, S<sub>44</sub></td><td>[GPa<sup>-1</sup>]</td><td>Table 3</td></tr><tr><td></td><td>G, E<sup>eca</sup></td><td>[GPa]</td><td>Table 4 & Eq. (30)</td></tr><tr><td>Thermal expansion coefficient</td><td>α<sub>T</sub></td><td>[K<sup>-1</sup>]</td><td>Fig. A.2.c</td></tr><tr><td>Thermally activable slip threshold</td><td>τ<sub>c,s</sub><sup>a</sup></td><td>[MPa]</td><td>Fig. 2.b</td></tr><tr><td>Strain rate sensitivity</td><td>m</td><td>[-]</td><td>Fig. 2.a</td></tr><tr><td>Asymptotic hardening coefficient</td><td>h<sub>asymp</sub></td><td>[MPa]</td><td>Fig. 2.c</td></tr><tr><td>Latent hardening coefficient</td><td>q<sup>lat</sup></td><td>[-]</td><td>1.4</td></tr><tr><td>Ratio coefficient</td><td>h<sub>ratio</sub></td><td>[-]</td><td>35</td></tr><tr><td>Strength of exponential coefficient</td><td>h<sub>exp</sub></td><td>[-]</td><td>2000</td></tr></tbody></table>

In Fig. 3.a, the calculated-averaged stress-strain curves for the random isotropic equiaxed polycrystal (PC), mentioned in Section 2.2.3, are manually calibrated and adjusted to match those of the isotropic specimen (IS) across a range of temperatures. To guarantee the precision, the fitting error is controlled within 2% with a mesh size $h = 1.5$ µm. Due to the inherent limitations of the power-law model with the low hardening behavior, the fitting error increases to 6% at 1100 °C. However, this level of error is considered acceptable, as the stress magnitude at 1100 °C is significantly lower compared to other temperature ranges. The results demonstrate that the crystal plasticity model for the polycrystal effectively captures the stress-strain behavior observed in the isotropic specimen.

In the same tensile test configuration at room temperature, three different orientations of the single-crystal with directions (100), (110), and (111) respectively aligned along the tensile direction are examined at room temperature (20 °C) (Fig. 3.b). Notably, beyond
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) The averaged stress-strain curves for the random isotropic equiaxed polycrystal (PC) compared to the stress-strain curves for the isotropic specimen (IS); (b) Tensile tests on the single-crystal cube at room temperature. Direction *z* is the tensile direction.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. (a) Transfer of Euler angles from CA grid to finite element of mechanical mesh in the reduced grain approach; (b) Transfer of a RVE grain structure to four finite element meshes for sensitivity test.</figcaption></figure>

strain $\epsilon = 0.005$, the stress in traction aligned with $\langle 100 \rangle$ is higher than that of $\langle 110 \rangle$, while the traction along $\langle 111 \rangle$ exhibits the highest stress. The observed stress-strain behavior of these single-crystals is consistent with the findings in [22,55].

### 2.2.6. Reduced grain approach

As shown in Fig. 1, the modeling framework integrates several computational meshes exchanging information on regions of the whole domain, and tailored for specific calculations: the CA grid for calculating grain structure, the thermal mesh for thermal simulation, and the mechanical mesh for mechanical simulation. Crystal plasticity calculations are conducted on the mechanical mesh. Ideally, the CA grid would be used directly for mechanical computations to retain the complete grain structure information. However, this approach is computationally prohibitive due to the large number of CA cells. As a sustainable alternative, the mechanical mesh with a mesh size larger than the grid size is employed, enabling reduced computational cost.

In the present study, a simple approach is used, named reduced grain approach, with consists in assigning to each element of the mechanical mesh the crystal orientation (Euler angles) from the CA cell containing the element barycenter, as illustrated in Fig. 4.a.
<figure><img src="image_7.png" /><figcaption>Fig. 5. The deviation between the reference mesh and the other meshes.</figcaption></figure>

<caption>**Table 6**</caption>

Relative errors $e_{h_i}$ associated with three different mesh sizes, as compared to reference mesh size $h_1$.

<table><thead><tr><th rowspan="2">Temperature</th><th colspan="3">Error[%]</th></tr><tr><th>e<sub>h<sub>2</sub></sub></th><th>e<sub>h<sub>3</sub></sub></th><th>e<sub>h<sub>4</sub></sub></th></tr></thead><tbody><tr><td>20 °C</td><td>1.40</td><td>3.47</td><td>8.52</td></tr><tr><td>800 °C</td><td>0.73</td><td>1.73</td><td>4.95</td></tr><tr><td>900 °C</td><td>0.62</td><td>1.61</td><td>4.83</td></tr><tr><td>1100 °C</td><td>0.89</td><td>1.97</td><td>5.61</td></tr></tbody></table>

Due to the larger mesh size, compared to the cell size, there is a corresponding loss of grain information from the CA grid in the direct assignment procedure. This grain information loss becomes more significant as the mesh size increases, as illustrated in Fig. 4.b.

In order to evaluate the influence of mesh size in the context of the reduced grain approach, four elementary cubic domains (edge size 100 µm, about 3000 grains, and averaged grain size 8.6 µm) are meshed differently and submitted to a virtual tensile test, performed in similar strain rate conditions and at four different temperatures: 20, 800, 900, and 1100 °C. The first cube, which serves as a reference, has a structured mesh with a mesh size which is the same as CA cells: $h_1 = 1.5$ µm, each CA cell being divided into six tetrahedra. The three other cubes have a non-structured mesh with mesh size respectively 3, 5, and 10 µm. The distribution of the first Euler angle is shown in Fig. 4.b for each of the four meshes.

Each grain of regular size is represented by a cluster of finite elements, enabling adjustments to the grain boundary and the volume fraction of the grain within the RVE by modifying the element size. The mechanical meshes ensure a random distribution of grains in the cube with a large number of grains. As a result, the cube exhibits an isotropic behavior for the four mechanical meshes. The different average stress-strain curves are plotted in Fig. 5. The errors between the averaged stress-strain curves (Table 6) obtained with the reference mesh $h_1$ and the three other meshes across the different temperatures are calculated as follows:

$$e_{h_i} = \int_0^t \frac{|\langle\sigma_{zz_{h_1}}\rangle - \langle\sigma_{zz_{h_i}}\rangle|}{\langle\sigma_{zz_{h_1}}\rangle} dt \quad i = 2, 3, 4 \qquad (31)$$

Table 6 presents a synthetic recap of the errors, showing that the error $e_{h_4}$ is considerably larger than $e_{h_2}$ and $e_{h_3}$. Table 6 demonstrates that the error $e_{h_1}$ is larger at lower temperature, reaching 8.5% at room temperature. While adjusting the mesh size can reduce the error, another significant portion of the error originates from the inherent properties of the material.

In Fig. 4.b, each grain in meshes consists of multiple elements. As the mesh size increases from the finest mesh $h_1$ to the coarsest mesh $h_4$, the grain boundaries in the reference mesh $h_1$ become noticeably less accurate in the coarser mesh $h_4$. Concurrently, the volume fraction of each grain changes significantly due to the expansion or contraction of neighboring grains within the RVE domain. As shown in Fig. 5, the stress-strain response derived from mesh $h_4$ deviates considerably from that of the finer meshes across the
<table><caption>**Table 7**
Computational time for tensile tests on the random isotropic equiaxed polycrystal with different mesh types.</caption><thead><tr><th>Mesh</th><th>Size [µm]</th><th>Number of elements</th><th>Computational time [s]</th></tr></thead><tbody><tr><td>h<sub>1</sub></td><td>1.5</td><td>1,804,578</td><td>4923</td></tr><tr><td>h<sub>2</sub></td><td>3</td><td>366,067</td><td>969</td></tr><tr><td>h<sub>3</sub></td><td>5</td><td>80,682</td><td>267</td></tr><tr><td>h<sub>4</sub></td><td>10</td><td>4,722</td><td>36</td></tr></tbody></table>

same temperature range. This suggests that variations in grain boundary and grain volume fraction, both induced by mesh coarsening, are primary contributors to discrepancies in stress-strain behavior. Conversely, when grain boundaries and volume fractions remain preserved across different mesh resolutions, the stress-strain response remains relatively consistent.

This stability condition is naturally satisfied in grain structures formed by AM, which often produce large, columnar grains with pronounced crystallographic texture. In such cases, mesh refinement, where the mesh element size is significantly smaller than the grain size, does not substantially alter the grain boundary or volume fractions. As a result, the mechanical response of these large columnar grains is much less sensitive to mesh resolution.

To demonstrate this, traction tests are conducted along three directions (x, y and z) for a grain structure with marked textures, referred to as RVE #1. This RVE, generated using the CAFE method from L-PBF process simulation, will be described in detail in section 3, where its full microstructural characteristics reflect AM grain structures. The deviations of stress-strain curves between mesh $h_3 = 5$ µm and mesh $h_5 = 30$ µm for RVE #1 (see Fig. D.3 in the Appendix D) under the three traction directions (x, y, and z) (Fig. D.4.b-d in the Appendix D) are 1.20%, 2.30%, and 1.89% respectively. These are significantly smaller than the deviation ($e_{h_4} - e_{h_3}$) = 5.05% (Table 6 at room temperature) between mesh $h_3 = 5$ µm and mesh $h_4 = 10$ µm observed in the polycrystalline cube with random grain orientations (Fig. D.4.a in the Appendix D), despite $h_5$ being three times larger than $h_4$. These results confirm that grain structures with large columnar grains, such as RVE #1, where the grain size exceeds the mesh element size, exhibit minimal sensitivity to mesh refinement.

Therefore, CP simulations benefit greatly from the application of the reduced grain approach for grain structures generated in AM, composed of large columnar grains.

Table 7 reports the computational times for different mesh sizes employed in the tensile tests on the PC illustrated in Fig. 5. The simulations were performed on workstation equipped with two ADM Roma 7502 processors, each featuring 32 cores of 2.5GHz. The number of elements in mesh $h_4$ is reduced by a factor of 380 compared to mesh $h_1$, resulting in a significant decrease in computational time from $h_1$ to $h_4$. Consequently, the reduced grain method offers significant advantages, including a substantial reduction in computational time, preservation of computational accuracy, and applicability to part-scale models.

### 3. Verification of the crystal plasticity model for RVEs

In this section, the sensitivity of the coupled methodology, combining crystal plasticity and the reduced grain approach, is evaluated using different grain structures that arise from additive manufacturing. This serves as a preliminary step toward applying the methodology to large-scale L-PBF simulations, as described in Section 4.

RVEs reflecting diverse grain structures are employed to investigate the sensitivity of mechanical response to microstructural morphology. Each RVE measures 300 × 300 × 300 µm³ and is extracted from a larger volume of 600 × 600 × 1100 µm³. Within this larger domain, the grain structure formed via epitaxial growth is generated using the hybrid approach presented in [3]. The simulations incorporate varying sets of L-PBF processing parameters, including scanning strategy, laser spot velocity, laser power, and hatch spacing, all of which are critical in determining the formation of columnar grains. These columnar grains, in turn, significantly influence the unique mechanical behavior in each RVE. Notably, the grain structures within the RVEs exhibit no periodicity along their boundaries; as a result, periodic boundary conditions are not applied.

Four cubic RVEs of same dimensions are utilized, as shown in Fig. 6. The first RVE, denoted as RVE #1 (Fig. 6.a), is generated using back and forth scan tracks along the x direction ($L_x$), with laser parameters specified in Table 8. This RVE exhibits a grain orientation with (110) along the y and z directions, and (100) along the x direction. RVE #2 and RVE #3 (Fig. 6.b-c) follow a similar scanning strategy as RVE #1 but differ by the hatch spacing and the laser power as detailed in Table 8. In contrast, RVE #4 is generated using scan tracks deposited in both the x and y directions subsequently ($L_{xy}$) on each layer, resulting in a more randomized grain orientation within the structure (Fig. 6.d). All RVE structures are generated on a CA grid with a resolution $h_{grid} = 1.5$ µm. This resolution enables to capture intricate microstructural details.

Virtual tensile tests are studied on each of the specified RVEs along the x, y, and z directions to investigate their anisotropic mechanical behavior (Fig. 7). Using the reduced grain approach, the grain orientations of RVEs, represented by Euler angles within the CA grid, are transferred to the mechanical mesh with a resolution $h_4 = 10$ µm. These tests are conducted at room temperature under a strain rate of $\dot{\epsilon} = 10^{-3}$ s⁻¹. In RVE #1, a concentration of randomly oriented grains appears at the center of each scan track, creating thin layers of random grain orientations between columnar grains along the y direction. The spacing between these layers corresponds to the hatch spacing, as shown in Fig. 6.a.

The stress response of RVE #1 observed in tensile tests conducted along the x, y, and z directions, as illustrated in (Fig. 7.a-c) and (Fig. 8), reveals significant anisotropic mechanical behavior. Specifically, the stress response ($\sigma_{xx}$) for a traction along x differs
<figure><img src="image_8.png" /><figcaption>Fig. 6. The inverse pole figures and volumetric pole figures of (a) RVE #1, (b) RVE #2, (c) RVE #3, (d) RVE #4.</figcaption></figure>

<caption>Table 8</caption>

Process parameters of the four cases: type of scan path, hatch spacing $h_s$, laser power $P_L$, and laser beam velocity $v_L$.

<table><thead><tr><th>RVE case</th><th>Scan</th><th>$h_s$ [µm]</th><th>$P_L$ [W]</th><th>$v_L$ [mm s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>1</td><td>$L_x$</td><td>76</td><td>120</td><td>400</td></tr><tr><td>2</td><td>$L_x$</td><td>50</td><td>120</td><td>400</td></tr><tr><td>3</td><td>$L_x$</td><td>76</td><td>80</td><td>400</td></tr><tr><td>4</td><td>$L_{xy}$</td><td>76</td><td>120</td><td>400</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7. RVE #1. Distribution of tensile stress component at strain $\epsilon = 0.08$ in tensile tests applied along x, y, and z directions: (a) $\sigma_{xx}$, (b) $\sigma_{yy}$, (c) $\sigma_{zz}$.</figcaption></figure>

markedly from the stress responses $\langle\sigma_{yy}\rangle$ and $\langle\sigma_{zz}\rangle$ in traction along y and z, respectively. In contrast, the stress responses $\langle\sigma_{yy}\rangle$ and $\langle\sigma_{zz}\rangle$ exhibit similar behavior. This shows that the primary grain orientations in RVE #1 contribute to an anisotropic mechanical behavior.

To further investigate the influence of grain orientation on the mechanical response, a single-crystal cube, labeled as SXCube, is generated to compare with RVE #1. In SXCube, a uniform single grain orientation with $\langle 110 \rangle$ along the y and z directions and $\langle 100 \rangle$ along the x direction are considered. Notably, under a tensile test along the z direction, $\sigma_{zz}$ of SXCube is identical to $\sigma_{yy}$ under traction along the y direction. Consequently, the stress response of RVE #1 is expected to resemble that of SXCube. The results in
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Comparison between the pole figures and average axial stress response in RVE #1 under tensile tests along x, y and z directions (continuous curves) respectively with the corresponding response of the single-crystal cube (SXCube dashed lines).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Distribution of stress component $\sigma_{xx}$ at strain $\epsilon = 0.08$ in different RVEs submitted to a tensile test along x direction (a) #1, (b) #2, (c) #3, and (d) #4.</figcaption></figure>

Fig. 8 indicate that the stress response $\langle\sigma_{xx}\rangle$ of RVE #1 aligns closely with $\sigma_{xx}$ of SXCube under tensile test along the x direction. However, notable differences are observed in the stress response $\langle\sigma_{yy}\rangle$ of RVE #1 compared to $\sigma_{yy}$ of the single-crystal cube in the tensile test along the y direction. Specifically, $\langle\sigma_{yy}\rangle$ is higher than $\langle\sigma_{xx}\rangle$ in RVE #1, whereas $\sigma_{yy}$ is lower than $\sigma_{xx}$ in SXCube. This shows that the stress response of RVE #1 is influenced by contributions from grains not aligned with the dominant grain orientation, and this influence is non-negligible.

RVEs #2, #3, and #4 (Fig. 6.b-d) are subjected to the same set of tensile tests in the three directions and at room temperature as RVE #1 to evaluate their mechanical response. As illustrated in (Fig. 9.a-c), the average tensile stresses of RVEs #2 and #3 exhibit similar trends compared to RVE #1. Additionally, the inverse pole figures and volumetric pole figures for RVEs #2 and #3 (Fig. 6.b-c) reveal a dominant grain orientation $\langle 100 \rangle$ along the x direction, although with increased noise levels compared to RVE #1. Consequently, RVE #2 and #3 exhibit a similar stress response of component $\sigma_{xx}$ to RVE #1.

Furthermore, due to the use of similar scanning strategies, the stress responses of $\langle\sigma_{xx}\rangle$, $\langle\sigma_{yy}\rangle$, and $\langle\sigma_{zz}\rangle$ under tensile tests along x, y and z directions, respectively, are consistent among RVEs #1, #2, and #3, as shown in (Fig. 10.a-c). In contrast, RVE #4, depicted in Fig. 6.d, has a less marked texture. Consequently, the stress-strain curves in Fig. 10.d for tensile tests in three different x, y, and z directions display nearly identical behaviors, indicating the isotropic mechanical behavior.

In the tensile test along x direction for RVE #4, as shown in Fig. 9.d, no significant stress accumulation in layers is observed, in contrast to RVEs #1, #2, and #3 (Fig. 9.a-c). In Fig. 9, the overall stress level in RVE #4 is noticeably higher than in RVEs #1, #2, #3. This difference arises from the dependence of grain orientation on the tensile direction. Specifically, the dominant grain orientation $\langle 100 \rangle$ in RVEs #1, #2, #3 is aligned with traction along x direction, facilitating slip system activation at lower stress levels. In contrast, the more random grain orientation in RVE #4 requires higher stress to activate all slip systems. These observations are consistent with
<figure><img src="image_12.png" /><figcaption>Fig. 10. Stress-strain curves for tensile tests in directions x, y, and z for (a) RVE #1, (b) RVE #2, (c) RVE #3, and (d) RVE #4.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. The propeller with laser scanning only in x direction for all of layers: (a) scanning strategy, (b) fully constructed part.</figcaption></figure>

findings from the previous section, further emphasizing the critical influence of the polycrystalline structure and the volume fraction of grains on the mechanical behavior.

## 4. Application to the L-PBF manufacturing of a small turbine propeller

The previously developed methodology is applied to the numerical simulation of the L-PBF construction of a small turbine propeller, as illustrated in Fig. 11, designed within a domain of 5.2 × 5.2 × 2.35 mm³. The propeller, made of IN718, is constructed layer by layer, with a scanning laser beam moving back and forth along the x direction for each layer. The laser trajectory, automatically
<table><caption>**Table 9**
Process parameters for the propeller manufacturing.</caption><thead><tr><th>Characteristics</th><th>Variables</th><th>Units</th><th>Values</th></tr></thead><tbody><tr><td>Laser power</td><td>P<sub>L</sub></td><td>[W]</td><td>120</td></tr><tr><td>Velocity</td><td>v<sub>L</sub></td><td>[mm s<sup>-1</sup>]</td><td>400</td></tr><tr><td>Hatch spacing</td><td>h<sub>s</sub></td><td>[µm]</td><td>60–82</td></tr><tr><td>Layer thickness</td><td>Δz<sup>l</sup></td><td>[µm]</td><td>25</td></tr><tr><td>Total layers</td><td>n</td><td>[-]</td><td>90</td></tr><tr><td>Interlayer dwell time</td><td>t<sup>l</sup><sub>dwell</sub></td><td>[s]</td><td>11.12</td></tr><tr><td>Diameter of the laser beam</td><td>φ<sub>L</sub></td><td>[µm]</td><td>70</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Computational domains of (a) the thermal mesh (TM) and (b) the mechanical mesh (M2).</figcaption></figure>

generated by the open source software Slic3r [56], accounts for the complex geometry of the propeller. As a result, the scan tracks of successive layers are not vertically aligned (Fig. 11.a), differing from the fabrication strategy employed for the RVEs in Section 3. Consequently, the dominant grain orientation in the propeller deviates from that observed in the RVEs. The process parameters are detailed in Table 9.

The grain structure formed at the end of the process has been simulated by the CAFE model using the hybrid method, on a CA grid with a cell size of $h_{grid} = 4.5$ µm. Compared to the previous RVEs, the significantly larger size of the propeller necessitates an increase in the cell size, resulting in the loss of some microstructural details. However, the overall texture of the grain structure, including the formation of columnar grains, is effectively captured as illustrated further. The grain structure information is subsequently transported from the CA grid to the mechanical mesh with a coarser resolution, having a mesh size of $h_{mech} = 50$ µm, ten times larger than that of the CA cell.

The thermo-mechanical simulation of the construction process is conducted as presented in Section 2, with the computational domains of the thermal and mechanical meshes shown in Fig. 12.a-b. The simulation of each layer encompasses different stages: heating and material deposition using the flash heating method [8], followed by an interlayer dwell time corresponding to the spreading of a new powder layer in the real process. The temperature peaks during the heating phase and then falls to a minimum after the dwell time. To track the temperature evolution throughout the process, a virtual sensor is positioned within the substrate at the coordinates (0.35, 0, -0.025). This sensor records temperature variations over the 90 printed layers. As shown in Fig. 13.a, the area of the horizontal cross-section A (at a height of 0.1 mm) is larger than that of section B (at a height of 1.45 mm) and section C (at a height of 2 mm). Consequently, the temperature evolution, depicted in Fig. 13.b, indicates that the temperature decreases significantly after the dwell time as the cross-sectional area of the printed layers reduces. As mentioned in Section 2, the thermo-mechanical framework integrates the computed temperature field from the thermal mesh (Fig. 13.c), which is mapped onto the mechanical mesh for thermo-mechanical simulation (Fig. 13.d).

## 4.1. Influence of grain texture on stress distribution

For the analysis and comparison, several distinct crystal structures are generated, each characterized by specific grain morphologies and volumetric pole figures, as shown in Fig. 14. The grain structure $G_{110}$ consists of a single-crystal configuration with (110) orientation along x and z directions and a (100) orientation along the y direction (Fig. 14.a). This is achieved by assigning a fixed grain orientation, represented by the Euler angles $(\phi_1, \phi, \phi_2) = (180^\circ, 90^\circ, 45^\circ)$, to every element of the mechanical mesh. The grain
<figure><img src="image_15.png" /><figcaption>Fig. 13. (a) Considered cross-sections of the propeller; (b) Temperature evolution during the construction extracted from the thermal sensor; (c) Temperature field in thermal mesh (encompassing the non-scanned powder and the gas); (d) Temperature field transported to the mechanical mesh.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. The grain morphology and volumetric pole figures of grain structures (a) $G_{110}$, (b) $G_{rand}$, and (c) $G_{L_x}$.</figcaption></figure>

structure $G_{rand}$ consists of a random distribution of grain orientations, defined by random Euler angles with the following conditions ($\phi_1 \in [-180^\circ, 180^\circ]$, $\phi \in \pm\text{acos}([-1, 1]) * 180^\circ/\pi$, $\phi_2 \in [-180^\circ, 180^\circ]$), applied to the mechanical mesh (Fig. 14.b). In contrast, the grain structure $G_{L_x}$ consists of the CAFE-predicted grain structure, with laser scanning applied in back and forth mode along the x direction for all layers [4] (Fig. 14.c). It is expected that these distinct grain distributions significantly influence the mechanical response of the propeller during its construction process, as well as the stress and distortion distribution at the end of the fabrication. This paper only presents the final state (stress and distortion).
<figure><img src="image_17.png" /><figcaption>Fig. 15. The mechanical response obtained for different grain structures: $G_{110}$, $G_{rand}$, and $G_{L_x}$, with (topline) $\sigma_{xx}$, (middle) $\sigma_{yy}$, and (bottom) $\sigma_{zz}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. (a) The location of lines 1 and 2 in the propeller, (b) Stress components $\sigma_{xx_{L1}}^{G_{110}}$, $\sigma_{yy_{L1}}^{G_{110}}$, and $\sigma_{zz_{L1}}^{G_{110}}$ extracted from line 1, and (c) The stress components $\sigma_{xx_{L1}}^{G_{110}}$, $\sigma_{xx_{L1}}^{G_{rand}}$, $\sigma_{xx_{L1}}^{L_x}$, and $\sigma_{xx_{L1}}^{ISO}$ with respect to grain structures $G_{110}$, $G_{rand}$, $G_{L_x}$, and isotropic structure $ISO$ extracted from line 1.</figcaption></figure>

In a geometry such as this propeller, isotropic mechanical behavior would induce stress components $\sigma_{xx}$ in cross-section E similar to $\sigma_{yy}$ in cross-section D (Fig. 13.a). Similarly, the stress component $\sigma_{yy}$ in cross-section E would mirror the stress component $\sigma_{xx}$ in cross-section D. Looking at Fig. 15, second column ($G_{rand}$), the two first lines $\sigma_{xx}$ and $\sigma_{yy}$, the expected result for a random equiaxed grain structure is effectively obtained. This is not the case for the grain structure $G_{110}$, which clearly shows an anisotropic response. For the grain structure $G_{L_x}$, the mirroring effect is somewhat retrieved, indicating that the x-y anisotropy is perceptible again, but
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** (a) Symmetrical axes x and y on the cross-section A for grain structure $G_{rand}$, (b) The stress distribution comparison of cases $G_{rand}$ and $G_{110}$, (c) The stress distribution comparison of cases $G_{L_x}$ and $G_{110}$.</figcaption></figure>

not as marked as for $G_{110}$. Moreover, for the $G_{L_x}$ grain structure, the effect of grain morphology leads to stress concentration aligned with the interfaces of the columnar grains.

To clarify the impact of the anisotropy arising from the microstructure, the stress and displacement components along two identified vertical lines 1 and 2 (Fig. 16.a) are analyzed and compared. Fig. 16.b illustrates the stress components along line 1 associated with the single-crystal grain structure $G_{110}$, showing variations in $\sigma_{xxL_1}^{G_{110}}$, $\sigma_{yyL_1}^{G_{110}}$, and $\sigma_{zzL_1}^{G_{110}}$ from the base to the top of the propeller. The higher stress levels are localised in the basement for the three components. However, $\sigma_{xxL_1}^{G_{110}}$ and $\sigma_{yyL_1}^{G_{110}}$ are clearly positive (i.e. tensile), while $\sigma_{zzL_1}^{G_{110}}$ is negative (i.e. compressive). By further comparing the stress component $\sigma_{xx}$ along line 1 with $G_{110}$, $G_{rand}$, and $G_{L_x}$, Fig. 16.c demonstrates a similar global trend, but with significant differences between the grain structures, particularly in the basement and in the lower region of the part, from height -0.1 to 1.5 mm. In the upper region of the part, from height 1.5 to 2.25 mm, regardless of the grain structure, the stress profiles are similar. Additionally, the stress profile obtained with the phenomenological isotropic elastic-viscoplastic approach described in [31], with experimental data from [28], is presented with the label *Iso*. This profile exhibits a similar global trend to the others, albeit with some differences. Note in particular that $G_{110}$ and *Iso* provide smoother profiles than $G_{rand}$ and $G_{L_x}$ calculated with grain-dependent constitutive equations. This reveals the impact of grain structure on stress distribution.

Due to the significant differences in the stress response at the base of the propeller, the stress components in the horizontal cross-section A, defined in Fig. 13.a, are analyzed. Assuming an isotropic mechanical behavior, the profile of $\sigma_{xx}$ along the x axis is expected to be similar to that of $\sigma_{yy}$ along the y axis in this cross-section. Fig. 17.a demonstrates that for the random equiaxed microstructure $G_{rand}$, this is effectively verified. This is confirmed, more quantitatively, by the plot of the two profiles $\sigma_{xx}^{G_{rand}}$ and $\sigma_{yy}^{G_{rand}}$ given in Fig. 17.b which are very close to each other. In contrast, the single-crystal orientation $G_{110}$ exhibits two very different stress profiles of $\sigma_{xx}^{G_{110}}$ and $\sigma_{yy}^{G_{110}}$, highlighting the pronounced anisotropy. Additionally, Fig. 17.c illustrates that the effective grain structure $G_{L_x}$ presents a mixed grain orientation, combining characteristics of the $G_{110}$ configuration with perturbation associated with grain orientation, as presented in the pole figure of Fig. 14.c.

## 4.2. Influence of grain texture on distortion

Using the data extracted along vertical lines 1 and 2 (Fig. 16.a), the displacement analysis in Fig. 18 highlights the relationship between grain structure and mechanical response. To examine the anisotropic behavior of the uniform $G_{110}$ grain structure, the displacement component $u_y$ along line 1, denoted as $u_{yL_1}^{G_{110}}$, is compared to the displacement component $u_x$ along line 2, denoted as $u_{xL_2}^{G_{110}}$. This comparison mirrors the one done previously for stress analysis.
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** (a) Displacement components calculated when considering grain structures ($G_{110}$, $G_{rand}$, $G_{L_x}$), and when using the phenomenological isotropic elastic-viscoplastic law (Iso), plotted along lines 1 and 2, (b) Inverse pole figures and pole figures of determined cross-sections of the grain structure $G_{L_x}$.</figcaption></figure>

As shown in Fig. 18.a, the two vertical lines cross through horizontal cross-sections of three different types:

- Type 1 cross-sections (TS1): massive cross-sections in the basement

- Type 2 cross-sections (TS2): cross-type cross-sections in the propeller body

- Type 3 cross-sections (TS3): reduced annular cross-sections at the top of the part.

In TS1 (from height -0.1 to 0.5 mm) and TS2 (from 0.5 to 1.3 mm), the two displacements $u_{yL1}^{G_{110}}$ and $u_{xL2}^{G_{110}}$ have the same sign, but the amplitude of $u_{xL2}^{G_{110}}$ is higher than $u_{yL1}^{G_{110}}$, whereas in TS3 (from 1.3 to 2.25 mm), they are nearly identical. This response, observed
with the single-crystal $G_{110}$ grain structure, demonstrates that the impact of the anisotropic behavior also depends on the geometry of the part.

The $G_{110}$ grain structure is then compared to $G_{L_x}$. In TS1, the largest gap between $u_{xL_2}^{G_{110}}$ and $u_{yL_1}^{G_{110}}$, or between $u_{xL_2}^{G_{L_x}}$ and $u_{yL_1}^{G_{L_x}}$, occurs at height 0.1 mm. At this height, the $G_{L_x}$ grain structure in the cross-section, as shown in the pole figure (Fig. 18.b), exhibits a random grain orientation with a slight concentration of grains clustered around the (110) orientation. Despite this slight concentration, $u_{xL_2}^{G_{L_x}}$ remains higher than $u_{yL_1}^{G_{L_x}}$. In TS2, at height 1 mm, where the gap between $u_{xL_2}^{G_{L_x}}$ and $u_{yL_1}^{G_{L_x}}$ is the largest, the cross-section texture reveals a strong dominant grain orientation (110). This explains why $u_{xL_2}^{L_x}$ closely matches $u_{yL_1}^{G_{110}}$, while $u_{yL_1}^{L_x}$ closely matches $u_{yL_1}^{G_{110}}$. In TS3, at height 2 mm, the cross-section displays a strong dominant grain orientation with additional noise in the pole figure. However, $u_{xL_2}^{G_{L_x}}$ and $u_{yL_1}^{G_{L_x}}$ are nearly identical due to the geometric effect mentioned earlier. Additionally, in Fig. 18.a, the random distribution of grains in the structure $G_{rand}$ results in nearly identical displacement components, $u_{xL_2}^{G_{rand}}$ and $u_{yL_1}^{G_{rand}}$. This result is consistent with the displacement responses of the isotropic structure Iso, where $u_{xL_2}^{G_{Iso}}$ is identical to $u_{yL_1}^{G_{Iso}}$. Furthermore, the displacement behavior of the $G_{rand}$ structure is very similar to that of the Iso structure. These observations indicate the isotropic mechanical response obtained with the $G_{rand}$ grain structure and confirm the correctness of the calibration of the crystal elastic-viscoplastic model based on experimental data. Finally, the crystal elastic-viscoplastic model required 24 hours of CPU time, whereas the isotropic structure model was completed in 16 hours.

In summary, the mechanical analysis conducted on the propeller demonstrates that both grain structure and part geometry have a significant impact on the mechanical response of the component.

## 5. Conclusion

A thermo-mechanical modeling framework has been developed using a purely numerical approach to simulate grain structure-based mechanical behavior in L-PBF at part-scale. The framework consists of the following key components:

- A bi-mesh approach was implemented to improve the computational efficiency. The thermal simulation utilizes an extended mesh encompassing the part under construction, the non-exposed powder, and the gas domain. In contrast, the mechanical simulation adopts a simplified mesh that excludes the non-exposed powder.

- A layer-by-layer method involving flash heating was applied to deposit material and energy onto the thermal mesh. The separately pre-calculated grain structure and the temperature field solved in the thermal mesh were mapped onto the mechanical mesh. Mechanical simulations incorporated a crystal plasticity model to account for the potential anisotropic behavior of grain structures.

- A reduced grain approach was implemented to transfer grain information from cellular automaton cells to mechanical finite element meshes. This enables efficient thermo-mechanical simulations at the part-scale.

- For the thermo-mechanical simulation, the parameters of the crystal plasticity law for the nickel-based superalloy IN718 were calibrated based on the phenomenological elastic-visoplastic model and the parameters from [28], covering a temperature range from room temperature to solidus.

Mesh sensitivity analyses were conducted on random grain RVEs across various mesh sizes, demonstrating the relationship between solution error and mesh refinement when applying the reduced grain approach. Furthermore, mesh refinement has a smaller influence on error in case of columnar grain structure or marked texture.

Mechanical simulations for RVEs and thermo-mechanical simulations of a small-size turbine propeller, whose grain structures were investigated in [3], were performed using the optimized mesh. Although the results have not been experimentally validated and cannot be used to draw allow quantitative conclusions, the simulations clearly demonstrate the influence of microstructure on mechanical response. Notably, variations in residual stresses and total displacements were observed between the as-predicted polycrystalline structure and the arbitrary single-crystal structures. This emphasises the critical influence of grain orientation on mechanical behavior at part-scale.

Overall, the present study showcases a path for numerical modeling of the relationship between process parameters, grain structure formed, and stress/strain distributions, for a part under construction by L-PBF.

For the perspectives, this study reveals a certain sensitivity of the mechanical response to the mesh size. This sensitivity, which is inherent in the reduced grain approach, could be mitigated by alternative homogenization procedures that leverage grain information more effectively. Therefore, methods for calculating the average stresses of multiple grains can be examined to reduce errors in the mechanical behaviors arising from the transport of grain information in the mechanical mesh. In addition, further refinement of the modeling assumptions and enabling of detailed experimental validation are tasks that remain for future studies.

## CRediT authorship contribution statement

**Trung-Chien Vo:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Data curation; **Zixuan Li:** Software, Methodology; **Michel Bellet:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Conceptualization; **Charles-André Gandin:** Writing – review & editing, Validation, Supervision, Investigation, Conceptualization; **Yancheng Zhang:** Writing – review & editing, Validation, Supervision, Software, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.
**Data availability**

Data will be made available on request.

**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.

**Acknowledgment**

The authors would like to 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. Material parameters for IN718**

The temperature dependence of the thermal conductivity in the dense state $κ_m$ and the powder state $\langle κ_p \rangle$ are depicted in Fig. A.1 following [57].

The temperature dependence of the consistency coefficient $K$, the sensitivity $m_{vp}$, and The thermal expansion $\alpha_T$ are given in Fig. A.2.

<figure><img src="image_21.png" /><figcaption>Fig. A.1. Temperature dependence of the thermal conductivity for IN718 material and powder.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. A.2. (a) Consistency $K$, (b) strain rate sensitivity $m_{vp}$, and (c) linear thermal expansion coefficient $\alpha_T$ as function of temperature.</figcaption></figure>
# Appendix B. Thermal resolution

## B.1. Weak form and FE discretization

The weak form of energy conservation Eq. (11) with test function $w^*$ in thermal domain $\Omega$ with boundary $\delta\Omega = \Gamma_D \cup \Gamma_N$ can be expressed by:

$$ \int_{\Omega} w^* \frac{\partial H}{\partial t} d\Omega + \int_{\Omega} \nabla w^* \cdot (\{\kappa\} \nabla T) d\Omega - \int_{\Omega} w^* \dot{Q}_L^t d\Omega + \int_{\Gamma_N} w^* h_c (T - T_{ext}) d\Gamma = 0 \quad (B.1) $$

where $\{\rho h\}$ is noted as $H$, with $w^* = 0$ on $\Gamma_D$.

Using Galerkin FE method under full-implicit Euler-backward scheme [40,58], Eq. (B.1) can be solved using the following discrete equation for temperature $T$ at time $t + \Delta t$.

$$ R(T) = M(H^{t+\Delta t} - H^t) + KT - F - Q = 0 \quad (B.2) $$

As mentioned in Table 2, this is a non-linear system due to the temperature-dependent thermal properties. The resolution of Eq. (B.2) can be obtained by using the Newton-Raphson method with the tangent matrix:

$$ \frac{\partial R}{\partial T} = M \frac{\partial H}{\partial T} + K \quad (B.3) $$

where $\frac{\partial H}{\partial T}$ is given numerically, as demonstrated in the references [40,59]. The non-linear Eq. (B.2), solved at the iteration $(v + 1)$, becomes:

$$ \frac{\partial R^{(v)}}{\partial T} = -R(T^{(v)}) \quad (B.4) $$

This allows for the calculation of a new estimate: $T^{(v+1)} = T^{(v)} + \delta T^{(v+1)}$

## B.2. Cond_Split method

To avoid thermal shock, the Cond_Split method has been adopted based on previous works [41]. The governing discrete equation (Eq. (B.2)) by splitting the diffusion terms as:

$$ R(T) = M(H^{t+\Delta t} - H^t) + K^*T - F - Q - S = 0 \quad (B.5) $$

where

$$ S = (K^* - K)T \quad (B.6) $$

with $K^*_j = \int_{\Omega} \{\kappa^*\} \nabla N_i \nabla N_j d\Omega + \int_{\Gamma_N} N_i N_j h_c d\Gamma$, where $N_i$ and $N_j$ are shape functions of node $i$ and $j$ ($i, j \in [1, 4]$).

By assuming $S$ to be a known source term, an augmented thermal conductivity $\kappa^*$ can be defined by satisfying the diffusion depth condition with the time increment $\Delta t_{ts} = h_c^2 \rho C_p / \kappa$:

$$ \kappa^* = \begin{cases} \kappa & \text{if } \Delta t_{ts} \le \Delta t \\ \kappa \frac{\Delta t_{ts}}{\Delta t} & \text{if } \Delta t_{ts} > \Delta t \end{cases} \quad (B.7) $$

where $h_c$ is the size of the local element. As the source term is considered as an explicit approximation by using the known $T^t$ from the previous time step, Thus:

$$ S \approx (K^* - K)T^t \quad (B.8) $$

It should be noted that this Cond_Split method implied by Eq. (B.8) is local and temporary.

# Appendix C. Mechanical resolution

## C.1. Weak form and FE discretization

The weak form of the equilibrium equation (Eq. (14)) with test function $v^*$, combined with the mass conservation equation using the test function $q^*$ in mechanical domain $\hat{\Omega}$ with boundary $\delta\hat{\Omega} = \hat{\Gamma}_D \cup \hat{\Gamma}_N$ is described as:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\hat{\Omega}} \dot{\epsilon}(v^*) : \{s\}(v) d\hat{\Omega} - \int_{\hat{\Omega}} p \nabla \cdot v^* d\hat{\Omega} - \int_{\hat{\Gamma}_N} v^* \cdot \tau_{imp} d\hat{\Gamma} = 0 \quad (a) \\ \\ \displaystyle \int_{\hat{\Omega}} q^* \left(-\nabla \cdot v - \frac{\dot{p}}{\chi} + 3\alpha_T \dot{T}\right) d\hat{\Omega} = 0 \quad (b) \end{array} \right. \quad (C.1) $$

where $\dot{\epsilon}$ is the total strain rate in the global frame, $\tau_{imp}$ is the imported stress applied to the surface $\hat{\Gamma}_N$ (here, the imposed stress is ignored), and $s = \sigma + pI$, where $s$ is the global deviatoric stress and $p$ is the hydrostatic pressure with $v^* = 0$ on $\hat{\Gamma}_D$. Eq. (C.1) is
then discretized on a finite element mesh $\Omega_h$, composed of linear tetrahedral elements, with a $P_1$ + $P_1$ formulation [42,60]. It is just reminded here that this formulation basically consists in adding three extra degrees of freedom for the velocity at the center of each element to ensure the Brezzi-Babuška condition [61,62]:

$$v_h = v_h^L + b_h \qquad (C.2)$$

Hence, the discretized velocity field $v_h$ is divided into two parts: the first part, $v_h^L$, is linear over the element and satisfies the velocity-imposed boundary conditions. The second part, $b_h$, is of “bubble-type”: linear in each of the four sub-tetrahedra constituting the element and taking zero values along the entire boundary of the element. The finite element form of Eq. (C.1) is given as follows:

$$\left\{
\begin{array}{l}
\displaystyle \int_{\hat{\Omega}_h} \dot{\epsilon}(v_h^{L*}) : \{s\}(v_h^L) d\hat{\Omega} - \int_{\hat{\Omega}} p_h \nabla \cdot v_h^{L*} d\hat{\Omega} - \int_{\hat{\Gamma}_N} v_h^{L*} \cdot \tau_{imp} d\hat{\Gamma} = 0 \quad (a) \\[2ex]
\displaystyle \int_{\hat{\Omega}_h} \dot{\epsilon}(b_h^*) : s(b_h) d\hat{\Omega} - \int_{\hat{\Omega}} p_h \nabla \cdot b_h^* d\hat{\Omega} = 0 \quad (b) \\[2ex]
\displaystyle \int_{\hat{\Omega}_h} q_h^* \left( -\nabla \cdot v - \frac{\dot{p}_h}{\chi} + 3\alpha_T \hat{T} \right) d\hat{\Omega} = 0 \quad (c)
\end{array}
\right.
\qquad (C.3)$$

To simplify the appendix, the details of the resolution of Eq. (C.3) are not given here. It is worth noting that the “bubble” extra-unknowns $b_h$ can be eliminated because they are internal to each finite element. This explains that the second equation in Eq. (C.3) (b) can in fact be injected into the two others, yielding only $(v_{L,h}, p_h)$ as principal unknowns. Details on “bubble elimination” can be found in [58,63] and references therein. The resulting set of equations is then solved with the Newton-Raphson iterative method.

### C.2. Local resolution for each element

It should be noted that different constitutive equations cannot be used concurrently in a given finite element, i.e. TVP at an integration point and TEVP at another one. This is inherent to the $P_1$ + $P_1$ character of elements.

Thus, the state of each element should first be determined, and then the stress $\{s\}$ and $\{\partial s/\partial \dot{e}\}$ should be resolved locally and fed back into Eq. (C.3).

The standard Heaviside function $H_e$, based on element, is introduced as an indicator of the utilization of the metal constitutive model based on $P_0$ field of level-set function $\psi_e$.

$$H_e(\psi_e) = \begin{cases} 1 & \text{if } \psi_e \ge 0 \\ 0 & \text{if } \psi_e < 0 \end{cases} \qquad (C.4)$$

where $H_e(\psi_e) = 1$ means that the element is in the gas (G) domain, conversely the element is located in the metal domain if $H_e(\psi_e) = 0$.

In the metal (M) domain, another Heaviside function, $I_{EVP}$, can be obtained from the element temperature, $T_e$, as given in Eq. (C.5). $I_{EVP} = 1$ corresponds to the EVP model in Eq. (15), and $I_{EVP} = 0$ represents the mushy or liquid state of the element (VP model in Eq. (24)).

$$I_{EVP} = H(T_c - T_e) = \begin{cases} 1 & \text{if } T_e < T_c \\ 0 & \text{if } T_e \ge T_c \end{cases} \qquad (C.5)$$

As a result, the stress field $\{s\}$ in the global frame can be expressed as:

$$\{s\} = H_e \{s^M\} + (1 - H_e) s^G, \quad \langle s_M \rangle = I_{EVP} s^{EV P} + (1 - I_{EVP}) s^{VP} \qquad (C.6)$$

This gives $\{s\} = \langle s^M \rangle$ or $s^G$, and $\langle s^M \rangle = s^{EV P}$ or $s^{VP}$. Similarly, $\{\partial s/\partial \dot{e}\}$ can be determined.

The discrete form in time of Eq. (C.3) can be obtained from the temporal discretization of the material constitutive equation in different states. Indeed, by considering two instants $t$ and $t + \Delta t$, an implicit finite difference scheme of the material constitutive equation can be obtained, where the EVP law is resolved in the crystal (local) reference frame.

#### C.2.1. Local resolution for the EVP behavior

In the crystal (local) reference frame, to satisfy Eq. (15), the total deviatoric strain rate is extracted with the Voigt notation:

$$\begin{aligned}
\tilde{\mathbf{e}} &= \tilde{\mathbf{e}}^{el} + \tilde{\mathbf{e}}^{vp} && \text{and} && \tilde{\mathbf{e}}^{el} &= \tilde{\mathbf{S}}^{el} \tilde{\mathbf{s}} \\
\tilde{\mathbf{e}}^{vp} &= \sum_{\alpha=1}^{n} \dot{\gamma}^{a} \tilde{\mathbf{P}}^{\alpha}
\end{aligned}
\qquad (C.7)$$

The Newton-Raphson method is applied to find the components of deviatoric tensor $\tilde{s}$ in Eq. (C.7) with the residual:

$$\tilde{\mathbf{R}}(\tilde{s}) = \tilde{\mathbf{e}} - \frac{1}{\Delta t} \tilde{\mathbf{S}}^{el} (\tilde{s} - \tilde{s}') - \sum_{\alpha=1}^{n} \dot{\gamma}^{a} \tilde{\mathbf{P}}^{\alpha} = \mathbf{0} \qquad (C.8)$$

At each iteration ($v + 1$), the updated estimation of $\tilde{s}$ is calculated by:
$$ \frac{\partial \tilde{\mathbf{R}}}{\partial \tilde{\mathbf{s}}} \delta \tilde{\mathbf{s}}^{(\nu+1)} = - \tilde{\mathbf{R}}(\tilde{\mathbf{s}}^{(\nu)}) \quad \text{and} \quad \tilde{\mathbf{s}}^{(\nu+1)} = \tilde{\mathbf{s}}^{(\nu)} + \delta \tilde{\mathbf{s}}^{(\nu+1)} \qquad (C.9) $$

where the stiffness tensor $\partial \tilde{\mathbf{R}}/\partial \tilde{\mathbf{s}}$ for Eq. (C.8) is:

$$ \frac{\partial \tilde{\mathbf{R}}}{\partial \tilde{\mathbf{s}}} = -\frac{1}{\Delta t} \tilde{\mathbf{s}}^{el} - \sum_{a=1}^{n} \frac{\dot{\gamma}_0}{(\tau_c^a)^{\frac{1}{m}}} \frac{1}{m} |\tau^a|^{\frac{1}{m}-1} \tilde{\mathbf{P}}^a \otimes \tilde{\mathbf{P}}^a \qquad (C.10) $$

The convergence criterion is defined as $\max_i (\text{abs}(\delta \tilde{s}_i^{(\nu+1)}/\tilde{s}_i^{(\nu+1)})) < 10^{-2}$, where i the component index of the $\tilde{s}$. Then the reference tangent modulus can be derived by:

$$ \frac{\partial \tilde{\mathbf{s}}}{\partial \tilde{\mathbf{e}}} = \left[ \frac{1}{\Delta t} \tilde{\mathbf{s}}^{el} + \tilde{\mathbf{G}} \right]^{-1} \mathbf{M}_{dev} \qquad (C.11) $$

where

$$ \tilde{\mathbf{G}} = \sum_{a=1}^{n} \frac{\dot{\gamma}_0}{(\tau_c^a)^{\frac{1}{m}}} \frac{1}{m} |\tau^a|^{\frac{1}{m}-1} (\tilde{\mathbf{P}}^a \otimes \tilde{\mathbf{P}}^a) $$

$$ \mathbf{M}_{dev} = \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} $$

where n is the number of slip systems. $M_{dev}$ is the deviatoric matrix, which converts the full strain rate (Voigt notation) to the deviatoric strain rate.

The bridge between the local crystal reference frame (X) and global reference frame (G) is built by the three Euler angles ($\phi_1, \phi, \phi_2$) of each local grain. According to Bunge's convention, the rotation matrix M is defined to convert the stress tensor from the crystal reference frame (X) and the global reference frame (G):

$$ \mathbf{M} = \begin{bmatrix} \cos \phi_1 \cos \phi_2 - \sin \phi_1 \cos \phi \sin \phi_2 & -\cos \phi_1 \sin \phi_2 - \sin \phi_1 \cos \phi \cos \phi_2 & \sin \phi_1 \sin \phi \\ \sin \phi_1 \cos \phi_2 + \cos \phi_1 \cos \phi \sin \phi_2 & -\sin \phi_1 \sin \phi_2 + \cos \phi_1 \cos \phi \cos \phi_2 & -\cos \phi_1 \sin \phi \\ \sin \phi_2 \sin \phi & \cos \phi_2 \sin \phi & \cos \phi \end{bmatrix} \qquad (C.12) $$

As a result, the stress tensor s in the global reference frame, $s_{|X}$, can be derived from the above obtained local stress, $s_{|X}$, by

$$ s_{|G} = M s_{|X} M^T \qquad (C.13) $$

Similar relationship can be found for the strain rate:

$$ \dot{\mathbf{e}}_{|G} = \mathbf{M} \dot{\mathbf{e}}_{|X} \mathbf{M}^T \qquad (C.14) $$

The transformation of the local tangent modulus $\partial s/\partial \dot{\mathbf{e}}$ obtained in the crystal reference frame (Eq. (C.11)), into its expression in the global reference frame is now detailed. This transformation is required to establish the consistent tangent stiffness of the global mechanical problem. From Eqs. (C.13) and (C.14), two matrices A and B are deduced, which define the transformation from the crystal to the global reference frame.

$$ \tilde{s}_{|G} = A \tilde{s}_{|X} \quad \text{and} \quad \tilde{\mathbf{e}}_{|G} = B \tilde{\mathbf{e}}_{|X} \qquad (C.15) $$

with

$$ \mathbf{A} = \begin{bmatrix} M_{11} M_{11} & M_{12} M_{12} & M_{13} M_{13} & 2 M_{12} M_{13} & 2 M_{13} M_{11} & 2 M_{11} M_{12} \\ M_{21} M_{21} & M_{22} M_{22} & M_{23} M_{23} & 2 M_{22} M_{23} & 2 M_{23} M_{21} & 2 M_{21} M_{22} \\ M_{31} M_{31} & M_{32} M_{32} & M_{33} M_{33} & 2 M_{32} M_{33} & 2 M_{33} M_{31} & 2 M_{31} M_{32} \\ M_{21} M_{31} & M_{22} M_{32} & M_{23} M_{33} & M_{22} M_{33} + M_{23} M_{32} & M_{23} M_{31} + M_{21} M_{33} & M_{21} M_{32} + M_{22} M_{31} \\ M_{31} M_{11} & M_{32} M_{12} & M_{33} M_{13} & M_{32} M_{13} + M_{33} M_{12} & M_{33} M_{11} + M_{31} M_{13} & M_{31} M_{12} + M_{32} M_{11} \\ M_{11} M_{21} & M_{12} M_{22} & M_{13} M_{23} & M_{12} M_{23} + M_{13} M_{22} & M_{13} M_{21} + M_{11} M_{23} & M_{11} M_{22} + M_{12} M_{21} \end{bmatrix} \qquad (C.16) $$

and

$$ \mathbf{B} = \begin{bmatrix} M_{11} M_{11} & M_{12} M_{12} & M_{13} M_{13} & M_{12} M_{13} & M_{13} M_{11} & M_{11} M_{12} \\ M_{21} M_{21} & M_{22} M_{22} & M_{23} M_{23} & M_{22} M_{23} & M_{23} M_{21} & M_{21} M_{22} \\ M_{31} M_{31} & M_{32} M_{32} & M_{33} M_{33} & M_{32} M_{33} + M_{33} M_{32} & M_{33} M_{31} + M_{21} M_{33} & M_{21} M_{32} + M_{22} M_{31} \\ M_{21} M_{31} & M_{22} M_{32} & M_{23} M_{33} & M_{22} M_{33} + M_{23} M_{32} & M_{23} M_{31} + M_{21} M_{33} & M_{21} M_{32} + M_{22} M_{31} \\ M_{31} M_{11} & M_{32} M_{12} & M_{33} M_{13} & M_{32} M_{13} + M_{33} M_{12} & M_{33} M_{11} + M_{31} M_{13} & M_{31} M_{12} + M_{32} M_{11} \\ M_{11} M_{21} & M_{12} M_{22} & M_{13} M_{23} & M_{12} M_{23} + M_{13} M_{22} & M_{13} M_{21} + M_{11} M_{23} & M_{11} M_{22} + M_{12} M_{21} \end{bmatrix} \qquad (C.17) $$

Therefore, we can write:

$$ \frac{\partial \tilde{s}_{|G_i}}{\partial \tilde{\mathbf{e}}_{|G_j}} = \frac{\partial \tilde{s}_{|G_i}}{\partial \tilde{\mathbf{e}}_{|X_i}} \cdot \frac{\partial \tilde{\mathbf{e}}_{|X_i}}{\partial \tilde{\mathbf{e}}_{|G_j}} = \frac{\partial}{\partial \tilde{\mathbf{e}}_{|X_i}} \left( A_{ik} \tilde{s}_{|X_k} \right) \cdot \frac{\partial}{\partial \tilde{\mathbf{e}}_{|G_j}} \left( (B^{-1})_{lm} \tilde{\mathbf{e}}_{|G_m} \right) \qquad (C.18) $$
$$ \frac{\partial \tilde{s}_{|G_i}}{\partial \tilde{\epsilon}_{|G_j}} = A_{ik} \cdot \frac{\partial \tilde{s}_{|X_k}}{\partial \tilde{\epsilon}_{|X_j}} \cdot (B^{-1})_{lm} \delta_{jm} = A_{ik} \cdot \frac{\partial \tilde{s}_{|X_k}}{\partial \tilde{\epsilon}_{|X_j}} \cdot (B^{-1})_{lj} \quad (C.19) $$

Finally, we have:

$$ \left( \frac{\partial \tilde{s}}{\partial \tilde{\epsilon}} \right)_{|G} = A \left( \frac{\partial \tilde{s}}{\partial \tilde{\epsilon}} \right)_{|X} B^{-1} \quad (C.20) $$

### C.2.2. Local resolution for the VP behavior

Based on Eqs. (24) and (25), the Norton-Hoff model is used for the purely viscoplastic constitutive law:

$$ s = 2K \left( \sqrt{3\tilde{\epsilon}} \right)^{m-1} \dot{\epsilon}^{vp} \quad (C.21) $$

when $m = 1$

$$ s = 2K \dot{\epsilon}^{vp} \quad (C.22) $$

The tangent modulus can be derived directly in the global reference frame:

$$ \frac{\partial \tilde{s}}{\partial \tilde{\epsilon}} = K_1 \dot{\epsilon}^{vp} \otimes \dot{\epsilon}^{vp} + 2K_2 (\mathbb{I} \otimes \mathbb{I} - \frac{1}{3} \mathbb{I} \otimes \mathbb{I}) \quad (C.23) $$

where $K_2$ and $K_1$ are derived by $K_2 = K(3\dot{\epsilon}^2)^{\frac{m-1}{2}}$, $K_1 = 4(m-1)K_2/3\dot{\epsilon}^2$.

Since anisotropies do not exist in VP states for mushy, liquid and gas domains (three Euler angles are set to zero), rotation operations are not necessary for either deviatoric stress or tangent modulus.

## Appendix D. Mesh sensitivity analysis for the RVE structure

The grain representations for RVE #1, using mesh $h_3 = 5$ µm and coarse mesh $h_5 = 30$ µm, are illustrated in Fig. D.3.

Fig. D.4 below illustrates the deviation between meshes $h_3 = 5$ µm and $h_4 = 10$ µm for the equiaxed polycrystalline cube (Fig. D.4.a), and between meshes $h_3 = 5$ µm and $h_5 = 30$ µm for RVE #1 (Fig. D.4.b-d).

<figure><img src="image_23.png" /><figcaption>Fig. D.3. Mesh sensitivity analysis is illustrated using IPFs and Euler angles for the columnar grains within RVE #1.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. D.4.** Traction test of an equiaxed polycrystalline cube with a random grain orientation in the z direction (a); and the traction tests for RVE#1: (b) x direction, (c) y direction, and (d) z direction.</figcaption></figure>

References

[1] P. Stavropoulos, P. Foteinopoulos, A. Papacharalampopoulos, H. Bikas, Addressing the challenges for the industrial application of additive manufacturing: towards a hybrid solution, Int. J. Lightweight Mater. Manuf. 1 (3) (2018) 157-168. https://doi.org/10.1016/j.ijlmm.2018.07.002

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

[3] 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 (101804) (2023). https://doi.org/10.1016/j.mtla.2023.101804

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

[5] D. Kats, Z. Wang, Z. Gan, W.K. Liu, G.J. Wagner, Y. Lian, A physics-informed machine learning method for predicting grain structure characteristics in directed energy deposition, Comput. Mater. Sci. 202 (2022) 110958. https://doi.org/10.1016/j.commatsci.2021.110958

[6] M. Chiumenti, E. Neiva, E. Salsi, M. Cervera, S. Badia, J. Moya, Z. Chen, C. Lee, C. Davies, Numerical modelling and experimental validation in selective laser melting, Addit. Manuf. 18 (2017) 171-185. https://doi.org/10.1016/j.addma.2017.09.002

[7] Y. Zhang, G. Guillemot, M. Bernacki, M. Bellet, Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process, Comput. Methods Appl. Mech. Eng. 331 (2018) 514-535. https://doi.org/10.1016/j.cma.2017.12.003

[8] Y. Zhang, C.-A. Gandin, M. Bellet, Finite element modeling of powder bed fusion at part scale by a super-layer deposition method based on level set and mesh adaptation, J. Manuf. Sci. Eng. 144 (2022) 051001-1-051001-15. https://doi.org/10.1115/1.4052386

[9] D. Raabe, P. Klose, B. Engl, K.P. Imlau, F. Friedel, F. Roters, Concepts for integrating plastic anisotropy into metal forming simulations, Adv. Eng. Mater. 4(4) (2002) 9169-9180. https://doi.org/10.1002/1527-2648(200204)4:4<169::AID-ADEM169>3.0.CO;2-G

[10] M. Sachtleber, Z. Zhao, D. Raabe, Experimental investigation of plastic grain interaction, Mater. Sci. Eng. 336(1-2) (2002) 81-87. https://doi.org/10.1016/S0921-5093(01)01974-8

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

[12] D. Peirce, R.J. Asaro, A. Needleman, Material rate dependence and localized deformation in crystalline solids, Acta Metall. 31(12) (1983) 1951-1976. https://doi.org/10.1016/0001-6160(83)90014-7

[13] G. Cailletaud, Une approche micromécanique du comportement des polycristaux, Revue de physique appliquée 23(4) (1988) 35-365. https://doi.org/10.1051/rphysap:01988002304035300

[14] M. Anahid, M.K. Samal, S. Ghosh, Dwell fatigue crack nucleation model based on crystal plasticity finite element simulations of polycrystalline titanium alloys, J. Mech. Phys. Solids 59(10) (2011) 2157-2176. https://doi.org/10.1016/j.jmps.2011.05.003

[15] U.F. Kocks, A.S. Argon, M.F. Ashby, Thermodynamics and Kinetics of Slip, Pergamon Press (1975).

[16] J. Cheng, S. Ghosh, A crystal plasticity FE model for deformation with twin nucleation in magnesium alloys, Int. J. Plast. 67 (2015) 148-170. https://doi.org/10.1016/j.ijplas.2014.10.005

[17] G.I. Taylor, Plastic strain in metals, J. Inst. Metals 62 (1938) 307-324.

[18] A. Molinari, G.R. Canova, S. Ahzi, A self consistent approach of the large deformation polycrystal viscoplasticity, Acta Metall. 35(12) (1987) 2983-2994. https://doi.org/10.1016/0001-6160(87)90297-5
[19] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. London Ser. A. Math. Phys. Sci. 241 (1226) (1957) 376-396. https://doi.org/10.1098/rspa.1957.0133

[20] R.A. Lebensohn, C.N. Tomé, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys, Acta Metall. Mater. 41 (9) (1993) 2611-2624. https://doi.org/10.1016/0956-7151(93)90130- K

[21] P. Hu, Y. Liu, Y. Zhu, L. Ying, Crystal plasticity extended models based on thermal mechanism and damage functions: application to multiscale modeling of aluminum alloy tensile behavior, Int. J. Plast. 86 (2016) 1-25. https://doi.org/10.1016/j.ijplas.2016.07.001

[22] J. Li, I. Romero, J. Segurado, Development of a thermo-mechanically coupled crystal plasticity modeling framework: application to polycrystalline homogenization, Int. J. Plast. 119 (2019) 313-330. https://doi.org/10.1016/j.ijplas.2019.04.008

[23] S. Ghorbanpour, M.E. Alam, N.C. Ferreri, A. Kumar, B.A. McWilliams, S.C. Vogel, J. Bicknell, I.J. Beyerlein, M. Knezevic, Experimental characterization and crystal plasticity modeling of anisotropy, tension-compression asymmetry, and texture evolution of additively manufactured inconel 718 at room and elevated temperatures, Int. J. Plast. 125 (2020) 63-79. https://doi.org/10.1016/j.ijplas.2019.09.002

[24] 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 (2022) 825-845. https://doi.org/10.1007/s00466-021-02116-z

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

[26] Y.M. Wang, T. Voisin, J.T. McKeown, J. Ye, N.P. Calza, Z. Li, T. Zhu, Additively manufactured hierarchical stainless steels with high strength and ductility, Nat. Mater. 17(1) (2018) 63-71. https://doi.org/10.1038/nmat5021

[27] M. Lindroos, T. Pinomaa, A. Antikainen, J. Lagerbom, J. Reijonen, T. Lindroos, 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. https://doi.org/10.1016/j.addma.2020.101819

[28] A. Pautrat, Étude expérimentale quantitative de la solidification de l'inconel 718 en fonderie, Ph.D. thesis, Ecole Nationale Supérieure des Mines de Paris, 2013. https://pastel.hal.science/pastel-00998532.

[29] Q. Chen, G. Guillemot, C.A. Gandin, M. Bellet, Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials, Addit. Manuf. 16 (2017) 124-137. https://doi.org/10.1016/j.addma.2017.02.005

[30] CimLib library, APP Deposit. 2007, Deposit No: IDDN.FR.001.460005.002.S.A.2007.000.10000.

[31] 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, Comptes Rendus. Mécanique 346 (11) (2018) 1055-1071. https://doi.org/10.1016/j.crme.2018.08.008

[32] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114(1) (1994) 146-159. https://doi.org/10.1006/jcph.1994.1155

[33] A. Aalilija, C.-A. Gandin, E. Hachem, A simple and efficient numerical model for thermal contact resistance based on diffuse interface immersed boundary method, Int. J. Therm. Sci. 166 (2021) 106817. https://doi.org/10.1016/j.ijthermalsci.2020.106817

[34] Y. Zhang, H. Behlal, C.-A. Gandin, O. Senninger, G. Guillemot, M. Bellet, Homogenization methods for thermal study of support structure in laser powder bed fusion (L-PBF) — application to process numerical modeling, Int. J. Numer. Methods Fluid Flow 35 (1), (2025) 358-381. https://doi.org/10.1108/HFF-09-2024-0683

[35] M. Shakoor, M. Bernacki, P.O. Bouchard, Three-dimensional numerical modeling of ductile fracture mechanisms at the microscale, 4th ECCOMAS Young Investigator Conference-YIC (2017). https://minesparis-psl.hal.science/hal-01626736v1.

[36] K.C. Mills, Recommended Values of Thermophysical Properties for Selected Commercial Alloys, Woodhead publishing (2002).

[37] E.W. Lemmon, R.T. Jacobsen, Viscosity and thermal conductivity equations for nitrogen, oxygen, argon, and air, Int. J. Thermophys. 25 (2004) 21-69. https://doi.org/10.1023/B:IJOT.0000022327.04529.f3

[38] 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. https://doi.org/10.1016/j.apmt.2017.08.006

[39] P. Zehner, E.U. Schlünder, Wärmeleitfähigkeit von schüttungen bei mäßigen temperaturen, Chem. Ing. Tech. 42 (14) (1970) 933-941. https://doi.org/10.1002/cite.330421408

[40] A. Saad, C.-A. Gandin, M. Bellet, Temperature-based energy solver coupled with tabulated thermodynamic properties — Application to the prediction of macrosegmentation in multicomponent alloys, Comput. Mater. Sci 99 (2015) 221-231. https://doi.org/10.1016/j.commatsci.2014.12.009.

[41] V.D. Fachinotti, M. Bellet, Linear tetrahedral finite elements for thermal shock problems, Int. J. Numer. Methods Fluid Flow 16 (2006) 590-601. https://doi.org/10.1108/09615530610669120.

[42] 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. Eng. 356 (2019) 294-324. https://doi.org/10.1016/j.cma.2019.07.006

[43] J.R. Rice, Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity, J. Mech. Phys. Solids 19(6) (1971) 433-455. https://doi.org/10.1016/0022-5096(71)90010-X

[44] J.W. Hutchinson, Creep and plasticity of hexagonal polycrystals as related to single crystal slip, Metall. Trans. A 8 (9) (1977) 1465-1469. https://doi.org/10.1007/BF02642860

[45] C. Kumara, D. Deng, J. Moverare, P. Nylen, Modelling of anisotropic elastic properties in alloy 718 built by electron beam melting, Mater. Sci. Technol. 34 (5) (2018) 529-537. https://doi.org/10.1080/02670836.2018.1426258

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

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

[48] M. Bellet, V.D. Fachinotti, ALE Method for solidification modelling, Comput. Methods Appl. Mech. Eng. 193 (2004) 4355-4381. https://doi.org/10.1016/j.cma.2003.11.016

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

[50] Y. Chen, Y. Wang, Y. Lu, B. Han, K. Huang, X. Fang, Q. Zhang, A multiscale model for efficiently simulating laser powder bed fusion process with detailed microstructure and mechanical performance results, Mater. Sci. Eng. A 918 (2024) 147435. https://doi.org/10.1016/j.msea.2024.147435

[51] P. Haldipur, F.J. Margetan, R.B. Thompson, Estimation of single-crystal elastic constants from ultrasonic measurements on polycrystalline specimens, AIP Conf. Proc. 700 (1) (2004) 1061-1068. https://doi.org/10.1063/1.1711735

[52] G. Martin, N. Ochoa, K. Sai, E. Hervé-Luanco, G. Cailletaud, A multiscale model for the elastoviscoplastic behavior of directionally solidified alloys: application to FE structural computations, Int. J. Solids Struct. 51(5) (2014) 1175-1187. https://doi.org/10.1016/j.ijsolstr.2013.12.013

[53] S. Jothi, S.V. Merzlikin, T.N. Croft, J. Andersson, S.G.R. Brown, An investigation of micro-mechanisms in hydrogen induced cracking in nickel-based superalloy 718, J. Alloys Compd. 664 (2016) 664-681. https://doi.org/10.1016/j.jallcom.2016.01.033

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

[55] M.R. Gotterbarm, M. Seifi, D. Melzer, J. Düzgan, A.A. Salem, Z.H. Liu, C. Koerner, Small scale testing of IN718 single crystals manufactured by EB-PBF, Addit. Manuf. 36 (101449) (2020). https://doi.org/10.1016/j.addma.2020.101449

[56] A. Ranellucci, Slic3r - Open source 3D printing toolbox (2023). https://slic3r.org/.

[57] A. Queva, G. Guillemot, C. Moriconi, C. Metton, M. Bellet, Numerical study of the impact of vaporisation on melt pool dynamics in laser powder bed fusion-application to IN718 and ti-6Al-4V, Addit. Manuf. 35 (2020) 101249. https://doi.org/10.1016/j.addma.2020.101249
[58] Q. Chen, Modélisation numérique thermomécanique de fabrication additive par fusion sélective de lit de poudre par laser: Application aux matériaux céramiques, Ph.D. thesis, Paris Sciences et Lettres (ComUE), 2018. https://pastel.hal.science/tel-01814209.

[59] K. Morgan, R.W. Lewis, O.C. Zienkiewicz, An improved algorithm for heat conduction problems with phase change, Int. J. Numer. Methods Eng. 12 (7) (1978) 1191–1195. https://doi.org/10.1002/nme.1620120710.

[60] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the stokes equations, CALCOLO 21 (4) (1984) 337–344. https://doi.org/10.1007/BF02576171

[61] I. Babuška, Error-bounds for finite element method, Numer. Math. 16 (4) (1971) 322–333. https://doi.org/10.1007/BF02165003.

[62] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers, Rev. Fr. Autom. Inf. Rech. Oper. 8 (1974) 129–151. https://doi.org/10.1051/m2an/197408R201291

[63] O. Jaouen, Modélisation tridimensionnelle par éléments finis pour l'analyse thermo-mécanique du refroidissement des pièces coulées, Ph.D. thesis, École Nationale Supérieure des Mines de Paris, 1998. https://pastel.hal.science/pastel-00688778/.