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


Table 1. Process parameters used in the simulation.

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 (r), specific heat (c), and temperature (T) according to Eq. (1):

$$\mathbf{H} = \int \rho \mathbf{c}(\mathbf{T}) \mathbf{d} \mathbf{T} \tag{1}$$

• 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

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

http://dx.doi.org/10.5772/intechopen.71876

131

While the conductivity value is taken from experimental data (or from literature), the enthalpy

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

<sup>2</sup> average specific heat (2)

specific heat for transition (3)

ð Þ TL � TS enthalpy at liquid temperature (5)

HS <sup>¼</sup> <sup>r</sup>CSð Þ TS � <sup>T</sup><sup>0</sup> enthalpy at solid temperature (4)

<sup>H</sup><sup>þ</sup> <sup>¼</sup> HL <sup>þ</sup> <sup>r</sup>CLð Þ <sup>T</sup> � TL enthalpy above liquid temperature (6)

Temperature (�C) Solid Powder 6.1874 0.6187 6.5660 0.6566 12.2620 1.2262 18.1490 1.8149 23.4590 2.3459 27.8010 2.7801 29.9134 2.9913 14.9567 14.9567 14.9567 14.9567 7.4784 7.4784 7.4784 7.4784

is calculated from Eq. (1). Variables needed for the equation are collected in Table 3.

high-temperature simplifications apply. Data are shown in Table 2.

and can be easily calculated from the equations below.

L TL � TS

Cavg <sup>¼</sup> Cs <sup>þ</sup> CL

<sup>C</sup><sup>∗</sup> <sup>¼</sup> Cavg <sup>þ</sup>

Thermal conductivity (W/mK)

HL <sup>¼</sup> HS <sup>þ</sup> <sup>r</sup>C<sup>∗</sup>

Table 2. Variation of thermal conductivity with temperature.

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:


• 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.

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.

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.

$$\mathbf{C\_{avg}} = \frac{\mathbf{C\_s} + \mathbf{C\_L}}{\mathbf{2}} \tag{2}$$

$$\mathbf{C}^\* = \mathbf{C}\_{\text{avg}} + \left(\frac{\mathbf{L}}{\mathbf{T}\_\text{L} - \mathbf{T}\_\text{S}}\right) \tag{3}$$

$$\mathbf{H}\_{\rm S} = \rho \mathbf{C}\_{\rm S} (\mathbf{T}\_{\rm S} - \mathbf{T}\_{0}) \tag{4}$$

$$\mathbf{H}\_{\rm L} = \mathbf{H}\_{\rm S} + \rho \mathbf{C}^{\*} (\mathbf{T}\_{\rm L} - \mathbf{T}\_{\rm S}) \tag{6} \tag{5}$$

$$\mathbf{H}\_{+} = \mathbf{H}\_{\mathrm{L}} + \rho \mathbf{C}\_{\mathrm{L}} (\mathrm{T} - \mathrm{T}\_{\mathrm{L}}) \qquad \text{(enthalpy above liquid temperature)} \tag{6}$$


Thermal conductivity (W/mK)

aluminum and 4% of Vanadium. It has excellent biocompatibility properties, as well as good

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

130 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

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 (r), specific heat (c), and temperature

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

• 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 con-

• 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.

rcð Þ T dT (1)

H ¼ ð

properties of the liquid phase, the following simplifying assumptions are taken:

mechanical properties [12].

Table 1. Process parameters used in the simulation.

Process parameters

(T) according to Eq. (1):

ductor than the solid.

Table 2. Variation of thermal conductivity with temperature.


Table 3. Thermal parameters.

where CL: specific heat of solid; r: density; TS: solidus temperature; TL: Liquidus temperature; T0: 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.

#### 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]:

$$\mathbf{I} = \frac{\mathbf{2AP}}{\pi \mathbf{r}\_{\text{max}}^2} \cdot \mathbf{e} \begin{pmatrix} -2\frac{r^2}{r\_{\text{max}}^2} \end{pmatrix} \tag{7}$$

