Finite Element Thermal Analysis of Metal Parts Additively Manufactured via Selective Laser Melting

Abstract

In this chapter, a three-dimensional finite element model is developed to simulate the thermal behavior of the molten pool in selective laser melting (SLM) process. Laser-based additive manufacturing (AM) is a near net shape manufacturing process able to produce 3D objects. They are layer-wise built through selective melting of a metal powder bed. The necessary energy is provided by a laser source. The interaction between laser and material occurs within a few microseconds, hence the transient thermal behavior must be taken into account. A calibration procedure is carried out to fit the numerical solution with the experimental data. Once the calibration has corrected the thermal parameters, a dynamic mesh refinement is applied to reduce the computational cost. The scanning strategy adopted by the laser is simulated by a path simulator built using MatLab®, while numerical analysis is carried out using ANSYS®, a commercial finite element software. To improve the performance of the simulation, the two codes interact each other to solve the analysis. Temperature distribution and geometrical feature of the molten pool under different process conditions are investigated. Results from the FE analysis provide guidance for setting up the optimization of process parameters and develop a base for further residual stress analysis.

1. Introduction

Additive manufacturing is a 3D manufacturing process able to produce prototypes directly from a CAD file. Powder bed fusion (PBF) processes were, among others, the first commercialized AM processes. Nowadays, the most important powder bed fusion process is selective laser melting (SLM). Its basic method of operation is schematically shown in Figure 1 and, generally, all other powder bed fusion processes follow the same basic approach.

Figure 1.

Schematic drawing of SLM process.

The process fuses thin layers (typically 30÷60 μm thick) of metallic powder, which has been spread across the build area using a powder depositor and a leveler (wiper). The building area is enclosed in a chamber filled with inert gas to minimize oxidation and degradation of the powdered material. The base plate is heated in order to maintain the powder and the fabricated component at high temperature. Sometimes, infrared heaters are placed above the build platform to increase the temperature around the part being formed. The temperature within the chamber must be controlled to minimize the laser power requirements of the process (when preheating, less laser energy is required for fusion) and to prevent warping of the part due to non-uniform thermal expansion (curling).

Once the layer has been deposited and adjusted by the wiper, a laser beam is directed onto the powder bed and is moved using galvanometers in such a way that it thermally fuses the material to form the slice cross-section. Surrounding powder remains loose and serves as a support for subsequent layer deposition. After completing a layer, the build platform (base plate) is lowered and a new layer is laid. The beam scans the subsequent slice cross-section and the process repeats until the complete part is built. A cool-down period is typically required to allow the parts to uniformly achieve an adequately low temperature to be handled and exposed to ambient atmosphere. Finally, the parts are removed from the powder bed, loose powder is cleaned off the parts, and further finishing operations, if necessary, are performed.

Since the introduction of PBF, each new technology developer has introduced competing terminology to describe the mechanism by which fusion occurs, with variants of sintering and melting being the most popular. However, the use of a single word to describe the powder fusion mechanism is problematic as multiple mechanisms are possible. There are four different fusion mechanisms that are present in PBF processes. These include solid-state sintering, chemically-induced binding, liquid-phase sintering, and full melting [1, 2].

Since the attention of this chapter is mainly focused on SLM, the only full melting mechanism will be considered in detail. Full melting is indeed the mechanism most commonly associated with PBF processing of engineering metal alloys. In SLM, the entire region of material subjected to the heat source is melted to a depth exceeding the layer thickness. The thermal energy of subsequent scans of a laser is typically sufficient to remelt a portion of the previously solidified material. As a consequence, this type of melting is very effective at creating well-bonded, high-density structures.

Use of optimum process parameters is extremely important for producing satisfactory parts using PBF processes. Among them, special attention must be paid to the scanning strategy, viz. the movements applied by the laser on the powder surface. The path followed by the laser greatly influences the surface heat distribution, as it is responsible for highly localized temperature peaks. Scanning often occurs in two modes, contour and fill mode, as shown in Figure 2.

Figure 2.

An example of scanning strategy: fill mode and contour mode.

In contour mode, the outline of the part cross-section is scanned. This is typically done for surface finishing reasons around the perimeter. The rest of the cross-section is then scanned using a raster technique, whereby one axis is incrementally moved across the part being formed and the other axis is continuously swept back and forth. The contour can be scanned either before or after the cross-section, depending on the surface characteristics of the part being done. Multiple scanning strategies are available for the fill section as shown in Figure 3.

Figure 3.

The four different scanning strategies: meander, stripes, chess, and contour.

Since high thermal gradients result in large residual stresses, a great effort must be made to well define the movement of the heat source on the surface. The choice of the path is not unique and it must be done trying to find the most suitable one depending on the cross-section characteristics. Moreover, the scan strategy rotates at each subsequent layer as it is explained in Figure 4.

Figure 4.

The meander path rotation and the oriented path progress.

The rotation helps to balance the temperature of the working area. Despite the rotation, the path increment in horizontal direction remains constant and opposite to the gas flux. Consequently, the working area is kept clean because the melting slags drop far from the loose powder. The user can set the rotation angle. The angle used in this study has been suggested by the machine manufacturer and it is equal to 67°, so the rotation scheme repeats every 180 layers.

Nowadays, SLM specifically permits to manufacture highly customized products that are almost ready to use rather than mere prototypes. An object with very complex shape, which is almost impossible to produce with traditional technologies, can be easily created saving cost and time. Despite the benefit, these manufacturing processes are very problematic to control because of the high number of involved parameters. The numerical simulation helps to reduce the number of experimental trial-and-error tests necessary to optimize the process, minimizing the time and cost of manufacture of the final product while maintaining its quality unmodified. The thermal behavior of the molten pool is one of the most critical factors influencing the reliability of the part, as it affects geometrical accuracy, material properties, and residual stresses. In addition, the heat exchange between laser and powder is a very complex process that involves a lot of variables. Despite direct measurements of the thermal field are available, this is not enough to fully understand the molten pool behavior. Finite element analysis (FEA) is a powerful tool to gather more information about the process [3].

