Microstructural modeling of anisotropic plasticity in large scale additively manufactured 316L stainless steel

In this paper, the anisotropic mechanical response triggered by specific microstructures encountered in wire + arc additively manufactured 316L stainless steels is analyzed. For this purpose, a representative volume element is generated, with an internal geometry that is based on an average periodic fusion zone shape. The large columnar grain geometry with different preferred directions within a single fusion zone is reconstructed with an anisotropic Voronoi algorithm. Grain orientations are extracted from experimental data of different regions within a single sample. Crystal plasticity finite element simulations are performed to predict the macroscopic yield behavior under various loading angles. It is observed that the anisotropic response varies significantly within a single sample, which is essentially caused by differences in the texture. A spatial correlation between the orientation of a grain and its location within the fusion zone is identified, allowing to distinguish clusters of grains with a similar orientation. As a result, strain localizes in a global shear band across the entire height of the fusion zone, which prominently influences the macroscopic response. A comparison of the results relative to a standard isotropic Voronoi grain geometry shows that the elongated geometry of the grains does not influence the macroscopic yield behavior significantly. The local crystallographic orientations are dominating instead, whereby the correlation between location and orientation of grains has to be taken into account.


Introduction
Additive manufacturing (AM) is a technique in which parts or products are manufactured by layer wise addition of raw material. This technology is also known as 3D printing or rapid prototyping. AM technology is increasingly attractive for industry, given its flexible geometric capabilities, its reduced design cycle time and the reduced number of spare parts required in stock (Ashley, 1995;Mazumder et al., 1999;Zhang et al., 2018). Many different AM processes have been developed in the past 30 years, which have been applied to various fields including aerospace (Angrish, 2014), electronics (Tavakoli et al., 2017), maritime (Kim et al., 2014), biomedical (Bose et al., 2018), and food (Lipton et al., 2015).
Various AM processes exist for the 3D printing of metallic components. Many of these techniques are categorized in either powder bed fusion or direct metal deposition processes, based on the feeding of raw material . Within the latter category, processes are further categorized based on their feedstock and energy source used to melt the material. This paper focuses on wire + arc additive manufacturing (WAAM), which is suitable for 3D printing of large-scale (∼ 1 to > 10 m) metallic components, due to its large metal feed rate (Ding et al., 2015;Williams et al., 2016). In this process, which is similar to multi-pass fusion welding, a metal wire is used as feedstock. Using an electric arc, locally, part of the workpiece is melted on which new metal material is deposited. While the electric torch moves, the melt pool solidifies, leaving behind a new bead of solid metal.
Another potential of metal AM stems from the layer-by-layer production of the product. At each time during the process, local heat input can be controlled. This can ultimately lead to controllable local microstructures and corresponding mechanical properties (DebRoy et al., 2018;Mazumder et al., 1999). To ensure structural integrity of the parts and to take optimal benefit from this local controllability, a thorough understanding of the microstructural characteristics, residual stresses and resulting mechanical properties is required. Multiple process-related factors control these aspects, including metallic composition and the complete thermal-mechanical history in each material point. Various empirical studies have been performed on the influence of process parameters on structure and accompanying mechanical properties (Ding et al., 2015;Lu et al., 2017;Shamsaei et al., 2015;Wang et al., 2013). A few studies have tried to quantitatively capture these effects in computational models (Ahmadi et al., 2016;Biswas et al., 2019;Taheri Andani et al., 2018). This work aims to quantitatively model the effect of a WAAM-specific microstructure on the mechanical properties.
In order to unravel the complex microstructure-property relationship of WAAM produced products, a micromechanical model is constructed to predict the mechanical response (Francois et al., 2017). Therefore, a microstructural volume element is generated that is representative for a WAAM microstructure. For that purpose, grain geometry and texture are taken into account. During the solidification of the melt pool, crystals tend to grow in the direction of the local temperature gradient, leading to large columnar grains (Martina et al., 2015;Wang et al., 2013). Not only the grain shape is driven by the thermal history, also the grain orientation is strongly affected. This work focuses on AM-processed 316L stainless steel, consisting of a face centered cubic crystal lattice. In these metals, the 100 family of crystal directions is the so-called easy-growth direction (DebRoy et al., 2018). This is an energetically favorable crystal orientation in which grains grow fastest. This crystal direction will therefore align with the local temperature gradient, leading to an anisotropic texture. The anisotropy due to both the grain geometry and texture are taken into account upon construction of the computational microstructure.
Plastic deformation in a metal is governed by the movement of dislocations in the crystallographic slip systems. In order to model the constitutive behavior of metals at the scale of individual grains, a crystal plasticity model is an appropriate choice (Asaro, 1983). In this framework, individual atoms are not modeled, but the effect of the crystal stacking on plastic deformation is taken into account via the slip systems. The crystal plasticity finite element framework is here exploited to predict the anisotropic yield behavior of WAAM produced 316L stainless steel parts.
The paper is organized as follows. In section 2, the WAAM microstructure of a 316L stainless steel sample is characterized. From this analysis, a microstructural model is constructed in section 3. Here, grain geometry and texture are taken into account. In section 4, the crystal plasticity constitutive model is presented. Next, simulations are carried out, from which the anisotropic yield behavior is characterized in section 5. The representativeness of the microstructural volume element and the influences of grain geometry and spatial arrangement of grain orientations are discussed. Finally, conclusions are drawn in section 6.
In the following, scalars, vectors, second-order tensors and fourth-order tensors are denoted by a (or A), a, A and 4 A, respectively. Using Einstein summation convention, standard tensor operations are denoted as follows: dyadic product, a b = a i b j e i e j ; dot product, A · B = A ik B kj e i e j ; double contraction, A : B = A ij B ji , where e l (with l = 1, 2, 3) are the Cartesian basis vectors. Furthermore, the transpose operation (•) T of a second-order tensor is defined as A = A ij e i e j , A T = A ji e i e j . Finally, the material time derivative is denoted as (•).

