**Mathematical Model for Prediction and Optimization of Weld Bead Geometry in All-Position Automatic Welding of Pipes**

### **Baoyi Liao 1,2, Yonghua Shi 1,2,\*, Yanxin Cui 1,2, Shuwan Cui 1,2, Zexin Jiang <sup>3</sup> and Yaoyong Yi <sup>4</sup>**


Received: 3 September 2018; Accepted: 21 September 2018; Published: 25 September 2018

**Abstract:** In this study all-position automatic tungsten inert gas (TIG) welding was exploited to enhance quality and efficiency in the welding of copper-nickel alloy pipes. The mathematical models of all-position automatic TIG weld bead shapes were conducted by the response surface method (RSM) on the foundation of central composition design (CCD). The statistical models were verified for their significance and adequacy by analysis of variance (ANOVA). In addition, the influences of welding peak current, welding velocity, welding duty ratio, and welding position on weld bead geometry were investigated. Finally, optimal welding parameters at the welding positions of 0◦ to 180◦ were determined by using RSM.

**Keywords:** all-position automatic tungsten inert gas (TIG) welding; optimal welding parameters; response surface method (RSM); lap joint; weld bead geometry

### **1. Introduction**

With the fast development of heavy industries, all-position automatic tungsten inert gas (TIG) welding of pipes had been widely applied to industries such as shipbuilding, nuclear power, chemical industries, and natural gas transportation. Welding plays an important role in joining pipes [1]. Traditional all-position welding of pipes is manual arc welding, but it has low welding efficiency, poor stability, and high cost. Moreover, the manual operation of pipe welding is a challenge to the welder. In order to overcome the aforementioned shortcomings of manual arc welding, an all-position automatic TIG welding process has been developed in recent years. At present, there are four types of automatic TIG butt-welding equipment for pipelines: closed welding heads, open welding heads, orbital welding trolleys, and thick-walled narrow gap TIG welding heads.

In all-position automatic TIG welding of pipes, especially for vertical and overhead position welding, it is a challenge to control the liquid metal flow down from the molten pool under the action of gravity. Therefore, there is a need to adjust welding parameters to keep the droplet transition smooth. In order to obtain high-quality weld beads, selecting optimal weld parameters and controlling the weld bead profile are important. As is well known, weld bead shape has a noteworthy influence on the mechanical properties of lap joints. In all-position automatic TIG welding of pipes, the weld bead profile is impacted not only by the weld parameters such as weld background current, weld peak current, weld speed, weld voltage, but also by the welding position. Therefore, it is hard to obtain several optimal parameters to receive perfect weld bead profiles. To realize this goal, the connection between weld parameters and weld bead shapes need to be established. A lot of research has been undertaken to establish mathematical models to optimize parameters and obtain the relationship [2,3]. Xu et al. [4] optimized the narrow gap all-position gas metal arc (GMA) welding parameters by employing the response surface method (RSM) and building regression models to predict the parameters of weld bead geometry. Rao et al. [5] had obtained the influence of welding parameters and statistical models for the prediction of weld bead shapes in pulsed GMA welding. Koleva [6] researched the connection with electron beam welding parameters and weld bead geometry; moreover, the mathematical model can optimize welding parameters. Karthikeyan and Balasubramanian used RSM to optimize friction stir spot welding procedure parameters and to obtain maximum lap shear strength of the lap joint [7]. The statistical model was built at an invariable position in inchoate research. However, in all-position automatic TIG welding of pipes, the welding position should be used as the input parameter of the statistical model because the weld includes flat, vertical and overhead positions [8], and welding position impacts the weld bead profile.

RSM can be used to predict the weld bead shape and mechanical properties in the welding process [9–12]. The ultimate goal of this paper is to build a new all-position automatic TIG pipes' welding procedure and a mathematical model of weld bead geometry using RSM.

### **2. Materials and Methods**

The schematic diagram of the system for all-position automatic TIG welding of pipelines is shown in Figure 1. The system includes an iOrbital 5000 welding machine, a TOA77 welding torch, a control panel, two argon tanks, a welding fixture and a cool water circulation system. The welding process was as follows:


**Figure 1.** The schematic diagram of the system for all-position automatic tungsten inert gas (TIG) welding of pipes.

The workpieces used in this study are UNS (Unified Numbering System for Metas and Alloys) C70600 with dimensions of Φ57 mm × 300 mm × 2 mm, but there is a flaring end of the pipe and the diameter is 62 mm, as shown in Figure 2. The chemical compositions of the parent metal are

listed in Table 1. In this study, experiments were executed with one-sided welding with double-sided shaping without an opening groove and welding wire, because the experiments melt the parent metal as filler metals. There is an angle of 30 degrees between the torch and the pipe vertical direction. Before welding, two pipes without a groove were fixed with a gap of 0.5 mm, and they were then lap welded. To ensure the stability of the welding process and avoid weld joints from being influenced by impurities and surface oil contamination, the pipes should be cleaned before welding. Pure argon (99.9% purity) was selected as a shielding gas with a flow rate of 15 L/min on the surface of the pipes and a flow rate of 20 L/min in the pipes to prevent the base metal from oxidizing in the air and enhancing the quality of the weld bead. The tungsten electrode diameter was 3.2 mm and the taper was 60◦.

**Figure 2.** The schematic diagram of the lap joint.

**Table 1.** Chemical compositions of parent metal (wt%).


In Figure 3, the specimens for weld bead shape observation were cut from the cross-section of the pipes, sealed with bake lite resin, then ground with a series of emery papers (grit size 400, 600, 800, 1200, 1500, 2000 and 2500), polished with a 1.0-μm diamond paste, and etched with an etchant of 5 g ferric chloride + 16 mL hydrochloric acid + 60 mL ethanol. Finally, a high-dynamic camera model NSC1003 (New Imaging Technologies, Paris, France) was used to photograph the weld bead shape and calculate the weld width, weld depth and weld thickness by image processing software. Reducing measurement error is important to the experiments, and therefore each value was measured three times. Then the average of the three measured values was calculated as experimental data.

**Figure 3.** The cross-section of the weld bead. (weld width (W), weld depth (D), weld thickness (T)).

In this investigation, due to the working range of individual factor being extensive, a central composition rotatable four-factor, five-level factorial design matrix was selected. The design software Design-Expert (V 10. 0. 7, Stat-Ease, MN, USA) was used to establish the design matrix and process the experimental data.

The steps of this investigation are as follows:


### **3. Creating Mathematical Model**

Based on the initial experimental results, four very important welding parameters affecting the weld bead geometry are the welding peak current (I), welding velocity (V), welding duty ratio (d) and welding position (P). The duty ratio is the ratio of welding peak current time to the pulse period. In order to create a mathematical model to describe and forecast the weld bead geometry in all-position automatic welding, I, V, d, and P were chosen as input parameters. The schematic diagrams of the weld bead shape and welding position are shown in Figures 3 and 4, respectively.

**Figure 4.** Schematic diagram of welding position.

Before this experiment, the control variable method was used to find out the working range of input parameters. The working range was determined by the steady welding procedure and invisible weld defects. In order to obtain optimized welding parameters, the establishment of a test parameter matrix adopts the center composed design (CCD) in the regression design method. In the CCD design, the upper and lower limit values of the input parameters were coded as ± *β*, the value of *β* is dependent on the number of input parameters and *β* = (2*k*) 1/4, *<sup>k</sup>* is the number of input parameters. The upper level and lower level were coded as ±1 respectively and the center point was coded as 0. In this research, there are 4 input parameters, so *β* is 2. Therefore, the upper and lower limit values of the input parameters were coded as ±2. The codes values for intermediate levels can be calculated by Equation (1):

$$X\_i = 2[2X - (X\_{\max} + X\_{\min})] / \left[X\_{\max} - X\_{\min}\right] \tag{1}$$

In Equation (1), *Xi* is the desired code value of a variable *X*, and range of values of *X* is from *Xmin* to *Xmax*; *Xmax* is the upper limits and *Xmin* is lower limits of the variable *X* [13,14].

