Finite Element Analysis of Thermal Stress and Thermal Deformation in Typical Part during SLM

Selective laser melting (SLM) constitutes an additive manufacturing (AM) method. Many issues such as thermal strain and macro-thermal deformation, which are caused by the thermal stress of different process parameters, are not clear. In this paper, an efficient and fast manufacturing simulation method was researched based on a moving heat source model and an elastoplastic theory of welding simulation, which was studied based on the thermodynamic coupling algorithm with a software-developed application for the SLM process. Subsequently, typical case results of thin and hollow plate part formation and the corresponding performances were simulated in detail. The results demonstrated that the effective thermal stress increased as the layer height increased from the surface layer to the substrate, while the thermal strain followed an approximate change rule. In addition, the stress was released from the underlying substrate when the support was removed. Moreover, the largest single axis plane stress was changed from tension to compression from the edge to the center, finally reaching equilibrium. In particular, maximum macro thermal deformation occurred at the printed support structure to the samples, displaying similar results in other locations such as the corners. Finally, the effectiveness of the simulation could be verified from the realistic printed part, which could provide proof for the quality prediction of the part that is actually forming.


Introduction
Due to the complex processing technology of thin and hollow plate parts, negative production and material loss occur through traditional machining. However, technology, along with the development of additive manufacturing (AM), have become the preferred forming process plan for small quantities of complex parts [1,2]. Selective laser melting (SLM) is a typical rapid forming technology [1,3]. However, in practice, the metal powder is spread layer by layer followed by laser melting in SLM, which can produce a molten pool with a low-sized area and high temperature. Additionally, the temperature gradient is produced from high heat as well as being simultaneously transferred to the surrounding powder or printed body formation quickly [3][4][5], which makes the molten pool rapidly cool within microseconds. In addition, the molten pool shrinks, and thermal stress with the surrounding solid is produced under the interaction of condensation. When the stress is significantly high, macroscopic geometry deformation is produced, which is caused by the warping and cracking of the printed parts due to thermal stress. Therefore, a high number of printed parts failed to be produced, resulting in high amounts of material waste. Consequently, in this paper, the thermal field, thermal stress, SLM technique through forming errors. Additionally, an error model evaluation simulation of parts, formed during melting and solidification shrinkage (0.1/0.2 mm), was developed.
According to the literature, the universal applicability of such results is so limited that the experimental research could not be regularly conducted during SLM. Simultaneously, low efficiency simulation has led to the negative overall effect of certain parts through poor predictions, especially during variable cross-section forming. As an example, a certain amount of stress concentration caused failure, but it was difficult to measure and simulate. Consequently, the efficient simulation of complex parts on forming quality has higher practical significance.
The research object of this paper was a typical thin and hollow plate part. The structure had a thin panel and a water cooling channel at the middle of the panel for heat dissipation, while a barbed water nozzle interface with the hose was provided with the panel. In traditional processes, it is also necessary to combine the milling with a traditional machining center, along with the welding of the upper and lower boards to manufacture this part. Therefore, quality control is difficult. Consequently, an enterprise will attempt to adopt an SLM process. However, due to the warping caused by thermal deformation, many problems still exist in actual production such as the high flatness errors of the panel, the imprecision of the cylindricity values of the two nozzles, and the variation in the center distance between the two nozzles. As the new process would still produce a high amount of waste, a thermodynamic coupling algorithm based on welding software was developed in this paper, while the forming process simulation of a part was conducted. The stress simulation results in the height direction and in the plane position of the sample were analyzed, along with the statistical deformation data. Finally, the validity of the algorithm was proven by the experimental sample results.

