Elasto-plastic residual stress analysis of selective laser sintered porous materials based on 3D-multilayer thermo-structural phase-field simulations

,


Introduction
In the recent decade, additive manufacturing (AM) emerged from a niche technology to a widely applicable means of production.Among a variety of AM techniques, powder bed fusion (PBF) stands out as a strong technique for producing intricate geometries while maintaining good structural integrity and superior material properties [1,2,3,4,5].In PBF, layers of powder get fused by a beam one after another to build up three-dimensional (3D) geometries.Compared to conventional ones, this manufacturing process offers several advantages, such as great design flexibility, reduction in waste, enabling rapid tooling and decreasing production cycles.
Selective laser sintering (SLS) is one prominent example of PBF, which utilizes the tuned incident beam (mostly continuous laser scan or laser pulse) to bind the free powders layer by layer.Due to its well-controlled powder bed temperature compared to other PBF techniques like selective laser melting (SLM), where the significant melting phenomenon can be captured, SLS enables only partial melting of particles and produces samples with relatively high porosity [6,7,8,9,10].In this regard, SLS has great potential in applications that require complex geometries with high porosity.For instance, SLS is applicable in manufacturing porous components for medication applications, especially medical scaffolds, and artificial bones [11,12,13].It is feasible in manufacturing functional materials with a low processing temperature, such as ferromagnetic materials [14,15,16,17,18].By precisely controlling the geometry and structure of the printed material, it may also be possible to create the desired strain gradients and electric polarization necessary to generate the flexoelectric effect [19,20,21].
Residual stress has been a critical issue since it can affect mechanical properties, dimensional accuracy, corrosion resistance, crack growth resistance, and performance of AM parts [22].There are already investigations focusing on the spatial distribution of the residual stress within the AM-built parts, which were carried out by experimental measurements and theoretical estimations [23,24,25,26].Mercelis et al. explained the origins of residual stress in parts produced by AM.They experimentally measured the residual stress distribution in a selective laser sintering produced part and compared it with analytical and numerical solutions [23].Pant et al. employed a layer-by-layer finite element approach to predict the residual stress and validated the model using measurements from neutron diffraction [25].Ibraheem et al. predicted the thermal profile and residual stress of SLS processed H13 tool steel using a finite element model.High porosity was observed in the range 25-40%, yet a homogenized powder bed was utilized to simulate the thermal evolution and subsequently the residual stress [27].
The calculation of residual stress in literature is often performed in two sequential steps: first the transient temperature is simulated in the entire domain and then the thermal history is used to calculate the stress and strain evolution in thermomechanical simulations.The accuracy of these calculations critically relies on the quality of the thermal history.The temperature gradient mechanism (TGM) model and the cool-down stage (CDS) model are two commonly used models to explain the development of plastic strain and residual stress formation mechanisms [23,28,29].In the TGM model, the deformation of the material in the molten pool/fusion zone is restrained by surrounding materials during both the heating and cooling stages due to a large temperature gradient inside and outside the overheated region (where the on-site temperature is beyond the melting point).In the heating stage, the plastic strain is developed in the surrounding material due to an expansion of the heated material.As the heated material cools down, it starts to shrink, which is, however, counteracted by the plastically deformed surrounding material.Thus, a tensile residual stress develops.The CDS model, on the other hand, was proposed to elucidate the residual stress due to the layerwise features of the AM process and the AM-manufactured parts.This occurs because, in both the heating and cooling stages, the deformation of the upper-layers is restrained by the fused lower-layers and/or substrate.Meanwhile, melted material in the lower-layers will undergo a remelting and re-solidification cycle.Both factors result in tensile residual stress in the top layer due to the shrinkage of the material during the cooling stage.These two models provide phenomenological aspects by employing idealized homogeneous layers for the understanding of the formation mechanism of residual stresses.They have been widely adopted in numerical analysis of residual stress at single layer scale models [30,31,32,33] and part-scale multilayer models [34,35,36].
In these aforementioned works, homogeneous layers are used for analyses of the thermal history and the mechanical response in an AM-manufactured part.However, due to the complex morphology of the powder bed on the powder scale and the resultant inhomogeneity of the local thermal properties, the thermal profile is subjected to a high level of non-uniformity as well.This implies heterogeneity of the thermal stress, the plastic deformation and the residual stress on the powder scale.In fact, mechanical properties of AM-built parts, such as fracture strength and hardening behavior, are significantly influenced by the local defects and local stress [37,38].The diversity in the powder bed structure and the packing density introduces thermal heterogeneity in the form of mesoscopic high-gradient temperature profile and asynchronous on-site thermal history.Such variability leads to varying degrees of thermal expansion and, therefore, thermal stresses.Additionally, the presence of stochastic inter-particle voids and lack-of-fusion pores contribute to evolving disparities in thermo-mechanical properties [7,39,40].The thermal gradient across the powder particles is very sharp as per the modeling results by Panwisawas et al. [41], which is in agreement with experimental fundings [42] and is also recaptured in our former numerical works [7,39,40].Furthermore, the layerwisely build-up process results in local interaction between the newly deposited and previously deposited layers with high surface roughness.Thus, the residual stress formation on the powder scale should be revisited for the SLS process.
For this purpose, the powder bed morphology, the heat transfer and the porous microstructure evolution during the SLS process should be primarily considered.One promising approach is the phase-field simulation.For instance, phase-field simulations of the AM process can reveal the in-process microstructure evolution, providing insights into key features such as porosity, surface morphology, temperature profile, and geometry evolution [43,7].Based on a non-isothermal phase-field multilayer simulation results of selective laser sintered porous structures, thermo-elastic simulations and homogenization of the elastic properties have been investigated in our previous work [39].In combination with the single layer phase-field thermo-microstructural simulations, the first powder-resolved elasto-plastic simulations of the residual stress in a SLS-processed magnetic alloy has also been performed [40].Results show that during the cooling stage, the partially melted and interconnected particles result in plastic deformation due to the shrinkage of the fusion zone, but mostly at the necking region because of stress concentration in porous microstructures.By combining the thermo-structural phase-field method and the thermo-mechanical calculations, this work demonstrates the promising capability of the multiphysics simulation scheme for the local residual stress analysis of AM-built materials on the powder scale.In the current article, we extend this multiphysics simulation scheme for elasto-plastic multilayer SLS process, providing systematic simulations of the local plastic deformation and the residual stress.These allow us to reveal the influential aspects related to the printing process and the powder bed.