The factor levels and coded values have been listed in Table 2. The experimental design matrix (Table 3) includes 30 sets of coded situations and constitutes a full replication four-factor factorial design of 16 points, 8 star points and 6 center points [7].


**Table 2.** Important weld parameters and their levels for copper–nickel alloy pipe.


**Table 3.** Trial design matrix and response of copper–nickel alloy lap joints.

In order to better reflect the weld bead geometry, reasonable response parameters need to be selected. The response parameters were named weld width (W), weld depth (D), and weld thickness (T). Figure 3 shows the cross section of weld bead.

On the foundation of central composite design matrix, the value of input parameters, response parameters and the regression model can be built. The connection between measured response and the input parameters could be shown as *y* = *f* (*x*1, *x*2, ... , *xi*) + ε, where *y* is the value of measured response, *xi* is the value of the input parameter, and ε is the systematic error. *Y* is a power transformation of *y* [11,15], so the second-order polynomial can be expressed as Equation (2):

$$Y = b\_0 + \sum\_{i=1}^{4} b\_i \mathbf{x}\_i + \sum\_{i=1}^{3} \sum\_{j=i+1}^{4} b\_{ij} \mathbf{x}\_i \mathbf{y}\_j + \sum\_{i=1}^{4} b\_{ii} \mathbf{x}\_i^2 \tag{2}$$

where, *b*<sup>0</sup> is the average of the measured response, the regression coefficients such as *bi*, *bij* and *bii* depend on linear, interaction and squared terms of factors, respectively.

Almost all response surface method problems can be approximated by these polynomials, and the regression coefficients could be obtained by the least squares method.

The regression model adopts the method of stepwise regression. Firstly, the regression model eliminates the insignificant terms and calculates the regression coefficients until the significant terms and the lack-of-fit terms of the regression model meet the requirements of the regression model. Finally, the relational expressions about all-position automatic TIG welding of the pipe within the range of 0–180◦ were obtained as Equations (3)–(8), which show the relationship between weld geometry shape and input parameters.

The final equations in terms of coded factors are given as follows:

$$\text{W} = 3.07 + 0.52 \text{ I} - 0.32 \text{ V} + 0.73 \text{ d} + 0.96 \text{ P} + 0.19 \text{ I} \times \text{d} + 0.27 \text{ I} \times \text{P} + 0.22 \text{ d}^2 - 0.16 \text{ P}^2 \tag{3}$$

$$\text{D} = 0.64 + 0.12\,\text{I} - 0.087\,\text{V} + 0.16\,\text{d} + 0.20\,\text{P} + 0.047\,\text{I} \times \text{d} + 0.065\,\text{I} \times \text{P} + 0.034\,\text{d}^2 - 0.032\,\text{P}^2 \tag{4}$$

$$\text{T} = 2.36 + 0.042\,\text{I} - 0.17\,\text{V} + 0.17\,\text{d} + 0.13\,\text{P} - 0.098\,\text{I} \times \text{d} - 0.12\,\text{V} \times \text{d} - 0.11\,\text{P}^2 \tag{5}$$

The final equations in terms of actual factors are given as follows:

$$\begin{array}{l}\text{W}=100.19346-0.64732\text{ I}-0.10788\text{ V}-2.29476\text{ d}-0.20122\text{ P}+10.012846\text{ I}\times\text{d}+0.00197176\text{ I}\times\text{P}+0.00898719\text{ d}^2-0.000783063\text{ P}^2\end{array} \tag{6}$$

$$\begin{aligned} \mathbf{D} &= 22.87088 - 0.15965 \,\mathbf{I} - 0.028847 \,\mathbf{V} - 0.47648 \,\mathbf{d} - 0.050128 \,\mathbf{P} + 0.00031250 \,\mathbf{P} + 0.00031250 \,\mathbf{I} \\ &- 0.00311250 \,\mathbf{I} \times \mathbf{d} + 0.00047824 \,\mathbf{I} \times \mathbf{P} + 0.00134656 \,\mathbf{d}^2 - 0.0000160301 \,\mathbf{P}^2 \end{aligned} \tag{7}$$

$$\begin{aligned} \text{T} &= -62.80003 + 0.34011 \text{ I} + 0.34775 \text{ V} + 1.32831 \text{ d} + 0.012728 \text{ P} \\ &- 0.00652083 \text{ I} \times \text{d} - 0.00811250 \text{ V} \times \text{d} - 0.0000546879 \text{ P}^2 \end{aligned} \tag{8}$$

Analysis of variance (ANOVA) was used to determine the significance and suitability of the regression model. Tables 4–6 show the ANOVA analysis of the weld width, the weld depth and the weld thickness model respectively.


**Table 4.** Results of analysis of variance (ANOVA) for model of weld width.

**\*** Degree of freedom (df), a concept in statics, indicates the number of unconstrained variables in calculating a statistical magnitude. According to the usual definition, df = *n* − *k*, *n* is the number of samples and *k* is the number of constrained variables or conditional number, while *k* is also the quantity of the other independent statistical magnitude in calculating one statistical magnitude.


**Table 5.** Results of ANOVA for model of weld depth.

**Table 6.** Results of ANOVA for model of weld thickness.


Using Design Expert V 10.0.7 Software can calculate the value of coefficient and the significance of each coefficient was confirmed by Student's t test and *p* values. The values of "Prob > F" less than 0.0500 indicate model terms are significant and values greater than 0.1000 indicate the model terms are not significant [16,17]. Tables 4–6 show the result of ANOVA for the W model, D model and T models, respectively, and the models' F values are 30.86, 31.63 and 7.89. The probability of F (prob > F) is less than 0.0001, in other words, these models are significant. Sometimes, these models show that the test results of lack-of-fit are insignificant relative to the pure error, insignificant lack-of-fit represents that the quadratic model is adequate. According to Table 4, I, d and P are the most important factors of the W model; V, (I × d), (I × P), d<sup>2</sup> and P<sup>2</sup> also could affect W. From Table 5, I, d and P are the most important factors of the D model; V, (I × d), (I × P), d<sup>2</sup> and P<sup>2</sup> also could affect D. I, V, d, P, (I × d), (V × d) and P2 could affect T as shown in Table 6.

### **4. Verification of Models**

To assure that the established model can predict and control the weld bead shape in actual application, it should test the accuracy of the mathematical model. The test experiments were executed by assigning diverse values for experimental variables within their working limits, but distinguishing them from the values of the design matrix. The values of input parameters, predicted response, actual response and percentage errors are listed in Table 7 respectively. It shows that the percentage errors for any models are less than 9%, and all the percentage errors are within the scope of industrial engineering requirements. Therefore, the statistical models can predict and optimize weld bead shape.


**Table 7.** Predicted values and actual values of the weld bead geometry.

\*\* Percentage error = *actual value*−*predicted value predicted value* × 100.

### **5. Results and Discussion**

According to the all models, the prime and interaction influences of input weld parameters on weld bead geometry can be found.

### *5.1. Influences of Welding Peak Current on Weld Width (W), Weld Depth (D), and Weld Thickness (T)*

Considering the influence of the single factor welding peak current on the weld bead geometry, it can be shown from Figure 5 that the welding peak current increases with the increase in weld width, weld depth and weld thickness. This is because with weld peak current increasing leads to heat input increase per unit time, which is good for the weld metal melted and enhancing the deposition efficiency. As the peak current increases, the weld width and depth increase significantly, while weld thickness increases less. The weld thickness increases gradually from 2.331 to 2.436 mm with the increase in welding peak current from 114 to 126 A. This is because as the weld peak current increases, the arc force also increases, and the liquid metal is blown to both sides of the molten pool under the action of arc force, so the weld thickness increases indistinctly.

**Figure 5.** Influence of weld peak current on weld bead. (Welding velocity (V) = 63 mm/min, welding duty ratio (d) = 50%, welding position (P) = 90◦).