Heat Source
The laser source was a Gaussian heat source, with the maximum heat flux density that decays exponentially far away from the center, which can be described according to [17] q = q 0 exp(− 2r 2 ω 2 ) (1) where q is the heat flux density (J/m 2 ·s); r is the radial distance from the beam center; ω is the radius of the beam; and q 0 is the intensity of the beam at r = 0: where P is the power of laser (W) and A is the absorptivity of the powder material, which can be calculated as 1 minus the reflectivity of the material and the transmissivity of light (usually, the absorptivity is considered 0.7 for 316L stainless steel). Next, the heat source model was given a picture, as shown in Figure 1a. The scanning paths of the heat source are drawn in Figure 1b, which is the moving path of the laser source by layer to layer in different directions. Each layer path changed by 30 degrees, which cycled from 0 to 180 (0, 30, 60, 90, 120, 0, 30, etc.). Appl. Sci. 2019, 9,

Heat-Transfer Equation
According to the definition of the heat flux density per unit cross-sectional area of the object, where k is the thermal conductivity (W/m K); T Δ is the temperature difference from high to low (K); dr is the thickness of the material in the current layer (m). When Equation (4) is combined with Equation (3), the temperature difference can be used to calculate: Furthermore, the heat transfer equation is where ρ is the material density (kg/m 3 ), c is the specific heat capacity (J/kg K), and t is the interaction time.
The porosity of the powder can be calculated as bulk power bulk ρ ρ φ ρ − = . (7) Consequently, power bulk The thermal equilibrium equation satisfies the following classical 3D heat conduction equation, given by [17] where Q = q*(s*t) l (x, y, z, t) is the volumetric heat generation (W/m 3 ).

Thermal-Mechanical Coupling
The thermal strain is described in the case where the temperature field distribution is known [17]:

Heat-Transfer Equation
According to the definition of the heat flux density per unit cross-sectional area of the object, where k is the thermal conductivity (W/m K); ∆T is the temperature difference from high to low (K); d r is the thickness of the material in the current layer (m). When Equation (4) is combined with Equation (3), the temperature difference can be used to calculate: Furthermore, the heat transfer equation is where ρ is the material density (kg/m 3 ), c is the specific heat capacity (J/kg K), and t is the interaction time.
The porosity of the powder can be calculated as Consequently, The thermal equilibrium equation satisfies the following classical 3D heat conduction equation, given by [17] ρc where Q = q*(s*t) l (x, y, z, t) is the volumetric heat generation (W/m 3 ).

Thermal-Mechanical Coupling
The thermal strain is described in the case where the temperature field distribution is known [17]: Appl. Sci. 2019, 9, 2231 5 of 16 where ε is the strain, α e is the coefficient of thermal expansion (1/K), and T ref is the reference temperature at the previous layer. For the overall elastic strain,

Stress-Strain Equations
The material is in an elastic state as each layer of powder begins to melt with the laser beam. For simplified isotropic materials, the thermal strain equation can be obtained through the constitutive equations of elasticity theory [17]: where γ is the shearing strain; E, µ, and G are the modulus of elasticity, Poisson's ratio, and shear modulus, respectively. The thermal stress equation can be obtained through the thermal strain equation solution: As the melting temperature increases, the stress-strain increases, and the material enters a plastic state. According to the Mises yield criterion and the Prandtl-Resus flow criterion, the incremental constitutive equation is used to calculate the elastic-plastic strain: dε ii = 1 2G dσ ii + dλσ ii , dε jj ..., dε kk ...; dγ ij = 1 G dτ i j + 2dλτ ij , dγ jk ..., dγ ik ... dγ ij = 1 G dτ i j + 2dλτ ij , dγ jk ..., dγ ik ... (16) (i, j, k ∈ x, y, z) where dλ is the positive instantaneous constant and the effective stress is calculated as follows:

