Numerical Modeling of the Additive Manufacturing (AM) Processes of Titanium Alloy

1. Introduction

It is easy to understand why industry and, especially, aerospace engineers love titanium. Titanium parts weigh roughly half as much as steel parts, but its strength is far greater than the strength of many alloy steels giving it an extremely high strength-to-weight ratio. Most titanium alloys are poor thermal conductors, thus heat generated during cutting does not dissipate through the part and machine structure, but concentrates in the cutting area. The high temperature generated during the cutting process also causes a work hardening phenomenon that affects the surface integrity of titanium, and could lead to geometric inaccuracies in the part and severe reduction in its fatigue strength [Benes, 2007]. On the contrary, additive manufacturing (AM) is an effective way to process titanium alloys as AM is principally thermal based, the effectiveness of AM processes depends on the material's thermal properties and its absorption of laser energy rather than on its mechanical properties. Therefore, brittle and hard materials can be processed easily if their thermal properties (e.g., conductivity, heat of fusion, etc.) are favorable, such as titanium. Cost effectiveness is also an important consideration for using additive manufacturing for titanium processing. Parts or products cast and/or machined from titanium and its alloys are very expensive, due to the processing difficulties and complexities during machining and casting. AM processes however, have been found to be very cost effective because they can produce near-net shape parts from these high performance metals with little or no machining [Liou & Kinsella, 2009]. In the aerospace industry, titanium and its alloys are used for many large structural components. When traditional machining/cast routines are adopted, conversion costs for these heavy section components can be prohibitive due to long lead time and low-yield material utilization [Eylon & Froes, 1984]. AM processes have the potential to shorten lead time and increase material utilization in these applications. The following sections 1.1, 1.2 and 1.3 summarize the fundamental knowledge for the modeling of additive manufacturing processes.

1.1. Additive manufacturing

Additive manufacturing can be achieved by powder-based spray (e.g., thermal spray or cold spray), sintering (e.g., selective laser sintering), or fusion-based processes (or direct metal deposition) which use a laser beam, an electron beam, a plasma beam, or an electric arc as an energy source and either metallic powder or wire as feedstock [Kobryn et al., 2006]. For the aerospace industry which is the biggest titanium market in the U.S. [Yu & Imam, 2007], fusion-based AM processes are more advantageous since they can produce 100% dense functional metal parts. This chapter will focus on fusion-based AM processes with application to titanium.

Numerical modeling and simulation is a very useful tool for assessing the impact of process parameters and predicting optimized conditions in AM processes. AM processes involve many process parameters, including total power and power intensity distribution of the energy source, travel speed, translation path, material feed rate and shielding gas pressure. These process parameters not only vary from part to part, but also frequently vary locally within a single part to attain the desired deposit shape [Kobryn et al., 2006]. Physical phenomena associated with AM processes are complex, including melting/solidification and vaporization phase changes, surface tension-dominated free-surface flow, heat and mass transfer, and moving heat source. The variable process parameters together with the interacting physical phenomena involved in AM complicate the development of process-property relationships and appropriate process control. Thus, an effective numerical modeling of the processing is very useful for assessing the impact of process parameters and predicting optimized conditions.

Currently process-scale modeling mainly addresses transport phenomena such as heat transfer and fluid dynamics, which are closely related to the mechanical properties of the final structure. For example, the buoyancy-driven flow due to temperature and species gradients in the melt pool strongly influences the microstructure and thus the mechanical properties of the final products. The surface tension-driven free-surface flow determines the shape and smoothness of the clad. In this chapter, numerical modeling of transport phenomena in fusion-based AM processes will be presented, using the laser metal deposition process as an example. Coaxial laser deposition systems with blown powder as shown in Fig. 1 are considered for simulations and experiments. The material studied is Ti-6Al-4V for both the substrate and powder. As the main challenges in modeling of fusion-based AM processes are related to melting/solidification phase change and free-surface flow in the melt pool, modeling approaches for these physical phenomena will be introduced in Sections 1.2 and 1.3.

1.2. Modeling of melting/solidification phase change

