Crystal plasticity model of residual stress in additive manufacturing

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


Introduction
Additive manufacturing (AM) is becoming increasingly popular in industries. Metal AM processes such as selective laser melting (SLM) [1], selective electron beam melting (SEBM) [2] and selective laser sintering (SLS) [3] are layered approaches fabricating the products with complex shapes [4]. One of the bottlenecks for metal AM is the thermal stress caused by the high temperature gradients during the high-frequency heating/cooling cycles. High residual stress may result in defects [5,6] such as the part distortion [7], delamination [8] and poor fracture resistance [9,10,11]. The investigation of temperature history and stress evolution is therefore necessary for optimising the manufacturing parameters during the AM process. Due to the large temperature gradients and high-frequency heating/cooling cycles during the manufacturing process, it is difficult to monitor the temperature and thermal stress in real time during experiments. For the residual stress measurement, the experimental approaches, such as the neutron diffraction [12,13], X-ray diffraction [14], contour method [15] and semi-destructive hole drilling method [16,17], are costly and sometimes limited to surface measurements. Numerical simulations are efficient to study the mechanical behaviour during the manufacturing process. The finite element method (FEM) has been widely used in simulating thermal stresses in AM [18]. One type of such models is based on the inherent strain method, where the deposited regions are activated with predefined strains to calculate the final deformation and residual stress of the part [19]. These models are computationally efficient, but the accuracy may be compromised due to the uncertainties of the inherent strains. Another type is thermo-mechanical coupling, where the temperature profile is calculated by solving the governing equations of heat transfer and then serves as the load for thermal stress calculation [20,21]. Thermo-mechanical models are mostly constructed based on many geometrical simplifications and physical assumptions as the heat transfer part does not incorporate the thermal-fluid flow effects [22]. Most of the FEM simulations failed to resolve the track morphology and molten pool evolution as all the tracks are assumed to be uniform. Therefore, those simulations are only able to predict the overall part distortion at the macroscopic scales, but cannot explain the thermal stress at the microscale. An emerging type of the thermal stress models is the coupling between the computational fluid dynamics (CFD) and FEM models [23,24,25]. These coupled CFD-FEM models analyse the thermal stress and deformation using the high-accuracy temperature profiles output from the CFD simulations, particularly the powder-scale thermal-fluid flow simulations [26], instead of those from simplified analytical models [27]. Cheon et al. [24] proposed a CFD-FEM framework for welding to analyse the thermal stress along a single welding line, but the model did not consider the residual stress reset due to melting. Bailey et al. [23] implemented the temperature and surface profiles from thermal-fluid flow simulations into an FEM residual stress model to simulate the multi-track manufacturing cases of the laser directed energy deposition (DED) process. The element activation and deactivation method in ABAQUS was used to simulate the melting and solidification, which is troublesome to input command lines to activate elements of the solidified materials in the input file and frequently encounters problems on numerical convergence and numerically-induced deformation of the newly-activated elements. In our previous work [25], we simulated the multi-layer multi-track SLM process by mapping the temperature profiles from the high-fidelity thermal-fluid flow simulation into the FEM model, enabling the high resolution of the molten pool evolution and the morphology (rough surfaces and inner voids) of the FEM model. The quiet element method was used to simulate the melting and solidification, where the melted materials are assigned with null material properties. There are two critical but common problems in the aforementioned thermo-mechanical models and CFD-FEM models: 1. How to represent the melting and solidification while avoiding numerically-induced deformation and numerical divergence. In the previous models, two methods are commonly used: first, the element activation method [28,29,30,31,32] can activate elements for solidified materials but cannot deactivate elements again due to remelting, while the surface nodal temperature obtained by the interpolation between active and inactive elements [33] can cause significant inaccuracy and divergence problems. Second, the quiet element method [32,33,34], which is essentially the same as the stiffness degradation method commonly used in fracture mechanics [35,36,37,38,39], can remove the effect of the liquid phase on the solidified region, however upon solidification the liquid phase becomes solid with initial deformation that is incompatible with the deformation of the surrounding solid region, which induces un-physical deformation and may also cause numerical divergence. 2. None of the previous models has considered the detailed grain structures and the resultant temperature-dependent anisotropic mechanical properties at the microscale. Using simplified homogenized constitutive laws, the previous models cannot resolve the micro-scale residual stress [40], which can result in unique mechanical properties of the AM parts significantly differing from those by the conventional manufacturing [41].
In this paper, a crystal plasticity finite element model (CPFEM) [42,43] is developed, incorporating the temperature profile from the thermal-fluid flow simulation [25] and grain structures from the grain growth simulation [44,45] and experimental data [40], as to be described in Section 2.
Particularly, an element elimination and reactivation method is proposed to model the melting and solidification. More importantly, the algorithm can reinitialise the state variables, as tested in a simple case in Section 3, such as the plastic deformation of the solidifying elements, which is critical to avoid the pre-existing plastic deformation just after solidification. Simulation cases of SLM of 316 stainless steel are conducted as described in Section 4. The proposed element elimination and reactivation method and the commonly used stiffness degradation method are compared in terms of prediction accuracy, convergence and computational efficiency. Based on the simulation results, the formation of residual stresses in different grains and the correlation between plastic deformation, depth and residual stress are discussed in Section 5.