Results
In this work, the simulation workflow is arranged in three stages, as shown in Fig. 1a.The powders were first deposited into a prime simulation domain based on the discrete element method (DEM) under a given gravitational force.Then, non-isothermal phase-field simulations were performed for the coupled thermo-structural evolution during the SLS process.Upon completion of a layer, the resulting microstructure is voxelized and re-imported into the DEM program for the deposition of the next powder layer until reaching the desired layer number or height.Meanwhile, the subsequent thermo-elasto-plastic calculations were performed to investigate the development of plastic strain and stress under the quasi-static thermo-structures by mapping the nodal values of the order parameters (OP, indicating the chronological-spatial distribution of the phases) and temperature to a subdomain.It is worth noting that the thermo-elasto-plastic calculation of the next layer was restarted from the former one with continuous history of the residual stress and accumulated plastic strain, though the OPs and temperature were reiterated in the prime domain.As it is usually done in the literature, we assume thereby that heat transfer is only strongly coupled with microstructure evolution (driven by diffusion and underlying grain growth), but mechanics impose negligible effects on both heat transder and microstructure evolution.
SLS processes with a constant beam spot diameter D L , and various combinations of beam power P as well as scan speed v were simulated.The processing window P ∈ [15,30] W and v ∈ [75, 150] mm s −1 was set in order to compare it with previous studies [7,39].D L = D FWE2 = 200 µm (FWE2 represents the full-width at 1/e 2 ) was adopted as the nominal diameter of the beam spot, within which around 86.5% of the power is concentrated.The full-width at half maximum intensity (FWHM) then equals to 0.588D FWE2 = 117.6 µm, characterizing 50% power concentration within the spot.The powder bed is preheated with a preheating temperature T 0 = 0.4T M = 680 K, which was embodied by temperature initial condition (IC) and boundary condition (BC).Both powder and substrate are SS316L.Powders are firstly deposited into a 250 × 500 × 400 µm 3 prime simulation domain under a given gravitational force based on the discrete element method (DEM).This domain contains the free space of 150 µm height for positioning powders and a substrate of 250 µm thickness.In-total four layers of SLS processes were performed for each pair of processing parameters (i.e., P and v, hereinafter as P-v pair).Each newly deposited layer maintains the same powder size distribution and a volume fraction of 12% with respect to the free space of 250 × 500 × 150 µm 3 .This guarantees each deposited layer sufficiently holds one to two particles in thickness, as illustrated in Fig. 1b.The powder size follows Gaussian distribution with a mean diameter of µ d = 20 µm, a standard deviation ς d = 5 µm, and a cut-off bandwidth d cut = 10 µm.Fig. 1c presents the powder size distribution of the first layer (common for all P-v combinations).After each layer's beam scan (i.e., the beam leaves the prime domain), the simulation continues for twice the duration of the scanning period to emulate a natural cooling stage.Details regarding the process as well as thermo-elasto-plastic modeling, simulation setup, and parameterization are explicitly given in the Methods section.

Thermo-structural evolution during multilayer SLS
Fig. 2a 1 -a 4 present the transient thermo-structural profiles for a single scan of P = 20 W and v = 100 mm s −1 where the beam spot is consistently positioned along scan direction (SD) across various building layers.In the overheated region (T ≥ T M ), particles can undergo complete or partial melting.This prompts molten materials to flow from convex to concave points, making the fusion of particles possible and, therefore, forging the fusion zone.In regions where the temperature remains below the melting point (T < T M ), however, the temperature is still high enough to cause diffusion as the diffusivity grows exponentially w.r.t.temperature.This is evidenced by the necking phenomenon among neighboring particles.Temperature profiles demonstrate a strong dependence on the local morphology.Notably, the concentrated temperature isolines can be observed around the neck region among particles, indicating an increased temperature gradient (up to around 50-100 K µm −1 , comparing to 1 K µm −1 in the densified region).It is worth noting that such thermal inhomogeneity induced by stochastic transient morphology can hardly be resolved by the simulation works employing homogenized powder bed [29,44,26].
The heat dissipation in the powder bed also changes considerably as upper-layers are continuously built and processed.This can be noticed by a varying shape and significance of the overheated region across layers.As a comparison, the overheated region at 4 th layer (Fig. 2a 1 ) is greater in comparison to the initial layer (Fig. 2a 4 ).This can simply be reasoned by increasing porosity in fused lower-layers that block the heat dissipation.At the initial layer, this dissipation is affected the least as there is no lower-layer but substrate.We further probed the temperature history at surface points located at the center of the scan path, as shown in Fig. 2c 1 , where the recurring peaks in every single-layer SLS stage (shaded section) present the thermal cycle.The first peaks of distinctive cycles (probed on different surface points) nearly match.In an identical cycle (e.g., C1), the peaks gradually decrease as the SLS processes advance.Once the spot leaves the probing point, the temperature drop is rapid at first (mostly due to the heat conduction driven by the high-gradient temperature profile) and then becomes moderate (due to heat convection and radiation) as the cooling stage begins, when the heat convection and radiation are effective.This steep temperature drop gradually disappears at points from fused lower-layers as the upper-layers are continuously built, comparing the thermal cycles at point C1.Meanwhile, the temperature drop at distinct points gradually coincides as the cooling stage proceeds, comparing the thermal cycles at points C1-C4.
The multilayer thermo-structural evolution is also significantly affected by processing parameters, as the profiles of varying beam powers and scan speeds presented in Fig. 2b 1 -b 4 .Increasing the beam power and/or decreasing the scanning speed is observed to intensify heat accumulation at the beam spot, resulting in an enlargement of the overheated region.When v is held constant at 100 mm s −1 , decreasing P from 30 W (Fig. 2b 1 ) to 15 W (Fig. 2b 2 ) results in a less pronounced overheated region.On the other hand, maintaining a constant P = 20 W and increasing v from 75 mm s −1 (Fig. 2b 3 ) to 150 mm s −1 (Fig. 2b 4 ) leads to a reduction in the overheated region as well.It is also evident that increasing P and/or decreasing v produces less porosity in fused lower-layers, which may further change the thermal conditions and improve the heat dissipation.Fig. 2c 2 presents the probed thermal cycle of point C1 under various P with v = 100 mm s −1 maintained.One can tell that the temperature of the case P = 30 W quickly drops from a higher peak to the value coincided with one of the case P = 20 W at the end of scan duration, implying an improved heat dissipation.It should also be noted that surface point C1, located at the initial layer, experienced three peaks of temperature that are beyond T M due to the extended depth of the fusion zone at higher beam power.In Sec.3.1 we will continue the discussion regarding the relation between porosity and processing parameters.