Microstructural characterization
A 316L stainless steel sample is printed on a substrate of the same material using a robotic gas metal arc welding process. An electric arc with a current of 160 A and a voltage of 20.5 V is thereby used. The sample consists of 10 layers with 10 beads each, with overall dimensions 150 × 50 × 25 mm. Electron backscatter diffraction (EBSD) experiments are performed on parts of the 3D printed sample in order to characterize the shape and orientation of the microstructural grains. In figure 1a, an image of the sample is shown. Throughout this paper, the coordinate system indicated in this figure will be used: the deposition direction corresponds to the movement direction of the electric arc, the transverse printing direction marks the direction in which neighboring beads in a single layer are located and the building direction is the direction in which layers are stacked upon each other. Two cross sections from the sample are employed to obtain the EBSD data; one in the building-deposition plane and one in the building-transverse plane. Notice that the former is cut through the center of the fusion zone. An image of the latter is shown in figure 1b. Here, the fusion interfaces can be identified easily. Clearly, a very repetitive fusion interface pattern is produced by the process. A single fusion zone is indicated by the black dashed line as well. An averaged shape of such fusion zone will be utilized later to generate a representative model of the microstructure. In figure 2, the results of both EBSD scans are depicted. First, the grain geometry is considered. Both figures reveal that the grains are columnar and relatively large (∼ 50 to 100 µm along the shortest axis). Consequently, grain size effects, including tension-compression asymmetry, are negligible (El Shawish and Cizelj, 2017) and a standard crystal plasticity model can be used to describe the mechanical behavior of the grains. In the EBSD scan parallel to the deposition direction (figure 2a), the preferred main shape direction of the grains is aligned with the building direction, which is in accordance with the temperature gradient for relatively large scan speeds (DebRoy et al., 2018). Also in the cross-section perpendicular to the deposition direction (figure 2b), a clear preferred grain geometry can be identified. In this figure, the dashed lines indicate the fusion zone boundaries. The long axes of the crystals are directed towards the center of this fusion zone, which is also in accordance with the local temperature gradient during solidification of the metal.
Besides grain geometry, texture is also one of the main global characteristics of a metal microstructure. 316L stainless steel consists of austenitic grains, with a face-centered cubic (FCC) crystal lattice. In figure 2, the colors represent the crystallographic directions that align with the building direction. In the EBSD data for the cross-section parallel to the deposition direction (figure 2a), the 100 family of crystallographic directions of the majority of the grains is aligned with the building direction. The 100 direction is the so-called easy-growth direction of the FCC crystal (DebRoy et al., 2018). During solidification, crystals grow from the bottom of the melt pool in the direction of the temperature gradient. These crystals may originate from new nuclei or from epitaxial growth. Crystals with their easy-growth direction aligned with the temperature gradient will overtake grains with a less preferable orientation, leading to a microstructure with almost all grains aligned with the building direction. Note that figure 2a is only a 2D slice of a complex 3D microstructure, which, in general, yields an incomplete view of the texture. This can indeed be observed in figure 2b, where the EBSD data for the cross-section perpendicular to the deposition direction is depicted. In this figure, the fusion zone boundaries have also been added for clarity. It can be observed that in the center of the fusion zone the 100 family of directions of the grains is again aligned with the building direction. At some distance from the fusion zone center, the grain orientation is different due to the difference in the local temperature gradient, which gives preference to grains with their easy-growth direction aligned with this local temperature gradient leading to spatial variations in the texture.

Microstructural model
In this section, the generation of the representative volume element (RVE) of the microstructure in the cross-section of a fusion zone is described. The average fusion zone shape that has been experimentally obtained by Palmeira Belotti et al. (in prep.) is used as a starting point. The grain geometry will be incorporated in the computational RVE and the acquisition of a representative set of orientations will be explained.