Crystal plasticity framework
A crystal plasticity material model including thermal response is used for this work. The deformation gradient F can be decomposed into elastic, thermal and plastic parts [46,47,48,49]: where F e is the elastic part including the stretching and rigid body rotation, F th is thermal part and F p is the irreversible plastic deformation, which evolves according to [50,51]: whereγ α is the shear rate of slip system α. The vectors m α and n α are unit vectors which describe the slip direction and normal of slip plane α, respectively. n is the number of slip systems in active state.

3
The plastic strain rate on slip system α of face-centered cubic metals can be calculated by a powerlaw equation [52,53]:γ where τ α is the resolved shear stress, τ α c is the temperature dependent critical resolved shear stress (CRSS),γ 0 and m are constants, which indicate the reference strain rate and rate sensitivity of slip, respectively. An exponential model calibrated by experiments can be used to calculate the CRSS change with temperature decrease of each slip system [54,55]: where τ α c (T ) is the CRSS of slip system α at temperature T , T 0 is the reference temperature, k A , k B and k C are constants. For slip system α, its hardening behaviour is influenced by other slip systems β as [56]: where h αβ is the hardening matrix given by [57,58]: where q αβ is a hardening coefficient matrix for latent hardening, h 0 is an initial hardening term, τ sat is the saturation slip resistance and a is a constant. The parameters of the crystal plasticity model are calibrated using the experimental data in [40]. For the thermo-mechanical response during the laser scan process, the Green-Lagrange strain tensor can be expressed by the equation composed of thermal and elastic deformation gradients as: where I is the identity matrix. To describe how the size of the substance changes with the temperature, the linear thermal expansion coefficient α l of 316L SS is introduced [59]. The volumetric thermal expansion coefficient can be expressed by α v = 3α l [60]. Since the thermal expansion coefficient of 316L stainless steel varies slightly with temperature [61,62], and the temperature dependence is given by fitting the experimental data into a linear equation [63]: where α 0 denotes the thermal expansion coefficient at room temperature and α 1 is the linear coefficient of the temperature dependence of α v . The infinitesimal volumetric change dV from T 0 to T can be described by [48]: where V is the volume of the region and T 0 is the reference temperature (room temperature). Substituting equation (8) into equation (9), the thermal eigenstrain tensor α can be expressed by [64]: The 2nd Piola-Kirchhoff stress can be calculated by: where C is the rank four elasticity tensor of 316L stainless steel [65,66]. The elasticity tensor of the solid is temperature dependent [54], which can be described by a linear function for the components: where C solid i j (T 0 ) is the elasticity tensor component at room temperature, and the derivative term indicates the change of the components caused by temperature change. Phase change can be modelled by changing the stiffness tensor C, as to be explained in section 2.2. Since positive strain represents tension and negative strain represents compression, the same applies to the stress tensor components. In the following, compression will appear as a negative stress component while tension will appear as a positive stress component. The temperature T is obtained from thermal-fluid flow simulations [25]. It should be mentioned that the temperature mapping from thermal-fluid to mechanical simulation is one-way, and the heat transfer analysis is merely conducted in the high-fidelity thermal-fluid flow simulation, from which the thermal information are extracted for the crystal plasticity calculation. Thus, deactivation of the melted parts will not affect the temperature distribution in the representative volume.

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

