Influence of Pore Characteristics on Anisotropic Mechanical Behavior of Laser Powder Bed Fusion–Manufactured Metal by Micromechanical Modeling

In recent times, additive manufacturing (AM) has proven to be an indispensable technique for processing complex 3D parts because of the versatility and ease of fabrication it offers. However, the generated microstructures show a high degree of complexity due to the complex solidification process of the melt pool. In this study, micromechanical modeling is applied to gain deeper insight into the influence of defects on plasticity and damage of 316L stainless steel specimens produced by a laser powder bed fusion (L‐PBF) process. With the statistical data obtained from microstructure characterization, the complex AM microstructures are modeled by a synthetic microstructure generation tool. A damage model in combination with an element deletion technique is implemented into a nonlocal crystal plasticity model to describe anisotropic mechanical behavior, including damage evolution. The element deletion technique is applied to effectively model the growth and coalescence of microstructural pores as described by a damage parameter. Numerical simulations show that the shape of the pores not only affects the yielding and hardening behavior but also influences the porosity evolution itself.


Introduction
The potential of additive manufacturing (AM) has rapidly grown as it allows for the fabrication of complex 3D parts with minimal effort. All the available AM techniques are based on an incremental layer-by-layer building approach, which gives them an edge over conventional manufacturing techniques. Several studies have focused on the AM processing of metallic components, and laser powder bed fusion (L-PBF) and electron beam powder bed fusion (E-PBF) have drawn increased attention in this regard. [1,2] A focused laser in L-PBF or an electron beam in E-PBF is used as the heating source to selectively melt the layer of powder deposited on the build platform taking into account a computer-generated geometry file for the structural information and a machines file in which the process parameters are listed. Regardless of the utilized heating source, the general approach of these methods is first, a 3D Computer Aided Design (CAD) model is created, and the model is then sliced into several layers of the required thickness for deposition. These data are used as input for melting and building the component using the heat source. The current work focuses on austenitic stainless steel (AISI 316L), a material that has excellent processing properties with the L-PBF process. [3][4][5][6] During L-PBF processing, a fast heat-up of the deposited 316L metallic powder layer and a rapid solidification of the melt pool occur. Also, depending on the laser intensity, constant reheating and recooling of the previously melted layers take place. Due to this complex solidification process, the resultant microstructures contain elongated as well as equiaxed grains with directiondependent material properties. Figure 1a shows the electron back scatter diffraction (EBSD) map and the corresponding orientation triangle for the L-PBF-manufactured 316L sample. The L-PBF process parameters such as laser power, exposure time, point distance, and hatch spacing can be changed to control the energy imparted to the build and thereby controlling the size of the melt area. They have been widely studied with regard to their influence on the generation of pores and microcracks within the microstructure. [7,8] The formation of pores is determined by the particle's flowability and the way in which the powder is fed. Figure 1b shows the pores generated within the microstructure of the L-PBF-manufactured 316L sample. As these microstructural defects determine the effective mechanical properties of the material, the influence of AM process parameters on its resultant microstructure and the microstructuremechanical property interdependency are topics of active research. [3,9] Many experimental studies [3,[10][11][12] have analyzed this important link between AM process parameters, AM microstructures, and the resultant material properties.
From this perspective, the micromechanical modeling approach, which predicts the mechanical behavior by using synthetic microstructure simulations, serves as an efficient tool to further establish the interdependency. This approach explicitly takes into account the important microstructural features such as grain size and shape distribution functions along with orientation and misorientation distribution functions (ODF and MDF) to model synthetic microstructures. Crystal plasticity (CP) models have been successfully used to account for important micromechanical features and capture complex material deformation behaviors. [13] In the context of AM, CP models have been used to understand the influence of microstructural level properties on the macroscopic mechanical response of the materials. A physicsbased CP model has been used by Motaman et al. [14] to investigate the impact of AM process-induced morphological and texture heterogeneities on the macroscopic strain hardening behaviors. By isolating the two heterogeneities, they have shown that grain morphology has a direct impact on the anisotropic plasticity, whereas texture independently influences the curvature of the strain hardening curve. In another similar work, Biswas et al. [15] have shown that the effect of grain shape on the macroscopic flow behavior is small in comparison to that of texture. And they emphasize that the crystallographic orientation has a significant impact on the macroscopic flow behavior, whereas the misorientation angle influences only the local mechanical events. Other works, such as Andani et al. [16] and Ahmadi et al., [17] model explicitly the complex AM melt pool geometries, overlaps, and layer orientations and use CP models for capturing the effect of processing defects and textures on the mechanical response. A CP model in combination with in situ synchrotron X-ray diffraction experiments has also been used to study the impact of microscale residual stresses on the tensioncompression asymmetry of AM 316L stainless steel. [18] In the literature, CP models are often extended to include the effects of pore growth and coalescence. [19][20][21] In the works of Asim et al., [22,23] they have investigated the void growth in a dual-phase Ti alloy at the interface of dissimilar α-β crystals. Chen et al. [24] have analyzed the void growth on grain boundaries in Ni-based superalloys. In the recent work of Diehl et al., [20] the damage around a void was comprehensively studied by a coupled CP-phase field fracture approach. Also, a multiscale investigation was performed in the work of Shang et al. [25] to study the effect of grain size on void evolution and fracture in polycrystalline austenite steel 316LN.
As most of the available studies consider single-crystal materials, the current work focuses on additively manufactured polycrystalline material to gain a more in-depth insight into the influence of defects on anisotropy in plasticity, including damage evolution. Emphasis is laid on modeling explicitly the complex microstructural features such as realistic grain size, morphology, and texture. To model the complex synthetic microstructures, a novel collision-driven particle-packing approach is used. Pores of varying volume fractions and shapes are modeled by removing mesh elements from the generated synthetic microstructures. To describe the mechanical response of the polycrystals, a nonlocal CP model is used; [26] a linear damage evolution model based on plastic strains and the element deletion technique is used to model the growth and coalescence of pores as the material softens. Finally, the effects of pore volume fraction and pore shapes are investigated to understand the local mechanical behavior.

