- Email: [email protected]
Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process
Accepted Manuscript Title: Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process Author: Manuela Galati Luca Iuliano Alessandro Salmi Eleonora Atzeni PII: DOI: Reference:
S2214-8604(16)30107-5 http://dx.doi.org/doi:10.1016/j.addma.2017.01.001 ADDMA 145
To appear in: Received date: Revised date: Accepted date:
20-5-2016 21-12-2016 2-1-2017
Please cite this article as: Manuela Galati, Luca Iuliano, Alessandro Salmi, Eleonora Atzeni, Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process, http://dx.doi.org/10.1016/j.addma.2017.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process Manuela Galati*, Luca Iuliano, Alessandro Salmi and Eleonora Atzeni Politecnico di Torino, Department of Management and Production Engineering (DIGEP), Torino, Italy Manuela Galati (Corresponding Author) Politecnico di Torino, Department of Management and Production Engineering (DIGEP) Corso Duca degli Abruzzi, 24 - 10129 Torino, Italy email: [email protected] phone: +39 011 090.7280 fax: +39 011 090.7299
Abstract The combination of additive manufacturing principles and electron beam (EB) technology allows complex metal parts, featuring excellent quality material, to be produced, whenever traditional methods are expensive or difficult to apply. Today, the optimization of process parameters, for a given metal powder, is generally attained through an empirical trial and error approach. Process simulation can be used as a tool for decision-making and process optimization, since a virtual analysis can help to facilitate the possibility of exploring ‘‘what if’’ scenarios. In this work, a new type of modelling has been introduced for energy source and powder material properties and it has been included in a thermal numerical model in order to improve the effectiveness and reliability of Electon Beam Melting (EBM) FE simulation. Several specific subroutines have been developed to automatically calculate the powder properties as temperature functions, and to consider the position of the beam during scanning as well as the material state changes from powder to liquid in the melting phase and from liquid to solid during cooling. A comparison of the numerical results and experimental data taken from literature has shown a good forecasting capability. The average deviations of the simulation from an experimental scan line width have been found to be below about 15%.
1
Keywords Electron Beam Melting (EBM); Additive Manufacturing; Metal powders; Simulation; Thermal modelling
1. Introduction There is a growing interest in industry in the possibility of increasing flexibility and reducing production costs and lead-time. These aims can be achieved by adopting more efficient and effective manufacturing processes. Additive manufacturing (AM) technologies offer the potential of responding to this request, as they can revolutionize the manufacturing approach [1-3]. Today, additive processes allow for the direct production of complex functional or end-usable metal parts [4]. The AM philosophy is that each object can be fabricated using a layer-by-layer process. This process involves three key phases: modelling the part using three-dimensional computer aided design (CAD) data, slicing the part into cross sections (layers) and building the part by overlapping layers, according to the CAD data. Numerous additive techniques are available for metals and they mainly differ on the basis of the construction method [3]. However, each AM process involves three main aspects: energy, material, and data or process control [5]. In particular, AM for metal can be classified as non-melting or melting processes, depending on which heat energy causes the partial or full melting of the raw material. The electron beam melting (EBM) process belongs to melting processes category, during which an electron beam (energy) melts metallic powders (material) to obtain a specify shape and solid massive part (data). The EBM process is currently being used in several applications in the aerospace and medical fields. The aim is to fabricate complex parts with excellent material quality [6]. Titanium superalloy components are examples of manufacturing by means of the EBM process, and for the use of traditional methods could be difficult and expensive [6]. Hence, EBM is not only a viable alternative to traditional processes, but also a process that offers exclusive benefits. The quality of electron beam melted components is affected to great extent by the position of the part in the job chamber and by the size of the component. Moreover, the EBM process must be properly calibrated for each new powder. Currently, a long and expensive empirical trial and error approach is usually adopted to identify an appropriate combination of process parameters for a given metal powder. Process simulation could be used as a suitable alternative tool for decision-making and process optimisation, since virtual analysis 2
facilitates the possibility of exploring “what if” scenarios [7]. In fact, temperature distribution and melt pool size simulations may help to better control of the quality of the built part. Therefore, in order to define a reliable simulation model and to optimise the quality of produced parts, the research community is currently studying and modelling the phenomena that occur during the EBM process. The simulation of the EBM process involves three fundamental aspects: a suitable analytical model of the material properties, a modelling of the thermal energy source, and as a consequence, a numerical implementation in order to solve the complexity of the phenomena, which takes into account the aforementioned aspects. Different approaches have been reported in the literature, and the main differences between them regard the modelling of powders. In the research developed by Körner, Attar and Heinl [8] and subsequent works [9-11], the powder was particle-toparticle modelled through a mesoscopic approach. In that study, researchers developed a tool that considers the numerical behaviour of each single powder particle. An algorithm was used to simulate the sequential addition of powder layers to generate a model composed of spherical particles. Despite the accuracy of the approach adopted to represent the EBM process, the simulation procedure is complex and expensive, in terms of time calculation. Instead, if the powder is considered as a continuum, it is possible to use the FE method as way of approximating the solution, and thus decreasing the calculation time significantly. The main challenges of this approach pertain to necessity of having to apply a suitable analytical model of the powder material properties and of having to take into account the material state change, because powder materials should have the same properties as the bulk material once the powders are molten and under cooling. Three dimensional FE numerical models have been developed by Qi, Yan, Lin and Zhang [12], Zäh, Lutzmann, Kahnert and Walchshäusl [13] and subsequent work [14], Shen and Chou [15] and subsequent works [16-21]. Qi, Yan, Lin and Zhang [12] defined the solid thermal conductivity as a linear function of porosity and neglected the emissivity of powder and the material state change. Zäh, Lutzmann, Kahnert and Walchshäusl [13] used the Zenher/Bauer/Schlünder model to model the density, specific heat and thermal conductivity of powders. Nevertheless, they did not consider the material state change and the emissivity as a temperature function [22]. In the works [15-21], the researchers assumed that the specific heat and the latent heat of fusion of the powders were the same as those of the solid bulk material, and modelled the thermal conductivity according to the studies of Tolochko, Arshinov, Gusarov, Titov, Laoui and Froyen [23]. Furthermore, they assumed that emissivity is independent of temperature and they used the Sih and Barlow’s model [22] to model powder emissivity. However, 3
conflicting hypotheses were introduced in this model because Tolochko’s model [23], which is used for thermal conductivity, considers the body-center cubic (BCC) packing of powders, while Sih and Barlow’s model [22], which is instead used for emissivity, considers the random packing of powders. In addition, the numerical implementation of the thermal behaviour of the material, after beam passage and the subsequent cooling phase, has not been described clearly in the literature. However, all of these studies ignored the heat dissipated by viscous forces and convection. In addition, in general, perfect wetting was assumed and the irradiation phenomenon was neglected [8, 12]. As far as the heat source modelling is concerned, when similarity between EBM and the selective laser melting (SLM) process is assumed, the Gaussian volumetric distribution is often chosen to model the heat source intensity in the same way as for laser energy [13-16, 18, 24]. Owing to the lack of knowledge on the EBM process, the simulation is tuned using coefficients to improve the match between experimental and numerical temperature results. These coefficients are introduced in the energy source model as efficiencies and reduce the amplitude of the Gaussian distribution [14, 19]. Qi, Yan, Lin and Zhang [12] considered the energy source as a uniform volumetric distribution, adjusted arbitrarily by a coefficient, which affects a single prismatic element. Many researchers have validated their numerical models considering EB welding (EBW) or laser processes [15, 18]. The models validated considering EBM experiments have shown an error of the simulated melt pool width of more than 50% [19]. Thus, the currently used FE models need to be improved to obtain a better understanding of what occurs during the EBM process. In this work, new analytical formulations of the energy source and of the properties of the powders are proposed for EBM process FE simulation. The study has focused on the energy source model going beyond the limit of the link with SLM, by evaluating the physical impact of the electron beam on the preheated powder bed. The aim has been to simplify the energy source model and to observe the main phenomena that occur during the interaction between electrons and powder bed. As far as the powder material is concerned, BCC lattice has been assumed for the packing of the powders and analytical formulations have been provided for thermal conductivity and emissivity. These new formulations have been implemented in a three-dimensional FE thermal model. The aim has been to develop an affordable FE numerical model that is able to predict three main quantities as well as their evolution over time: 1) the temperature distribution; 2) the conductive effects and 3) the melt pool size (width, length and depth). Several subroutines have been developed to automatically compute the EBM
4
phenomena, including the powder properties, the material state change and the energy source. Finally, in order to verify the applicability of this approach, the numerical model has been validated against literature data on the EBM process.
2. Modelling Properties of Powder Material The powders used in the EBM process are usually produced by gas atomization, and the obtained particle shape is spherical [25]. The main parameters that characterize the powders particles are size distribution and relative density [26]. The density of the powder bed is usually assumed to be equal to 68% of the corresponding solid density [16]. Hence, it can be assumed that the powder can be approximated as a packed bed of equal spheres with a mean radius, while the relative packing density corresponds to the BCC atomic packing factor (Figure 1). Therefore, porosity 𝜑 is equal to 0.32 [27]. Preheating before the melting phase sinters the powder particles that form the circular necks [19] as can be seen in Figure 1.
If the powder bed is modelled as a continuum, it is possible to define the material properties which are independent of any particular coordinate system. Hence, the powder density 𝜌𝑝𝑜𝑤𝑑𝑒𝑟 can be defined using the bulk density 𝜌𝑏𝑢𝑙𝑘 as: 𝜌𝑝𝑜𝑤𝑑𝑒𝑟 = (1 − 𝜑)𝜌𝑏𝑢𝑙𝑘
(1)
where the density values are measured in kilograms per millimeter cubed (kgmm-3). As far as the thermal properties of powder material are concerned, many researchers have indicated that the specific heat and the latent heat of powders can be considered the same those of solid material [14, 15, 22, 23]. The thermal conductivity and the emissivity were modelled as explained in the subsequent sections.
2.1. Thermal conductivity Since the hypothesis of the BCC structure and the vacuum in the built chamber are satisfied, the Tolochko model can be applied to estimate the thermal conductivity of the powder [23]. This model 5
takes into account radiation heat through the pores and conduction heat between adjacent particles. The thermal conductivity of the powder 𝜆𝑝𝑜𝑤𝑑𝑒𝑟 can be written as: 𝜆𝑝𝑜𝑤𝑑𝑒𝑟 = 𝜆𝑟 + 𝜆𝑐
(2)
where 𝜆𝑟 is the effective thermal conductivity due to radiation while 𝜆𝑐 is due to the heat transfer through the necks in the BCC lattice, and all conductivity values are expressed in WK-1mm-1. 𝜆𝑟 and 𝜆𝑐 can be computed as: 𝜆𝑟 =
16 3
𝑙𝜎𝑇 3
𝜆𝑐 = √3𝜆𝑏𝑢𝑙𝑘 𝑥
(3) (4)
where 𝑙 , expressed in millimetres, is the mean free path of the photon between the scattering events, 𝜎
is the Stefan-Boltzmann constant whose value is 5.67·10–14 Wmm-2·K-4, 𝑇 is the absolute
temperature, 𝜆𝑏𝑢𝑙𝑘 is the thermal conductivity of the bulk material and 𝑥 is the relative size of the necks, which is assumed equal to 0.09, according to the experimental measurements of Tolochko, Arshinov, Gusarov, Titov, Laoui and Froyen [23].
2.2. Emissivity According to Sih and Barlow [22], the emissivity of the powder bed depends on the emission from the powder particles and the emission from the holes in the powder bed. Consequently, emissivity can be written as: 𝜀𝑝𝑜𝑤𝑑𝑒𝑟 = 𝐴𝐻 𝜀𝐻 + (1 − 𝐴𝐻 )𝜀𝑏𝑢𝑙𝑘
(5)
where 𝐴𝐻 is expressed in mm2 and is the surface fraction that is occupied by radiation emitting holes, and 𝜀𝐻 is the emissivity of the holes, which is defined as: 𝜀𝐻 = 𝜀
𝜀𝑏𝑢𝑙𝑘 𝑏𝑢𝑙𝑘 +𝑓(1−𝜀𝑏𝑢𝑙𝑘 )
(6)
where 𝜀𝑏𝑢𝑙𝑘 is the emissivity of the bulk material and 𝑓 is the fraction of the total cavity surface that is cut away by the emitting hole.
6
Sih and Barlow [22] provided the expressions of the parameters 𝐴𝐻 and 𝑓 for a randomly packed bed. In this paper, a formulation of the parameters 𝐴𝐻 and 𝑓 has been obtained, according to the BCC structure hypothesis, that is: 𝐴𝐻
=
𝑓=
ℎ𝑜𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑎𝑐𝑒
(7)
𝑡𝑜𝑡𝑎𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 ℎ𝑜𝑙𝑒 𝑠𝑒𝑐𝑡𝑖𝑜𝑛
(8)
𝑡𝑜𝑡𝑎𝑙 ℎ𝑜𝑙𝑒 𝑠𝑢𝑟𝑓𝑎𝑐𝑒
The hole surface is a function of the number of holes per unit of cell. Considering spherical holes with diameter 𝑑, expressed in millimetres, it is possible to estimate the number of holes 𝑛ℎ𝑜𝑙𝑒𝑠 per volume of unit cell 𝑉𝐶𝐸𝐿𝐿 . The number of holes of the BCC structure can be determined from the following definition, where 𝑟𝐵𝐶𝐶 is the ratio between the hole diameter and particle diameter. 𝑛ℎ𝑜𝑙𝑒𝑠 =
𝜑𝑉𝐶𝐸𝐿𝐿 𝑉ℎ𝑜𝑙𝑒
2.94𝜑
=𝑟
𝐵𝐶𝐶
(9)
3
Assuming that 𝑟𝐵𝐶𝐶 is equal to 0.4142 [27], it results that there are approximately 13.24 holes per unit cell or 6.62 holes per sphere. Therefore, the area fraction, 𝐴𝐻 , occupied by the pores is 𝐴𝐻 =
𝑛ℎ𝑜𝑙𝑒𝑠
𝜋𝑑2 4
0.735𝜑
𝜋𝑑2 𝜋𝐷2 𝑛ℎ𝑜𝑙𝑒𝑠 +𝑛𝑠𝑝ℎ𝑒𝑟𝑒𝑠 4 4
= 0.735𝜑+0.5𝑟
𝐵𝐶𝐶
(10)
Similarly, 𝑓 can be calculated from the definition as the ratio between the cross section and the total surface of the hole. Therefore, the parameters 𝐴𝐻 and 𝑓 are equal to about 0.53 and 0.25, respectively, for BCC structure.
3. Heat Source Modelling In EBM processes, the electron beam scans the top surface and melts the powder material in order to build a layer. The use of an electron beam to melt metallic powder distinguishes EBM from SLM because the nature of the energy of the electrons is different from that of the photons. In fact, in the EBM process, the thermal energy to melt, sinter, heat, or vaporize material is obtained through the conversion of the kinetic energy of the electrons into thermal energy [28]. In order to obtain a proper modelling of the energy source, the impact of the electron beam on the preheated powder bed was studied by means of simulation with the Monte Carlo method. This method 7
estimates the trajectories of the electrons within a sample, taking into account the mean free paths of the electrons and the probability of interaction of the phenomena that occur during collision. The factors that play key roles are the material density, the alloy elements and their atomic number, the energy or accelerating voltage of the beam, and the focus beam diameter. The simulations were run for several materials and several focus beam diameters that are generally used in the EBM process, assuming the beam had a perpendicular impact with the top surface. Figure 2 shows the trajectories of the electrons in a cross section of a 316L sample of the preheated powder layer, with the density calculated according to (7) and a mean diameter of particle equal to 0.075 mm. The simulation is run with an acceleration voltage of 60 kV and a beam diameter of 0.400 mm. The back-scattered electrons are represented by red lines, whereas the blue lines refer to the not back-scattered electrons. The grey lines in Figure 2 show the distribution of the electrons along the radius and along the layer depth. The simulation results show two significant phenomena. The first one is that the electron energy rapidly decreases as a function of the depth. The depth of the interaction volume is much smaller than the thickness of the layer, according to [29]. Above all, the penetration depth is negligible compared to the typical powder particle size. Thus, the heat conduction effects along the depth between the particles are more important than the energy transferred by electrons. The second phenomenon is a lateral spread of the beam, which, according to [30], results in a greater effective diameter, 𝐷𝐸 , than the nominal focus beam diameter. Hence, 𝐷𝐸 is the diameter that characterises the total circle surface where the electrons interact with the powders. Furthermore, if the section of energy distribution per area is observed, it can be seen that the energy inside the focus beam diameter is similar to a flattened Gaussian distribution, except for the presence of a single central peak. In Figure 2, it is possible to notice that 𝐷𝐸 is approximately 30% larger than the focus diameter. Consequently, the heat source model in the simulation is described by means of a bi-dimensional distribution, because the kinetic energy transfer of the electrons, which only affects the top surface of the layer and the penetration depth, is negligible. In addition, in order to simplify the model, the real distribution is flattened, considering a uniform distribution on a circle portion of the top surface, which is described by focus beam diameter 𝐷. However, the surface heat flux applied in the simulation is multiplied by a corrective coefficient of global flux 𝜂 so as to take into account power conservation. For instance, considering a mono-dimensional Gaussian distribution, Figure 3 graphically shows how the model of the energy source has been simplified, according to the aforementioned considerations. The intensity of heat source 𝑞, expressed in Wmm-2, is calculated from: 8
𝑞=𝜂
𝑈𝐼
(11)
𝑆
where 𝑈 is the acceleration voltage measured in kilovolts, 𝐼 is the beam current, expressed in milliamps, and 𝑆 = 𝜋𝐷2 ⁄4 is the area of 𝛤, and is measured in millimetres square.
The 𝜂 formulation is provided to ensure power conservation that considers the real distribution of electrons by 𝐷𝐸 . This coefficient is expressed in Eq. 12 and calculated by means of numeric integration. 1
(𝑥1 −𝑥10 )2 +(𝑥2 −𝑥20 )2
𝜂 = ∫𝛤 2𝜋𝜎2 𝑒𝑥𝑝 (
2𝜎2
) dΓ
(12)
where, using the millimetre as the unit of length, 𝑥1 , 𝑥2 are the generic Cartesian coordinates, 𝑥10 , 𝑥20 are the coordinates of the starting point of the centre of the focal spot of the electron beam and 𝜎 = 𝐷𝐸 ⁄4 is the standard deviation necessary to consider the effective beam diameter and lateral spread. When this approach is adopted, the 𝜂 coefficient becomes a function of 𝐷𝐸 , and thus of the material, of the acceleration voltage and of the focus beam diameter and no adjustment needs to be applied to the flux value. For instance, for the previous example of the simulation of the stainless steel 316L, 𝜂 is equal to 0.22. If the focus beam diameter and the acceleration voltage are changed to 0.200 mm and 50 kV, respectively, the 𝜂 coefficient becomes equal to 0.13. In addition, in order to evaluate the difference, in terms of time calculation, the same thermal simulation was run considering the Gaussian distribution taken from literature and the uniform distribution proposed in this study. The calculation total time was considerably less when the uniform distribution is applied. In fact, when the uniform distribution was used, it was observed that the time increment for convergence was of a greater order of magnitude than that obtained when the Gaussian function was applied. Furthermore, the total heat flux applied to the surface using a Gaussian distribution was found to be dependent to great extent on the mesh size.
9
4. FE Thermal model Several complex physical phenomena occur during the EBM process. After the application of the beam energy, heat transfer takes place through three different mechanisms: through heat conduction between the powder particles, through heat conduction between the powder bed and the bulk substrate, and through irradiation from the powder bed to the chamber. During the process, the material first undergoes a phase change from powder to molten metal and then a second phase change in which the metal solidifies again as bulk material [15]. The analysis of the process is complex, because the scan speed is high and the phase change occurs quasi-instantaneously. Moreover, the temperature distribution is not homogeneous within the part/layer. In order to obtain fully dense parts without delamination, the molten pool must wet the previously consolidated material and the surrounding powder. In addition, after the melt pool solidifies, the top surface should be smooth enough to enable the powder of the next layer to spread, and a good surface finish to be obtained [31]. Therefore, the wetting of the solid phase by the liquid phase and capillary forces are important for the process [28]. In order to ensure a correct wetting, the energy beam must be able to melt both the top material and a portion of the bulk substrate. Moreover, with the formation of the melt pool, the heat transfer increases and produces strong thermal gradients. The temperature distribution throughout the depth causes the generation of flotation forces, while the surface temperature distribution leads to turbulent flows, such as Marangoni convection flows. Furthermore, the phase change, i.e. melting or solidification of the material during the process, also affects the heat transfer. In order to simplify the analysis, under the assumption of perfect wetting and negligible capillarity forces and evaporation, only the main phenomena that cause the heat transfer have been considered in the model, such as conduction between the powder particles and between the powder bed and the bulk substrate, as well as irradiation from the powder bed to the chamber. According to [32], the effects of Marangoni convection flows are irrelevant for higher scan speeds than 100 mm/s, which are common in EBM processes. Fundamental physics laws can be applied when the powder bed and the bulk substrate are modelled as a continuum. Thus, conservation of mass, conservation of momentum, and conservation of energy were used to obtain differential equations that describe the thermal behaviour during the process. Energy conservation was considered to calculate the temperature distribution. The absolute temperature distribution 𝑇(𝑥1 , 𝑥2 , 𝑥3 ) in a homogeneous body, within which there is no matter 10
exchange, is expressed in Cartesian coordinates. In this model, the electron beam is perpendicular to the 𝑥1 𝑥2 -plane. The energy balance per unit volume within an infinitesimal control volume of dimensions 𝑑𝑥1 𝑑𝑥2 𝑑𝑥3 and in differential form at a given time instant can be expressed as: 𝐷𝑒
−∇ ∙ 𝒒 = 𝜌 𝐷𝑡
(13) 𝐷𝑒
where 𝒒 is the surface heat flux vector, 𝜌 = 𝜌(𝑇) is the density and 𝐷𝑡 is the material derivate of the thermal energy density. The thermal energy density 𝑒 can be written as: 𝑒 = c𝑇 + ∆ℎ
(14)
where c is the specific heat, 𝑇 = 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 𝑡) is the temperature, as a function of both space and time 𝑡, and ∆ℎ is the latent enthalpy, which is defined as: 𝐿
𝑇 ≥ 𝑇𝑙
𝑇−𝑇𝑠
∆ℎ = {𝑓𝑠 𝐿 = 𝑇 −𝑇 𝐿 𝑙
𝑠
𝑇𝑠 < 𝑇 < 𝑇𝑙
0
(15)
𝑇 ≤ 𝑇𝑠
where 𝑇𝑠 and 𝑇𝑙 are the solidus and liquidus temperatures, respectively, 𝑓𝑠 is the solid fraction and 𝐿 is the latent heat of fusion, measured in Jkg-1. The surface heat flux vector 𝒒 is described by Fourier's law as: 𝒒 = −𝜆∇𝑇
(16)
where 𝜆 = 𝜆(𝑇) is the thermal conductivity, expressed in WK-1mm-1. The starting and boundary conditions are: 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 0) = 𝑇𝑝𝑟𝑒ℎ𝑒𝑎𝑡 with (𝑥1 , 𝑥2 , 𝑥3 ) ∈ 𝐷 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 0) = 𝑇𝑟 with (𝑥1 , 𝑥2 , 𝑥3 ) ∉ 𝐷 𝑇(𝑥1 , 𝑥2 , 𝑥3 , ∞) = 𝑇𝑟 with (𝑥1 , 𝑥2 , 𝑥3 ) ∉ 𝐷 −𝜆
𝜕𝑇
| = 𝑞(𝑥1 , 𝑥2 , 𝑥3 , 𝑣, 𝑡) − 𝑞𝑟𝑎𝑑
𝜕𝑛 𝛤
−𝜆
𝜕𝑇
|
𝜕𝑛 Ω−𝛤
= −𝑞𝑟𝑎𝑑
where 𝐷 is the union between the substrate and the layer domains and Ω as the top surface of the layer. The heat flux is a function of the coordinates, beam speed 𝑣 and time 𝑡. 𝑇𝑝𝑟𝑒ℎ𝑒𝑎𝑡 and 𝑇𝑟 are the 11
preheating temperature of the build chamber temperature, respectively. Because of radiation, heat loss 𝑞𝑟𝑎𝑑 ,
expressed in Wmm-2, can be written as: 𝑞𝑟𝑎𝑑 = 𝜀𝜎(𝑇 4 − 𝑇𝑟 4 )
(17)
where 𝜀 = 𝜀(𝑇) is the emissivity and 𝜎 is the Stefan-Boltzmann constant, whose value is 5.67·10– 14
Wmm-2K-4. The heat loss due to convection is negligible because the EBM process is performed
under a high vacuum. The finite element technique is used to find the solution to the heat transfer analysis. A thermal simulation was implemented to predict the temperature, size and shape of the melt pool and scan line during the EBM process. The three-dimensional FE model was developed in Abaqus/Standard.
4.1. Geometry and Mesh The geometry of FE consists of one substrate and one thin layer above the substrate. The substrate represents the solid bulk that has already been processed. A single thick layer is modelled on the top of the substrate as unsintered powder. The mesh consists of 8-node linear heat transfer DC3D bricks according to Abaqus Documentation [33]. A specific finer mesh was used within a portion of the powder layer of where and in the proximity the heat flux was applied. This mesh strategy was used to obtain detailed results on the inside of the powder layer close to the incident electron beam. However, the size of the mesh elements was also chosen to avoid long running times and to ensure the absence of spurious oscillations in the solution, caused by a numerical relationship between the minimum usable time increment, the element size and the thermophysical properties [33]. The electron beam flux moves along the 𝑥1 axis to simulate a single track. Figure 4 shows the model configuration and the boundary conditions.
4.2. Heat Source An Abaqus DFLUX user subroutine has been used to apply the heat flux to the top surface of the layer. The location of this flux changes with time due to the movement of the beam. The DFLUX subroutine is used to define the flux distribution as a function of the position and time. The user subroutine DFLUX is called at the beginning of each time increment and at each flux integration point. It then reads the simulation time and calculates the current position of the centre of the electron beam. The user code includes the motion law of the beam (scanning mode), the analytical formulation for the 12
energy source and the equation to select the Γ surface where the heat flux is applied. For example, equation (18) simulates the movement of the beam spot along the 𝑥1 axis 𝐷 2
(𝑥1 − 𝑥10 − 𝑥1̇ 𝑡)2 + (𝑥2 − 𝑥20 )2 ≤ ( ) 2
(18)
where 𝑥1̇ is the beam scan speed, expressed in mms-1, 𝐷 is the beam diameter and 𝑥10 and 𝑥20 are the coordinates of the starting point of the beam centre along the 𝑥1 and 𝑥2 axis, respectively, and the millimetre is used as the unit of length.
4.3. Material Properties Owing to the large temperature range involved in the EBM process, the thermophysical properties should be expressed as functions of the temperature. However, the thermal behaviour of the powder and of the bulk material is significantly different. Consequently, when the temperature is higher than the melting point, the thermo-physical properties of the powder are the ones that correspond to the properties of the liquid metal. After melting and during cooling, the thermophysical properties become those of the bulk. The melting phase of the powders causes the material state change from powder properties to the bulk properties. Thus, the thermal behaviour is a function of the temperature and of the material state. The temperature dependent values of the bulk material have been extracted from technical databases. The properties of the corresponding powder material have been calculated from the respective bulk material by applying the specific models described in Section 2. 4.3.1. Material properties, phase change and material state change The user subroutine UMATHT has been developed to define the thermal behaviour as a function of the temperature and the material state (bulk or powder). A solution-dependent state variable, named MAT_ID, works like an index of the material state for each node of each element. Initially, MAT_ID is set to 0 and the powder properties are assigned to the nodes of the elements. When the powder melts, MAT_ID changes to 1 and the bulk properties are assigned to the material calculation points. The subroutine UMATHT reads the MAT_ID value and the temperature at each material calculation point at the start of the increment, and then calculates and updates the material properties expressed in the previous section. In addition, according to (15), the subroutine UMATHT modifies the enthalpy function in order to take into account latent heat effects and the solid fraction during the 13
melting and cooling phases. Finally, the subroutine UMATHT updates the MAT_ID values, if required, according to the temperature solution and the current value of MAT_ID. The MAT_ID value can no longer change from 1 back to 0. Thus, UMATHT is used to compute the specific heat, the thermal conductivity, and the material state change and the phase change. The required inputs are listed below:
Thermal conductivity of the bulk material (𝜆𝑏𝑢𝑙𝑘 ) as a function of the temperature
Specific heat (𝑐) as a function of the temperature
Latent heat
Solidus temperature
Liquidus Temperature
Powder porosity
Neck ratio
Particle size
Stefan-Boltzmann constant
4.3.2. Emissivity No user subroutine to define the emissivity as a function of the temperature and of the material state is available in the Abaqus Standard. In order to model the heat radiation loss, a pseudo heat transfer coefficient ℎ𝑟𝑎𝑑 has been defined in equation (19) and introduced into a user subroutine called FILM. ℎ𝑟𝑎𝑑 = 𝜀(𝑇)𝜎(𝑇 2 + 𝑇𝑟 2 )(𝑇 + 𝑇𝑟 )
(19)
The pseudo heat transfer coefficient has units of measurement of the film coefficient and it includes emissivity as a function of the temperature. The required inputs are the emissivity 𝜀(𝑇) of the bulk material and the chamber temperature. 4.3.3. Density The user subroutine USDFLD is used to update the density value from the powder density to the bulk density after melting. This subroutine is linked to the subroutine GETVRM which reads the MAT_ID information and calculates the corresponding density.
4.4. How the FE model works 14
The flowchart in Figure 5a shows the model structure after the preheating phase. Initially, the MAT_ID and temperature values are read by the user subroutines UMATHT and FILM, at the start of the time increment for each material calculation point. The temperature dependent properties are computed according to the obtained MAT_ID value and the density is updated. The user subroutine DFLUX then reads the current total time and applies the heat flux at each flux integration point according to the motion law of the beam. The software then solves the heat transfer problem. When the temperature value at the material calculation points is higher than the solidus temperature and the material is powder, the MAT_ID value changes from 0 to 1. The MAT_ID variable is updated, if necessary, in the user subroutine UMATH. The cycle is repeated for each time increment until the simulation ends.
The flowchart in Figure 5b shows the calculation method for when the material change is not considered. The difference is significant. In fact, in this case, the properties of the material are only temperature dependent. Consequently, the powder thermo-physical properties are the ones that correspond to the properties of the liquid metal when the temperature is higher than the melting point. However, when the temperature drops, the properties again become those of the powder, even through the material is actually in bulk form. This approach leads to unreliable numerical results because the powder properties are significantly different from the corresponding bulk material, for example, those of the density and the thermal conductivity (Figure 6). It is possible to explain unreliable numerical results by considering the thermal diffusivity of the material. The thermal diffusivity value describes the heat transfer rate of a material, which is equal to the thermal conductivity divided by the density and the specific heat. If the thermal diffusivity is calculated from the relative properties of the bulk and of the powder material (Figure 6), the thermal diffusivity value of the bulk is higher than that of the powder material for each temperature value. For this reason, the heat moves rapidly through the bulk material, as it is able to conduct heat quickly, while the heat transfer through the powder material is slow, that is, close to zero, because the powder has low thermal diffusivity. Thus, if the material change procedure is not considered in the simulation after passage of the beam, the heat transfer is calculated between the melt pool and the already scanned material that has the powder properties, instead of the bulk material. In this case, the temperature gradient is high and the heat transfer is low. Thus, the dimensions of the melt pool, in terms of width, length and depth, are reduced compared to those in which the surrounding material has changed
15
properties in the bulk material. In fact, in the case in which the material change procedure is deactivated, it has been observed, for instance, that the width of the scan line is approximately equal to the beam diameter. Instead, as will be shown in the next section, the width of the scan line is several times greater than the beam diameter, because the previously scanned material, which is cooling as bulk material, transfers the heat rapidly because of the conduction effects, and thus increases the heat affected area.
5. Model Validation The proposed model has been validated against the results of experiments available in the literature in terms of width of the scan line. The validation was based on the experiments conducted by Qi, Yan, Lin, He and Zhang [25]. In their experiments, several stainless steel 316L single scan lines were fabricated, varying the beam current and the scan speed. Material Properties Density at 293 K, bulk [kg/mm3]
7.97
Apparent density of the powder [kg/mm3]*
5.42
Solidus Temperature, TS [K]
1658
Liquidus Temperature, TL [K]
1699
Latent heat, L [kJ/kg]
273
Average particle size of the powders [mm]
0.075
Table 1 and Figure 6 report the thermophysical properties of the material as implemented in the model. The bulk material properties were extracted from [12, 14], whereas the powder properties were computed according to the models outlined in Section 2. Table 2 summarises the details of the experiments, related to the simulation inputs introduced to reproduce the real experiments. In order to resemble the test setup, the entire workpiece was modelled as powder material. The dimensions of the model are 10 × 10 × 10 mm3 and the finer mesh size is 0.05 mm.
The width of the melt pool is measured by the MAT_ID state variable. As can be observed in Figure 7 and Figure 8, it is possible to distinguish in the simulation results: the bulk material (soft grey, MAT_ID 0), the powder material (intense grey, MAT_ID 1) and the boundary line (red line). The 16
boundary line separates the powder area from the area that has undergone the material change. Each single powder particle near the boundary line can physically be in any of these positions (Figure 7): (a) tangent to the boundary line in the powder material part; (b) melted and solidified in the bulk material part, and thus theoretically tangent to the boundary line in the bulk material; (c) on the boundary line. The powder particles in position (a) are sintered, due to heat conduction, and as a consequence, they are part of the solidified scan line [34]. Similarly, the powder particles in position (c) belong to the scan line, but they are melted. However, the scan line width obtained from the simulation results does not take into account the particles in position (a) and/or (c), therefore the maximum dispersion of the simulation results is considered equal to twice the particle size. Figure 9 compares numerical and experimental measurements of the width of each single track. The numerical measurement of the width of the scan line refers to the maximum width, Wnumeric max. This approach allows the width of the scan line to be defined, even when the boundary line is irregular, as shown in Figure 8. The interval bars show the deviations from the numerical measurements. The extremes of the interval bar are: Wmax, which is the maximum scan width when at least one powder particle is in position (a) for each boundary line, and the minimum width of scan line Wnumeric
min
measured
considering the minimum distance between the boundary lines. According to the experimental results, the width of the scan line, for each combination of process parameters, is greater than of the focus beam diameter. The comparison shows a good agreement between the numerical results and the experimental data. The deviation between the experimental data and the corresponding Wnumeric max values is less than 25%. The maximum deviation drops to 13% when the beam current is higher than 0.5 mA. Furthermore, the experimental width is always included within the error bar of the numerical results, except for when the beam current and the scan speed are equal to 0.5 mA and 2 mm/s, respectively. According to main effects plot, Figure 10 highlights the average outcome for each value of each scan speed and beam current for the experimental and numerical results. The model appears to be slightly sensitive to low values of the scan speed. However, except for low values of the scan speed, the numeric response follows the experimental results when the scan speed and the beam current change. In fact, the numeric width of the scan line decreases as the scanning speed increase and increases as the beam current increases. The deviation between the experimental and numerical results may be attributed to:
17
the effect of the Marangoni convection flows that has been neglected. This hypothesis affects the numerical results for low scan speeds; if the scan speed were increased, the deviation between the numerical and experimental results would be reduced;
a balling effect that occurs for a low heat flux value [14]; the scan line contour is more blurred for a low scan speed [25]. Thus, the sides of the track are not clearly recognizable and the measurement of the width becomes more difficult. The simulation results partially show this behaviour, because the width of the scan line is indented for a low beam current while a constant width is shown for a higher beam current (Figure 8);
the melting line that soaks the powders. When the scanning speed is high, the melting line, in fact, obtains less energy from the electron beam and soaks less powder [25]. This phenomenon could be analysed by introducing the mechanical behaviour of the powders, but this is beyond the scope of this work.
However, the deviation is generally acceptable for any processes simulation. The numerical outcomes show that the model is able to correctly represent the phenomena that occur during the EBM process. Hence, the proposed modelling of a uniform distribution seems to be able to correctly model the energy source, as is the model that is used to model the thermophysical properties of the powder material.
6. Conclusions In this work, a three-dimensional thermal FE model of the EBM process has been developed in order to implement a new modelling of the energy source and powder material properties. First, the impact of the electron beam on a preheated powder bed was studied. This analysis led to the energy source being modelled as a uniform distribution, and a corrective coefficient of the global flux being introduced as a function of the material, the acceleration voltage and the focus beam diameter. Because of the high temperature gradient and the presence of different material states, i.e. powder and bulk, the thermophysical properties were modelled as a function of the temperature and material state. The properties of the bulk material were extracted from technical literature, whereas the material properties of the gas-atomized powder were calculated assuming a BCC packing structure and introducing a formulation for the emissivity. Several user subroutines were used to solve the heat transfer problem in Abaqus/Standard. The user subroutine DFLUX was implemented to apply a heat 18
flux and to move the electron beam spot. The user subroutines UMATHT, USDFLD, and FILM were used to calculate the properties of the powder material. Moreover, conduction effects, radiation losses, melting and cooling phases with a latent heat effect, and a coupled phase change and material change were considered. Experimental data taken from the literature were used as a reference for validation purposes to test the response of the numerical model to conductive effects and to different widths of the scan lines. However, the deviation was found to be less than 25% when the beam current was equal to 0.5 mA, and dropped to 13% when the beam current values were higher. These deviations are acceptable for a general response of the simulation processes. This limited deviation may be ascribed to the fact that certain phenomena had not been taken into consideration. The trend of the numerical results was found to be similar to that of the experimental data. This behaviour indicates that the model is able to correctly represent the main phenomena that occur during the EBM process. Hence, the proposed modelling of the material properties as well as the uniform distribution of the energy source seem to constitute the correct approach to simulate the EBM process. The model could be used as a tool to identify the most appropriate set of parameters, in terms of beam characteristics, powder properties and process control, and to consider their interactions. Improvements to the existing model will include a temperature analysis, in the hatching mode, and the implementation of the mechanical properties of the material in order to have a complete coupled thermo-mechanical model for the evaluation of stresses and deformations induced by beam passage.
19
References [1] E. Atzeni, L. Iuliano, P. Minetola, A. Salmi, Redesign and cost estimation of rapid manufactured plastic parts, Rapid Prototyping J 16(5) (2010) 308-317. [2] N. Hopkinson, P. Dickens, Rapid prototyping for direct manufacture, Rapid Prototyping J 7(4) (2001) 197-202. [3] T. Wohlers, Additive Manufacturing State of the Industry (Wohlers Report 2014), Fort Collins, USA2014. [4] E. Atzeni, A. Salmi, Economics of additive manufacturing for end-usable metal parts, The International Journal of Advanced Manufacturing Technology 62(9-12) (2012) 1147-1155. [5] J. Milberg, M. Sigl, Electron beam sintering of metal powder, Prod. Eng. Res. Devel. 2(2) (2008) 117-122. [6] S. Biamino, A. Penna, U. Ackelid, S. Sabbadini, O. Tassa, P. Fino, M. Pavese, P. Gennaro, C. Badini, Electron beam melting of Ti-48Al-2Cr-2Nb alloy: Microstructure and mechanical properties investigation, Intermetallics 19(6) (2011) 776-781. [7] M. Chiumenti, M. Cervera, A. Salmi, C. Agelet de Saracibar, N. Dialami, K. Matsui, Finite element modeling of multipass welding and shaped metal deposition processes, Computer Methods in Applied Mechanics and Engineering 199(3740) (2010) 2343-2359. [8] C. Körner, E. Attar, P. Heinl, Mesoscopic simulation of selective beam melting processes, Journal of Materials Processing Technology 211(6) (2011) 978-987. [9] R. Ammer, M. Markl, U. Ljungblad, C. Körner, U. Rüde, Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method, Computers & Mathematics with Applications 67(2) (2014) 318-330. [10] R. Ammer, U. Rüde, M. Markl, C. Körner, Modeling of Thermodynamic Phenomena with Lattice Boltzmann Method for Additive Manufacturing Processes. [11] T. Scharowsky, A. Bauereiß, R. Singer, C. Körner, Observation and numerical simulation of melt pool dynamic and beam powder interaction during selective electron beam melting, Proceedings of the Solid Freeform Fabrication Symposium, 2012, pp. 815-820. [12] H. Qi, Y. Yan, F. Lin, R. Zhang, Scanning method of filling lines in electron beam selective melting, Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 221(12) (2007) 1685-1694. [13] M. Zäh, S. Lutzmann, M. Kahnert, F. Walchshäusl, Determination of Process Parameters for Electron Beam Sintering (EBS), Excerpt from the Proceedings of the COMSOL Conference Hannover, 2008. [14] M.F. Zäh, S. Lutzmann, Modelling and simulation of electron beam melting, Production Engineering 4(1) (2010) 1523. [15] N. Shen, K. Chou, Thermal Modeling of Electron Beam Additive Manufacturing Process: Powder Sintering Effects, ASME 2012 International Manufacturing Science and Engineering Conference collocated with the 40th North American Manufacturing Research Conference and in participation with the International Conference on Tribology Materials and Processing, American Society of Mechanical Engineers, 2012, pp. 287-295. [16] N. Shen, Y. Chou, Numerical thermal analysis in electron beam additive manufacturing with preheating effects, Proceedings of the 23rd Solid Freeform Fabrication Symposium, Austin, TX, 2012, pp. 774-784. [17] N.G. Shen, K. Chou, Thermal Modeling of Electron Beam Additive Manufacturing Process Powder Sintering Effects, Proceedings of the Asme International Manufacturing Science and Engineering Conference, 2012 (2012) 287-295. [18] N. Shen, K. Chou, Simulations of Thermo-Mechanical Characteristics in Electron Beam Additive Manufacturing, ASME 2012 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2012, pp. 67-74. [19] B. Cheng, S. Price, J. Lydon, K. Cooper, K. Chou, On Process Temperature in Powder-Bed Electron Beam Additive Manufacturing: Model Development and Validation, J Manuf Sci E-T Asme 136(6) (2014). [20] S. Price, B. Cheng, J. Lydon, K. Cooper, K. Chou, On Process Temperature in Powder-Bed Electron Beam Additive Manufacturing: Process Parameter Effects, J Manuf Sci E-T Asme 136(6) (2014). [21] B. Cheng, S. Price, X. Gong, J. Lydon, K. Cooper, K. Chou, Speed Function Effects in Electron Beam Additive Manufacturing, ASME 2014 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2014, pp. V02AT02A003-V02AT02A003. [22] S.S. Sih, J.W. Barlow, Emissivity of powder beds, Solid Freeform Fabrication Symposium Proceedings, Center for Materials Science and Engineering, Mechanical Engineering Department and Chemical Engineering Department, the University of Texas at Austin, 1995, pp. 7-9. [23] N.K. Tolochko, M.K. Arshinov, A.V. Gusarov, V.I. Titov, T. Laoui, L. Froyen, Mechanisms of selective laser sintering and heat transfer in Ti powder, Rapid Prototyping Journal 9(5) (2003) 314-326. [24] E. Soylemez, J.L. Beuth, K. Taminger, Controlling Melt Pool Dimensions over a Wide Range of Material Deposition Rates in Electron Beam Additive Manufacturing, Solid Freeform Fabrication Proceedings, 2010, pp. 571-582.
20
[25] H.B. Qi, Y.N. Yan, F. Lin, W. He, R.J. Zhang, Direct metal part forming of 316L stainless steel powder by electron beam selective melting, P I Mech Eng B-J Eng 220(11) (2006) 1845-1853. [26] M. Krishnan, E. Atzeni, R. Canali, F. Calignano, D. Manfredi, E.P. Ambrosio, L. Iuliano, On the effect of process parameters on properties of AlSi10Mg parts produced by DMLS, Rapid Prototyping J 20(6) (2014). [27] R.M. German, Particle Packing Characteristics, Metal Powder Industries Federation1989. [28] E. Attar, Simulation of selective electron beam melting processes, Dr.-Ing., University of Erlangen, Nuremberg, Germany (2011). [29] K. Alexander, B. Andreas, K. Carolin, Modelling of electron beam absorption in complex geometries, Journal of Physics D: Applied Physics 47(6) (2014) 065307. [30] T.R. Mahale, Electron Beam Melting of Advanced Materials and Structures, Industrial Engineering, North Carolina State University, 2009. [31] M. Agarwala, D. Bourell, J. Beaman, H. Marcus, J. Barlow, Direct selective laser sintering of metals, Rapid Prototyping Journal 1(1) (1995) 26-36. [32] A. Gusarov, I. Yadroitsev, P. Bertrand, I. Smurov, Heat transfer modelling and stability analysis of selective laser melting, Applied Surface Science 254(4) (2007) 975-979. [33] Abaqus, ABAQUS Documentation, Dassault Systèmes, Providence (USA), 2014. [34] O. Cansizoglu, O. Harrysson, D. Cormier, H. West, T. Mahale, Properties of Ti–6Al–4V non-stochastic lattice structures fabricated via electron beam melting, Materials Science and Engineering: A 492(1) (2008) 468-474.
21
Figure 1 – (a) the elementary model of a powder (BCC structure) before the preheating phase, (b) the elementary model after the preheating phase.
Figure 2 – Electron trajectories in a cross section of a 316L solid sample. The electron trajectory simulation has been run by means of CASINO software [27]. The red lines refer to back-scattered electron trajectories, while the blue lines refer to the paths of non-back-scattered electrons. The grey lines illustrate the distribution of the electrons along the radius and along the layer depth. An acceleration voltage of 60 kV and a beam diameter of 0.400 mm have been adopted. The effective diameter 𝐷𝐸 is larger than the nominal focus beam diameter.
22
Figure 3 – Heat flux distribution along the radius. The grey dashed line is the Gaussian distribution function that is used in literature [13-21] that will be adjusted by the efficiencies; the black dashed line is the Gaussian distribution corrected by means of 𝐷𝐸 ; the continuous black line is the approximation of the Gaussian distribution corrected by means of 𝐷𝐸 .
Figure 4 - Model configuration and boundary condition
23
Figure 5 – FE model structure (the subroutines UMATHT and DFLUX) with the material change procedure (a) and without the material change procedure (b)
24
Figure 6 – Specific heat of the bulk material [12]. The specific heat of the powder material can be considered the same as that of the bulk material [14, 15, 22, 23]. Thermal conductivity and specific heat of the bulk [12] and powder material according to the model presented in Section 2.1. Emissivity value for bulk material [14] and powder material according to the model presented in Section 2.2.
25
Figure 7 – Possible position of a single powder particle: (a) tangent to the boundary line in the powder material part; (b) melted and solidified in the bulk material part and thus theoretically tangent to the boundary line in the bulk material; (c) over the boundary line.
Figure 8 – Scan line for 0.5 mA of beam current and 3.3 mm/s of scan speed. The indented scan line can be observed. The numeric maximum width of scan line Wnumeric max is approximately 0.55 mm. The minimum width Wnumeric min is equal to 0.45 mm. Considering as the powder particles near the boundary line could be sintered, the maximum width of the scan line, Wmax, is equal to the sum of Wnumeric max and twice the particle diameter (75µm).
26
Figure 9 – Comparison between the experimental and numerical results
27
Figure 10 – Plot of the main effects Table 1 Material properties of 316L Stainless Steel [12, 25]. Material Properties Density at 293 K, bulk [kg/mm3]
7.97
Apparent density of the powder [kg/mm3]*
5.42
Solidus Temperature, TS [K]
1658
Liquidus Temperature, TL [K]
1699
Latent heat, L [kJ/kg]
273
Average particle size of the powders [mm] * value calculated according to Eq. 1
0.075
Table 2 Simulation inputs on the basis of experimental parameters [25]. Simulations parameters Beam diameter, 𝐷 [mm]
0.200
Accelerating voltage, 𝑈 [kV]
50
𝜂*
0.13
Beam current, 𝐼 [mA] Scan speed, 𝑥1̇ [mm/s] Surface flux, 𝑞 [W/mm2]*
0.5 2.3
3.3
0.8 5.0
2.3
3.3
96.32 154.11 * values calculated according to Section 4.2
28
1.0 5.0
2.3
3.3 192.64
5.0
S2214-8604(16)30107-5 http://dx.doi.org/doi:10.1016/j.addma.2017.01.001 ADDMA 145
To appear in: Received date: Revised date: Accepted date:
20-5-2016 21-12-2016 2-1-2017
Please cite this article as: Manuela Galati, Luca Iuliano, Alessandro Salmi, Eleonora Atzeni, Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process, http://dx.doi.org/10.1016/j.addma.2017.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Modelling energy source and powder properties for the development of a thermal FE model of the EBM additive manufacturing process Manuela Galati*, Luca Iuliano, Alessandro Salmi and Eleonora Atzeni Politecnico di Torino, Department of Management and Production Engineering (DIGEP), Torino, Italy Manuela Galati (Corresponding Author) Politecnico di Torino, Department of Management and Production Engineering (DIGEP) Corso Duca degli Abruzzi, 24 - 10129 Torino, Italy email: [email protected] phone: +39 011 090.7280 fax: +39 011 090.7299
Abstract The combination of additive manufacturing principles and electron beam (EB) technology allows complex metal parts, featuring excellent quality material, to be produced, whenever traditional methods are expensive or difficult to apply. Today, the optimization of process parameters, for a given metal powder, is generally attained through an empirical trial and error approach. Process simulation can be used as a tool for decision-making and process optimization, since a virtual analysis can help to facilitate the possibility of exploring ‘‘what if’’ scenarios. In this work, a new type of modelling has been introduced for energy source and powder material properties and it has been included in a thermal numerical model in order to improve the effectiveness and reliability of Electon Beam Melting (EBM) FE simulation. Several specific subroutines have been developed to automatically calculate the powder properties as temperature functions, and to consider the position of the beam during scanning as well as the material state changes from powder to liquid in the melting phase and from liquid to solid during cooling. A comparison of the numerical results and experimental data taken from literature has shown a good forecasting capability. The average deviations of the simulation from an experimental scan line width have been found to be below about 15%.
1
Keywords Electron Beam Melting (EBM); Additive Manufacturing; Metal powders; Simulation; Thermal modelling
1. Introduction There is a growing interest in industry in the possibility of increasing flexibility and reducing production costs and lead-time. These aims can be achieved by adopting more efficient and effective manufacturing processes. Additive manufacturing (AM) technologies offer the potential of responding to this request, as they can revolutionize the manufacturing approach [1-3]. Today, additive processes allow for the direct production of complex functional or end-usable metal parts [4]. The AM philosophy is that each object can be fabricated using a layer-by-layer process. This process involves three key phases: modelling the part using three-dimensional computer aided design (CAD) data, slicing the part into cross sections (layers) and building the part by overlapping layers, according to the CAD data. Numerous additive techniques are available for metals and they mainly differ on the basis of the construction method [3]. However, each AM process involves three main aspects: energy, material, and data or process control [5]. In particular, AM for metal can be classified as non-melting or melting processes, depending on which heat energy causes the partial or full melting of the raw material. The electron beam melting (EBM) process belongs to melting processes category, during which an electron beam (energy) melts metallic powders (material) to obtain a specify shape and solid massive part (data). The EBM process is currently being used in several applications in the aerospace and medical fields. The aim is to fabricate complex parts with excellent material quality [6]. Titanium superalloy components are examples of manufacturing by means of the EBM process, and for the use of traditional methods could be difficult and expensive [6]. Hence, EBM is not only a viable alternative to traditional processes, but also a process that offers exclusive benefits. The quality of electron beam melted components is affected to great extent by the position of the part in the job chamber and by the size of the component. Moreover, the EBM process must be properly calibrated for each new powder. Currently, a long and expensive empirical trial and error approach is usually adopted to identify an appropriate combination of process parameters for a given metal powder. Process simulation could be used as a suitable alternative tool for decision-making and process optimisation, since virtual analysis 2
facilitates the possibility of exploring “what if” scenarios [7]. In fact, temperature distribution and melt pool size simulations may help to better control of the quality of the built part. Therefore, in order to define a reliable simulation model and to optimise the quality of produced parts, the research community is currently studying and modelling the phenomena that occur during the EBM process. The simulation of the EBM process involves three fundamental aspects: a suitable analytical model of the material properties, a modelling of the thermal energy source, and as a consequence, a numerical implementation in order to solve the complexity of the phenomena, which takes into account the aforementioned aspects. Different approaches have been reported in the literature, and the main differences between them regard the modelling of powders. In the research developed by Körner, Attar and Heinl [8] and subsequent works [9-11], the powder was particle-toparticle modelled through a mesoscopic approach. In that study, researchers developed a tool that considers the numerical behaviour of each single powder particle. An algorithm was used to simulate the sequential addition of powder layers to generate a model composed of spherical particles. Despite the accuracy of the approach adopted to represent the EBM process, the simulation procedure is complex and expensive, in terms of time calculation. Instead, if the powder is considered as a continuum, it is possible to use the FE method as way of approximating the solution, and thus decreasing the calculation time significantly. The main challenges of this approach pertain to necessity of having to apply a suitable analytical model of the powder material properties and of having to take into account the material state change, because powder materials should have the same properties as the bulk material once the powders are molten and under cooling. Three dimensional FE numerical models have been developed by Qi, Yan, Lin and Zhang [12], Zäh, Lutzmann, Kahnert and Walchshäusl [13] and subsequent work [14], Shen and Chou [15] and subsequent works [16-21]. Qi, Yan, Lin and Zhang [12] defined the solid thermal conductivity as a linear function of porosity and neglected the emissivity of powder and the material state change. Zäh, Lutzmann, Kahnert and Walchshäusl [13] used the Zenher/Bauer/Schlünder model to model the density, specific heat and thermal conductivity of powders. Nevertheless, they did not consider the material state change and the emissivity as a temperature function [22]. In the works [15-21], the researchers assumed that the specific heat and the latent heat of fusion of the powders were the same as those of the solid bulk material, and modelled the thermal conductivity according to the studies of Tolochko, Arshinov, Gusarov, Titov, Laoui and Froyen [23]. Furthermore, they assumed that emissivity is independent of temperature and they used the Sih and Barlow’s model [22] to model powder emissivity. However, 3
conflicting hypotheses were introduced in this model because Tolochko’s model [23], which is used for thermal conductivity, considers the body-center cubic (BCC) packing of powders, while Sih and Barlow’s model [22], which is instead used for emissivity, considers the random packing of powders. In addition, the numerical implementation of the thermal behaviour of the material, after beam passage and the subsequent cooling phase, has not been described clearly in the literature. However, all of these studies ignored the heat dissipated by viscous forces and convection. In addition, in general, perfect wetting was assumed and the irradiation phenomenon was neglected [8, 12]. As far as the heat source modelling is concerned, when similarity between EBM and the selective laser melting (SLM) process is assumed, the Gaussian volumetric distribution is often chosen to model the heat source intensity in the same way as for laser energy [13-16, 18, 24]. Owing to the lack of knowledge on the EBM process, the simulation is tuned using coefficients to improve the match between experimental and numerical temperature results. These coefficients are introduced in the energy source model as efficiencies and reduce the amplitude of the Gaussian distribution [14, 19]. Qi, Yan, Lin and Zhang [12] considered the energy source as a uniform volumetric distribution, adjusted arbitrarily by a coefficient, which affects a single prismatic element. Many researchers have validated their numerical models considering EB welding (EBW) or laser processes [15, 18]. The models validated considering EBM experiments have shown an error of the simulated melt pool width of more than 50% [19]. Thus, the currently used FE models need to be improved to obtain a better understanding of what occurs during the EBM process. In this work, new analytical formulations of the energy source and of the properties of the powders are proposed for EBM process FE simulation. The study has focused on the energy source model going beyond the limit of the link with SLM, by evaluating the physical impact of the electron beam on the preheated powder bed. The aim has been to simplify the energy source model and to observe the main phenomena that occur during the interaction between electrons and powder bed. As far as the powder material is concerned, BCC lattice has been assumed for the packing of the powders and analytical formulations have been provided for thermal conductivity and emissivity. These new formulations have been implemented in a three-dimensional FE thermal model. The aim has been to develop an affordable FE numerical model that is able to predict three main quantities as well as their evolution over time: 1) the temperature distribution; 2) the conductive effects and 3) the melt pool size (width, length and depth). Several subroutines have been developed to automatically compute the EBM
4
phenomena, including the powder properties, the material state change and the energy source. Finally, in order to verify the applicability of this approach, the numerical model has been validated against literature data on the EBM process.
2. Modelling Properties of Powder Material The powders used in the EBM process are usually produced by gas atomization, and the obtained particle shape is spherical [25]. The main parameters that characterize the powders particles are size distribution and relative density [26]. The density of the powder bed is usually assumed to be equal to 68% of the corresponding solid density [16]. Hence, it can be assumed that the powder can be approximated as a packed bed of equal spheres with a mean radius, while the relative packing density corresponds to the BCC atomic packing factor (Figure 1). Therefore, porosity 𝜑 is equal to 0.32 [27]. Preheating before the melting phase sinters the powder particles that form the circular necks [19] as can be seen in Figure 1.
If the powder bed is modelled as a continuum, it is possible to define the material properties which are independent of any particular coordinate system. Hence, the powder density 𝜌𝑝𝑜𝑤𝑑𝑒𝑟 can be defined using the bulk density 𝜌𝑏𝑢𝑙𝑘 as: 𝜌𝑝𝑜𝑤𝑑𝑒𝑟 = (1 − 𝜑)𝜌𝑏𝑢𝑙𝑘
(1)
where the density values are measured in kilograms per millimeter cubed (kgmm-3). As far as the thermal properties of powder material are concerned, many researchers have indicated that the specific heat and the latent heat of powders can be considered the same those of solid material [14, 15, 22, 23]. The thermal conductivity and the emissivity were modelled as explained in the subsequent sections.
2.1. Thermal conductivity Since the hypothesis of the BCC structure and the vacuum in the built chamber are satisfied, the Tolochko model can be applied to estimate the thermal conductivity of the powder [23]. This model 5
takes into account radiation heat through the pores and conduction heat between adjacent particles. The thermal conductivity of the powder 𝜆𝑝𝑜𝑤𝑑𝑒𝑟 can be written as: 𝜆𝑝𝑜𝑤𝑑𝑒𝑟 = 𝜆𝑟 + 𝜆𝑐
(2)
where 𝜆𝑟 is the effective thermal conductivity due to radiation while 𝜆𝑐 is due to the heat transfer through the necks in the BCC lattice, and all conductivity values are expressed in WK-1mm-1. 𝜆𝑟 and 𝜆𝑐 can be computed as: 𝜆𝑟 =
16 3
𝑙𝜎𝑇 3
𝜆𝑐 = √3𝜆𝑏𝑢𝑙𝑘 𝑥
(3) (4)
where 𝑙 , expressed in millimetres, is the mean free path of the photon between the scattering events, 𝜎
is the Stefan-Boltzmann constant whose value is 5.67·10–14 Wmm-2·K-4, 𝑇 is the absolute
temperature, 𝜆𝑏𝑢𝑙𝑘 is the thermal conductivity of the bulk material and 𝑥 is the relative size of the necks, which is assumed equal to 0.09, according to the experimental measurements of Tolochko, Arshinov, Gusarov, Titov, Laoui and Froyen [23].
2.2. Emissivity According to Sih and Barlow [22], the emissivity of the powder bed depends on the emission from the powder particles and the emission from the holes in the powder bed. Consequently, emissivity can be written as: 𝜀𝑝𝑜𝑤𝑑𝑒𝑟 = 𝐴𝐻 𝜀𝐻 + (1 − 𝐴𝐻 )𝜀𝑏𝑢𝑙𝑘
(5)
where 𝐴𝐻 is expressed in mm2 and is the surface fraction that is occupied by radiation emitting holes, and 𝜀𝐻 is the emissivity of the holes, which is defined as: 𝜀𝐻 = 𝜀
𝜀𝑏𝑢𝑙𝑘 𝑏𝑢𝑙𝑘 +𝑓(1−𝜀𝑏𝑢𝑙𝑘 )
(6)
where 𝜀𝑏𝑢𝑙𝑘 is the emissivity of the bulk material and 𝑓 is the fraction of the total cavity surface that is cut away by the emitting hole.
6
Sih and Barlow [22] provided the expressions of the parameters 𝐴𝐻 and 𝑓 for a randomly packed bed. In this paper, a formulation of the parameters 𝐴𝐻 and 𝑓 has been obtained, according to the BCC structure hypothesis, that is: 𝐴𝐻
=
𝑓=
ℎ𝑜𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑎𝑐𝑒
(7)
𝑡𝑜𝑡𝑎𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 ℎ𝑜𝑙𝑒 𝑠𝑒𝑐𝑡𝑖𝑜𝑛
(8)
𝑡𝑜𝑡𝑎𝑙 ℎ𝑜𝑙𝑒 𝑠𝑢𝑟𝑓𝑎𝑐𝑒
The hole surface is a function of the number of holes per unit of cell. Considering spherical holes with diameter 𝑑, expressed in millimetres, it is possible to estimate the number of holes 𝑛ℎ𝑜𝑙𝑒𝑠 per volume of unit cell 𝑉𝐶𝐸𝐿𝐿 . The number of holes of the BCC structure can be determined from the following definition, where 𝑟𝐵𝐶𝐶 is the ratio between the hole diameter and particle diameter. 𝑛ℎ𝑜𝑙𝑒𝑠 =
𝜑𝑉𝐶𝐸𝐿𝐿 𝑉ℎ𝑜𝑙𝑒
2.94𝜑
=𝑟
𝐵𝐶𝐶
(9)
3
Assuming that 𝑟𝐵𝐶𝐶 is equal to 0.4142 [27], it results that there are approximately 13.24 holes per unit cell or 6.62 holes per sphere. Therefore, the area fraction, 𝐴𝐻 , occupied by the pores is 𝐴𝐻 =
𝑛ℎ𝑜𝑙𝑒𝑠
𝜋𝑑2 4
0.735𝜑
𝜋𝑑2 𝜋𝐷2 𝑛ℎ𝑜𝑙𝑒𝑠 +𝑛𝑠𝑝ℎ𝑒𝑟𝑒𝑠 4 4
= 0.735𝜑+0.5𝑟
𝐵𝐶𝐶
(10)
Similarly, 𝑓 can be calculated from the definition as the ratio between the cross section and the total surface of the hole. Therefore, the parameters 𝐴𝐻 and 𝑓 are equal to about 0.53 and 0.25, respectively, for BCC structure.
3. Heat Source Modelling In EBM processes, the electron beam scans the top surface and melts the powder material in order to build a layer. The use of an electron beam to melt metallic powder distinguishes EBM from SLM because the nature of the energy of the electrons is different from that of the photons. In fact, in the EBM process, the thermal energy to melt, sinter, heat, or vaporize material is obtained through the conversion of the kinetic energy of the electrons into thermal energy [28]. In order to obtain a proper modelling of the energy source, the impact of the electron beam on the preheated powder bed was studied by means of simulation with the Monte Carlo method. This method 7
estimates the trajectories of the electrons within a sample, taking into account the mean free paths of the electrons and the probability of interaction of the phenomena that occur during collision. The factors that play key roles are the material density, the alloy elements and their atomic number, the energy or accelerating voltage of the beam, and the focus beam diameter. The simulations were run for several materials and several focus beam diameters that are generally used in the EBM process, assuming the beam had a perpendicular impact with the top surface. Figure 2 shows the trajectories of the electrons in a cross section of a 316L sample of the preheated powder layer, with the density calculated according to (7) and a mean diameter of particle equal to 0.075 mm. The simulation is run with an acceleration voltage of 60 kV and a beam diameter of 0.400 mm. The back-scattered electrons are represented by red lines, whereas the blue lines refer to the not back-scattered electrons. The grey lines in Figure 2 show the distribution of the electrons along the radius and along the layer depth. The simulation results show two significant phenomena. The first one is that the electron energy rapidly decreases as a function of the depth. The depth of the interaction volume is much smaller than the thickness of the layer, according to [29]. Above all, the penetration depth is negligible compared to the typical powder particle size. Thus, the heat conduction effects along the depth between the particles are more important than the energy transferred by electrons. The second phenomenon is a lateral spread of the beam, which, according to [30], results in a greater effective diameter, 𝐷𝐸 , than the nominal focus beam diameter. Hence, 𝐷𝐸 is the diameter that characterises the total circle surface where the electrons interact with the powders. Furthermore, if the section of energy distribution per area is observed, it can be seen that the energy inside the focus beam diameter is similar to a flattened Gaussian distribution, except for the presence of a single central peak. In Figure 2, it is possible to notice that 𝐷𝐸 is approximately 30% larger than the focus diameter. Consequently, the heat source model in the simulation is described by means of a bi-dimensional distribution, because the kinetic energy transfer of the electrons, which only affects the top surface of the layer and the penetration depth, is negligible. In addition, in order to simplify the model, the real distribution is flattened, considering a uniform distribution on a circle portion of the top surface, which is described by focus beam diameter 𝐷. However, the surface heat flux applied in the simulation is multiplied by a corrective coefficient of global flux 𝜂 so as to take into account power conservation. For instance, considering a mono-dimensional Gaussian distribution, Figure 3 graphically shows how the model of the energy source has been simplified, according to the aforementioned considerations. The intensity of heat source 𝑞, expressed in Wmm-2, is calculated from: 8
𝑞=𝜂
𝑈𝐼
(11)
𝑆
where 𝑈 is the acceleration voltage measured in kilovolts, 𝐼 is the beam current, expressed in milliamps, and 𝑆 = 𝜋𝐷2 ⁄4 is the area of 𝛤, and is measured in millimetres square.
The 𝜂 formulation is provided to ensure power conservation that considers the real distribution of electrons by 𝐷𝐸 . This coefficient is expressed in Eq. 12 and calculated by means of numeric integration. 1
(𝑥1 −𝑥10 )2 +(𝑥2 −𝑥20 )2
𝜂 = ∫𝛤 2𝜋𝜎2 𝑒𝑥𝑝 (
2𝜎2
) dΓ
(12)
where, using the millimetre as the unit of length, 𝑥1 , 𝑥2 are the generic Cartesian coordinates, 𝑥10 , 𝑥20 are the coordinates of the starting point of the centre of the focal spot of the electron beam and 𝜎 = 𝐷𝐸 ⁄4 is the standard deviation necessary to consider the effective beam diameter and lateral spread. When this approach is adopted, the 𝜂 coefficient becomes a function of 𝐷𝐸 , and thus of the material, of the acceleration voltage and of the focus beam diameter and no adjustment needs to be applied to the flux value. For instance, for the previous example of the simulation of the stainless steel 316L, 𝜂 is equal to 0.22. If the focus beam diameter and the acceleration voltage are changed to 0.200 mm and 50 kV, respectively, the 𝜂 coefficient becomes equal to 0.13. In addition, in order to evaluate the difference, in terms of time calculation, the same thermal simulation was run considering the Gaussian distribution taken from literature and the uniform distribution proposed in this study. The calculation total time was considerably less when the uniform distribution is applied. In fact, when the uniform distribution was used, it was observed that the time increment for convergence was of a greater order of magnitude than that obtained when the Gaussian function was applied. Furthermore, the total heat flux applied to the surface using a Gaussian distribution was found to be dependent to great extent on the mesh size.
9
4. FE Thermal model Several complex physical phenomena occur during the EBM process. After the application of the beam energy, heat transfer takes place through three different mechanisms: through heat conduction between the powder particles, through heat conduction between the powder bed and the bulk substrate, and through irradiation from the powder bed to the chamber. During the process, the material first undergoes a phase change from powder to molten metal and then a second phase change in which the metal solidifies again as bulk material [15]. The analysis of the process is complex, because the scan speed is high and the phase change occurs quasi-instantaneously. Moreover, the temperature distribution is not homogeneous within the part/layer. In order to obtain fully dense parts without delamination, the molten pool must wet the previously consolidated material and the surrounding powder. In addition, after the melt pool solidifies, the top surface should be smooth enough to enable the powder of the next layer to spread, and a good surface finish to be obtained [31]. Therefore, the wetting of the solid phase by the liquid phase and capillary forces are important for the process [28]. In order to ensure a correct wetting, the energy beam must be able to melt both the top material and a portion of the bulk substrate. Moreover, with the formation of the melt pool, the heat transfer increases and produces strong thermal gradients. The temperature distribution throughout the depth causes the generation of flotation forces, while the surface temperature distribution leads to turbulent flows, such as Marangoni convection flows. Furthermore, the phase change, i.e. melting or solidification of the material during the process, also affects the heat transfer. In order to simplify the analysis, under the assumption of perfect wetting and negligible capillarity forces and evaporation, only the main phenomena that cause the heat transfer have been considered in the model, such as conduction between the powder particles and between the powder bed and the bulk substrate, as well as irradiation from the powder bed to the chamber. According to [32], the effects of Marangoni convection flows are irrelevant for higher scan speeds than 100 mm/s, which are common in EBM processes. Fundamental physics laws can be applied when the powder bed and the bulk substrate are modelled as a continuum. Thus, conservation of mass, conservation of momentum, and conservation of energy were used to obtain differential equations that describe the thermal behaviour during the process. Energy conservation was considered to calculate the temperature distribution. The absolute temperature distribution 𝑇(𝑥1 , 𝑥2 , 𝑥3 ) in a homogeneous body, within which there is no matter 10
exchange, is expressed in Cartesian coordinates. In this model, the electron beam is perpendicular to the 𝑥1 𝑥2 -plane. The energy balance per unit volume within an infinitesimal control volume of dimensions 𝑑𝑥1 𝑑𝑥2 𝑑𝑥3 and in differential form at a given time instant can be expressed as: 𝐷𝑒
−∇ ∙ 𝒒 = 𝜌 𝐷𝑡
(13) 𝐷𝑒
where 𝒒 is the surface heat flux vector, 𝜌 = 𝜌(𝑇) is the density and 𝐷𝑡 is the material derivate of the thermal energy density. The thermal energy density 𝑒 can be written as: 𝑒 = c𝑇 + ∆ℎ
(14)
where c is the specific heat, 𝑇 = 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 𝑡) is the temperature, as a function of both space and time 𝑡, and ∆ℎ is the latent enthalpy, which is defined as: 𝐿
𝑇 ≥ 𝑇𝑙
𝑇−𝑇𝑠
∆ℎ = {𝑓𝑠 𝐿 = 𝑇 −𝑇 𝐿 𝑙
𝑠
𝑇𝑠 < 𝑇 < 𝑇𝑙
0
(15)
𝑇 ≤ 𝑇𝑠
where 𝑇𝑠 and 𝑇𝑙 are the solidus and liquidus temperatures, respectively, 𝑓𝑠 is the solid fraction and 𝐿 is the latent heat of fusion, measured in Jkg-1. The surface heat flux vector 𝒒 is described by Fourier's law as: 𝒒 = −𝜆∇𝑇
(16)
where 𝜆 = 𝜆(𝑇) is the thermal conductivity, expressed in WK-1mm-1. The starting and boundary conditions are: 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 0) = 𝑇𝑝𝑟𝑒ℎ𝑒𝑎𝑡 with (𝑥1 , 𝑥2 , 𝑥3 ) ∈ 𝐷 𝑇(𝑥1 , 𝑥2 , 𝑥3 , 0) = 𝑇𝑟 with (𝑥1 , 𝑥2 , 𝑥3 ) ∉ 𝐷 𝑇(𝑥1 , 𝑥2 , 𝑥3 , ∞) = 𝑇𝑟 with (𝑥1 , 𝑥2 , 𝑥3 ) ∉ 𝐷 −𝜆
𝜕𝑇
| = 𝑞(𝑥1 , 𝑥2 , 𝑥3 , 𝑣, 𝑡) − 𝑞𝑟𝑎𝑑
𝜕𝑛 𝛤
−𝜆
𝜕𝑇
|
𝜕𝑛 Ω−𝛤
= −𝑞𝑟𝑎𝑑
where 𝐷 is the union between the substrate and the layer domains and Ω as the top surface of the layer. The heat flux is a function of the coordinates, beam speed 𝑣 and time 𝑡. 𝑇𝑝𝑟𝑒ℎ𝑒𝑎𝑡 and 𝑇𝑟 are the 11
preheating temperature of the build chamber temperature, respectively. Because of radiation, heat loss 𝑞𝑟𝑎𝑑 ,
expressed in Wmm-2, can be written as: 𝑞𝑟𝑎𝑑 = 𝜀𝜎(𝑇 4 − 𝑇𝑟 4 )
(17)
where 𝜀 = 𝜀(𝑇) is the emissivity and 𝜎 is the Stefan-Boltzmann constant, whose value is 5.67·10– 14
Wmm-2K-4. The heat loss due to convection is negligible because the EBM process is performed
under a high vacuum. The finite element technique is used to find the solution to the heat transfer analysis. A thermal simulation was implemented to predict the temperature, size and shape of the melt pool and scan line during the EBM process. The three-dimensional FE model was developed in Abaqus/Standard.
4.1. Geometry and Mesh The geometry of FE consists of one substrate and one thin layer above the substrate. The substrate represents the solid bulk that has already been processed. A single thick layer is modelled on the top of the substrate as unsintered powder. The mesh consists of 8-node linear heat transfer DC3D bricks according to Abaqus Documentation [33]. A specific finer mesh was used within a portion of the powder layer of where and in the proximity the heat flux was applied. This mesh strategy was used to obtain detailed results on the inside of the powder layer close to the incident electron beam. However, the size of the mesh elements was also chosen to avoid long running times and to ensure the absence of spurious oscillations in the solution, caused by a numerical relationship between the minimum usable time increment, the element size and the thermophysical properties [33]. The electron beam flux moves along the 𝑥1 axis to simulate a single track. Figure 4 shows the model configuration and the boundary conditions.
4.2. Heat Source An Abaqus DFLUX user subroutine has been used to apply the heat flux to the top surface of the layer. The location of this flux changes with time due to the movement of the beam. The DFLUX subroutine is used to define the flux distribution as a function of the position and time. The user subroutine DFLUX is called at the beginning of each time increment and at each flux integration point. It then reads the simulation time and calculates the current position of the centre of the electron beam. The user code includes the motion law of the beam (scanning mode), the analytical formulation for the 12
energy source and the equation to select the Γ surface where the heat flux is applied. For example, equation (18) simulates the movement of the beam spot along the 𝑥1 axis 𝐷 2
(𝑥1 − 𝑥10 − 𝑥1̇ 𝑡)2 + (𝑥2 − 𝑥20 )2 ≤ ( ) 2
(18)
where 𝑥1̇ is the beam scan speed, expressed in mms-1, 𝐷 is the beam diameter and 𝑥10 and 𝑥20 are the coordinates of the starting point of the beam centre along the 𝑥1 and 𝑥2 axis, respectively, and the millimetre is used as the unit of length.
4.3. Material Properties Owing to the large temperature range involved in the EBM process, the thermophysical properties should be expressed as functions of the temperature. However, the thermal behaviour of the powder and of the bulk material is significantly different. Consequently, when the temperature is higher than the melting point, the thermo-physical properties of the powder are the ones that correspond to the properties of the liquid metal. After melting and during cooling, the thermophysical properties become those of the bulk. The melting phase of the powders causes the material state change from powder properties to the bulk properties. Thus, the thermal behaviour is a function of the temperature and of the material state. The temperature dependent values of the bulk material have been extracted from technical databases. The properties of the corresponding powder material have been calculated from the respective bulk material by applying the specific models described in Section 2. 4.3.1. Material properties, phase change and material state change The user subroutine UMATHT has been developed to define the thermal behaviour as a function of the temperature and the material state (bulk or powder). A solution-dependent state variable, named MAT_ID, works like an index of the material state for each node of each element. Initially, MAT_ID is set to 0 and the powder properties are assigned to the nodes of the elements. When the powder melts, MAT_ID changes to 1 and the bulk properties are assigned to the material calculation points. The subroutine UMATHT reads the MAT_ID value and the temperature at each material calculation point at the start of the increment, and then calculates and updates the material properties expressed in the previous section. In addition, according to (15), the subroutine UMATHT modifies the enthalpy function in order to take into account latent heat effects and the solid fraction during the 13
melting and cooling phases. Finally, the subroutine UMATHT updates the MAT_ID values, if required, according to the temperature solution and the current value of MAT_ID. The MAT_ID value can no longer change from 1 back to 0. Thus, UMATHT is used to compute the specific heat, the thermal conductivity, and the material state change and the phase change. The required inputs are listed below:
Thermal conductivity of the bulk material (𝜆𝑏𝑢𝑙𝑘 ) as a function of the temperature
Specific heat (𝑐) as a function of the temperature
Latent heat
Solidus temperature
Liquidus Temperature
Powder porosity
Neck ratio
Particle size
Stefan-Boltzmann constant
4.3.2. Emissivity No user subroutine to define the emissivity as a function of the temperature and of the material state is available in the Abaqus Standard. In order to model the heat radiation loss, a pseudo heat transfer coefficient ℎ𝑟𝑎𝑑 has been defined in equation (19) and introduced into a user subroutine called FILM. ℎ𝑟𝑎𝑑 = 𝜀(𝑇)𝜎(𝑇 2 + 𝑇𝑟 2 )(𝑇 + 𝑇𝑟 )
(19)
The pseudo heat transfer coefficient has units of measurement of the film coefficient and it includes emissivity as a function of the temperature. The required inputs are the emissivity 𝜀(𝑇) of the bulk material and the chamber temperature. 4.3.3. Density The user subroutine USDFLD is used to update the density value from the powder density to the bulk density after melting. This subroutine is linked to the subroutine GETVRM which reads the MAT_ID information and calculates the corresponding density.
4.4. How the FE model works 14
The flowchart in Figure 5a shows the model structure after the preheating phase. Initially, the MAT_ID and temperature values are read by the user subroutines UMATHT and FILM, at the start of the time increment for each material calculation point. The temperature dependent properties are computed according to the obtained MAT_ID value and the density is updated. The user subroutine DFLUX then reads the current total time and applies the heat flux at each flux integration point according to the motion law of the beam. The software then solves the heat transfer problem. When the temperature value at the material calculation points is higher than the solidus temperature and the material is powder, the MAT_ID value changes from 0 to 1. The MAT_ID variable is updated, if necessary, in the user subroutine UMATH. The cycle is repeated for each time increment until the simulation ends.
The flowchart in Figure 5b shows the calculation method for when the material change is not considered. The difference is significant. In fact, in this case, the properties of the material are only temperature dependent. Consequently, the powder thermo-physical properties are the ones that correspond to the properties of the liquid metal when the temperature is higher than the melting point. However, when the temperature drops, the properties again become those of the powder, even through the material is actually in bulk form. This approach leads to unreliable numerical results because the powder properties are significantly different from the corresponding bulk material, for example, those of the density and the thermal conductivity (Figure 6). It is possible to explain unreliable numerical results by considering the thermal diffusivity of the material. The thermal diffusivity value describes the heat transfer rate of a material, which is equal to the thermal conductivity divided by the density and the specific heat. If the thermal diffusivity is calculated from the relative properties of the bulk and of the powder material (Figure 6), the thermal diffusivity value of the bulk is higher than that of the powder material for each temperature value. For this reason, the heat moves rapidly through the bulk material, as it is able to conduct heat quickly, while the heat transfer through the powder material is slow, that is, close to zero, because the powder has low thermal diffusivity. Thus, if the material change procedure is not considered in the simulation after passage of the beam, the heat transfer is calculated between the melt pool and the already scanned material that has the powder properties, instead of the bulk material. In this case, the temperature gradient is high and the heat transfer is low. Thus, the dimensions of the melt pool, in terms of width, length and depth, are reduced compared to those in which the surrounding material has changed
15
properties in the bulk material. In fact, in the case in which the material change procedure is deactivated, it has been observed, for instance, that the width of the scan line is approximately equal to the beam diameter. Instead, as will be shown in the next section, the width of the scan line is several times greater than the beam diameter, because the previously scanned material, which is cooling as bulk material, transfers the heat rapidly because of the conduction effects, and thus increases the heat affected area.
5. Model Validation The proposed model has been validated against the results of experiments available in the literature in terms of width of the scan line. The validation was based on the experiments conducted by Qi, Yan, Lin, He and Zhang [25]. In their experiments, several stainless steel 316L single scan lines were fabricated, varying the beam current and the scan speed. Material Properties Density at 293 K, bulk [kg/mm3]
7.97
Apparent density of the powder [kg/mm3]*
5.42
Solidus Temperature, TS [K]
1658
Liquidus Temperature, TL [K]
1699
Latent heat, L [kJ/kg]
273
Average particle size of the powders [mm]
0.075
Table 1 and Figure 6 report the thermophysical properties of the material as implemented in the model. The bulk material properties were extracted from [12, 14], whereas the powder properties were computed according to the models outlined in Section 2. Table 2 summarises the details of the experiments, related to the simulation inputs introduced to reproduce the real experiments. In order to resemble the test setup, the entire workpiece was modelled as powder material. The dimensions of the model are 10 × 10 × 10 mm3 and the finer mesh size is 0.05 mm.
The width of the melt pool is measured by the MAT_ID state variable. As can be observed in Figure 7 and Figure 8, it is possible to distinguish in the simulation results: the bulk material (soft grey, MAT_ID 0), the powder material (intense grey, MAT_ID 1) and the boundary line (red line). The 16
boundary line separates the powder area from the area that has undergone the material change. Each single powder particle near the boundary line can physically be in any of these positions (Figure 7): (a) tangent to the boundary line in the powder material part; (b) melted and solidified in the bulk material part, and thus theoretically tangent to the boundary line in the bulk material; (c) on the boundary line. The powder particles in position (a) are sintered, due to heat conduction, and as a consequence, they are part of the solidified scan line [34]. Similarly, the powder particles in position (c) belong to the scan line, but they are melted. However, the scan line width obtained from the simulation results does not take into account the particles in position (a) and/or (c), therefore the maximum dispersion of the simulation results is considered equal to twice the particle size. Figure 9 compares numerical and experimental measurements of the width of each single track. The numerical measurement of the width of the scan line refers to the maximum width, Wnumeric max. This approach allows the width of the scan line to be defined, even when the boundary line is irregular, as shown in Figure 8. The interval bars show the deviations from the numerical measurements. The extremes of the interval bar are: Wmax, which is the maximum scan width when at least one powder particle is in position (a) for each boundary line, and the minimum width of scan line Wnumeric
min
measured
considering the minimum distance between the boundary lines. According to the experimental results, the width of the scan line, for each combination of process parameters, is greater than of the focus beam diameter. The comparison shows a good agreement between the numerical results and the experimental data. The deviation between the experimental data and the corresponding Wnumeric max values is less than 25%. The maximum deviation drops to 13% when the beam current is higher than 0.5 mA. Furthermore, the experimental width is always included within the error bar of the numerical results, except for when the beam current and the scan speed are equal to 0.5 mA and 2 mm/s, respectively. According to main effects plot, Figure 10 highlights the average outcome for each value of each scan speed and beam current for the experimental and numerical results. The model appears to be slightly sensitive to low values of the scan speed. However, except for low values of the scan speed, the numeric response follows the experimental results when the scan speed and the beam current change. In fact, the numeric width of the scan line decreases as the scanning speed increase and increases as the beam current increases. The deviation between the experimental and numerical results may be attributed to:
17
the effect of the Marangoni convection flows that has been neglected. This hypothesis affects the numerical results for low scan speeds; if the scan speed were increased, the deviation between the numerical and experimental results would be reduced;
a balling effect that occurs for a low heat flux value [14]; the scan line contour is more blurred for a low scan speed [25]. Thus, the sides of the track are not clearly recognizable and the measurement of the width becomes more difficult. The simulation results partially show this behaviour, because the width of the scan line is indented for a low beam current while a constant width is shown for a higher beam current (Figure 8);
the melting line that soaks the powders. When the scanning speed is high, the melting line, in fact, obtains less energy from the electron beam and soaks less powder [25]. This phenomenon could be analysed by introducing the mechanical behaviour of the powders, but this is beyond the scope of this work.
However, the deviation is generally acceptable for any processes simulation. The numerical outcomes show that the model is able to correctly represent the phenomena that occur during the EBM process. Hence, the proposed modelling of a uniform distribution seems to be able to correctly model the energy source, as is the model that is used to model the thermophysical properties of the powder material.
6. Conclusions In this work, a three-dimensional thermal FE model of the EBM process has been developed in order to implement a new modelling of the energy source and powder material properties. First, the impact of the electron beam on a preheated powder bed was studied. This analysis led to the energy source being modelled as a uniform distribution, and a corrective coefficient of the global flux being introduced as a function of the material, the acceleration voltage and the focus beam diameter. Because of the high temperature gradient and the presence of different material states, i.e. powder and bulk, the thermophysical properties were modelled as a function of the temperature and material state. The properties of the bulk material were extracted from technical literature, whereas the material properties of the gas-atomized powder were calculated assuming a BCC packing structure and introducing a formulation for the emissivity. Several user subroutines were used to solve the heat transfer problem in Abaqus/Standard. The user subroutine DFLUX was implemented to apply a heat 18
flux and to move the electron beam spot. The user subroutines UMATHT, USDFLD, and FILM were used to calculate the properties of the powder material. Moreover, conduction effects, radiation losses, melting and cooling phases with a latent heat effect, and a coupled phase change and material change were considered. Experimental data taken from the literature were used as a reference for validation purposes to test the response of the numerical model to conductive effects and to different widths of the scan lines. However, the deviation was found to be less than 25% when the beam current was equal to 0.5 mA, and dropped to 13% when the beam current values were higher. These deviations are acceptable for a general response of the simulation processes. This limited deviation may be ascribed to the fact that certain phenomena had not been taken into consideration. The trend of the numerical results was found to be similar to that of the experimental data. This behaviour indicates that the model is able to correctly represent the main phenomena that occur during the EBM process. Hence, the proposed modelling of the material properties as well as the uniform distribution of the energy source seem to constitute the correct approach to simulate the EBM process. The model could be used as a tool to identify the most appropriate set of parameters, in terms of beam characteristics, powder properties and process control, and to consider their interactions. Improvements to the existing model will include a temperature analysis, in the hatching mode, and the implementation of the mechanical properties of the material in order to have a complete coupled thermo-mechanical model for the evaluation of stresses and deformations induced by beam passage.
19
References [1] E. Atzeni, L. Iuliano, P. Minetola, A. Salmi, Redesign and cost estimation of rapid manufactured plastic parts, Rapid Prototyping J 16(5) (2010) 308-317. [2] N. Hopkinson, P. Dickens, Rapid prototyping for direct manufacture, Rapid Prototyping J 7(4) (2001) 197-202. [3] T. Wohlers, Additive Manufacturing State of the Industry (Wohlers Report 2014), Fort Collins, USA2014. [4] E. Atzeni, A. Salmi, Economics of additive manufacturing for end-usable metal parts, The International Journal of Advanced Manufacturing Technology 62(9-12) (2012) 1147-1155. [5] J. Milberg, M. Sigl, Electron beam sintering of metal powder, Prod. Eng. Res. Devel. 2(2) (2008) 117-122. [6] S. Biamino, A. Penna, U. Ackelid, S. Sabbadini, O. Tassa, P. Fino, M. Pavese, P. Gennaro, C. Badini, Electron beam melting of Ti-48Al-2Cr-2Nb alloy: Microstructure and mechanical properties investigation, Intermetallics 19(6) (2011) 776-781. [7] M. Chiumenti, M. Cervera, A. Salmi, C. Agelet de Saracibar, N. Dialami, K. Matsui, Finite element modeling of multipass welding and shaped metal deposition processes, Computer Methods in Applied Mechanics and Engineering 199(3740) (2010) 2343-2359. [8] C. Körner, E. Attar, P. Heinl, Mesoscopic simulation of selective beam melting processes, Journal of Materials Processing Technology 211(6) (2011) 978-987. [9] R. Ammer, M. Markl, U. Ljungblad, C. Körner, U. Rüde, Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method, Computers & Mathematics with Applications 67(2) (2014) 318-330. [10] R. Ammer, U. Rüde, M. Markl, C. Körner, Modeling of Thermodynamic Phenomena with Lattice Boltzmann Method for Additive Manufacturing Processes. [11] T. Scharowsky, A. Bauereiß, R. Singer, C. Körner, Observation and numerical simulation of melt pool dynamic and beam powder interaction during selective electron beam melting, Proceedings of the Solid Freeform Fabrication Symposium, 2012, pp. 815-820. [12] H. Qi, Y. Yan, F. Lin, R. Zhang, Scanning method of filling lines in electron beam selective melting, Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 221(12) (2007) 1685-1694. [13] M. Zäh, S. Lutzmann, M. Kahnert, F. Walchshäusl, Determination of Process Parameters for Electron Beam Sintering (EBS), Excerpt from the Proceedings of the COMSOL Conference Hannover, 2008. [14] M.F. Zäh, S. Lutzmann, Modelling and simulation of electron beam melting, Production Engineering 4(1) (2010) 1523. [15] N. Shen, K. Chou, Thermal Modeling of Electron Beam Additive Manufacturing Process: Powder Sintering Effects, ASME 2012 International Manufacturing Science and Engineering Conference collocated with the 40th North American Manufacturing Research Conference and in participation with the International Conference on Tribology Materials and Processing, American Society of Mechanical Engineers, 2012, pp. 287-295. [16] N. Shen, Y. Chou, Numerical thermal analysis in electron beam additive manufacturing with preheating effects, Proceedings of the 23rd Solid Freeform Fabrication Symposium, Austin, TX, 2012, pp. 774-784. [17] N.G. Shen, K. Chou, Thermal Modeling of Electron Beam Additive Manufacturing Process Powder Sintering Effects, Proceedings of the Asme International Manufacturing Science and Engineering Conference, 2012 (2012) 287-295. [18] N. Shen, K. Chou, Simulations of Thermo-Mechanical Characteristics in Electron Beam Additive Manufacturing, ASME 2012 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2012, pp. 67-74. [19] B. Cheng, S. Price, J. Lydon, K. Cooper, K. Chou, On Process Temperature in Powder-Bed Electron Beam Additive Manufacturing: Model Development and Validation, J Manuf Sci E-T Asme 136(6) (2014). [20] S. Price, B. Cheng, J. Lydon, K. Cooper, K. Chou, On Process Temperature in Powder-Bed Electron Beam Additive Manufacturing: Process Parameter Effects, J Manuf Sci E-T Asme 136(6) (2014). [21] B. Cheng, S. Price, X. Gong, J. Lydon, K. Cooper, K. Chou, Speed Function Effects in Electron Beam Additive Manufacturing, ASME 2014 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2014, pp. V02AT02A003-V02AT02A003. [22] S.S. Sih, J.W. Barlow, Emissivity of powder beds, Solid Freeform Fabrication Symposium Proceedings, Center for Materials Science and Engineering, Mechanical Engineering Department and Chemical Engineering Department, the University of Texas at Austin, 1995, pp. 7-9. [23] N.K. Tolochko, M.K. Arshinov, A.V. Gusarov, V.I. Titov, T. Laoui, L. Froyen, Mechanisms of selective laser sintering and heat transfer in Ti powder, Rapid Prototyping Journal 9(5) (2003) 314-326. [24] E. Soylemez, J.L. Beuth, K. Taminger, Controlling Melt Pool Dimensions over a Wide Range of Material Deposition Rates in Electron Beam Additive Manufacturing, Solid Freeform Fabrication Proceedings, 2010, pp. 571-582.
20
[25] H.B. Qi, Y.N. Yan, F. Lin, W. He, R.J. Zhang, Direct metal part forming of 316L stainless steel powder by electron beam selective melting, P I Mech Eng B-J Eng 220(11) (2006) 1845-1853. [26] M. Krishnan, E. Atzeni, R. Canali, F. Calignano, D. Manfredi, E.P. Ambrosio, L. Iuliano, On the effect of process parameters on properties of AlSi10Mg parts produced by DMLS, Rapid Prototyping J 20(6) (2014). [27] R.M. German, Particle Packing Characteristics, Metal Powder Industries Federation1989. [28] E. Attar, Simulation of selective electron beam melting processes, Dr.-Ing., University of Erlangen, Nuremberg, Germany (2011). [29] K. Alexander, B. Andreas, K. Carolin, Modelling of electron beam absorption in complex geometries, Journal of Physics D: Applied Physics 47(6) (2014) 065307. [30] T.R. Mahale, Electron Beam Melting of Advanced Materials and Structures, Industrial Engineering, North Carolina State University, 2009. [31] M. Agarwala, D. Bourell, J. Beaman, H. Marcus, J. Barlow, Direct selective laser sintering of metals, Rapid Prototyping Journal 1(1) (1995) 26-36. [32] A. Gusarov, I. Yadroitsev, P. Bertrand, I. Smurov, Heat transfer modelling and stability analysis of selective laser melting, Applied Surface Science 254(4) (2007) 975-979. [33] Abaqus, ABAQUS Documentation, Dassault Systèmes, Providence (USA), 2014. [34] O. Cansizoglu, O. Harrysson, D. Cormier, H. West, T. Mahale, Properties of Ti–6Al–4V non-stochastic lattice structures fabricated via electron beam melting, Materials Science and Engineering: A 492(1) (2008) 468-474.
21
Figure 1 – (a) the elementary model of a powder (BCC structure) before the preheating phase, (b) the elementary model after the preheating phase.
Figure 2 – Electron trajectories in a cross section of a 316L solid sample. The electron trajectory simulation has been run by means of CASINO software [27]. The red lines refer to back-scattered electron trajectories, while the blue lines refer to the paths of non-back-scattered electrons. The grey lines illustrate the distribution of the electrons along the radius and along the layer depth. An acceleration voltage of 60 kV and a beam diameter of 0.400 mm have been adopted. The effective diameter 𝐷𝐸 is larger than the nominal focus beam diameter.
22
Figure 3 – Heat flux distribution along the radius. The grey dashed line is the Gaussian distribution function that is used in literature [13-21] that will be adjusted by the efficiencies; the black dashed line is the Gaussian distribution corrected by means of 𝐷𝐸 ; the continuous black line is the approximation of the Gaussian distribution corrected by means of 𝐷𝐸 .
Figure 4 - Model configuration and boundary condition
23
Figure 5 – FE model structure (the subroutines UMATHT and DFLUX) with the material change procedure (a) and without the material change procedure (b)
24
Figure 6 – Specific heat of the bulk material [12]. The specific heat of the powder material can be considered the same as that of the bulk material [14, 15, 22, 23]. Thermal conductivity and specific heat of the bulk [12] and powder material according to the model presented in Section 2.1. Emissivity value for bulk material [14] and powder material according to the model presented in Section 2.2.
25
Figure 7 – Possible position of a single powder particle: (a) tangent to the boundary line in the powder material part; (b) melted and solidified in the bulk material part and thus theoretically tangent to the boundary line in the bulk material; (c) over the boundary line.
Figure 8 – Scan line for 0.5 mA of beam current and 3.3 mm/s of scan speed. The indented scan line can be observed. The numeric maximum width of scan line Wnumeric max is approximately 0.55 mm. The minimum width Wnumeric min is equal to 0.45 mm. Considering as the powder particles near the boundary line could be sintered, the maximum width of the scan line, Wmax, is equal to the sum of Wnumeric max and twice the particle diameter (75µm).
26
Figure 9 – Comparison between the experimental and numerical results
27
Figure 10 – Plot of the main effects Table 1 Material properties of 316L Stainless Steel [12, 25]. Material Properties Density at 293 K, bulk [kg/mm3]
7.97
Apparent density of the powder [kg/mm3]*
5.42
Solidus Temperature, TS [K]
1658
Liquidus Temperature, TL [K]
1699
Latent heat, L [kJ/kg]
273
Average particle size of the powders [mm] * value calculated according to Eq. 1
0.075
Table 2 Simulation inputs on the basis of experimental parameters [25]. Simulations parameters Beam diameter, 𝐷 [mm]
0.200
Accelerating voltage, 𝑈 [kV]
50
𝜂*
0.13
Beam current, 𝐼 [mA] Scan speed, 𝑥1̇ [mm/s] Surface flux, 𝑞 [W/mm2]*
0.5 2.3
3.3
0.8 5.0
2.3
3.3
96.32 154.11 * values calculated according to Section 4.2
28
1.0 5.0
2.3
3.3 192.64
5.0