Phase 1 Phase 2 Stiffness
Solid Solid Newton-Raphson approach to find the Cauchy stress tensor σ is applied. To obtain the 2nd Piola-Kirchhoff stress in equation (11), a function ψ = S − C(E e − α) is defined, the absolute value of which must be minimised. An iterative return mapping algorithm, in which the stress S is the variable is implemented. The stress increment during the return mapping algorithm can be written as [69,70]: where i and i + 1 indicate two subsequent iterations. The Cauchy stress σ is obtained by [71]: The Cauchy stress is used by the finite element solver to find the displacement field u that leads to stress equilibrium. The parameters used in the simulations and the elastic constants of 316L stainless steel in Voigt notation are listed in Table 2.

Implementation of the element elimination and reactivation method
To model the melting and solidification of material in the SLM process, an element elimination and reactivation approach is employed. This approach is based on a moving subdomain paradigm that is originally developed [74], which is based on MOOSE's subdomain restricted system [68]. Specifically, the physical domain is divided into an active subdomain (Ω a ) and an inactive subdomain (Ω i ). The thermo-mechanical finite element computation is defined and carried out in Ω a , while Ω i only carries the geometric information. In other words, the elements in the inactive domain, or inactive elements, do not contribute to the equation system that calculates stress equilibrium. The active and inactive subdomains are determined only by the prescribed temperature field. At a typical FEM time step, any active element T ∈ Ω a with an average temperature that is above the melting temperature (i.e., T avg > T m ) or below the gas temperature (i.e., T avg < T g ) is moved from the active subdomain to the inactive one, i.e.: This process mimics the solid material being eliminated from the physical domain, which happens when the material is melted around the laser beam, and is therefore referred to as element elimination. Similarly, when the average temperature of one inactive element T ∈ Ω i falls to the solid temperature range (i.e., T g ≤ T avg ≤ T m ), it is moved to the active subdomain, i.e. : This mimics the process of material being added back to the solid physical domain, which happens when the material solidifies behind the laser beam, thus is referred to as element reactivation. After elements being transferred between the two subdomains, several tasks need to be accomplished before the subsequent computation. First, the boundary information is updated based on the coverage of the new subdomains. Here, the subdomain boundaries are updated for both Ω a and Ω i . The boundary conditions are projected to the updated boundaries correspondingly. Second, the initial conditions of the displacement variables and the material properties that define the crystal plasticity behaviour, such as the plastic deformation (see Section 2.1), are projected to the quadrature points of the reactivated elements. Projection of the initial condition indicates that the displacement, strain, and stress calculations are reset to zero. This allows for the restart of the plastic deformation calculation during the simulation, which cannot be implemented easily in the stiffness degradation method [75]. Indeed, if the elements are not eliminated and reactivated, they can still carry a certain amount of load. In the stiffness degradation method, a sudden change of the plastic deformation due to melting and solidification causes convergence problems because an elastic deformation jump, and consequent stress discontinuity, occur. Finally, the element elimination and reactivation method resolves the issue that is brought by large deformation from the reactivated elements if they are set to carry the displacement history upon reactivation.

Test case: element elimination during tension
The element elimination and reactivation method is tested on a simple example geometry, as shown in Figure 2. The aim of this test is to demonstrate that inactive elements do not contribute to the stress equilibrium calculations, i.e. they do not exert stress on the active elements, and to prove 8 that boundary conditions are implemented correctly on the continuously changing free surface. A cubic representative volume is constituted of 10 × 10 × 10 elements and has a size of 1 µm × 1 µm × 1 µm. Boundary condition u x = 0 is applied on the surface x = 0, u y = 0 on y = 0 and y = 1 µm, u z = 0 on z = 0 and z = 1 µm. Displacement boundary condition is imposed on the surface x = 1 µm. u x is increased up to 0.001 µm on that surface before the element elimination algorithm starts. The element elimination algorithm is based on a prescribed temperature field, which is chosen in such a way that elements from the surface z = 1.0 µm are progressively removed, down to the surface z = 0.5 µm. The displacement u x on the surface x = 1 µm is kept constant during element elimination, as shown in Figure 2(a)-(c), therefore ε xx remains constant at the value of 0.001. Plasticity and thermal eigenstrain are not activated in this simulation case, and the elastic constants are set as C 11 = 1.5 GPa and C 12 = 0.75 GPa. ε yy and the off diagonal strain components are zero. When no element elimination has taken place, the boundary condition u z = 0 on the surface z = 1 µm holds, therefore ε zz = 0 and u z = 0 everywhere in the representative volume, as shown in Figure 2(d). This lateral constraint is reflected by the large positive value of the stress component σ zz in Figure 2(g). Therefore, a tensile stress is present on that surface. After the element elimination algorithm starts, the top surface perpendicular to the z axis becomes unconstrained. Therefore, lateral contraction takes place because of the constant displacement along the x axis. The lateral contraction can be calculated as follows: ε zz = − (C 12 /C 11 ) ε xx = −0.0005. For example, in Figure 2

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