Fusion-based AM processes involve a melting/solidification phase change. Numerical modeling of the solidification of metal alloys is very challenging because a general solidification of metal alloys involves a so-called “mushy region” over which both solid and liquid coexist and the transport phenomena occur across a wide range of time and length scales [Voller, 2006]. A rapidly developing approach that tries to resolve the smallest scales of the solid-liquid interface can be thought of as direct microstructure simulation. In order to simulate the microstructure development directly, the evolution of the interface between different phases or different microstructure constituents has to be calculated, coupled with the physical fields such as temperature and concentration [Pavlyk & Dilthey, 2004]. To this approach belong phase-field [Beckermann et al., 1999; Boettinger et al., 2002; Caginalp, 1989; Karma & Rappel, 1996,1998; Kobayashi,1993; Provatas et al., 1998; Steinbach et al., 1996; Warren & Boettinger, 1995; Wheeler et al., 1992], cellular-automaton [Boettinger et al., 2000; Fan et al., 2007a; Gandin & Rappaz, 1994; Grujicic et al. 2001; Rappaz & Gandin, 1993; Zhu et al., 2004], front tracking [Juric & Tryggvason, 1996; Sullivan et al., 1987; Tryggvason et al., 2001], immersed boundary [Udaykumar et al., 1999, 2003] and level set [Gibou et al., 2003; Kim et al. 2000] methods. Due to the limits of current computing power, the above methods only apply to small domains on a continuum scale from about 0.1 μm to 10 mm.

Figure 1.

Schematics of a coaxial laser metal deposition system with powder injection

To treat the effects of transport phenomena at the process-scale (~ 1 m), a macroscopic model needs to be adopted, where a representative elementary volume (REV) is selected to include a representative and uniform sampling of the mushy region such that local scale solidification processes can be described by variables averaged over the REV [Voller et al., 2004]. Based on the REV concept, governing equations for the mass, momentum, energy and species conservation at the process scale are developed and solved. Two main approaches have been used for the derivation and solution of the macroscopic conservation equations. One approach is the two-phase model [Beckermann & Viskanta, 1988; Ganesan & Poirier, 1990; Ni & Beckermann, 1991], in which the two phases are treated as separate and separate volume-averaged conservation equations are derived for solid and liquid phases using a volume averaging technique. This approach gives the complete mathematical models for solidification developed today, which have the potential to build a strong linkage between physical phenomena occurring on macroscopic and microscopic scales [Ni & Incropera, 1995]. However, the numerical procedures of this model are fairly involved since two separate sets of conservation equations need to be solved and the interface between the two phases must be determined for each time step [Jaluria, 2006]. This places a great demand on computational capabilities. In addition, the lack of information about the microscopic configuration at the solid-liquid interface is still a serious obstacle in the implementation of this model for practical applications [Stefanescu, 2002]. An alternative approach to the development of macroscopic conservation equations is the continuum model [Bennon & Incropera, 1987; Hills et al., 1983; Prantil & Dawson, 1983; Prescott et al., 1991; Voller & Prakash, 1987; Voller et al., 1989]. This model uses the classical mixture theory [Muller, 1968] to develop a single set of mass, momentum, energy and species conservation equations, which concurrently apply to the solid, liquid and mushy regions. The numerical procedures for this model are much simpler since the same equations are employed over the entire computational domain, thereby facilitating use of standard, single-phase CFD procedures. In this study, the continuum model is adopted to develop the governing equations.

1.3. Modeling of free-surface flow

In fusion-based AM processes, the melt pool created by the energy source on the substrate is usually modelled as a free-surface flow, in which the pressure of the lighter fluid is not dependent on space, and viscous stresses in the lighter fluid is negligible. The techniques to find the shape of the free surface can be classified into two major groups: Lagrangian (or moving grid) methods and Eulerian (or fixed grid) methods. In Lagrangian methods [Hansbo, 2000; Idelsohn et al., 2001; Ramaswany& Kahawara, 1987; Takizawa et al., 1992], every point of the liquid domain is moved with the liquid velocity. A continuous re-meshing of the domain or part of it is required at each time step so as to follow the interface movement. A special procedure is needed to enforce volume conservation in the moving cells. All of this can lead to complex algorithms. They are mainly used if the deformation of the interface is small, for example, in fluid-structure interactions or small amplitude waves [Caboussat, 2005]. In Eulerian methods, the interface is moving within a fixed grid, and no re-meshing is needed. The interface is determined from a field variable, for example, a volume fraction [DeBar, 1974; Hirt & Nichols, 1981; Noh & Woodward, 1976], a level-set [ Sethian, 1996, 1999] or a phase-field [Boettinger et al., 2002; Jacqmin, 1999]. While Lagrangian techniques are superior for small deformations of the interfaces, Eulerian techniques are usually preferred for highly distorted, complex interfaces, which is the case for fusion-based additive manufacturing processes. For example, in AM processes with metallic powder as feedstock, powder injection causes intermittent mergers and breakups at the interface of the melt pool, which needs a robust Eulerian technique to handle.