Materials and Used Powder
In this work, the austenitic stainless steel AISI 316 (face-centered cubic crystal structure) was processed using L-PBF. The powder was produced by gas atomization and has a particle size of 20-63 μm in the classified state. The chemical composition was measured by means of spark emission spectrometry on compressed samples and is listed in Table 1 together with the nominal composition. Due to the spherical morphology and the low satellite density, the powder has a good flowability of 13.93 s/50 g (measured using a Hall Flowmeter according to DIN EN ISO 4490, t ¼ 22.4 C, air humidity ¼ 50.3%) and a bulk density of 4.34 kg dm À3 (measured according to DIN EN ISO 697).

Sample Preparation
The L-PBF system Realizer SLM 100 was used for the sample setup. The system has a continuously operating ytterbium laser with a wavelength of 1065-1075 nm, a fixed laser diameter of 90 μm, and a maximum laser power of 100 W (effectively measured laser energy on the building platform level P eff ¼ 73.5 W). The samples with a size of 5 Â 5 Â 2.5 mm were produced on a construction platform of the same material with a diameter of 120 mm. During the construction job, the laser was moved in a meandering manner (x-y interlayer stagger approach) over the powder bed (h ¼ 80 μm). After the densification of one slice, a further powder layer (30 μm) was deposited, and the laser was moved in a perpendicular course. The powder was processed in an inert Ar atmosphere with an O 2 content of less than 0.3% by volume. Based on previous parameter studies from Röttger, [3] samples with a targeted porosity of 2 vol% were prepared for this study with a maximum laser power of 100 W, a point distance ¼ hatch distance of 100 μm and a exposure time of 250 μs, leading to a specific energy input of 0.234 J mm À3 . For tensile tests, corresponding tensile samples were made in the horizontal direction on the building platform with the same parameters.
The porosity of the samples was determined using the quantitative image analysis. For this purpose, micrographs were previously binarized with a light microscope and evaluated using the Aquinto A4i Analysis software. The EBSD technique was used to investigate the grain size and the crystallographic orientations of the grains. An EBSD detector type EDAX/TSL from TexSem Laboratories was used, which was connected to a LEO Gemini 1530 VP scanning electron microscope. For the EBSD investigations, an acceleration voltage of 21 keV and a working distance of 13 mm were used. During the EBSD measurements, the samples were tilted by 70 to the pole piece. The colors in the EBSD mappings represent the crystal directions in the surface normal direction (ND).