The first studies about numerical simulation applied to PBF were developed during the 1990s of the last century. In 1998, Williams and Deckard [4] were the first that setup a framework for the numerical analysis applied to PBF. The numerical model was able to study the effects of process parameters on the selective laser sintering (SLS). Afterward, Shiomi and Yoshidome proposed a model where melting and solidifying process could be studied, using FE analysis [5]. In this chapter, the behavior of a molten pool due to a single laser spot is considered. The laser works in pulse mode and the melted part of the powder are assumed to change into a sphere that increases the dimensions after each pulsed irradiation. Fisher and Romano [6] investigated the variation in sintering process using pulse and continuous heat source. They noticed that during continuous wave interaction, the grains are homogeneously heated (the authors use an analogy to distinguish between different heating: they say that with homogeneously heating, the grains are cooked), while during pulsed laser interaction, the heating is no longer homogeneous (the grains are roasted). As mentioned before, the interaction of laser radiation with powder bed is a very complex mechanism greatly responsible for the mechanical characteristics of the final product. Gusarov et al. made a great effort in developing a numerical model that is able to simulate the heat transfer between laser and powder [7, 8]. It is worth mentioning that all the new models developed for the SLM simulation always refer to the models that have just been developed for the welding process [9] inasmuch SLM can be considered as a series of micro welding processes. Among all the numerical frameworks proposed in the literature, it seems that the most powerful tool for the process simulation is the finite element method. The first example of 3D finite element analysis is the paper of Contuzzi and Campanelli [10]. The aim of this study is to evaluate the temperature evolution in a 3D part being formed with SLM process. The molten liquid phase is considered by introducing the specific latent heat while the thermal properties are kept constant, aiming to reduce the computational cost of the analysis. The model proposed by Kolossov and Boillat [11] instead, allows for the nonlinear behavior of thermal conductivity and specific heat due to the temperature change and phase transformation. Introducing nonlinear behavior for the material increases the reliability of the simulation. However, the measurements of the thermal properties are greatly affected by uncertainty. The thermal properties can be directly measured [12, 13, 14] or calculated with thermal models [15]. A more precise thermal conductivity of the powder bed, named effective thermal conductivity, can be calculated using the relation proposed by Dai and Shaw [16]. Their model encompasses the effect of temperature- and porosity-dependent thermal conduction and radiation as well as the temperature-dependent natural convection. Not only material properties but also heat source definitions greatly affect the numerical results. Hashemzadeh and Chen [17] collected a variety of heat distributions that can be used to model the heat exchange. Articles [11] and [18] are an example where the heat source can be easily modeled as a 2D flux applied on the powder surface, neglecting the effect of the laser penetration into the powder. Li and Wang [19] explained how a 3D heat flux can be modeled in order to improve the numerical results related to the molten pool depth. Despite the great accuracy of all the previous models, results coming from the numerical simulations are slightly different from the experimental evidence. This is due to the high complexity of the physics involved in SLM and indeed to a large number of simplifications needed to simulate the process. Hu and Kovacevic [18] tried to reduce the mismatch between numerical and experimental data with a calibration procedure that adjusts the process parameters in order to fit the molten pool dimension retrieved from the simulation with the experimental measurements. It can be noticed from the previous papers that the discretization of the domain is always kept constant so that all the analysis must be applied to microscopic scale in order to keep the computational time low. Patil, Pal et al. [20, 21] explained how the dynamic mesh refinement can be applied to reduce the computational cost needed to solve such problems (like SLM), where different levels of mesh density are required to capture the localized phenomena.

Despite the strong reduction in computational time due to the dynamic mesh refinement, it is not feasible to simulate the thermal field evolution of an entire object, also for very small components. This is due to the extremely large number of spots that are required to melt the cross-section of a single layer. This problem can be partially solved with an analytical solution of the thermal field [22, 23]. Nevertheless, this solution returns the global temperature of a single layer and it is not possible to distinguish between different heat distributions due to the multiple scanning strategies. Therefore, the simulation presented in this chapter is mainly focused on the study of the thermal evolution of a microscale domain. The FEA domain presented here includes only one layer of powder and its dimensions are chosen to reduce as much as possible the number of elements. A model that is able to simulate the thermal evolution of a complete part, taking also into consideration the effect of the scanning strategy, is still an open challenge.

SLM technology is widely used in different domains of the industry, such as aerospace, automotive, and consumer goods. However, the most important industrial application is the medical and surgical field. In this context, SLM is acting a major transformation of the traditional production techniques, more and more surgical implants are fabricated with PBF technologies. Since the industrial applications of SLM are continuously increasing, numerical simulations of the process become fundamental to predict the mechanical properties of the parts starting from the behavior of the thermal field. ANSYS® and ABAQUS® are an example of general purpose software that is widely used for the FE simulation of SLM process. Moreover, the industries are so interested in the new manufacturing possibilities offered by the SLM, so that numerical simulations were no more confined to academic research but became a powerful tool for the companies willing to develop new commercial products. As a result, specific software has launched onto the marketplace, for example, 3DSIM (based on dynamic adaptive mesh refinement [24]), SIMUFACT (based on the inherent strain theory [25, 26]), and NETFABB.

The framework presented in this chapter involves both ANSYS® and MatLab® to perform the FE analysis of a SLM process on a titanium alloy powder. At the beginning, all the process parameters needed for the simulation are collected. Then, the laser path is simulated using MatLab®. This permits the definition of each point of heat application. Their spatial coordinates are imported into ANSYS® to apply the thermal load into the FE model. The performances of the simulation are improved through a tight interaction between MatLab® and ANSYS®. The main analysis environment is ruled by MatLab® and ANSYS®, which is launched in batch mode only when it is needed. Moreover, a dynamic mesh refinement is applied to the model to reduce the computational cost of the simulation. Finally, a calibration procedure is proposed to correct the titanium alloy properties.

2. Process parameters

Thermal residual stresses are greatly affected by the temperature distribution. Moreover, the thermal field is greatly affected by laser characteristics and scanning strategies. Therefore, a proper setup of machines parameters is a key issue in reducing thermal residual stresses. Even if different scanning strategies can be adopted, the meandering path only will be considered in the following as it seems to be the most reliable one, ensuring a uniform thermal field and a low machining time. In this configuration, the laser scans the powder with the trajectory shown in Figure 2. The most important process parameters that are related to the different laser scanning strategies are listed below:

• Laser power [150÷500 W]

• Hatch distance [50÷105 μm]

• Point distance [20÷75 μm]

• Path rotation [0÷90°]

• Layer thickness [30÷60 μm]

• Time exposure [30÷70 μs]

• Powder absorbance [0.3÷0.5]

Changing these parameters can lead to different material characteristics, for example, material density, surface roughness, and porosity. All these characteristics are directly responsible for the mechanical properties of the sample. It is important, indeed, to have a powerful tool that is able to predict the thermal field with respect to different laser parameters.