Among the Eulerian methods, VOF (for Volume-Of-Fluid) [Hirt & Nichols, 1981] is probably the most widely used. It has been adopted by many in-house codes and built into commercial codes (SOLA-VOF [Nichols et al, 1980], NASA-VOF2D [Torrey et al 1985], NASA-VOF3D [Torrey et al 1987], RIPPLE [Kothe & Mjolsness 1991], and FLOW3D [Hirt & Nichols 1988], ANSYS Fluent, to name a few). In this method a scalar indicator function, F, is defined on the grid to indicate the liquid-volume fraction in each computational cell. Volume fraction values between zero and unity indicate the presence of the interface. The VOF method consists of an interface reconstruction algorithm and a volume fraction advection scheme. The features of these two steps are used to distinguish different VOF versions. For modeling of AM processes, an advantage of VOF is that it can be readily integrated with the techniques for simulation of the melting /solidification phase change. VOF methods have gone through a continuous process of development and improvement. Reviews of the historical development of VOF can be found in [Benson, 2002; Rider & Kothe, 1998; Rudman, 1997; Tang et al., 2004]. In earlier versions of VOF [Chorin, 1980; Debar, 1974; Hirt & Nichols, 1981; Noh & Woodward, 1976], reconstruction algorithms are based on a piecewise-constant or “stair-stepped” representation of the interface and advection schemes are at best first-order accurate. These first-order VOF methods are numerically unstable in the absence of surface tension, leading to the deterioration of the interface in the form of flotsam and jetsam [Scardovelli & Zaleski, 1999]. The current generation of VOF methods approximate the interface as a plane within a computational cell, and are commonly referred to as piecewise linear interface construction (PLIC) methods [Gueyffier et al., 1999; Rider & Kothe, 1998; Youngs, 1982, 1984]. PLIC-VOF is more accurate and avoids the numerical instability [Scardovelli & Zaleski, 1999].

2. Mathematical model

2.1. Governing equations

In this study the calculation domain for a laser deposition system includes the substrate, melt pool, remelted zone, deposited layer and part of the gas region, as shown in Fig.2. The continuum model Bennon & Incropera, 1987; Prescott et al., 1991 is adopted to derive the governing equations for melting and solidification with the mushy zone. Some important terms for the melt pool have been added in the momentum equations, including the buoyancy force term and surface tension force term, while some minor terms in the original derivation in [Prescott et al., 1991 have been neglected. The molten metal is assumed to be Newtonian fluid, and the melt pool is assumed to be an incompressible, laminar flow. The laminar flow assumption can be relaxed if turbulence is considered by an appropriate turbulence model, such as a low-Reynolds-number k-ε model [Jones & Launder, 1973]. The solid and liquid phases in the mushy zone are assumed to be in local thermal equilibrium.

Figure 2.

Schematic diagram of the calculation domain for laser metal deposition process

For the system of interest, the conservation equations are summarized as follows:

Mass conservation:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho V)=0$$ E1

Momentum conservation:

$$\frac{\partial }{\partial t}(\rho u)+\nabla \cdot (\rho Vu)=\nabla \cdot ({\mu }_l\frac{\rho }{{\rho }_l}\nabla u)-\frac{\partial p}{\partial x}-\frac{{\mu }_l}{K_x}\frac{\rho }{{\rho }_l}(u-u_s)+\rho g_x[1-\alpha (T-T_0)]+F_{Sx}$$ E2

$$\frac{\partial }{\partial t}(\rho v)+\nabla \cdot (\rho Vv)=\nabla \cdot ({\mu }_l\frac{\rho }{{\rho }_l}\nabla v)-\frac{\partial p}{\partial y}-\frac{{\mu }_l}{K_y}\frac{\rho }{{\rho }_l}(v-v_s)+\rho g_y[1-\alpha (T-T_0)]+F_{Sy}$$ E3

Energy conservation:

$$\frac{\partial }{\partial t}(\rho h)+\nabla \cdot (\rho Vh)=\nabla \cdot (k\nabla T)-\nabla \cdot [\rho (h_l-h)(V-V_s)]+S$$ E4

In equations (1)-(4), the subscripts s and l stand for solid and liquid phase, respectively. t, μ, and T are time, dynamics viscosity and temperature, respectively. u and v are x-direction and y-direction velocity components. The continuum density ρ, vector velocity V, enthalpy h, and thermal conductivity k are defined as follows:

$$\rho =g_s{\rho }_s+g_l{\rho }_l$$ E5