Micromechanical Modeling
A micromechanical model is developed to predict the plastic deformation behavior of realistic microstructures by means of homogenization techniques. The complex polycrystalline microstructure of the L-PBF-processed material is described by a representative volume element (RVE), which considers the key microstructural features such as grain morphologies, grain size distributions, and crystallographic orientations. Defects such as pores are modeled by deleting mesh elements of the RVE in the finite element calculation by using the Abaqus-Python scripting environment. [27] The nonlocal CP model and the damage model are implemented in a user-defined material subroutine (UMAT) to capture the stress-strain response of RVE. In the following subsections, different aspects of the modeling approach such as the RVE generation, the pore modeling, the damage model, and the nonlocal CP model parameterization are detailed. For the description of the nonlocal CP model, the reader is kindly referred to the Appendix.

Microstructure Modeling
To efficiently model the complex microstructures of L-PBFprocessed metals, detailed information regarding their grain and pore morphologies must be considered. For this, the statistical data of the microstructure are utilized, which are characterized by scanning electron microscopy (SEM) and light optical microscopy (LOM). The EBSD technique is used to quantify the grain size and the crystallographic orientations, whereas the porosity information is quantified by LOM micrographs. Both these techniques are described in Section 2.3.

RVE Generation
As observed in Figure 1, the morphologies of the grains are elongated along the building direction (BD). Here, BD refers to the direction along which the material deposition is done in AM. Furthermore, the two directions orthogonal to the BD are referred to as normal (ND) and transverse (TD) directions, respectively. To construct an RVE that consists of these complexly shaped grains, the synthetic microstructure generation tool Kanapy [28] is used. Kanapy approximates grains as ellipsoidal particles following a particular size distribution and packs them into a pre-defined simulation box representing the RVE. The general framework used in the packing process is the collision detection [29] and response system for particles under random motion in the simulation box. The interaction between particles can be modeled by breaking it down into stages of collision detection and response, and the interaction between the particles and the simulation box can be modeled by evaluating whether the particles cross the boundaries of the box. If periodicity is enabled, periodic images on the opposite boundaries of the box are created; otherwise, the particle position and velocity vectors have to be updated to mimic the bouncing-back effect. For collision detection, the algebraic separation condition developed by Wang et al. [30] is used to determine whether two static ellipsoids overlap. Consider two ellipsoids in their standard form A∶X T AX ¼ 0 and B∶X T BX ¼ 0 in R 3 , where X ¼ ½x, y, z, 1 T . According to Wang et al., [30] the characteristic equation is given as Equation (1) can be solved, and, depending on the nature of the roots obtained, the overlap or separation condition can be determined. Wang et al. [30] have established that Equation (1) has at least two negative roots and that, depending on the nature of the remaining two roots, the separation conditions are given: 1) A and B are separated if f ðλÞ ¼ 0 has two distinct positive roots; 2) A and B touch externally if f ðλÞ ¼ 0 has a positive double root; 3) A and B overlap for all other cases.
If collision is established, the particles should bounce back and continue their motion in opposite directions. The position vectors r A , r B are used to determine the angles ϕ and θ, which represent the in-plane and out-of-plane directions for the bounce-back response. The velocity vectors v A , v B are updated accordingly.
As each particle has to be checked for collision with every other particle in the domain, the approach can be computationally expensive. To overcome this, a two-layer collision detection scheme is implemented wherein the outer layer utilizes an octree spatial partitioning data structure [31,32] to estimate which particles should be checked for collision. The inner layer consists of a bounding spheres hierarchy, [33] which carries out the collision detection only if the bounding spheres between two particles overlap. Thus, the method is able to handle a large number of collision checks in an efficient way. For a more detailed description regarding the collision detection scheme and its implementation, refer to Prasad et al. [28] The user-defined simulation box size and the ellipsoid size distribution are used for creating the simulation box and ellipsoids. The simulation begins by randomly placing ellipsoids of null volume inside the box, and each ellipsoid is given a random velocity vector for movement. As the simulation proceeds, the ellipsoids grow in size along their axes and also collide with one another, thus updating their position and velocities. The simulation terminates as soon as all the ellipsoids have reached their defined volumes; the process is shown in Figure 2.
Once the ellipsoids are tightly packed with minimal acceptable overlaps, they are processed further for meshing. As Voronoi tessellation, which was used in the previous work for 2D smooth-boundary RVE generation, [34] cannot be applied to anisotropic particles such as ellipsoids, a voxel-based mesh-generating routine that utilizes a convex hull for discretizing the RVE is used now. An example of ellipsoidal packing and the corresponding voxelized RVE are shown in Figure 3. The generated RVE with hexahedral elements can be exported to Abaqus. [27] 3.1.2. Modeling Pores within RVE As observed from Figure 1, the pore morphologies can be approximated by spheres and ellipsoids. The pores are placed inside the RVE by deleting the mesh elements, and their placement is controlled by using Abaqus-Python scripts. The AM laser power and the scanning speed control the specific energy input, and it is reported in Röttger et al. [3] that with increase in specific energy input the porosity of L-PBF-manufactured specimens decreases. These parameters also influence the size and the shape of the generated pores. Thus, in this work, pore volume fraction and pore shape are the two aspects that are varied to generate different combinations of test cases. Table 2 lists the test cases considered in this regard. To study the influence of the pore volume fraction, three test cases (represented as "A" in Table 2) with 2%, 5%, and 8% volume fraction of pores are realized. These values are within the porosity range observed experimentally. [3] The pores in all these test cases are modeled with a mean pore size of 6.2 μm, as this is the mean pore size observed through LOM micrographs of light microscopy. To study the influence of the pore shape, three test cases (represented as "B" in Table 2) with 9.3 μm mean pore size at 2% and 8% porosity volume fraction are realized. As the pore shape is a difficult parameter to quantify from LOM micrographs, the equivalent diameter (d) is used for computing the major axis (a) and the www.advancedsciencenews.com www.aem-journal.com two minor axes (b, c) lengths of the elongated pores . The elongated pores approximated by ellipsoids have an aspect ratio of 3:1:1 (a∶b∶c) and are placed inside the RVE with their major axis oriented along the sample's BD. Figure 4 depicts the 3D view of RVEs containing equiaxed and elongated pores. A bigger mean pore size (d ¼ 9.3 μm) is chosen for case study B, because the computed minor axes lengths of the elongated pores (for d ¼ 6.2 μm) are smaller than the mesh size (3 μm).