Initial Conditions and Boundary Conditions
The initial condition of room temperature and some attribute parameters distributed throughout at time t = 0 can be applied as where u, v, and w is the displacement in three directions-x, y, and z, respectively. Meanwhile, as the natural convection and surrounding radiation with the gases were considered, the boundary conditions of the simulation model can be defined. The heat flux at the boundary is equal to the laser heat flux at the boundary unit.
where n is the normal direction of the model surface, and w is the model border.
where h is the heat transfer coefficient at the model surface; T w is the temperature of the model surface; T f is the temperature of the surrounding gas. ε is the radiation coefficient of the black body, and σ is the Stefan-Boltzmann constant.
Next, according to the input laser parameters such as P, r, A, and v as well as the property parameters of the material at different temperatures T, such as k, c, ρ, α, E, µ, and G, the software established the stiffness matrix to solve the stress balance equation (Navier equation), which can be obtained through the displacement of any unit [u i , v i , w i ]. Through the Cauchy equation, the stress-strain value of any unit can be solved by Equations (13)-(18).

Material Properties of 316L
The basic properties and the thermal properties of 316L stainless steel powder and solid 316L at different temperatures are shown in Table 1. Meanwhile, the mechanical properties of printed 316L are anisotropic in practice [38,39]. However, it is very difficult to obtain the whole basic parameter data of the anisotropic materials of 316L and was conducted on simplification processing as isotropic for the thermodynamic simulation in this model.

Simulation Flow
The simulation methodology is presented in Figure 2. This simulation included two aspects. On one hand, the heat source definition, the heat transfer mode, and the thermal-mechanical coupling calculations were carried out, based on the finite element theory. On this basis, the simulation model was utilized for a layered slice, whereas the layer thickness data were added to the processing technology file. In addition to the scanning path development, a common path striped pattern was input to the processing technology file. On the other hand, process parameters such as the material properties, printing support, process equipment, and laser parameter settings could be set within the corresponding parameter window. Moreover, the reasonable value could be modified on the basis of previous simulation results.
calculations were carried out, based on the finite element theory. On this basis, the simulation model was utilized for a layered slice, whereas the layer thickness data were added to the processing technology file. In addition to the scanning path development, a common path striped pattern was input to the processing technology file. On the other hand, process parameters such as the material properties, printing support, process equipment, and laser parameter settings could be set within the corresponding parameter window. Moreover, the reasonable value could be modified on the basis of previous simulation results.

Simulation Parameters
According to the current resources, the simulation was performed with a PC (8 cores, 3.60 GHz CPU, and 32GB memory). The surface grid of the model was 0.5 mm. Additionally, the main laser parameters and process parameters are presented in Table 2.

Simulation Parameters
According to the current resources, the simulation was performed with a PC (8 cores, 3.60 GHz CPU, and 32GB memory). The surface grid of the model was 0.5 mm. Additionally, the main laser parameters and process parameters are presented in Table 2.

Simulation Runs
The original 3D model of the component was designed as presented in Figure 3a. It included a panel with a relatively flat surface, a water-cooling channel, an internal connection with the water nozzle as well as many external cooling fins on the plate. Based on the above thermal-mechanical algorithm and software platform simulation, the simulation produced a diagram, layer-by-layer, as presented in Figure 3b. This diagram contained a 3 mm support, and certain points were marked on the graph.

