4.1. Numerical model

The governing heat transfer equation can be written as:

$$-\left(\frac{\partial q\_{\text{x}}}{\partial x} + \frac{\partial q\_{\text{y}}}{\partial y} + \frac{\partial q\_{\text{x}}}{\partial z}\right) + \mathbf{Q} = \rho c \frac{\partial T}{\partial t} \tag{8}$$

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

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

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

$$\mathbf{q}\_x = -\mathbf{k}\_x \frac{\partial T}{\partial z}$$

Notice that for isotropic, material thermal conductivity is the same: kx ¼ ky ¼ kz.

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

$$\text{div}([\mathbf{k}] \cdot \nabla T) + \mathbf{Q} = \rho \mathbf{c} \frac{\partial T}{\partial t} \tag{9}$$

is valid only when the convection boundary conditions apply. Tf is the temperature of the

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

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

139

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

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

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

environment into the working chamber.

numerical solution, while preserving results accuracy.

where k½ �¼ k I½ � 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:

$$\begin{aligned} \left[\mathbf{C}\_{\text{thermal}}\right] \left\{ \dot{\mathbf{T}}(\mathbf{t}) \right\} + \left[\mathbf{k}\_{\text{thermal}}\right] \left\{ \mathbf{T}(\mathbf{t}) \right\} &= \left\{ \mathbf{F}(\mathbf{t}) \right\} \qquad \qquad \qquad \mathbf{t} \in \{\mathbf{0}, \mathbf{t}'\} \end{aligned} \tag{10}$$

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

i. The thermal stiffness matrix kthermal ½ � is expressed as follows:

$$\left[\mathbf{k}\_{\text{thermal}}\right] = \int\_{\text{V}} \left[\mathbf{B}\_{\text{thermal}}\right]^{\text{T}} \left[\mathbf{k}\right] \left[\mathbf{B}\_{\text{thermal}}\right] \text{dV} + \int\_{\text{S}} \mathbf{h} \left[\mathbf{N}\right]^{\text{T}} \mathbf{N} \,\text{ds} \tag{11}$$

where Bthermal ½ � 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, kthermal ½ � has dimension 8 � 8. 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:5 w=m<sup>2</sup>K for Argon [20].

ii. The thermal specific heat matrix Cthermal ½ � is expressed as follows:

$$\left[\mathbf{C}\_{\text{thermal}}\right] = \int\_{\text{V}} \left[\mathbf{N}\right]^{\text{T}} \left[\mathbf{N}\right] \,\rho c \,dV\tag{12}$$

where N½ �¼ ½ � N1N2N3…N8 are the three-dimensional nodal shape functions of size 1 � 8. Once computed, Cthermal ½ � has dimension 8 � 8. r is the mass density and c is the specific heat.

iii. The thermal flux vector f g F tð Þ is expressed as follow:

$$\cdot \{ \mathbf{F}(\mathbf{t}) \} = \int\_{S} (\mathbf{q}(\mathbf{x}, \mathbf{t}) \cdot \hat{\mathbf{n}}) \cdot [\mathbf{N}]^{\mathrm{T}} d\mathbf{S} + \int\_{S} \mathbf{h}[\mathbf{N}]^{\mathrm{T}} \mathbf{T}\_{\mathrm{f}} \, d\mathbf{S} \tag{13}$$

where <sup>q</sup> is the input heat flux depending on boundary conditions and <sup>x</sup> <sup>¼</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> � �<sup>T</sup> is the position vector. dS denotes the surface area of the element. The second surface integral in Eq. (13) is valid only when the convection boundary conditions apply. Tf is the temperature of the environment into the working chamber.

qz ¼ �kz

General formulation of governing differential equation can be obtained substituting Fourier's

div kð Þþ ½ �∙∇T Q ¼ rc

Eq. (9) can be decomposed into a weak formulation [30] as the following system of differential

\_ð Þ � � <sup>þ</sup> kthermal ½ �f g T tð Þ <sup>¼</sup> f g F tð Þ <sup>t</sup><sup>∈</sup> <sup>0</sup>; <sup>t</sup>

½ � k Bthermal ½ �dV þ

where Bthermal ½ � 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, kthermal ½ � has dimension 8 � 8. 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

ð S

Notice that for isotropic, material thermal conductivity is the same: kx ¼ ky ¼ kz.

law component in Eq. (8):

equation of first order in t:

specific heat.

where k½ �¼ k I½ � is the thermal conductivity matrix.

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

Cthermal ½ � T t

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

kthermal ½ �¼

i. The thermal stiffness matrix kthermal ½ � is expressed as follows:

ð V

which has been assumed to be 12:5 w=m<sup>2</sup>K for Argon [20].

ii. The thermal specific heat matrix Cthermal ½ � is expressed as follows:

iii. The thermal flux vector f g F tð Þ is expressed as follow:

f g F tð Þ ¼

ð S

Cthermal ½ �¼

ð V

ð Þ q xð Þ ; <sup>t</sup> <sup>∙</sup>n<sup>b</sup> <sup>∙</sup>½ � <sup>N</sup> <sup>T</sup>

where <sup>q</sup> is the input heat flux depending on boundary conditions and <sup>x</sup> <sup>¼</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> � �<sup>T</sup> is the position vector. dS denotes the surface area of the element. The second surface integral in Eq. (13)

where N½ �¼ ½ � N1N2N3…N8 are the three-dimensional nodal shape functions of size 1 � 8. Once computed, Cthermal ½ � has dimension 8 � 8. r is the mass density and c is the

> dS þ ð S h N½ �<sup>T</sup>

Bthermal ½ �<sup>T</sup>

∂T ∂z

∂T

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

<sup>0</sup> f g (10)

h N½ �TNds (11)

Tf dS (13)

½ � <sup>N</sup> <sup>T</sup>½ � <sup>N</sup> <sup>r</sup>c dV (12)

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

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

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

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:


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.