Stress and plastic strain evolution during multilayer SLS process
To analyze the overall stress and plastic strain evolution during SLS, the average quantities (incl.von Mises stress σ e , effective plastic strain p e , and temperature T) in the powder bed are defined as σP e = Ω ′ ρσ e dΩ Ω ′ ρdΩ , pP e = P Ω ′ ρp e dΩ Ω ′ ρdΩ where Ω ′ is the volume of the simulation domain without the substrate meshes, and ρ is the OP indicating the substance with ρ = 1.They are hereinafter termed as PB-averaged quantities.Fig. 3a presents that σP e develops as the SLS process advances on distinctive layers.When the beam spot enters the domain, accompanied by the appearance of the overheated region, σP e begins to decrease while T P continues to rise to its peak.This is attributed to losing structural integrity caused by full/partial melting within the overheated region, resulting in zero stress as shown in Fig. 3b 1 .The areas surrounding the overheated region also exhibit relatively low stress due to a significant reduction in stiffness at elevated temperatures.Stress accumulates faster around concave features, such as surface depressions and particle sintering necks, compared with the stress around convex features as well as unfused powders away from the fusion zone.This is due to localized temperature gradients, as explained in Sec.2.1.As cooling progresses, stress decreases in convex features and unfused powders, yet remains concentrated around concave features (see Figs. 3b 2 -3b 3 ).The stress at the end of each cooling stage (Figs.3b 3 -3b 6 ) serves as the residual stress of corresponding processed layers.After the deposition of new layers, σP e has an instant drop due to the addition of stress-free substances (powders), then a new cycle begins.It is worth noting that the peak values of the stress cycle (which are also residual stresses of each processed layer) are almost identical, implying that each single-layer SLS generates nearly the same amount of residual stress by average in the powder bed.For upper-layers, σ e experiences an additional reduction in the cooling stage, which may be attributed to the more delicate variations in the on-site stress components and will be discussed in Fig. 4.
The accumulated plastic strain p e , on the other hand, presents an overall increasing tendency vs time during the single-layer SLS and cooling stage, as shown in Fig. 3c.This distinguishes the plastic strain history from the stress cycle, which suffers a reduction during SLS due to loss of structural integrity.The continuous accumulation of p e leads to a locally concentrated profile in the vicinity of the melt zone (see Fig. 3d 2 -3d 3 ).The high temperature gradient at the front and bottom of the overheated region during single layer SLS also results in the localized p e surrounding the fusion zone (see Fig. 3d 3 -3d 6 ), where the existing high thermal stress locally activates the plastic deformation.A similar explanation applies to the concentrated p e around pores and concave features such as sintering necks near the fusion zone.In contrast, unfused powders and the substrate away from the melt zone show minimal p e , since the thermal stress at these locations is insufficient to induce material plastification due to lower local temperatures.
We probed the history of the stress components at six points, where L 1 0-L 4 0 locates at the center of 1 st -4 th layers, L 1 1 locates at the boundary of the fusion zone, and L 1 2 locates outside the fusion zone in the 1 st layer, as shown in the inset of Fig. 4. The results mainly show the difference of the stress field inside and outside the fusion zone.Since the temperature cools down from above the melting temperature to pre-heating temperature, the rise of the elasticity leads to a significant increase of σ e on point L 1 0 for temperature below 0.65T M .As a comparison, the increament of the σ e becomes smaller for the point L 1 1 and the point L 1 2, as shown in Figs.4a 1 -4a 3 .It can also be seen that σ e for the mentioned three points (L 1 0-L 1 2) during SLS of 2 nd -and 3 rd layers are very similar when they are all outside the fusion zone.Similar results can be observed for the development of normal stress components σ ii (i = x, y, z).The main differences in normal stress for the three points are mostly for the SLS process of the 1 st layer.σ ii are tensile in the fusion zone while compressive outside the fusion zone.This is due to the thermal stress in the fusion zone being tensile while it is compressive outside the fusion zone.Curves for the 2 nd and 3 rd layers are similar because the three points are all outside the fusion zone.Shear stress components σ ij (i, j = x, y, z; i ̸ = j) remain small inside and outside the fusion zone.The variation of the shear stress is comparably more visible at the boundary of the fusion zone (L 1 1), as shown in Fig. 4a 2 .This is due to a large gradient of the thermal strain and material mechanical properties across the boundary of the fusion zone.With a larger variation of the shear stress components, variation of σ e at the highest temperature during SLS process of each layer is also more visible.The saturation of the stress field can be observed after the 2 nd layer.It means that the upper layer has limited influence on the plastic deformation of previous layers, which experience lower temperatures and lower temperature gradients outside the fusion zone.
For points in the layer-wise direction (L 1 0 to L 4 0 at the boundary of each layer), the results in Fig. 4b show a direct comparison between the development of the stress field of a point in the microstructure during SLS process of different layers.It is worth noting that the development of the σ e of each point is very similar at the beginning.The saturation of the σ e for each point can be observed in the layer-wise direction when the points are outside the fusion zone.However, variation of the σ e is much more pronounced when the overheated region is moving close to the point.This indicates that the point is located around the boundary of the fusion zone, which encounters a large variation of mechanical properties of the material w.r.t. the high-gradient temperature and thus undergoes a large variation of the shear stress.

Powder-resolved mesoscopic residual stress formation mechanism
For SLS-processed porous microstructures, the residual stress is generally related to stress concentration under thermal loading in the microstructure.In Fig. 5a 3 , and 5a 5 , the particles undergo a low degree of fusion (due to low level of overheating) due to relatively low beam power or high scan speed.The stress concentration mostly occurs at the necking region, leading to accumulated plastic strain and residual stress.In Fig. 5a 1 and 5a 4 , with a higher degree of fusion due to a relatively high beam power or low scan speed, accumulated plastic strain and residual stress tend to concentrate on both necking region and inter-layer region (as indicated by the black lines in Fig. 5c and 5d).Note that the inter-layer boundary has a high surface roughness which also causes stress concentrations.In Fig. 5a 2 , with the combination of highest P and lowest v within the processing window, residual stress can be observed in the whole fusion zone, which makes it difficult to identify the residual stress concentration.Besides, accumulated plastic strain and residual stress can be observed on the boundary of the substrate and the 1 st layer, especially at high beam power.This is because a high beam power generates a much larger fusion zone, which even penetrates into the substrate.When the substrate with continuum material is melted, it generates a large plastic deformation and thus the residual stress, which can be significantly larger than those in the porous microstructure, as shown in Fig. 5a 2 .These results show that the residual stress is directly related to the degree of fusion and thus related to SLS processing parameters.
To gain a better understanding of the distribution of residual stress, we analyzed the stress state of a few representative positions in the microstructure, as shown in Fig. 5c and 5d.The former surfaces of fused layers are marked by black lines, and the bottom of the fused strut, a.k.a., the fusion zone boundary (FZB) of the 1 st layer, is indicated by white lines.Fig. 5c 1 and 5d 1 show that the residual stress and the accumulated plastic strain are solely concentrated in the necking region and inter-layer region for the case with a low degree of fusion.In Fig. 5c 2 , with the high degree of fusion (the smallest porosity acheived in the selected processing window), the residual stress and the accumulated plastic strain still tend to be solely concentrated in the inter-layer region, yet less distinguishable compared to the former case.By comparing Lamé's stress ellipsoids at six points on the boundaries of different layers, it is worth noting that the stress states close to the boundary between layers are very similar, e.g. the stress states of P 2 and P 3 are similar to that of P 5 and P 6 , no matter how high the degree of fusion is.It means that the residual stress is still formulated due to stress concentration at the inter-layer region.On the other hand, the stress states of P 1 and P 4 at the top layer are quite different.This is because the residual stress of the top layer is directly determined by the thermal loading which depends on the surface morphology and morphology-induced thermal inhomogeneity.
Based on the above observations, we propose a powder-resolved mesoscopic residual stress formation mechanism, which is summarized in the schematic illustrated in Fig. 6.The residual stress formulated in porous microstructures manufactured by a multilayer SLS process contains two primary contributions: (i) Residual stress is directly caused by the stress concentration at the necking region of partially melted particles under thermal loading.The partially melted particles are inter-connected via necking regions, which is the weak link of the microstructure.Therefore, both the thermal expansion in the overheated region and the thermal contraction during the cooling stage cause severe plastic deformation at the necking region, which is one primary source of the accumulated plastic strain and residual stress.
(ii) Residual stress due to interaction between the upper-and lower-/fused layers in the layerwisely build-up process.In the cooling stage, the shrinkage of the upper-layer after overheating results in tensile stress on itself and compressive stress on the lower-/fused layer, which causes plastic deformation of the porous structure, especially at the inter-layer region with stress concentration due to a high surface roughness, which is the other primary source of the accumulated plastic strain and residual stress, as indicated by the white lines in Fig. 6a.
The proposed mechanism is based on the detailed simulation results using the powder-resolved thermo-mechanical model, which evidently demonstrates the formation and distribution of residual stress in the porous structure.The accumulated plastic strain results in structural distortion of the fused strut, which is schematically illustrated in Fig. 6a and demonstrated by the contour plot of deformation component u z in Fig. 6b.Compared to the aforementioned TGM and CDS models proposed by Mercelis and Kruth [23] correspondingly, which ignored the difference between parts manufactured by SLM and SLS, the partial melting as the main fusion mechanism plays an important role in the residual stress in the microstructure.There are some similarities, the concept of the TGM model also applies to the proposed mechanism because the temperature gradient indeed directly leads to residual stress.However, in the SLS process, since the power bed is fully relaxed before fusion, particles provide very weak constrain on the fusion zone and thus have limited plastic deformation.The plastic deformation is majorly accumulated during the cooling stage due to thermal contraction and stress concentration at the necking region.For the CDS model, the proposed model also considered the shrinkage of the top layer resulting in a tensile residual stress in the top layer and a compressive residual stress in the previous layers.However, the top layer in the simulated SLS process has very complex boundaries with previous layers.These boundaries are the source of the stress concentration which leads to the dominant residual stress at the inter-layer regions.
Based on the TGM model and the CDS model, Mercelis and Kruth [23] also introduced a theoretical model to predict the relationship between residual stress and the part height.Similarly, we can use the proposed powder-resolved model and mesoscopic residual stress formation mechanism to predict the dependence of the residual stress on the microstructure porosity, which is also a characteristic geometric parameter of a porous microstructure.Since the porous structure induces complex inhomogeneity of both material properties and thermal and mechanical loading, and we also need to consider the porosity is directly related to SLS processing parameters, it is more straightforward to propose phenomenological models based on our simulation results to predict the relationship between the residual stress and SLS processing parameters, as will be discussed in the next section.