Microstructural geometry
The computational geometry is generated using an anisotropic Voronoi algorithm (van Nuland et al., 2021). Standard Voronoi algorithms cannot be applied due to the large columnar grains which are oriented in different directions within a single fusion zone. The anisotropic Voronoi algorithm considers elliptical (instead of circular) growth from seed points. This preferred growth direction is specified for each individual seed point. The size of the growth vector determines the aspect ratio of the growth ellipse, which increases with the distance from the center of the fusion zone as observed in the experimental data. Inserting N g seed points, resulting in N g grains in the region of interest and specifying the corresponding preferred directions aligned towards the center of the fusion zone (see figure 3a), a representative geometry is obtained as depicted in figure 3b.
The shape of the fusion zone is periodic with matching boundaries as indicated by colors in figure 3. A periodic RVE emerges by considering a rectangular domain effectively containing one fusion zone. The rectangular periodic RVE, as shown in figure 3c, has dimensions 3.85 × 2.34 mm. Note that due to this choice of RVE, epitaxial growth is not taken into account in terms of the grain geometry, i.e. the fusion interfaces always act as grain boundaries. The depicted microstructural geometry consists of N g = 60 grains which are representative for the experimental data, as clarified next.

Representative orientation set
Next to the grain geometry, grain orientations that are representative for the microstructure need to be assigned as well. The complete set of orientations obtained from the EBSD scan (figure 2b) is depicted in the equal area pole figures shown in figure 4a. A set of N g orientations is extracted from the experimental data that adequately represents the complete texture of the 3D printed microstructure. This is accomplished by making use of an iterative procedure. First, grains are identified from the EBSD scan using a misorientation threshold. The average orientation per grain is determined and this list of orientations is ordered based on the size of the corresponding grain. The average orientations of the largest N g grains is used as an initial guess. Then the next orientation (with number N g + 1) is added to this set. Thereafter, the error of all possible sets of N g orientations out of the pool of N g + 1 orientations is evaluated with respect to the original EBSD data, and is defined as the L2 norm of the difference in orientation density at 1000 uniformly distributed random orientations. The set of N g orientations with the smallest error is used in the next iteration where the orientation with number N g + 2 is added. This process is considered to have reached convergence when an orientation set is obtained that has not changed during the last 20 iterations. The pole figures of the resulting representative set are shown in figure 4b. By construction, the representative set shows a good match with the experimental data.

Assignment of orientations to grains
Once the representative set of orientations is acquired, the individual orientations are assigned to specific grains within the computational microstructure. In the experimental data (see figure 2b), there is a clear correlation between the main shape direction of each grain (the geometrical long axis) and the easy growth direction of the crystal lattice (which is the 100 family of directions). The match between an individual orientation and a specific grain is quantified based on the projection of the easy growth directions of that crystal orientation onto the main grain shape direction. The orientation with the best match is assigned to the grain. Priority is given to the grains with largest size, since they have the largest contribution to the final texture. The resulting RVE is visualized in the top row of figure 5. In the bottom row of this figure, the corresponding experimental EBSD data is visualized. Again, a good agreement between the RVE and experimental EBSD data has been achieved, both in terms of the geometry and the texture.
In order to study the representativeness of the generated microstructural volume element (section 5.1), three other EBSD scans of different locations within the same sample are also used for the generation of the representative texture. Here, the same iterative procedure for the generation of the orientation set and the same algorithm to assign specific orientations to individual grains is adopted.

Microstructural model variants
To unravel the main microstructural features that affect the anisotropic plastic response, multiple cases are studied in which the spatial arrangement of grain orientations and the geometry of the grains are varied. To examine the influence of the correlation between the main shape direction of each grain and its orientation, results are compared to the case in which orientations are randomly assigned to grains. Also the influence of the geometry of the grains is investigated by comparing results to a microstructure with grains not having a preferred growth direction. For the latter, a standard Voronoi microstructure is employed as a representative geometry. The same seed point locations as in the fusion zone RVE are used, but the internal fusion interfaces are not considered. The microstructure is obtained with an isotropic Voronoi algorithm and is imposed to be periodic by periodically arranging the Voronoi seed points. Additionally, the orientations are assigned to grains in a specific manner. Since the same seed point locations are used as in the fusion zone geometry, the RVE is generated in the following manner: the orientation that belongs to a specific grain in the fusion zone RVE is assigned to the grain which originates from the same seed point in the Voronoi geometry. In this manner, solely the effect of the shape of the grains can be studied, without altering the spatial locations of the grains within the RVE.