Damage Model
It is well known that damage due to growth and coalescence of pores leads to ductile fracture in metals. Generally, CP models are coupled with damage models to investigate material failure.
Other works in the literature couple CP with phase-field method to model damage around a void [20] and fatigue crack initiation, and propagation. [35] CP-based crack-nucleation models for polycrystalline microstructures have been investigated in previous studies. [36,37] A criterion for the onset of damage can be used in combination with the damage evolution models for damage modelling. [38,39] The initiation criterion can be specified by providing the equivalent plastic strain at the onset of ductile damage as a function of stress triaxiality, Lode angle, and strain rate. [40,41] However, such a modeling approach, in general, features a macroscopic phenomenological nature, owning great predictive capability in application but also requiring a significant amount of effort for parameter identification by experiments at various stress states and loading conditions. Although such criteria can be used in combination with CP models, the data required are difficult to obtain, as experiments covering a range of stress triaxialities must be conducted at the microscopic scale. With the current state of the art in micromechanical characterization and testing, it is still a challenge to conduct all these tests, and therefore no systematic data sets in this regard have been published. In addition, access to such testing devices is still very limited to only a few labs. Assumptions can be made to upscale the damage initiation by a localization criterion, [42,43] but it is so far only valid for dual-phase materials and not yet generalized. Therefore, in this work, the damage   www.advancedsciencenews.com www.aem-journal.com initiation criterion is modeled by defining a lower limit for the equivalent plastic strain. The equivalent plastic strain (p) is the integral of the equivalent plastic strain rate (ṗ), which is computed using the Frobenius norm aṡ whereĖ p is the plastic strain rate tensor. The plastic strain tensor (E p ) is deviatoric and is computed using the plastic deformation gradient (F p ). The plastic deformation gradient is computed using the plastic velocity gradient (L p ), which depends on the plastic shear rate (γ α ) and the Schmidt tensor (M α ), given by Equation (7) in the Appendix. For values of the equivalent plastic strain below the lower limit (p l ), the damage variable equals zero (D ¼ 0). Once the value exceeds the lower limit, the damage variable evolves linearly, as shown in Figure 5. Materials undergoing damage show pronounced loss in yield strength and degradation of their stiffness. In numerical modeling, as the damage variable evolves, the stiffness components are deteriorated proportional to ð1 À DÞ and the stress tensor in the material is given by σ ¼ ð1 À DÞσ. The linear damage variable evolves as where p u is the upper limit of equivalent plastic strain that defines the state of complete damage. Once this value is reached, the damage variable is set to a maximum value (D ¼ 0.99 for numerical stability). And once the value of D at all the integration points of a severely deformed element reaches its maximum, that particular element is deleted from the FE mesh in Abaqus by using a status variable. [27] It is noted that the microscopic damage modeling formulation is inspired by a successful macroscopic concept, as demonstrated with various steels. [39,44] The damage initiation and linear evolution law is a simplification of the nonlinear damage process in physics. However, it is good practice to simulate the process and give an assessment indicator to distinguish damaged materials from damage-free ones.