### *5.2. Influences of Weld Velocity on W, D and T*

Figure 6 indicates that the increase of welding velocity leads to the decrease in the heat input and volume of liquid filler metal per unit time of the weld bead, so the weld width, the weld depth and the weld thickness are correspondingly reduced.

**Figure 6.** Influence of welding velocity on weld bead. (Welding peak current (I) = 120 A, d = 50%, P = 90◦).

### *5.3. Effects of Duty Ratio on W, D and T*

As shown in Figures 7–9, as the duty ratio increases, the weld width, depth and thickness increase correspondingly. This is due to increases in the heat input and volume of the liquid filler metal per unit time of the weld bead. The weld width and depth have a quadratic parabolic relationship with the duty ratio, and the slope of the curve growth is getting bigger and bigger. Because the duty ratio is increased, the weld heat input to the weld bead is larger and the more heat that is accumulated per unit time, the more favorable is the melting of the weld metal so the rate of weld width and weld depth growth is increasing. It can be seen from Figure 9 that the weld thickness increases linearly with increasing duty ratio, although the increase in duty ratio and the accumulation of welding heat contribute to the weld thickness, but excessive heat leads to liquid metal loss on the weld surface.

**Figure 7.** Influence of duty ratio on weld width.

**Figure 8.** Influence of duty ratio on weld depth.

**Figure 9.** Influence of duty ratio on weld thickness.

### *5.4. Effects of Welding Position on W, D and T*

As shown in Figures 10 and 11, as the degree of welding position increases, the weld width and weld depth increase correspondingly in the welding interval of 0–180◦. The weld width and weld depth have a quadratic parabolic relationship to the degree of welding position, and the slope of the curve is getting smaller and smaller due to the gravity of the molten pool and the flow of molten metal along the weld bead during welding. Figure 12 shows the weld thickness increases first and then decreases. The gravity of the molten pool at different welding positions was shown in Figure 13. It was decomposed into a tangential force Gt and a radial force Gr. When welding in the 0–180◦ interval, the molten pool is subjected to the tangential force Gt, which leads to the molten liquid metal flowing down along the weld bead, and the flowing liquid metal can preheat the remaining weld bead and fill the weld bead, so the shape parameters of the weld width, weld depth and the weld thickness will increase with the increasing of the welding position degree. But when welding is in the 0–90◦interval, the molten pool is subjected to the radial force Gr, and Gr = G cosθ. As the degree of the welding position increases, the radial force becomes smaller and smaller, and the direction of the radial force points to the center of the pipe; In the 90–180◦ interval welding, the molten pool is subjected to the direction of the radial force back to the center of the pipe, and Gr = G sin(θ − 90); with the welding position degree increasing, Gr becomes larger and larger, hindering the increase in weld width, depth and thickness. Therefore, the slope of the curve growth is getting smaller and smaller. Figure 12 shows the peak value of weld thickness for all values of P is received when P is 116◦.

**Figure 10.** Influence of welding position on weld width.

**Figure 11.** Influence of welding position on weld depth.

:HOGLQJSRVLWLRQ3GHJUHH

**Figure 12.** Influence of welding position on weld thickness.

**Figure 13.** Gravity of molten pool at different welding positions.

### *5.5. Effects of Two-Factor Interaction on W, D and T*

The interaction of I and d are shown in Figures 14–16. The speed of the increase in W and D with the increases in I increases as d increases, and the speed of the increase in W and D with d increases as I increases. This is because that I and d have an active influence on heat input of the weld width and depth. In Figure 16, T increases as I increases when d is less than 50%, while it decreases as I increases when d is more than 50%.

The interaction of I and P are shown in Figures 17 and 18. The rate of the increase in W and D with the increases in I increases gradually as P increases, and the speed of increase in W and D with P increases as I increases.

According to Figure 19, when the welding velocity is small, the T increases with the increase of the duty ratio. As the welding velocity increases, the influence of the duty ratio on T becomes smaller and smaller, because the welding velocity is large and the time of the welding process is short, so the welding heat input and the accumulated welding heat is small. When the duty ratio is small, the welding velocity has little effect on T. When the duty ratio is greater than the critical value, T decreases as the welding velocity increases.

**Figure 14.** Interaction of peak current and duty ratio on weld width.

**Figure 15.** Interaction of peak current and duty ratio on weld depth.

**Figure 16.** Interaction of peak current and duty ratio on weld thickness.

**Figure 17.** Interaction of peak current and welding position on weld width.

**Figure 18.** Interaction of peak current and welding position on weld depth.

**Figure 19.** Interaction of weld velocity and duty ratio on weld thickness.

### *5.6. Optimization of the Welding Parameters by Numerical Method*

In order to obtain the ideal weld bead shape and avoid welding defects such as incomplete penetration and weld collapse, numerical analysis was used to optimize the welding parameters. The goal, lower, upper limits and importance for every input and response parameters of the standard are shown in Table 8. Table 9 lists optimal welding parameters at 0, 45◦, 90◦, 135◦ and 180◦ welding position. The Figure 20 has shown the cross sections of weld bead with the optimal parameters at disparate welding position in 0–180◦, it shows that the gap between tubes is different for each case; there are two factors lead to this phenomenon. One factor is welding deformation: the welding process generates a lot of heat, which lead to pipes deformation. The other factor is welding sequence: welding experiments weld the bead at the 0◦ position firstly, so the weld bead was solidified first at the 0◦ welding position, which results in other gaps between tubes being confirmed.




**Table 9.** Optimal parameters.

**Figure 20.** Cross sections of weld bead at different welding position from 0–180◦. (**a**) 0◦ welding position, (**b**) 45◦ welding position, (**c**) 90◦ welding position, (**d**) 135◦ welding position.

### **6. Conclusions**

In this research, the all-position automatic welding of pipes has been researched and statistically analysed. The main conclusions drawn from this research are as follows:


**Author Contributions:** B.L. performed the experiments, analyzed the data, and wrote the original manuscript. Y.S. conceived and designed the experiments, and contributed significantly to analysis and manuscript improvement. Y.C. and S.C. helped perform the analysis with constructive discussions. Z.J. and Y.Y. helped perform the experiments and analyse the data.

**Funding:** This research was funded by Science and Technology Planning Project of Guangzhou City (Grant No. 201604046026), Science and Technology Planning Project of Guangdong Province (Grant No. 2015B010919005), and National Natural Science Foundation of China (Grant No. 51374111).

**Acknowledgments:** The authors gratefully acknowledge Guangzhou Shipyard International Company Limited for providing experimental materials.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Volumetric Heat Source Model for Thermal Modeling of Additive Manufacturing of Metals**

### **Yabin Yang 1,\* and Xin Zhou <sup>2</sup>**


Received: 27 September 2020; Accepted: 15 October 2020; Published: 22 October 2020

**Abstract:** In additive manufacturing of metallic materials, an accurate description of the thermal histories of the built part is important for further analysis of the distortions and residual stresses, which is a big issue for additively manufactured metal products. In the present paper, a computationally volumetric heat source model based on a semianalytical thermal modeling approach is proposed. The proposed model is applied to model the thermal response during a selective laser melting (SLM) process. The interaction between the laser and the material is described using a moving volumetric heat source. High computational efficiency can be achieved with considerable accuracy. Several case studies are conducted to examine the accuracy of the proposed model. By comparing with the experimentally measured melt-pool dimensions, it is found that the error between the predictions obtained by the proposed model and the experimental results can be controlled to less than 10%. High computational efficiency can also be achieved for the proposed model. It is shown that for simulating the thermal process of scanning a single layer with the dimension of 2 mm × 2 mm, the calculation can be finished in around 110 s.

**Keywords:** additive manufacturing; thermal modeling; volumetric heat source; computational efficiency

### **1. Introduction**