Simulation setup
Simulations are performed on the microstructural RVE by means of finite element analyses. Therefore, the 2D geometry is discretized by 9279 (quadratic) serendipity elements with reduced integration (i.e. 4 Gauss points) and 76 quadratic triangular elements with full integration (i.e. 3 Gauss points). In the outof-plane direction (the deposition direction), a plane strain assumption is adopted that reflects the large dimensions of the fusion zone in the deposition direction. Due to the periodicity of the microstructure, periodic boundary conditions are applied on the RVE in the two in-plane directions. Therefore, for every pair of nodes on the top-bottom and right-left boundaries, the periodic boundary conditions can be written as Here, u T , u B , u R and u L are the displacements of nodes on the top (∂Ω T 0 ), bottom (∂Ω B 0 ), right (∂Ω R 0 ) and left (∂Ω L 0 ) boundaries, respectively, as depicted in figure 6. Likewise, u (p) are the displacements of the prescribed corner nodes p = 1, 2, 4, see figure 6. A macroscopic tensile strain of ε = 2.5 % is applied, representing an average in-plane uniaxial macroscopic stress state. Notice that the stress state in the out-of-plane (deposition) direction will be non-zero due to the plane strain condition. A strain rate oḟ ε = 5 · 10 −5 s −1 is used, corresponding to the strain rate of the single crystal experiments from which the material parameters have been identified (Karaman et al., 2001), see also section 4.3. Multiple loadings under different angles are considered in order to investigate the plastic anisotropy. The algorithm for applying macroscopic uniaxial loading conditions under an arbitrary angle, in combination with periodic boundary conditions, is elaborated in appendix A. The coupling between the microscale and the macroscopic quantities is obtained by volume averaging according to classical homogenization theory (Kouznetsova et al., 2001). The macroscopic first Piola-Kirchoff stress tensor P is obtained according to in which V 0 is the RVE microscopic volume in the undeformed configuration, f (p) are the external forces on the prescribed corner nodes p and x (p) 0 are the corresponding position vectors in the undeformed configuration.

Crystal plasticity
Numerical simulations are performed on the generated RVE. A crystal plasticity model is used to describe the constitutive behavior of individual grains at the scale of interest. In this section, the governing equations will be discussed (Asaro, 1983;Kalidindi et al., 1992). This framework has been implemented in the hypela2 subroutine for MSC.Marc/Mentat.

Kinematics
First, the total deformation gradient tensor F is multiplicatively split into an elastic and a plastic part, marked by the subscripts "e" and "p", respectively, reading In the crystal plasticity model, plastic deformation occurs via slip on the individual slip systems. Consider slip system α that is characterized by its slip plane normal n α 0 and slip direction s α 0 , defined with respect to the initial configuration. Both vectors have unit length. The total plastic velocity gradient tensor L p due to the presence of N s slip systems, each having a shear rateγ α , is then given by with the (non-symmetric) Schmid tensor P α 0 = s α 0 n α 0 defined with respect to the intermediate configuration. In this case of an FCC crystal lattice, the {111} 110 family of slip systems (i.e. N s = 12) is considered.

Constitutive relations
The elastic response of a crystal is described in terms of the 2nd Piola-Kirchhoff stress tensor in the intermediate stress-free configuration S e = JF −1 e · σ · F −T e , with J = det (F ) and σ the Cauchy stress tensor, as in which E e = 1 2 (C e − I) is the elastic Green-Lagrange strain tensor and C e = F T e · F e the elastic right Cauchy-Green deformation tensor and I the unit tensor. Furthermore, cubic symmetry is used for the fourth-order elasticity tensor 4 C.
The plastic deformation is governed by the resolved shear stress τ α on each slip system: The relation between the shear rateγ α and the resolved shear stress τ α is given by a power law relation: Here,γ 0 is the initial slip rate, g α is the slip resistance and m is the strain-rate sensitivity which is often taken very small, such that it closely resembles a rate-independent yield criterion. Hardening of the material is taken into account by defining an evolution equation for the slip resistance according tȯ The PAN hardening law (Peirce et al., 1983) is employed to model the strain hardening behavior of 316L stainless steel (Amelirad and Assempour, 2019;El Shawish and Cizelj, 2017). Therefore, the hardening coefficients h αβ are defined as in which δ αβ is the Kronecker delta, h 0 is the hardening modulus, g 0 the initial slip resistance, g ∞ the saturation slip resistance, q the latent hardening ratio differentiating between self-hardening and latent hardening and Γ the total amount of accumulated slip defined by

Crystal plasticity parameters
The parameters of this constitutive model for 316L stainless steel are obtained from literature (El Shawish and Cizelj, 2017;Gonzalez et al., 2012), and are listed in table 1. El Shawish and Cizelj (2017) identified the plasticity parameters from fitting crystal plasticity simulations to single crystal experimental data (Karaman et al., 2001;Kashyap and Tangri, 1995).

Results & Discussion
In order to investigate the anisotropy of the yield stress of the bulk material of a 3D printed 316L stainless steel, multiple loadings under different angles are considered. The results of these simulations are depicted in figure 7a. Here, the first Piola-Kirchhoff normal stress (engineering stress) component P is plotted against the macroscopic engineering strain ε. Both quantities are taken in the direction of the loading. The macroscopic yield stress P Y has been extracted from the results as the macroscopic stress P at 0.2 % macroscopic plastic strain. These values have been marked in the polar figure 7b. In both figures, the yield stress anisotropy can be identified, where the yield stress is maximum for loading under 30 • and minimum for loading under 45 • with respect to the transverse direction. This reveals a difference in yield stress of 12 MPa.