Simulation convergence
The convergence speed of the simulations using the stiffness degradation method depends mainly on the residual stiffness parameter q r . A value smaller than 0.001 significantly increases the simulation time, but there is very small difference in the plastic deformation between the simulations with q r = 0.01 and q r = 0.001, as will be shown in the following. Simulations with both q r = 0.01 and q r = 0.001 are carried out. If the element elimination method is used, the time step must be decreased immediately before and after the temperature steps at which element elimination takes place, as described in section 2. In this way, good convergence can be obtained. Specifically, the time step is reduced to 0.1 µs in an interval of 2 µs around each temperature step, i.e. every 45 µs. For instance, considering the first temperature step, the time step is reduced to 0.1 µs at t = 44 µs and it is increased again to 1 µs after t = 46 µs. This time step reduction applies only during the laser scan, while time step is kept constant at 1 µs during the subsequent cooling stage, which is the second part of the simulation, carried out after the laser beam has exited the geometry. By contrast, in the stiffness degradation method, the value of the time step is always 1 µs. Since the eliminated elements do not contribute to the residual or Jacobian calculation of the solver, the computation of each time step is slightly faster if the element elimination method is used. However, a reduction of the time step size before and after the element elimination steps is necessary in that case. Therefore, the simulations carried out with the element elimination method are not necessarily faster, unless the number of eliminated elements is a significant fraction of the total number of elements. This is shown by the computational costs reported in table 3. As stated before, the element elimination steps require a reduction of the time step, however,   Figure 5(a). Normally, reactivated elements have zero displacement because the displacement variable is reinitialised. In the MOOSE framework, the displacement in the reinitialised nodes is assigned in such a way as to get an av-erage zero displacement at the eight integration points, as shown in Figure 5(b). Therefore, when an element is reactivated with zero displacement, the excessive distortion leads to divergence, as shown in Figure 5(c). Several attempts have been carried out in which the initial conditions on the displacement after element reactivation are changed as a function of the coordinates. Surprisingly, a simple solution is to reinitialise u z = 3 µm when elements are reactivated, while u x and u y are reinitialised to zero everywhere. This displacement value is closer to the average displacement of the molten pool region on the top surface in the simulation using the residual stiffness method, while it also leads to good convergence. Other reinitialisation values for u z , such as 1 µm, 2 µm, 5 µm and 6 µm have been used, but they have lead to excessive element distortion just after reinitialisation. Only values between 3 µm and 4 µm remove the element distortion problem, as shown in Figure 5(d). This implementation is made possible by the MOOSE framework because the initial conditions can be set as a function of time. At t = 0, the displacement is initialised to zero, while for t > 0, the components u x , u y are reinitialised to zero and u z is reinitialised to 3 µm everywhere in the geometry. This strategy can be applied to any additive manufacturing simulation and can be extended by setting the initial conditions as a more complex function of space and time.