Additive manufacturing (AM), defined by ASTM international [1] as "the process of joining materials to make parts from 3D model data, usually layer upon layer, as opposed to subtractive manufacturing and formative manufacturing methodologies", opens up a new era to design novel structural materials with complex geometries. Direct energy deposition and powder-bed fusion are the two main AM approaches for metallic materials [2]. For both approaches, a 3D object is usually first sliced into thousands of 2D layers. Focused thermal energy—such as the laser, electron-beam, or plasma-arc—is then employed on a working plane to melt and fuse the material (coming from wire or powder) locally with the designed path in accordance with the corresponding 2D cross-sectional layout. After one layer is finished, the working plane is lowered for a small distance and the heat source scans the subsequent slice. This process repeats until the complete 3D object is built. Support structures may sometimes be needed to eliminate overhanging. Finally, the support structures are removed and the built part is cut from the baseplate.

Residual stresses and distortions are major issues for the AM process of metals due to the complex set of heat-cooling cycles during the process [3,4]. Part distortions can disqualify precision components with tight dimensional tolerances, and high values of residual stresses can lead to part failure while being built. It is therefore of great interest to optimize the process parameters and part design in order to minimize the undesired residual stresses and distortions. One approach is experimental trial-and-error. A more systematic way to enhance the understanding of the interplay between the processes parameters, design, and final part quality is computational process modeling. Since the transient temperature response is the root cause of part distortions and residual stresses, a computationally efficient thermal model for the AM process of metals is the key step to increase the build quality and repeatability leading to products with superior mechanical properties.

Although the mechanisms for the two AM processes of direct energy deposition and powder-bed fusion are different, the modeling techniques for predicting the thermal transients during theses two AM processes in part scale share many common features [5]. In order to obtain certain computational efficiency, the powder or the wire is usually represented as a continuum having effective thermal properties [2], and the interaction between the focused thermal energy and the material can be simplified as a surface or volumetric heat source applied on the material [6]. The temperature transient is usually obtained by solving the Fourier's law of heat conduction with Neumann and Dirichlet boundary conditions (BCs) using the finite element method (e.g., [7–9]). Phase transitions within the melt pool can also be neglected as a second-order effect to improve the computational efficiency [10,11]. Even with these simplifications, the challenge of the thermal modeling of AM process in part scale is still how to accurately predict the thermal transient within a reasonable amount of time. The modeled part is usually orders of magnitude larger than the heat source, thus resulting in very fine discretizations in both space and time domains to describe the movement of the heat source [2]. It should be noted that even the heat source is modeled as a dimensionless moving point source, the multiscale nature of the problem and the numerical requirements still pose a computational challenge.

One way to address the separation of the scales in the problem is to use an adaptive mesh refinement in the vicinity of the heat source. However, this requires an additional remeshing step within time integration. Zhang et al. [12] developed an adaptive remeshing technique to reduce the computational cost for modeling the heat-transfer process in selective laser melting (SLM, belonging to the family of powder-bed fusion technique). Although the computational efficiency is certainly improved, the influence of the heat source (i.e., the laser) moving path is not fully considered. Instead, the powder-bed deposition is simplified by the scale of an entire layer or fractions of each layer, and each fraction is heated entirely for a certain effective time interval and then cools down. As the thermomechanical response in the metal AM process is very sensitive to the heat source moving path [13], an adaptive remeshing technique which takes into account the influence of the real heat source moving strategy may need to be further developed.

Yang et al. [14,15] proposed a semianalytical thermal approach based on the superposition principle and applied it in modeling the thermal response of the SLM process. The total temperature is decoupled as the superposition of an analytical field and a numerical field. The analytical field corresponds to the temperature caused by the moving heat source in a semi-infinite space, for which a closed-form expression exits. The numerical field is employed to account for the BCs and solved numerically. In the semianalytical approach, there is no need to give fine discretization for the heat source. The mesh size in solving the numerical field scales with the dimension of the modeled part and hence, the computational efficiency can be improved. In the model developed by Yang et al. [14,15], the moving heat source is discretized as a number of dimensionless point sources. Considerable accuracy of the thermal predictions can be achieved at a distance of 100 μm away from the point heat source [14]. However, predictions of the near-field temperature evolution, especially in the vicinity of the center of the melting zone, are less accurate [16,17]. Furthermore, in AM process, the nonaxisymmetric melt-pool observations with respect to the moving direction of the heat source suggest the heat energy distributes nonuniformly over a certain volume, and this cannot be captured by the dimensionless point heat source, for which the energy diffuses isotropically in the space domain. Therefore, it would be more appropriate to model the heat source as volumetric. Li et al. [18] compared the effect of using surface and volumetric heat sources in modeling the laser melting of ceramic materials and found the model incorporating the volumetric heat source increased the accuracy of melt-pool predictions.

In the present paper, taking the SLM process as an example, a volumetric heat source model based on the semianalytical thermal approach proposed by Yang et al. [14,15] is developed to predict the thermal histories of the built part in SLM. Compared to the point heat source model [14,15], the proposed volumetric heat source model is able to capture the nonaxisymmetric nature of the melt pool and give more accurate temperature predictions in the vicinity of the heat source. Meanwhile, the computational efficiency of the proposed volumetric heat source model is not impaired. The predicted melt-pool dimensions by modeling the heat source in SLM as point and volumetric sources are compared. The accuracy of the proposed volumetric heat source model is evaluated by comparing the corresponding simulation results with experiments. Section 2 introduces the volumetric heat source model based on the semianalytical thermal approach. Three case studies are investigated in Section 3 to evaluate the accuracy of the proposed model. The article concludes with a reiteration of the most salient points of the study.

### **2. Model Description**

The semianalytical thermal model developed by Yang et al. [14,15] is briefly introduced in Section 2.1 The SLM process is taken as an example to explain how the semianalytical approach is utilized to model the thermal transients. The volumetric heat source model is then detailed in Section 2.2.

### *2.1. Semianalytical Thermal Model*

Consider a 3D body *V* that has already been built on the baseplate, as shown in Figure 1a. The top, lateral, and bottom surfaces of the body *V* are represented by *∂V*top, *∂V*lat, and *∂V*bot, respectively. In SLM, a thin layer is laid on the top surface *∂V*top. The bottom surface *∂V*bot is bonded to the baseplate and the lateral surface *∂V*lat is in contact with the powder. Since the mean conductivity of the solid body *V* is much larger than that of the powder [19], it can be assumed that there is no heat transfer between body *V* and the powder, and hence, the lateral surface *∂V*lat is thermally insulated. The uppermost layer of powder is also neglected since its overall heat capacity is also negligible. The moving laser is modeled as a moving heat source enforced on the top surface *∂V*top. During SLM, the baseplate is usually preheated and keeps a relatively constant temperature [20,21]. This is considered as prescribing a fixed temperature *Tc* on the bottom surface *∂V*bot. Consequently, the temperature of body *V* caused by the moving heat source on *∂V*top is governed by the heat equation

$$\frac{\partial T}{\partial t} = a\nabla^2 T + \frac{\dot{Q}\_v}{\rho c\_p} \tag{1}$$

with the BCs

$$
\frac{\partial T}{\partial \mathbf{x}\_i} \mathbf{n}\_i = \mathbf{0}, \quad \text{on} \quad \partial V\_{\text{top}} \quad \text{and} \quad \partial V\_{\text{lat}}, \quad i = 1, 2, 3 \tag{2}
$$

$$T = T\_{\mathfrak{c}\prime} \quad \text{on} \quad \partial V\_{\text{bot}}.\tag{3}$$

The initial condition is

$$T = T\_{\text{ini}\prime} \quad \text{at} \quad t = 0. \tag{4}$$