Boundary Conditions and Homogenization
The RVE is modelled such that it is periodic in all directions. The periodic boundary conditions methodology introduced by Smit et al. [45] is implemented on the RVE to capture the deformations in a realistic manner. The displacements of the nodes on the opposite surfaces of the RVE boundary are coupled, and thus the strain is captured. With this coupling, the complete transfer of stress across the opposite surfaces is also achieved. For a detailed description of the implementation of periodic boundary conditions, please refer to Boeff. [46] The stress and strain tensors of the RVE under periodic boundary conditions are homogenized by fulfilling the macrohomogeneity condition. An overview of such averaging theorems is provided in Nemat-Nasser. [47] The RVE is simulated under the uniaxial tensile loading condition to obtain the flow curve.

RVE Convergence Study
In addition to modeling the complex grain morphologies, the size of the RVE has a vital influence on the mechanical behavior. Therefore, an RVE convergence study must be performed to evaluate the ideal RVE size for further simulations. Based on grain size and shape statistics obtained from EBSD, four RVEs with increasing number of grains 298, 520, 702, and 902 are generated by using Kanapy. The size of the RVE correlates directly with the number of grains in them, so the four RVEs shown in Figure 6 are of increasing side lengths. Their mechanical response is evaluated through numerical simulations, and they are analyzed for convergence based on the yield stress values (calculated using the 0.2% offset method). Table 3 provides the mesh element and the yield stress values for the four considered RVEs. Selecting the appropriate RVE depends not only on the convergence of the yield strength, but also on the number of mesh elements as this influences the computational costs. An RVE with 520 grains is appropriate for the current study as the yield stress value converges and the computational cost with respect to FEM simulation time is not high.

Parameterization of Nonlocal CP Model
The experimental tensile test curve obtained from Röttger et al. [3] as shown in Figure 7 is used for parameterizing the nonlocal CP model. As the experimental curve belongs to the sample with a negligible amount of pores, a nonporous RVE with 520 grains is used for the parameterization. The pole figures corresponding to the EBSD map (Figure 1a) are shown in Figure 8. It can be seen that for the L-PBF-manufactured sample, there is no preferred crystallographic orientation. Thus, the crystallographic orientations for the grains in the RVE are chosen such that the texture index is close to "1", representing a texture close to random texture. Nonlocal CP parameters are adapted to match the experimental result. From the flow curve comparison in Figure 5. Evolution of the damage parameter as a function of equivalent plastic strain. p l and p u represents the lower and upper limits for damage initiation and complete damage, respectively.
www.advancedsciencenews.com www.aem-journal.com Figure 7, the parameters summarized in Table 4 yield a good agreement with the experiment and are used in the rest of this study.