The SLM machine available for this study is a Renishaw^® model AM250. The main characteristic of this machine is the pulsed laser technology. The process parameters used in the simulation are collected in Table 1 and refers to the machine setup.

Process parameters
Laser power (W)	160
Hatch distance (μm)	50
Point distance (μm)	20
Layer thickness (μm)	60
Path rotation (°)	67
Time exposure (μm)	30
Powder absorbance	0.5

Table 1.

Process parameters used in the simulation.

2.1. Material properties

The most suitable metallic material for surgical implants is Ti-6Al-4 V and the data collected in the subsequent tables always refer to this alloy. Ti-6Al-4 V is a titanium alloy with 6% of aluminum and 4% of Vanadium. It has excellent biocompatibility properties, as well as good mechanical properties [12].

Since the melting and cooling process is governed by nonlinear phenomena, the material properties used in the simulation must be temperature dependent. The parameters needed for the FE thermal analysis are thermal conductivity, density, and specific heat. In addition, for the processes involving phase changes, enthalpy is requested as well to account for the latent heat of the material. Specifically, the enthalpy (H) is related to density (ρ), specific heat (c), and temperature (T) according to Eq. (1):

$\mathrm{H}=\int \mathrm{\rho }\mathit{c}\left(\mathrm{T}\right)\mathrm{dT}$ E1

The thermal properties of Ti-6Al-4 V can be easily found in the technical literature [12, 13]. However, these properties are defined only for the solid bulk state and are not suitable to represent the material behavior in powder form or above the melting point. Regarding the properties of the liquid phase, the following simplifying assumptions are taken:

• The apparent powder density is assumed to be 60% of that of the bulk material [27, 28].

• The powder’s thermal conductivity might be calculated from the modified Zehner-Schlunder’s equation [15] or from the Dai and Shaw’s model [16]. Nevertheless, constants involved in the equation are given in the literature with a high uncertainty; therefore, accurate results are not straightforward. Despite the presence of those analytical models, thermal conductivity is simply supposed to be one order of magnitude less than the corresponding parameter for bulk material. In general, the powder is a worse heat conductor than the solid.

• The numerical model requires the definition of thermal parameters for the whole temperature domain, even though thermal conductivity has no real meaning in liquid and vapor phase. Starting from the consideration that thermal gradients are hindered in the liquid phase (and even more in vapor) since crystal lattices are randomly oriented, a first trial conductivity can be half of the value at the transition point. Moreover, the thermal conductivity is supposed to remain constant through liquid and vapor phase. The lack of reliability caused by these assumptions will be reduced by the calibration procedure.

• A distinction between powder and bulk thermal properties holds only in the solid-phase. Above the melting point, the thermal behavior is the same both for powder and bulk.

As a result, thermal conductivity below melting point is experimentally measured while for high-temperature simplifications apply. Data are shown in Table 2.

Thermal conductivity (W/mK)
Temperature (°C)	Solid	Powder
0	6.1874	0.6187
19	6.5660	0.6566
399	12.2620	1.2262
800	18.1490	1.8149
1198	23.4590	2.3459
1499	27.8010	2.7801
1605	29.9134	2.9913
1660	14.9567	14.9567
3265	14.9567	14.9567
3295	7.4784	7.4784
5000	7.4784	7.4784

Table 2.

Variation of thermal conductivity with temperature.

While the conductivity value is taken from experimental data (or from literature), the enthalpy is calculated from Eq. (1). Variables needed for the equation are collected in Table 3.

Thermal parameters	Solid	Powder
Density (kg/m^3)	4220	2532
Solidus temp (°C)	1605	1605
Liquidus temp (°C)	1660	1660
Vapor temp (°C)	3265	3265
Saturation temp (°C)	3295	3295
Specific heat for solid (J/kgK)	708.8	708.8
Specific heat for liquid (J/kgK)	1000	1000
Specific heat for vapor (J/kgK)	1500	1500
Latent heat of fusion (J/kg)	365,000	365,000
Latent heat of vapor (J/kg)	9,376,200	9,376,200

Table 3.

Thermal parameters.

The integration domain from Eq. (1) is divided into steps, from reference temperature (0°C) to limit temperature (5000°C). Each step corresponds to a different alloy phase: solid, liquid, and vapor. In order to reduce the computational cost, the density is kept constant throughout the temperature domain. The specific heat changes as a function of the temperature, but it is considered constant at each integration step. As a consequence, the enthalpy behavior is linear and can be easily calculated from the equations below.

${\mathrm{C}}_{\mathrm{avg}}=\frac{{\mathrm{C}}_{\mathrm{s}}+{\mathrm{C}}_{\mathrm{L}}}{2}\left(\text{average specific heat}\right)$ E2

${\mathrm{C}}^{\ast }={\mathrm{C}}_{\mathrm{avg}}+\left(\frac{\mathrm{L}}{{\mathrm{T}}_{\mathrm{L}}-{\mathrm{T}}_{\mathrm{S}}}\right)\left(\text{specific heat for transition}\right)$ E3

${\mathrm{H}}_{\mathrm{S}}=\mathit{\rho }{\mathrm{C}}_{\mathrm{S}}\left({\mathrm{T}}_{\mathrm{S}}-{\mathrm{T}}_0\right)\left(\text{enthalpy}\mathrm{at}\text{solid temperature}\right)$ E4

${\mathrm{H}}_{\mathrm{L}}={\mathrm{H}}_{\mathrm{S}}+\mathit{\rho }{\mathrm{C}}^{\ast }\left({\mathrm{T}}_{\mathrm{L}}-{\mathrm{T}}_{\mathrm{S}}\right)\left(\text{enthalpy}\mathrm{at}\text{liquid temperature}\right)$ E5

${\mathrm{H}}_{+}={\mathrm{H}}_{\mathrm{L}}+\mathit{\rho }{\mathrm{C}}_{\mathrm{L}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{L}}\right)\left(\text{enthalpy above liquid temperature}\right)$ E6

where ${\mathrm{C}}_{\mathrm{L}}$ : specific heat of solid; $\mathrm{\rho }$ : density; T_S: solidus temperature; T_L: Liquidus temperature; T₀: reference temperature; T: saturation temperature; and L: latent heat.

Table 4 shows the values of enthalpy calculated with the previous equations. The column named equation refers to previously numbered equations used for the calculation.

