Zhibo Luo
email: zhibo.luo@mail.mcgill.ca

Yaoyao Zhao
email: yaoyao.zhao@mcgill.ca

Numerical simulation of part-level temperature fields during selective laser melting of stainless steel 316L

Keywords: finite element model, selective laser melting (SLM), adaptive mesh, thermal analysis

Numerical transient thermal analysis of a part during the selective laser melting process is a great challenge due to huge computational cost. An efficient simulation scheme is proposed to estimate the temperature distribution and history in part-level. The building domain is coarsely meshed at the beginning. Then deposited elements are activated layer by layer with desired material properties. The Gaussian distributed heat flux is applied to the top surface and travels according to the predefined scanning path. With the motion of input heat flux, the mesh is continuously refined and coarsened as the laser beam moves in and out. To improve the solution accuracy, temperature-dependent thermal properties in powder and solid state are incorporated in the finite element system. To reduce the computation cost, Gaussian line heat source with a large time step length is used to replace the traditional moving Gaussian point heat source. To achieve a compromise between the computational efficiency and the solution accuracy, the hybrid of Gaussian line heat source and Gaussian point heat source is adopted as heat input and results are compared with that of Gaussian point heat source. The proposed simulation scheme is validated by comparing the simulated geometry of the melt pool with experimental results. Simulation results show that the proposed simulation scheme has the ability to efficiently and accurately predict the temperature field of a part during the selective laser melting process and can be easily extended to other powder bed fusion process.

Introduction

Additive manufacturing (AM) is a newly arisen manufacturing process compared to the traditional subtractive manufacturing processes such as milling, drilling, or machining. It is easy to build a part with complicated internal or external geometric features with acceptable cost through the layer-by-layer method. Dalgarno and Stewart [START_REF] Dalgarno | Manufacture of production injection mould tooling incorporating conformal cooling channels via indirect selective laser sintering[END_REF] conclude that this layer-based manufacturing of complex internal conformal cooling channels is economically competitive compared to conventional methods. Selective laser melting (SLM) process is one of the most commonly used powder bed fusion AM processes. The SLM fabricated part is formed by bonding the powder material in the manner of layer-by-layer through the laser beam. Since the SLM process almost imposes no design constraints to designers and manufacturers, it has been extensively applied in injection mold, aerospace, automobile, and medical industry. However, high-temperature gradients during the thermal cycles of the SLM process tend to induce high thermal stresses which result in the distortion of deposited layers. Defects induced by thermal stress such as cracks, warping and distortion, fatigue failure, and delamination from the build plate [START_REF] Mazur | Deformation and failure behaviour of Ti-6Al-4V lattice structures manufactured by selective laser melting (SLM)[END_REF][START_REF] Wang | Research on the fabricating quality optimization of the overhanging surface in SLM process[END_REF][START_REF] Cherry | Investigation into the effect of process parameters on microstructural and physical properties of 316L stainless steel parts by selective laser melting[END_REF][START_REF] He | Rigid-flexible coupling virtual prototyping-based approach to the failure mode, effects, and criticality analysis[END_REF], are the main reasons for the dimension error and the failure of the printing process. These issues significantly preclude the wide adoption of AM in large-scale industrial production [START_REF] Mani | Measurement Science Needs for Real-time Control of Additive Manufacturing Powder Bed Fusion Processes[END_REF]. Many research attempts have been conducted to experimentally or numerically investigate these defects during the SLM process. Edwards and Ramulu [START_REF] Edwards | Effect of build direction on the fracture toughness and fatigue crack growth in selective laser melted Ti-6Al-4 V[END_REF] experimentally studied the fatigue crack growth characteristics of Ti-6Al-4V in SLM process and concluded that residual stress is the main factor in the growth of crack as the formation of microstructure is directional to the build orientation. Liu, et al. [START_REF] Liu | A study on the residual stress during selective laser melting (SLM) of metallic powder[END_REF] experimentally investigated the origin of residual stress in terms of temperature gradients. The effects of temperature gradient on residual stress and microstructure are experimentally explained. The process parameters such as energy density and scanning track length are studied to optimize the support structure. Ma and Bin [START_REF] Ma | Temperature and stress analysis and simulation in fractal scanningbased laser sintering[END_REF] numerically analyzed the influence of scanning pattern on the evolution of temperature, distortion, and thermal stress and concluded that 'S' scanning pattern is easier to be warp compared with the fractal scanning pattern due to the symmetric temperature distribution. Foroozmehr and Kovacevic [START_REF] Foroozmehr | Effect of path planning on the laser powder deposition process: thermal and structural evaluation[END_REF] also numerically investigated the effect of material deposition pattern on the residual stress and concluded that short scanning track has less temperature variation thus results in reduced residual stress in the deposited domain. Usually, extra support structures are needed for building overhang structures and prevent these features from warpage. However, the additional support structure is not easy to remove in the post-process and requires extra printing time and material cost. Vora, et al. [START_REF] Vora | AlSi12 in-situ alloy formation and residual stress reduction using anchorless selective laser melting[END_REF] proposed a novel method to relieve the stress during the building process without support structures by using two types of powder. This method does not guarantee the elimination of support structures due to complex thermal cycles among the melt pool. As the thermal history of the building domain is an accumulation process from the melt pool, Bauereiß, et al. [START_REF] Bauereiß | Defect generation and propagation mechanism during additive manufacturing by selective beam melting[END_REF] developed a mesoscopic numerical model to investigate the melting process in the melt pool. This model considered hydrodynamics such as capillary and wetting angles to predict the mechanism of fault generation. Cheng and Chou [START_REF] Cheng | Thermal Stresses Associated with Part Overhang Geometry in Electron Beam Additive Manufacturing: Process Parameter Effects[END_REF] and Zaeh and Branner [START_REF] Zaeh | Investigations on residual stresses and deformations in selective laser melting[END_REF] demonstrated the predominant reason for macro-level residual stress in parts fabricated via SLM process is the significant thermal gradient during rapid thermal cycles in different regions.

To overcome the high-temperature gradient during the SLM process, Simson, et al. [START_REF] Simson | Residual stress measurements on AISI 316L samples manufactured by selective laser melting[END_REF] and Hussein, et al. [START_REF] Hussein | Advanced lattice support structures for metal additive manufacturing[END_REF] conducted a large number of experiments to minimize the temperature gradient by optimizing the process parameters such as the laser power, scanning speed, scanning path, etc. This experimental method is not suitable for arbitrary part geometry due to high cost. In addition, it cannot obtain the in-situ temperature fields and thermal stress fields inside of the build domain through in-process measurement or monitoring. Finite element (FE) method has been widely applied to model the material thermal and mechanical behaviors during and after the building process. FE analysis of temperature distribution, thermal stress, and deformation in the SLM process is similar to that of laser multi-pass welding process [START_REF] Yilbas | Laser welding of AISI 316 steel: microstructural and stress analysis[END_REF][START_REF] Denlinger | Residual stress and distortion modeling of electron beam direct manufacturing Ti-6Al-4V[END_REF], which have been developed and modeled for four decades. Therefore, it can be extended to simulate the SLM process by increasing the welding tracks. FE method is suitable for simulating the temperature and thermal stress evolution of a single welding track or several tracks, but it is not suitable for simulating the whole history of material thermal-mechanical behavior of a large and complex part in the SLM process due to the huge computation cost. Safronov, et al. [START_REF] Safronov | Distortions and Residual Stresses at Layer-by-Layer Additive Manufacturing by Fusion[END_REF] concluded that the principal difficulty lays on the small size of the melt pool and large size of the building part compared with the welding process in which only have limited number and relatively large size of welding tracks. For example, Cheng, et al. [START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF] developed an FE model using Abaqus to simulate the temperature field and residual stress field of 8 mm × 8 mm × 1 mm part using 0.4 mm laser beam diameter. This model took more than 90 hours for the weak coupled thermal-mechanical analysis. Denlinger, et al. [START_REF] Denlinger | Thermomechanical modeling of additive manufacturing large parts[END_REF] developed a model based on Pan Computing code [START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF], which can refine and coarse the mesh with the motion of heat source, to simulate the thermal and corresponding mechanical behaviors of a large part that scales to meters. Recently, Jayanath and Achuthan [START_REF] Jayanath | A Computationally Efficient Finite Element Framework to Simulate Additive Manufacturing Processes[END_REF] proposed an efficient numerical model with adaptive mesh to estimate the temperature and thermal stress for multiple layers. The computational time of this FE framework is reduced by about 50% compared to that of the conventional FE framework. Luo and Zhao [START_REF] Luo | A survey of finite element analysis of temperature and thermal stress fields in powder bed fusion Additive Manufacturing[END_REF] reviewed existing FE methods in the modeling of the SLM process and summarized following drawbacks: First, the existing multiscale modeling approach does not incorporate the influence of crucial process parameters such as scanning path, scanning speed, layer thickness, and hatching space, etc., all of which have a very important influence on the temperature distribution [START_REF] Ma | Temperature and stress analysis and simulation in fractal scanningbased laser sintering[END_REF]; Second, the spreading of powder material on top of deposited solid material is not easy to implement using available FE models due to quite different solid and powder material properties; Third, the moving Gaussian point heat source requires such small timestep that leads to unacceptable simulation time when modeling the SLM process in part level; Fourth, traditional static mesh scheme cannot trace the motion of laser beam and thus is not appropriate for a part level simulation.