Results and Discussion
Microstructural features such as grain size, grain shape, and crystallographic orientation of the RVE are kept constant to isolate and to study the influence of different aspects of pores on the plastic deformation behavior. From the RVE convergence study, the nonporous RVE with 520 grains (refer Figure 6) is chosen, and the Abaqus-Python script is used to place pores inside the RVE according to the test cases described in Table 2. Note that the damage model described in Section 3.2 is used only for the pore shape study to understand the influence of pore shapes on local events such as damage.   Figure 7. Equivalent stress-equivalent plastic strain curves from experiment [3] and simulation (fitted). The experimental curve is used as reference for fitting the model parameters. A nonporous RVE with 520 grains is shown as inset, where BD, TD, and ND stand for building, transverse, and normal directions.

Influence of Pore Volume Fraction
The Abaqus-Python script removes mesh elements from the RVE such that three variants of the RVE with increasing pore volume fractions of 2%, 5%, and 8% are generated. The size of the pores in all three RVEs correspond to 6.2 μm (refer Table 2). Uniaxial tensile test simulations with loading parallel to the BD are performed by using Abaqus, [27] and the results are plotted in Figure 9.
The simulation flow curves of the porous RVEs are compared with the nonporous one, and it is evident that, with increasing pore volume fraction, the strength of the material reduces. To better quantify this loss of strength, the modified Voce law, [48] which is an empirical isotropic hardening law, is fitted to the flow curves. The modified Voce law is given by where Y 0 , R 0 , R inf , and b are model parameters, and ε p is the equivalent plastic strain. The Voce model parameters are fitted using the nonlinear least squares fitting method, [49] and the flow curves thus generated are shown in dashed line style in Figure 9.
The four sets of Voce model parameters corresponding to the nonporous and three porous test cases are plotted in Figure 10. The secondary y-axis is used to show the percentage change in their values with respect to the nonporous test case. Y 0 corresponds to the initial yield stress, and its value decreases linearly, indicating loss in strength with increasing pore volume fraction. Quantitatively, the percentage loss in yield strength for the extreme case of 8% pore volume fraction with respect to the nonporous RVE is %9%. R 0 and R inf correspond to the slope and the difference between the saturation and the initial yield stress, respectively. Their values also decrease linearly similar to that of Y 0 with the maximum variation within 10%. The hardening parameter b governs the rate of saturation of the exponential term in Equation (4). The variation in values with respect to the nonporous case is the smallest for this parameter, indicating that the saturation rate remains approximately the same although the porosity increases in volume fraction. With a linear variational trend observed for Y 0 , R 0 , R inf values and an almost negligible variation in the b value, it can be concluded that the porosity has marginal influence on the hardening behavior.
In addition to the flow curves, the distribution of stresses in the deformed RVEs is further investigated, as shown in Figure 11. Compared to the nonporous RVE, the RVEs with pores show a localized concentration of stresses around the pores. Such localized stress concentrations are vital for the onset of local events such as damage. With the quantification of the Voce parameters, a pronounced effect on the hardening behavior is not observed as the pore volume fraction increases. However, the presence of pores strongly influences the localized stress concentration, possibly leading to other important local events.
Value 247 [59] 106 [59] 71 [59] 147 250 20 2.25 160 0.03 1.4 Figure 9. Equivalent stress-equivalent plastic strain curves for RVEs with increasing pore volume fraction. Uniaxial tensile test simulation results are represented by markers and the corresponding the modified Voce law fitted values are represented by curves.
www.advancedsciencenews.com www.aem-journal.com Figure 10. Voce model parameters fitted to the flow curves in Figure 9 are plotted against the pore volume fraction. Percentage changes in the parameters with respect to their corresponding maximum value are plotted as secondary y-axis for each subplot. To study the influence of local events on the mechanical behavior, the damage model described in Section 3.2 is used in the pore shape study discussed next.

