Numerical Modeling of Metal-Based Additive Manufacturing Using Level Set Methods

The additive manufacturing technology produces parts by adding material in layers [1]. Compared with conventional manufacturing methods which cut off the excess material, additive manufacturing uses an energy source to melt the material and deposit it layer by layer [2]. AM allows fabricating products with great freedom in the shape, configuration, and function.

Many materials can be included in the AM process, such as metals, ceramics, or composites. In this paper, we focus on metal-based AM process since metal products are extensively employed in the industry. Metal-based AM products can be obtained with either a direct or an indirect procedure [3]. In the direct procedure, the metal particles are melted in the AM process. Selective laser melting (SLM) and electron beam melting (EBM) are two standard direct methods. In SLM or EBM, the powder is spread in a layer and melted by a laser or electron beam. When the temperature drops, the molten material solidifies and forms a layer of the part. This procedure is repeated with each layer adhered to the last one until the part is completed. In the indirect way, the layers of the products are bonded by a binder. The Binder Jetting-based AM is one of the commonly used indirect methods. The metal particles are first sintered, and after infiltration, the metal particles are bonded together. Though AM techniques have been widely developed these years, the process is not completely explored, which lays the groundwork for the monitoring and control of AM. As in Fig. 1 [4], both the direct and the indirect ways involve complex multiphysics phenomena, for example, melting and solidification, fluid flow, heat transfer, vaporization, or radiation. Also, in the SLM or EBM processes, the wetting mechanism should be studied as well [5]. Those physical phenomena are closely tied up with the AM process parameters which directly affect the quality of the products [6]. However, those critical parameters, such as laser power, scanning speed, scanning space, particle size or packing density, are mostly determined empirically. Those methods resting on empirical evidence are often unreliable and hard to manage when the manufacturing process varies.

Nowadays, predictive science-based simulation technologies have arisen as an effective tool to investigate the impact of the AM process parameters. Physical phenomena, such as the heat transfer in solid, heat absorption from the heat source, or heat radiation between solid with the surrounding gas, have been successfully simulated [7,8]. Fan and Liou [9] simulated the AM process of titanium alloy and compared the simulation results with experiment, which showed a good agreement. However, some phenomena are not well captured in current computational models due to the complexity of physics, for instance, the Marangoni convection due to surface tension gradients, the motion of liquid, and the radiation losses at the fluid/gas interface.

The objective of this paper is to simulate the metal-based AM process at the mesoscale using level set methods and multiphysics finite element simulation. In this research, the simulation model is built for the SLM process, where the material is spread in a layer and then melted by a laser beam. This research focuses on titanium alloy powders with a size between 10 $μm$ and 150 $μm$⁠. The heat transfer phenomena are modeled by solving the partial differential equations of physics transportation. By applying the level set method, the complex changes of particle geometries and the liquid/gas interface can be captured at high fidelity.

This paper is organized as follows: Section 2 describes the mathematical method and simulation techniques for modeling the metal-based AM process. In Sec. 3, two numerical examples simulating the power-based SLM melting process are presented. The computational results and future work are discussed in Sec. 4.

As introduced in Sec. 1, AM contains many physical processes. To reproduce the behavior of the AM system with a computer model, we first need set up the physical models before carrying out the numerical simulation. In this section, two principal physical models are introduced.

One is modeling the phase transition process. One way to model this is to simulate the solid–liquid interface at small scale. This approach is straightforward, but demands a lot of computing resources, so it only applies to small domains (0.1 μm to 10 mm) [9]. For the macroscopic transport problem, the representative elementary volume (REV) is introduced in dealing with the size difference between the fully solid or liquid region and the mushy zone. This model selects a zone to include a representative and uniform sampling of the mushy region, so the local scale solidification can be characterized by variables averaged over the REV [10]. Based on the concept, the transportation equations for the conservation in the process can be developed and solved.