Phenomenological relation for porosity control
Controlling the porosity of processed sample through the tuning of the processing parameters plays a central role in tailoring the end-up performance of the porous materials in AM, as many properties are determined by (or related to) porosity, including but not limited to mechanical strength, permeability, acoustic/optical absorptivity and various effective conductivities.We start the discussion with phenomenologically relating the porosity to the processing parameters (in this work P and v).In this work, a nominal porosity is defined as with the substance OP ρ.Ω ′′ represents the volume of a post-processed simulation domain with surface and unfused powders sufficiently removed (termed as "virtual polishing").Fig. 7a presents the local porosity that evaluated segment-wise along z-(BD) and y-direction (perpendicular to SD) to help identify the representative domain with a width W R and a height H R for porosity calculation, while the complete length (L R = 250 µm) along x-direction (SD) is included due to exhibited relatively minor variations in local porosity [39].We then conduct the virtual polishing on simulated multilayer SLS-processed microstructures with H R = 60 µm (containing the coordinate range z ∈ [−135, −75] µm) and W R = D FWHM , and proceed calculating their porosity using Eq.(2) with the post-processed domain size The resulting P-v map of porosity is in shown Fig. 7c.
Selected microstructures illustrated in Figs.7b 2 -7b 4 for the cases varying P with constant v; Figs.7b 5 -7b 7 for the cases varying v with constant P. Fig. 7b 1 is the microstructure under a reference processing parameter (P = 20 W, v = 100 mm s −1 ), while Figs.7b 8 and 7b 9 are respectively the one with minimum and maximum porosity.It demonstrates that improved densification can be achieved by either increasing P or decreasing v, where smoothed surface morphology implies an enhanced partial melting.Porosity drops from 31% to 18% for the increase of P from 15 to 30 W (v = 100 mm s −1 ) and from 29% to 21% for the decrease of v from 150 to 75 mm s −1 (P = 20 W).The presented tendencies imply a possible allometric relation φ ∝ P −m v n with positive indices m and n.On the other hand, combining beam power and scan speed as one characteristic quantity, P/v is widely adapted to define a specific energy density, notably the volumetric energy input as where H is the powder layer thickness (here H = 60 µm) and W is the scan track width (here W takes D FWHM ).Evidently, φ exhibits an overall decrease tendency with rising U, meaning lower porosity can be achieved as higher specific energy input.
We then examined the simulation results with a proposed phenomenological relation following Refs.[45,7], read as with a densification factor defined as ϱ indicates the ratio between a reduced porosity w.r.t. the maximum porosity reduction in the chosen processing window, as φ 0 and φ min represent the initial and minimum achieved porosity, respectively.φ min usually varies between 3% to 30% for metallic materials [45].In pursuit of uniformity, we have chosen φ min = 3% in this study, aligning with our previous single-layer SLS simulation [7].φ 0 , however, is difficult to be determined for multilayer SLS simulation since there is sorely one or two particles in thickness for every new layer while the old layers have been fused.In this sense, we evaluate φ 0 by three routes: (i) Assuming that the porosity of the powder bed region, which is away from the beam spot (i.e., without significant densification), is the initial one of the powder bed, φ 0 can be read as the converged value from segment-wise porosity evaluation along y-direction (Fig. 7a 2 ), which is φ 0 = 27.8%∼ 33.6%.
(ii) As the particle size distribution in this work is assumed to be Gaussian-type, φ 0 is then estimated by statistic randomclose-packing model of spherical particle as φ 0 = 0.366 − 0.0257 (ς d /µ d ) with µ d and ς d the mean and standard deviation of the particle diameter [46,44].Taking µ d = 20 µm and ς d = 5 µm, it can be calculated as φ 0 = 36.0%.
(iii) We also piled multiple densification-free powder stacks using DEM method with the same domain volume fraction of the particles (48.0%) as the overall one of the multilayer SLS simulation, in which each layer deposits powders of 12% domain volume fraction.The measured nominal porosity is thereby regarded as the initial one, as φ 0 = 40.0%∼ 40.5%.φ 0 from route (i) directly reflects the porosity of on-site unfused powders, but the influences from the potential necked particles due to thermal processes cannot be well eliminated.φ 0 from route (ii) is based on the statistics of random-closed-packed particles, which is rather ideal compared to the practical powder bed created by powder spreading.Since route (iii) creates particle stacks via simulating the deposition in powder spreading, the calculated φ 0 may be still inflated as the morphological variability of the deposition surface is missing.Considering all of these factors, we have selected φ 0 = 36.0%,which stands at the midpoint among all the assessed values.In Fig. 7d, linear regression between − ln(1 − ϱ), calculated from simulated φ, and U is presented.For comparison, regressions on data by single-layer SLS simulation and experiment are also illustrated [45,7].The regressed line together with the 95% confidence interval (CI 95% ) of the multilayer SLS simulations is right in between those of the single-layer SLS simulation and experiments.As the coefficient K ϱ is related to the material and the size distribution of the powders, the multilayer result K ϱ = 0.016 ± 0.001 mm 3 J −1 demonstrates an improved coherence with the experimental one K ϱ = 0.013 mm 3 J −1 compared the single-layer one K ϱ = 0.019 mm 3 J −1 , which can suffer from inflated porosity mostly due to insufficient volume microstructure for porosity calculation.It is also worth noting that this relation between − ln(1 − ϱ) and U is examined on the microstructures both with and without fusion zone.To characterize the formation of the fusion zone, a historical indicator ξ was added in the system to emulate the phenomenological fusion of the strut during the multilayer SLS process following our former work [44,40], which is initialized as zero and irreversibly transitions as one once T ≥ T M .In Figs.7b 1 -7b 9 the fusion zones are denoted for selected P and v.The map can be thereby divided into characteristic regions where a continuous/discontinuous fusion zone or no fusion zone is formed.Increasing U from 17.01 to 22.68 J mm −3 , the transition from microstructures without fusion zone to with continuous fusion zone also implies the switch of densification mechanism from pure solid-state sintering to partial-melting-induced fusion (similar to the liquid-state sintering).This may explain the change of the tendency of − ln(1 − ϱ) when U < 22.68 J mm −3 in Fig. 7d.