Enthalpy (J/m^3)
Temperature (°C)	Solid	Powder	Equation
0	0	0
1605	4.8008e + 9	2.8805e + 9	Eq. 3
1660	6.5398e + 9	6.5398e + 9	Eq. 4
3265	1.3312e + 10	1.3312e + 10	Eq. 5
3295	1.4210e + 10	1.4210e + 10
5000	6.5237e + 10	6.5237e + 10

Table 4.

Variation of enthalpy with temperature.

2.2. Heat source

The heat source can be modeled both as a surface or a volumetric thermal load [17]. In order to keep the computational cost as low as possible, the heat source is considered as a 2D heat flux applied on the surface of the powder bed. The thermal load transferred by the laser is called laser irradiance and can be represented as a Gaussian distribution [18]:

$\mathrm{I}=\frac{2\mathrm{AP}}{\mathrm{\pi }{\mathrm{r}}_{\max }^2}\bullet {\mathrm{e}}^{\left(-2\frac{{\mathrm{r}}^2}{{\mathrm{r}}_{\max }^2}\right)}$ E7

where A: powder absorbance; P: laser powder; and ${\mathrm{r}}_{\max }$ : laser beam radius is defined as the radius in which the power density is reduced from the peak value by a factor of ${\mathrm{e}}^2$ [29].

Figure 5 shows irradiance values of the laser heat source. The laser radius is 35 μm.

Figure 5.

Laser beam irradiance with a Gaussian distribution.

3. Path simulation

A good knowledge of the heat source movements is fundamental to develop a reliable simulation that is able to predict the temperature distribution into the working area. The simulation is developed using MatLab^® and takes a CAD file storing the geometry of the part in STereo Lithography interface format (STL) as input. This file extension describes volumes through raw unstructured triangulated surfaces (tessellation). For each triangle, the unit normal and vertices (ordered by the right-hand rule) are collected using a three-dimensional Cartesian coordinate system. Figure 6 shows an example of tessellation.

Figure 6.

An example of CAD model tessellation.

Even if the simulation can analyze any kind of CAD geometry, a simple square block was chosen for the sake of simplicity. Since a non-uniform heat distribution is greatly responsible for thermal residual stresses, different scanning strategies can be used to reduce temperature peaks on the powder bed. Among them, the most suitable paths are: meandering, stripes, chess, and contour (Figure 3). The strategy chosen for the simulator is the meandering path, since it is a good trade-off between high deposition rate and low-temperature gradients. Moreover, for further reduction in local temperature peaks, the laser path changes orientation at each layer. Path simulation carried out with MatLab^® takes into account all these characteristics. At the beginning, a slicing algorithm is formulated to slice the 3D-geometry into a given number of layers as requested by the real process. Once having the spatial location and orientation of the triangles (from STL files), it is possible to intersect their triangular surfaces with horizontal planes (layers) using geometrical properties. Figure 7 explains the working principle of the slicing algorithm.

Figure 7.

A 2D example of a triangles tessellation.

Instead of a three-dimensional example, the triangles are sketched in 2D-space. The dashed line represents the slicing plane. Vertices are grouped with respect to Z-coordinates as listed in Table 5.

Vertices below the slicing plane	Vertices above the slicing plane
1 → A	9 → A
2 → E	10 → B
4 → B	12 → E
7 → C	14 → D
8 → D	15 → C

Table 5.

Vertices are grouped with respect to the slicing plane.

The triangles having the Z-coordinate higher than the height of the slicing plane belong to the upper vertices group, while the remaining belongs to the other group. The triangles that share vertices in both groups are those involved in the slicing procedure. Contour related to each cross-section is calculated from the intersection between the plane and the triangle surfaces as it is shown in Figure 8. Coordinates of the intersection point are called Z_p and X_p.

Figure 8.

The intersection points are determined through geometrical relations.

As a result, being the contour related to all cross-sections, the meandering path is applied to each layer. Since the path rotates at each layer, meandering slope changes and also the intersection between the path and the contour. Instead of a straightforward rotation of the path, it seems preferable to apply first a geometrical transformation to the contour and then apply the path, keeping its slope horizontal. This indeed allows a simple evaluation of the intersection points. Figure 9 shows an example that can help to understand the procedure.

Figure 9.

The geometrical transformation applied to cross-section for two different path rotations.

Referring to Figure 9, a square contour (blue) must be filled with a meandering contour tilted first with a slope of 45° and then 67°. Instead of tilting the path, the contour is rotated (cyan) by an angle corresponding to the desired slope of the path. Then, the contour is stretched so that the meander path is held constant and the hatch distance is kept unitary. Notice that the number of hatch spaces has been previously calculated, and keeping a unitary hatch distance will result in a more reliable intersection procedure.

The result of the slicing procedure is illustrated in Figure 10.

Figure 10.

Results from the path simulator.

As shown in Figure 10, the scanning strategy is defined for each layer with a sequential rotation (multiple of 67°). The dots represent the laser spot locations and the line used to connect them is indicated to emphasize the meandering path. The coordinates of the laser spot are stored into a file and can be used by ANSYS^® for the heat source application. Notice that distances between laser spots in Figure 10 are not in true scale. Meandering paths are shown just for illustrating how the laser moves on the powder surface.

4. Modeling approach

In SLM, the energy needed to melt the powder bed is provided by a laser source. The portion of material under the laser is heated because of the interaction between electromagnetic waves and powder grains. This type of heat transfer occurs in a very short time interval (microseconds) and provokes material modifications due to both phase (liquid and solid) and aggregation (powder and bulk) state changes. When the laser heats the surface, powder grains undergo very rapid heating that melts the material in the localized region surrounding the irradiated spot. After that, the laser is moved forth and the molten pool starts cooling and solidifying. At the end, the material has changed its aggregation state from powder to bulk. Since the path meanders on the surface, the material undergoes multiple reheating processes, sometimes above the melting point.

All previous characteristics would lead to a very complex and cumbersome simulation unless some simplifications are applied to the numerical model. The simulation is formulated taking into consideration all the SLM features, even if they are applied in a simpler way. For example, with the aim of reproducing both phase (solid-liquid-vapor) and aggregation state (bulk-powder) transformations, only material properties are defined, rather than complex thermodynamic models. Applications involving phase change can be approached using ANSYS^® through elements with enthalpy property capabilities.

The thermal transient analysis is necessary to take into account the high heating and cooling rate. Moreover, since the laser works in pulsed mode, the analysis is fully solved for each application point. The iterative algorithm forces the analysis to be solved, deleted, and restarted at each step. Consequently, nodal results must be continuously saved and uploaded through a mapping procedure as will be extensively explained later.

4.1. Numerical model