Simulation Runs
The original 3D model of the component was designed as presented in Figure 3a. It included a panel with a relatively flat surface, a water-cooling channel, an internal connection with the water nozzle as well as many external cooling fins on the plate. Based on the above thermal-mechanical algorithm and software platform simulation, the simulation produced a diagram, layer-by-layer, as presented in Figure 3b. This diagram contained a 3 mm support, and certain points were marked on the graph. Subsequent to the finite element thermal-mechanical coupling solution, the stress nephograms of several forming stages are presented in Figure 4. Following laser melting, the melt solidification was achieved, while thermal stress regions were macroscopically generated for different temperature gradients as the cooling progressed. The reason for stress occurrence was that different molecular motion rates under different temperature gradients can cause different orientations and spacing of Subsequent to the finite element thermal-mechanical coupling solution, the stress nephograms of several forming stages are presented in Figure 4. Following laser melting, the melt solidification was achieved, while thermal stress regions were macroscopically generated for different temperature gradients as the cooling progressed. The reason for stress occurrence was that different molecular motion rates under different temperature gradients can cause different orientations and spacing of the crystals, resulting in stress formation with the slowest heat transfer point at the center, along with tension at the side and pressure in the center. Through the effective stress cloud diagram, it was observed that the plane stress was evenly distributed throughout the material. Furthermore, the four stages of the effective stress nephogram and its ruler are demonstrated in Figure 4a-d. From the first to the fourth stage, the stress value gradually increased, and the layer-by-layer cyclical superposition effect of the thermal stress from the base layer to the upper layer was produced as a result. This could be obtained because the lower layer was continuously heated and cooled by the upper layer. Following the formation of several hundred layers, the normal stress distribution of the model could be observed at the fourth stage. The stress value of the support part exceeded the solid stress value of the part from the cloud image color distribution. In addition, it can be observed in Figure 4e that the stress was reduced after the substrate was removed. Furthermore, this stress, following the thermal stress release, occurred as residual stress. Specifically, the existence of the above stresses would inevitably lead to strain generation and macroscopic thermal deformation. The simulation results are thoroughly analyzed in Section 4.

Emulation Results and Discussion
According to the thermo-mechanical coupling algorithm and stress-strain elasticity theory simulations, relevant data statistics were produced. The main conclusions are as follows.

Emulation Results and Discussion
According to the thermo-mechanical coupling algorithm and stress-strain elasticity theory simulations, relevant data statistics were produced. The main conclusions are as follows.

Effective Stress and Strain
In the formation analysis of heat stress subsequent to different levels of height, each floor included four locations (Points 1-4, as shown in Figure 3b). As presented in Figure 5a, the printing started from the support layers to the height position of the plate layers at 4 mm, and the thermal stress gradually increased from the upper plate surface to the lower layer near the substrate as the layers accumulated. Since the heat transfer temperature gradient was the highest, it led to faster solidification and higher thermal stress at the base plate layer due to "metal solidification theory" [40,41]. As the other powders were formed through laser exposure, which caused the layer-by-layer circulation heating, the current printing layer produced a "micro annealing effect" of stress release before each layer. In contrast, its influence was lower than that of the stress-induced rapid solidification. Consequently, the bottom stress was the highest and the upward stress became slightly lower with the forming height. Moreover, water nozzles with barbs existed at the height of 4 mm. These were added with the new support, but the heat transfer medium was a panel of considerable temperature, which led to a reduction in the part cooling rate and the stress decrease in the curves. Following this, the curve entered the normal stress accumulation phase as before. Simultaneously, the thermal stress caused the thermal strain of the material, which constituted microcosmic thermal deformation [35]. From Figure 5a,b, it could be observed that the thermal plastic strain and thermal stress variations had the same tendency. Furthermore, the effective stress extension from 0 to 4 mm in the large panel structure that was printed from the substrate was analyzed (single support layer was accumulated on a single basis), while the increase occurred as follows: Through the least square method cubic polynomial fitting function curve, Through MATLAB calculations, P = 0.0196, −0.9584, 7.2315, and 289.6667 because Therefore, the thermal stress curve monotonically decreases from the substrate to the upper layer, and vice versa.

Stress Release in Removed Substrate
After the substrate processing was removed from the simulation, the Z negative direction constraint of the substrate was relieved. Under the action of stress, the inner bound crystal alignment and its grain boundary could be adjusted. Consequently, the stress could be released and the residual stresses of the parts could be formed. The effective stress distribution prior to and following the removal of the substrate could be observed from Point 1 to Point 2 in Figure 6a. It was quite clear that Furthermore, the effective stress extension from 0 to 4 mm in the large panel structure that was printed from the substrate was analyzed (single support layer was accumulated on a single basis), while the increase occurred as follows: Through the least square method cubic polynomial fitting function curve, Through MATLAB calculations, P = 0.0196, −0.9584, 7.2315, and 289.6667 because Therefore, the thermal stress curve monotonically decreases from the substrate to the upper layer, and vice versa.