Phenomenological relations for residual stress and plastic strain control
Taking the PB-averaged residual von Mises stress σP e and effective plastic strain pP e (defined in Eq. (1)) at the end of the simulations, Fig. 8 presents the distributions of σP e and pP e w.r.t. the P and v in the chosen processing window.Generally, the rise in σP e corresponds to an increase in the specific energy input U (i.e., increase P with maintained v or increase v with maintained P), resulting in further concentrated residual stress around concave features (surface depressions and particle sintering necks) and across layers, as already presented in Fig. 5a and 5b.Meanwhile, as the locally concentrated residual stresses and plastic strain are evident in the processed powder bed, understanding how these local quantities are distributed and how much effect the extreme values (here the local maximums ) impose on the evaluated averages σP e and pP e becomes essential.In that sense, we sampled all nodal points in the processed powder bed and performed statistics on σ e and p e of these points.The probability distribution of local σ e and p e are presented in Supplementary Fig. 1a.The third quartile (Q3) is employed to characterize the concentration of the local quantities.Then, P-v regions with Q3(σ e ) and Q3(p e ) correspondingly beyond the referencing values, which are PB-averages σ P e0 = 247 MPa and p P e0 = 0.012 from the case P = 20 W and v = 100 mm s −1 , are denoted in Fig. 8a and 8b (the complete P-v maps for Q3(σ e ) and Q3(p e ) are shown in Supplementary Fig. 2).These can be statistically interpreted as that the processed powder bed with P and v selected within the regions have at least 25% local substances obtaining σ e and p e beyond the referencing σ P e0 and p P e0 , respectively, while the rest of local substances have adequate σ e and p e up to these referencing values.In other words, the SLS processes with P and v selected from these regions tend to have processed powder bed with relatively higher localized residual stresses and accumulated plastic strains.
To understand the relationship between σP e and U, we interpreted σP e as the stored mechanical energy density after the multilayer SLS processes, which can be regarded as the residue of the energy density imported via the beam-induced thermal effect in the powder bed after all types of in-process dissipation.In this regard, the nonlinear regression analysis of an energy conversion law σP ]} was conducted on the simulation data.The parameters U P th and σ P ∞ adopt the physical meanings as the volumetric energy input at stress-free state and the saturated stress at infinity energy input, respectively, and σ P ∞ /U P σ is the increasing rate (slope) of σP e vs. U at stress-free state, as shown in the inset of Fig. 8c.The result in Fig. 8c presents a high correlation (R 2 = 98.41%) between simulated σP e vs. U with narrow confidence interval, demonstrating the applicability of proposed energy conversion law in predicting the residual stress in the powder bed.
Unlike the residual stress, the effective plastic strain p e is a measure of cumulative plastic deformation at any given moment during the process and lacks a clear physical picture of its relationship with volumetric energy input U. Therefore, the nonlinear regression analysis of an allometric scaling law pP e = C P p (U) I P p with parameters C P p and I P p was conducted to phenomenologically relate pP e to U, as shown in Fig. 8d.The analysis gives a relatively low correlation coefficient R 2 = 87.4% of the pP e (U) relation compared with one of σP e (U), with an expanded confidence interval at high U range.Regressed I P p = 1.17 ± 0.13 indicates an almost linear scalability of accumulated plastic strain in the processed powder bed on energy input via beam scan.It is worth noting that the proposed scaling law can be challenged as U appears to not being able to uniquely identify pP e , in other words, U of certain value can correspond to multiple p e value, as also depicted in Fig. 8d.Nonetheless, this scaling law can find it feasible in estimating strength of in-process plastification of the microstructure at the given specific energy input.
Since both pP e and σP e consider a complete powder bed with unfused particles, we further defined two average quantities that only take the residual stress and effective plastic strain inside the fused strut into account, denoted respectively σS e and pS e , which are calculated as σS e = Ω ξσ e dΩ Ω ξdΩ , pS e = Ω ξ p e dΩ Ω ξdΩ (5) with the fusion zone indicator ξ.Fig. 9 presents the distribution of σS e and pS e w.r.t.P and v in the chosen processing window.Contrasting with those shown in Fig. 8, a notable distinction in the maps presented in Fig. 9c and 9d is the emergence of regions where the selected P and v fail to form continuous fused strut.Meanwhile, the profiles of σS e and pS e receive significant influences from the size of the strut.Comparing Fig. 9a 1 with 9a 2 -9a 3 and 9a 4 -9a 5 , the increase of σS e follows the direction of increasing U as well, which simultaneously leads to the enlargement of the fused strut.As depicted in Figs.5c and 9a, concentrated residual stress is primarily found within the upper-layers of the strut.When the strut size is enlarged, it encompasses a larger volume with elevated maximum σ e and high-σ e region, resulting in an increased σS e with the rise in U. pS e exhibits similar tendency w.r.t.P and v as σS e , with an increase in U resulting in an increase in pS e .Nonetheless, regions with highly accumulated p e locates at the bottom of each layer's fusion zone, as illustrated in Figs.5d and 9b.At high U, such accumulation intensifies, especially at the bottom of the strut (which is also the bottom of the initial layer's fusion zone).Simultaneously, enlarged depth of a fusion zone indicates an extended remelting in the fused lower-layers, which removes some accumulated p e by the process or former layer within the strut, as shown in Fig. 6a.It also concentrates the high-p e region further to the bottom of the strut, notably the profile in Fig. 9b 4 .Eventually, the rise of pS e w.r.t.U is relatively less "rapid" than the one of σS e , evident in slightly sparser contours in Fig. 9d compared to 9c.Statistics of the local σ e and p e in the fused strut are also performed (Supplementary Fig. 1b), and the P-v regions characterizing Q3(σ e ) and Q3(p e ) beyond the referencing values, which are strut averages σ S e0 = 215 MPa and p S e0 = 0.016 from the same case P = 20 W and v = 100 mm s −1 , are denoted in the Fig. 9c and 9d (the complete P-v maps for Q3(σ e ) and Q3(p e ) are shown in Supplementary Fig. 3).As many local points outside the fusion zone are excluded from the statistics, extended regions of Q3(σ e ) > σ S e0 and Q3(p e ) > p S e0 are evident, demonstrating that more P-v combinations can lead to the fused strut containing at least 25% local σ e and p e are beyond the averages σ S e0 and p S e0 of the same chosen reference case.One can also conclude that most formed continuous struts more or less have localized residual stresses and accumulated plastic strains, based on the observation of narrow un-marked P-v regions between one with the discontinuous strut and one with Q3 of local quantities beyond the chosen references.In this sense, secondary thermal treatment should be positively considered to anneal the residual stresses inside the SLS-processed struts.
Nonlinear regression analyses were also conducted on σS e and pS e employing the proposed energy conversion law for residual stress and scaling law for effective plastic strain.Results are correspondingly presented in Fig. 9e and 9f.Notably, comparing to those of σP e (U) and pP e (U), the correlation of σS e (U) and pS e (U) present decline, respectively, with enlarged confidence interval at both low and high U ranges.This attributes to the removal of the influences from the unfused particles.Moreover, for σS e (U), regressed threshold energy input presents a negative value as U S th = −11.11± 17.43 J mm −3 , indicating a required energy output to achieve stress-free state.This is explainable as the σ e already exists in a just-formed strut.In other word, primary SLS process (i.e., the SLS without subsequent thermal post-process to release residual stress) cannot achieve samples with stress-free strut.Regressed saturated stress σ S ∞ = 260.05± 24.5 MPa also presents decline comparing to the regressed one on PB-averaged data.For pS e (U), a reduced regressed innndex of the scaling law pS e (U), i.e., I S p = 0.30 ± 0.06 in Fig. 9f, demonstrates a sublinear scalability of pS e (U) compared with one of pP e (U), which is almost linear.It also reflects a reduced growth rate of pS e at high U, and the underlying intensified remelting, which removes some accumulated p e within the strut while further concentrates p e around the strut bottom, and the comparably faster increase in strut size shall be the reason.Nonetheless, information conveyed by σS e (U) and pS e (U) is more relevant for practical application, as unfused particles shall be removed during post-process of a practical SLS.Further examination and validation of the proposed laws for residual stress and plastic strain w.r.t.specific energy input are expected in the future numerical and experimental studies.