The governing heat transfer equation can be written as:

$-\left(\frac{\partial {\mathit{q}}_{\mathrm{x}}}{\partial \mathit{x}}+\frac{\partial {\mathit{q}}_{\mathrm{y}}}{\partial \mathit{y}}+\frac{\partial {\mathit{q}}_{\mathrm{z}}}{\partial \mathit{z}}\right)+\mathrm{Q}=\mathit{\rho c}\frac{\partial \mathit{T}}{\partial \mathit{t}}$ E8

where $q_x$ , $q_y$ and $q_z$ are components of heat flow through unit area. According to Fourier’s law:

${\mathrm{q}}_{\mathrm{x}}=-{\mathrm{k}}_{\mathrm{x}}\frac{\partial \mathit{T}}{\partial \mathit{x}}$

${\mathrm{q}}_{\mathrm{y}}=-{\mathrm{k}}_{\mathrm{y}}\frac{\partial \mathit{T}}{\partial \mathit{y}}$

${\mathrm{q}}_{\mathrm{z}}=-{\mathrm{k}}_{\mathrm{z}}\frac{\partial \mathit{T}}{\partial \mathit{z}}$

Notice that for isotropic, material thermal conductivity is the same: ${\mathrm{k}}_{\mathrm{x}}={\mathrm{k}}_{\mathrm{y}}={\mathrm{k}}_{\mathrm{z}}$ .

General formulation of governing differential equation can be obtained substituting Fourier’s law component in Eq. (8):

$\text{div}\left(\left[\mathrm{k}\right]\bullet \nabla \mathit{T}\right)+\mathrm{Q}=\mathit{\rho c}\frac{\partial \mathit{T}}{\partial \mathit{t}}$ E9

where $\left[\mathrm{k}\right]=\mathrm{k}\left[\mathrm{I}\right]$ is the thermal conductivity matrix.

Eq. (9) can be decomposed into a weak formulation [30] as the following system of differential equation of first order in $t$ :

$\left[{\mathrm{C}}_{\text{thermal}}\right]\left\{\dot{\mathrm{T}}\left(\mathrm{t}\right)\right\}+\left[{\mathrm{k}}_{\text{thermal}}\right]\left\{\mathrm{T}\left(\mathrm{t}\right)\right\}=\left\{\mathrm{F}\left(\mathrm{t}\right)\right\}\mathrm{t}\in \left\{0,{\mathrm{t}}^{\prime }\right\}$ E10

A description of the three matrices in Eq. (10) is given below:

1. The thermal stiffness matrix $\left[k_{\text{thermal}}\right]$ is expressed as follows:

   $\left[{\mathrm{k}}_{\text{thermal}}\right]={\int }_{\mathrm{V}}{\left[{\mathrm{B}}_{\text{thermal}}\right]}^{\mathrm{T}}\left[\mathrm{k}\right]\left[{\mathrm{B}}_{\text{thermal}}\right]\mathrm{dV}+{\int }_{\mathrm{S}}\mathrm{h}{\left[\mathrm{N}\right]}^{\mathrm{T}}\mathrm{Nds}$ E11

   where $\left[B_{\mathit{\text{thermal}}}\right]$ is the matrix containing the first derivatives of shape functions. The size of this matrix related to a brick element comprised of eight integration points is 3 × 8. Once computed, $\left[k_{\mathit{\text{thermal}}}\right]$ has dimension 8 × 8. $\mathit{dV}$ denotes the volume of the element. The surface integral is valid when the bulk is exposed to convection boundary conditions (i.e., boundary condition with imposed flux). $h$ is the convective heat transfer coefficient which has been assumed to be $12.5w/m^{2}K$ for Argon [20].

2. The thermal specific heat matrix $\left[C_{\text{thermal}}\right]$ is expressed as follows:

   $\left[{\mathrm{C}}_{\text{thermal}}\right]={\int }_{\mathrm{V}}{\left[\mathrm{N}\right]}^{\mathrm{T}}\left[\mathrm{N}\right]\mathit{\rho c}\mathit{dV}$ E12

   where $\left[\mathrm{N}\right]=\left[{\mathrm{N}}_1{\mathrm{N}}_2{\mathrm{N}}_3\dots {\mathrm{N}}_8\right]$ are the three-dimensional nodal shape functions of size 1 × 8. Once computed, $\left[C_{\mathit{\text{thermal}}}\right]$ has dimension 8 × 8. $\rho$ is the mass density and $c$ is the specific heat.

3. The thermal flux vector $\left\{F\left(t\right)\right\}$ is expressed as follow:

$\left\{\mathrm{F}\left(\mathrm{t}\right)\right\}={\int }_{\mathrm{S}}\left(\mathbf{q}\left(\mathbf{x},\mathbf{t}\right)\bullet \widehat{\mathbf{n}}\right)\bullet {\left[\mathrm{N}\right]}^{\mathrm{T}}\mathrm{dS}+{\int }_{\mathrm{S}}\mathrm{h}{\left[\mathrm{N}\right]}^{\mathrm{T}}{\mathrm{T}}_{\mathrm{f}}\mathrm{dS}$ E13

where q is the input heat flux depending on boundary conditions and $\mathbf{x}={\left\{\mathrm{x},\mathrm{y},\mathrm{z}\right\}}^{\mathrm{T}}$ is the position vector. $\mathit{dS}$ denotes the surface area of the element. The second surface integral in Eq. (13) is valid only when the convection boundary conditions apply. $T_f$ is the temperature of the environment into the working chamber.

Simulation reliability is subjected to the accuracy of the numerical model. Although geometrical features try to reproduce as much as possible the real SLM environment, simplifications must be taken regarding the boundary and loading conditions in order to reduce the computational cost. The attention is mainly focused on the melting process of one single powder layer, even if the algorithm allows for the simulation of multiple layers. Adding more layers, the number of elements grows as well as the number of spots; hence, the computational cost dramatically increases. As it will be presented in Section 5, dynamic mesh refinement is a powerful way to reduce the number of elements and, therefore, the time needed for the numerical solution, while preserving results accuracy.

The FE model is based on SLM real workspace. The geometry, as it is shown in Figure 11, is divided into two regular and constant meshes. The first one is at the surface and it is finer because it simulates the powder bed (1 mm x 0.2 mm x 0.06 mm) where the temperature variations are more important. To reduce the computation time, the base plate (1 mm x 0.2 mm x 0.3 mm) is discretized with a coarser mesh. The height must ensure that the bottom border will not interfere with the surface temperature. Mapped mesh guarantees nodal consistency at the interface between powder and base. Convergence analysis has been done to validate the mesh size. Powder elements affected by the heat source, are assigned with different material properties, as it is shown in Figure 11.