Influence of Pore Shape
A nonporous RVE with 492 grains (91 125 mesh elements) along with a customized Abaqus-Python script is used to generate RVEs with equiaxed pores of 8% porosity volume fraction and elongated pores of 2% and 8% porosity volume fraction (see Table 2). The long axes of the elongated pores are oriented along the BD of the sample. Uniaxial tensile test simulations with loads parallel to the BD and TD are performed, and the results are plotted in Figure 12. The damage model outlined in Section 3.2 is used to describe the local damage behavior. In this regard, the lower limit (p l ¼ 0.15) and the upper limit (p u ¼ 0.22) for the equivalent plastic strain are chosen such that the resulting model, with elongated pores (8%) aligned along the BD, reaches its uniaxial tensile strength at around 8% total strain. [50] Figure 12a shows that the pore shape in combination with the loading direction has significant influence on the flow behavior. It can be seen that the flow curves for the RVE with elongated pores show anisotropy between the two loading directions, and the anisotropy gets amplified with increasing pore volume fraction. The difference in the yield points and the corresponding hardening behavior between the curves can be attributed to the cross-sectional shape of the pores orthogonal to the particular loading direction. The softening observed is a consequence of the material stiffness degradation included in the damage model. The shape of the pores and the crystallographic orientation of the grains surrounding the pore directly contribute to the loss of material strength.
The strain at which plastic instability occurs can be computed using the Considère criterion. [51] Based on this definition, a tangent line of the true stress (σ t )nominal strain (ϵ n ) curve defines the point of instability. The nominal strain value thus computed is converted into a true plastic strain (ϵ pl t ), and the points of instability (σ t , ϵ pl t ) for each of the flow curves in Figure 12a are represented by markers and listed in Table 5. Comparing the strain values between the two loading directions, the plastic instability can be seen to occur earlier for the TD loading for RVE with elongated pores, whereas for the equiaxed pored RVE the plastic instability along BD and TD occur at approximately the same strain value.
The evolution of pore volume fraction is plotted in Figure 12b as a function of equivalent plastic strain. Elongated pores loaded along the TD tend to grow significantly faster than the ones loaded along the BD, thus contributing to the pronounced damage observed in Figure 12a. In contrast, growth of the equiaxed pores for the two loading cases is similar and lies in between that of the 8% elongated loading cases. The crosssectional views of RVEs with elongated pores (8%) loaded along BD and TD at various values of homogenized equivalent plastic strains are shown in Figure 13a,b. The severely deformed elements are removed from the FE mesh by Abaqus once they satisfy the damage criterion (i.e., D ≥ 0.99). Thus, with this element deletion technique, the growth and coalescence of pores are modeled effectively in a mean field approach. It is evident from these figures that pore coalescence in the TD loading occurs earlier at smaller equivalent plastic strains and, furthermore, signifies the effect of the loading direction on the onset of plastic instability with respect to the pore shape. Irrespective of the loading direction, the material softens as pores start to coalesce (compare Figure 12a with Figure 13a,b), and with further applied load, the coalescence intensifies. The intensities of the plastic strain contour around the pores signify that the pore shape in combination with the crystallographic grain orientations around the pores contributes to the local equivalent plastic strains.  In contrast, the RVEs with equiaxed pores show no such anisotropic behavior until the onset of damage. This clearly indicates the influence the pore shape has on the mechanical behavior. Figure 14a,b shows the cross-sectional views of RVEs with equiaxed pores (8%) loaded along the BD and TD at various values of homogenized equivalent plastic strains. Similar to Figure 13a,b, the pores coalesce along the directions perpendicular to the loading. But, in contrast to the elongated pores, the onset of coalescence for the two loading directions occurs at similar homogenized equivalent plastic strains. Although the two loading directions show similar yield strength and hardening behavior, the material softening due to damage shows an anisotropic trend. As the pores are equiaxed, there exists no shape effect that contributes to this behavior. It is rather the crystallographic orientations of the grains surrounding the pores that are responsible for this anisotropy. These grains deform differently for the two loading directions because their orientations influence the local stress state.
To qualitatively determine the influence of the crystallographic orientations on the softening anisotropy, two additional  www.advancedsciencenews.com www.aem-journal.com configurations of RVEs with equiaxed pores are generated (equiaxed pore configuration 2 and 3). As the pore placement is random, the flow curves for the two new configurations when loaded along the BD and TD are expected to be different. Figure 15 shows the softening anisotropy among various configurations of RVE with equiaxed pores. The influence of the crystallographic orientations on the softening anisotropy is evident, as the anisotropic difference varies with each configuration. The simulation results presented in this work are in general in line with the results of Carlton et al. [52] and Choo et al., [53] who have performed experimental investigations on an L-PBFfabricated 316L alloy to understand and quantify the effect of various pore characteristics on the mechanical behavior. In situ tensile tests using synchrotron radiation microtomography were performed by Carlton et al. [52] to track the damage evolution, and their results show that pore size, morphology, and distribution significantly affect the fracture mechanism. They observed that pore growth and coalescence, which eventually leads to material failure, was defect driven with cracks originating from pre-existing pores. And, specimens with low porosity displayed high amounts of ductility in contrast to high-porosity specimens. These results correlate well with our simulation results (see Figure 12a).
In the investigation of Choo et al., [53] in situ tensile tests using high-energy synchrotron X-ray diffraction were performed to study the effects of pore type, density, and orientation with respect to the loading direction on the local stress state of the material. According to their results, specimens with elongated pores (%3%) aligned perpendicular to the loading direction exhibited lower yield strength and limited ductility in comparison to specimens with pores aligned parallel to the loading direction. This behavior has been attributed to the difference in the local stress concentrations that arise near the curvatures of the pores when loaded along different directions. Similar observations have been made from our simulation results (see Figure 12a, Figure 13a,b, and Table 5), which demonstrates the capability of the model in describing real material behavior. Based on these observations, the simulation framework can be extended to include other microstructural features (such as melt pool geometry, crystallographic texture, pore size distribution, and cracks) that play a major role in affecting the mechanical properties of AM metals.