Conclusion
In this work, we proposed a powder-resolved multilayer simulation scheme for producing porous materials using selective laser sintering, combing FEM-based non-isothermal phase-field simulation and thermo-elasto-plastic simulation with temperaturedependent material properties.This work has presented the mesoscopic evolution of stress and plastic strain on a transient thermo-microstructure under various beam power (P) and scan speed (v).Process-property relationships between porosity, residual stress and effective plastic strain and the volumetric energy input (U ∝ P/v) have also been demonstrated and discussed.The following conclusions can be drawn from the present work: (i) We proposed in this work a new powder-resolved mesoscopic residual stress formation mechanism for porous materials manufactured by the SLS process, which lead to the structural distortion appeared in fused strut.It was demonstrated with simulation results that the stress concentration at the necking region of the partially melted particles and inter-layer region between different layers provide dominant accumulated plastic strain and residual stress in the porous material.
(ii) Based on the proposed residual stress formation mechanism, we examined the proposed phenomenological relation between the porosity (densification) and the volumetric energy input U. Regression analysis on the resulting porosity from multilayer SLS simulations suggested an improved coherence with the experimental data, as the regressed densification coefficient K ϱ = 0.016 ± 0.001 mm 3 J −1 in this work is compared to the experimental K ϱ = 0.013 mm 3 J −1 and the one from our former single-layer simulation K ϱ = 0.019 mm 3 J −1 .
(iii) Two types of average quantities, namely PB-averaged ( σP e and pP e ) and strut averaged ones ( σS e and pS e ), were defined to characterize the residual stress and plastic strain within the powder bed and fused strut, correspondingly.The relationships between these quantities and volumetric energy input (U) are unveiled by conducting nonlinear regression analyses.The average residual stress ( σP e and σS e ) relates to U by the energy conversion law, while the average effective plastic strain ( pP e and pS e ) relates to U by the allometric scaling law.(iv) Attributing to the removal of influences from unfused particles, the correlation of the relations σS e (U) and pS e (U) present drops compared with one of σP e (U) and pP e (U), respectively.Saturation behavior is observed on both σP e (U) and σS e (U), while the linear scalability in pP e (U) degenerates into sublinear one in σS e (U), demonstrating a reduced growth rate of pS e at high U.
Despite the feasibility of the multilayer simulation scheme in recapitulating mesoscopic formation of porosity, residual stress and plastic strain under given processing parameters; and the proposed mechanism in explaining the structural distortion of SLS-produced samples, several points should be further examined and discussed in future works: (i) The present work omits the consideration of chronological-spatial distribution of the thermal-elasto-plastic properties among polycrystals, as the properties such as elasticity and crystal plasticity vary from grain to grain with distinct orientations.
(ii) The present findings are examined at relatively low specific energy input U and, correspondingly, low generated residual stress and accumulated plastic strain.Pore formation is also limited to lack-of-fusion mechanism.It is anticipated to conduct further simulations with relatively high energy input, where the mechanisms like keyholing co-exist with high residual stress and plastic strain.Extension of the proposed mechanism into the high-U range together with the following examination and validation are also expected.