The time is represented by *t* and the thermal diffusivity *α* = *k*/*ρcp*, where *k* is the conductivity, *ρ* is the density and *cp* is the specific heat. The heat generation rate is represented by *Q*˙ *v*. The Cartesian coordinate system is illustrated in Figure 1a. The normal to the boundary surfaces of *∂V*top and *∂V*top are denoted as *ni*, and *T*ini is the initial temperature.

It should be noted that Equation (1) is linear by assuming the thermal properties *k*, *ρ*, and *cp* are temperature independent. Although for most materials, the specific heat *cp* and conductivity *k* are both temperature dependent. However, it has been demonstrated by many studies [15,22,23] that by choosing appropriate effective thermal constants *cp* and *k*, the linear governing Equation (1) can still result in accurate prediction of the temperature field. The energy losses due to convection, radiation, and phase transitions are not explicitly accounted for and instead an effective laser power is later employed to implicitly characterize these energy losses.

**Figure 1.** (**a**) A laser scanning is applied on the top surface of the body *V*. (**b**) The scanning line is discretized by a finite number of heat sources. Image sources by mirroring the original heat sources are added with respect to the boundaries. In a semi-infinite space, the temperature caused by the original discretized heat sources is *T*˜, and the temperature caused by the image sources is *T*˘. (**c**) Temperature filed *T*ˆ is employed to account for the boundary conditions (BCs). The total temperature *T* can be decomposed as the superposition of *T*˜, *T*˘, and *T*ˆ.

Since the governing Equation (1) and the BCs expressed by Equations (2) and (3) are all linear, the temperature *T* can be decomposed as

$$T = \mathcal{T} + \mathcal{T} + \mathcal{T},\tag{5}$$

where *T*˜ is the analytical temperature field caused by the moving heat source in a semi-infinite space, for which the boundary surface coincides with the top surface *∂V*top of the body *V*, as shown in Figure 1b. By discretizing the moving heat source as a finite number of individual heat sources, as shown in Figure 1b, the analytical field *T*˜ can be expressed as

$$\mathcal{T} = \sum\_{I=1}^{N} \mathcal{T}^{(I)},\tag{6}$$

where *T*˜(*I*) is the temperature field caused by the *I*th heat source in the semi-infinite space and *N* is the total number of heat sources. The source density, which represents the number of heat sources per unit length, is given by

$$
\rho\_s = \frac{1}{v\Delta t'} \tag{7}
$$

where *v* is the heat source moving speed and Δ*t* is the duration between the activation of two consecutive sources.

In Equation (6), the *T*˘ and *T*ˆ fields are employed to account for the BCs. The *T*˘ is the image source field, which is the temperature field caused by the image sources in the semi-infinite space (see Figure 1b) and can also be obtained analytically, and *T*ˆ is the complementary field being solved numerically (see Figure 1c). The image sources are added by mirroring the original heat sources with respect to the boundary surfaces. As illustrated in Figure 2, image source *J* (1) <sup>1</sup> is added by mirroring the heat source *I* with respect to boundary *∂B*(1). Consider that the power associated with heat source *I* is *P*. If the power of image source *J* (1) <sup>1</sup> is *<sup>P</sup>*, the heat flux on boundary *<sup>∂</sup>B*(1) caused by heat source

*I* and image source *J* (1) <sup>1</sup> would be 0; while if the power of image source *J* (1) <sup>1</sup> is −*P*, the temperature on boundary *∂B*(1) caused by heat source *I* and image source *J* (1) <sup>1</sup> would be 0. Therefore, the no heat flux and prescribed temperature BCs shown in Equations (2) and (3) can be satisfied by adding the image sources. However, take Figure 2 as an example, image source *J* (1) <sup>1</sup> would affect the heat flux and temperature on boundary *∂B*(2), and thus, a second-order image source *J* (2) <sup>2</sup> (the subscript denotes the order of the image source and the superscript represents the boundary by which the image source is mirrored) needs to be added. It is obvious that for the parallel boundaries shown in Figure 2, an infinite number of image sources will finally be needed. Therefore, to reduce the computational cost, only a finite number of image sources are added and the BCs are finally accounted for by the complementary field *T*ˆ. The *T*˘ is expressed as

$$\mathcal{T} = \sum\_{J=1}^{M} \mathcal{T}^{(J)},$$

where *M* is the total number of image sources. The analytical expression of *T*˘(*J*) is the same as *T*˜(*I*), which is detailed discussed in Section 2.2.

**Figure 2.** The strategy of adding image sources. The original heat source is represented by the red source and the image source is represented by the blue source. Image source *J* (1) <sup>1</sup> is added by mirroring the original source with respect to boundary *∂B*(1), and image source *J* (2) <sup>1</sup> is added by mirroring the original source with respect to boundary *∂B*(2). Image source *J* (2) <sup>2</sup> is added by mirroring image source *J* (1) <sup>1</sup> with respect to boundary *<sup>∂</sup>B*(2). The subscript of *<sup>J</sup>* (2) <sup>2</sup> denotes the order of the image source, and the superscript represents the boundary by which the image source is mirrored. For the parallel boundaries *∂B*(1) and *∂B*(2), an infinite number of image sources need to be added.

The complementary temperature field *T*ˆ is obtained by solving

$$
\frac{\partial \hat{T}}{\partial t} = \kappa \nabla^2 \hat{T} \tag{9}
$$

with the BCs

$$\frac{\partial \vec{T}}{\partial \mathbf{x}\_{i}} \mathbf{n}\_{i} = -\frac{\partial \vec{T}}{\partial \mathbf{x}\_{i}} \mathbf{n}\_{i} - \frac{\partial \vec{T}}{\partial \mathbf{x}\_{i}} \mathbf{n}\_{i\prime} \quad \text{on} \quad \partial V\_{\text{top}} \quad \text{and} \quad \partial V\_{\text{lat}\prime} \tag{10}$$

$$
\dot{T} = -\mathcal{T} - \mathcal{T} + T\_{\text{c}} \quad \text{on} \quad \partial V\_{\text{b}\text{ot}\prime} \tag{11}
$$

and initial condition

$$
\hat{T} = -\hat{T} - \check{T} + T\_{\text{irr}\nu} \quad \text{at} \quad t = 0,\tag{12}
$$

where Equations (10) to (12) are a direct consequence of Equations (2) to (4).

The analytical field *T*˜ serves to partially capture the steep temperature gradient in the vicinity of the heat source, and the image source field *T*˘ is employed to describe the steep temperature gradient when the heat source is close to the boundary. As a result, the temperature gradient in *T*ˆ field is relatively smooth, and thus, a relatively coarse mesh can be applied in solving *T*ˆ.

It should be noted that in the semianalytical approach, as there is the need to add image sources, body *V* needs to be a convex polyhedron (where any two points within the polyhedron can be connected by a line). For arbitrary complex geometries, the semianalytical approach can still be utilized by using the *T*ˆ field only to account for the BCs. Then, without the image source field and since the temperature gradient will become very steep when the heat source is close to the boundary, a fine discretization in the vicinity of the boundary may be necessary for solving the *T*ˆ field to obtain certain accuracy. The computational efficiency and accuracy for using the *T*ˆ field only to account for BCs to model the nonconvex part are detailed discussed in [15]. In the present paper, for the purpose of validating the accuracy of the proposed volumetric heat source model, only convex body *V* is considered so that the accuracy of the model is not very sensitive to the mesh size with the assistance of the image sources.

### *2.2. Volumetric Heat Source*

In the model developed by Yang et al. [14,15], sources *I* = 1 to *N* and the corresponding images sources are modeled as dimensionless point sources. In the present paper, for individual source *I*, it is described as a half-ellipsoidal volumetric source (see Figure 3), and the corresponding analytical solutions for *T*˜ and the gradient of *T*˜ are developed.

**Figure 3.** The discretized heat source is modeled as a half-ellipsoidal volumetric source.

Consider the energy distribution for the heat source *I*, shown in Figure 3, expressed as

$$q(\mathbf{x}\_i) = q\_0 \exp(-2\frac{(\mathbf{x}\_1 - \mathbf{S}\_1^{(I)})^2 + (\mathbf{x}\_2 - \mathbf{S}\_2^{(I)})^2}{r\_l^2} - 2\frac{(\mathbf{x}\_3 - \mathbf{S}\_3^{(I)})^2}{c^2}), \quad i = 1, 2, 3,\tag{13}$$

where *<sup>q</sup>*<sup>0</sup> is the maximum energy density at the position of *xi* <sup>=</sup> *<sup>S</sup>*(*I*) *<sup>i</sup>* . The parameters *rl* and *c* are illustrated in Figure 3. In SLM, the *rl* can be considered as the laser spot radius and the parameter *c* can be seen as the optical penetration depth (OPD) of the laser. Equation (13) satisfies that the energy density reduces to *q*0/e2 at the surface of the ellipsoidal source.

The total energy of the heat source is *Q*, and conservation energy requires that

$$Q = \int\_0^\infty \int\_{-\infty}^\infty \int\_{-\infty}^\infty q(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \mathbf{dx}\_1 \mathbf{dx}\_2 \mathbf{dx}\_3 \,. \tag{14}$$

Evaluating Equation (14) yields

$$q\_0 = \frac{4\sqrt{2}Q}{\pi\sqrt{\pi}cr\_l^2}.\tag{15}$$

*Metals* **2020**, *10*, 1406

If heat source *I* is considered as a single dimensionless point source, the corresponding temperature *T*˜(*I*) *<sup>p</sup>* in a semi-infinite space is given by [24]

$$\mathcal{T}\_p^{(I)}(\mathbf{x}\_i^{(P)}, t) = 2Q\mathbf{G}(\mathbf{x}\_i^{(P)}, t),\tag{16}$$

where

$$G(\mathbf{x}\_i^{(P)}, t) = \frac{1}{\rho c\_p (\pi a \eta)^{3/2}} \exp(-\frac{(\mathbf{x}\_1 - \mathbf{x}\_1^{(P)})^2 + (\mathbf{x}\_2 - \mathbf{x}\_2^{(P)})^2 + (\mathbf{x}\_3 - \mathbf{x}\_3^{(P)})^2}{\eta}).\tag{17}$$

The *η* = 4*α*(*t* − *t* (*I*) <sup>0</sup> ) and *x* (*P*) *<sup>i</sup>* is the coordinates for the point of interest. Source *I* is activated at *t* = *t* (*I*) <sup>0</sup> . Function *G*(*x* (*P*) *<sup>i</sup>* , *t*) is the Green function, which is a fundamental solution of Equation (1). Due to the boundary effect caused by the boundary surface of the semi-infinite space, there is a factor 2 in Equation (16).

Due to the linearity of Equation (1), the temperature for point *x* (*P*) *<sup>i</sup>* at moment *t* induced by a heat source with the energy density of *q*(*xi*) is equal to *q*(*xi*)*G*(*x* (*P*) *<sup>i</sup>* , *t*). In addition, we need to consider the boundary effect caused by the boundary surface of the semi-infinite space. This can be solved by adding an image source with respect to the boundary surface of the semi-infinite space. Consequently, the temperature *T*˜(*I*) caused by a volumetric heat source shown in Figure 3 is expressed as

$$\tilde{T}^{(I)} = \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} q(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) G(\mathbf{x}\_1^{(P)}, \mathbf{x}\_2^{(P)}, \mathbf{x}\_3^{(P)}, t) d\mathbf{x}\_1 d\mathbf{x}\_2 d\mathbf{x}\_3. \tag{18}$$

Finally, the temperature *T*˜(*I*) is given by

$$\mathcal{T}^{(I)} = B \frac{\exp\left(-\left(I\_{xy} + I\_z\right)\right)}{D},\tag{19}$$

where

$$\begin{aligned} B &= \frac{2\sqrt{2}Q}{\rho c\_p r^{3/2}}, \\ I\_{xy} &= -\frac{2[(\mathbf{x}\_1^{(P)} - \mathbf{S}\_1^{(I)})^2 + (\mathbf{x}\_2^{(P)} - \mathbf{S}\_2^{(I)})^2]}{2\eta + r\_I^2}, \\ I\_z &= -\frac{2(\mathbf{x}\_3^{(P)} - \mathbf{S}\_3^{(I)})^2}{2\eta + c^2}, \\ D &= (2\eta + r\_I^2)\sqrt{2\eta + c^2}. \end{aligned} \tag{20}$$

It can be seen that if *rl* and *c* both become 0, source *I* would become a dimensionless point source, and Equation (19) would be equivalent to the temperature caused by a point source as shown in Equation (16). For any given heat source, the total energy *Q* = *PA*Δ*t*, where *P* is the power and *A* is a nondimensional coefficient which implicitly accounts for the energy losses during the SLM process, such as due to laser absorptivity, convection, radiation, and phase transitions. The expression for any given image source *J* is the same as Equation (19) but with different source positions *S*(*J*) *<sup>i</sup>* . Therefore, with Equation (19), the analytical field *T*˜ and image source field *T*˘ can be finally obtained using Equations (6) and (8), respectively.

The gradients of *T*˜ and *T*˘ used for solving the *T*ˆ field are given by

$$\begin{split} \frac{\partial \mathcal{T}}{\partial \mathbf{x}\_{i}} &= -\frac{4\mathbf{x}\_{i}^{(P)}}{2\eta + \lambda^{2}} \mathcal{T}, \\ \frac{\partial \mathcal{T}}{\partial \mathbf{x}\_{i}} &= -\frac{4\mathbf{x}\_{i}^{(P)}}{2\eta + \lambda^{2}} \mathcal{T}, \end{split} \tag{21}$$

where *λ* = *rl* for *i* = 1, 2 and *λ* = *c* for *i* = 3.

### **3. Results and Discussions**

Three numerical examples are investigated in the present study. For the first example, a single laser scan is applied on a very large baseplate to mimic the scanning of the first layer. The baseplate can be assumed to be a semi-infinite space and thus, the total temperature field *T* simply becomes *T*˜ in the absence of any BCs. The predicted melt-pool width and depth by the proposed volumetric heat source model are compared with the experiments reported in [25]. The material used is stainless steel (SS) 316L. For the second example, a single laser scan is applied on a baseplate with the cubic profile and thus, the BCs need to be considered. The material Ti6Al4V is considered in this example to demonstrate that the proposed model can be applied to different materials. The accuracy of the proposed model is further validated by comparing the predicted melt-pool width and depth with the experiments [8]. A third example is presented to demonstrate the ability of the proposed model to simulate the multiple laser scanning process with various heat source moving paths in the AM process.

The material properties of SS316L and Ti6Al4V are tabulated in Table 1. As suggested by Yang et al. [14], the values quoted in Table 1 are representative for a temperature which is slightly higher than the melting point (*Tm*) of SS 316L and Ti6Al4V, respectively. The thermal properties for SS316L and Ti6Al4V in Table 1 correspond to the temperatures of 1727 oC [26] and 2227 oC [27], respectively. The temperature fields *T*˜, *T*˘, and *T*ˆ are all calculated using an in-house Matlab code. The *T*˜ and *T*˘ fields are directly calculated using Equations (6) and (8), and the *T*ˆ is solved with an explicit finite difference scheme, centered in space and forward in time.

**Table 1.** Material properties.


### *3.1. A Single Laser Scan on a Semi-Infinite Space*

A single laser scan is applied on a semi-infinite space as schematically illustrated in Figure 3 to mimic scanning the first SS 316L powder layer on a very large baseplate [25]. The initial temperature is set to be 25 oC. The laser spot radius *rl* is 27 μm and the OPD *c* is set as 60 μm, which is close to the diameter of the powder (54 μm).

A convergence study with respect to the source density *ρ<sup>s</sup>* is first investigated. The maximum temperature experienced by a point with a distance of 20 μm to the laser scanning line with various source density *ρ<sup>s</sup>* is calculated. This distance is close to the laser spot radius. The source density *<sup>ρ</sup><sup>s</sup>* <sup>=</sup> <sup>100</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> Wms−<sup>1</sup> is taken as the reference. As shown in Figure 4, the error <sup>|</sup>*T*˜(*ρs*) <sup>−</sup> *<sup>T</sup>*˜ ref|/*T*˜ ref for two different *P*/*v* is plotted as a function the source density *ρs*. It can be seen that the error quickly reduces to a value close to 0 and converges at *<sup>ρ</sup><sup>s</sup>* = <sup>10</sup> × 104 Wms<sup>−</sup>1. Hence, for the following calculations, the source density *<sup>ρ</sup><sup>s</sup>* = <sup>10</sup> × <sup>10</sup><sup>4</sup> Wms−<sup>1</sup> is used.