Stress Release in Removed Substrate
After the substrate processing was removed from the simulation, the Z negative direction constraint of the substrate was relieved. Under the action of stress, the inner bound crystal alignment and its grain boundary could be adjusted. Consequently, the stress could be released and the residual stresses of the parts could be formed. The effective stress distribution prior to and following the removal of the substrate could be observed from Point 1 to Point 2 in Figure 6a. It was quite clear that the closer the substrate, the higher the stress release, and the closer the proximity to the stress-free state at the bottom. When the parts reached the nozzle, the new support stresses could be preserved. As a result, both the residual stress and the thermal stress displayed an opposite trend prior to and following the direct removal of the substrate. Given that the internal molecules of the matrix severely changed the state of the crystal due to high thermal stress, the crystal was the first to make changes when the new equilibrium mechanism was obtained. This phenomenon conformed to the theory of the microstructure phase transition [41]. Additionally, the residual strain change trend was similar to the residual stress trend presented in Figure 6b. It could be observed that the simulation was more successful for the actual thermodynamic coupling.

X/Y/Z Uniaxial Plane Stress
Furthermore, the plane state of the uniaxial stresses of X, Y, and Z axes were analyzed at the panel height of 0.5 mm. As presented in Figure 7a, the uniaxial stress of X alternated between the two sides, which was in accordance with the relationship between the acting and reacting forces. Therefore, the plane stress tension was offset and the plane stress balance was reached. Figure 7b presents a similarity to the thermal stress distribution. The reason for this stress alternation was that following the melting of a layer of powder, the edge position was cooled first, thus becoming the stress distribution boundary. The subsequent solidification position caused the tension to increase as the central formation was symmetrical to the tensile stress on both sides of the boundary from the edge to the internal cooling center. Moreover, the internal solidification contraction produced contraction pressure, while as the cooling was extended to the center, the equilibrium compressive stress was produced. As presented in Figures 7a,b, from the center of geometry, the tension was drawn to the pressure zone and subsequently to the edge. Specifically, as presented in Figure 7c, the Z direction thermal stress cloud edge stress was high and the center tension was low, which led to the parts' susceptibility to warping deformation.

X/Y/Z Uniaxial Plane Stress
Furthermore, the plane state of the uniaxial stresses of X, Y, and Z axes were analyzed at the panel height of 0.5 mm. As presented in Figure 7a, the uniaxial stress of X alternated between the two sides, which was in accordance with the relationship between the acting and reacting forces. Therefore, the plane stress tension was offset and the plane stress balance was reached. Figure 7b presents a similarity to the thermal stress distribution. The reason for this stress alternation was that following the melting of a layer of powder, the edge position was cooled first, thus becoming the stress distribution boundary. The subsequent solidification position caused the tension to increase as the central formation was symmetrical to the tensile stress on both sides of the boundary from the edge to the internal cooling center. Moreover, the internal solidification contraction produced contraction pressure, while as the cooling was extended to the center, the equilibrium compressive stress was produced. As presented in Figure 7a,b, from the center of geometry, the tension was drawn to the pressure zone and subsequently to the edge. Specifically, as presented in Figure 7c, the Z direction thermal stress cloud edge stress was high and the center tension was low, which led to the parts' susceptibility to warping deformation. The distributions of axial stress are presented in Figures 7d-f. σx expresses the stress variation law of the extended molten pool along the scanning direction, σy expresses the change law of thermal stress in the trajectory direction, and σz expresses the stress change principle along the layer thickness accumulation. Each group could be divided into three areas, as presented in Figures 7d-f. The stress in Area I was positive and the edge was cooled first, therefore producing high residual tensile stress. The stress in Area II was the maximum compressive stress of the opposite force. The stress in Area III was the pull force from the edge tensile stress to the center. From the graph of the σx/σy/σz comparison, the distribution curve regularity was the same, but the amplitude varied highly. Therefore, the maximum residual stress of 316L stainless steel under this process parameter was σxMax = 295 MPa, σyMax = 283 MPa, and σZMax = 305 MPa, while the residual stress was lower than the 316L yield strength of 310 Mp < σs < 500 MPa. Additionally, these values were lower than the tensile strength 600 MPa < σp < 750 MPa [42].