The molten metal flow in the AM processes is a free-surface flow, which means the viscous stresses in the molten metal flow are assumed to be zero. Both the Lagrangian method and the Eulerian method have been applied to model the shape of the free-surface flow [9]. However, the Lagrangian method becomes more difficult in dealing with fluid problems with topological changes. The Eulerian method is more suitable to handle moving boundary problems. For instance, the Eulerian method can be used in modeling the powder-based AM process for tracking the mergers and splits at the interface of the melt pool [9].

where $u$ is the velocity at the interface. The application of level set methods to fluid dynamics problem has been carried out, for instance, Chung and Das [12] used the level set method to solve the Stefan problem and successfully tracked the interface between the liquid and solid phases. In this paper, the free-surface motion of the melt pools can be captured by solving the Hamilton–Jacobi partial differential equation [13,14]. The later motion is determined by the velocity field $u$⁠, which rests on the external physical or geometric information of the interface. Once the initial level set $ϕ$ (x, y, t = 0) is identified, a proper velocity field needs to be calculated to update the level set function and drive the boundary.

The computational domain for the SLM process modeling comprises the melt pools, deposited metal powder, and the surrounding air, as shown in Fig. 1. The continuum model, which describes the liquid in collective variables, is employed to estimate the properties of the solid–liquid coexisting area. The molten metal is assumed to be a Newtonian fluid, with a linear relationship between the rate of its deformation and the applied shear stress. The melt pool is considered as an incompressible laminar flow, the density of which is constant with the changes in pressure and temperature.

The conservation equations for the whole computational domain are formulated as below.

where $ρ,t,V,μ,k,T,p,A, and h$ are density, time, velocity, molten fluid dynamic viscosity, heat conductivity, temperature, pressure, permeability, and enthalpy, respectively. $Fst$ represents the capillary forces and $S$ represents the source term which related to the gravity, and $Fc$ refers to the thermocapillary force. The subscription s represents solid, and l represents liquid.

where $K0$ is the permeability coefficient, $gl$ is the mass fraction of liquid, and $ϵ0$ is a small constant. From the equation, the permeability A vanishes in the liquid region and approaches to infinity in the solid region. By adding this damping term, the momentum equation is valid for both solid and liquid phases.

In this study, the level set method tracks the interface on a fixed grid. The source terms in the momentum equations are applied as interfacial forces in the fluid flow, such as capillary force, buoyancy force, or thermocapillary force [16].

where $α,g,and Tref$ are the thermal expansion coefficient, the acceleration of gravity, and the reference temperature, correspondingly.

where h is a parameter for the smooth Dirac delta.

The phase change during the melting and solidification process is very complicated. It is a nonlinear problem due to the absorption and releases of latent energy at the melting point. The phase change is a temperature-dependent process. The enthalpy–porosity method [20] is used to simulate this phenomenon. It is assumed that the phase change occurs only when the temperature reaches the melting temperature $Tm$⁠. According to the temperature field, the computational domain can be decomposed into three regions: the solid region, the liquid region, and the mushy zone, as shown in Eq. (16). We use one and zero to indicate the liquid and solid phase, respectively. The region with a fraction between zero and one is the mushy area. This area is considered as a porous medium, where both the porosity and the fluid velocity decrease to zero in the solid area.

where $ϵ$ denotes the half width of the transition zone across the liquid and solid phases. The temperatures of the solid and liquid phases, as well as the enthalpy describing the total heat content of each phase, are represented by the following equations:

where $ρ,t,V,μ,k,T,p,K,h, and cp$ are density, time, velocity, molten fluid dynamic viscosity, heat conductivity, temperature, pressure, permeability, enthalpy, and heat capacity, respectively. The subscription s represents solid; l represents liquid; $Lm$ represents the latent heat for melting, which is the absorbed energy when phase changes from solid to liquid.

where $fs$ and $fl$⁠, $gs$ and $gl$ refer to mass fractions and volume fractions for solid and liquid phases, respectively; $cp$ is the heat capacity; $cps and cpl$ are the specific heat for solid metal and liquid metal.