**Figure 4.** The temperature convergence of *T*˜ as a function of the source density *ρs*.

The predicted melt-pools by the volumetric heat source model within the *x*<sup>1</sup> − *x*<sup>2</sup> plane and *x*<sup>2</sup> − *x*<sup>3</sup> plane are shown in Figure 5a,b, respectively. The power is 150 W and the speed is 0.8 m/s. The melt-pool is shown as the white area in Figure 5 and determined as the material points heated above the melting point. A regular hexahedral grid with the size of 10 μm is employed to plot Figure 5. The melt-pool width and depth along the scanning line at *P* = 150 W and *v* = 0.8 m/s are shown in Figure 6. By examining the temperature of the material points in the vicinity of the scanning line (along the width and depth directions, respectively), the boundary of the melt-pool is determined to be the first material point with the temperature below the melting point. The resolution for determining the melt-pool width and depth in this example is 1 μm. It can be seen that the melt-pool width and depth are both small at the beginning of the laser scan, but quickly reach a higher steady-state value. The melt-pool width and depth then decrease towards the completion of the laser scan. Hence, the calculated width and depth of the melt-pool shown in the following are determined by the steady-state value along the laser scanning line.

**Figure 5.** The temperature profile within the (**a**) *x*<sup>1</sup> − *x*<sup>2</sup> plane and (**b**) *x*<sup>2</sup> − *x*<sup>3</sup> plane for the power of 150 W and the speed of 0.8 m/s. The melt-pool is shown as the white area.