Numerical simulation of the SLM process can be classified into three levels based on the dimension of simulated domain: micro-level simulation of heat, mass, and fluid flow among the melt pool using CFD model in the range of several hundred microns; meso-level simulation of temperature and displacement fields of several scanning tracks or layers; and macro-level simulation of the whole SLM process in part scale. The numerical simulation of micro-level and meso-level have been widely studied and revealed many physical phenomena among the melt pool such as balling effects induced by the surface tension between the melted and unmelted powder [START_REF] Körner | Fundamental consolidation mechanisms during selective beam melting of powders[END_REF], influence of Plateau-Rayleigh capillary instability on fragmentation of the scanning track [START_REF] Gusarov | Heat transfer modelling and stability analysis of selective laser melting[END_REF], and material transportation between several layers caused by material accumulation in the melt pool [START_REF] Vandenbroucke | Selective laser melting of biocompatible metals for rapid manufacturing of medical parts[END_REF], etc. However, the FE modeling of SLM process in macro-level has much more practical use and has many degrees of improvement. Though there already exist other methods to simulate the SLM building process such as inherent strain method [START_REF] Keller | New method for fast predictions of residual stress and distortion of AM parts[END_REF] and multiscale method [START_REF] Markl | Multi-Scale Modeling of Powder-Bed-Based Additive Manufacturing[END_REF], the inherent strain method only conducts a static mechanical analysis without thermal load for each layer. And the multiscale method conducts a series of simulations in different levels and then reduces the scale of the problem by applying equivalent heat flux for one or several layers. These attempts cannot incorporate the key parameters that significantly influence the material thermal and mechanical behaviors thus cannot accurately predict the quality of the final product.

Very limited studies have been conducted to investigate the simulation scalability at macro level. This research aims to efficiently and accurately simulate the temperature field during the SLM process in part level. Compared with other numerical models that attempt to reduce the scale of the problem through layer model or multiscale model, this research detailly simulates the SLM process track-by-track and layer-by-layer with an efficient simulation scheme. This new simulation scheme has more freedom to estimate the effect of key process parameters on the temperature field of a part. Two ways are proposed to improve the efficiency and accuracy: Kelly error based adaptive mesh scheme is applied to reduce the cell number; line heat source model is developed to increase the time step length. The adaptive mesh scheme can generate the best FE mesh for SLM process due to highly non-uniform distributed input heat flux in space and time. In addition, the adaptive mesh scheme in this research is based on the Kelly error estimation of temperature gradients among each cell, which can minimize the error induced by mesh density. A new hybrid heat source which combines Gaussian distributed line heat flux and point heat flux is developed to improve the simulation efficiency by increasing time step size while keeping the solution accuracy around the contours in each layer. As thermal stresses are mainly caused by temperature gradients, then the obtained temperature fields can be used to identify the possible thermal concentration region which may lead to high thermal stress and deformation in the real fabrication process. Also, the estimated temperature history of each point can be used to feed the evolution of microstructure. Next, the transient thermal analysis for the SLM process is elucidated comprehensively in section 2. Then, the systematic design of the simulation scheme based on an open source FE library is discussed in section 3. In section 4, the efficiency and accuracy of the proposed simulation scheme are analyzed and compared with traditional simulation scheme. The predicted size and shape of the melt pool using the proposed new heat source models are compared with experimental results. Finally, the effect of different heat sources on temperature distribution and history of a part with overhang features will be studied and analyzed in section 5.

Finite element model for thermal analysis

In the SLM process, if a reference coordinate is fixed to the laser spot, this reference framework denotes a Eulerian frame. In the Eulerian frame, if the laser spot is sufficiently far from the border, the evolution of temperature distribution in powder bed does not change and thus can be modeled as a steady-state problem. Thermal analysis using the Eulerian reference frame is mainly used to calculate the temperature distribution and geometry of the melt pool [START_REF] Mishra | A heat-transfer and fluid-flow-based model to obtain a specific weld geometry using various combinations of welding variables[END_REF][START_REF] Zacharia | Weld pool development during GTA and laser beam welding of Type 304 stainless steel; Part I-theoretical analysis[END_REF][START_REF] Kim | Modeling of temperature field and solidified surface profile during gas-metal arc fillet welding[END_REF], or simulate other physical phenomena such as buoyancy, powder melt pool interaction, surface tension, and magneto-hydro-dynamic effects [START_REF] Zacharia | Three-dimensional transient model for arc welding process[END_REF][START_REF] Han | Modeling of laser deposition and repair process[END_REF]. However, it cannot capture the transient thermal response of the whole build domain. This paper will use the Lagrangian reference frame, where a reference coordinate is fixed to the powder bed, to establish the transient FE formulations for thermal analysis of the SLM process.

Governing equation

For the heat transferred in powder material, the scattered particle material is considered as uniform and continuum isotropic solid material in the SLM process, which is the fundamental assumption in many research [START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF][START_REF] Hussein | Finite element simulation of the temperature and stress fields in single layers built without-support in selective laser melting[END_REF]. For continuum isotropic solid material, the governing equation of heat transfer process can be derived as,

𝜌𝐶 ! "#(𝒙,') "' = ∇ • [𝑘∇𝑇(𝒙, 𝑡)] + 𝑄(𝒙, 𝑡) ∀𝒙 ∈ Ω, 𝑡 ∈ (0, 𝑡 )*+ ) (1)
Where 𝐶 ! is the specific heat capacity, 𝑘 is the thermal conductivity tensor, 𝜌 is the density, 𝑄 is the internal volume heat generation. Note that 𝑘 and 𝐶 ! are temperature-dependent material properties.

In the SLM process, Fischer, et al. [START_REF] Fischer | Sintering of commercially pure titanium powder with a Nd: YAG laser source[END_REF], Ma and Bin [START_REF] Ma | Temperature and stress analysis and simulation in fractal scanningbased laser sintering[END_REF] claimed that the thermal heat is mainly absorbed through laser irradiation, which can only penetrate the surface of the solid material within a few microns. This penetration depth is much smaller compared with the element size in the FE mesh. Therefore, the laser beam in the SLM process is considered as a surface-related heat flux.

