Material properties for Ti6Al4V and main process parameters used in simulations. ^{a}Value for commercially pure titanium was used.
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 strengthtoweight 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 nearnet 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 lowyield 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 powderbased spray (e.g., thermal spray or cold spray), sintering (e.g., selective laser sintering), or fusionbased 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], fusionbased AM processes are more advantageous since they can produce 100% dense functional metal parts. This chapter will focus on fusionbased 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 tensiondominated freesurface 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 processproperty 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 processscale 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 buoyancydriven 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 tensiondriven freesurface flow determines the shape and smoothness of the clad. In this chapter, numerical modeling of transport phenomena in fusionbased 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 Ti6Al4V for both the substrate and powder. As the main challenges in modeling of fusionbased AM processes are related to melting/solidification phase change and freesurface 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
Fusionbased 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 socalled “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 solidliquid 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 phasefield [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], cellularautomaton [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.
To treat the effects of transport phenomena at the processscale (~ 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 twophase model [Beckermann & Viskanta, 1988; Ganesan & Poirier, 1990; Ni & Beckermann, 1991], in which the two phases are treated as separate and separate volumeaveraged 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 solidliquid 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, singlephase CFD procedures. In this study, the continuum model is adopted to develop the governing equations.
1.3. Modeling of freesurface flow
In fusionbased AM processes, the melt pool created by the energy source on the substrate is usually modelled as a freesurface 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 remeshing 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 fluidstructure interactions or small amplitude waves [Caboussat, 2005]. In Eulerian methods, the interface is moving within a fixed grid, and no remeshing is needed. The interface is determined from a field variable, for example, a volume fraction [DeBar, 1974; Hirt & Nichols, 1981; Noh & Woodward, 1976], a levelset [ Sethian, 1996, 1999] or a phasefield [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 fusionbased 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 VolumeOfFluid) [Hirt & Nichols, 1981] is probably the most widely used. It has been adopted by many inhouse codes and built into commercial codes (SOLAVOF [Nichols et al, 1980], NASAVOF2D [Torrey et al 1985], NASAVOF3D [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 liquidvolume 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 piecewiseconstant or “stairstepped” representation of the interface and advection schemes are at best firstorder accurate. These firstorder 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]. PLICVOF 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 lowReynoldsnumber
For the system of interest, the conservation equations are summarized as follows:
Mass conservation:
Momentum conservation:
Energy conservation:
In equations (1)(4), the subscripts
Here the subscripts
The volume fraction of liquid
The phase enthalpy for the solid and the liquid can be expressed as:
where
In Eqs. (2) and (3), the third terms on the righthand side are the drag interaction terms, and
2.2. Surface tension
The surface tension force,
where
The conventional approach when dealing with surface tension is to use finite difference schemes to apply a pressure jump at a freesurface 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 freesurface VOF data [Brackbill et al., 1992].
A wellknown problem with VOF (and other Eularian methods) modeling of surface tension is socalled “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 breakup [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 forcebalance 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 pressuregradient terms can be avoided. Within a VOF framework, such forcebalance 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 secondorder 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 ydirection component of the interface normal vector is larger than the xdirection component, interface curvature (in 2D) is given by
where
2.3. Tracking of the free surface
The free surface of the melt pool is tracked using the PLICVOF [Gueyffier et al., 1999; Scardovelli & Zaleski, 2000, 2003]. The Volume of Fluid function, F, satisfies the following conservation equation:
The PLICVOF 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:
where
2.4. Boundary conditions
Energy balance at the free surface satisfies the following equation:
where terms on the righthand side are laser irradiation, convective heat loss, radiation heat loss and evaporation heat loss, respectively.
On the bottom surface and side surfaces, boundary conditions are given by
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 fusionbased additive manufacturing processes, the material experiences a large variation from room temperature to above the melting temperature. For Ti6Al4V, 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^{1}kg^{1} at room temperature to 831 JK^{1}kg^{1} at liquidus temperature. Thermal conductivity varies from 7 to 33.4 Wm^{1}K^{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 thermophysical 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:
Eqs. (1)  (3) and the related boundary conditions are solved iteratively using a twostep projection method [Chorin, 1968] to obtain velocities and pressures. Thermophysical 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 BiCGSTAB (for BiConjugate Gradient Stabilized) method [Barrett et al., 1994] is used to solve the linear system of equations.
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 thermophysical properties are updated.
Equation (15) is solved using the PLICVOF [Gueyffier et al., 1999; Scardovelli & Zaleski, 2000, 2003] to obtain the updated free surface and geometry of the melt pool.
. 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}) 

[Mills, 2002] 
Liquid specific heat (J kg^{1} K^{1})  831.0  [Mills, 2002] 
Thermal conductivity (W m^{1} K^{1}) 

[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 kg1) 
9.83 10^{6}  [Mills, 2002] 
Dynamic viscosity (N m1 s1)  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 m1)  1.525 – 0.28×10^{3}(T – 1941K)^{a}  [Mills, 2002] 
Thermal expansion coefficient (K1) 
1.1 × 10^{5}  [Mills, 2002] 
Laser absorption coefficient  0.4  
Ambient temperature (K)  300  
Convective coefficient (W m2 K1) 
10 
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. Ti6AlV4 plates with a thickness of 0.25 inch were selected as substrates. Ti6AlV4 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, 5axis CNC machine, and monitoring subsystem. The laser system used was a Nuvonyx ISL1000M 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 Ti6AlV4, the system is covered in an environmental chamber to supply argon gas. The melt pool peak temperature is measured by a noncontact 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 realtime at 500 or 1000 Hz using a National Instruments realtime control system. A 420 mA signal is sent to the realtime 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 Zaxis 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 twodimensional 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.
The samples were crosssectioned using a WireEDM machine to measure dilution depth. An SEM (Scanning Electron Microscope) line trace was used to determine the dilution of the clad layer. The deposited Ti6Al4V is of Widmansttaten structure. The substrate has a rolled equiaxed 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 Ti6Al4V. The small quantity makes sure that it does not drastically change the deposit features of a 100% Ti6Al4V 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.
Good agreements between measured and simulated dilution depths can be found in Figs. 810. 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 lackoffusion 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 Ti6Al4V specimens were inspected at Quality Testing Services Co. using ultrasonic and radiographic inspections to determine the extent of lackoffusion 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 lackoffusion was detected in specimens deposited using 1,200 watts and higher energy levels. However, lackoffusion 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.
4. Conclusion
This chapter has outlined the approach for mathematical and numerical modeling of fusionbased 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 freesurface flow in the melt pool.
The laser deposition process is selected as an example of fusionbased additive manufacturing processes. The governing equations, auxiliary relationships, and boundary conditions for the solidification system and freesurface flow are presented. The main challenge for modeling of the surface tensiondominant 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 fusionbased 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 BoeingSt Louis and Dr. James Sears at South Dakota School of Mines and Technology are also acknowledged.
References
 1.
Aulisa E. Manservisi S. Scardovelli R. Zaleski S. 2007 Interface reconstruction with leastsquares fit and split advection in threedimensional Cartesian geometry ,  2.
Barrett R. Berry M. Chan T. F. Demmel J. Donato J. Dongarra J. Eijkhout V. Pozo R. Romine C. Van der Vorst H. 1994  3.
Beckermann C. Viskanta R. 1988 Doublediffusive convection during dendritic solidification of a binary mixture.  4.
Beckermann C. Diepers H. J. Steinbach I. Karma A. Tong X. 1999 Modeling Melt Convection in PhaseField Simulations of Solidification ,  5.
Benes J. 2007 Cool Tips for Cutting Titanium, In:  6.
Bennon W. D. Incropera F. P. 1987 A continuum model for momentum, heat and species transport in binary solidliquid phase change systemsI. Model formulation .  7.
Benson D. J. 2002 Volume of fluid interface reconstruction methods for multimaterial problems .  8.
Bhat M. S. Poirier D. R. Heinrich J. C. 1995 Permeability for cross flow through columnardendritic alloys .  9.
Boettinger W. J. Coriell S. R. Greer A. L. Karma A. Kurz W. Rappaz M. Trivedi R. 2000 Solidification microstructures: recent developments, future directions ,  10.
Boettinger W. J. Warren J. A. Beckermann C. Karma A. 2002 PhaseField Simulation of Solidification ,  11.
Boyer R. Welsch G. Collings E. W. 1994 Materials properties handbook: titanium alloys, ASM International,9780871704818 Materials Park, OH  12.
Brackbill J. U. Kothe D. B. Zemach C. 1992 A continuum method for modeling surface tension,  13.
Caboussat A. 2005 Numerical Simulation of TwoPhase Free Surface Flows .  14.
Carman P. C. 1937 Fluid flow through granular beds .  15.
Caginalp G. 1989 Stefan and HeleShaw type models as asymptotic limits of the phasefield equations .  16.
Choi M. Greif R. Salcudean M. 1987 A study of the heat transfer during arc welding with applications to pure metals or alloys and low or high boiling temperature materials ,  17.
Chorin A. J. 1968 Numerical solution of the NavierStokes equations.  18.
Chorin A. J. 1980 Flame advection and propagation algorithms,  19.
Clyne T. W. Kurz W. 1981 Solute redisribution during solidification with rapid solid state diffusion.  20.
Cummins S. J. Francois M. M. Kothe D. B. 2005 Estimating curvature from volume fractions .  21.
Debar R. 1974  22.
Drummond J. E. Tahir M. I. 1984 Laminar viscous flow through regular arrays of parallel solid cylinders .  23.
Eylon D. Froes F. H. 1984 Titanium NetShape Technologies.  24.
Fan Z. Sparks T. E. Liou F. Jambunathan A. Bao Y. Ruan J. Newkirk J. W. 2007 Numerical Simulation of the Evolution of Solidification Microstructure in Laser Deposition. Proceedings of the 18th Annual Solid Freeform Fabrication Symposium, Austin, TX, USA, August 2007  25.
Fan Z. Stroble J. K. Ruan J. Sparks T. E. Liou F. 2007 Numerical and Analytical Modeling of Laser Deposition with Preheating. Proceeding s of ASME 2007 International Manufacturing Science and Engineering Conference,0791842908 Georgia, USA, October, 2007,37 51  26.
Fan Z. 2007 Numerical and Experimental Study of Parameter Sensitivity and Dilution in Laser Deposition, Master Thesis, University of MissouriRolla, 2007  27.
Fan Z. Jambunathan A. Sparks T. E. Ruan J. Yang Y. Bao Y. Liou F. 2006 Numerical simulation and prediction of dilution during laser deposition.  28.
Ferziger J. H. Peric M. 2002  29.
Francois M. M. Cummins S. J. Dendy E. D. Kothe D. B. Sicilian J. M. Williams M. W. 2006 A balancedforce algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework .  30.
Frenk A. Vandyoussefi M. Wagniere J. Zryd A. Kurz W. 1997 Analysis of the lasercladding process for stellite on steel .  31.
Fuster D. Agbaglah G. Josserand C. Popinet S. Zaleski S. 2009 Numerical simulation of droplets, bubbles and waves: state of the art .  32.
Gandin Ch.A Rappaz M. 1994 A coupled finite elementcellular automaton model for the prediction of dendritic grain structures in solidification processes ,  33.
Ganesan S. Poirier D. R. 1990 Conservation of mass and momentum for the flow of interdendritic liquid during solidification ,  34.
Ganesan S. Chan C. L. Poirier D. R. 1992 Permeability for flow parallel to primary dendrite arms .  35.
Gerlach D. Tomar G. Biswas G. Durst F. 2006 Comparison of volumeoffluid methods for surface tensiondominant twophase flows .  36.
Gibou F. Fedkiw R. Caflisch R. Osher S. 2003 A Level Set Approach for the Numerical Simulation of Dendritic Growth .  37.
Grujicic M. Cao G. Figliola R. S. 2001 Computer simulations of the evolution of solidification microstructure in the LENS^{TM} rapid fabrication process,  38.
Gueyffier D. Li J. Nadim A. Scardovelli R. Zaleski S. 1999 Volume of Fluid interface tracking with smoothed surface stress methods for threedimensional flows.  39.
Hansbo P. 2000 A FreeLagrange Finite Element Method using SpaceTime Elements .  40.
Hills R. N. Loper D. E. Roberts P. H. 1983 A thermodynamically consistent model of a mushy zone .  41.
Hirt C. W. Nichols B. D. 1981 Volumeoffluid (VOF) for the dynamics of free boundaries. Journ al of Computational Physics,39 1 (January 1981),201 225 ,00219991  42.
Hirt C. W. Nichols B. D. 1988  43.
Idelsohn S. R. Storti M. A. Onate E. 2001 Lagrangian Formulations to Solve Free  44.
Surface Incompressible Inviscid Fluid Flows.  45.
Jacqmin D. 1999 Calculation of twophase NavierStokes flows using phasefield modeling .  46.
Jaluria Y. 2006 Numerical Modeling of Manufacturing Processes. In:  47.
Jones W. P. Launder B. E. 1973 The calculation of lowReynoldsnumber phenomena with a twoequation model of turbulence.  48.
Juric D. Tryggvason G. 1996 A FrontTracking Method for Dendritic Solidification ,  49.
Karma A. Rappel W. J. 1996 Phasefield method for computationally efficient modeling of solidification with arbitrary interface kinetics ,  50.
Karma A. Rappel W. J. 1998 Quantitative phasefield modeling of dendritic growth in two and three dimensions .  51.
Kim Y. T. Goldenfeld N. Dantzig J. 2000 Computation of dendritic microstructures using a level set method .  52.
Knoll D. A. Kothe D. B. Lally B. 1999 A new nonlinear solution method for phase change problems .  53.
Kobayashi R. 1993 Modeling and numerical simulations of dendritic crystal growth .  54.
Kobryn P. A. Ontko N. R. Perkins L. P. Tiley J. S. 2006 Additive Manufacturing of Aerospace Alloys for Aircraft Structures. In:  55.
Kothe D. B. Mjolsness R. C. Torrey M. D. 1991 RIPPLE: A Computer Program for Incompressible Flows with Free surfaces, Technical Report, LA12007 MS, Los Alamos National Lab  56.
Lafaurie B. Nardone C. Scardovelli R. Zaleski S. Zanetti G. 1994 Modelling Merging and Fragmentation in Multiphase Flows with SURFER .  57.
Liou F. Kinsella M. 2009 A Rapid Manufacturing Process for High Performance Precision Metal Parts .  58.
Liovic P. Francois M. Rudman M. Manasseh R. 2010 Efficient simulation of surface tensiondominated flows through enhanced interface geometry interrogation .  59.
Lips T. Fritsche B. 2005 A comparison of commonly used reentry analysis tools ,  60.
López J. Hernández J. 2010 On reducing interface curvature computation errors in the height function technique .  61.
Meier M. Yadigaroglu G. Smith B. L. 2002 A novel technique for including surface tension in PLICVOF methods .  62.
Mills K. C. 2002 Recommended Values of Thermophysical Properties for Selected Commercial Alloys , Woodhead Publishing Ltd,9781855735699 Cambridge  63.
Muller I. A. 1968 A Thermodynamic theory of mixtures of fluids,  64.
Ni J. Beckermann C. 1991 A volumeaveraged twophase model for transport phenomena during solidification ,  65.
Ni J. Incropera F. P. 1995 Extension of the Continuum Model for Transport Phenomena Occurring during Metal Alloy Solidification, Part I: The Conservation Equations.  66.
Nichols B. D. Hirt C. W. Hotchkiss R. S. 1980  67.
Noh W. F. Woodward P. R. 1976 SLIC (simple line interface calculation),  68.
Pavlyk V. Dilthey U. 2004 2004). Numerical Simulation of Solidification Structures During Fusion Welding. In:  69.
Pilliod Jr J. E. Puckett E. G. 2004 Secondorder accurate volumeoffluid algorithms for tracking material interfaces.  70.
Poirier D. R. 1987 Permeability for flow of interdendritic liquid in columnardendritic alloys.  71.
Prantil V. C. Dawson P. R. 1983 Application of a mixture theory to continuous casting. In:  72.
Prescott P. J. Incropera F. P. Bennon W. D. 1991 Modeling of Dendritic Solidification Systems: Reassessment of the Continuum Momentum Equation.  73.
Provatas N. Goldenfeld N. Dantzig J. 1998 Efficient Computation of Dendritic Microstructures using Adaptive Mesh Refinement.  74.
Ramaswany B. Kahawara M. 1987 Lagrangian finite element analysis applied to viscous free surface fluid flow.  75.
Rappaz M. Gandin Ch. A. 1993 Probabilistic modelling of microstructure formation in solidification processes.  76.
Renardy Y. Renardy M. 2002 PROST: A parabolic reconstruction of surface tension for the volumeoffluid method.  77.
Rider W. J. Kothe D. B. 1998 Reconstructing volume tracking.  78.
Rudman M. 1997 Volumetracking methods for interfacial flow calculations.  79.
Scardovelli R. Zaleski S. 1999 Direct numerical simulation of freesurface and interfacial flow,  80.
Scardovelli R. Zaleski S. 2000 Analytical relations connecting linear interfaces and volume fractions in rectangular grids,  81.
Scardovelli R. Zaleski S. 2003 Interface Reconstruction with LeastSquare Fit and Split EulerianLagrangian Advection.  82.
Scheil E. 1942  83.
Sethian J. A. 1996  84.
Sethian J. A. 1999  85.
Shirani E. Ashgriz N. Mostaghimi J. 2005 Interface pressure calculation based on conservation of momentum for front capturing methods.  86.
Stefanescu D. M. 2002  87.
Steinbach I. Pezzolla F. Nestler B. Seeszelberg M. Prieler R. Schmitz G. J. Rezende J. L. L. 1996 A phase field concept for multiphase systems.  88.
Sullivan .,Jr J. M.; Lynch ,D.R. O’Neill ,K. 1987 ). Finite element simulation of planar instabilities during solidification of an undercooled melt. Journal of Computational Physics, Vol. 69, No. 1, (March 1987), pp. 81111, ISSN 00219991  89.
Swaminathan C. R. Voller V. R. 1992 A general enthalpy method for modeling solidification processes.  90.
Takizawa A. Koshizuka S. Kondo S. 1992 Generalization of physical component boundary fitted coordinate method for the analysis of free surface flow.  91.
Tang H. Wrobel L. C. Fan Z. 2004 Tracking of immiscible interfaces in multiplematerial mixing processes,  92.
Torrey M. D. Cloutman L. D. Mjolsness R. C. Hirt C. W. 1985  93.
Torrey M. D. Mjolsness R. C. Stein L. R. 1987  94.
Tryggvason G. Bunner B. Esmaeeli A. Juric D. AlRawahi N. Tauber W. Han J. Nas S. Jan Y. J. 2001 A FrontTracking Method for the Computations of Multiphase Flow.  95.
Udaykumar H. S. Mittal R. Shyy W. 1999 Computation of SolidLiquid Phase Fronts in the Sharp Interface Limit on Fixed Grids.  96.
Udaykumar H. S. Marella S. Krishnan S. 2003 Sharpinterface simulation of dendritic growth with convection: benchmarks.  97.
Voller V. R. Prakash C. 1987 A Fixed Grid Numerical Modeling Methodology for ConvectionDiffusion Mushy Region Phase Change Problems.  98.
Voller V. R. Brent A. Prakash C. 1989 The Modeling of Heat, Mass and Solute Transport in Solidification Systems,  99.
Voller V. R. Mouchmov A. Cross M. 2004 An explicit scheme for coupling temperature and concentration fields in solidification models.  100.
Voller V. R. 2006 Numerical Methods for PhaseChange Problems. In:  101.
Warren J. A. Boettinger W. J. 1995 Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phasefield method.  102.
West R. 1985 On the permeability of the twophase zone during solidification of alloys.  103.
Wheeler A. A. Boettinger W. J. Mc Fadden G. B. 1992 Phasefield model for isothermal phase transitions in binary alloys.  104.
Youngs D. L. 1982 Timedependent multimaterial flow with large fluid distortion. In: Numerical Methods for Fluid Dynamics, Morton, K.W. & Baines, M. J. (Ed.),273 285 , Academic Press,9780125083607 UK  105.
Youngs D. L. 1984  106.
Yu K. O. Imam M. A. 2007 Development of Titanium Processing Technology in the USA.  107.
Zhu M. F. Lee S. Y. Hong C. P. 2004 Modified cellular automaton model for the prediction of dendritic growth with melt convection.