Representativeness of the microstructure
Up to this point, only data from a single EBSD scan has been used as input for the representative texture. In order to assess the representativeness of this microstructure, three other EBSD scans of the  figure 8a, labeled as texture I, II, III and IV. From these, four representative orientation sets have been extracted and these are again assigned to specific grains by shape correlation. The resulting generated RVEs are depicted in figure 8b. Note that texture I corresponds to the experimental data and RVE which have been considered before. As can also be observed, the computational geometry of the grains is not altered, solely the texture is different. In addition to the previous analysis, simulations are also performed for loading directions between 90 • -180 • . The results are shown in figure 8c. Note that all line colors correspond with the arrow color in between figures 8a, 8b and 8c. It can be seen that for all cases, loading under 0 • -90 • gives a similar response as loading under 90 • -180 • . The macroscopic yield stress is thus rotation symmetric over an angle of 90 • . The geometry of the grains does not possess this symmetry, suggesting that the geometry of the grains may not have a significant influence on the anisotropy of the yield stress. Also, textures I, II and IV have similar trends in the macroscopic yield stress; a maximum is around 30 • loading and a dip around 45 • . Texture III results in a different macroscopic response with a larger anisotropy; the difference between the maximum and minimum yield stress is 41 MPa. Since in all four cases the geometry of the grains is equal, it can be concluded that the texture is of great importance for the macroscopic yield behavior. It should also be noticed that these textures are based on experimental EBSD data from a single 3D printed sample, indicating that the macroscopic response can vary significantly from fusion zone to fusion zone. This characteristic is confirmed to be present in more WAAM printed samples, as for example described by Gordon et al. (2019). When the process can be controlled at such level that these variations are eliminated and different specific textures are induced, a large controlled spatial variation in macroscopic anisotropy can be reached. However, first a more thorough understanding of the cause of these variations is required.
In the following, focus is put on textures I and III (as defined in figure 8). In figure 9, pole figures of the RVE with texture I and III are presented. Here, each marker represents a grain direction from the 100 family (in figure 9a) or from the 111 family (in figure 9b). The color of the marker indicates the corresponding grain averaged Green Lagrange strain E = 1 2 (F T · F − I), from which the normal strain E in the direction of loading is taken. First, the results of texture I under 45 • loading are considered. The corresponding pole figure in figure 9b reveals that some of the grains remain nearly undeformed and other grains undergo large strains. This depends on the orientation of the grain. Two groups of orientations can be identified; the two horizontal dark blue bands and the two tilted dark red bands. These correspond to grains which have their preferred growth direction (a direction from the 100 family) aligned with the building direction or aligned under an angle of approximately 45 • with respect to the transverse direction, respectively, as can also be identified in the corresponding pole figure in figure 9a.
In figure 9, similar pole figures are also presented for various loading angles between 0 • -90 • . In texture I, upon rotation of the loading, the two different groups of grains will undergo a similar amount of deformation (e.g. around 15 • ) or one of the two groups will deform more than the other (e.g. around 0 • and 45 • ) due to preferential orientation of the slip systems. However in all cases, the amount of deformation for grains belonging to the same group is similar.
In texture III, a classification between grains can also be made. However, as evidenced by figure 9b when comparing the pole figures corresponding to 45 • loading for textures I and III, texture III contains far fewer grains that belong to the group that has its preferential growth direction under approximately 45 • with respect to the transverse direction. Grains that belong to this part of the texture are actually deformed largely under loading of 45 • . When these grains do not exist, other grains, with less preferred slip systems will plastically deform at higher stress levels, leading to a higher yield stress. On the other hand, around 0 • and 90 • loading, in texture III, there are more favorably oriented grains (i.e. grains with the 100 direction aligned with the building direction), leading to a lower yield stress under these angles, which is in correspondence with the results depicted in figure 8c. Thus, there is more preferred texture present in III than in I. Consequently, a different macroscopic yield behavior and a larger anisotropy result.
As discussed before, there is a correlation between the orientation of the grains and their location within the fusion zone. The two different groups of grain orientations originate from different temperature gradient directions within the melt pool. On the left side of the RVE (which represents the center of the fusion zone), grains are growing in the building direction. However, on the right side of the RVE (which represents the left side of the fusion zone), grains are growing more to the right. Here most of the grains with their preferred growth direction under approximately 45 • are located. Thus there is a spatial correlation of preferred orientations inside the RVE, which will be referred to as a grain cluster. Since there is also a relation between the orientation of the grain and its deformation, localization is prone to happen. In order to further elaborate on this, in figure 10, the equivalent Green Lagrange strain E eq = 2 /3 E : E in the microstructure is visualized at the end of the simulation. Various loading angles are also indicated. loading, strain localizes in a slip band across multiple grains over the full height of the right side of the RVE. This is the side in which most of the grains have their preferential growth direction under approximately 45 • with respect to the transverse direction. As a consequence of the formation of a large slip band across the full height of the RVE, a large part of the microstructure remains undeformed. Therefore, a lower stress is needed to reach the applied macroscopic deformation. This causes the dip in yield stress around 45 • loading.