where 𝜕𝑇/𝜕𝐧 is the temperature gradient along the normal direction, 𝑇 7 is the prescribed temperature on the boundary 𝛤 , , 𝑞 ./*0 and 𝑞 12+ are the prescribed thermal convection term and irradiation term on the boundary 𝛤 # , 𝑞 3456 is the input heat flux.

F 𝑞 ./*0 = ℎ • [𝑇(𝒙, 𝑡) -𝑇 29: ] ∀𝒙 ∈ 𝛤 ./*0 ⊆ 𝛤 # , 𝑡 ∈ (0, 𝑡 )*+ ) 𝑞 12+ = 𝜖𝜎 • [𝑇 ; (𝒙, 𝑡) -𝑇 29: ; ] ∀𝒙 ∈ 𝛤 12+ ⊆ 𝛤 # , 𝑡 ∈ (0, 𝑡 )*+ ) (3)
where ℎ, 𝜖 is the temperature-dependent heat convection coefficient and surface emissivity, respectively. 𝜎 is the Stefan-Boltzmann constant.

During the printing process, for each layer, the input heat flux is applied to the top surface.

-𝑘 • "#(𝒙,')

"𝐧 = 𝑞 3456 ∀𝒙 ∈ 𝛤 3456 , 𝑡 ∈ (0, 𝑡 )*+ ) (4)
The initial condition with constant temperature distribution can be written as, 𝑇(𝒙, 0) = 𝑇 29: ∀𝒙 ∈ Ω, 𝑡 = 0 (5)

Finite element formulation of the problem

This section describes the space and time discretization for the actual building process. First, the Galerkin discretization method is applied to convert the continuum problem with infinite number of unknowns into a discrete one with finite number of unknowns. This includes two major steps: to obtain the weak formulation of the problem in a Hilbert space 𝑉 which has continuous form; and to approximate the continuous Hilbert space with a subspace of finite dimension 𝑉 < through a linear combination of the basis vectors in 𝑉 < . Second, approximate the continuous space with a subspace of finite dimension and obtain the finite element formulation.

Weak formulation of the problem

For the transient heat transfer problem, the strong form equation should be discretized in both time and space. First, equation ( 1) is discretized in time by applying the 𝜃-scheme time discretization, then the semi-discrete version of the strong form can be obtained,

𝜌𝐶 ! " # ! $# !"# 𝜏 -𝜃∇(𝑘 " • ∇𝑇 " ) -(1 -𝜃)∇(𝑘 "$% • ∇𝑇 "$% ) = 𝜃𝑄 " + (1 -𝜃)𝑄 "$% ( 6
)
Where 𝑛 is the time step number, 𝑇 * is used to approximate 𝑇(𝒙, 𝑡 * ) at time 𝑡 * . 𝐶 ! * , 𝑘 * and 𝑄 * is the specific heat capacity, thermal conductivity, and internal volume heat generation at time 𝑡 * , respectively. 𝜏 = 𝑡 * -𝑡 *>? is the time step size. 𝜃 is the parameter that controls the time discretization method. This research chooses 𝜃 = 0.5, which generates the Crank-Nicolson time discretization and has second order convergence.

Second, the whole domain Ω is discretized in space. After discretization of the problem in time, then by multiplying a test function 𝑣, integrating by parts, and applying boundaries, it gives the following discrete equations.

𝜌𝐶 𝑝 𝑛 (𝑇 * , 𝑣) @ -𝜌𝐶 𝑝 𝑛 (𝑇 *>? , 𝑣) @ + 𝜏𝑘 𝑛 𝜃(∇𝑇 * , ∇𝑣) @ + 𝜏𝑘 𝑛-1 (1 -𝜃)(∇𝑇 *>? , ∇𝑣) @ = 𝜏𝜃(𝑄 * , 𝑣) @ + 𝜏(1 -𝜃)(𝑄 *>? , 𝑣) @ + 𝜏𝜃(𝐧 • 𝑘 𝑛 ∇𝑇 * , 𝑣) "@ + 𝜏(1 -𝜃)(𝐧 • 𝑘 𝑛 ∇𝑇 *>? , 𝑣) "@ (7)
On the boundary 𝛤 , ,

𝐧 • 𝑘∇𝑇 = 0 (8)
On the boundary 𝛤 # (note that

𝛤 # = 𝛤 ./*0 = 𝛤 12+ = 𝛤 3456 ), R 𝐧 • 𝑘∇𝑇 = -ℎ(𝑇) • (𝑇 -𝑇 29: ) 𝑜𝑛 𝛤 ./*0 𝐧 • 𝑘∇𝑇 = -𝜖(𝑇)𝜎 • (𝑇 ; -𝑇 29: ; ) 𝑜𝑛 𝛤 12+ 𝐧 • 𝑘∇𝑇 = 𝑞 3456 (𝑇) 𝑜𝑛 𝛤 3456 (9)

Discretized finite element formulation

The continuous space 𝑉 is approximated with a linear combination of finite number of subspace 𝑉 < . Therefore, the continuous temperature field 𝑇 of a part can be approximated with a linear combination of node temperature in each cell, and

𝑇 = ∑ 𝑇 A ) 𝑣 A ) * ! AB?
. Substitution of this into equation ( 9) and then get the algebraic form:

[𝑀 * + 𝜏𝜃𝐾 * + 𝜏𝜃𝑀 ./*0 * + 𝜏𝜃𝑀 12+ * ]𝑇 * = 𝑀 * 𝑇 *>? -𝜏(1 -𝜃)𝐾 *>? 𝑇 *>? + 𝜏𝜃𝐹 * + 𝜏(1 -𝜃)𝐹 *>? -𝜏(1 -𝜃)𝑀 ./*0 *>? 𝑇 *>? -𝜏(1 -𝜃)𝑀
Where:

𝑀 * is the mass matrix:

𝑀 AC * = 𝜌𝐶 𝑝 𝑛 (𝑣 A , 𝑣 C ) @ ;
𝐾 * is the stiffness matrix:

𝐾 AC * = 𝑘 𝑛 (∇𝑣 A , ∇𝑣 C ) @ ; 𝑀 ./*0 * is the convection boundary matrix: 𝑀 ./*0 * AC = ℎ 𝑛 (𝑣 A , 𝑣 C ) D "#$% ; 𝑀 12+ * is the radiation boundary matrix: 𝑀 12+ * AC = 𝜖 * 𝜎(𝑇 * + 𝑇 29: )(𝑇 *E + 𝑇 29: E )(𝑣 A , 𝑣 C ) D &'( ; 𝑉 ./*0 * is the convection boundary vector: 𝑉 ./*0 * A = ℎ 𝑛 (1, 𝑣 A ) D "#$% ; 𝑉 12+
* is the radiation boundary vector:

𝑉 12+ * A = 𝜖 * 𝜎(𝑇 * + 𝑇 29: )(𝑇 *E + 𝑇 29: E )(1, 𝑣 A ) D &'( ;
𝐹 * is the right-hand-side vector: 𝐹 A * = (𝑄 * , 𝑣 A ) @ ; 𝑓 * is the heat flux boundary vector:

𝑓 A * = (𝑞 3456 * , 𝑣 A ) D )*+, .
Clearly, equation ( 10) is a set of fully discrete nonlinear equations. Next step is to linearize the above equation. Then the linear equation can be solved using general solver. An iterative procedure will be needed to find the converged solution.

Systematic design of the simulation scheme

This research develops a new efficient simulation scheme for SML process by using the open source FE library named Deal.II, which is mainly developed by Bangerth, et al. [START_REF] Bangerth | deal. II-a general-purpose object-oriented finite element library[END_REF]. Figure 2 presents the simulation scheme for the SLM process. Initially, the whole domain is coarsely meshed without assigning the degree of freedom (DoF) to each node. At the same time, a scanning path file is generated based on the CAD model and process parameters. The scanning path file includes process parameters such as laser power, scan orientation, speed, and time step for each scan track and layer. After the coarse mesh and scanning path file is generated, combined with the material property file, the simulation process starts to simulate the building process track-by-track and layer-by-layer. These scan tracks and layers are predefined in the scanning path file. For example, in the simulation of the current layer, cells belong to current layer are activated by assigning corresponding DoFs to each cell. Then, boundary conditions and powder material phases are updated and applied to newly activated cells. Once the activation of new layer cells and the updated domain is cooling for dwell time, the simulation process continues until all layer nodes are visited. Next, key modules which have been highlighted in Figure 2 will be discussed.

Modeling of the scanning path

As laser spot changes rapidly with respect to the time and space, the input heat flux is not easy to apply to the building part. An XML format file, which can identify the position of the laser spot at any time, will be used to drive the simulation process layer by layer at different time and position. Figure 3 illustrates the structure of the scanning path file which is generated based on the input CAD file and process parameters. The machine node mainly consists of process parameters for all layers. The layer node includes the process parameters for the current layer. Similarly, the scan track node includes necessary parameters for each scan track. This file drives the simulation process at different time steps and can trace the position of the laser spot according to the current time.

Modeling of material property

Thermal conduction is the main heat transfer mechanism in the SLM process. Thermal conductivity of powder bed is a function of density or porosity [START_REF] Kundakciolu | Thermal and molten pool model in selective laser melting process of Inconel 625[END_REF]. Yagi et al. [START_REF] Yagi | Studies on effective thermal conductivities in packed beds[END_REF] derived the effective thermal conductivity of the packed powder bed using this equation

𝑘 )33 = F & G - ?HIG -/G . ( 11
)
Where 𝑘 7 and 𝑘 K is the solid thermal conductivity and surrounding gas thermal conductivity, respectively. 𝜌 1 is the relative density of powder bed. Φ is the empirical coefficient taken as equal to = 0.02 × 10 E(L.N>F & ) . Usually, the thermal conductivity of the solid material is experimentally determined at different temperatures.

Another important thermal property is the specific heat, which is highly nonlinear with the change of temperature and phase transformation. During the melting and solidification process, the phase change is modeled based on the apparent heat capacity method [START_REF] Hu | Mathematical modelling of solidification and melting: a review[END_REF]. In this method, the specific heat is artificially increased around the solidus and liquidus point to capture the absorbed/released latent heat during the melting/solidification process. For example, the solid material requires more energy to maintain the temperature around the melting point and becomes liquid material. This process is equivalent to include the latent heat during the melting process. The effective specific heat capacity w.r.t temperature is modeled to incorporate the latent heat during melting/solidification [START_REF] Wolff | A framework to link localized cooling and properties of directed energy deposition (DED)-processed Ti-6Al-4V[END_REF].

𝐶 ! )33 (𝑇) = R 𝐶 ! (𝑇) 𝑇 < 𝑇 7 𝐶 ! + ∆𝐻 -/(𝑇 4 -𝑇 7 ) 𝑇 7 < 𝑇 < 𝑇 4 𝐶 ! (𝑇) 𝑇 4 < 𝑇 (12)
Where 𝐶 ! is the specific heat capacity, 𝑇 7 and 𝑇 4 is the solidus temperature and liquidus temperature, respectively. ∆𝐻 -is the change in enthalpy during the phase change and the enthalpy of the solid or powder is defined as,

𝐻 -= ∫ 𝜌𝐶 ! 𝑑𝑇 (13)
In this research, temperature-dependent material properties are used to capture the rapid change of material properties with the variation of temperature. In addition, each cell will experience the same material phase change from powder, to melted liquid, and finally to solid during the melting and solidification process as shown in Figure 4. Initially, the material phase type of a cell owns powder phase. Then, the material phase type transits to melted liquid when the temperature rises to the melting point. Finally, the liquid material solidifies when the temperature cools down to the room temperature. After that, the solidified material can only transform between remelted liquid and solid when the laser beam scans neighboring scan tracks or subsequent layers. It is not convenient to illustrate the relationship between porosity and temperature during the building process. Therefore, material properties are modeling as a piece function of temperature.

For each active cell, the material property keeps the same which is a function of the average cell temperature. A material property file is defined to tabulate the desired material property with prescribed temperature. For a given temperature 𝑇, the property 𝑀 can be obtained as,

𝑀 ! O (𝑇) = R 𝑓 ? (𝑇) 𝑇 ≤ 𝑇 ? 𝑓 E (𝑇) 𝑇 ? < 𝑇 < 𝑇 E 𝑓 P (𝑇) … (14)
Where 𝑝 indicates powder, solid, or liquid material type. 𝑞 indicates thermal conductivity, specific heat, convection coefficient, radiation coefficient, surface emissivity, and so on. 𝑓 ? , 𝑓 E , 𝑓 P presents the linear piecewise function under different temperature ranges respectively.

Powder Liquid Solid

T solid > T melt T liquid < T melt

Modeling of material deposition

It has been known that hundreds of layers are required to build a complex part, which results in a huge amount of computational cost. To simulate the layer-by-layer fabrication process, two different strategies are widely used to mimic the material spreading process in the SLM: the quiet cell method and the inactive cell method [START_REF] Debroy | Additive manufacturing of metallic components -Process, structure and properties[END_REF][START_REF] Afazov | Distortion prediction and compensation in selective laser melting[END_REF]. For the quiet cell method, the whole continuous domain meshes into finite element domain at the beginning. However, material properties of the non-deposited cells are prescribed to a very small value so that the non-deposition material will not affect the thermal analysis. For example, the conductivity and specific heat are assigned with a small value which will prevent the thermal heat from transferring into quiet cells. For the inactive cell method, after the whole domain meshes, only nodal DoF corresponding to active cells are assigned while non-active cells have no DoF. Therefore, the total cell number is changed with the activation of each new layer of material and the assemble algebraic system is smaller compared with the quiet cell method. Both methods have merits and drawbacks. For example, in the quiet cell method, the number of equations does not change and no need to renumber cells. However, if the value of material property of the quiet cell is too small, the final assemble mass matrix or stiffness matrix may be ill-conditioned and results in non-convergence of the finite element system. If it is too large, it will lose accuracy due to the transfer of heat into the quiet cell and the accumulated errors in each time step. In the inactive method, only active cells have DoF. Therefore, the finite element system is smaller at a certain layer compare with the quiet cell method. In addition, no ill-condition of the assemble matrix is introduced by inactive cells. However, extra efforts are required to renumber the cell and node as the number of DoF always change.

In this work, a new cell activation strategy is proposed by combining the advantages of quiet cell method and inactive cell method. Figure 5 illustrates the simulation of the material deposition process in detail. First, the geometric domain is meshed with the same size of hexahedron cells in 3D without DoFs. Second, a new layer of cells is activated with DoFs in each cell and an initial powder material phase is assigned to newly activated cells. Then the simulation process in the current layer proceeds without changing the size of the assembled matrix until the next layer of cells are activated. This process continues until the last layer of material is activated and simulated.

Once the temperature of powder material exceeds melting temperature, the phase of the material transfers into liquid and then changes to solid as it is cooled down. When a new layer of material is activated, a thermal analysis of the cooldown process will be executed to simulate the powder spreading process. The bottom surface is in contact with the substrate where Dirichlet boundary condition is applied. The lateral surface is immersed in the powder therefore thermal isolation condition is applied. The thermal irradiation and convection boundaries, as well as the surfacerelated heat source, are applied to the top external free surface.

Modeling of the heat source

Gaussian point heat source (GPHS) model has been widely for FE modeling of the SLM process. Hussein, et al. [START_REF] Hussein | Finite element simulation of the temperature and stress fields in single layers built without-support in selective laser melting[END_REF] and Ren, et al. [START_REF] Ren | Simulation of transient temperature field in the selective laser sintering process of W/Ni powder mixture[END_REF] assumed that if 85% of total thermal energy is captured within a circle projected by the laser beam, the GPHS can be derived as,

𝑞(𝑥, 𝑦, 𝑡) = EQR ST / exp[- EU6 / HV / W T / ] ( 15
)
Where 𝜂 is the efficiency of the laser power, 𝑃 is the laser power, 𝜔 is the radius of the laser beam.

Equation [START_REF] Simson | Residual stress measurements on AISI 316L samples manufactured by selective laser melting[END_REF] is the Gaussian point heat flux distribution in Cartesian coordinates. In the Lagrangian frame, the position of the laser beam is related to time, which can be written as

F 𝑥 = 𝑥 L + 𝑣 6 𝑡 y = 𝑦 L + 𝑣 V 𝑡 (16)
where, 𝑣 6 , 𝑣 V is the scan speed in x and y direction respectively.

Usually, the laser beam scans the powder bed very quickly. Therefore, Irwin and Michaleris [START_REF] Irwin | A line heat input model for additive manufacturing[END_REF] made an assumption that the moving of the GPHS forms Gaussian line heat source (GLHS). For the GPHS, the timestep ∆𝑡 should satisfy the relation of ∆𝑡 ≤ 𝜔/𝑣. Otherwise, the GPHS will skip over some cells and induce incorrect solutions. Therefore, a large number of timesteps is mandatory to simulate the moving GPHS along scan tracks. To get rid of these limitations, equation ( 15) is averaged over the period of a timestep ∆𝑡 and yields the GLHS model, 𝑞 l(𝑥, 𝑦, 𝑡) = 1/∆𝑡 ∫ 𝑞(𝑥, 𝑦, 𝑡)𝑑𝑡

' 0 H∆' ' 0 (17)
Substitution of equation ( 15) into equation ( 17) yields the GLHS model,

𝑞 l(𝑥, 𝑦, 𝑡) = QR ∆'•0•T•√ES exp m- EV 1/ T / n erf ( √E6 1 T )| 'B' 0 ' 0 H∆' (18)
where, Where 𝛼 is the scan orientation of a scan track.

r 𝑥 [ = 𝑥𝑐𝑜𝑠𝛼 + 𝑦𝑠𝑖𝑛𝛼 = (𝑥 -𝑥 L -𝑣 6 𝑡)𝑐𝑜𝑠𝛼 + (y -𝑦 L -𝑣 V 𝑡)𝑠𝑖𝑛𝛼 𝑦 [ = -𝑥𝑠𝑖𝑛𝛼 + 𝑦𝑐𝑜𝑠𝛼 = -(𝑥 -𝑥 L -𝑣 6 𝑡)𝑠𝑖𝑛𝛼 + (y -𝑦 L -𝑣 V 𝑡)𝑐𝑜𝑠𝛼 (19)
Equation [START_REF] Denlinger | Residual stress and distortion modeling of electron beam direct manufacturing Ti-6Al-4V[END_REF] shows that the timestep can be enlarged without jumping any cells. Also, as the timestep ∆𝑡 → 0, the GLHS model becomes the moving GPHS model. Figure 6 (a) and (b) illustrate the GPHS and the GLHS along a single scanning track. The application of GLHS instead of GPHS can significantly reduce the computation cost due to the increased timestep. Obviously, the replacement of moving GPHS model with the GLHS model would lead to some extent of accuracy problem. To overcome this issue, a hybrid of Gaussian point and Gaussian line heat source model is proposed in this research. Specifically, for each scan track, the GPHS is applied at the start and end of the scanning track. While the GLHS is applied to the middle section as shown in Figure 6 (c). However, length of the line heat input cannot be too large as a significant amount of error will be introduced. The previous research [START_REF] Luo | Numerical simulation of temperature fields in powder bed fusion process by using hybrid heat source model[END_REF] has detailly investigated the influence of line heat input length and position of line heat input on temperature fields among the powder bed. It concludes that when the length of each line heat input is 1 mm and moving GPHS is applied for 1 mm around the boundaries for each track, the accuracy and simulation speed can reach a good balance. In this research, for a given scanning track vector of 𝑺𝑬 zzzzz⃗ that starts and ends at point S and E respectively, and the coordinate of current laser spot of 𝑷(𝑥, 𝑦), the hybrid heat source (HHS) 𝑞 }(𝑥, 𝑦, 𝑡) is modeled using the following rules (Figure 7).

Adaptive mesh refinement scheme

The small diameter of the laser beam necessitates high dense mesh around it. Otherwise, there are not enough quadrature points to capture the input heat flux. On the other hand, the density of cells far away from the heat source should be decreased to reduce the scale of the FE system. In addition, with the motion of the laser beam, the total mesh should be adaptively refined and coarsen as the laser beam travels. Figure 8 (a) shows that the coarse mesh cannot capture the maximum heat input due to lack of enough quadrature points. Figure 8 (b) shows that the global refined mesh has enough quadrature points to capture the Gaussian distributed heat input by increasing the scale of the FE mesh. However, Figure 8 (c) shows that the adaptively refined mesh can accurately capture the Gaussian distributed heat input without significantly increasing the scale of the FE system.

There are mainly three different adaptive mesh methods: h-adaptive, r-adaptive, and p-adaptive [START_REF] Ghorashi | Goal-oriented error estimation and mesh adaptivity in 3d elastoplasticity problems[END_REF]. The h-adaptive method discretizes an element into four quadrilaterals in 2D and 8 hexahedra in 3D without altering shape functions. The r-adaptive method finds the optimum mesh topology for a given discretization. The p-adaptive method increases the polynomial order of shape functions of an element without discretizing it. Due to the complex geometry of the AM part, the r-adaptive method is not suitable. The mesh close to the laser spot is typically refined 3 to 5 levels, the p-adaptive method in this case will introduce more hanging nodes between two adjacent elements in which extra efforts are needed to constrain these hanging nodes. Therefore, the hadaptive method is considered in this research to refine and coarsen the mesh based on the estimated element-wise errors.

The Deal.II FE library has the ability to easily implement the mesh refinement and coarsen operation [START_REF] Riedlbauer | Thermomechanical finite element simulations of selective electron beam melting processes: performance considerations[END_REF]. The mesh refinement criterion is based on the Kelly error estimation, which is developed by Kelly, et al. [START_REF] Kelly | A posteriori error analysis and adaptive processes in the finite element method: Part I-error analysis[END_REF] for the Laplace equation. This error estimator utilizes the jump of solution gradient at cell faces and scales it by the size of the cell. The implementation of this error estimator for a cell is defined as

𝑒 \ E = ∑ 𝐶 3 ∫ m 2"5 2 "𝐧 n E "\ ) 𝑑𝜊 3∈"\ (20)
Where 𝑎𝜕𝑢 < /𝜕𝐧 denotes the solution gradient's jump at cell faces. 𝑓 is the cell face. 𝐶 3 is the scale of the cell size and usually can be taken as the greatest length of the diagonals of cell face.

The mesh will be refined in the region where has relatively large Kelly error and coarsen in the region where has relatively small Kelly error. To maintain the geometry boundary, another coarsen criterion is considered: if 8 adjacent and fine cells have the same material phase, and the Kelly error of each cell is small, these 8 cells can be coarsened as one cell. Otherwise, these 8 cells are forbidden to be coarsened. This criterion guarantees that the boundary of the part is distinguishable from the surrounding powder bed. The adaptive mesh will introduce hanging nodes which lie on the edge of two neighboring cells as shown in Figure 8 (c). For example, the solution of node N5 or N9 is half of the sum of neighboring nodes to make sure a continuous solution. This section presents a simple case to validate the temperature field in the melt pool under the GPHS and GLHS. In addition, simulation efficiency of different simulation schemes is compared.

A scanning track travels along the center of the top surface is simulated. Figure 9 describes the initial coarse mesh and the initial material phase of each cell. This mesh will be adaptively refined and coarsen based on the estimated Kelly error of each cell [START_REF] Kelly | A posteriori error analysis and adaptive processes in the finite element method: Part I-error analysis[END_REF]. Figure 10 and Table 1 give the temperature-dependent material properties of Stainless Steel 316L and process parameters modified from Gusarov, et al. [START_REF] Gusarov | Model of radiation and heat transfer in laser-powder interaction zone at selective laser melting[END_REF] and Tran and Lo [START_REF] Tran | Heat transfer simulations of selective laser melting process based on volumetric heat source with powder size consideration[END_REF]. Note that the specific heat is artificially increased around the melting point to encompass the latent heat.

Efficiency analysis

This research utilizes two approaches to improve the efficiency of proposed simulation scheme. The first approach is the Kelly error based adaptive mesh scheme which can greatly reduce the cell number while keeping the solution accuracy; The second approach is the development of new GLHS model which can significantly increase the step size and thus reduce the simulation time. This section mainly compares the computation time and solution accuracy of different simulation schemes.

Temperature fields at different time steps and under different simulation schemes are illustrated in Figure 11. Figure 11 (a)(b) and (c)(d) applies static mesh and adaptive mesh respectively. The number of active cells can be tremendously decreased by using the adaptive mesh strategy compared with the static mesh strategy which has been adopted by many researches [START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF][START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF][START_REF] Jayanath | A Computationally Efficient Finite Element Framework to Simulate Additive Manufacturing Processes[END_REF].

In addition, the simulation process can start with a coarse mesh and then activate the material layer by layer without a dense mesh at the beginning. Therefore, the simulation cost can be reduced. Figure 11 (a)(c) and (b)(d) applies the moving GPHS and the GLHS respectively. A relatively large time step size can be used for the GLHS compared to the time step size used for the GPHS. In this validation case, the time step is set to 0.001s for the GPHS and 0.01s for the GLHS. Only three segments of line heat input can be used to simulate this 5 mm scanning track. Therefore, the application of the GLHS model can greatly reduce the simulation cost. Figure 12 shows the temperature profile along the center of the scanning track at different time. It can be noted that the temperature distributions under different heat sources around the melt pool are quite different. The maximum temperature of the GLHS is much lower than that of the GPHS as the total input heat flux is uniformly distributed along the line according to equation ( 17) and [START_REF] Denlinger | Residual stress and distortion modeling of electron beam direct manufacturing Ti-6Al-4V[END_REF]. However, the temperature distribution far from the melt pool is very close to the temperature distribution under the GPHS. This indicates that for the material under the GLHS far from the melt pool experience the same thermal cycles with the material under the GPHS. When the laser beam travels to the end of the first line heat input (t = 10.4 ms as shown in Figure 12 (a)), the material around the line heat input is melted. As the laser beam moves to the end of second line heat input (t=20.8 ms as shown in Figure 12 (b)), the temperature of the material around the first line heat input is almost the same to that of the GPHS. Similarly, both of the temperature distributions around the first and second line heat inputs are very close to that of the GPHS when the laser beam moves to the end of third line heat input (t=31.2 ms as shown in Figure 12 (c)). In summary, as the powder material will be melted under the GLHS and the size of the melt pool is much smaller compared to the part dimension, all material will experience the same thermal cycle after the melted material is cooling from the melted pool. As shown in Figure 12, the temperature profiles of the melt pool under the GPHS and the GLHS are quite different. To evaluate the temperature difference between the GPHS and the GLHS, the following equation is used to evaluate the average temperature difference over the scanning track between the GPHS and the GLHS:

𝜖(𝑡 * ) = ' ? 4 ∑ m #(6 0 )># 1 (6 0 ) #(6 0 ) n E ∆𝑙 4 6 0 BL × 100%, 𝑡 * ∈ [0, 𝑙 𝑣 ⁄ ] (21)
where 𝑇 and 𝑇 [ is the temperature obtained from the simulation scheme using the GPHS and the new GLHS, respectively. 𝑡 * is the time step number. 𝑥 L is the position of the sampling point. 𝑙 is the length of this single scanning track and ∆𝑙 = 0.1 𝑚𝑚. 𝜖(𝑡 * ) is the averaged percent error of temperature difference along a scanning track at time step 𝑡 * . Then the corresponding 𝜖 at time step 10.4 ms, 20.8 ms, and 31.2 ms is 11.23%, 4.38%, and 8.76%, respectively. Note that the temperature difference around the start and end of a scanning track is relatively large. Therefore, hybrid application of the GPHS to the start and end of a scanning track and the GLHS to the middle of a scanning track is a potential way to decrease the temperature difference.

Based on the simulation results as shown in Figure 11 and Figure 12, Table 2 compares the simulation efficiency and temperature difference of a single scanning track between different simulation schemes. The solution of the moving GPHS with static mesh is taken as reference.

Compared to the GPHS, the time step number for the GLHS can be significantly decreased due to increased time step size. Compared to the static mesh strategy, the active cell number of adaptive mesh strategy can be reduced by more than half in this case. However, the reduced computational time is not proportional to the decreased cell number. This can be attributed to the fact that extra operations are needed to transfer solution between two different meshes. For the GLHS, the temperature difference between the GPHS and the GLHS is confined below 10% while the efficiency can be improved by more than 2 times. Therefore, the adaptive mesh scheme combines with the GLHS is a good balance between the solution accuracy and computation cost.

Melt pool validation using the GPHS

Based on the process parameters from Table 1, scanning tracks with different scanning speeds are experimentally printed on the substrate and the experimental results will be used to validate the simulated results. Figure 13 (a) and (b) show the experimental images of the melt tracks and corresponding cross-sections conducted by Gusarov, et al. [START_REF] Gusarov | Model of radiation and heat transfer in laser-powder interaction zone at selective laser melting[END_REF]. Figure 13 (c) ~ (e) present the material phase, cross-section, and top view of the simulated melt pool respectively using moving GPHS. Both experiment and simulation results show that melt pool width and contact width are reduced with the increase of scanning speeds. In addition, as the scanning speed is over 200 mm/s, experimental results show no contact between the melt pool and the substrate. Though the simulated contact width is not zero at the scanning speed of 200 mm/s, the depth is almost zero. This is the reason that leads to a broken melt track which can be observed in the top right of Figure 13 (a). Table 3 summarized the experiment and simulation results for melt pool width and contact width. It can be noted that the dimension of the simulated melt pool and experimental melt track obtained a good agreement with 11% ~ 25% error. Therefore, the proposed simulation scheme can predict the dimension of the melt pool with acceptable error in the SLM process.

(a) (b) (c) (d) (e) 165um 140um
Contact width in present study (um) 137 88 30 12 Experimental contact width of Gusarov, et al. [START_REF] Gusarov | Model of radiation and heat transfer in laser-powder interaction zone at selective laser melting[END_REF] (um) 120 100 --

Melt pool validation using the GLHS

This section compares the temperature field of melt pool under the GLHS and the GPHS. The same process parameters are used. In addition, the length of line heat input is set to 1 mm, the time step is set to 0.01s, and scanning speed is set to 160 mm/s. The simulated temperature distributions are presented and compared in Error! Reference source not found. (a) and (b). Cross-sections of the two simulated scanning tracks obtained a good agreement except that the length of melt pool under the GLHS is significantly different. Considering the dimension of a part in the real printing process, the influence of the melt pool in length direction on the final temperature distribution of a part can be ignored, which can be proven in Error! Reference source not found.. In Error! Reference source not found., the temperature history curve of point p1 within the GLHS tends to be the same with the curve of the GPHS after 0.03 second. Table 4 summarizes the simulation results by using two heat sources. It can be noted that both of point p1 and p2 can achieve melting point within the moving GPHS and GLHS. However, the simulation time by using GLHS is much lower than the simulation time by using moving GPHS as shown in Table 2. Therefore, GLHS is suitable to replace the GPHS when modeling the temperature history and material transition among the melt pool. The simulation time can be reduced greatly without inducing too much temperature deviation.

Numerical application

The validation of the melt pool and simulation of several scanning tracks cannot show the ability of the proposed simulation scheme. Therefore, a numerical simulation of the SLM process in partlevel is done as a case study to investigate the thermal response of a "T" part with overhang features. The purpose of this case study is to indicate whether the proposed FE model can predict the temperature distribution and history in part-level. The obtained temperature fields of each scanning track and each layer are very useful. For designers, the predicted temperature fields of a part can be used to guide the optimization of geometrical modeling during the design stage; for manufacturers, the predicted heat concentration area may need more material to conduct the heat or prevent from warpage so manufacturers have more confidence in the manufacturability; for academic users, the predicted temperature data can be used for further simulation of displacement/stress fields and microstructure evolution in part level and optimize the process parameters without conducting a real experiment.

Simulation design

A thin-wall part with two different overhang features is simulated using the proposed simulation scheme. As shown in Figure 16, the dimension of the whole domain is 14.4 × 14.4 × 12.8 𝑚𝑚 P . The substrate domain is 6.4 mm in height. The thin wall is 1 mm in thick and 5 mm in height, and it is immersed in the build domain (14.4 × 1.6 × 6.4 𝑚𝑚 P ). The two overhang features have the same size except that the left one is inclined 15 ° w.r.t. the horizontal plane. The layer profile of the part is not visible until the coarse mesh is refined at the boundary of the part. The origin of coordinates is in the center of the top substrate surface as shown in Figure 16. Scanning tracks are driven by an input XML file which consists of all scanning path and heat source information. A simple zigzag scanning pattern is used in this application and there are 100 layers in total with a layer thickness of 50 um.

Figure 16. The dimension of the build domain and part domain

In this application, thermal convection, thermal radiation, and surface-related heat flux are applied on the top surface. Note that the top surface is continuously updated when a new layer material is activated. The height of deposited material varies with the increase of new layer of material during the building process. The lateral surface is considered as adiabatic since it is surrounded by powder material which has a very low thermal conductivity compared to the solidified material. The substrate temperature of the bottom surface is fixed to 300 K. The powder material is the stainless steel 316L and its temperature-dependent thermal properties are described in Figure 10. Other process parameters are listed in Table 5. The maximum size of the coarse mesh is 1.6 mm. The layer thickness is 50 um. Therefore, the critical mesh size close to the heat affected zone is 50 um, which is too small to mesh the whole domain in part level. To minimize the number of cells, adaptive mesh scheme and material activation scheme are applied to reduce the scale of the FE system. Figure 17 presents the result of the adaptive mesh and material deposition scheme in the middle/end of the printing process.

The gray inactive elements denote non-deposited powder material while the red active elements denote deposited material. As shown in Figure 17 (a) and (b), the coarse mesh is found in the substrate domain and far from the current construction layer, while the region close to the laser beam is continuously to capture the detailed thermal diffusion. In addition, as shown in Figure 17 (c), the mesh close to the boundaries/edges keeps refined from coarsening to maintain the part shape, and the coarse mesh can be observed within the part. Utilizing the presented mesh scheme with the deposition of new layers, the number of total active cells in the simulation system can be kept between 116k and 125k after 100 layers material is deposited when using moving GPHS.

Simulation results and analysis

After deposition of 100 layers and cooling down for 1 s in each layer, the simulated temperature distributions in the entire domain under three different heat sources are given in the upper of Figure 18. The part domain (consists of solid elements) is shown in the bottom of Figure 18 after removing powder elements. For the part domain, when using GPHS, the highest temperature after cooling down for 1 second is around 1060 K close to the end heating point of the previous layer. The lowest temperature is around 430 K in the corner of the bottom part connected to the substrate. As a comparison, the highest temperature by using the GLHS and the HHS is 1180 K and 1115 K respectively, and the lowest temperature is around 440 K and 435 K respectively.

Table 6 lists the number of cells and simulation time for different heat sources. All the three simulations are run on an Intel Core i7-6700K @ 4GHz processor. It can be noted that the simulation time can be greatly reduced by using the GLHS or HHS when it is compared to the GPHS. As shown in Figure 19, temperature distributions along a certain curve under different heat sources are very similar. Though a difference between the GPHS and the GLHS can be observed, the trend of the temperature gradient is almost the same. Considering the temperature distribution along HT1, the temperature tends to be uniform among the center solidified region while non-uniform among the surrounding powder material. Considering the temperature distribution along HT2 and HT3, temperature gradients tend to be lower in both left and right overhang than that of the center solidified region. This can be attributed to the large solidified material in the center region, which can be served as a heat sink. As the powder material resides below overhangs, the thermal heat in the overhang region cannot be dissipated as fast as that in the center region and higher temperature can be observed in the overhang regions. In addition, the temperature in the right inclined overhang is higher than the temperature in the left overhang due to less solidified material in the right overhang. For the temperature distribution along vertical lines, the distinct transition can be found from the powder material to solidified material in the line VT2 and VT3. It also can be observed that temperature difference among the powder material (VT2/VT3) is higher than the temperature difference among the solidified material (VT1).

Temperature history of testing points:

To evaluate the temperature evolution of the substrate and part, the temperature history of point A ([0,0,0] mm), point B ([-0.8,-0.8,0] mm), point C ([3,0,3] mm), and point D ([-3,0,3] mm) are recorded after each time step and drawn in Figure 20. It can be noted that each scanning track will be re-melted by adjacent track and influence by several layers. As shown in Figure 20 (a), when the laser beam approaches point A, the temperature of point A is changed rapidly for consecutive three tracks. And this influence decreases as the increase of deposition layer. For the GPHS, the temperature of point A and point B oscillate around 515 K and 440 K with the increase of deposited material, respectively. The average temperature of point B is lower than that of point A is because point B is in the powder domain and less thermal heat is conducted to the powder domain. The temperature of point C and point D oscillate around 1120 K and 1100 K with the increase of deposited material, respectively. Though the point C and point D is symmetric about the z-axis, the relatively high temperature of point C can be attributed to the low thermal conductivity of powder material which accumulates the heat for a longer time compared to the solid material. Similar temperature history trends for those points can be found in the HHS and the GLHS except that there is a small difference about the converge temperature. Therefore, it can be concluded that each point experiences the same thermal cycles in the HHS and GLHS compared to the thermal cycles of corresponding point in the GPHS.

Conclusions

Numerical modeling of the SLM process in part level has been a challenge for decades due to the unacceptable computation cost for simulating a real part that using a moving Gaussian point heat source. To solve this problem, this research proposes an efficient simulation scheme for the SLM process. A new Gaussian line heat source and hybrid heat source are used to accelerate the simulation process by adopting a large acceptable time step. A material deposition scheme, which can reduce the active cell number, is applied to mimic the powder spread process. In addition, an adaptive mesh scheme is applied to refine the mesh around the region where heat input is applied and coarsen the mesh where has a low-temperature gradient. To simulate the printing process accurately, temperature-dependent material properties and material phase transition are The simulation efficiency of this new simulation scheme can be improved by more than 4 times using the Kelly error based adaptive mesh and the GLHS compared to the traditional scheme that utilizes static mesh and the moving GPHS. While the average temperature difference along a scanning track is below 10%. In addition, the difference of the melt pool dimension between the numerical result and experiment result can be limited to an acceptable error of 11% ~ 25%. More importantly, this simulation scheme for part level simulation of the SLM can be run on a general workstation rather than the HPC. Numerical experiment results of the 'T' part show that temperature distribution tends to be more uniform in the solidified region compared to the powder region. Therefore, larger temperature gradients can be found among the overhang feature where it is surrounded by less solid material compared to the region that has more solid material. The simulation result reveals that hybrid of Gaussian point heat source and line heat source can be used to replace the Gaussian point heat source without losing too much accuracy while greatly reduce the simulation time. In addition, for a certain point, it will experience the same thermal cycle under different heat sources.

The developed FE model can be used to analyze and optimize the process parameters such as scanning speed, layer thickness, and scanning path without setting up a real printing experiment.

The estimated temperature history of each point also can be used for further material thermomechanical behavior and microstructure evolution in part level. This new simulation scheme for the SLM process can be easily extended to other powder-bed-fusion AM process by making some small modifications. Future work should be focused on the reduction of cell elements and assembly work as well as iteration numbers when solving nonlinear heat transfer equations.

Figure 2 .

 2 Figure 2. Systematic simulation scheme for SLM process

Figure 4 .

 4 Figure 4. Material transition relationship

Figure 5 .

 5 Figure 5. Material deposition strategy

Figure 6 .

 6 Figure 6. Heat sources and scanning orientation

Figure 7 .

 7 Figure 7. Rules for applying the hybrid heat source of a single scanning track

Figure 8 .

 8 Figure 8. Locally refined adaptive mesh and constraint rules for hanging nodes

Figure 9 .Figure 10 .

 910 Figure 9. Description of initial FE Mesh and material phase

Figure 11 .

 11 Figure 11. Temperature distributions under different simulation schemes

Figure 12 .

 12 Figure 12. Temperature profiles along the center scanning track under different time steps

Figure 13 .

 13 Figure 13. Comparison of experiment results for melt track and temperature distributions of melt pool at various scan

Figure 14 .Figure 15 .

 1415 Figure 14.Comparison of temperature distributions in different heat source

Figure 17 .

 17 Figure 17. Adaptive mesh and material deposition at the middle/end of the building process

Figure 18 .

 18 Figure 18. Temperature distribution at the end of the building process after cooling for 1 second 5.2.1 Temperature distribution along testing curves as the build domain is very thin in the thickness direction, only the temperature distribution along horizontal and vertical directions are presented here. The temperature distribution along selected curves are presented in Figure 19. The vertical line VT1 ([0,0,0] ~ [0,0,5] mm) crosses from the center of the top substrate to the center of the top part. The vertical line VT2 ([-4,0,0] ~ [-4,0,5] mm) and VT3 ([4,0,0] ~ [4,0,5] mm) cross from the top substrate to the left and right overhang, respectively. The horizontal line HT1([-7.2,0,1.6] ~ [7.2,0,1.6] mm) crosses the bottom cooled region. The horizontal line HT2 ([-7.2,0,3.2] ~ [7.2,0,3.2] mm) crosses the two overhang features. The horizontal line HT3 ([-7.2,0,4.8] ~ [7.2,0,4.8] mm) crosses the top hot region. The transition between different materials is denoted by small circles as shown in Figure 19 (b) and Figure 19 (c). As the heat source will melt some powder material around boundaries, the final solidified domain will be a little larger compared to the designed part domain, which can be noted in the transition point from powder to solid in Figure 19 (b).

Figure 19 .

 19 Figure 19. Temperature distribution along different curves through the part domain after a dwell time

Figure 20 .

 20 Figure 20. Temperature history of testing points under different heat sources

  history of point A, B, C, and D with the GPHS (b) Temperature history of point A, B, C, and D with the HHS (c) Temperature history of point A, B, C, and D with the GLHSincorporated in this model. The shape and temperature distributions of the melt pool under different heat sources are comparable to other research and experiment results.

  Figure1illustrates the simulation domain and boundaries during the SLM process. The boundary 𝛤 , is the bottom surface of the substrate where Dirichlet boundary is applied. The boundary 𝛤 -is the lateral surface of the whole domain where is surrounded by the loose powder and thermal insulation is applied there. The boundary 𝛤 # is the top surface of each layer indicates a Robin boundary. Thermal convection (𝛤 ./*0 ), radiation (𝛤 12+ ), and input heat flux (𝛤 3456 ) boundaries will be applied on top of each layer. These three boundaries satisfy the relation of ∂𝛺 = 𝛤 , ∪ 𝛤 -∪ 𝛤 # , and 𝛤 # = 𝛤 ./*0 = 𝛤 12+ = 𝛤 3456 .

	Figure 1. Build domain and boundary conditions
	Boundaries in this research are listed below
	⎩ ⎧ ⎨ -𝑘 • 𝑇(𝒙, 𝑡) = 𝑇 7 ∀𝒙 ∈ 𝛤 , , 𝑡 ∈ (0, 𝑡 )*+ ) -𝑘 • "#(𝒙,') "𝐧 = 𝑞 ./*0 + 𝑞 12+ + 𝑞 3456 "#(𝒙,') "𝐧 = 0 ∀𝒙 ∈ 𝛤 -, 𝑡 ∈ (0, 𝑡 )*+ ) ∀𝒙 ∈ 𝛤 # , 𝑡 ∈ (0, 𝑡 )*+ )

Scan track 2 …… …… Layer 1 Layer 2 Layer 3 Scan path Scan path file

	Root node
	Machine node
	total layers
	velocity
	beam diameter
	Layer 1
	Layer node
	total scan tracks
	scan orientation
	thickness
	hatching distance
	part height
	layer start point
	layer end point
	layer start time
	layer end time
	layer cooling time
	time step
	Scan track 1
	heat source type
	scan track length
	start point
	end point
	start time
	end time
	Layer 2

Table 1 .

 1 Process parameters used in the validation of melt pool

	Process parameters	Values
	Effective laser power	45 W
	Scanning speed, 𝑣	120, 160, 200, 240 mm/s
	Hatch spacing	50 um
	Spot size	50 um
	Layer thickness	50 um
	Material	Stainless Steel 316L
	Melting point	1700 K
	Latent heat	300 kJ/kg

Table 2 .

 2 Efficiency comparison and temperature difference between different simulation schemes

	Simulation scheme	Time step	Cell number	Wall CPU	Temperature	Efficiency
		size (s)		time (s)	difference (%)	improved (%)
	Static mesh + GPHS	0.001	129200	511	-	-

Table 3 .

 3 Comparison of experimental results with simulated results for melt pool width and contact width

	Scanning speed (mm/s)	120	160	200	240
	Simulated width in present study (um)	237	212	187	175
	Experimental width of Gusarov, et al. [50] (um)	200	190	165	140

Table 4 .

 4 Simulation time and melt pool dimensions within the moving GPHS and GLHS

	Heat source	Step size (s)	Melt pool dimension (um)	Dimension difference (%)
	GPHS	0.001	Width: 220
			Depth: 70	Width: 2.3
	GLHS	0.01	Width: 215	Depth: 7.1
			Depth: 65

Table 5 .

 5 Simulation process parameters

	Process parameters	Values
	Effective laser power	60 W
	Scanning speed, 𝑣	150 mm/s

Table 6 .

 6 Comparison of the simulation cost for different heat sources

		GPHS	HHS	GLHS
	Number of active cells	<125000	<167000	<175000

Acknowledgments

This research work is supported by National Sciences and Engineering Research Council of Canada Discovery Grant RGPIN 436055-2013 and China Scholarship Council.


Electrotechnical Commission (IEC). Power Performance Measurements of Electricity Producing Wind

 Turbines;IEC 61400-12-1 Ed. 2;International Electrotechnical Commission: Geneva, Switzerland, 20172. Albers, A., 2010.Turbulence and shear normalisation of wind turbine power curve. In: European Wind Energy Conference andExhibition 2010, EWEC 2010. Vol. 6. pp. 4116-4123. 3. Electrotechnical Commission (IEC). Power performance of electricity-producing wind turbines based on nacelle anemometry; IEC 61400-12-2;International Electrotechnical Commission: Geneva, Switzerland, 2013 4. Milan, Patrick, et al. "18 -Power Curves for Wind Turbines ." Wind Power Generation andWind Turbine Design, WIT Press, Billerica, MA, 2010. The global aim of the turbulence normalisation, or normalisation, is to correct the power signal and to reduce the bias induced, in part, by the mathematical averaging process of the dataset over a 10min period. The process can't be fully described here but is publicly accessible [2], [1]. The correction is applied as follow:

𝑃 𝐼 𝑟𝑒𝑓 𝑣 = 𝑃 𝑣 -𝑃 𝑠𝑖𝑚,𝐼 𝑣 + 𝑃 𝑠𝑖𝑚,𝐼 𝑟𝑒𝑓 𝑣
Where 𝑃 𝐼 𝑟𝑒𝑓 𝑣 is the normalised power output, 𝑃 𝑣 is the mean measured power output, 𝑃 𝑠𝑖𝑚,𝐼 𝑣 and 𝑃 𝑠𝑖𝑚,𝐼 𝑟𝑒𝑓 𝑣 are the power outputs for, respectively, the measured TI and the reference TI. These 𝑃 𝑠𝑖𝑚 are derived using:

P sim v = න 𝑣=0 ∞ 𝑃 I=0 𝑣 . f 𝑣 𝑑𝑣
Where 𝑃 I=0 𝑣 is the Zero-Turbulence PC function and f 𝑣 is the wind distribution within the 10min period. As the latter is never stored, f 𝑣 is almost always approximated by g 𝑣 , a standard gaussian function.

The determination of the 𝑃 I=0 𝑣 itself is an iterative process that involves the computation of an initial PC, its binned wind speed and TI as well as prior knowledge about the studied wind turbine. The TI correction is applied after the air density normalisation on each 10min point of the dataset.

Turbulence Intensity (TI or 𝐼) affects wind turbine power curve analysis when performed using 10-minute datasets. This influence has been well documented, its mathematical component is understood and its influence on the Power Curve (PC) characterised.

Over the years, methods have emerged to counteract the resulting bias on PC and, since its update in 2017, the IEC 61400-12-1 [1] has also added a method for turbulence intensity correction, annex M. Although informative, this annex provides a common ground to deal with this issue that should be used by the wind industry.

Unlike air density normalisation of the wind speed, the approach for TI correction can feel more convoluted. This perceived high complexity can be seen as a barrier to entry by data analysts and could explain its low adoption by the community.

In this study, we assess the effect of the IEC annex M methodology when applied to wind measurement from SCADA, rather than undisturbed met mast or lidar measurements.

Abstract Fundamentals: Turbulence Intensity Normalisation

Implementation

Assessing the effect on the Power Curve

When implementing the algorithm with the intention of performing monthly analysis on a >3000 WTGs portfolio, computational efficiency is key. To solve the integration of the TI correction, Monte-Carlo or quadrature implementations are simply too costly when deployed at this scale.

Fortunately, as long as the f 𝑣 is simplified by g 𝑣 or another defined fitted function, and 𝑃 I=0 𝑣 remains a function defined by part, an analytical solution to the process can be found. This closed solution is several orders of magnitude faster to run than any other implementation.

When analysing the monthly PC, the Power_Sdv is slightly but consistently reduced (-6.9% in the analysis above). When deploying that metric on whole portfolio, we found an average reduction of the binned Power_Sdv of about -5%. The main expected effect of the normalisation on the dataset is obviously to correct the power signal as if the TI was constantly equal to 𝐼 𝑟𝑒𝑓 . The method typically results in lower power values in the lower windspeed range, and higher power values at higher windspeeds close to the "knee" of the curve. This is because TI tends to be higher than the standard reference value (𝐼 𝑟𝑒𝑓 ), often set to 10% (or 0.1). The well-known visualisation of this effect can be seen o the left [4]:

Conclusions

• When undisturbed wind measurements are not available, the pragmatic solution for fleet scale performance analysis is to use nacelle-based anemometry from SCADA, despite the known influence of the blades.

• Given the results of this pragmatic approach, it was deemed suitable for immediate deployment and used for regular, automated analysis of the performance of over 3000 WTGs.

• The inclusion of TI correction has reduced the perceived deviation from expected performance, therefore reducing the number of "false alarms" raised by automated analysis.

• The quantity of monthly analysis performed forced us to implement the methodology in the most efficient way, in the form of its analytic solution.

Computation time and power required for Monte-Carlo or integration algorithm implementations was shown to be unfeasible for a large-scale deployment.

As stated above, we are here testing the TI normalisation using wind measurement from SCADA rather than the preferred met mast or lidar sources.

An additional benefit of the method is also the reduction of the scatter of the PC when observed across a period where TI had varied. Checking that this reduction did occur is an effective way to verify if the process had a positive impact or not . If that scatter remains the same or increases, it implies that the method had either no effect or a negative one on the dataset. On the other hand, any significant reduction on the scatter would indicate that the normalisation reduces uncertainty.

The graphs below compares the standards deviation (Power_Sdv) of the raw power and the corrected power for one of the WTGs in the portfolio.

Although deviation from the standards [2, 3], we perform careful filtering of the dataset to remove wake periods, icing, curtailment and faults, and we assume that TI derived from nacelle anemometry correlates with undisturbed TI. Since met masts are not available for all WTGs in our portfolio, we cannot automatically check the quality of data provided by the nacelle anemometry of each WTGs. However, we can check that the correction behaves as the theory dictates. An example of that verification can be seen on the right:

Turbulence Intensity normalisation of power curves and large-scale fleet validation Julien TISSOT, Alexander SCHLACHER, Christopher GRAY

	Replace	Replace
	with QR code	with QR code

i4SEE TECH GmbH, Austria

PO.246

Visit our homepage at i4see.com