Thermal Deformation
The total deformation of thermal deformation was analyzed, following the height interval of 1 mm. At the four edge points (Edge 1-4), the middle of the water nozzle (Point 1-4) and the midpoint of the boundary line (Middle 1-2) were considered, respectively. Consequently, the data of each point The distributions of axial stress are presented in Figure 7d-f. σ x expresses the stress variation law of the extended molten pool along the scanning direction, σ y expresses the change law of thermal stress in the trajectory direction, and σ z expresses the stress change principle along the layer thickness accumulation. Each group could be divided into three areas, as presented in Figure 7d-f. The stress in Area I was positive and the edge was cooled first, therefore producing high residual tensile stress. The stress in Area II was the maximum compressive stress of the opposite force. The stress in Area III was the pull force from the edge tensile stress to the center. From the graph of the σ x /σ y /σ z comparison, the distribution curve regularity was the same, but the amplitude varied highly. Therefore, the maximum residual stress of 316L stainless steel under this process parameter was σ xMax = 295 MPa, σ yMax = 283 MPa, and σ ZMax = 305 MPa, while the residual stress was lower than the 316L yield strength of 310 Mp < σ s < 500 MPa. Additionally, these values were lower than the tensile strength 600 MPa < σ p < 750 MPa [42].

Thermal Deformation
The total deformation of thermal deformation was analyzed, following the height interval of 1 mm. At the four edge points (Edge 1-4), the middle of the water nozzle (Point 1-4) and the midpoint of the boundary line (Middle 1-2) were considered, respectively. Consequently, the data of each point deformation with height are presented in Figure 8. It could be observed that the deformation amount of the four edges was the highest, whereas the maximum deformation occurred at the height of 3 mm from the supporting layer to the solid layer δ Max = 0.068 mm (which could meet the tolerance of ±0.1 mm). The second deformation was at the surrounding edges where the minimum deformation was at the center of the part. Even if the part with the water nozzle increased in height, no large deformation occurred. It was clear that the deformation of the parts mainly occurred to handle the transition position of the large plane edge point, especially from the supporting layer to the solid layer.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 17 deformation with height are presented in Figure 8. It could be observed that the deformation amount of the four edges was the highest, whereas the maximum deformation occurred at the height of 3 mm from the supporting layer to the solid layer δMax = 0.068 mm (which could meet the tolerance of ±0.1 mm). The second deformation was at the surrounding edges where the minimum deformation was at the center of the part. Even if the part with the water nozzle increased in height, no large deformation occurred. It was clear that the deformation of the parts mainly occurred to handle the transition position of the large plane edge point, especially from the supporting layer to the solid layer.

Experimental Results and Analysis
In order to verify the simulation results, typical parts were formed with a concept machine in SLM, as presented in Figure 9a. Additionally, a sample with the same process parameters was printed that only reflected the panel structure, while the size was 8 × 10 × 5 (mm). Subsequently, on the third day after forming, the residual stress was tested with a PROTO MG2000 XRD detector (Mn target, β from 20° to −20°, and 2θ ≈ 152.80°), as presented in Figure 9b. In the experiments, the electrolysis corrosion method was utilized to measure the different heights, while the thickness of the striped layer was approximately 0.5 mm.  Figure 10 presents a comparison analysis diagram of the sample with the simulation and experimental data subsequently after forming. In order to unify the comparison range, the simulation