Figure 11.

Boundary condition and different material properties: molten pool (red elements circular area on the top), powder bed (blue elements, fine-meshed volume on the top), and base plate (cyan elements, coarse-meshed volume on the bottom).

A mapped mesh employing hexahedral 8-node elements is adopted to reduce the computational cost while maintaining high thermal field resolution. Specifically, thermal brick elements called SOLID70 with the following characteristics are used:

• Conduction and enthalpy capabilities

• Eight nodes (no mid-edge node capability)

• Applicable to a 3-D transient thermal analysis

• Mapped mesh

Different thermal properties can be associated with the same element type, hence it makes it possible to distinguish the different behavior of powder, base plate, and molten pool. Referring to Figure 11, the base plate elements (cyan) are associated with constant thermal properties, as it is supposed that the base plate is not affected by the thermal field. This assumption helps to reduce the non-linearity. Different element properties are also associated with the layer in order to distinguish between the inert powder bed (blue elements) and the grains that undergo the melting process. These elements are depicted in red and represent the dimension of the molten pool. Different thermal properties have just been explained in Section 2.1.

As mentioned before, the algorithm is iterative and the system must be solved at each laser spot application. The diagram in Figure 12 shows how the solving algorithm is carried out.

Figure 12.

A flow chart representing the main algorithm.

4.2. Boundary condition

The system of equations resulting from Eq. (10) can be solved once the prescribed boundary conditions (BCs) have been substituted. In FE thermal analysis, the possible BCs are:

1. Imposed temperature

2. Imposed heat flux

3. Flux due to convection ruled by the temperature difference

4. Flux due to radiation ruled by the fourth power of the absolute temperature

Boundary conditions depend on how the system interacts with the external environment:

1. Top surface:

   The working chamber is filled with Argon to reduce the alloy powder oxidation. Natural convection applies overall on the top surface, apart from the localized area where the laser heat flux is imposed. Since the emitting radiation flux makes the analysis highly nonlinear, its effect is not considered here. To solve this problem, an empirical relationship has been proposed [9, 18], which combines the effect of radiation and convection into a lumped heat transfer coefficient.

2. Lateral surfaces:

   Since the powder conductivity is very low, lateral surfaces can be considered as adiabatic, hence the heat flux imposed is equal to zero ( $\mathbf{q}\left(\mathbf{x},\mathrm{t}\right)=0$ ).

3. Bottom surface:

   In SLM machines, the base plate is heated between 80°C and 130° C, depending on the machine model. Bottom nodes are constrained with imposed temperature or with convection conditions. In this work, the bottom surface is constrained with convection boundary condition. As a consequence, a convection coefficient must be chosen in order to reproduce the convective exchange conditions into the base plate.

BC applied to the numerical model is summarized in Figure 11.

Not only boundary conditions, but also initial conditions (ICs) are requested to solve the numerical model. Initial conditions can be imposed setting up a starting temperature for all the nodes. These temperatures are used in transient solutions as the first step temperatures, hence at a time equal to zero:

$\mathrm{T}\left(\mathbf{x},\mathrm{t}=0\dots {\mathrm{t}}^{\prime }\right)={\mathrm{T}}_0\left(\mathbf{x}\right)$ E14

Moreover, since transient solution occurs at each cycle, initial condition must also be set at the beginning of each load step. It follows that initial condition applied to load step n are the nodal temperature obtained from the solution at step n-1.

4.3. Mapping procedure

Element undergoing phase change must be continuously updated with different material properties to simulate the melting and cooling process. When the average temperature of an element is higher than the melting point, the element is provided with different material properties that allow tracking the molten pool behavior. ANSYS^® cannot easily change material properties while the transient solution is running, not even using restart options. Therefore, it is mandatory to solve the analysis before modifying material properties. During post-processing, the temperatures of each element are analyzed and the material properties are changed accordingly.

The iterative algorithm helps to keep the analysis simple, even if it requires the element properties to be deleted at the end of each iteration so that they must be continuously saved and resumed at the beginning of the next iteration. Moreover, the mesh and, hence, the element spatial location are not constant throughout the iteration due to the dynamic mesh refinement (see Section 5). In fact, the FE environment is rebuilt iteratively with different element densities depending on the laser location, as it is shown in Figure 13.

Figure 13.

Comparison between OLD_MESH and NEW_MESH.

OLD_MESH and NEW_MESH refer to listed mesh entities. Each row of the list contains elements and nodes tracking number, their spatial coordinates, and the related properties.

The procedure able to assign correctly the temperature and material properties between two different mesh environments is called mapping procedure and is carried out in sequence by MatLab^® and ANSYS^®. A mapping algorithm is a useful tool that is able to save nodal temperatures from the previous load step, evaluate and assign material type with respect to the element average temperature, and finally, restore the data in the subsequent iteration as initial conditions. To avoid misunderstanding, it is worth noticing the difference between nodal and element properties: temperatures are the values assigned to nodes, while the material number is assigned to the elements. Due to ANSYS^® programming language, different thermal behaviors can be assigned to the same element type using material numbers. This is the reason why in this work the expression material properties has been used with the same meaning as thermal properties.

The flowchart presented in Figure 14 helps to understand the mapping algorithm.

Figure 14.

A flowchart illustrating the mapping procedure between two subsequent spots.

At the beginning, the elements and nodes are listed by ANSYS^® in a file with the related material number and temperature. This occurs in the post-procedure step related to the n cycle (NEW_MESH). The file is imported into MatLab^® and compared with the previous mesh file, just saved before from the n-1 cycle (OLD_MESH). Referring to Figure 13, elements and nodes are compared with respect to their spatial location and divided into two groups: common and uncommon entities. The dashed squares in Figure 13 highlight the difference between common and uncommon mesh.

Data coming from the previous analysis (step n-1) are assigned to the next one (step n) regardless of the grouped entities. Since the common elements and nodes share the same spatial location, the properties are simply transferred from the OLD_MESH to the NEW_MESH. The mapping algorithm takes the (element and node) spatial coordinates from OLD_MESH and searches the corresponding location in NEW_MESH (notice that the reference point for the element localization is the centroid). The mapping based on spatial coordinates is needed as the common entities share the same location but not the same tracking number because of the different meshes. Consequently, the temperatures and material numbers are transferred to NEW_MESH and can be simply assigned as the initial condition with respect to the entities number.