Comparison between the computational methods
The two computational methods presented in sections 2.2 and 2.3 are compared. The plastic deformation is a very important quantity because it determines the residual stress after the laser scan. The diagonal components of F p − I are averaged over the central part of the representative volume, constituted of the highlighted elements in Figure 6(d). This procedure is carried out for simulations using the element elimination method and the stiffness degradation method. Figures 6(a)-(c) show the comparison of simulation results by the two methods. During laser scan, the hot region undergoes thermal expansion and applies compression on the neighbouring regions. Therefore, plastic compression along the x and y axes takes place in the central region. Since plastic deformation is isochoric, the elements expand along the z axis, as shown in Figure 6(c). The plastic deformation is always larger using the element elimination and reactivation method. If the stiffness degradation method is used, the residual plastic deformation component along the y axis, after laser scan and cooling, is underestimated by about 20%. The difference is particularly important at time 400-500 µs, when the laser just passes the representative volume but the melting pool is still deep. In this time interval, the stiffness degradation model underestimates the plastic deformation by a factor 2.
Overall, in the stiffness degradation method, the presence of low stiffness elements prevents the elements near the melting pool to accommodate part of the plastic deformation. Simulations with both q r = 0.01 and q r = 0.001 are carried out and the plastic deformation along the y axis is shown in Figure 7 when the laser beam is approximately in the centre of the representative volume. The strong similarity between the plastic deformation in the two simulations confirms that the range of values chosen for the parameter q r does not affect the simulation results. The other components of the plastic deformation gradient are also very similar, therefore a value q r = 0.001 is suitable for this kind of simulations.

Grain orientation and plastic deformation
The element elimination and reactivation method will be used in the following simulation results to understand the correlation between plastic deformation, residual stress and grain orientation. Firstly, the correlation between residual plastic deformation (after the laser scan and cooling) and grain orientation is analysed. The scatter plots in Figure 8(a)-(c) are obtained as follows: the plastic deformation components are averaged over each element in the representative volume and one element is selected every 16 µm and represented as a point in the scatter plots. The distribution of these elements is uniform in the representative volume and this algorithm allows to select a representative set of elements that do not fill completely the scatter plot. Given the Euler angles in each selected element, the maximum Schmid factor among the slip systems is calculated for load along a particular direction (x, y or z axes respectively) [79]. This is motivated by the fact that, despite of the complex load undergone by different regions, the stress component along one direction should in principle induce plastic deformation along that particular direction. The colour of the points in Figure 8(a)-(c) corresponds to the position of the element along the z axis. As shown in Figure 8, the maximum plastic deformations occur on some points that are at an intermediate depth (about 40 µm to 80 µm). This is related to the solidification and thermal stress around the molten pool, as also revealed in Figure 11. Moreover, it is clear that a higher number of grains with larger maximum Schmid factors along a specific direction (x, y or z) tends to have larger plastic deformation along that direction. However, there are grains with large maximum Schmid factors that show little plastic deformation. Many of this kind of points are dark blue in Figure 8(a)-(c), which means they are near the substrate. Because the substrate is fixed during simulation, the plastic deformation near the substrate is limited. Nevertheless, there are still some points at intermediate depth that do not have much plastic deformation, this is quite surprising and reflects the complexity of the thermal load in different regions. This behaviour will be clarified in the following, when the plastic deformation during laser scan is shown near the melting pool. There are also some grains with larger Schmid factor that show little plastic deformation and are located near the top stress-free surface (dark red points), as shown in Figure 8(a)-(c). To investigate the directionality of the plastic deformation, the range of the plastic deformation components is investigated, which shows a greater correlation with the maximum Schmid factor [80], as shown in Figure 9(a)-(c). The range is calculated as the difference between maximum and minimum values of the plastic strain components for a given maximum Schmid factor. This correlation does not depend strongly on the depth z. The majority of the grains undergo plastic compression along the x and y axes, while plastic tension is observed along the z axis. This is