Experimental Results and Analysis
In order to verify the simulation results, typical parts were formed with a concept machine in SLM, as presented in Figure 9a. Additionally, a sample with the same process parameters was printed that only reflected the panel structure, while the size was 8 × 10 × 5 (mm). Subsequently, on the third day after forming, the residual stress was tested with a PROTO MG2000 XRD detector (Mn target, β from 20 • to −20 • , and 2θ ≈ 152.80 • ), as presented in Figure 9b. In the experiments, the electrolysis corrosion method was utilized to measure the different heights, while the thickness of the striped layer was approximately 0.5 mm.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 17 deformation with height are presented in Figure 8. It could be observed that the deformation amount of the four edges was the highest, whereas the maximum deformation occurred at the height of 3 mm from the supporting layer to the solid layer δMax = 0.068 mm (which could meet the tolerance of ±0.1 mm). The second deformation was at the surrounding edges where the minimum deformation was at the center of the part. Even if the part with the water nozzle increased in height, no large deformation occurred. It was clear that the deformation of the parts mainly occurred to handle the transition position of the large plane edge point, especially from the supporting layer to the solid layer.

Experimental Results and Analysis
In order to verify the simulation results, typical parts were formed with a concept machine in SLM, as presented in Figure 9a. Additionally, a sample with the same process parameters was printed that only reflected the panel structure, while the size was 8 × 10 × 5 (mm). Subsequently, on the third day after forming, the residual stress was tested with a PROTO MG2000 XRD detector (Mn target, β from 20° to −20°, and 2θ ≈ 152.80°), as presented in Figure 9b. In the experiments, the electrolysis corrosion method was utilized to measure the different heights, while the thickness of the striped layer was approximately 0.5 mm.  Figure 10 presents a comparison analysis diagram of the sample with the simulation and experimental data subsequently after forming. In order to unify the comparison range, the simulation  Figure 10 presents a comparison analysis diagram of the sample with the simulation and experimental data subsequently after forming. In order to unify the comparison range, the simulation data were only focused on 3-8 mm. The samples were tested on the substrates directly (substrate present) as well as subsequently with wire-electrode cutting (substrate removed), and the test data were drawn as a black curve, as presented in Figure 10a. The black solid line shows that the forming stress of the unremoved substrate was identical to the thermal stress trend in the simulation. Simultaneously, the residual stress trend of the stress release after the removal of the substrate was also consistent with the black dotted line trend. It can also be seen in Figure 10b that the values of the two sets of emulation stress deviation before and after the removal of the substrate were very small, and measurement errors of experimental stress might exist. However, the trend of stress release was consistent before and after the removal of the substrate. This proved the simulation algorithm's validity.

Residual Stress
Appl. Sci. 2019, 9, x FOR PEER REVIEW 14 of 17 data were only focused on 3-8 mm. The samples were tested on the substrates directly (substrate present) as well as subsequently with wire-electrode cutting (substrate removed), and the test data were drawn as a black curve, as presented in Figure 10a. The black solid line shows that the forming stress of the unremoved substrate was identical to the thermal stress trend in the simulation. Simultaneously, the residual stress trend of the stress release after the removal of the substrate was also consistent with the black dotted line trend. It can also be seen in Figure 10b that the values of the two sets of emulation stress deviation before and after the removal of the substrate were very small, and measurement errors of experimental stress might exist. However, the trend of stress release was consistent before and after the removal of the substrate. This proved the simulation algorithm's validity.  Figure 11 presents the macroscopic deformation of the forming object. Figure 11a presents the physical image of a set of parameters that were not simulated. Moreover, the thermal deformation was quite high (δMax ≈ 2.3 mm). In particular, severe warping occurred at the edge and four corners, causing forming failure. However, it could be validated that the higher deformation position was at the corner and the maximum deformation amount was at the support of the solidified parts. On the other hand, the thermal deformation data were obtained through the selection of a better set of parameters from the simulation that was within the yield strength range of the material. The physical part of the corresponding process parameters is presented in Figure 11b. The design requirements of the parts were basically attained, and the best result could be optimized through the simulation to obtain improved part formation.