Conclusions
In this work, a micromechanical model has been used to study the anisotropic mechanical behavior of L-PBF-manufactured metal in the presence of pores. The complex porous microstructures used in this study have been synthesized using statistical data obtained from EBSD and LOM micrographs. The nonlocal CP model, [26] originally developed to describe the hardening caused by strain gradients in polycrystal deformation, has been extended to include a damage evolution law to simulate material softening. Furthermore, by coupling the damage model with the element deletion technique, the growth and coalescence of pores have been modeled effectively in a mean field approach. In this way, the evolution of a pore volume fraction is obtained that can be used to investigate the influence of pore shape on the mechanical behavior. It has been demonstrated that the anisotropy in the hardening and damage behaviour between the different loading directions is influenced by the shape of the pores. Microstructures with elongated pores show significant anisotropy with increasing pore volume fractions and earlier onset of damage when loaded orthogonally to the orientation of the long axes of the elongated pores. In contrast, microstructures with equiaxed pores do not exhibit such anisotropy in mechanical behavior. However, with ongoing evolution of the damage a preferred direction evolves, and the material becomes anisotropic in nature. This is caused by the crystallographic orientations of the grains surrounding the pores, which results in an anisotropic plastic deformation. Although the CP models are computationally demanding, to account for the anisotropic material behavior due to complex local events they are the ideal choice. Modeling pores that follow a size distribution and investigating the effects of pore sizes on the mechanical behavior can be incorporated into the model in a straightforward manner. From a future perspective, an in-depth investigation of the effect of texture on porosity evolution is expected to increase our understanding of local damage evolution.

Nonlocal CP Model
The nonlocal CP model proposed by Ma and Hartmaier [26] has been used in this work, the fundamental aspects of which are detailed in the works of Rice, Peirce et al., and Kalidindi. [54][55][56] This section describes only the nonlocal aspects of the model, which has been validated experimentally in Engels. [57] Starting from the kinematics of deformation, the multiplicative decomposition of the total deformation gradient F into its elastic F e and plastic parts F p is given by Figure 15. Comparison of equivalent stress-equivalent plastic strain curves for RVEs with different configurations of equiaxed pores. The placement of pores in the three configurations is randomized; however, the porosity volume fraction remains 8%.