Plastic deformation and residual stress
The residual stress, after laser scan and cooling, is due to the heterogeneity of the plastic strain in different grains. In order to maintain strain compatibility, the elastic deformation must compensate the difference in plastic deformation, therefore residual stresses are generated, as modelled by equations (1) and (11). The correlation between residual stress and plastic deformation can be visualised using the scatter plots in Figure 10. Each point represents an element in the representative volume, selected every 16 µm, with its corresponding residual stress and plastic strain, averaged over the integration points. The colour of each point represents the z coordinate of the corresponding element. It is certainly true that there is a correlation between residual stress and plastic strain. Compressive plastic strain along the x and y axes corresponds to tensile residual stress components along those directions, while tensile plastic strain along the z axis corresponds to compressive residual stress along z. However, few points show the opposite behaviour: for instance, a tensile residual stress σ yy is associated with a small tensile plastic strain along the y axis in some regions. Overall, the tensile residual stress components σ xx and σ yy are similar and large, therefore these pre-tensions can induce fracture along both the x and y directions. The plastic deformation depends also on the depth. In Figure 10, it is evident that the regions at the top of the representative volume, represented by dark red points, are not the ones accommodating the largest plastic strain. The points with the largest plastic deformation are the ones with depth in the interval 30 µm < z < 80 µm. These are the points that are closer to the bottom of the melting pool during laser scan. However, the points at the top surface can have residual stresses that are similar to points that undergo larger plastic deformations. This shows that the influence of the plastic deformation in the substrate is very important; it will be discussed in more details in section 5. It is important to note that some elements show values of the residual stress components that are larger than the critical resolved shear stress to induce plastic deformation. This is due to large volumetric stress in those elements. Moreover, some dark red points in Figure 10(c) have residual stress σ zz that is different from zero. These points are picked just below the surface, in the coordinate interval 90 µm < z < 100 µm. Indeed the stress component σ zz on the free top surface is very close to zero, as shown in Figure 12(c). In order to understand the plastic deformation in the representative volume induced by the laser scan, the components of F p are shown in Figure 11(a)-(c) when the laser beam is approximately in the centre of the representative volume. The thermal expansion ahead of the laser beam leads to plastic compression along the x axis. This is visible on the two sides and at the bottom of the melting pool in Figure 11(a), where blue areas are present. By contrast, the plastic deformation along the y axis is compressive at the bottom of the melting pool but it is tensile on the two sides, as shown by the red regions in Figure 11(b). Since plastic deformation is isochoric, the region at the bottom of the melting pool expands along the z axis, as shown in Figure 11(c). Because of the subsequent solidification, the melting pool becomes more shallow and the plastic deformation concentration shown in Figure 11(a)-(c) spreads in the central region of the representative volume. This is the reason why the regions with the largest plastic deformation are concentrated at depth 30 µm < z < 80 µm, as shown in Figure 10. At the top free surface, the plastic deformation is slightly smaller than in the centre of the representative volume, as shown in Figure 8(a)-(c). This is the reason why not all grains with larger Schmid factor at the top surface can accommodate large plastic deformations, as stated in section 4.4. The diagonal components of the residual Cauchy stress tensor are shown in Figures 12(a)-(c). σ xx and σ yy are heterogeneous at the top surface and in the centre of the representative volume, as shown in Figures 12(a) and (b). This shows that they are strongly affected by the grain orientation and that the substrate is applying an important constraint along x and y on the top surface. In fact, even if the top surface is free, large values of σ xx and σ yy are present. The compressive stress σ zz is concentrated in the substrate but it is not very large at the top surface because of the free boundary condition.

Discussion
The simulations suggest the following mechanism for the build up of the residual stresses, as shown in Figure 13. The thermal expansion around the molten pool induces plastic compression along the horizontal directions, x and y, and plastic tension along the out-of-plane z direction. This is particularly obvious at the bottom of the melting pool, as shown in Figure 11. After the laser scan is over, tensile residual stresses build up along x and y directions to compensate the plastic compression. This effect is more relevant in the centre of the representative volume, however constraints are present between the regions close to the top surface and the centre. This induces residual stresses along the x and y directions also close to the top surface, as shown in Figure 12. The different grain orientations cause the heterogeneity of the residual stress. Grains with slip systems favourably oriented for plastic slip during loads along the x, y and z directions are more likely to develop plastic deformation along the corresponding directions. The simulations show that the maximum Schmid factor along the three axes is a suitable quantity to predict the corre- sponding plastic deformation and consequent residual stress, as shown in Figures 8 and 10. The simulations show that, in the stiffness degradation method, the presence of regions with a small residual stiffness imposes a constraint on the regions near the melting pool, preventing more plastic deformation to take place. This is the origin of the underestimation of the plastic deformation and of the residual stress using the stiffness degradation method. By contrast, the element elimination and reactivation method predicts higher plastic deformation in selective laser melting simulations. This difference between the two methods must be taken into account in future studies. Regions with larger tensile residual stresses along specific directions are more likely to undergo fracture perpendicular to those directions. Therefore, the present simulations suggest that fracture perpendicular to the x and y axes is more probable because of the magnitude of the corresponding residual stresses σ xx and σ yy . These two directions have similar residual stresses, as shown in Figure 10. This is consistent with laser-induced microcracking [81], in which fracture is observed at the top surface perpendicular to the laser scan direction.

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