Regarding the uncommon entities, their properties cannot be directly transferred to the model because there is a spatial mismatch between the two meshes. The solution is to perform an interpolation. The elements and nodes locations from OLD_MESH are the input values, while the NEW_MESH entities are the target. The interpolation scheme scans all the elements (nodes) in the NEW_MESH, searching for the location that suits better the elements (nodes) belonging to the OLD_MESH. When the best solution is found, the properties can be interpolated between the two meshes. A different interpolation scheme is applied for temperature and material properties. The material number is assigned to the target with respect to the nearest OLD_MESH element and no interpolation is needed. However, the temperatures are assigned to target nodes with a more complicated scheme: not only the nearest node from OLD_MESH is chosen, but also a group of surrounding nodes that properly fit the target. Therefore, the temperature is assigned by means of an interpolation scheme that can be performed on a surface (2D interpolation) or on a volume (3D interpolation) regardless of the precision requested in the analysis. Finally, the interpolated value can be transferred to NEW_MESH and used for the initial condition.

The framework for the mapping procedure requests a lot of time for being set up, because it needs a strong interaction between ANSYS^® and MatLab^®. MatLab^® is used for grouping the entities and for mapping the common entities. ANSYS^® is chosen for the uncommon nodes taking advantage of the in-built powerful interpolation algorithm. Despite the complexity, a mapping method based on common and uncommon entities guarantees a strong reduction in the computational cost. The bottleneck of a traditional mapping procedure is the time-consuming interpolation algorithm. With this solution, the interpolation is applied only to a limited number of elements and not to the entire domain.

5. Mesh refinement

FEA is a useful tool to return an approximation of physics variables. Obviously, the computation time is a decisive factor to make numerical simulations competitive with respect to trial-and-error experiments. Traditionally, the bottleneck of a transient FEA is the time requested to compute the temperature field at each laser beam position. This gets even worse considering nonlinear material properties, the high amount of load steps, and elements number. In fact, every load step is divided into sub-steps to satisfy transient time integration rules. Consequently, the factor mainly responsible for the prolonged simulation time is the amount of elements and sub-steps; thus, in order to minimize the computational cost, this must be reduced as much as possible.

Suppose that the model has been built applying a uniform mapped mesh to the entire domain. Moreover, the mesh density has been increased as much as possible since FEM can predict more accurate results when the number of elements is high. Generally, it is often recommended to increase the elements density in the neighborhood of a certain zone where the results are requested to be more accurate. A typical example is when the stress concentration factor of a shouldered shaft subjected to an axial force needs to be determined. In such a case, the elements are concentrated in the vicinity of the fillet in order to obtain more reliable results.

However, the simulation time is proportional to the elements number and expected to increase enormously. One possible solution to overcome this long computing time is to use the dynamic mesh refinement (DMR) approach. It involves an independent mesh refinement of multiple sub-domains. This strategy allows to further refine, independently, the meshes in a hierarchic manner to reach a higher resolution.

It is mainly composed of two parts:

• Mesh refinement: increase element density in a region while having a coarse mesh in the remaining domain.

• Dynamic: the mesh is dynamically adapted according to the problem’s nature, e.g., boundary conditions, constitutive laws or geometry.

We will adopt a dynamic mesh refinement so as to adjust the local mesh refinement to the position of the laser spot. This has the great advantage of solving load steps with much fewer elements, hence the simulation time will benefit too. The dynamic part is a priority in this case, since there is a need to iteratively update the mesh at the end of each load step according to the laser spot position. The mesh can be rearranged in many ways; one option is to divide the mesh into different levels of refinement. As shown in Figure 15, three levels of increasing refinement degree have been implemented, namely level 1,2, and 3.

Figure 15.

An example of dynamic mesh refinement approach.

There are essentially two methods in order to build the mesh with different refinement levels, viz. bottom-up and top-down approaches. The former builds the entire domain starting with the coarsest mesh (level 1) and subsequently digs and removes the elements in order to generate a mesh of level 2 and so on for all the refinement levels. The latter method differs inasmuch as the finest mesh is built at the beginning and the remaining meshes are generated accordingly. ANSYS^® does not let the user modify the mesh once it is created, thus the second method was preferred over the first one.

5.1. Bonded contact technique

In order to impose discontinuous mesh levels to work properly, there is a need to connect the two parts to restore the continuity of the field variable. As it can be seen from Figure 15, mesh level 2 and 3 do not share the same nodes; hence, there is no mesh continuity between the two parts. Mesh compatibility was intentionally lost in order to further reduce the number of elements outside the heat-affected zone (HAZ). As a matter of fact, there are two main techniques to ensure continuity between incompatible meshes: bonded contact and constraint equations. Since the latter introduces additional constraint equations thus increasing the computational cost and the memory request, DMR based on bonded contact is preferred. The state of these contact elements never changes throughout the simulation, whereby not introducing additional sources of nonlinearity.

DMR requires an additional routine that permits data transfer from the previous mesh to the newly created and adapted mesh, as it has been explained in Section 4. Figure 16 shows the flowchart related to the DMR procedure. Moreover, it helps to understand how the mapping procedure is matched to well-fit the DMR requirements.

Figure 16.

A flow chart for the mesh refinement procedure.

At the beginning, the ANSYS^® simulation of the first laser spot is solved and data including mesh and nodal temperature are stored in the external file OLD_MESH. Subsequently, the spot position is moved and the mesh is updated and saved as NEW_MESH. At this point, ANSYS^® stops working and MatLab^® ad-hoc procedure will take OLD_MESH and NEW_MESH as inputs. The routine will sort nodes as common and uncommon and will return the information to ANSYS^®. At this point, the temperature field of the previous mesh can be applied to the new mesh as an initial condition in the following way:

• Common new nodes will have the same temperature as the old ones

• Uncommon new nodes will have an interpolated temperature from the old ones due to the *MOPER APDL command.

As a result, the simulation time dropped from 15 minutes/spot to 71 seconds/spot reducing the calculation time by 92%. As already mentioned before, the main parameter, which greatly affects the simulation time, is the number of substeps. The optimal value thereof was found through a convergence analysis based on the plot shown in Figure 17. It can be noted that an increment in the time step size $∆t$ has a much more pronounced effect on the solution time rather than on maximum nodal temperature (a measure of the solution accuracy). A reasonable trade-off between accuracy and simulation time is a time step size of 3 μs, allowing for 1% error in maximum temperature estimation and a computation time reduction of 80%. It is worth noting that the overall time reduction, due to mesh refinement and time step size reduction, is equal to 98.5%.