**Figure 6.** The melt-pool width and depth along the laser scanning line at the power of 150 W and the speed is 0.8 m/s.

The melt-pool width and depth by the proposed volumetric heat source model for various laser powers and speeds, and the corresponding experimental measurements are reported in [25] are plotted in Figure 7a,b. The coefficient *A* is set to be 0.43 to result in the best agreement for the volumetric heat source model. The simulation results by the point heat source model are also plotted in Figure 7.

The uncertainty of the experimental measurements is 5 μm. As explained in [25], the experimental uncertainty is mainly because that the plasma/metal vapor plume during the SLM process can change the laser absorptivity, which causes fluctuations of the effective power. Incorporating a variable coefficient *A* in the proposed model to accurately account for the variance of the effective power during the SLM would be quite difficult and thus, a more common approach [25] to take the coefficient *A* as a constant is employed. It can be observed from Figure 7 that the melt-pool width and depth increase with the increasing *P*/*v*. The value of *P*/*v* for the power of 200 W and speed of 1.2 m/s is the same as that for the power of 300 W and speed of 1.8 m/s, and thus, close values of melt-pool width and depth can be observed in Figure 7 both for the experimental and simulation results.

**Figure 7.** The predicted melt-pool (**a**) width and (**b**) depth and the experimental measurements under various laser powers and scanning speeds.

In Figure 7, it can also be seen that the predicted melt-pool width and depth by the volumetric heat source model agree well with the experimental measurements for all the four sets of laser power and speed. In comparison, the predicted melt-pool width by the point heat source model is overestimated (see Figure 7a) and the depth is underestimated (Figure 7b). The predicted melt-pool depth by the point heat source model is actually half of the predicted width as the heat energy diffuses isotropically in the space domain for the point heat source. If we tune the coefficient *A* to make the width calculated by the point heat source model agree well with the experiments, the calculated depth will also reduce accordingly. Therefore, the difference between the half width and depth of the melt-pool cannot be captured by the point heat source model. In the volumetric heat source model, by adjusting the parameter *c* to characterize the OPD of the laser, the melt-pool with shallow or deep depth can be conveniently captured.

The error of the predicted results in Figure 7 with respect to the average experimental measurements is evaluated by

$$\varepsilon = \frac{|\mathcal{Y} - \mathcal{Y}\_{\text{exp}}|}{\mathcal{Y}\_{\text{exp}}},\tag{22}$$

where *Y* represents the predicted melt-pool width or depth, and *Y*exp represents the averaged experimental melt-pool width or depth. Figure 8 shows the errors of the predicted melt-pool width and depth for both the point and volumetric heat source models. The rectangular bars represent the average error *ε*¯ for the four sets of laser power and speed, and the error bar corresponds to the maximum and minimum error, respectively, in the four sets of laser power and speed. It can be seen from Figure 8 that the maximum errors of the predicted melt-pool width and depth for the point heat source model are close to 20%, while the corresponding maximum errors for the volumetric heat source model are less than 10%. The average errors of the predicted melt-pool width and depth for the point heat source model are around 15%, and the corresponding average errors for the volumetric heat source model are around 5%.

**Figure 8.** The errors *ε* = |*Y* − *Y*exp|/*Y*exp of the predicted width and depth with respect to the experiments.

### *3.2. A Single Laser Scan on a Finite Space*

In this example, a single scan experiment reported in literature [8] is simulated using the proposed model. As shown in Figure 9, the baseplate has a dimension of 4 mm × 2 mm × 0.5 mm and the material is Ti6Al4V. The bottom surface of the baseplate is set to be at 20 oC during the whole process. The laser spot radius is 26 μm and the OPD *c* is set to be 40 μm, which is also close to the powder diameter (35 μm). The coefficient *A* is set to be 0.77 as suggested in [8]. The laser speed is 0.2 m/s and various laser powers are investigated. A total number of 64 cubic finite difference cells is used to

discretize the baseplate, resulting in a cell size of 1 mm <sup>×</sup> 0.5 mm <sup>×</sup> 0.125 mm. The *<sup>T</sup>*˜ and *<sup>T</sup>*˘ fields can be readily obtained for any point of interest *x* (*P*) *<sup>i</sup>* within the body analytically, while the *<sup>T</sup>*<sup>ˆ</sup> field can only be calculated at the nodes of the finite difference cells. Therefore, a linear interpolation of known *T*ˆ values is utilized to estimate the value of *T*ˆ(*x* (*P*) *<sup>i</sup>* ) on any given material point, which is given by

$$
\hat{T}(\mathbf{x}\_i^{(P)}) = \sum\_{q=1}^{N\_q} \Phi^{(q)} \hat{T}(\mathbf{x}\_i^{(q)}),\tag{23}
$$

where *q* = 1, 2, ... , *N* is the index for the grid point with *T*ˆ(*x* (*q*) *<sup>i</sup>* ) and the total number of grid points of the finite difference cell is denoted as *Nq*. For an 8-node hexahedral finite difference cell (*Nq* = 8) used in this example, the function Φ(*q*) is expressed as

$$\Phi^{(q)} = \frac{1}{8} (\mathbf{1} + \mathfrak{J}\_i \mathfrak{k}^{(q)}),\tag{24}$$

where *ξ<sup>i</sup>* is the normalized position of *x* (*P*) *<sup>i</sup>* in a right-handed local coordinate system having an origin located at the center of the cell and *ξ*(*q*) is the normalized position of *x* (*q*) *<sup>i</sup>* given in the same coordinate system. The resolution for determining the melt-pool width and depth in this example is 1 μm.

**Figure 9.** A single laser scan is applied on a finite space.

The predicted melt-pool width and depth obtained by the point and volumetric heat source models as well as the experimental measurements are shown in Figure 10a,b, respectively. The experimental results were measured for different places along the scanning line by optical microscopy based on the solidified microstructure [8]. As explained in [8], because the boundaries of the melt-pool were difficult to characterize, some large error bars can be observed in Figure 10. It can be seen from Figure 10 that the trend of the predicted melt-pool width and depth obtained by the point and volumetric heat source models agree well with that of the experimentally obtained results. The melt-pool width and depth increase with the increasing *P*/*v*. The melt-pool widths predicted by the point and volumetric heat source models both agree reasonably well with the experiments. However, in Figure 10b, it can be seen that the melt-pool depth is underestimated by the point heat source model at the powers of 20 W and 40 W. In comparison, the melt-pool depths predicted by the volumetric source model have good agreements with the experimentally measured values at all the four sets of laser power. This is because the different heat energy distributions within the laser scanning plane (i.e., the *x*<sup>1</sup> − *x*<sup>2</sup> plane) and along the depth direction can be well captured by the volumetric heat source model. In comparison, for the point heat source model, both the heat energy diffusion within the laser scanning plane and the the phenomenon of penetration along the depth direction cannot be accurately described. For the powers of 60 W and 80 W, the melt-pool depths happen to be approximately half of the melt-pool width, and thus, the predicted depths obtained by the point heat source model agree well with the experimental results.