Non-isothermal phase-field model for multilayer SLS processes
The non-isothermal phase-field model employed in this work to simulate the multilayer is based on our former works [7,47,39,40].For the sake of completeness, in this section we summarize the essentials of the employed model.For clarity, vectors are denoted by italic symbols with accented arrow (e.g., ⃗ v).Bold symbols (e.g., K, M, σ, ε) are for 2 nd -order tensors, and blackboard bold symbols (e.g., C) are for 4 th -order ones.
The model adopts both conserved and non-conserved OPs for representing the microstructure evolution of polycrystalline material.The conserved OP ρ indicates the substance, including the unfused and partially melted regions, while the nonconserved OPs {η i }, i = 1, 2, . . .distinguish particles with different crystallographic orientations.To guarantee that the orientation fields ({η i }) get valued only in the substance (ρ = 1), a numerical constraint ∑ i η i + (1 − ρ) = 1 is imposed to the simulation domain.This constraint also applies to the substrate where ρ = 1 always holds.
The thermo-structural evolution is governed by where the local free energy density is formulated as Here f (T, ρ, {η i }) contains multiple local minima, representing the thermodynamic stability of the states, i.e., substance, atmosphere/pore and grains with distinct orientations, with the coefficients W sf (T) and W gb (T) related to the barrier heights among minima, varying with temperature [7,40].The heat contribution f ht thereby manifests the stability of states by shifting the minima according to the local temperature field.c r is a relative specific heat landscape formulated by the volumetric specific heats in the substance c ss and pore c at , i.e., c r = h ss c ss + h at c at with corresponding interpolation functions h ss , h at indicating the substance (including solid and liquid) and the atmosphere/pore, respectively.When the material reached melting point T M , the extra contribution due to the latent heat L M is mapped by the interpolation function h M , which approaches unity when T → T M [7,40].It is worth noting that W sf (T), κ sf (T), W gb (T), and κ gb (T) inherit their temperature dependencies from the surface energy γ sf (T) and grain boundary energy γ gb (T), respectively (see Sec. 4.3).At the equilibrium, they can be related to γ sf (T), γ gb (T), and a diffuse-interface width of the grain boundary l gb , which have been explained in our former work [7,40].The Cahn-Hilliard mobility employed in Eq. ( 6) equation adopts the anisotropic form.Contributions considered in this work contain not only the mass transfer through substance, atmosphere, surface and grain boundary, but also the diffusion enhancement due to possible partial melting [47,40,48], i.e.
where D path is the effective diffusivity of the mass transport via path = ss, at, sf, gb, and h ss , h at , h sf and h gb are again interpolation functions to indicate substance, atmosphere/pore, surface and grain boundary, respectively, which obtain unity only in the corresponding region [47,40,48].Mass transport along surface and grain boundary is also regulated by corresponding projection tensor T gb and T sf .It is worth noting that the partial melting contribution M M is treated as an enhanced surface diffusion when T → T M due to the assumption of a limited melting phenomenon around the surface of particles.In this sense, formulations in Eqs. ( 6) and ( 10) disregard the contributions from melt flow dynamics as well as from inter-coupling effects between mass and heat transfer, i.e., Soret and Dufour effects [47].The isotropic Allen-Cahn mobility is explicitly formulated by the the temperature-dependent GB mobility G gb , GB energy γ gb and gradient coefficient κ η as [49,50] The phase-dependent thermal conductivity tensor takes the continuity of the thermal flux along both the normal and tangential directions of the surface into account, and is formulated as [51,52] where K ss and K at are the thermal conductivity of the substance and pore/atmosphere.N sf is the normal tensor of the surface [51,40,48].It is worth noting that radiation contribution via pore/atmosphere is considered in K at as K at = K 0 + 4T 3 σ B ℓ rad /3 with K 0 the thermal conductivity of the Argon gas and σ B the Stefan-Boltzmann constant.ℓ rad is the effective radiation path between particles, which usually takes the average diameter of the powders [53].The beam-incited thermal effect is equivalently treated as a volumetric heat source with its distribution along the depth direction formulated in a radiation penetration fashion, as the powder bed is regarded as an effective homogenized optical medium, i.e.
in which the in-plane Gaussian distribution p xy with a moving center ⃗ r O (⃗ v, t).P is the beam power and ⃗ v is the scan velocity with its magnitude v = |⃗ v| as the scan speed, which are two major processing parameters in this work.The absorptivity profile function along depth da/dz which is calculated based on Refs.[7,54].It is obvious that when one ignores the effects of those latent heat induced by microstructure evolution (i.e. the evolution of the pore/substance as well as unique grains), Eq. ( 8) can be degenerated to the conventional Fourier's equation for heat conduction with an internal heat source.

Elasto-plastic model for thermo-mechanical analysis
The transient temperature field and substance field ρ, resulting from the non-isothermal phase-field simulations, are then imported into the quasi-static elasto-plastic model as the thermal load as well as the phase indicator for the interpolation of mechanical properties.The linear momentum equation reads where σ is the 2 nd -order stress tensor.Taking the Voigt-Taylor interpolation scheme (VTS), where the total strain is assumed to be identical among pores/atmosphere and substance, i.e. ε = ε ss = ε at .[55,56,57].In this regard, the stress can be eventually formulated by the linear constitutive equation where the 4 th -order elastic tensor is interpolated from the substance one C ss and the pores/atmosphere one C at , i.e., In this work, C ss is calculated from the temperature-dependent Youngs' modulus E(T) and poison ratio ν(T), while C at is assigned with a sufficiently small value.The thermal eigenstrain ε th is calculated using interpolated coefficient of thermal expansion, i.e., Here I is the 2 nd -order identity tensor.Meanwhile, for plastic strain ε pl and isotropic hardening model with the von Mises yield criterion is employed.The yield condition is determined as with where σ e is the von Mises stress.s is the deviatoric stress, and σ y is the isotropic yield stress when no plastic strain is present.The isotropic plastic modulus H can be calculated from the isotropic tangent modulus E t and the Young's modulus E as H = EE t /(E − E t ).p e is the effective (accumulated) plastic strain, which is integrated implicitly by the plastic strain increment ε pl employing radial return method and will be elaborated in Sec.4.4.

Boundary conditions and model parameters
Together with the governing equations shown in Eqs. ( 6)-( 8), the following BCs are also employed for the process simulations with the convectivity h, Stefan-Boltzmann constant σ B , the hemispherical emissivity ε, and the pre-heating (environmental) temperature T 0 .Γ T and Γ B are correspondingly the top and bottom boundaries of the simulation domain, and Γ S is the set of all surrounded boundaries.Γ B ∪ Γ S ∪ Γ T = Γ, as the schematic shown in the inset of Fig. 1b.⃗ n is the normal vector of the boundary.Eq. (18) physically shows the close condition for the mass transfer, restricting no net mass exchange with the environment.Eq. (19) shows the heat convection and radiation, which are only allowed via the pore/atmosphere at the boundary (masked by h at ). Eq. (20) emulates a semi-infinite heat reservoir under the bottom of the substrate with constant temperature T 0 , consistently draining heat from the simulation domain and reducing its temperature back to T 0 .As summarized in Sec.4.1, this simulation requests following parameters/properties: the thermodynamic parameters W sf , W gb , κ ρ , and κ η ; the kinetic properties (diffusivities and GB mobility) D path (path = ss, at, sf, gb) and G gb ; and the thermal properties K ss and K at .Among them, W sf , W gb , κ ρ , and κ η are parameterized by the temperature-dependent surface as well as GB energies, and a nominal diffuse interface width l gb as (21) with normalized tendencies τ sf (T) and τ gb (T) that reach unity at and κ η = κ η τ gb (T).Constants W sf , W gb , κ ρ , κ η also satisfy a constraint (W sf + W gb )/κ ρ = 6W gb /κ η derived from the coherent diffuse-interface profile at equilibrium (Supplementary Note 1 of Ref. [7]).These parameter/properties are collectively shown in Table 1.
For the thermo-elasto-plastic calculations, a 250 × 250 × 250 µm 3 subdomain was selected from the center of the prime domain to eliminate the boundary effects (Fig. 1a).The above momentum balance is subjected to the following rigid support BC as restricting the displacement ⃗ u along the normal direction of all boundaries except the top (Γ T ), which is traction free.As summarized in Sec.4.2, temperature-dependent E, ν, α, E t , and σ y are as listed in Table 2, where piecewise linear interpolation was employed to implement their temperature dependence.

