Zhibo Luo

Yaoyao Zhao
email: yaoyao.zhao@mcgill.ca

Efficient thermal finite element modeling of selective laser melting of Inconel 718

Keywords: finite element method, adaptive mesh, selective laser melting, heat transfer, MPI

In the powder bed fusion process, an accurate prediction of the transient temperature field of a part is essential to calculate the subsequent thermal stress evolution and microstructure propagation in that part. The experimental method is time-consuming and expensive since the temperature field is controlled by many process parameters.

Numerical heat transfer models can be used to estimate the temperature field at any time point. However, traditional numerical simulation schemes are not suitable for the layerwised fabrication process due to the extremely high computational cost. The computational cost mainly relies on the element number and time step size. This research provides a new efficient and part-level simulation scheme based on an open-source finite element library, which is able to adaptively refine and coarsen the mesh and solve finite element equations with multiple processors in a parallel way. Here, a new mesh strategy that aims to reduce the element number while keeping the solution accuracy is developed.

The simulation speed is 12x -18x faster compared with the traditional simulation scheme depending on the scale of the simulated domain and number of processors. Simulation results have been compared with the experimental results of an Inconel 718 component.

It is shown that the testing point in the simulation experiences the same thermal cycles of the same point in the experiment. This simulation scheme can also be used to optimize the process parameters such as scanning pattern, scan velocity, and layer thickness and can be easily extended to other additive manufacturing processes.

Introduction

Selective laser melting (SLM), one of the widely used powder bed fusion (PBF) additive manufacturing (AM) process, is capable of selectively melting the powder material with prescribed geometry and forming functional metal part with complex structures layer-bylayer [START_REF] Debroy | Additive manufacturing of metallic components -Process, structure and properties[END_REF][START_REF] Technologies | Terminology, Standard terminology for additive manufacturing technologies[END_REF]. Due to the ability to fabricate quality parts, the SLM process has attracted significant attention in the academic field and has been applied extensively in industries such as health-care [START_REF] Averyanova | Manufacture of Co-Cr dental crowns and bridges by selective laser Melting technology: This paper presents the successful application of the selective laser melting technology in dental frameworks manufacturing from Co-Cr alloy using Phenix PM 100T Dental Machine over a production period of 14 months[END_REF][START_REF] Vandenbroucke | Selective laser melting of biocompatible metals for rapid manufacturing of medical parts[END_REF], aerospace [START_REF] Rochus | New applications of rapid prototyping and rapid manufacturing (RP/RM) technologies for space instrumentation[END_REF], and injection modeling [START_REF] Dalgarno | Manufacture of production injection mould tooling incorporating conformal cooling channels via indirect selective laser sintering[END_REF]. However, high thermal gradient induced issues, such as distortion, residual stress, pores, delamination, and fatigue fracture of the final part, preclude the wide adoption of the SLM process in industries [START_REF] Mani | Measurement Science Needs for Real-time Control of Additive Manufacturing Powder Bed Fusion Processes[END_REF]. These defects are mainly controlled by more than 100 process parameters and at least 20 of them lead to first-order influence on these defects, such as scanning pattern, speed, powder properties, and laser beam properties [START_REF] Yan | Data-driven multi-scale multi-physics models to derive process-structureproperty relationships for additive manufacturing[END_REF]. As the experimental trial-and-error method is time-consuming and costly, numerical simulation has been widely investigated to predict the in-process thermal and mechanical behaviors to optimize those parameters. Traditional simulation approach such as the finite element method (FEM) origins from the welding process with single or several tracks [START_REF] Denlinger | Residual stress and distortion modeling of electron beam direct manufacturing Ti-6Al-4V[END_REF][START_REF] Mackwood | Thermal modelling of laser welding and related processes: a literature review[END_REF].

However, the SLM process consists of thousands of scanning tracks that lead to unacceptable computational time. Therefore, scaling the simulation from microscale to macroscale is a promising way to reduce the scale of a simulation model [START_REF] Lavery | A review of Computational Modelling of Additive Layer Manufacturing-multi-scale and multi-physics[END_REF][START_REF] Li | On the Simulation Scalability of Predicting Residual Stress and Distortion in Selective Laser Melting[END_REF][START_REF] Papadakis | A computational reduction model for appraising structural effects in selective laser melting manufacturing: A methodical model reduction proposed for time-efficient finite element analysis of larger components in Selective Laser Melting[END_REF][START_REF] Markl | Multi-Scale Modeling of Powder-Bed-Based Additive Manufacturing[END_REF]. For example, Li. C et al. [START_REF] Li | On the Simulation Scalability of Predicting Residual Stress and Distortion in Selective Laser Melting[END_REF] developed a multiscale simulation model using ABAQUS. In this method, microscale thermal analysis of a point of a single scanning track was simulated to determine the heating rate, heating time, and melt pool depth. Then, these parameters were used to determine the equivalent body heat flux of a scanning track in the mesoscale hatch model. Finally, in the macroscale part model, a block-layer consisting of several layers was heated by using the equivalent heat flux. Similarly, Papadakis [START_REF] Papadakis | A computational reduction model for appraising structural effects in selective laser melting manufacturing: A methodical model reduction proposed for time-efficient finite element analysis of larger components in Selective Laser Melting[END_REF] developed a computational reduction model to accelerate the thermal-mechanical analysis of the SLM process. Heat inputs were used as thermal loads to estimate the temperature field and distortion by using SYSWELD in three different levels: a detailed moving heat source which could mimic the scanning path of the laser beam, a temperature profile of unit length that was defined on the scanning track, and a temperature profile representing the total heat input per unit volume that was defined on a whole layer or several layers. The simulation result indicates a good correlation between the reduced heat input approach and the moving heat input approach. This multiscale simulation method is faster compared with the detailed thermal-mechanical simulation model as the transient simulation on millions of scanning tracks has beyond the ability of a powerful computer even computer clusters [START_REF] Li | A multiscale modeling approach for fast prediction of part distortion in selective laser melting[END_REF]. However, mapping from the microscale model to the macroscale model will lose some key information such as heating and cooling rates, maximum temperature, melting and solidification of the melt pool, and re-melting of bottom layers. In addition, critical process parameters, such as scanning pattern, scan speed, and layer thickness, cannot be included in a macroscale model. During the modeling of a process-structure-property framework for AM processes [START_REF] Yan | Data-driven multi-scale multi-physics models to derive process-structureproperty relationships for additive manufacturing[END_REF][START_REF] Smith | Linking process, structure, property, and performance for metal-based additive manufacturing: computational approaches with experimental support[END_REF], the thermal distribution and evolution are the basis for other physical processes. Mechanical responses, such as thermal stress and distortion [START_REF] Hodge | Implementation of a thermomechanical model for the simulation of selective laser melting[END_REF], fatigue performance [START_REF] Song | Differences in microstructure and properties between selective laser melting and traditional manufacturing for fabrication of metal parts: A review[END_REF], material properties, and microstructure evolution [START_REF] Chen | A hybrid finite-element and cellular-automaton framework for modeling 3D microstructure of Ti-6Al-4V alloy during solid-solid phase transformation in additive manufacturing[END_REF], depending on the thermal history during the SLM process. Although the inherent strain method [START_REF] Keller | New method for fast predictions of residual stress and distortion of AM parts[END_REF], which is originally used in structural analysis of welding components [START_REF] Ueda | A predicting method of welding residual stress using source of residual stress (report I): characteristics of inherent strain (source of residual stress)(mechanics, strength & structural design)[END_REF], has been used to estimate the residual stress and distortion of final parts [START_REF] Afazov | Distortion prediction and compensation in selective laser melting[END_REF], it does not incorporate the accumulation of thermal stress during fabrication process and cannot include the influence of layer thickness, powder property, and scanning pattern into the quality of final parts. Therefore, considering the framework of process-structure-property-performance for AM, detailed prediction of thermal history is essential for subsequent analysis of thermal stress and microstructure evolutions.