The laser beam is modeled as a Gaussian beam, which means the intensity I takes a Gaussian distribution. The highest intensity exists in the center of the laser beam and keeps decreasing while away from the center. The intensity distribution of laser power along the distance of the laser beam center is shown as in Fig. 2 

where $η$ is the absorptivity coefficient; R and r are the beam radius and the distance from the calculated point to the beam center; and $Plaser$ stands for the power of laser beam. The laser beam moves starting from the center on the top surface. In the 2D model, we use a sine function along the x direction to simulate the motion of laser beam.

where the subscript g denotes the gas and the subscript m stands for the metal. The area with the negative sign of the level set function describes the metallic phases. The area with the positive sign designates the gas.

When applied to solve the multiphase flow problem, the level set function must be reinitialized every several iterations to maintain the numerical stability [22] and conserve mass of the phases [23]. In this research, we implement the model with the multiphysics finite element software comsol 5.0. The phase-injection method proposed by Zimmerman [24] is used to enforce the mass conservation. The basic idea is to add a constraint into the level set function to match with the calculated mass loss.

where $φ$ represents the mass change before and after the level set function is updated. The constraint is only applied to the interface, and the level set function need not be a signed distance function.

In Eq. (31), $ε$ represents the half width of the transition zone along the interface, where the smoothed Heaviside function changes from zero to one across the boundary [25].

where $ρ,cp$⁠, and k represent the density and heat capacity, heat conductivity, respectively. The subscripts m and g represent metal and gas.

where $hc$ is the heat transfer coefficient, $σ$ is the Stefan–Boltzmann's constant, $ε$ is the surface radiation emissivity, and $T∞$ is the atmosphere temperature.

Also, the thermal contact between melt pool and the walls is neglected.

In this experiment, the algorithm for free surface tracking and surface tension is tested. In the verification studies, the material properties of Ti–6Al–4V are considered constant in both solid and liquid phases. In this part, the model is built to simulate melting and spreading of a single particle onto a substrate. The particle is assumed to be spherical, which in this 2D case is a circle with radius 0.3 mm. The size of the study domain is 1.5 mm long and 0.7 mm wide. The particle is represented by a level set function $ϕ$⁠. The interface between the particle and gas locates where $ϕ$ equals zero. The objective of this model is to study the case involving both the free surface flow and the melting process of a metal particle. The effects of convection are neglected. The method can be extended to a more complex situation.

where $x,y,r$ refer to the x, y coordinates and the radius of metal particle.

The melting and spreading of a single particle problem is solved for a 10 ms period of time with a time step 0.01 ms. Experimental results are shown in Fig. 3. At each time step, the figures in the first row display the changes of the particle. Figures in the second row refer to the changes of temperature field as time increases.

From the observation of computational results, the process can be broken down into three parts:

1. (1)

   Heating (0–0.39 ms): during the heating procedure, the temperature in the domain keeps increasing. Because the heat conductivity is higher in metal, the heat transfers faster in the metal particle.

2. (2)

   Phase change (0.40–0.42 ms): during this period, the particle melts. The latent heat absorption is observed within the range of transition temperatures in Fig. 4. The contours represent temperature, and the legend illustrates the heat capacity caused by latent heat.

3. (3)

   Spreading (0.43–10 ms): The temperature of the whole computational domain is higher than the melting temperature. The phase change has completed in the particle. There are only two phases in the domain: liquid and gas. The liquid spreads under the effects of surface tension and gravity.

The level set function at zero level can clearly show the deformation of the surface as illustrated in Fig. 5.

In Eq. (42), the first term on the right side $qlaser$ stands for the heat power of the laser beam; the second and third terms in the right side represent the heat loss between the melt pool and ambient gas; $σ, ε, and hc$ are the Stefan–Boltzmann constant, radiation emissivity, and heat transfer coefficient, respectively; $T∞$ refers to the ambient temperature, which equals to 500 K in this study. By applying the delta function, we only apply the heat loss to the L/G interface.