Numerical implementation and parallel computing
The theoretical models are numerically implemented via the finite element method within the program NIsoS, developed by authors based on MOOSE framework (Idaho National Laboratory, ID, USA) [58].8-node hexahedron Lagrangian elements are chosen to mesh the geometry.A transient solver with preconditioned Jacobian-Free Newton-Krylov method (PJFNK) and backward Euler algorithm was employed for both problems.
For non-isothermal phase-field simulations, the Cahn-Hilliard equation in Eq. (6) was solved in a split way, The constraint of the order parameters was fulfilled using the penalty method.To reduce computation costs, h-adaptive meshing and timestepping schemes are used.The additive Schwarz method (ASM) preconditioner with incomplete LU-decomposition subpreconditioner was also employed for parallel computation of the vast linear system, seeking the balance between memory consumption per core and computation speed [59].Due to the usage of adaptive meshes, the computational costs vary from case to case.The peak DOF number is on the order of 10,000,000 for both the nonlinear system and the auxiliary system.The peak computational consumption is on the order of 10,000 core-hour.More details about the FEM implementation are shown in the Supplementary Note 4 of Ref. [7].
For thermo-elasto-plastic simulations, a static mesh was utilized to avoid the hanging nodes generated from h-adaptive meshing scheme.In that sense, the transient fields T and ρ of each calculation step were uni-directionally mapped from the non-isothermal phase-field results (with h-adaptive meshes) into the static meshes, assuming a weak coupling between thermostructural and mechanical problems in this work.This is achieved by the MOOSE-embedded SolutionUserObject class and associated functions.The parallel algebraic multigrid preconditioner BoomerAMG was utilized, where the Eisenstat-Walker (EW) method was employed for determining linear system convergence.It is worth noting that a vibrating residual of non-linear iterations would show without the employment of EW method for this work.The DOF number of each simulation is on the order of 1,000,000 for the nonlinear system and 10,000,000 for the auxiliary system.The computational consumption is approximately 1,000 CPU core-hour.
A modified radial return method was employed to calculate the integral of the plasticity as well as to determine the yield condition during the process with the temperature-dependent elasto-plastic properties at any time t with the time increment ∆t.This method employs the trial stress σ ⋆ calculated assuming an elastic new strain increment ∆ε where the elasticity tensor is obtained under the current temperature field T t .Once the trial stress state is outside the yield condition (Eq.( 17)), i.e., the plastic flow exists, the stress is then projected onto the closet point of the expanded yield surface with the normal direction determined as n y = 3s ⋆ /2σ ⋆ e with the von Mises σ ⋆ e and deviatoric trial stress s ⋆ [60,61] .Meanwhile, assuming isotropic linear hardening under T t of every timestep, the amount of the effective plastic strain increment ∆p for returning the stress state back to the yield surface is calculated in an iterative fashion adopting Newton's method where p t and ∆p t are updated as p (t+∆t) and ∆p (t+∆t) at the end of the timestep (t + ∆t).G(T t ) and H(T t ) are shear and isotropic plastic modulus calculated under the local temperature at the current timestep (T t ).Here G(T t ) = E(T t )/[2 + 2ν(T t )] with E(T t ) and ν(T t ) correspondingly the temperature-dependent Young's modulus and Poisson ratio.Knowing that the plastic strain increment ∆ε pl = ∆pn y following the normality hypothesis of plasticity [60], the updated stress and plastic strain at the end of timestep are thereby obtained by in which the vanishing of the stress as well as the plastic strain beyond a stress-free temperature T A (in this work T A = T M ) is also considered.Eqs. ( 23)-( 27) are sequentially executed on every timestep and update the quantities under T t .It is worth noting that the linear-interpolated H is implemented for both pand T-dependence with where f ( Ti , p t ) is a piecewise function with the grids Ti (i = 1, 2, ...) and T i is inside the section bound by Ti and Ti+1 , i.e., T t ∈ ( Ti , Ti+1 ].Considering reduction of the non-linearity, the p-independent H(T) was practically employed in the simulations, as formulated in Eq. (24).          ) effective plastic strain p e in the distorted fused strut with varying beam power P and scan speed v, which are marked as points in the P-v maps of (c) strut-averaged residual stress σS e and (d) strut-averaged plastic strain pS e , respectively.The dotted lines maps represent the volumetric energy input U isolines.P-v regions with third quantiles Q3(σ e ) and Q3(p e ) beyond the referencing strut-averages σ S e0 = 215 MPa and p P e0 = 0.016 of the the processed microstructure under P = 20 W and v = 100 mm s −1 are also denoted along with ones indicating different geometries of fused struts.Nonlinear regressions of (e) σS e and (f) pS e on U with the regression parameters indicated correspondingly.

Fig. 7 .
Fig. 7. Dependence of microstructure on processing parameters.(a) Segment-wise porosity of the sample produced under various processing parameters along (a 1 ) building direction (BD) and (a 2 ) perpendicular to scan direction (SD), where the range of the substrate as well as the beam spot size (D FWHM ) are denoted.Representative height (H R ) and width (W R ) for porosity calculation are also selected.(b 1 )-(b 9 ) SLS-processed microstructures (with fusion zone indicated in red) under varying beam power and scan speed, which are marked as points in the porosity P-v map (c).The dotted lines represent the volumetric energy input U isolines and the dash-dotted line represents the median porosity isoline (24.1%).(d) Phenomenological relation between densification factor ϱ (calculated from porosity) and U.

Fig. 8 .
Fig. 8. Dependence of average residual stress and effective plastic strain microstructure on processing parameters inside the powder bed (PB).P-v maps of (a) PB-averaged residual stress σP e and (b) PB-averaged plastic strain pP e .The dotted lines maps represent the volumetric energy input U isolines.P-v pairs for simulated profiles in Fig. 5a and 5b are denoted correspondingly.P-v regions with third quantiles Q3(σ e ) and Q3(p e ) beyond the referencing PB-averages σ P e0 = 247 MPa and p P e0 = 0.012 of the processed microstructure under P = 20 W and v = 100 mm s −1 are also denoted.Nonlinear regressions of (c) σP e and (d) pP e on U with the regression parameters indicated correspondingly.

ScanFig. 9 .
Fig.9.Dependence of average residual stress and effective plastic strain microstructure on processing parameters inside the fused strut.Simulated profiles of (a 1 )-(a 5 ) residual von Mises stress σ e and (b 1 )-(b 5 ) effective plastic strain p e in the distorted fused strut with varying beam power P and scan speed v, which are marked as points in the P-v maps of (c) strut-averaged residual stress σS e and (d) strut-averaged plastic strain pS e , respectively.The dotted lines maps represent the volumetric energy input U isolines.P-v regions with third quantiles Q3(σ e ) and Q3(p e ) beyond the referencing strut-averages σ S e0 = 215 MPa and p P e0 = 0.016 of the the processed microstructure under P = 20 W and v = 100 mm s −1 are also denoted along with ones indicating different geometries of fused struts.Nonlinear regressions of (e) σS e and (f) pS e on U with the regression parameters indicated correspondingly.

Table 1 .
[68]rial properties of the bulk SS316L used in the non-isothermal phase-field simulations.Here R represents the ideal gas constant.Activation energy is obtained from[68]while the prefix factor is estimated as unity at T M after normalization.** Estimated as 100D sf /2(W sf + W gb ). *