where A: powder absorbance; P: laser powder; and rmax: laser beam radius is defined as the radius in which the power density is reduced from the peak value by a factor of e2 [29].

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

http://dx.doi.org/10.5772/intechopen.71876

133

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

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

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

coordinate system. Figure 6 shows an example of tessellation.

Figure 5. Laser beam irradiance with a Gaussian distribution.

3. Path simulation

principle of the slicing algorithm.

Enthalpy (J/m<sup>3</sup> )


Table 4. Variation of enthalpy with temperature.

Finite Element Thermal Analysis of Metal Parts Additively Manufactured via Selective Laser Melting http://dx.doi.org/10.5772/intechopen.71876 133

Figure 5. Laser beam irradiance with a Gaussian distribution.

where A: powder absorbance; P: laser powder; and rmax: laser beam radius is defined as the radius in which the power density is reduced from the peak value by a factor of e2 [29].

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

## 3. Path simulation

where CL: specific heat of solid; r: density; TS: solidus temperature; TL: Liquidus temperature;

Thermal parameters Solid Powder

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

) 4220 2532

Table 4 shows the values of enthalpy calculated with the previous equations. The column

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

Temperature (�C) Solid Powder Equation

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

<sup>∙</sup><sup>e</sup> �<sup>2</sup> <sup>r</sup><sup>2</sup> r2 max 

(7)

T0: reference temperature; T: saturation temperature; and L: latent heat.

132 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

laser irradiance and can be represented as a Gaussian distribution [18]:

0 00

3295 1.4210e + 10 1.4210e + 10 5000 6.5237e + 10 6.5237e + 10

2.2. Heat source

Enthalpy (J/m<sup>3</sup>

)

Table 4. Variation of enthalpy with temperature.

Table 3. Thermal parameters.

Density (kg/m<sup>3</sup>

named equation refers to previously numbered equations used for the calculation.

<sup>I</sup> <sup>¼</sup> <sup>2</sup>AP πr<sup>2</sup> max 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.

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.

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

Vertices below the slicing plane Vertices above the slicing plane

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

http://dx.doi.org/10.5772/intechopen.71876

135

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

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

it is shown in Figure 8. Coordinates of the intersection point are called Zp and Xp.

 ! A 9 ! A ! E 10 ! B ! B 12 ! E ! C 14 ! D ! D 15 ! C

Table 5. Vertices are grouped with respect to the slicing plane.

points. Figure 9 shows an example that can help to understand the procedure.

Figure 8. The intersection points are determined through geometrical relations.

Figure 6. An example of CAD model tessellation.

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.

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


Table 5. Vertices are grouped with respect to the slicing plane.

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 Zp and Xp.

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.

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

Figure 8. The intersection points are determined through geometrical relations.

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

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

Table 5.

Figure 6. An example of CAD model tessellation.

134 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 7. A 2D example of a triangles tessellation.

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

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

http://dx.doi.org/10.5772/intechopen.71876

137

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 (bulkpowder) transformations, only material properties are defined, rather than complex thermodynamic models. Applications involving phase change can be approached using ANSYS®

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 proce-

just for illustrating how the laser moves on the powder surface.

through elements with enthalpy property capabilities.

The governing heat transfer equation can be written as:

� <sup>∂</sup>q<sup>x</sup> ∂x þ ∂q<sup>y</sup> ∂y þ ∂q<sup>z</sup> ∂z

where qx, qy and qz are components of heat flow through unit area. According to Fourier's law:

∂T ∂x

∂T ∂y

qx ¼ �kx

qy ¼ �ky

þ Q ¼ rc

∂T

<sup>∂</sup><sup>t</sup> (8)

dure as will be extensively explained later.

4.1. Numerical model

4. Modeling approach

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

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.