Residual Deformation
(a) (b) Figure 11. Physical thermal deformation of forming: (a) serious deformation; (b) ordinary deformation.  Figure 11 presents the macroscopic deformation of the forming object. Figure 11a presents the physical image of a set of parameters that were not simulated. Moreover, the thermal deformation was quite high (δ Max ≈ 2.3 mm). In particular, severe warping occurred at the edge and four corners, causing forming failure. However, it could be validated that the higher deformation position was at the corner and the maximum deformation amount was at the support of the solidified parts. On the other hand, the thermal deformation data were obtained through the selection of a better set of parameters from the simulation that was within the yield strength range of the material. The physical part of the corresponding process parameters is presented in Figure 11b. The design requirements of the parts were basically attained, and the best result could be optimized through the simulation to obtain improved part formation.

Residual Deformation
Appl. Sci. 2019, 9, x FOR PEER REVIEW 14 of 17 data were only focused on 3-8 mm. The samples were tested on the substrates directly (substrate present) as well as subsequently with wire-electrode cutting (substrate removed), and the test data were drawn as a black curve, as presented in Figure 10a. The black solid line shows that the forming stress of the unremoved substrate was identical to the thermal stress trend in the simulation. Simultaneously, the residual stress trend of the stress release after the removal of the substrate was also consistent with the black dotted line trend. It can also be seen in Figure 10b that the values of the two sets of emulation stress deviation before and after the removal of the substrate were very small, and measurement errors of experimental stress might exist. However, the trend of stress release was consistent before and after the removal of the substrate. This proved the simulation algorithm's validity.  Figure 11 presents the macroscopic deformation of the forming object. Figure 11a presents the physical image of a set of parameters that were not simulated. Moreover, the thermal deformation was quite high (δMax ≈ 2.3 mm). In particular, severe warping occurred at the edge and four corners, causing forming failure. However, it could be validated that the higher deformation position was at the corner and the maximum deformation amount was at the support of the solidified parts. On the other hand, the thermal deformation data were obtained through the selection of a better set of parameters from the simulation that was within the yield strength range of the material. The physical part of the corresponding process parameters is presented in Figure 11b. The design requirements of the parts were basically attained, and the best result could be optimized through the simulation to obtain improved part formation.

Conclusions
In this paper, the simulation of laser additive manufacturing was studied through thermodynamic coupling and elastoplastic theory. The simulation results demonstrated that the thermal stress and its strain corresponded to the geometrical deformation, when the part was quickly heated and the material rapidly solidified. Through simulation case analysis, the main conclusions were as follows: 1.
According to the reasonable design simulation algorithm and process parameters, the simulation process of the SLM parts was completed. Consequently, the thermodynamic coupling algorithm was simple and effective, while the thermal stress and thermal strain data were calculated to draw the curve.

2.
For the simulated sample, the thermal stress increased from the substrate to the upper surface layer by layer. While the stress was reduced when the substrate was removed, the stress decreased from the substrate to the upper surface.

3.
For the forming plane, the uniaxial stress reached the tension and pressure equilibrium state.
Through the σ x /σ y /σ z comparison, the same three areas-Areas I, II, and III (in Figure 7d-f)-showed tensile and compressive stress, but the amplitude and breadth were highly different. Furthermore, the maximum residual stresses of the 316L stainless steel under these process parameters were σ xMax = 295 MPa, σ yMax = 283 MPa, and σ zMax = 305 MPa.

4.
The changes of strain and thermal stress trends were relative, while the macroscopic geometric deformation occurred. From the simulation results, it could be deduced that the thermal deformation and thermal stress were in positive proportion. The higher deformation occurred at the corner of the same plane and the maximum deformation amount was located at the support toward the solid parts.

5.
After the simulation, a set of optimizing technological parameters in the experiment was printed in the case model. The characteristics of the experimentally formed parts were consistent with the simulation results. It was observed that the simulation analysis could effectively predict the quality of the forming parts, which proved instructive to the practical production.