Here, instead of using constant thermal properties, we consider temperature-dependent properties in both solid and liquid phases in Table 1, the value of the constant in this model is given in Table 2. The metal powder is distributed in a rectangular domain with 1.5 mm long and 1 mm wide. The laser beam is set at the point (0.25 mm, 0.6 mm) with a radius of 3 mm. Initially, the geometry of the powder bed is shown in Fig. 6. The laser beam intensity is 1000 W. The constants of physical properties and the process parameters are detailed in Tables 1 and 2.

Finally, the computational domain is discretized with a quadrilateral mesh with the maximum element size equal to $1.8×10−5$ mm. The capillary force, buoyancy force, and thermos-capillary force are applied to the L/G interface.

Figure 6 shows the simulation results. At each time spot, the figure on the left side shows the melt pool geometry, the temperature, and velocity field. The numbered contours refer to the temperatures with the SI unit of Kelvin (K); the black contours are the interface between metal and gas. The arrows refer to the velocity field of the molten metal flow. It is observed that the highest temperature happens in the area nearest to the center of the laser beam. By the end of the simulation, the highest temperature reaches to 2750 K. The particles change from the solid phase to the liquid phase. The deformation of molten particles can be observed. Figures on the right side demonstrate the liquid faction evolution during the phase change. From the results, the phase transition begins in the area closest to the laser beam. Only regions higher than the melting point are simulated as liquid. This is in accordance with the fraction function discussed in Sec. 2.2. The phase change is an energy absorbing or releasing process. In this case, the absorbed energy while melting is called the absorbing latent heat assimilated only within the range of transition temperatures.

As discussed in Sec. 1, the problem is a highly nonlinear multiphysics one. There are three models to be solved: the heat transfer, the fluid flow, and level set. The heat transfer and fluid flow are physics-related models, and the geometric level set method is only used to track the interface. To find an efficient way to solve the problem, both the fully coupled solver and the segregated solver are utilized in this study.

The fully coupled solver solves the heat transfer, the fluid flow partial differential equations, and the level set equation simultaneously. This approach is straightforward but is memory costing and time consuming. The segregated approach first solves the physics governing equations and then the level set equation at each step. The computer solves one specific model at a time and inherits other state variables from the previous step [28]. The comparison of the convergence plots of each solver is shown in Fig. 7.

From the observation, by using the segregated solver, the number of total computational time steps reduces from 180 to 155, which means the approach is faster than the fully coupled approach. Also, the convergence plot for segregated solver is smoother than the full direct couple solver. So, we can conclude that the segregated solver shows greater efficiency in solving this multiphysics problem.

In this study, a 2D numerical model has been built to simulate the melting process of metal-based SLM. The metal melting process is studied in different aspects: the temperature field, the velocity field, and the interface between liquid and gas phases. The model couples heat transfer and fluid dynamics with level set methods to trace the boundary of melt pools. A segregated finite element solver is adopted to solve the multiphysics models in a sequential way. Compared with the fully coupled solver, the segregated solver shows higher efficiency.

There are still some issues needing to be examined and addressed in the future research. In the current implementation, the effects of evaporation are neglected, but in the real AM process, the temperature in the melt pool may exceed the evaporating point. So, the mass loss during evaporation should be considered later. The wetting effects between molten particles with unmolten particles and the thermal contacts between melt pool and the substrate are not studied. The above issues will be further explored in our future research. Also, the reinitialization of the level set function must be further improved to avoid the possible interface shifting. The current work focuses on 2D particles. The extension of the proposed method to a more complex 3D model will be further investigated. For our next step, the simulation of the residual stress, which leads distortion and deformation in the parts [29], will be predicted using the stated variables calculated in our model [30]. The deformation can be presented by level set function as well. In addition, the experimental validation is needed in future work [1].