Generally, the detailed finite element (FE) based thermal modeling method can be classified into three categories: heat conduction model, thermal-fluid flow model, and thermal-fluid-vaporization model [START_REF] Gan | Benchmark Study of Thermal Behavior, Surface Topography, and Dendritic Microstructure in Selective Laser Melting of Inconel 625[END_REF]. The heat conduction model considers the heat diffusion but neglects the fluid flow convection enhanced by the Marangoni effect and heat loss caused by vaporization inside the melt pool. This method is fast as it does not need to solve the fluid flow speed and other physical quantities within the melt pool, but the simulation accuracy is not as high as the other two models. The neglection of Marangoni driven flow and vaporization leads to an overestimation of the peak temperature and an underestimation of the solid cooling rate as well as an overestimation of the solidification cooling rate within the melt pool [START_REF] Gan | Benchmark Study of Thermal Behavior, Surface Topography, and Dendritic Microstructure in Selective Laser Melting of Inconel 625[END_REF][START_REF] Manvatkar | Heat transfer and material flow during laser assisted multi-layer additive manufacturing[END_REF]. The thermal-fluid flow model not only considers the heat conduction but also the liquid flow inside the melt pool driven by the Marangoni effect which increases the heat convection by the thermo-capillary flow. This method is more accurate but requires a higher computational cost. The thermal-fluidvaporization model further considers the heat loss caused by vaporization inside the melt pool. Among these three models, despite the thermal-fluid-vaporization method is the most accurate one, it is not suitable for part-level simulation due to more time steps. In addition, more element number is required and more physical quantities need to be solved during the whole simulation process. Therefore, a transient heat conduction model is used in this research for predicting the transient thermal fields and thermal history to reduce computational costs. This research has developed a new efficient and accurate FE model for predicting transient material thermal behavior during the SLM process in part-level. It mainly provides a simulation tool to comprehensively record the temperature field and history of each time step and each node for further mechanical analysis and microstructure evolution. In addition, this simulation tool can be used to optimize process parameters in a more practical way due to a user-defined scanning path file indicating the parameters for each scanning track in each layer. Similar simulation models for thermal analysis in SLM or other PBF processes exist [START_REF] Denlinger | Thermomechanical modeling of additive manufacturing large parts[END_REF][START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF][START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF][START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF][START_REF] Tran | Heat transfer simulations of selective laser melting process based on volumetric heat source with powder size consideration[END_REF][START_REF] Mozaffar | Acceleration strategies for explicit finite element analysis of metal powder-based additive manufacturing processes using graphical processing units[END_REF]. For example, Denlinger and Michaleris [START_REF] Denlinger | Thermomechanical modeling of additive manufacturing large parts[END_REF] developed a part-scale thermal-mechanical analysis using a quiet inactive element activation strategy on the platform of Pan computing [START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF]. A large part deposited for 107 layers was simulated and the simulated distortion has a 29% error compared with the experimental result. Zhang et al. [START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF] developed a 3D macroscopic FE transient thermal model using adaptive mesh technique to preserve a conformal mesh along with the interface between workpiece and powder. The temperature field of a complex Inconel 718 part is simulated for 17 hours using 60 processors. Bo Cheng et al. [START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF] developed a 3D multiple-layer FE model for the SLM using ABAQUS, where volumetric heat source with certain penetration was accounted for heat flow in the melt pool. This coupled thermalmechanical simulation of 64 mm ! domain took 92 hours. 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] developed a comprehensive heat transfer model for a single scanning track using the COMSOL Multiphysics module. A volumetric heat source considering the particle size was used and the simulated peak temperature among the melt pool was validated with experimental results. This model can be run on a common workstation to predict the melt pool size and the contact width between the melt pool and deposition region. However, the aforementioned models mainly use either commercial codes, such as Netfabb, 3DSim, GeonX, Ansys, Abaqus, Additiveworks, and Simufact additive, or existing FE packages that cannot be easily modified with desired process parameters.

The developed thermal FE model in this research is based on the open-source FE library named Deal.II [START_REF] Bangerth | deal. II-a general-purpose object-oriented finite element library[END_REF]. This research developed several methods to improve simulation efficiency while keeping solution accuracy. First, a new adaptive mesh scheme is developed to control the total active element involving in the simulation domain while accurately to capture the input heat flux. This adaptive mesh scheme can maintain dense mesh around the heat source and coarse mesh away from the heat source. It can also reduce the number of mesh refinements that take extra time to interpolate solutions between two different FE meshes. In addition, as the layer number increases, this mesh scheme guarantees the total element number does not increase tremendously. Another similar adaptive mesh scheme is developed by Mukherjee et al. [START_REF] Mukherjee | Heat and fluid flow in additive manufacturing-Part I: Modeling of powder bed fusion[END_REF]. This scheme changes the grid system every step with the same amount of elements to guarantee the dense grid region moved along with a laser spot. The grid size decreased ahead of the laser spot and increased at the back of the laser spot as the laser spot traveled through each scanning track. However, this mesh scheme is only suitable for a part with regular shape as the total element number is not constant for an irregular part domain, i.e., complex geometry parts. Second, a temperature-dependent material property file is designed and incorporated in the simulation. This material file specifies thermal properties of powder, liquid, and solid material at different temperatures so that the heat is transferred accurately in the simulated domain. The material file also includes the influence of latent heat on temperature fields using the apparent heat capacity method.

Third, the assembly of global matrices is optimized. And the setup of degree of freedom (DoF), the mesh refinement, assembly of the FE system, and solving of nonlinear equations are computed parallelly. These can significantly reduce the simulation time for each time step. Fourth, this research has also designed a scanning path file that includes the laser energy density, scanning patterns, speed, and time step for each layer. The scanning path file provides users more freedom to control process parameters so it can fit the real fabrication process as much as possible. For example, during a real fabrication process, the laser power may change for different layers. For each layer, the laser energy density is lower in boundary regions than in internal regions. Layer thicknesses are different in a real printing process due to different shrinkage rates. This may influence the heat transfer in deposited domains. These problems can be easily overcome by using the predefined scanning path file.

Transient thermal FE model for the SLM process

Transient thermal simulation of the SLM process in part-level requires an efficient simulation scheme that can balance the computational cost and solution accuracy. This research has developed an efficient simulation scheme that can dramatically reduce the simulation time with satisfactory accuracy. This is achieved through analyzing the fabrication process from a track to a layer and finally to a part, optimizing the mesh scheme strategy, modeling the material deposition and material phase transition mechanism, and improving the efficiency of solving nonlinear equations. This section presents the computational heat conduction model for estimating temperature fields/history during the whole SLM process in part-level.

Figure 1. Heat transfer process in the SLM process