Figure 17.

Time consumption and maximum temperature trend.

The uniform mesh model and the one implementing DMR was tested applying a laser beam on a straight line and the results are shown in Figure 21. It can be seen that the temperature scale has only some little negligible variations. The DMR model well represents the physical phenomena and is a trade-off between result accuracy and computation time.

6. Calibration

It has been proven that the computational performances can be strongly enhanced using a simplified numerical model. Nevertheless, results are greatly affected by the lack of accuracy due to several simplifications applied to the model. A calibration procedure is necessary to reduce this issue: the material properties and the boundary conditions can be modified trying to fit numerical results with experimental data. The comparison is based on the molten pool dimension measured from a solidified seam. Moreover, since only a single seam is needed for the calibration, the geometrical domain can be halved along the symmetry plane, saving computational cost. The main reason for the calibration is that ANSYS^® does not consider elements which behave as liquid elements. The only way to simulate melting and cooling is to change the material properties and in particular thermal conduction and enthalpy, even though conduction does not apply for liquids.

As a consequence, the thermal properties are not well-defined for temperatures above the melting point. In this situation, a convective parameter should be used. This consideration permits the change of the parameters as needed, trying to simulate the convective behavior with a fictitious conduction parameter. The calibration procedure aims to adjust the thermal properties (only above the melting point) whereby fitting the molten pool size to the experimental data. The parameters adjusted by the calibration are enthalpy and thermal conductivity: the enthalpy is modified to control the temperature field while the conductivity is mainly responsible for the size of the molten pool. Specifically, metallographic inspections, as shown in Figure 19, are used to estimate width and depth of the molten pool.

Figure 18 shows a single molten seam obtained overlapping multiple layers. The seam was melted using the process parameters listed in Table 1. It has been cut and analyzed in order to gather information about the width and depth of the molten pool. Measured values are:

• Width = 183 ± 38 μm

• Depth = 107 ± 38 μm

Figure 18.

Comparison between the constant mesh and the dynamic mesh refinement.

The deviation related to measurements is mainly related to the narrow geometry and tiny dimensions of the object. Its width is only 4-5 times larger than the dimensions of the metal powder particles; therefore, the profile is not regular. It represents the minimum thickness that can be obtained with a single scan of the laser beam on the powder bed. Notice from Figure 19 that the molten seam undergoes the re-melting process with the application of successive layers and therefore the depth is not a reliable parameter. Nevertheless, the object is helpful in order to evaluate the real width of the molten pool.

Figure 19.

Micrograph of a single molten seam.

6.1. Calibration results

Due to the uncertainty related to the depth, only the molten pool width is taken into account, while the former issue will be addressed in future work. At the beginning, a trial simulation is carried out to check how the temperature field is sensitive to the parameters change. A directly measured thermal field is not available for this work; hence, the comparison is done with respect to results retrieved from the literature [20]. The enthalpy is indeed modified to keep the thermal field under control. New values for enthalpy are shown in Table 6. Only the last enthalpy value is modified increasing the specific heat for vapor by a factor of 10. This helps to decrease the maximum nodal temperature.

Enthalpy (J/m^3)
Powder	Solid	Powder
5000	5.2448e + 11	5.2448e + 11

Table 6.

Calibrated value for enthalpy.

The thermal behavior of molten pool is shown in Figure 20, which gives an idea about how elevated is the temperature of the zone irradiated by the laser.

Figure 20.

Thermal behavior of the molten pool.

This is due to the fact that the molten pool width is narrower than the experimental data and the conductivity needs to be increased. The calibration is, in a nutshell, an iterative algorithm that changes the conductivity with a trial factor as long as the numerical data well reproduce the experimental measurement. The algorithm involves MatLab^® and ANSYS^® as it is shown in the diagram presented in Figure 21.

Figure 21.

A flow chart for the calibration procedure.

The correction factor for conductivity is shown in Table 7. Notice that the correction is applied only to those values above the melting temperature.

Thermal conductivity (W/mK)
Temperature (°C)	Solid	Powder
1660	14.9567*5	14.9567*5
3265	14.9567*5	14.9567*5
3295	7.4784*5	7.4784*5
5000	7.4784*5	7.4784*5

Table 7.

Calibrated value for thermal conductivity.

Results coming from calibration are shown in Figure 22.

Figure 22.

Molten pool behavior with calibrated parameters.

7. Results

After parameters calibration, the simulation of a SLM process can be carried out. Because of the high number of laser spots, the simulation must be applied only to a small portion of the powder bed. Only one layer is considered and the adopted laser scanning strategy is the meander path. In order to gather information about the thermal field evolution into the bed powder, the time evolution of the temperature field is sampled on a spot selected as a temperature probe. The corresponding results are shown in Figure 23. The small window shows powder bed and meandering path. The black point along the meandering path represents the probe, which the temperature graph refers to.

Figure 23.

Temperature development of a point at the surface of the powder bed.

A very narrow temperature peak can be noticed in Figure 23. The highest peak is due to the heat source applied directly onto the probe. The other peaks are related to the reheating of the solidified area as the heat source is applied on the surrounding areas. In this example, the scan element is subject to remelting only once. The melting and cooling process occurs with very high gradients and this is the main source of the thermal residual stresses affecting the as-built parts.

8. Conclusions

A three-dimensional FE model is developed using ANSYS^® to study the thermal behavior of the molten pool in building a single layer via SLM process. At the beginning, the scanning strategy adopted by the laser is simulated by a path simulator built using MatLab. Then, the FE analysis framework is extensively explained with special regard to thermal properties applied to the model. Dynamic mesh refinement is used to reduce the computational cost of the simulation. Special care is taken in devising a mapped mesh discretization scheme, ensuring that the traveling subdomain centered on the laser spot changes as less as possible the mesh of the remaining subdomain. Finally, a calibration procedure is applied to fit the numerical results with the experimental measurements. The simulation results agree reasonably well with experimental and literature results and give some insight into the mutual interaction among the process parameters. Useful indications can be gained to optimize the process parameters, to estimate the adhesion between the layers, and to identify the best building strategy. This model can be further developed by incorporating the nodal temperature field into a structural analysis for predicting the resulting stress and strain field.

Acknowledgments

This work is part of the FAMAC Research Project, co-sponsored by Eurocoating S.p.A. and Provincia Autonoma di Trento (Regional Public Authority).