A nonlinear thermal analysis using the finite element method was also conducted in [8], in which the nonlinear thermal properties of Ti6Al4V were employed. The temperature distributions along the width direction (see line 1 in Figure 9) at the powers of 40 W and 80 W at a certain time calculated by the proposed volumetric heat source model were compared with the results obtained by the nonlinear model in [8], as shown in Figure 11. Good agreements can be observed between the nonlinear model and the proposed approach. The disparities of the temperature obtained by the two models shown in Figure 11 are characterized by *κ* = |*T*vol − *T*non|/*T*non, where *T*vol is the temperature obtained by the proposed volumetric heat source model and *T*non is the temperature obtained by the nonlinear model proposed in [8]. It is found the *κ* is less 10%. The small disparities are expected, as a nonlinear set of thermal properties which is temperature dependent is used in the nonlinear model while constant thermal properties are employed in the proposed model.

**Figure 10.** The predicted melt-pool (**a**) width and (**b**) depth with the experimental measurements.

**Figure 11.** The temperature along the width direction (see line 1 in Figure 9) obtained by the proposed volumetric heat source model and the nonlinear model reported in [8].

### *3.3. Building a New Layer with Multiple Laser Scans*

Finally, a numerical example for scanning a new layer with multiple scans is conducted. Two scanning strategies are considered in the proposed model. The material of Ti6Al4V is employed

in this example. As shown in Figure 12a, a cube with the dimension of *l* = 2 mm is assumed already built, and two scanning strategies (unidirectional, see Figure 12b; and alternating, see Figure 12c) are considered to build the new layer. A total number of 64 cubic finite difference cells is used to discretize the body. The temperature of the bottom surface of the cube is set to be 20 oC during the whole process. An arbitrary set of process parameters is chosen to investigate the temperature histories of point G1 and G2 (see Figure 12a). The laser spot radius is set to be 50 μm and the OPD is 60 μm. The power of the laser is set to be 40 W with a speed of 0.8 m/s. The hatching distance between adjacent tracks is 90 μm and a total number of 21 tracks are applied on the top surface *∂V*top. The thermal response is calculated by the volumetric heat source model.

Figure 13a shows the temperature histories of point G1 under the unidirectional and alternating scanning strategies. The temperature evolutions under the two scanning strategies are well captured. The first peak temperature for the alternating scanning strategy appears earlier than that for the unidirectional scanning strategy. This is because compared with the unidirectional scanning strategy, the first track in the alternating scanning strategy arrives earlier at point G1. The second laser track is the same in both alternating and unidirectional scanning strategies, and thus, the second peak temperatures of point G1 for both scanning strategies occur at the same moment. It also demonstrates that the temperature history of point G1 is sensitive to the scanning strategy. The temperature field *T*˜ without considering the BCs under the unidirectional scanning strategy is also plotted in Figure 13a. The peak temperature of *T*˜ is only 1439 oC, while the peak temperature of *T* is 4493 oC. It can be observed that after including the BCs, the total temperature *T* is reasonably much higher than the *T*˜, which indicates it is essential to consider the BCs in the thermal modeling of the AM process.

**Figure 12.** (**a**) A total of 21 laser tracks are applied on the top surface *∂V*top of the cube, which is discretized by 64 finite difference cells. The temperature histories of point G1 (1.5, 0, 0) mm and G2 (1.5, −0.5, 0) mm are investigated. Two scanning strategies, (**b**) unidirectional and (**c**) alternating, are employed.

The temperature histories of point G2 under the two different scanning strategies are compared in Figure 13b. It can be seen that the peak temperature of point G2 is reasonably lower than that of point G1 because point G1 is at the boundary. The peak temperature in Figure 13b is around 1000 oC, which is lower than the melting point 1650 oC. This indicates that the set of process parameters employed may cause the lack of fusion. It is also important to note that scanning strategies can modify the temperature histories which, in turn, have an impact on the resulting residual stresses.

The calculation time for this example by the volumetric heat source model is around 110 s. The calculation associated with *T*ˆ field is only 0.45 s. The high efficiency is because the steep temperature gradient is captured by the *T*˜ and *T*˘ fields, which can be calculated analytically, and thus, coarse mesh can be applied in solving *T*ˆ. With this high computational efficiency, it would be very efficient to investigate the thermal histories of an AM part and optimize the process parameters. The calculation is performed using a single-core Intel i7 6600U quad-core processor with a clock speed of 2.60 GHz and 8 GB RAM.

**Figure 13.** Temperature histories of (**a**) point G1 and (**b**) point G2.

### **4. Limitations and Future Extensions**

In the proposed model, a half-ellipsoid heat source is employed to describe the interaction between the laser and the materials in SLM. Since the laser spot radius is usually small in SLM, it is appropriate to assume a circle heat source profile within the laser scanning plane (*x*<sup>1</sup> − *x*<sup>2</sup> plane). When the size of heat source becomes large—such as in wire arc additive manufacturing, for which the characteristic size of the heat source is usually in mm scale—a more complex heat source profile, such as the Goldak double ellipsoidal heat source [16], may need to be applied. Therefore, suitable heat source profiles for different metal AM processes with the corresponding analytical solutions of *T*˜ and Grad (*T*˜) need to be further developed. Moreover, in the half-ellipsoid heat source, the OPD of the laser is characterized by the constant parameter *c*, which is independent of the laser power and speed. However, the OPD may vary for different processing parameters, and this cannot be accounted for by the present model. Hence, a thermal model which is capable of characterizing the variance of OPD as a function of the process parameters is also of interest to investigate.

### **5. Conclusions**

A computationally efficient volumetric heat source model based on the semianalytical thermal approach is proposed to simulate the thermal response of the produced part in the AM process. It has shown that the asymmetry of the melt pool in width and depth directions can be well captured by the proposed volumetric heat source model. Two materials of SS 316L and Ti6Al4V are employed to validate the accuracy of the proposed model. As the steep temperature gradients are mainly captured by the analytical field *T*˜ and the image source field *T*˘, which can be obtained analytically, high computational efficiency can be achieved in solving *T*ˆ field. Compared with the point heat source model, the volumetric heat source model has higher accuracy and in the meantime, the computational cost of the proposed volumetric heat source model is not impaired because of the concise expressions of *T*˜ and Grad(*T*˜). Moreover, the computational efficiency of the proposed model is independent of the laser spot radius and OPD, indicating that various laser profiles can be easily modeled. The high computational efficiency of the proposed model suggests it is promising to be incorporated in an optimization process (optimizing the topology of the built part or process parameters) to reduce the residual stresses and distortions.

**Author Contributions:** Conceptualization, Y.Y.; methodology, Y.Y.; validation, Y.Y.; resources, X.Z.; writing–review and editing, Y.Y., X.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Science and Technology on Plasma Dynamics Laboratory, Air Force Engineering University, Xi'an, China (Project No. 614220206021808).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