The simplified heat transfer process for the SLM process is shown in Figure 1. The build domain ( Ω ), boundary conditions, and input heat flux are indicated. Specifically, a prescribed temperature is applied to the bottom substrate (𝛤 " ). Thermal insulation is applied to the lateral of the whole domain due to the weak heat conduction ability of powder material (𝛤 # ) far from the laser scanning area [START_REF] Dong | Finite element simulations of temperature distribution and of densification of a titanium powder during metal laser sintering[END_REF]. The heat loss caused by radiation and convection from the lateral surface is ignored in the FE domain in order to reduce the computational cost. Thermal radiation and convection conditions are applied to the top surface of the current deposited layer (𝛤 $ ). These boundaries satisfy this relation:

∂Ω = 𝛤 " ∪ 𝛤 # ∪ 𝛤 $ .

Assumptions

The following assumptions are made to guarantee an acceptable computational cost.

-Densities of powder, solid, and liquid material are assumed to be constant in each phase throughout the whole simulation process [START_REF] Gan | Benchmark Study of Thermal Behavior, Surface Topography, and Dendritic Microstructure in Selective Laser Melting of Inconel 625[END_REF] to guarantee the conservation of mass; -The top layer surface of the deposited domain is assumed to be flat throughout the whole simulation process [START_REF] Mukherjee | Heat and fluid flow in additive manufacturing-Part I: Modeling of powder bed fusion[END_REF];

-The fluid flow inside the melt pool due to the Marangoni effect is ignored [START_REF] Wolff | A framework to link localized cooling and properties of directed energy deposition (DED)-processed Ti-6Al-4V[END_REF];

-The heat loss and mass loss inside the melt pool due to vaporization are not considered in the calculations [START_REF] Gan | Benchmark Study of Thermal Behavior, Surface Topography, and Dendritic Microstructure in Selective Laser Melting of Inconel 625[END_REF][START_REF] Gan | Modeling of thermal behavior and mass transport in multi-layer laser additive manufacturing of Ni-based alloy on cast iron[END_REF]; -Powder material is considered as a continuous solid material and its thermal properties are the same before its melt point. After the melting point, it will be considered as liquid material;

-Thermal properties of powder and solid material are assumed to be temperaturedependent [START_REF] Gan | Numerical simulation of thermal behavior and multicomponent mass transfer in direct laser deposition of Co-base alloy on steel[END_REF];

-The input heat is assumed to be a Gaussian distribution [START_REF] Gan | Numerical simulation of thermal behavior and multicomponent mass transfer in direct laser deposition of Co-base alloy on steel[END_REF] and treated as a volume heat source using the double ellipsoid model [START_REF] Denlinger | Thermomechanical model development and in situ experimental validation of the Laser Powder-Bed Fusion process[END_REF];

Governing equations

During the SLM process, if the powder material is considered as continuous material and the coordinate system is fixed to the build plate, the transient heat transfer equation in Lagrangian reference frame can be written as

𝜌𝐶 % (𝑇) &$(𝒙,*) &* = ∇ • [𝑘(𝑇)∇𝑇(𝒙, 𝑡)] + 𝑄(𝒙, 𝑡), ∀𝒙 ∈ Ω, 𝑡 ∈ (0, 𝑡 ,-. ) (1)
where 𝑘(𝑇) and 𝐶 % (𝑇) are temperature-dependent thermal conductivity and specific heat capacity, respectively. 𝜌 is the density and 𝑄 is the internal volume heat source.

The low thermal conduction ability of powder material makes it possible to consider the boundary of powder bed as thermal insulation. The heat flux constitutive function is given as

𝑞(𝒙, 𝑡) = 𝑘(𝑇)∇𝑇(𝒙, 𝑡) (2)
Then, three boundaries are used in this research:

a. Dirichlet boundary with a given temperature at the bottom substrate

𝑇(𝒙, 𝑡) = 𝑇 / ∀𝒙 ∈ 𝛤 " , 𝑡 ∈ (0, 𝑡 ,-. ) (3)
b. Thermal insulation of lateral powder bed surface

𝑞(𝒙, 𝑡) = 0 ∀𝒙 ∈ 𝛤 # , 𝑡 ∈ (0, 𝑡 ,-. ) (4)
Equation ( 1) is a transient nonlinear problem. The first step is to discretize the equation in time, then discretize the domain in space through the finite element method. The 𝜃scheme time discretization method will be used to approximate the continuous solution 𝑇 at some time 𝑡 -with discretized solution 𝑇 -.

𝜌𝐶 % -$ ! 8$ !"# 9 ! -𝜃𝑘 -∇ • (∇𝑇 -) -(1 -𝜃)𝑘 -8: ∇ • (∇𝑇 -8: ) = 𝜃𝑄 -+ (1 -𝜃)𝑄 -8: (7)
where 𝜏 -= 𝑡 --𝑡 -8: is the time step size. 𝐶 % -, 𝑘 -, and 𝑄 -are temperature-dependent specific heat, thermal conductivity, and heat source at temperature 𝑇(𝒙, 𝑡 -), respectively.

𝑇 -is used to approximate 𝑇(𝒙, 𝑡 -) at time step 𝑡 -. Here, 𝜃 = 0.5 generates the Crank-Nicolson time discretization.

Given the 𝜃-scheme time discretization, the space discretization happens by multiplying equation [START_REF] Mani | Measurement Science Needs for Real-time Control of Additive Manufacturing Powder Bed Fusion Processes[END_REF] with test function 𝜑, integrating by part, and then applying boundaries. This yields the following FE system.

[𝑀 -+ 𝜏 -𝜃𝐾 -+ 𝜏 -𝜃𝑀 01-2 - + 𝜏 -𝜃𝑀 63. -]𝑇 - = 𝑀 -𝑇 -8: -𝜏 -(1 -𝜃)𝐾 -8: 𝑇 -8: + 𝜏 -𝜃𝐹 -+ 𝜏 -(1 -𝜃)𝐹 -8:
-

;
𝐹 -is the right-hand-side vector: 𝐹 ; -= (𝑄 -, 𝜑 ; ) = ;

𝑓 -is the heat flux boundary vector: 𝑓 ; -= (𝑞 @ABC -, 𝜑 ; ) > *+,-.

The first mass term in the right-hand side of equation ( 8) is exactly written as:

𝑀 -𝑇 -8: = 𝜌𝐶 % -∑ (𝜑 ; -, 𝜑 < -8: )𝑇 < -8: D !"# <EF , 𝑖 = 0,1,2, … 𝑁 - (9)
where 𝑁 -8: and 𝜑 -8: are the total DoFs and shape function at step 𝑛 -1, respectively.

When the FE meshes used in two sequential time steps are different, i.e., the number of DoFs in two sequential time steps are different ( 𝑁 -≠ 𝑁 -8: ), (𝜑 ; -, 𝜑 < -8: ) yields a rectangular mass matrix. It is difficult to compute these integrals because if we loop over the cells at time step 𝑛 and evaluate 𝜑 < -at quadrature points of these cells, we need to evaluate 𝜑 < -8: at the same quadrature points which may not necessarily correspond to the cells of the mesh at time step 𝑛 -1. What is worse, 𝜑 < -8: may not even be defined via these cells. These integrals can be realized by matching the cells in two different meshes at the same place and by integrating shape functions on two meshes. But this method costs too much time if the number of cells is large, so it is not practical to be implemented in codes. This research utilizes an interpolation operation that can interpolate the solution from the old mesh into the new mesh each time the mesh is changed. Then, equation ( 9) can be rewritten as

𝜌𝐶 % -∑ (𝜑 ; -, 𝜑 < -8: )𝑇 < -8: D !"# <EF = S𝜑 < -, 𝐼 G -𝑇 G -8: U , 𝑖 = 0, 1, 2, … 𝑁 - (10)
where 𝐼 G -is the interpolation operator onto the mesh field used in the 𝑛 th time step. -8: is the approximated solution.

𝑇 G -8: = ∑ 𝜑 < -8: D !"# F 𝑇 <
This method is easy to be implemented in the FE code for time-dependent adaptive meshes. Similarly, the stiffness matrix 𝐾 -can also use this method to coordinate the shape functions between two sequential time steps with different FE meshes.

Equation ( 8) is a set of fully discrete nonlinear equations. The first step is to linearize these equations. Then, the linearized equations can be solved using the Conjugate Gradient method. This nonlinear algebraic equation system is a special structure with the form

𝐴(𝑢)𝑢 = 𝑏(𝑢) (11)
where 𝐴(𝑢) is a (𝑁 + 1) × (𝑁 + 1) matrix function of 𝑢 and 𝑏 is a vector function of 𝑢: 𝑏 = (𝑏 F , … , 𝑏 D ) . Iteration method will be needed to calculate the unknown solution 𝑢 = (𝑢 F , … , 𝑢 D ).

In this research, a relaxed Picard iteration method will be adopted to linearize the above nonlinear equation system. For a given nonlinear equation 𝐴(𝑢)𝑢 = 𝑏(𝑢) and an initial guess 𝑢, iterate step 1 to step 3 until 𝑢 convergence:

• Step 1, solve the linear equation system for 𝑢 H , 𝐴(𝑢)𝑢 H = 𝑏(𝑢);

• Step 2, 𝑢 = 𝜔𝑢 H + (1 -𝜔)𝑢, 0 < 𝜔 < 1; • Step 3, update 𝑢 = 𝑢.
The solution will converge when the following equation is satisfied ‖𝑢 -𝑢‖ < 10 8? ~10 87 (12)

Heat source

The heat source model is critical to simulate the heat input in the SLM process. The powder material is continuously melted when the laser spot moves along a scanning track and solidified after the laser spot completely passes. In this circumstance, the heat source is a function of space and time, which is termed as moving Gaussian point heat source.

Depending on the penetration depth of the heat flux, the moving Gaussian point heat source can be modelled as surface-related where heat flux is only applied on top of the surface; and volume-related where heat flux is generated within a small volume [START_REF] Luo | A survey of finite element analysis of temperature and thermal stress fields in powder bed fusion Additive Manufacturing[END_REF]. The laser radiation can penetrate the powder bed into several micrometers and undergoes multiple relflections inside the powder bed [START_REF] Mukherjee | Heat and fluid flow in additive manufacturing-Part I: Modeling of powder bed fusion[END_REF], for example, it can penetrate maximum 20 𝑢𝑚 in nickel [START_REF] Ma | Temperature and stress analysis and simulation in fractal scanningbased laser sintering[END_REF] and 65 𝑢𝑚 in titanium [START_REF] Fischer | Sintering of commercially pure titanium powder with a Nd: YAG laser source[END_REF]. When the penetration depth is smaller than the cell size, the input heat flux can be modelled as surface-related Gaussian point heat source. Otherwise, the input heat flux can be modelled as volume-related Gaussian point heat source. Most of the research model the laser beam as surface-related heat source when low energy density is used [START_REF] Li | A multiscale modeling approach for fast prediction of part distortion in selective laser melting[END_REF][START_REF] Hodge | Implementation of a thermomechanical model for the simulation of selective laser melting[END_REF][START_REF] Dong | Finite element simulations of temperature distribution and of densification of a titanium powder during metal laser sintering[END_REF]. Other research model the laser beam as volume-related heat source when high energy density is used [START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF][START_REF] Cheng | Stress and deformation evaluations of scanning strategy effect in selective laser melting[END_REF][START_REF] Denlinger | Effect of stress relaxation on distortion in additive manufacturing process modeling[END_REF][START_REF] Bruna-Rosso | Selective laser melting finite element modeling: Validation with high-speed imaging and lack of fusion defects prediction[END_REF].

For the surface-related heat source, the heat flux can be applied in 𝑥, 𝑦 directions at the top surface of powder bed. Then, the moving Gaussian point heat source can be written as:

𝑞(𝑥, 𝑦, 𝑧, 𝑡)| IEJKLMJN PQ NPR STUQVWK = 2𝜂𝑃 𝜋𝜔 ? exp h- 2𝑥 ? + 2𝑦 ? 𝜔 ? i ( 13
)
where 𝜂 is the laser absorptivity of the powder material, 𝑃 is the laser power, 𝜔 is the maximum radius of the laser spot.

For the volume-related heat source, the heat flux can be applied in 𝑥, 𝑦, 𝑧 directions in the powder bed. Then, the moving Gaussian point heat source can be written as:

𝑞(𝑥, 𝑦, 𝑧, 𝑡) = 6√3𝜂𝑃 𝑎𝑏𝑐𝜋√𝜋 exp h- 3𝑥 ? 𝑎 ? - 3𝑦 ? 𝑏 ? - 3𝑧 ? 𝑐 ? i (14)
where 𝑎, 𝑏, 𝑐 presents the semi-axes of the ellipsoid in 𝑥, 𝑦, 𝑧 directions, respectively.

Equation ( 13) and ( 14) are the Gaussian point heat flux distribution in Cartesian coordinates. In the Lagrangian reference frame, the position of the laser spot is related to time, which can be written as:

p 𝑥 = 𝑥 F + 𝑣 C 𝑡 𝑦 = 𝑦 F + 𝑣 X 𝑡 𝑧 = ℎ I (𝑡) (15)
where (𝑥 F , 𝑦 F , 𝑧 F ) is the start point of a scanning track. 𝑣 C = 𝑣cos𝜃, 𝑣 X = 𝑣sin𝜃, and 𝜃 is the scanning orientation. ℎ I (𝑡) is the height of current deposited part and increases with the increase of layer number.

It should be noted that when equation ( 13) is used as a surface-related heat source, it is applied as a boundary condition upon the top deposition layer. When equation ( 14) is used as a volume-related heat source, it is applied as an internal volumetric heat source.

In this research, as the laser beam has a high energy density that can reflect several times inside the powder material, thus the internal volumetric heat source 𝑄 is modeled using the moving double ellipsoid model [START_REF] Goldak | A new finite element model for welding heat sources[END_REF].

Simulation scheme

Figure 2. Simulation scheme for simulating part-level SLM process

The AM process usually consists of hundreds of layers. During the building of each layer, the laser spot travels back and forth hundreds of times. Therefore, the basic simulation unit is a scanning track. This research has designed a data structure to drive the simulation process track-by-track and layer-by-layer. The systematic simulation scheme for the SLM process is shown in Figure 2. After obtaining the necessary input data such as the process parameters and scanning path, the transient thermal analysis starts. For the simulation of each layer, the thermal analysis is divided into two stages, the heating process and cooling process. During the heating process, the laser spot scans the prescribed layer geometry track-by-track. After all the scanning tracks of this layer are simulated, the powder bed is cooled down for dwell time, then a new layer of powder is spread upon the deposited FE domain. This process continues until the last layer is simulated. The data structure of the scanning path file is illustrated in Figure 3. This file can be generated from the scanning path by tracking the motion of the laser spot. The simulation process will follow the sequence of the layer number and scanning track number defined in this file. This research improves the simulation efficiency from a single scanning track to a single layer with multiple tracks and finally a part with multiple layers.

As it is known that the simulation time is directly determined by the total element number, the focus of this research is to reduce the element number while keeping the accuracy of the final solution.

Efficient modeling approaches

This section presents modeling approaches from a single scanning track to a single layer and finally to a full part. For simplicity and consistency, the three simulations are performed in the same domain, and the same process parameters and constant material properties are used.

Simulation of a single scanning track

The transient temperature evolution of a 6.4 mm scanning track on top of a substrate is simulated. The FE model and initial coarse mesh are shown in Figure 4. A linear hexahedra element that possesses 8 nodes is used for this simulation. The maximum cell size, which is 0.8 mm, is defined as 0-level refined. The minimum cell size, which is 0.05 mm, is defined as 4-levels refined. The refining of a hexahedra element generates 8 children hexahedra elements that are only half size of the parent element. Similarly, the coarsening of 8 adjacent elements forms a large element. The constant material properties of Inconel 718 and process parameters for the computational analysis are listed in Table 1.

Minimum element size and time step size

The diameter of the laser spot is extremely small compared with a part. The ratio of a part dimension to a laser diameter ranges from 50 to more than 1000. Therefore, the mesh should be very fine around the laser spot to capture the effective input heat flux and the high thermal gradient in the heat-affected zone. However, the number of elements increases enormously if the part or even just several layers have meshed with very fine elements. Melt pool temperature distributions after 2.05 ms in different refinement levels are shown in Figure 5(a). The maximum temperature is plotted against varying mesh refinement levels in Figure 5(b). Since the maximum temperature stays constant as the increased mesh density, the refined mesh yields a convergent temperature distribution [START_REF] Romano | Laser Additive Melting and Solidification of Inconel 718: Finite Element Simulation and Experiment[END_REF]. In this case, the laser diameter is set to 110 𝜇𝑚. Considering the solution accuracy

Maximum cell size: 0.8 mm Minimum cell size: 0.05 mm and FE system scale, 5 levels of mesh refinement, which has a mesh density of 64000 elements/mm 3 , is a good compromise in this research. The maximum temperature difference between the 5 levels refinement and 7 levels refinement is 2.4%. For this transient heat transfer problem, the heat source is moving with a certain velocity 𝑣. The heat source should not skip any element when it travels. Then the maximum time step is restrained by

∆𝑡 = . 2 (16)
where 𝑑 is the diameter of the laser spot.

Temperature-dependent material properties

To accurately predict the temperature distribution of the SLM processed part, temperature-dependent material properties should be incorporated in solving this heat transfer problem. Most researches treat the density and thermal conductivity as a function of porosity of the powder bed using the Maxwell model [START_REF] Cernuschi | Modelling of thermal conductivity of porous materials: application to thick thermal barrier coatings[END_REF]. Then the effective thermal conductivity and density of the powder bed are derived as [START_REF] Foteinopoulos | On thermal modeling of Additive Manufacturing processes[END_REF] 𝑘 ,@@ = Y .

:8F.[\@ , 0 ≤ 𝑓 < 1 (17)

𝜌 ,@@ = (1 -𝑓)𝜌 / (18)
where 𝑘 / is the thermal conductivity of solid material, 𝑓 is the volumetric fraction of porosity.

However, it is not convenient to illustrate the relationship between porosity and temperature during the building process. In this research, thermal conductivity, thermal convective coefficient, radiation coefficient, and surface emissivity are modelled as a function of temperature. Then, a material file is defined to tabulate the desired material properties with prescribed temperatures. For a given temperature 𝑇, the property value (𝑉) can be obtained as,

𝑉 % ] (𝑇) = { 𝑓 : (𝑇) 𝑇 ≤ 𝑇 : 𝑓 ? (𝑇) 𝑇 : < 𝑇 < 𝑇 ? 𝑓 ! (𝑇) … (19)
where 𝑝 indicates powder, solid, or liquid material type. 𝑞 indicates thermal conductivity, specific heat, convection coefficient, radiation coefficient, or surface emissivity.

𝑓 : , 𝑓 ? , and 𝑓 ! present the function at different temperature ranges, respectively.

Material phase transition mechanism

The powder material's property is fairly different from the solid material at the same temperature. Therefore, a mechanism should be established to indicate the material phase of an element at a certain temperature. In this research, the material undergoes three different phases from powder to liquid and finally to solid type. During the melting and solidification process, phase change induced latent heat is modelled based on the apparent heat capacity method [START_REF] Hu | Mathematical modelling of solidification and melting: a review[END_REF]. When material phase changes from solid to liquid, the specific heat is artificially increased around the liquidus to encompass the absorbed latent heat during the melting process. This artificially increased heat capacity is called apparent heat capacity. The temperature range between solidus and liquidus is referred to "mushy zone". The apparent heat capacity can be a function of temperature. The latent heat is equal to the integration of the apparent heat capacity over the temperature in the mushy zone.

Based on the apparent heat capacity, the effective specific heat capacity with respect to the temperature is modelled 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].

𝐶 % ,@@ (𝑇) = { 𝐶 % (𝑇) 𝑇 < 𝑇 / 𝐶 % + ∆𝐻 # /(𝑇 A -𝑇 / ) 𝑇 / < 𝑇 < 𝑇 A 𝐶 % (𝑇) 𝑇 A < 𝑇 (20)
where 𝐶 % is the specific heat capacity, 𝑇 / and 𝑇 A are the solidus temperature and liquidus temperature, respectively. ∆𝐻 # is the enthalpy change during the phase change. The enthalpy of the solid or powder is defined as

𝐻 # = ∫ 𝜌𝐶 % 𝑑𝑇 (21)
The material phase transition rule is illustrated in Figure 6. As this is a simplified model, powder splashing, Marangoni effect, and fluid flow are not considered.

Analysis of computational cost for static mesh scheme and adaptive mesh scheme

This single scanning track is simulated using a static mesh scheme and adaptive mesh scheme, respectively. For the adaptive mesh scheme, the criterion for refining and coarsening is based on the Kelly error estimation [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: [START_REF] Afazov | Distortion prediction and compensation in selective laser melting[END_REF] where 𝑎 @ 𝜕𝑢 G /𝜕𝒏 denotes the solution gradient's jump at cell faces, 𝑓 is the cell face, 𝐶 @ is the scale of the cell size and usually can be taken as the greatest length of the diagonals of cell faces.

Powder Liquid Solid T powder > T melt T solid > T melt T liquid < T melt 𝑒 ? = ∑ 𝐶 @ ∫ ƒ 3 * &B / &𝒏 " ? 𝑑𝑜 &^* @∈&^
By using this refinement criterion, the mesh is only refined in the region where the temperature gradient is relatively large and is coarsened in the region where the temperature gradient is relatively small. As the temperature gradient around the melt pool is always large, this refinement criterion is based on the location of the laser spot. Table 2 summarizes the computational cost of these simulations. Each simulation is completed on the Intel Xeon® Gold 5122 CPU @3.6GHz with 64GB RAM. With the increase of element number, the computational cost increases as well. The adaptive mesh scheme can greatly reduce the element number and thus reduce the simulation time. As the laser spot moves with prescribed scan speed, the mesh should be adaptively refined and coarsened for each time step. However, extra efforts are needed to refine the mesh, renumber the nodes, and interpolate shape functions between two different meshes.

Therefore, refinement frequency cannot be very high. In terms of the build domain dimension and solution accuracy, adaptively refining the mesh every 2~5 steps is a better choice compared with the static mesh scheme. Then, the simulation time of the adaptive mesh scheme can be reduced to 44%~49% of the simulation time of the static mesh scheme. The simulated temperature distribution of this scanning track after 2.91 ms (step-51) using two mesh schemes is plotted in Figure 7. The temperature difference between the adaptive mesh refined every 2 steps and the static mesh is around 2.1% using the following temperature distribution error for the current time step. It can be concluded from Table 2 and Figure 7 that adaptively refining of the elements close to the laser spot every 2 steps is a much more efficient mesh scheme with enough accuracy.

Optimization of the assembly of global matrix

To solve the FE system, the matrices on the left-hand side and vectors on the right-hand side of equation ( 8) need to be assembled for each time step. As shown in Table 3, the assembly of matrix and vector occupies 50%~70% of wall CPU time in each time step.

Besides, each time step needs 10~20 iterations to solve the nonlinear algebraic equations. In each time step, as the FE mesh stays the same in each iteration, the global matrices only need to be assembled once if the thermal properties are constant. However, when temperature-dependent material properties are used, thermal properties of each node are changed as the temperature of each node varies in different iterations. Therefore, the assembly of each local element and distribution of local elements' contribution to the global matrix are needed in each iteration. This consumes a large amount of computational time as each time step may need to iterate 10~20 times to get a converged solution. To reduce the matrix assembling work, a reduced-assembly matrix is developed.

The basic idea comes from the fact that the temperature difference of an element far from the laser spot between two sequential time steps is very small and can be ignored. In this circumstance, the local contribution of these elements far from the laser spot is always the same between two sequential time steps. In other words, for each time step, the full global matrix only needs to be assembled once in the first iteration. In the following iterations, only the local contribution of the elements whose temperature is quite different from the previous iteration will be reassembled and their local contributions will be distributed to the global matrix. By doing this, the time cost for assembling the full global matrix can be greatly reduced. For the adaptive mesh used in this research, the solution of the hanging node is obtained by averaging the solutions of neighboring nodes [START_REF] Luo | Numerical simulation of part-level temperature fields during selective laser melting of stainless steel 316L[END_REF]. For example, if node 5 in Figure 8(b) is a hanging node, then 𝑇 \ = 1/2(𝑇 ? + 𝑇 a ). The laser spot only influences a very small region between two sequential time steps, shown in Figure 8(a). For elements that are outside of the heat-affected region, the temperature variation can be ignored. The local and global numbering rules of a node that determines the distribution of local matrix entries to the global matrix entries is presented in Figure 8(b). In the matrix assembly process of these four elements, shown in Figure 8 where 𝐴 ;< , = (∇𝜑 ; , , ∇𝜑 < , ) = 3 is the local stiffness contribution between node 𝑖 and node 𝑗 from element 𝑒. 𝑘 , is the temperature-dependent thermal conductivity of element 𝑒. If the average temperature of element 2 is significantly changed and the average temperatures of elements 1, 3, and 4 stay the same compared with the previous iteration, only the gray regions in the global stiffness matrix need to be updated. Therefore, the computational cost for assembling the global matrix can be significantly reduced when the scale of the FE mesh becomes very large. Table 4 indicates that by adopting this matrix assembly strategy, the assembling time can be decreased to 55.3% of the time from the regular assembly strategy for a single time step.

MPI method to further reduce assembling workload

Up to now, only one processor is involved in the calculation of solving the heat transfer problem. With the increase of element number, this is not an efficient way to refine the mesh, assemble matrix, solve nonlinear equations, and so on. The mesh can be segmented into several pieces and each piece can be calculated by one processor.

Therefore, multiple processors can be used to parallelly compute the nonlinear algebraic equations. This research uses the Message Passing Interface (MPI) to solve the transient heat transfer problem in a parallel way. Deal.II library provides interfaces for PETSc/Trilinos, which is a commonly used package for parallel linear algebra, and package METIS can be used to partition meshes parallelly and efficiently.

The computational time of simulating the 6.4 mm scanning track with static mesh and adaptive mesh scheme is summarized in Table 5. For static mesh, the more processors involved in the simulation process, the faster the simulation process is. For adaptive mesh, the simulation time does not linearly decrease with the increase of processors numbers.

It can be noted that the static mesh scheme is faster than the adaptive mesh scheme when more processors are involved. This is attributed to the fact that for the refinement process, extra time is needed for communication and synchronization between these processors to form a complete mesh. In addition, since constant material properties are used for static mesh schemes, the total FE system only needs to be assembled once at the beginning. While for adaptive mesh scheme with temperature-dependent material properties, the FE system needs to be reassembled and solved for each time step.

Considering the real printing process, the adaptive mesh scheme is much more efficient than the static mesh scheme. In this case, the simulation efficiency can be improved by maximum 5.5 times for static mesh scheme with constant material properties; The simulation efficiency of the adaptive mesh scheme can be improved by 1.9~2.5 times depending on the number of CPUs and the refinement frequency. As the scale of the FE system becomes large, more processors lead to a much faster simulation because the communication time becomes less important compared to the assembly, setup updates, and process solving.

Simulation of a single layer with multiple scanning tracks

Simulation of a single layer repeats the previous simulation of one scanning track for hundreds of times, which is rather time-consuming. The previous adaptive mesh scheme is not an efficient way for simulating multiple scanning tracks due to tremendous mesh refinement operations. This section focuses on optimizing the adaptive mesh scheme based on the location of the scanning track.

Optimization of adaptive mesh scheme

As the laser spot moves, the mesh close to the heat source should be dynamically refined while the mesh far from the heat source should be coarsened. Though the mesh refinement and coarsening can well control the element number, this operation is timeconsuming and leads to the reassembly of the nonlinear FE system. For the transient thermal simulation of one layer or multiple layers, this research has developed another adaptive mesh scheme to compromise the element number and heat input. Elements close to the current scanning track rather than around the laser spot are refined. This may lead to a little increase in element number for each step but can greatly reduce the matrix assembling and mesh refinement operations. In addition, the solution accuracy can be guaranteed as the fine mesh can always be found around the laser spot. The comparison of time costs between these two adaptive mesh schemes on simulating a single layer with 58 scanning tracks using a constant zigzag scanning pattern is shown in Table 6. Eight processors are used and the simulation time is significantly reduced by using the scanning track based adaptive mesh scheme. After several numerical experiments, it has been proved that the simulation efficiency can be improved by 1.83 times compared with the laser spot based adaptive mesh scheme. The temperature distribution and material phase transition of a layer are presented in Figure 9, where the 30 th scanning track is selected as a representative. The mesh is clearly refined based on the current scanning track. then the maximum coarsen level is 4 -log ? 𝑛 A3X,6 . In summary, for a certain layer number 𝑛 A3X,6 , it can be written as

𝐻 %36* = 𝑛 A3X,6 • 𝑇 *G;0Y = 2 Y • 𝑇 *G;0Y , (𝑘 = 0,1,2, … ) (24)
𝑛 A3X,6 = 2 F + 2 : + 2 ? … + 2 * , 𝑡 = 0,1,2, … (25)
Solve the above equation for 𝑡, then the maximum coarsen level is

𝑛 0136/,- 43C = ' 4 -𝑡, 𝑡 < 4 0, 𝑡 ≥ 4 (26)
For example, as shown in Figure 10

Simulation of a part with multiple layers

Most existing FE models can only simulate the SLM process for a single layer or several layers. This is mainly hindered by the huge element number due to inappropriate mesh scheme and material deposition strategy. This section presents a new material deposition strategy to simulate the powder spread process and develops the criterion for the mesh refinement scheme between layers. In addition, as the activation of a new layer of elements renumbers the DoFs, a solution transfer strategy is developed to transfer solutions between the old mesh and new mesh.

Material deposition strategy

To model the powder spreading process for finite element analysis, metal particles are assumed to act as solid-like material with powder material properties. The geometry model is sliced into an FE model layer-by-layer at the beginning and only the elements located among deposited layers are involved in the heat transfer analysis. For a layer that

is not yet deposited, no DoF will be assigned to its elements. The mechanism of material deposition and material phase transition during the simulation is shown in Figure 11. A modified inactive element method is used in this research. The basic idea of the inactive element method is to mesh the whole continuous domain into FE meshes at the beginning, then only deposited elements are assigned with DoFs while non-deposited elements have no DoFs. In addition, the material phase of the new active element is a powder type. Then the powder element is converted to liquid or solid type as the laser spot goes in and out.

Mesh refinement and coarsening scheme between layers

Figure 12. Adaptive mesh refinement and coarsening scheme between layers

After cooling for dwell time, a new layer of material will be activated and DoFs will be assigned to newly activated elements. When a new layer of elements is being activated, the size of the elements on the top layer should not exceed layer thickness. Otherwise, the deposited surface becomes uneven, which influences the following mesh refinement for that layer. Figure 12 illustrates the adaptive mesh refinement scheme when a new layer of elements is activated upon the deposited FE domain. First, elements close to the top layer are refined until the maximum element size reaches the layer thickness. Then, a new layer of elements is activated by using newly assigned DoFs with initial powder material type. After that, the mesh is refined based on the location of scanning tracks. table, an initial temperature condition will be assigned to those nodes. The solution transfer between layer 1 and layer 2 after the FE domain are cooled for 0.01 second is illustrated in Figure 14. In the simulation, newly activated elements are assigned with powder material type and initial room temperature. Then, a full part is simulated using the aforementioned techniques and algorithms to estimate the temperature history and temperature fields during the SLM process of Inconel 718. nodes. The coarse mesh is created using Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand postprocessing facilities[END_REF]. The refined mesh of the first layer is shown in Figure 15(c). This mesh has 50,000 hexahedral elements and 58,484 nodes.

The maximum element size is 0.

Material properties and process parameters

The temperature-dependent thermal properties of Inconel 718 are listed in Table 7 [START_REF] Romano | Laser Additive Melting and Solidification of Inconel 718: Finite Element Simulation and Experiment[END_REF][START_REF] Ahn | Thermal Characteristics in the Cutting of Inconel 718 Superalloy Using CW Nd: YAG Laser[END_REF]. The process parameters for the simulation are listed in Table 8. It has been proved that a variation of 50% of convection coefficient value only results in less than 1% difference of final results [START_REF] Dunbar | Experimental validation of finite element modeling for laser powder bed fusion deformation[END_REF]. Therefore, a constant convection coefficient is applied to the top layer.

Adaptive mesh control and material deposition

FE mesh in different tracks and layers is shown in Figure 16 based on the adaptive mesh scheme in track level and layer level under a constant zigzag scanning pattern. Figure and the mesh at the top deposited layers is refined. Therefore, the scale of the FE system does not increase significantly, the computational cost is decreased, and the solution accuracy is guaranteed.

Validation

In this section, experimental in-situ measurements of temperature history [START_REF] Dunbar | Development of experimental method for in situ distortion and temperature measurements during the laser powder bed fusion additive manufacturing process[END_REF] will be used to validate the proposed simulation scheme for modelling the SLM process in part level.

Test point (0, 0, 0)

The Inconel 718 powder is used and its temperature-dependent material properties are listed in Table 7. The process parameters are listed in Table 8. The experimental setup is shown in Figure 19(a). A K-type thermocouple was used to capture the temperature history of the testing point at (12, 0, 0) mm of the thin metal plate. Four supporting screws were used to hold the thin metal plate to the bottom substrate in order to allow space for inserting thermocouples. The printed part is shown in Figure 19

Comparison of temperature results

The comparison of the experimental and simulated temperature history of testing points under different scanning patterns is shown in Figure 20. In case 1, the scanning direction rotated 67° for each layer. In case 2, a constant zigzag scanning pattern was applied for each layer. As the probe of the thermocouple could be damaged by the high temperature of the melt pool, it was attached far from the deposited area and could only record the temperature history of the thin plate.

As shown in Figure 20, the in-situ recorded temperatures increased with the increase of printing layers and finally reached a plateau at around 380 seconds when the last layer was printed. After that, the temperatures rapidly went down to room temperature. The temperature difference between case 1 and case 2 was very small at the beginning and became more perceivable as the building progressed. The simulated results followed the same trend. However, the simulated temperature history at the beginning did not fit well with the experimental results. This can be attributed to the simplified boundary conditions and the physical thermal process of the FE model. The neglection of the Marangoni effect may lead to an overestimation of the peak temperature among the melt pool, which could delay the thermal heat for conduction. After 60 seconds, the temperature differences between the simulation and experiment results were about 4 K and 6 K for case 1 and case 2, respectively. The simulated results overpredicted the experimental results by an acceptable 5% at the testing points. This can also be attributed to the neglection of the Marangoni effect, which could increase the heat convection by thermocapillary flow inside the melt pool. The cooling rates were almost the same for each case after the parts were completely built. According to Wolff et al. [START_REF] Wolff | Experimentally validated predictions of thermal history and microhardness in laser-deposited Inconel 718 on carbon steel[END_REF], the cooling rate should be significantly underestimated for the heat conduction model. This can be attributed to the fact that the Marangoni effect significantly influenced the melt pool while the testing point was fixed and located far from the melt pool. Therefore, this FE heat conduction model can accurately predict the thermal fields and history except for the regions close to the melt pool. For both scanning patterns, the differences in scanning times were very small, but perceivable temperature differences were observed. The temperature from the constant scanning pattern was lower than the corresponding temperature from the rotating scanning pattern. In summary, the thermal cycles from the numerical simulation are comparable to those from experimental results. The proposed simulation scheme for estimating transient temperature fields is able to predict the trend of the temperature evolution during the SLM process. number has a large influence on the computational cost, the proposed simulation scheme utilizes an adaptive mesh strategy based on the input scanning path file in order to reduce the scale of the FE system and the computational cost.

To clarify the efficiency of the proposed simulation model, the simulation is analyzed from a single scanning track to a layer with multiple scanning tracks, and finally to a full part with multiple layers. The simulation time can be reduced by 44%~49% using the proposed adaptive mesh scheme compared with the traditional static mesh scheme for simulating a single scanning track. In addition, the matrix assembly time can be reduced by 55.3% using the optimized matrix assembly method compared with the regular matrix assembly method for a single time step. The MPI method can further reduce the computational time by 1.9~2.5 times depending on the number of CPUs used and the refinement frequency.

Finally, in the simulation of multiple scanning tracks and multiple layers, an optimized adaptive mesh scheme based on the location of scanning tracks is developed and the simulation efficiency can be improved by 1.83 times compared with the laser spot based adaptive mesh scheme. In conclusion, the overall simulation efficiency can be improved by 12.6~18.6 times using the algorithms and techniques proposed in this research compared with the traditional static mesh scheme using a single processor. The ΔT case1 = 4 K ΔT case2 = 6 K experimental results show that the temperature differences can be controlled within 5% by using this simulation scheme compared with the static mesh scheme. In summary, this research proposes the following methods to improve the simulation efficiency:

• Design a scanning path file which contains the necessary process parameters for each track and each layer, and use it to drive the simulation process track-by-track;

• Develop a new material deposition model which can reduce the active element number, mimic the powder spreading process, and incorporate material phase transition among elements;

• Develop a new adaptive mesh scheme which can automatically refine and coarsen elements based on the current scanning track, and reduce the time for mesh change;

• Develop a new temperature-dependent material property distribution strategy which designates the properties of solid, liquid, and powder material at different temperatures;

• Optimize the global matrix assembly by distributing the changed local element contributions and updating the corresponding global matrix regions;

• Utilize a parallel computing method that can segment the domain into several subdomains and assign each subdomain to a processor. All the subdomains can be processed simultaneously to reduce the computational time.

Numerical results show that the melt pool will influence the neighboring scanning tracks for at least three layers, so the layer thickness is a key process parameter in optimizing the mechanical response of the final part. The scanning pattern has a large influence on the temperature distribution. The temperature difference between the simulation and experimental results varies 4~6 K by using constant zigzag scanning pattern and rotating scanning pattern.

This research provides an efficient simulation tool for modeling thermal response in the SLM process from a single track level to a single layer level, and to a part level. As thermal stresses are mainly caused by temperature gradients, the obtained temperature fields can be used to identify the possible concentrated thermal regions in order to avoid high thermal stress and deformation in the real fabrication process. The estimated temperature history of each point can also be used to feed the evolution of microstructure. In addition, this simulation tool can be used for other AM processes such as electron beam melting and directed energy deposition process with little modification.

Figure 3 .

 3 Figure 3. The data structure of a scanning path file

Figure 4 .

 4 Figure 4. Initial FE model for simulation analysis

Figure 5 .

 5 Figure 5. Comparison of temperature distribution in different mesh density

Figure 6 .

 6 Figure 6. Material phase transition rules

  𝑥 ; , 𝑦 ; , 𝑧 ; ) ∈ 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑠𝑐𝑎𝑛 𝑡𝑟𝑎𝑐𝑘[START_REF] Gan | Benchmark Study of Thermal Behavior, Surface Topography, and Dendritic Microstructure in Selective Laser Melting of Inconel 625[END_REF] where 𝑇 /*3*;0 ; and 𝑇 3.3%*;2, ; are the temperatures of node-𝑖 belonging to the current scanning track in static mesh and adaptive mesh, respectively. 𝑙 is the length of the scanning track. ∆𝑙 is the cell size.

Figure 7 .

 7 Figure 7. Temperature distributions of a single scanning track in two different mesh schemes

(a) 5

 5 levels refinement static mesh (b) Adaptively refined mesh every 2 time step

Figure 8 .

 8 Figure 8. Optimization of matrix assemble

Figure 9 .

 9 Figure 9. Temperature distribution and material phase transition of 30 th scanning track using scanning track based adaptive mesh scheme 3.2.2 Scanning track based adaptive mesh scheme in simulating a layer The melt pool generated by high energy laser beam influences neighboring scanning tracks and beneath layers [49]. Therefore, a bounding box is defined as the region for mesh refinement. The width of the bounding box should be at least three times large of the laser diameter and the depth should cover three layers of elements. As shown in Figure 10(a), fine elements can be found around the melt pool shown in a red elliptical shape. Elements outside of the bounding box are attempted to be coarsened. To keep a flat deposition surface, the maximum coarsening level in the top layer is varied for different layers. Figure 10(b) illustrates that if elements in the top layer are over coarsened, the deposition surface will become uneven. Then, in the following mesh refinement operations, the number of active elements will be changed and the mesh refinement process, solution interpolation between two meshes, element activation, and matrix assembly will become complicated. Therefore, a rule is devised to calculate the maximum coarsening level of the top layer. The minimum and maximum element size are set to 0.05 mm and 0.8 mm, respectively. The minimum size is exactly the layer thickness (0.05 mm). Then a coarse element with a size of 0.8 mm needs to be refined 4 times to reach the minimum element size. If the part height 𝐻 %36* (layer number times layer thickness) is 2 Y times of the layer thickness (𝑇 *G;0Y ),

  (a) Temperature distribution of 30th scanning track (b) Material phase transition of 30th scanning track

is 4 .

 4 (c)~(f), top elements in layer 1 have a size of 0.05 mm (layer thickness) and 𝑛 0136/,- 43C When simulating layer 2, 𝑛 0136/,- 43C is 3. Similarly, according to the above equation (26), as the simulation process proceeds to simulate layer 3 and layer 4, the maximum coarsen level of top elements becomes 4 and 2, respectively.

Figure 10 .

 10 Figure 10. Scanning track based adaptive mesh scheme in simulating a layer

  elements belongs to the top of layer 3: refined 4 levels substrate

Figure 11 .

 11 Figure 11. The mechanism of element activation and material phase transitionTo successfully activate a new layer of elements, the hp finite element method is used to indicate the element type of each cell in the FE domain. The hp finite element method was initially proposed by Babuska and Guo in the 1980s[START_REF] Guo | The hp version of the finite element method[END_REF] as an alternative to reduce the approximation error of the shape function solutions by either decreasing the cell size h or increasing the polynomial degree p. In this research, a combination of two different finite elements is used to represent the active domain or inactive domain. As shown in Figure11, a special FE container (FECollection) including the regular finite element with a certain polynomial degree (FE_Q) and the special finite element with zero degree of freedom (FE_Nothing) is created. Then this FEColletion is distributed to each cell based on the indicator of the deposition domain (indicator 0) or void domain (indicator 1). In the

Figure 13 .

 13 Figure 13. Node solution storage and transfer between two sequential layers The total element number only changes when a mesh refinement and coarsening operation is applied. When a new layer of elements is activated, the total element number stays the same and the total DoF increases. According to the node numbering rules, the global node number changes as well when a new layer is activated. This is explained in Figure 13 (a) in detail. For each element, the mapping relation between local node number and global node number changes. A new mapping relation should be rebuilt to link the local node number and the global node number. This will lead to another problem: How to map the existing node solution to the new DoFs. As the coordinates of each vertex does not change when a new layer is activated, it is possible to establish a mapping rule between the node coordinates and the node solution. In this research, each node solution is mapped to its corresponding vertex coordinates. When a new layer of elements is activated and a new node numbering is rebuilt, the current node solution will try to find the previously saved solution through the mapping function.Figure 13(b) illustrates the mapping relation between node solutions and node coordinates. The nodal old solution is saved together with corresponding node coordinates before assigning DoFs to the newly activated elements. After that, try to transfer the old nodal solution to the new nodal solution through the mapping relation. For those nodes that are not found in the mapping

Figure 14 . 4 . 1

 1441 Figure 14. Node solution transfer between layer 1 and layer 2

8

 8 mm and the minimum element size is 0.05 mm. Once the elements of the first layer are activated and refined, elements close to the first scanning track are refined for an additional time. After the first track is simulated, the elements close to the second track are refined and the elements located in the first track are coarsened. This process continues until all scanning tracks in the current layer are simulated. The total simulation time for this domain is 22.2 hours by using 8 cores of Intel Xeon® Gold 5122 CPU @3.6GHz with 64 GB RAM.

Figure 15 . Part geometry and initial coarse mesh 4 . 2

 1542 Figure 15. Part geometry and initial coarse mesh 4.2 Numerical simulation

  (a) part dimension (b) initial coarse mesh (c) refined mesh of first layer

  16(a)(b) and Figure 16(c)(d) show the FE mesh in the simulation of different scanning tracks in layer 22 and layer 23, respectively. The total active element number of layer 22 is 41,761 with 44,687 DoFs. The total active element number of layer 23 is 75,809 with 73,101 DoFs. According to equation (26), the minimum refinement level of layer 22 is 3 and the maximum element size at the top of layer 22 is 0.1 mm. When the last scanning track is finished, these elements need to be refined for one more time in order to activate the elements in layer 23 and guarantee a flat surface. Then, the minimum refinement level of layer 23 is 4 and the maximum element size at the top of layer 23 is 0.05 mm.

Figure 16 . Adaptive mesh based on track level and layer level 4 . 3 results 4 . 3 . 1

 1643431 Figure 16. Adaptive mesh based on track level and layer level4.3 Simulation results4.3.1 Temperature fields of a layerThe temperature fields and material phase transitions during the printing process is shown in Figure17. The temperature field of layer 18 cooled down for 12 seconds after the printing is shown in Figure17(a). The temperature fields of layer 19 and layer 21 are shown in Figure 17(b) and (c), respectively. In Figure 17(c), the melt pool is formed

Figure 17 . 4 . 3 . 2

 17432 Figure 17. Temperature fields and material transitions during the SLM process4.3.2 Temperature history of a testing pointThe temperature history of the point located at (0, 0, 0) mm, which is the center of the top surface of the substrate, is shown in Figure18. When the laser spot approaches the point, the melt pool can influence the neighboring 2~3 scanning tracks. The thermal heat

Figure 18 .

 18 Figure 18. Temperature history of a point at the top of the substrate

  (b) and Figure 19(c).

Figure 19 .

 19 Figure 19. Experiment design and the printed part [54]

  Built green part (c) Post-build substrate attached to the baseplate through support set scew

Figure 20 .

 20 Figure 20. Comparison of in-situ measured temperature history (solid curves) and simulated temperature history (dashed curves) under rotating and constant scanning pattern during the SLM process

Table 1 . Constant thermal properties and process parameters for the thermal analysis Parameter Value

 1

	Absorptivity of laser spot, 𝜂	0.6
	Layer thickness	50 𝜇𝑚
	Scan speed, 𝑣	960 mm/s

Ambient Temperature 300 K Conductivity, 𝑘 12 W/m K Surface convection coefficient, ℎ 25 W/m 2 K Density, 𝜌 8146 kg/m 3 Laser power, 𝑃 180 W Diameter of laser, 𝑑 0.11 mm 𝑎, 𝑏, and 𝑐 of ellipsoid heat source model 110, 55, and 55 𝜇𝑚

Table 2 . Comparison of computational cost using static mesh scheme and adaptive mesh scheme for a single scanning track

 2

	Mesh scheme	Element number	Node number	Wall CPU time (s)
	Static mesh -5 level refinement	152,896	220,472	577
	Adaptive mesh -every 1 step			353
	Adaptive mesh -every 2 step	40,756 ~ 67,916	59,106 ~ 93,088	284
	Adaptive mesh -every 5 step			253

Table 3 . Wall CPU time-consuming percentage for consecutive 60 time steps

 3

	Matrix/Vector assembly	Mesh refinement and coarsening	Solve linear FE system	Output solution	Setup system	other
	50%~70%	5%~20%	5%~10%	10%~20%	5%~10% 5%~15%

Table 4 . Comparison of computational cost for a single time step using regular and optimized matrix assembly methods

 4

		Number of elements	Number of DoFs	Wall CPU time
	Regular matrix assembly method Improved matrix assembly method	40,756	59,106	1295 ms 717 ms

Table 5 . Comparison of the simulation time using multiple processors

 5

	Processor		Wall CPU time (second)
	number	static mesh	refine every 1 step	refine every 2 steps
	1	577	353	284
	2	327	270	193
	4	175	204	137
	6	130	198	124
	8	105	187	115

Table 6 . Comparison of single layer simulation time between laser spot based and scanning track based adaptive mesh schemes Adaptive mesh scheme Wall CPU time (s) Efficiency improved

 6

	laser spot based	6939	-
	scanning track based	2450	183%

Table 7 . Temperature-dependent thermal properties of Inconel 718

 7

	Temperatur e (K)	Thermal conductivity (W/mK) Powder Solid	Specific heat(J/K kg)	Enthalpy (KJ/m 3 )	Emissivity
	298	1.96	8.9	435	0	0.539
	373	2.35	10.8	455	2.69e5	0.533
	473	2.78	12.9	479	6.49e5	0.533
	573	3.24	15.2	497	1.04e6	0.534
	673	3.69	17.4	515	1.45e6	0.534
	773	3.96	18.7	427	1.86e6	0.535
	873	4.38	20.8	558	2.27e6	0.535
	973	4.61	21.9	568	2.72e6	0.536
	1073	5.62	26.9	680	3.19e6	0.536
	1173	5.40	25.8	640	3.71e6	0.537
	1273	5.58	26.7	620	4.18e6	0.537
	1373	5.91	28.3	640	4.66e6	0.538
	1443	6.11	29.3	650	4.98e6	0.538
	1609	6.17	444	720	7.22e6	0.329
	1673	6.17	444	720	7.43e6	0.332
	1773	6.17	444	720	7.86e6	0.337
	1873	6.17	444	720	8.28e6	0.341

Table 8 . Process parameters and constant material properties for the simulation

 8

	Process parameters	Value
	Laser power, 𝑃	280 W
	Laser radius, 𝜔	55 𝜇𝑚
	Thermal absorptivity, 𝜂	0.4
	Scanning speed, 𝑣	960 mm/s
	Scanning pattern	constant/rotate 67° zigzag
	Dwell time after each layer	12 s
	Hatch spacing	110 µm
	Layer thickness	50 µm
	Ambient temperature	300 K
	Material	Inconel 718 powder
	Convection coefficient, ℎ	25 W/m 2 /K [37]
	Density	8146 kg/m 3
	Solidus temperature	1533 K
	Liquidus temperature	1609 K

Acknowledgments: This research work is supported by the National Sciences and Engineering Research Council of Canada (NSERC) Strategic Network for Holistic Innovation In Additive Manufacturing (HI-AM) with NSERC Project Number: NETGP 494158-16.