$$V=f_s{V}_s+f_l{V}_l$$ E6

$$h=f_s{h}_s+f_l{h}_l$$ E7

$$k=g_s{k}_s+g_l{k}_l$$ E8

Here the subscripts s and l stand for solid and liquid phase, respectively. f _s and f _l refer to mass fractions of solid and liquid phases, and g _s and g _l are volume fractions of solid and liquid phases. To calculate these four quantities, a general practice is that g _l (or g _s ) is calculated first and then the other three quantities are obtained according to the following relationships:

$$f_l=\frac{g_l{\rho }_l}{\rho }\cdot f_s=\frac{g_s{\rho }_s}{\rho }\cdot g_s+g_l=1\cdot f_s+f_l=1$$ E9

The volume fraction of liquid g _l can be found using different models, such as the level rule, the Scheil model [Scheil, 1942], or the Clyne and Kurz model [Clyne & Kurz, 1981]. For the target material Ti-6Al-4V, it is assumed that g _l is only dependent on temperature. The g _l (T) function is given by Swaminathan & Voller, 1992:

$$g_l=\{\begin{array}{c}0\text{if}T{T}_s^{}\\ \frac{T-T_s}{T_l-T_s}\text{if}{T}_s\le T\le T_l\\ 1\text{if}T{T}_l\end{array}$$ E10

The phase enthalpy for the solid and the liquid can be expressed as:

$$h_s={\displaystyle {\int }_0^T{c}_s(T)dT}$$ E11

$$h_l={\displaystyle {\int }_0^{T_s}{c}_s(T)dT}+{\displaystyle {\int }_{T_s}^{T_{}}{c}_l(T)dT+L_m}$$ E12

where L _m is the latent heat of melting. c _s and c _l are specific heat of solid and liquid phases.

In Eqs. (2) and (3), the third terms on the right-hand side are the drag interaction terms, and K _x and K _y are the permeability of the two-phase mushy zone in x- and y- directions, which can be calculated from various models Bhat et al., 1995; Carman, 1937; Drummond & Tahir, 1984; Ganesan et al., 1992; Poirier, 1987; West, 1985. Here the mushy zone is considered as rigid (i.e. a porous media). If the mushy zone is modeled as a slurry region, these two terms can be treated as in [Ni & Incropera, 1995]. In Eqs. (2) and (3), the fourth terms on the right-hand side are the buoyancy force components due to temperature gradients. Here Boussinesq approximation is applied. α is the thermal expansion coefficient. The fifth terms on the right-hand side of Eqs. (2) and (3) are surface tension force components, which will be described in Section 2.2 below. The term S in Eq. (4) is the heat source.

2.2. Surface tension

The surface tension force, F _S , is given by:

$$F_S=\gamma \kappa \widehat{n}+{\nabla }_S\gamma$$ E13

where $\gamma$ is surface tension coefficient, κ the curvature of the interface, $\widehat{n}$ the unit normal to the local surface, and ${\nabla }_S$ the surface gradient operator. The term $\gamma \kappa \widehat{n}$ is the normal component of the surface tension force. The term ${\nabla }_S\gamma$ represents the Marangoni effect caused by spatial variations in the surface tension coefficient along the interface due to temperature and/or species gradients. It causes the fluid flow from regions of lower to higher surface tension coefficient.

The conventional approach when dealing with surface tension is to use finite difference schemes to apply a pressure jump at a free-surface discontinuity. More recently, a general practice is to model surface tension as a volume force using a continuum model, either the Continuum Surface Force (CSF) model [Brackbill et al., 1992] or the Continuum Surface Stress (CSS) model [Lafaurie et al., 1994]. The volume force acts everywhere within a finite transition region between the two phases. In this study, the CSF model is adopted, which has been shown to make more accurate use of the free-surface VOF data [Brackbill et al., 1992].

A well-known problem with VOF (and other Eularian methods) modeling of surface tension is so-called “parasitic currents” or “spurious currents”, which is a flow induced solely by the discretization and by a lack of convergence with mesh refinement. Under some circumstances, this artificial flow can be strong enough to dominate the solution, and the resulting strong vortices at the interface may lead to catastrophic instability of the interface and may even break-up [Fuster et al., 2009, Gerlach et al. 2006]. Two measures can be taken to relieve or even resolve this problem. One measure is to use a force-balance flow algorithm in which the CSF model is applied in a way that is consistent with the calculation of the pressure gradient field. Thus, imbalance between discrete surface tension and pressure-gradient terms can be avoided. Within a VOF framework, such force-balance flow algorithms can be found in [Francois et al., 2006; Y.Renardy & M. Renardy, 2002; Shirani et al., 2005]. In this study, the algorithm in [Shirani et al., 2005] is followed. The other measure is to get an accurate calculation of surface tension by accurately calculating interface normals and curvatures from volume fractions. For this purpose, many methods have been developed, such as those in [Cummins et al, 2005; Francois et al., 2006; López & Hernández, 2010; Meier et al., 2002; Pilliod Jr. & Puckett, 2004; Y.Renardy & M. Renardy, 2002]. The method we use here is the height function (HF) technique, which has been shown to be second-order accurate, and superior to those based on kernel derivatives of volume fractions or RDF distributions [Cummins et al, 2005; Francois et al., 2006; Liovic et al., 2010]. Specifically, we adopt the HF technique in [López & Hernández, 2010] that has many improvements over earlier versions (such as that in [Torrey et al., 1985]) of HF, including using an error correction procedure to minimize estimation error. Within the HF framework, suppose the absolute value of the y-direction component of the interface normal vector is larger than the x-direction component, interface curvature (in 2D) is given by

$$\kappa =\frac{H_{xx}}{{(1+H_x{}^2)}^{3/2}}$$ E14

where H is the height function, H _x and H _xx are first-order and second-order derivatives of H, respectively. H _x and H _xx are obtained by using a finite difference formula. Interface normals are calculated based on the Least-Squares Fit method from [Aulisa et al., 2007].

2.3. Tracking of the free surface

The free surface of the melt pool is tracked using the PLIC-VOF [Gueyffier et al., 1999; Scardovelli & Zaleski, 2000, 2003]. The Volume of Fluid function, F, satisfies the following conservation equation:

$$\frac{\partial F}{\partial t}+(V\cdot \nabla )F=0$$ E15

The PLIC-VOF method consists of two steps: interface reconstruction and interface advection. In 2D calculation, a reconstructed planar surface becomes a straight line which satisfies the following equation:

$$n_{x}x+n_{y}y=d$$ E16

where n _x and n _y are x and y components of the interface normal vector. d is a parameter related to the distance between the line and the coordinate origin of the reference cell. In the interface reconstruction step, n _x and n _y of each cell are calculated based on volume fraction data, using the Least-Squares Fit method from [Aulisa et al., 2007]. Then the parameter d is determined to match the given volume fraction. Finally given the velocity field, the reconstructed interface is advected according to the combined Eulerian-Lagrangian scheme in [Aulisa et al., 2007].

2.4. Boundary conditions

Energy balance at the free surface satisfies the following equation:

$$k\frac{\partial T}{\partial n}=\frac{\eta (P_{laser}-P_{atten})}{\pi R^2}-h_c(T-T_{\infty })-\epsilon \sigma (T^4-T_{\infty }^4)-{\dot{m}}_e^{}{L}_v$$ E17

(17)

where terms on the right-hand side are laser irradiation, convective heat loss, radiation heat loss and evaporation heat loss, respectively. P _laser is the power of laser beam, P _atten the power attenuated by the powder cloud, R the radius of laser beam spot, the laser absorption coefficient, $${\dot{m}}_e^{}$$ the evaporation mass flux, L _v the latent heat of evaporation, h _c the heat convective coefficient, ε emissivity, σ the Stefan-Boltzmann constant, and $n$ the normal vector at the local interface. $${\dot{m}}_e^{}$$ can be evaluated according to the “overall evaporation model” in [Choi et al., 1987], and P_atten can be calculated according to Frenk et al., 1997 with a minor modification.

On the bottom surface and side surfaces, boundary conditions are given by

$$k\frac{\partial T}{\partial n}+h_c(T-T_{\infty })=0$$ E18

(

$$V=0$$ E19

Note that the radiation heat loss at these surfaces is neglected due to the fact that the temperature differences at these surfaces are not large.

2.5. Numerical Implementation

Finite difference and finite volume methods are used for spatial discretization of the governing equations. Staggered grids are employed where the temperatures, pressures and VOF function are located at the cell center and the velocities at the walls. In the numerical implementation, material properties play an important role. The material properties are generally dependent on temperature, concentration, and pressure. For fusion-based additive manufacturing processes, the material experiences a large variation from room temperature to above the melting temperature. For Ti-6Al-4V, many material properties experience large variations over this wide temperature range, as shown in Table 1. For example, the value of specific heat varies from 546 JK^-1kg^-1 at room temperature to 831 JK^-1kg^-1 at liquidus temperature. Thermal conductivity varies from 7 to 33.4 Wm^-1K^-1 over the same temperature range. Thus, the temperature dependence of the properties dominates, which necessitates a coupling of the momentum equations with the energy equation and gives rise to strong nonlinearity in the conservation equations.

The variable properties have two effects on the numerical solution procedure [Ferziger & Peric, 2002]. First, although an incompressible flow assumption is made, the thermo-physical properties need to be kept inside the differential operators. Thus, solution methods for incompressible flow can be used. Second, the momentum and energy conservation equations have to be solved in a coupled way. In this study, the coupling between momentum and energy equations is achieved by the following iterative scheme:

1. Eqs. (1) - (3) and the related boundary conditions are solved iteratively using a two-step projection method [Chorin, 1968] to obtain velocities and pressures. Thermo-physical properties used in this step are computed from the old temperature field. At each time step, the discretized momentum equations calculate new velocities in terms of an estimated pressure field. Then the pressure field is iteratively adjusted and velocity changes induced by each pressure correction are added to the previous velocities. This iterative process is repeated until the continuity equation is satisfied under an imposed tolerance by the newly computed velocities. This imposes a requirement for solving a linear system of equations. The preconditioned Bi-CGSTAB (for Bi-Conjugate Gradient Stabilized) method [Barrett et al., 1994] is used to solve the linear system of equations.

2. Eq. (4) is solved by a method [Knoll et al., 1999] based on a finite volume discretization of the enthalpy formulation of Eq. (4). The finite volume approach ensures that the numerical scheme is locally and globally conservative, while the enthalpy formulation can treat phase change in a straightforward and unified manner. Once new temperature field is obtained, the thermo-physical properties are updated.

3. Equation (15) is solved using the PLIC-VOF [Gueyffier et al., 1999; Scardovelli & Zaleski, 2000, 2003] to obtain the updated free surface and geometry of the melt pool.

4. . Advance to the next time step and back to step 1 until the desired process time is reached.

Physical Properties	Value	Reference
Liquidus temperature (K)	1923.0	[Mills, 2002]
Solidus temperature (K)	1877.0	[Boyer et al., 1994]
Evaporation temperature (K)	3533.0	[Boyer et al., 1994]
Solid specific heat (J kg^-1 K^-1)	$$\{\begin{array}{l}483.04+0.215TT\le 1268K\hfill \\ 412.7+0.1801T1268T\le 1923\hfill \end{array}$$	[Mills, 2002]
Liquid specific heat (J kg^-1 K^-1)	831.0	[Mills, 2002]
Thermal conductivity; (W m^-1 K^-1)	$$\{\begin{array}{l}1.2595+0.0157TT\le 1268K\hfill \\ 3.5127+0.0127T1268T\le 1923\hfill \\ -12.752+0.024TT1923\hfill \end{array}$$	[Mills, 2002]
Solid density (kg m^-3)	4420 – 0.154 (T – 298 K)	[Mills, 2002]
Liquid density (kg m^-3)	3920 – 0.68 (T – 1923 K)	[Mills, 2002]
Latent heat of fusion (J kg^-1)	2.86 10^5	[Mills, 2002]
Latent heat of evaporation; (J kg-1)	9.83 10^6	[Mills, 2002]
Dynamic viscosity (N m-1 s-1)	3.25 10^-3 (1923K) 3.03 10^-3 (1973K); 2.66 10^-3 (2073K) 2.36 10^-3 (2173K)	[Mills, 2002]
Radiation emissivity	0.1536+1.837710^-4 (T -300.0 K)	[Lips&Fritsche, 2005]
Surface tension (N m-1)	1.525 – 0.28×10^-3(T – 1941K)^a	[Mills, 2002]
Thermal expansion coefficient; (K-1)	1.1 × 10^-5	[Mills, 2002]
Laser absorption coefficient	0.4
Ambient temperature (K)	300
Convective coefficient; (W m-2 K-1)	10

Table 1.

Material properties for Ti-6Al-4V and main process parameters used in simulations. ^aValue for commercially pure titanium was used.

The time step is taken at the level of 10^-6 s initially and adapted subsequently according to the convergence and stability requirements of the Courant–Friedrichs–Lewy (CFL) condition, the explicit differencing of the Newtonian viscous stress tensor, and the explicit treatment of the surface tension force.

3. Simulation results and model validation

The parameters for the simulation were chosen based on the capability of our experimental facilities to compare the simulation results with the experimental measurements. A diode laser deposition system (the LAMP system of Missouri S&T) and a YAG laser deposition system at South Dakota School of Mines and Technology (SDSMT) were used for simulations and experiments. Ti-6Al-V4 plates with a thickness of 0.25 inch were selected as substrates. Ti-6Al-V4 powder particles with a diameter from 40 to 140 m were used as deposit material. Fig. 3 shows the typical simulation results for temperature, velocity and VOF function.

The numerical model was validated from different aspects. First, it was validated in terms of melt pool peak temperature and melt pool length. The experiments were performed on the LAMP system as shown in Fig. 4. The system consists of a diode laser, powder delivery unit, 5-axis CNC machine, and monitoring subsystem. The laser system used was a Nuvonyx ISL-1000M Diode Laser that is rated for 1 kW of output power. The laser emits at 808 nm and operates in the continuous wave (CW) mode. The laser spot size is 2.5 mm. To protect oxidization of Ti-6Al-V4, the system is covered in an environmental chamber to supply argon gas. The melt pool peak temperature is measured by a non-contact optical pyrometer that is designed for rough conditions, such as high ambient temperatures or electromagnetic interferences. A laser sight within the pyrometer allows for perfect alignment and focal length positioning; the spot size is 2.6 mm which encompasses the melt pool. The pyrometer senses the maximum temperature between 400 and 2500 (degrees C) and correlates the emissivity of the object to the resulting measurement. Temperature measurements are taken in real-time at 500 or 1000 Hz using a National Instruments real-time control system. A 4-20 mA signal is sent to the real-time system which is converted to degrees Celsius, displayed to the user and simultaneously recorded to be analyzed at a later date. Due to the collimator, the pyrometer is mounted to the Z-axis of the CNC at 42 (degrees) and is aligned with the center of the nozzle. Temperature measurements recorded the rise and steady state temperatures and the cooling rates of the melt pool. A complementary metal oxide semiconductor (CMOS) camera was installed right above the nozzle head for a better view in dynamically acquiring the melt pool image. The melt pool dimensions can be calculated from the image by the image process software.

Fig. 5 and Fig. 6 show the measured and predicted melt pool peak temperatures at different laser power levels and at different travel speeds, respectively. It can be seen from the plot that the general trend between simulation and experiment is consistent. At different power intensity level, there is a different error from 10 K (about 0.5%) to 121 K (about 5%). Fig. 7 shows measured and predicted melt pool length at different laser power levels. The biggest disagreement between measured and simulated values is about 7%. It can be seen that the differences between measured and predicted values at higher power intensities (higher power levels or slower travel speeds) are generally bigger than those at lower power intensities. This can be explained by the two-dimensional nature of the numerical model. A 2D model does not consider the heat transfer in the third direction. At a higher power level, heat transfer in the third dimension is more significant.

Figure 3.

Simulation results of laser deposition of Ti-6Al-4V

Figure 4.

Schematic of experimental setup

Figure 5.

Melt pool peak temperature comparison between simulation and experiment at different laser power levels (powder mass flow rate = 4.68g/min., travel speed = 20 inch/min.)

Figure 6.

Melt pool peak temperature comparison between simulation and experiment at different travel speeds (powder mass flow rate = 4.68g/min., laser power = 700 W)

Figure 7.

Melt pool length comparisons between simulation and experiment at different power levels (powder mass flow rate = 4. 68 g/min., travel speed = 20 inch/min.)

The samples were cross-sectioned using a Wire-EDM machine to measure dilution depth. An SEM (Scanning Electron Microscope) line trace was used to determine the dilution of the clad layer. The deposited Ti-6Al-4V is of Widmansttaten structure. The substrate has a rolled equi-axed alpha plus beta structure. Even though these two structures are are easily distinguishable, the HAZ is large and has a martensitic structure that can be associated with it. Hence a small quantity of tool steel in the order of 5% was mixed with Ti-6Al-4V. The small quantity makes sure that it does not drastically change the deposit features of a 100% Ti-6Al-4V deposit. At the same time, the presence of Cr in tool steel makes it easily identifiable by means of EDS scans using SEM. Simulation and experimental results of dilution depth are shown in Figs. 8 - 10.

Figure 8.

Comparison of dilution depth between simulation and experiment at different power levels (powder mass flow rate = 4.68 g/min., travel speed = 20 inch/min.)

Figure 9.

Comparison of dilution depth between simulation and experiment at different travel speeds and different laser power levels (powder mass flow rate = 4.68 g/min.)

Figure 10.

Comparison of dilution depth between simulation and experiment at different powder mass flow rates and different laser power levels (travel speed = 20 inch/min.)

Good agreements between measured and simulated dilution depths can be found in Figs. 8-10. The differences are from about 4.8% to 15.1%. It can be seen that an increase in the laser power will increase the dilution depth. An increase in the laser travel speed will decrease the dilution depth. It is clear that the dilution depth has a linear dependence on the laser power and the laser travel speed. This is easy to understand. As the laser power increases, more power is available for melting the substrate. As travel speed decreases, the laser material interaction time is extended. From Fig. 10, it can be seen that an increase in powder mass flow rate will decrease the dilution depth. But this effect is more significant at a higher level of laser power. It is likely that at a lower level of laser power, a significant portion of laser energy is consumed to melt the powder. Hence the energy available is barely enough to melt the substrate. Detailed discussion can be found in [Fan et al, 2006, 2007b; Fan, 2007].

Finally, the numerical model was validated in terms of its capability for predicting the lack-of-fusion defect. The test was performed using the YAG laser deposition system at South Dakota School of Mines and Technology (SDSMT). The simulation model determined that 1,200 watts would be the nominal energy level for the test. This means that based on the model, lack of fusion should occur when the laser power is below 1200 W. In accordance with the test matrix, seven energy levels were tested: nominal, nominal ± 10%, nominal ± 20%, and nominal ± 30%. Based on the predicted nominal value of 1,200 watts, the seven energy levels in the test matrix are 840, 960, 1080, 1200, 1320, 1440, and 1540 watts. The deposited Ti-6Al-4V specimens were inspected at Quality Testing Services Co. using ultrasonic and radiographic inspections to determine the extent of lack-of-fusion in the specimens. The determination of whether or not there exists lack of fusion in a deposited specimen can be explained using Figs. 11 - 13. First a substrate without deposit on it was inspected as shown in Fig. 11. Notice that the distance between two peaks are the thickness of the substrate. Then laser deposited specimens were inspected. If there is lack of fusion in a deposited specimen, some form of peaks can be found between the two high peaks in the ultrasonic graph, the distance of which is the height of the deposition and the thickness of the substrate. Fig. 12 shows an ultrasonic graph of a deposited specimen with a very good deposition. The ultrasonic result indicates there is not lack of fusion occurring between layers and the interface. The distance between two peaks is the height of the deposition and the thickness of the substrate. For the deposition as shown in Fig. 13, the lack of fusion occurs as the small peak (in circle) appears between two high peaks. The results revealed that no lack-of-fusion was detected in specimens deposited using 1,200 watts and higher energy levels. However, lack-of-fusion was detected in specimens deposited from lower energy levels (minus 10% up to minus 30% of 1,200 watts.). The test results validated the simulation model.

Figure 11.

Ultrasonic graph of a Ti-6Al-4V substrate

Figure 12.

Ultrasonic graph of a laser deposited Ti-6Al-4V specimen without lack of fusion

Figure 13.

Ultrasonic graph of a laser deposited Ti-6Al-4V specimen with lack of fusion

4. Conclusion

This chapter has outlined the approach for mathematical and numerical modeling of fusion-based additive manufacturing of titanium. The emphasis is put on modeling of transport phenomena associated with the process, including heat transfer and fluid flow dynamics. Of particular interest are the modeling approaches for solidification and free surface flow with surface tension. The advantages and disadvantages of the main modeling approaches are briefly discussed. Based on the comparisons, the continuum model is adopted for modeling of melting/solidification phase change, and the VOF method for modeling of free-surface flow in the melt pool.

The laser deposition process is selected as an example of fusion-based additive manufacturing processes. The governing equations, auxiliary relationships, and boundary conditions for the solidification system and free-surface flow are presented. The main challenge for modeling of the surface tension-dominant free surface flow is discussed and the measures to overcome the challenge are given. The numerical implementation procedures are outlined, with a focus on the effects of variable material property on the discretization schemes and solution algorithms. Finally the simulation results are presented and compared with experimental measurements. A good agreement has been obtained and thus the numerical model is validated. The modeling approach can be applied to other fusion-based manufacturing processes, such as casting and welding.

Acknowledgments

This research was partially supported by the National Aeronautics and Space Administration Grant Number NNX11AI73A, the grant from the U.S. Air Force Research Laboratory, and Missouri S&T’s Intelligent Systems Center and Manufacturing Engineering program. Their support is greatly appreciated. The help from Dr. Kevin Slattery and Mr. Hsin Nan Chou at Boeing-St Louis and Dr. James Sears at South Dakota School of Mines and Technology are also acknowledged.