Influential factors
In order to identify the main microstructural characteristics that govern the anisotropy of the macroscopic yield stress of 3D printed austenitic stainless steel, multiple cases are examined in which the spatial arrangement of grain orientations and the geometry of the grains are varied. This is done for textures I and III. All cases are depicted in figure 11a and the results of this analysis are shown in figure 11b for texture I and in figure 11c for texture III. In both figures, the blue solid line corresponds to the results of the RVE considered until now, i.e. the microstructural geometry described by the specific fusion zone shape with grains having a geometry obtained from their preferred growth direction. The orientations are obtained from a representative set and the locations of the orientations are correlated to the corresponding grain shape, as described in section 3.3. This case is referred to as "correlated, fusion zone", since the spatial arrangement of grain orientations is correlated to the shape of the grains and the geometry is the fusion zone specific grain geometry.

Effect of spatial arrangement of grain orientations
In order to identify the effect of the preferred location of all orientations on the anisotropy of the macroscopic yield stress, the following test case is performed: all 60 grains are classified based on their preferred growth direction in the two clusters that correspond to the two clusters of orientations that have been identified before. Then, all orientations are shuffled within their own cluster and assigned to a random grain corresponding to the same cluster. This preserves the clustering of similar orientations, but the direct neighbors of the grains are different. This test case is referred to as "cluster shuffled, fusion zone", since the orientations are only shuffled within their own clusters and the geometry is still fusion zone specific. Note that due to the fact that all grains do not have exactly the same size, a slightly different texture results. However, as shown later, this does not influence the macroscopic yield stress significantly. The results of the "cluster shuffled, fusion zone" analysis are depicted by the red solid line in figures 11b and 11c. For both cases, there is an influence on the exact values of the yield stress indicating that the misorientation distribution between neighboring grains has a visible effect on the anisotropy of the macroscopic yield stress. However, all features and trends that have been identified in the "correlated, fusion zone" case are also present in the "cluster shuffled, fusion zone" case, including the large anisotropy in the case of texture III. Furthermore, the dip in the yield stress for loading under 45 • in the case of texture I can also be recognized. This is due to the fact that the clustering of orientations is still present, which allows for a large shear band to appear over the full height of the RVE around this specific loading direction.
Next, the clustering of the orientations is completely eliminated by assigning the same set of 60 orientations to the grains randomly. Here, the clusters of orientations are not taken into account. This case is referred to as "fully shuffled, fusion zone", since the texture is fully shuffled and the geometry is still fusion zone specific. This analysis is performed for both textures I and III and the results are depicted by the yellow solid line in figures 11b and 11c, respectively. Figure 11b reveals that the dip around 45 • loading has disappeared. It is no longer possible to induce a shear band through multiple grains in a large area of the RVE. Furthermore, figure 11b shows that the yield stress decreased for loadings under approximately 0 • and 90 • with respect to the "correlated, fusion zone" case. In figure 10, for loading under these angles, a large global shear band under approximately 45 • is active. In order to get such a band in a periodic setting, it needs to pass through grains which have a less favorable orientation (i.e. the grains from the group which have their preferred growth direction around 45 • with respect to the transverse direction). This effectively increases the macroscopic yield stress. Upon shuffling all orientations, there are many small shear bands arising as depicted in figure 12, where the equivalent Green Lagrange strain E eq is shown at the end of the simulation of the "fully shuffled, fusion zone" RVE with texture I under 0 • loading. These smaller shear bands do not have to cross large areas of grains with less preferred orientations, which effectively decreases the yield stress. Here, again, the clustering of the orientations plays an important role. The spatial location of the oriented grains in combination with the preferred texture that has been identified in wire + arc additively manufactured 316L stainless steel thus has a prominent influence on the anisotropy of the macroscopic yield stress.

Effect of geometry
Also the influence of the geometry of the grains is investigated by comparing the results to a microstructure with grains having no preferred growth direction. Therefore, RVEs are generated using a standard Voronoi microstructure. As discussed in section 3.4, the same orientations and "shuffling" are applied in all cases; only the grain geometry is different. This is done for all three microstructures and for both textures I and III. This leads to the six Voronoi RVEs depicted in figure 11a.
First, the results of the "correlated, Voronoi" case are considered. Here, the orientations of the grains are correlated to their location and the Voronoi grain geometry is used. The results of this analysis are depicted by the dashed blue line in figures 11b and 11c for textures I and III, respectively. For both textures, there is no significant influence of the grain geometry on the macroscopic yield behavior, compared to the most realistic "correlated, fusion zone" case. Note that the size of the grains in the "correlated, Voronoi" case are not exactly equal to the size of the corresponding grains in the "correlated, fusion zone" case, leading to a slightly different texture. This also indicates that this effect does not influence the macroscopic yield stress significantly. In practice, the shape of grains cannot be influenced without simultaneously altering the texture, due to the aforementioned correlation between the main shape direction of the grain and the easy-growth direction of the crystal lattice. The insensitivity to the grain geometry is therefore particularly interesting for microstructural modeling, since a simple grain geometry suffices in predicting plastic anisotropy.
In order to further investigate the influence of the geometry of the grains, the results of the "cluster shuffled, Voronoi" and "fully shuffled, Voronoi" case are depicted by the red dashed line and yellow dashed line, respectively, in figures 11b (for texture I) and 11c (for texture III). In all cases, the macroscopic anisotropy predicted by the simple Voronoi structure is similar to the corresponding fusion zone grain geometry. Some differences can be recognized when comparing the "fully shuffled, Voronoi" case with the "fully shuffled, fusion zone" case for texture I. Under a loading of approximately 45 • , a somewhat larger yield stress is predicted using the Voronoi grain geometry. This is due to the fact that in the fusion zone geometry, there are locations in which there are only 2-3 grains across the height of the RVE. In this case, a full height shear band can be formed more easily, see figure 13a, where the equivalent Green Lagrange strain E eq is depicted at the end of the "fully shuffled, fusion zone" simulation with a loading under 45 • for the RVE with texture I. This slightly lowers the yield stress. In the Voronoi geometry, there are 5-8 grains across the height. Under 45 • loading there is no such pronounced slip band present, as can also be seen in figure 13b, which effectively increases the macroscopic yield stress. When the strain maps in figure 13 are compared to the strain maps for the "correlated, fusion zone" depicted in figure 10, it can be observed again that global strain patterns across multiple grains or even the entire RVE can be recognized in the response of the 3D printed 316L stainless steel samples. This is a result of the aforementioned correlation between the orientation and the location of the grains within the fusion zone. As discussed before, when this relation is taken into account in the modeling of the microstructure, the geometry of the grains does not have a significant influence on the macroscopic yield behavior. A simple grain geometry thus suffices for predicting the plastic anisotropy. However, it has been shown that the clustering of orientations and the misorientation distribution between neighboring grains has a visible effect on the anisotropy of the macroscopic yield stress. Therefore, if a simple grain geometry is used it remains necessary to correlate the location of oriented grains and the misorientation distribution with the experimental data, by using for example an algorithm as described by Biswas et al. (2019).

Conclusions
In this paper, the key features of the microstructure of wire + arc additively manufactured 316L stainless steel metals have been modeled. For this purpose, the microstructure has been characterized by means of electron backscatter diffraction (EBSD) analysis from which the grain geometry and texture has been extracted. This provides the required input for the generation of a representative volume element (RVE), which has the size of a single fusion zone due to microstructural periodicity. Internal fusion zone boundaries have been taken into account in the geometry of the RVE by making use of an average periodic fusion zone shape. The grain geometry has been obtained by employing an anisotropic Voronoi algorithm in order to take into account the large columnar grains with different preferred directions within a single fusion zone. The orientations of the grains have been obtained from a representative orientation set which has been extracted from the EBSD data. The orientations are assigned to individual grains based on a correlation between the main shape direction of each grain and the easy-growth direction of the crystal.
Numerical simulations have been performed in a plane perpendicular to the deposition direction of the electric arc. The mechanical anisotropic behavior of individual grains has been modeled by making use of a crystal plasticity framework. Periodic boundary conditions have been applied in order to obtain the bulk mechanical response of the material. Therefore, boundary and edge effects have not been considered. EBSD data from four different fusion zones within a single sample have been used as input for the generation of the RVEs. From the results of the numerical simulations, the macroscopic yield stress has been extracted under various loading angles in order to characterize the anisotropic response of the macroscopic yield behavior. Based on the response of individual grains, two groups of orientations have been identified corresponding to grains having their preferred growth direction aligned with the building direction or under an angle of approximately 45 • with respect to this direction. One out of the four aforementioned fusion zones is missing one of these groups of orientations, resulting in a significantly larger anisotropic behavior. This also indicates that the texture is of great importance for the macroscopic response, which can vary significantly within a single product. This also implies that when the process can be controlled at a level that eliminates these variations and specific textures can be obtained, a large controlled spatial variation in macroscopic anisotropy may be at reach. Furthermore, modeling multiple fusion zones with different textures is necessary to predict the full macroscopic response if these variations cannot be avoided. The two different groups of orientations originate from different local temperature gradients within a single fusion zone. Therefore, a correlation between the orientation of a grain and its location within the fusion zone has also been identified, leading to clustering of orientations from the same group. As a result, strain can localize in a global slip band across the entire height of the RVE. Depending on the local texture variations from bead to bead, this can in turn prominently influence the macroscopic response. The location of the orientations in combination with the texture has a significant influence on the plastic anisotropy. In the current work, the anisotropic yield surface was not yet complete for a 3D setting. This may be done in the future by extending the model to three dimensions and including simulations with loading in the deposition direction. Furthermore, experimental validation of the model is not possible at present due to the lack of the corresponding data.
The influence of the specific grain geometry has also been studied by comparing results to RVEs with a simple Voronoi grain geometry. From this analysis it is concluded that the elongated geometry of the grains does not influence the macroscopic yield behavior significantly. In practice, the shape of the grains cannot be influenced without simultaneously altering the texture, due to a preferred growth correlation of the grains. This conclusion is therefore mostly interesting for microstructural modeling since a simple grain geometry suffices in predicting the macroscopic yield response. However, as discussed before, it should be taken into account that the location of the orientations does have a pronounced effect on the macroscopic anisotropy.

A. Macroscopic uniaxial loading under an angle in combination with periodic boundary conditions
According to classical computational homogenization theories (Kouznetsova et al., 2001), when applying periodic boundary conditions (PBCs) at the microscopic scale, the prescribed displacement u (p) of corner node p at the microscopic scale (as indicated in figure 14a) can be related to the macroscopic deformation gradient F M according to in which x (p) 0 are the coordinates of the corner nodes p in the undeformed configuration and I is the unit tensor. The corresponding displacement components are then applied as prescribed displacement boundary conditions (BCs) on the individual degrees of freedom (DOFs) related to the finite element (FE) Cartesian coordinate system { e x , e y }, as depicted in figure 14a.
Equation (12) can be used to prescribe a fully known deformation gradient tensor. However, in case of uniaxial loading conditions, the RVE is required to be deformation controlled only in the loading direction and stress-free in the direction perpendicular to the loading direction. Assume that uniaxial loading conditions are applied in a different Cartesian coordinate system e x , e y and the uniaxial loading direction is aligned with e x . The required BCs are then written as e x · F M · e x = λ x = 1 + ε, e y · σ M · e y = e x · σ M · e y = e y · σ M · e x = 0, in which σ M is the macroscopic (i.e. volume averaged) Cauchy stress tensor. In order to meet these BCs, first, on corner node 1, prescribed displacement BCs are applied to both DOFs in order to suppress rigid body translations. Then, the displacement of corner nodes 2 and 4 can be written as u (p) = u x (p) e x + u y (p) e y , p = 2, 4, with u y (p) unknown. Moreover, under the requirement that corner node 1 is located at the origin (i.e. x (1) 0 = 0) and the macroscopic deformation gradient tensor F M is of the form F M = λ x e x e x + λ y e y e y + λ yx e y e x , where λ y and λ yx are unknown, the applied displacement in the direction of the loading u x (p) is obtained by combining equations (12) and (15) as = e x · (F M − I) · x (p) 0 , p = 2, 4.
Furthermore, in order to fulfill condition (14), traction free BCs are applied to corner nodes 2 and 4 in the direction of e y . Now, when the uniaxial loading direction is aligned with the FE coordinate system, i.e. e x = e x and e y = e y , the prescribed displacement BC (19) and traction-free BC are applied directly to the corresponding DOFs. The resulting BCs are also indicated in figure 14b.
In case macroscopic uniaxial loading conditions are applied in a direction that is not aligned with the FE coordinate system, a more elaborate application of the BCs is needed. A schematic representation of this situation is depicted in figure 14c. As a starting point, the relation between loading direction e x and the FE coordinate system { e x , e y } is given as e x = a x e x + a y e y .
Note that a x and a y are dependent variables. For example, when e x , e y is rotated over an angle θ with respect to { e x , e y }, then a x = cos (θ) and a y = sin (θ). The displacement of corner nodes 2 and 4 can also be written with respect to the FE coordinate system as u (p) = u (p) x e x + u (p) y e y , p = 2, 4, in which u (p) x and u (p) y are the corresponding corner node DOFs. The prescribed displacement BCs (19) need to be applied with respect to the FE coordinate system, which are related to the corresponding DOFs u (p) x and u (p) y of the corner nodes by substitution of relations (20) and (21) into equation (17) as u x (p) = a x u (p) x + a y u (p) y , p = 2, 4.
Relations (19) and (22) describe the prescribed displacement BCs that are to be applied. These can be implemented in a FE package by, for example, assigning u x (p) to the DOFs of an auxiliary node which is not connected to any element in the mesh. The position of this auxiliary node is arbitrary. Using this method, relation (19) is directly applied as a standard prescribed displacement BC to this auxiliary node. Furthermore, (22) is rewritten in the form of tying relations as also used for the application of PBCs.