Modeling and Simulation of Metal AM

Abstract

Additive manufacturing with metal components is a complex, and currently cyclic, process due to the physical phenomena that are occurring. These phenomena can be mathematically modeled in order to predict the outcome of a specific aspect of the build. Coupling the mathematical models can then be used to develop a complete simulation, which can produce estimates for a range of characteristics for a part built using additive manufacturing techniques. This chapter will investigate the main models used in the simulation of metal AM. These models will include the modeling of thermal behavior, fluid dynamics, stress, and a selection of other auxiliary models, which are necessary to complete the simulations. For each of the models investigated, the various modeling techniques that have been developed will be presented along with their limitations, validation techniques, and parameters necessary to model the process correctly.

1. Introduction

1.1. Metal AM

Additive manufacturing (AM) has its roots in the development of stereolithography (SLA) in 1951 by Munz [1]. The development of SLA slowly escalated into a powder bed process and eventually into all of the various forms of AM that can be seen on the market today. In its infancy, the material and process controls only allowed AM to make prototypes and casting inserts. This was mainly due to the very low mechanical strength of the plastics being used and the process controls that produced very poor surface finishes [2]. This lacking drove researchers to find new and better materials and processes in order to use these methods for manufacturing as opposed to just prototyping. Over the last several decades, this has been accomplished and the American Society for Testing Materials (ASTM) has defined seven unique AM processes (binder jetting, directed energy deposition, material extrusion, material jetting, powder bed fusion, sheet lamination, and vat photopolymerization) [3]. The transition from rapid prototyping (RP) to AM resulted predominantly through the innovation of using metals in AM techniques. Of the seven standard AM process, three have been heavily investigated for their application to metal AM, powder bed fusion, directed energy deposition, and sheet lamination. In the process of sheet lamination, ultrasonic vibrations are used to fuse the sheets together. Due to its extreme dissimilarity to the other processes, namely the metal is never melted, it will be excluded from this body of work. The other processes, powder bed fusion and directed energy deposition, are so similar that many of the mathematical models used can be applied to both processes. To understand how these models are used it is best to understand the processes they simulate.

The first method to be developed for metal AM was powder bed fusion. In this process, shown in Figure 1a, the powder is delivered to the work area by a powder spreader of some type. The spreader is typically a roller, a pliable wiper, or a rigid wiper. This spreader will move the powder from the delivery system, either another bed or a hopper, to the work area. Once the powder is spread evenly across the surface of the work area, an energy source, typically a laser or electron beam, is used to selectively melt the powder into the desired solid shape. This process is repeated by lowering the work area and applying more powder with the spreader [4].

Figure 1.

Schematics for most common metal AM processes.

The process of directed energy deposition, shown in Figure 1b, is radically different from the powder bed process. In this process an energy sources, again typically laser or electron beam, is used to melt a small portion of a substrate, sometimes called the build plate. Once a melt pool is generated, the material is fed into the melt pool. This material can take the form of either wire feedstock or powder, which is blown into the melt pool. The laser and feedstock are then traced along the surface of the substrate until the first layer is built. The machine is incremented in the z direction and the process is repeated. This is repeated until the entire part is built. As can be inferred, this process is not just limited to a 3 axis machine and more degrees of freedom can allow for more agility when building the part [5].

1.2. Modeling basics

In modeling a physical process, it is necessary to divide the space being modeled so that the equation can be computed. For AM processes, space representation has been divided into two main domain: Eulerian and Lagrangian. The Eulerian domain, which can also be called a field coordinate system, tracks specific locations in space to determine what is happening at a given location. This space, which can be seen in Figure 2a, is typically the space for finite element analysis (FEA). FEA can be used to simulate most physics processes. This process breaks the material down using a mesh structure and then the nodes are used to calculate the physics, which are being modeled. This mesh facilitates a quick and easy identification of neighbors using the connections of nodes found with the mesh. This technique is very beneficial when attempting to solve differential equations and has led to several solution methods being developed. However, since properties are typically saved as floating point values in the model, round off errors can lead to an arbitrary reduction or inflation of there properties.

Figure 2.

Depiction of domains used in metal AM modeling.

On the contrary, the Lagrangian domain, which can also be called a material coordinate system, tracks matter as it moves and reacts within the simulation space. In this domain, the material that is used will be broken up into discrete representative particles. Each of these particles is then tracked through space and time, calculating the various physical processes that apply to the particles. The Lagrangian method can be very beneficial for a couple of reasons. For starters, matter is never created nor destroyed by rounding of floating point numbers, due to the fact that the particle is given a fixed mass at the beginning of the simulation. In addition, the forces that the fluid is subject to are inherent within the simulation setup. Therefore, the pressure the fluid undergoes can be recorded without any additional calculations. And lastly, since particles are tracked there is no limit to where they are allowed to move. This can result in a much more accurate method of simulating fluid flow. This technique is not without its struggles. Since the material is divided into many little particles it can be very memory intensive and require substantial computing power. Additionally, this technique can make solving the differential equation that defines a system challenging due to the lack of inherent connectivity [6].

2. Thermal modeling

2.1. Modeling techniques

The heat transfer that occurs within a part can be broken down into three main types: conduction, convection, and radiation. Each of these can be calculated individually or the three can be combined to garner a full representation of the physical processes, which are taking place. For the sake of simplicity, each will be discussed on its own to highlight the necessary inputs and responses and then a complete picture will be drawn at the end.

The first physical phenomenon to be discussed is conduction. The heat conduction equation (or the heat diffusion equation) was first proposed by Jean Baptiste Joseph Fourier in 1807, resulting in it being one of the most mature modeling methods in AM. The full representation, in modern notation, can be seen in Eq. (1),

$$\nabla \cdot k\nabla T=c\frac{\mathit{\delta T}}{\mathit{\delta t}}$$ E1

where $k$ is the thermal conductivity, $T$ is the temperature, $c$ is the heat capacity of the material and $t$ is time [7]. To understand this equation, it can be helpful to look at a simplified model, namely the 1-dimensional steady state example, Eq. (2),

$${q^{\prime }}^{\prime }=-k\frac{\mathit{dT}}{\mathit{dx}}=k\frac{T_1-T_2}{L}$$ E2

where ${q^{\prime }}^{\prime }$ is the heat flux transferred within the part. From this representation of the equation, it can be seen that conduction is simply the material attempting to reach its lowest energy state, namely constant temperature. Though helpful for understanding, this approach applies to virtually no actual applications in AM modeling, for this, it is necessary to create a 3-dimensional transient model, shown in Eq. (3) [8].

$$\frac{{\delta }^{2}T}{\delta x^2}+\frac{{\delta }^{2}T}{\delta y^2}+\frac{{\delta }^{2}T}{\delta z^2}=\frac{\mathit{\rho c}}{k}\frac{\mathit{\delta T}}{\mathit{\delta t}}$$ E3

With the heat conduction within the part being handled by conduction, the other two phenomena only affect the surfaces of the material and are referred to as boundary conditions. The first, convection, is like conduction in the fact that utilizes contact to transfer heat from one object to another. Unlike conduction, which is between solid materials, convection is between a solid and a liquid, where air can be treated as a fluid. Newton, in 1701, saw that the relationship could be written as Eq. (4),

$$\frac{\mathit{dT}}{\mathit{dt}}\propto -\left(T-T_a\right)$$ E4

where $T$ is the temperature of the material, $T_a$ is the temperature of the fluid (i.e. ambient temperature), and $t$ is time. The proportionally was later found for various situations and became known as the heat transfer coefficient (h) [9].

The last phenomena is a contact-less heat transfer through radiation. The law of radiation was first experimentally suggested by Josef Stefan in 1879 and later analytically shown by Ludwig Boltzmann in 1884, together these men derived the Stefan-Boltzmann Law, shown in Eq. (5),

$${q^{\prime }}^{\prime }=\mathit{\epsilon \sigma }{T}_s^4$$ E5

where ${q^{\prime }}^{\prime }$ is the heat flux due to radiation, $\epsilon$ is the emissivity of the surface, and $\sigma$ is the Stefan-Boltzmann constant [10]. In order to implement these into a model, the process of conduction must be applied to the entirety of the material, whereas, the convection and radiation must only be applied to the skin as boundary conditions to ensure that the part cools properly based on the environment where it is located.

As can be seen from these equations, there are only five material properties, outline in Table 1, that are needed in order properly model the thermal history of an AM process. Most of these values must be found experimentally and therefore can be found in several locations including literature, material handbooks, and supplied by material vendors.

Thermal conductivity ( $W/\mathit{mK}$ )

Specific heat ( $J/\mathit{kgK}$ )

Density ( $\mathit{kg}/m^3$ )

Convective coefficient ( $W/m^{2}K$ )

Emissivity

Table 1.

Material properties required for thermal modeling.

2.2. Validation techniques

Once the mathematical models have been developed it is necessary to show that they are a true representation of reality. For thermal modeling in AM, this has been done through two main methods, direct measurement of temperature with an instrument or by measuring a more easily attained physical characteristic and relating that to the temperature. As can be seen from a representative set of papers in Table 2, more effort has been placed on using an instrument to validate the model as opposed to measuring another characteristic. This is due to the fact that the physical characteristic comparison relies on the correctness of all models that are involved and does not have a one to one correlation between the material temperature and the model. Each of these techniques has a unique set of applications, strengths and concerns.

Instrument validated		Physical Char. validated
IR/CCD camera	[11, 12, 13, 14]	Melt pool depth	[15, 16, 17]
Pyrometer	[18, 19]
Thermal couple	[19, 20, 21]

Table 2.

Breakdown of validation techniques.

The use of a camera, be it IR or CCD, is a very appealing technique for many researchers for several reasons. The first key feature is that a camera is a non-contact measurement. This has several implications, namely, the application to all methods of metal AM developed to date and ability to be used in-situ. Since temperature is time and spatially dependent in nature, in order to accurately measure the temperature of a material both are critical, which in a camera is related to its frame rate and resolution. Researchers have reported frame rates of 800 frames/sec [11] and resolutions of 0.1 × 0.1 mm [12]. Though a powerful tool, cameras are not without their faults. One fault, related to the frame rate, is that a camera measures the average temperature of the skin during the time the shutter is open. This means that any extremely abrupt changes in the temperature can be smoothed by the camera. This problem only materializes when a pulsed laser is used, which has a different interaction with the material that includes an extremely fast and short spike in temperature [22]. The most critical consideration, however, is that cameras are highly sensitive to distance and angle to the substrate [23]. It is critical to ensure that both of these are controlled to ensure that an accurate measurement is made.

The next instrument used, a pyrometer, is less accurate than the camera but is simple which has allowed for it to be used as a feedback sensor where a mathematical model allows for the estimation of the expected sensor reading. A pyrometer is a non-contact spot instrument that is capable of measuring the average temperature of an area of the material. It measures the thermal radiation that is emitted from the area of the material where it is aimed. From Planck’s radiation law and knowing the spectral distribution of the blackbody emissive power it is possible to derive Eq. (6) that is the effective temperature a pyrometer reads, assuming a Gaussian laser and a 1 mm radius for the pyrometer,

$$T_{\mathit{eff}}=\frac{\mathit{hc}}{\mathit{\lambda \sigma }\mathit{ln}\left(1+n/{\sum }_{i=1}^n{\left(e^{\frac{\mathit{hc}}{\mathit{\lambda \sigma }{T}_i}}-1\right)}^{-1}\right)}$$ E6

where $h$ is Planck’s constant, $c$ is the speed of light, $\lambda$ is the wavelength of the emitted radiation, $\sigma$ is Stefan-Boltzmann constant, $n$ is the number of small sampling areas within the pyrometer viewing area, and $T_i$ is surface temperature within the small $n$ areas [24]. This equation can then be included in a mathematical model to simulate the effect of a feedback system, which will give a much more accurate simulation.

The last instrument used for validation is a thermocouple, which is a contact spot measurement of the temperature. Thermocouples are limited to DED application because they must be fixed to the material. It also can not be used close to the melt pool or it will become detached and give inaccurate results. Lastly, it only is able to measure the average temperature of a specific location of the material, similar to the pyrometer. However, unlike the pyrometer that can be moved to follow the melt pool, thermocouples are fixed. This results in needing an array of thermocouples in order to garner a full understanding of the heat distribution within the material. For these reasons, thermocouples are typically used as a secondary validation.

The final validation technique is to use a more readily measured physical characteristic, namely melted track size, to show the correctness of the modeling approach. This method takes either a laser heating of the substrate or a single layer build and compares the results with the simulations findings. This is possible since the melting temperature of the material is known, therefore, the region of the simulation that reached a higher temperature can be flagged for analysis when the simulation is complete. The experiment is then sliced and the width and depth of the melted region are measured and then compared to the simulation. This method of validation is preferred by some due to the lack of needing any special equipment. However, the methods rely on the correctness of the fluid motion model, discussed in Section 3, and the thermal model previously discussed. This can lead to false rejection or validation of the model due to outside variables.

3. Fluid motion modeling

Fluid modeling can be done in either the Eulerian domain or the Lagrangian domain. If the Eulerian domain is used the Volume of Fluid (VOF) technique is implemented. However, if a Lagrangian domain is used then either a Smoothed Particle Hydrodynamics (SPH) or Position Based Dynamics (PBD) can be implemented.

3.1. Volume of fluid

VOF is a method that is based on a Eulerian domain and tracks the amount of fluid that is within a specific region as a percentage of a fully dense region. In this method, the region being investigated is broken down into a mesh, like all Eulerian domain problems. Each mesh grid is then given an F value based on the amount of fluid that is within the cell. This F value denotes the percentage of the cell that is occupied by fluid; therefore, a value of 1 represents a cell that is full of fluid, a value of 0.5 represents a cell that is half full of fluid, and a value of 0 represents a cell that contains no fluid.

The definition of the F term allows for the fluid boundary to be readily determined by finding all cells with an F value between, but not equaling, 0 and 1. Calculating the gradient of F within the domain also allows for the normal of the fluid surface to be defined as the direction of most rapid change. It is possible to step the F value through time with Eq. (7).

$$\frac{\mathit{\delta F}}{\mathit{\delta t}}+u\frac{\mathit{\delta F}}{\mathit{\delta x}}+v\frac{\mathit{\delta F}}{\mathit{\delta y}}=0$$ E7

This is coupled with the Navier–Stokes equations, Eq. (8) and Eq. (9), and the continuity equation, Eq. (10), to find the motion that the fluid endures [25],

$$\frac{\delta \left(\mathit{\rho u}\right)}{\mathit{\delta t}}+\frac{\delta \left(\rho u^2\right)}{\mathit{\delta x}}+\frac{\delta \left(\mathit{\rho uv}\right)}{\mathit{\delta y}}=-\frac{\mathit{\delta p}}{\mathit{\delta x}}+g_x+\nu \left[\frac{{\delta }^{2}u}{\delta x^2}+\frac{{\delta }^{2}u}{\delta y^2}+\xi \left(\frac{1}{x}\frac{\mathit{\delta u}}{\mathit{\delta x}}-\frac{u}{x}\right)\right]$$ E8

$$\frac{\delta \left(\mathit{\rho v}\right)}{\mathit{\delta t}}+\frac{\delta \left(\mathit{\rho vu}\right)}{\mathit{\delta x}}+\frac{\delta \left(\rho v^2\right)}{\mathit{\delta y}}=-\frac{\mathit{\delta p}}{\mathit{\delta y}}+g_y+\nu \left[\frac{{\delta }^{2}v}{\delta x^2}+\frac{{\delta }^{2}v}{\delta y^2}+\frac{\xi }{x}\frac{\mathit{\delta v}}{\mathit{\delta x}}\right]$$ E9

$$\frac{1}{c^2}\frac{\mathit{\delta p}}{\mathit{\delta t}}+\frac{\mathit{\delta u}}{\mathit{\delta x}}+\frac{\mathit{\delta v}}{\mathit{\delta y}}+\frac{\mathit{\xi u}}{x}=0$$ E10

where $\left(u,v\right)$ is the velocity in the $\left(x,y\right)$ directions for Cartesian coordinate system or $\left(r,z\right)$ directions for cylindrical coordinate systems, $\rho$ is the density, $\xi$ is a coordinate specific variable 0 for Cartesian and 1 for cylindrical, $\left(g_x,g_y\right)$ is the body acceleration, $\nu$ is the kinematic viscosity, $p$ is the pressure, $t$ is time, and $c$ is the adiabatic speed of sound in the fluid.

In order to step through time, the following steps must be followed:

1. Compute the first guess of the new velocities with explicit approximations of Eq. (8) and Eq. (9) using previous time step values as the boundary conditions.

2. Adjust the pressure in each cell to satisfy Eq. (10) and change velocities of adjacent cells accordingly. Iterate until pressures and velocities are found such that all equations are balanced.

3. Update F with Eq. (7) from the solved values

In order to model with the VOF techniques, the material properties outlined in Table 3 need to be known. Similar to heat transfer material properties, these properties must be found experimentally and can be found throughout literature and material properties handbooks.

Density ( $\mathit{kg}/m^3$ )

Kinematic viscosity ( $\mathit{Pa}.s$ )

Table 3.

Material properties required for VOF fluid modeling.

Though a mature and robust method, one of the main disadvantages of this method is its reliance on a Eulerian domain. This results in the fluid mass in a cell containing a free surface being represented by a floating point number. In a computer system, the precision is limited, which will lead to round off errors. In a technique as VOF, this can lead to loss or creation of mass based on the rounding of these numbers. This, in turn, breaks the law of conversation of mass. Another problem that arises from the use of VOF techniques is that the fluid is only allowed to move within predetermined space becuase the technique relies on the eularian domain. Even if some of the fluid physically would move outside this domain, this model will not allow this resulting in the potential for inaccurate results.

3.2. Smoothed particle hydrodynamics

The SPH method utilizes a Lagrangian domain and approximates the fluid as a set of particles that are tracked through space. This method is based on the fact that a function $f\left(x\right)$ can be written with an integral representation, Eq. (11),

$$f\left(x\right)={\int }_{\mathrm{\Omega }}f\left(x^{\prime }\right)\delta \left(x-x^{\prime }\right)d{x}^{\prime }$$ E11

where $\mathrm{\Omega }$ is the domain and $\delta$ is a Dirac function. For the SPH technique, it is then necessary to convert from a continuous function to a discrete function because the fluid is represented by a set of particles and not a continuous fluid. Additionally, the Dirac function can be replaced with a weighting function to result in Eq. (12),

$$f\left(x\right)=\sum _{b=1}^{N}f\left(x_b\right)W\left(x-x_b,h\right)\mathrm{\Delta }{V}_b$$ E12

where $b$ are the neighboring particles, $W$ is the weighting function, and $h$ is the radius of influence of the weighting kernel.

The weighting function acts as a smoothing function that gives a weight to the influence that the particles have on each other. This weighting function results in particles that are closer together having more influence than those that are farther away. The various weighting kernels have been proposed in literature including Gaussian Kernel Eq. (13), Bell shaped kernel Eq. (14), or Cubic spline kernel Eq. (15),

$$W\left(r,h\right)={\alpha }_d{e}^{-{\left(r/h\right)}^2}$$ E13

$$W\left(r,h\right)={\alpha }_d\left\{\begin{array}{ll}\left(1+3\frac{r}{h}\right){\left(1-\frac{r}{h}\right)}^2& \frac{r}{h}\le 1\\ 0& \frac{r}{h}>1\end{array}\right.$$ E14

$$W\left(r,h\right)={\alpha }_d\left\{\begin{array}{ll}\frac{3}{2}-{\left(\frac{r}{h}\right)}^2+\frac{1}{2}{\left(\frac{r}{h}\right)}^3& 0\le \frac{r}{h}<1\\ \frac{1}{6}{\left(1-\frac{r}{h}\right)}^3& 1\le \frac{r}{h}<2\\ 0& \frac{r}{h}\ge 2\end{array}\right.$$ E15

where $r$ is the distance between the particles begin investigated and ${\alpha }_d$ is a normalizing scalar.

The conservation of mass, Eq. (16), within this system is the first equation that needs to be written in this approximates form, Eq. (17)

$$\frac{\mathit{D\rho }}{\mathit{Dt}}=u\nabla \rho -\nabla \left(\mathit{\rho u}\right)$$ E16

$$\frac{D{\rho }_a}{\mathit{Dt}}=\sum _{b=1}^N{m}_b{u}_{\mathit{ab}}{\nabla }_a{W}_{\mathit{ab}}$$ E17

Next, conservation of momentum, Eq. (18), must be written in it’s approximate form, Eq. (19).

$$\frac{\mathit{Du}}{\mathit{Dt}}=-\frac{1}{\rho }\nabla p$$ E18

$$\frac{{\mathit{Du}}_a}{\mathit{Dt}}=-\sum _{b=1}^N{m}_b\left(\frac{p_b}{{\rho }_b^2}+\frac{p_a}{{\rho }_a^2}+{\mathrm{\Pi }}_{\mathit{ab}}\right){\nabla }_a{W}_{\mathit{ab}}+F$$ E19

where the artificial viscosity ${\mathrm{\Pi }}_{\mathit{ab}}$ is defined by Eq. (20),

$${\mathrm{\Pi }}_{\mathit{ab}}=\left\{\begin{array}{ll}\frac{-\alpha c_{\mathit{ab}}{\mu }_{\mathit{ab}}+{\mathit{\beta \mu }}_{\mathit{ab}}^2}{{\rho }_{\mathit{ab}}}& u_{\mathit{ab}}\left(x_a-x_b\right)<0\\ 0& u_{\mathit{ab}}\left(x_a-x_b\right)\ge 0\end{array}\right.$$ E20

where $c_{\mathit{ab}}$ is the average speed of sound of the particles, ${\rho }_{\mathit{ab}}$ is the average density of the particles, and $\mu$ is defined in Eq. (21).

$${\mu }_{\mathit{ab}}=\frac{{\mathit{hu}}_{\mathit{ab}}\left(x_a-x_b\right)}{{\left(x_a-x_b\right)}^2+{\left(0.01h^2\right)}^2}$$ E21

The volume forces, $F$ , is any force that is applied to the entire fluid, such as gravity. One final equation that is needed is the equation of state, which links the pressure and density, for an ideal gas it can be defined as Eq. (22) or for a quasi-incompressible fluid it can be defined as Eq. (23),

$$p=\frac{\mathit{RT}}{M}\left(\frac{\rho }{{\rho }_0}-1\right)$$ E22

$$p=c_0^2{\rho }_0\left(\frac{\rho }{{\rho }_0}-1\right)$$ E23

where $R$ is the ideal gas constant, $T$ is the temperature, $M$ molar mass, $c_0$ is the speed of sound in the material [26].

In order to step through time, the following steps must be followed for each particle:

1. Find all neighbors within smoothing kernel length

2. Solve Momentum Equation by solving for Artificial Viscosity (Eq. (20)) and Equation of State (Eq. (22) or Eq. (23)) and applying Volume Forces

3. Solve Continuity Equation (Eq. (17))

4. Update $\rho$ , x, u

Even though it has some very appealing properties, SPH is not without its faults. The main fault applying to AM modeling is that the SPH technique is very sensitive to density fluctuations resulting from a lack of neighbors. In AM simulation, this problem can arise when a simulation is beginning to develop a melt pool, when a melt pool is in the end stages of solidification, or if the energy source is moving so fast that only a small amount of material is melted. This fluctuation in density can be combated with extremely small time steps, however, this makes the computation time extremely long [27].

3.3. Position based dynamics

PBD is also a Lagrangian particle based method, similar to SPH. This method approximates the domain with particles that are subject to the standard Newtonian forces. These forces are applied to all of the particles and their position in space is tracked to determine the motion of the fluids.

In order to solve for the position of the particle sets, several functions must be applied to the particles. The first function that needs to be defined is the scheme for updating the velocities of the set of particles. This is done by adjusting the velocities based on any external forces, Eq. (24),

$$v_i^{\ast }=v_i+\mathrm{\Delta }{\mathit{tf}}_{\mathit{ext}}$$ E24

where $v_i^{\ast }$ is the new velocity of the particle, $v_i$ is the current velocity of the particle, $\mathrm{\Delta }t$ is the time step, and $f_{\mathit{ext}}$ is the sum of the external forces that are applied to the particle. From this, a guess is made as to the new positions of the particles, Eq. (25).

$$x_i^{\ast }=x_i+\mathrm{\Delta }{\mathit{tv}}_i$$ E25

Next, it is necessary to define the equations that enforce incompressibility of the fluid. The first equation is a function of the position ( $p_1\dots p_n$ ) of all of the neighbors, Eq. (26),

$$C_i\left(p_1\dots p_n\right)=\frac{{\rho }_i}{{\rho }_0}-1$$ E26

where ${\rho }_0$ is the rest material density and ${\rho }_i$ is the standard SPH density Eq. (19). Since the density of the fluid is assumed to be constant, it is possible to find the pressure change each particle experiences, Eq. (27),

$$\mathrm{\Delta }p\approx \nabla C\left(p\right)\lambda$$ E27

where $\lambda$ is a scalar. From this, it is possible to solve for $\lambda$ , Eq. (28)

$${\lambda }_i=-\frac{C_i\left(p_1\dots P_n\right)}{\sum _k{\left|{\nabla }_{p_k}{C}_i\right|}^2}$$ E28

Once $\lambda$ is known, it is then possible to determine the change in pressure that each particle experiences with Eq. (29)

$$\mathrm{\Delta }{p}_i=\frac{1}{{\rho }_0}\sum _j\left({\lambda }_i+{\lambda }_j\right)\nabla W\left(p_i-p_j,h\right)$$ E29

Once the pressure change is known, it is then possible to generate the actual positions ( $x_i^{\ast }$ ) of the particles with Eq. (30).

$$x_i^{\ast }=x_i+\mathrm{\Delta }{p}_i$$ E30

And lastly, with the actual motion of the particles known, it is possible to determine the final velocity of the particles from Eq. (31).

$$v_i^{\ast }=\frac{1}{\mathrm{\Delta }t}\left(x_i^{\ast }-x_i\right)$$ E31

To step through time each of these steps is to be applied to each particle in the simulation before moving on to the next step [27]:

1. Apply external forces (Eq. (24)) and predict the position (Eq. (25))

2. Find the nearest neighbors

3. Calculate the ${\lambda }_i$ (Eq. (28))

4. Calculate change in pressure (Eq. (29)) and collision detection

5. Update positions (Eq. (30)) and velocities (Eq. (31))

Like SPH, this method is appealing due to its inherent conservation of mass and Lagrangian representation does not limit the flow to a specific area. However, this technique relies on arbitrary values to generate an accurate appearing simulation. For this reason, tuning is needed to develop a set of inputs, as opposed to running an experiment to gather the necessary material properties.

4. Stress modeling

Once the evolution of heat through the process has been found, it is then possible to take this data and estimate the stresses and deformations that are the result of the cyclic process.

4.1. Modeling techniques

The total stress that accumulates during an AM process can be written as the sum of the individual contributors of the strain, Eq. (32),

$${\epsilon }_{\mathit{ij}}={\epsilon }_{\mathit{ij}}^E+{\epsilon }_{\mathit{ij}}^P+{\epsilon }_{\mathit{ij}}^T+{\epsilon }_{\mathit{ij}}^{\mathrm{\Delta }V}+{\epsilon }_{\mathit{ij}}^{\mathit{TrP}}$$ E32

where ${\epsilon }_{\mathit{ij}}$ is the total strain, ${\epsilon }_{\mathit{ij}}^E$ is elastic strain (Eq. (33)), ${\epsilon }_{\mathit{ij}}^P$ is the plastic strain (Eq. (34)), ${\epsilon }_{\mathit{ij}}^T$ is thermal strain (Eq. (35)), ${\epsilon }_{\mathit{ij}}^{\mathrm{\Delta }V}$ is strain from volumetric dilation strain (Eq. (36)), and ${\epsilon }_{\mathit{ij}}^{\mathit{TrP}}$ is the strain from phase transformation (Eq. (37)).

$${\epsilon }_{\mathit{ij}}^E=\frac{1+\nu }{E}{\sigma }_{\mathit{ij}}-\frac{\nu }{E}{\sigma }_{\mathit{kk}}{\delta }_{\mathit{ij}}$$ E33

$${\epsilon }_{\mathit{ij}}^P=\widehat{G}\left(\frac{\mathit{\delta F}}{{\mathit{\delta \sigma }}_{\mathit{kl}}}+\frac{\mathit{\delta F}}{\mathit{\delta T}}T+\sum \frac{\mathit{\delta F}}{\delta X_i}{X}_i\right)\frac{\mathit{\delta F}}{{\mathit{\delta \sigma }}_{\mathit{ij}}}$$ E34

$${\epsilon }_{\mathit{ij}}^T=\alpha \left(T-T_0\right){\delta }_{\mathit{ij}}$$ E35

$${\epsilon }_{\mathit{ij}}^{\mathrm{\Delta }V}=\frac{\mathrm{\Delta }V}{3V}\mathrm{\Delta }X$$ E36

$${\epsilon }_{\mathit{ij}}^{\mathit{TrP}}=\frac{5S_{\mathit{ij}}\mathrm{\Delta }V}{4\mathit{YV}}\left(2-2X_n-\mathrm{\Delta }X\right)\mathrm{\Delta }X$$ E37

where $E$ is Young’s modulus, $\nu$ is Poisson’s ratio, $\alpha$ is the thermal expansion coefficient, $Y$ is the yield stress of the weaker phase, F is the yield function (Eq. (38)) and G is the hardening function (Eq. (39)) [28].

$$F=F\left({\sigma }_{\mathit{ij}},{\epsilon }^P,\kappa ,T,X_i\right)$$ E38

$$\frac{1}{{\displaystyle \widehat{G}}}=-\left(\frac{\mathit{\delta F}}{{\mathit{\delta \sigma }}_{\mathit{mn}}^P}+\frac{\mathit{\delta F}}{\mathit{\delta \kappa }}{\sigma }_{\mathit{mn}}\right)\frac{\mathit{\delta F}}{{\mathit{\delta \sigma }}_{\mathit{mn}}}$$ E39

In order to solve for the stress that is accompanied by the AM process, each of the individual components of strain need to be solved for and added together. This can be a tedious and computationally expensive process, based on the parameters in Table 4.

Young’s modulus

Poisson’s ratio

Thermal expansion coefficient

Yield stress of the weaker phase

Table 4.

Material properties required for stress modeling.

To reduce the simulation time, the inherent strain approach has been developed. In this method the strains are lumped together, as seen in Eq. (40), to simplify the calculations.

$${\epsilon }^{\ast }={\epsilon }_{\mathit{ij}}-{\epsilon }_{\mathit{ij}}^E={\epsilon }_{\mathit{ij}}^P+{\epsilon }_{\mathit{ij}}^T+{\epsilon }_{\mathit{ij}}^{\mathrm{\Delta }V}+{\epsilon }_{\mathit{ij}}^{\mathit{TrP}}$$ E40

This technique will run a small scale simulation to determine the value of this inherent strain. These values are then superimposed as needed when the full-scale simulation is being modeled. This results in a simulation that is orders of magnitude faster than the previously described technique [29].

4.2. Validation techniques

Once the models have been developed it is necessary to validate the technique. There are several methods that have been used, which are reported in Table 5.

Presence of cracks	[13, 30]
Final part displacement	[31, 32, 21]
Digital image correlation (DIC)	[33, 34]
Neutron diffraction	[35, 20]
X-ray diffraction	[36, 37]

Table 5.

Frequency of stress analysis techniques.

The simplest of the methods that have been used is to observe the locations of any cracks that have developed during the build. The experimental crack locations are then compared to the simulation to determine if a high stress is predicted. This method is appealing due to the lack of specialized equipment required. However, the lack of detail results in this method only being able to be used as a method of qualitative validation and not a precise comparison.

A more precise comparison is needed in order to quantitatively compare the simulation and the experiment, and this has led to the following techniques. The first technique is to measure the distortion that occurs within the part during the build process. This has been done with either a laser displacement sensor or a 3-D scanner. If a laser displacement sensor is to be used then the substrate must be setup as a cantilever, as in Figure 3, which will allow one end to move as the material is heated and cools. The fluctuations that are observed in the end of the substrate are then compared to the simulation to determine the accuracy of the simulation. Another method of measuring the distortion of the part is to scan the part with a 3-D scanner after the part has been removed from the fixturing and allowed to flex. Each of these methods has its advantages and problems. The first is the accuracy of the validation technique. For a higher accuracy, the laser displacement sensor should be used, and it has been reported to have an accuracy of $\pm 1 \mathit{\mu m}$ [21], whereas 3D scanners have reported to have an accuracy of $\pm 500 \mathit{\mu m}$ [32]. Coupled with the accuracy is the area that is measured to determine the distortion. The laser displacement sensor only measures a single spot and tracks that through space and time, whereas the 3D scanner is capable of capturing the distortion that occurred in a larger section of the build. An additional advantage that the 3-D scanner has over the laser displacement sensor is the effect of misalignment of the initial deposit is less when compared to the laser displacement sensor. However, the 3D scanner cannot be used in-situ and the laser displacement sensor can. These methods allow for a quantitative analysis of the distortion which can be used to determine the stress within the part.

Figure 3.

Experimental setup using laser displacement sensor to measure distortion.

Another method of tracking the part’s distortion is to use digital image correlation (DIC). This method uses a camera to observe the part and track any distortion by comparing frames. Points are selected on the material and tracked through space as shown in Figure 4.

Figure 4.

Schematic of mini-tensile specimen with tracking points for DIC.

This method will inherently return the distortion that a part undergoes and it then it is possible to precisely determine any surface level stresses. Measuring the stresses is possible by starting with a stressed specimen, for example, one that has undergone the deposition process. The part is then selectively stress relieved to determine the relaxation that the material undergoes. The relaxation that the material experiences will be directly linked to the stress [33]. This method is highly appealing because the motion of the material is measured in pixels of the camera. Because of this, the resolution of this method is limited by only the camera, which can be increased by using a microscope in conjunction with a camera. The major drawback of this technique is that it is a destructive method, resulting in it not being an acceptable method for some applications.

If the exact strain, the displacement of the atoms from their rest positions, is needed then it is necessary to measure the actual motion of the atoms. This is done using Bragg’s Law and the scattering of either X-Rays or neutrons. The initial material is placed in an apparatus that allows for the initial diffraction pattern to be generated. This is the initial location of the atoms before the thermal processes have been applied to the material. The material is then subjected to the thermal processes. The material is then placed back into the apparatus and the diffraction pattern is again measured. This new pattern is the final locations of the atoms. From these two sets, it is possible to find how far the atoms moved from their base locations, this value is the strain. This strain value can then be related to the simulation to determine the accuracy of the simulation [38]. Though an accurate method, the use diffraction patterns are not without its faults. The first is the access to researchers, it is necessary to have a dedicated setup for measuring the diffraction, which can make it difficult for some to gain access. Additionally, this method can only be used at the end of a build and cannot be used to measure in-situ strain.

5. Miscellaneous models

In order to completely model the AM process it can be necessary to include an energy source, laser or e-beam, and if powders are used it can be necessary to track their position within a powder bed or blown powder.

5.1. Laser modeling

The most common laser distributions used in metal AM process are top hat, Figure 5b, and Gaussian, Figure 5a. Just as the material must be divided into the domain, a laser must be divided into smaller pieces for proper modeling. This is done by dividing the projected area into segments, which can be thought of as rays from the laser. The flux that is generated within each of these divisions of the laser is found by taking the integral of the functions over one division. This was done for the Gaussian laser distribution where Figure 5a was generated by integrating Eq. (41),

$$\varphi \left(x,y\right)={\varphi }_0\sqrt{1-\frac{x^2}{r_0^2}-\frac{y^2}{r_0^2}}$$ E41

where $\varphi \left(x,y\right)$ is the intensity as a specific $\left(x,y\right)$ location, ${\varphi }_0$ is the initial laser power, and $r_0$ is the laser radius [19]. The top hat beam, Figure 5b, is much easier to model. In this model, the laser power is simply divided by the number of divisions that can be placed within the laser boundaries. The heat flux that is found is then applied to a specific location of the domain based on the current laser location. These intensities are then multiplied by $\alpha$ , the material absorptivity, to determine the actual amount of energy that is absorbed by the material. This process can be applied to any distribution that can be imagined.

Figure 5.

Example of heat flux from example distributions.

There are two methods of applying the laser. The first is to define the surface and apply the flux directly to this region of the material and not attempting to determine if there is anything blocking the laser projection. The other, and more realistic method, is to perform ray tracing and apply the heat flux to the first location that is struck by the ray. Both method work but ray tracing is preferred for the sake of reality but can be costly to compute.

5.2. E-beam modeling

The other choice of a heat source for metal AM is an electron beam. Typically an electron beam is modeled as a Gaussian heat source, just as the previously mentioned laser. However, an electron beam will also penetrate the surface heating the material for a given depth. The intensity in the z-direction can be expressed as Eq. (42).

$${\varphi }_z=\frac{1}{0.75}\left(-2.25{\left(\frac{z}{S}\right)}^2+1.5\frac{z}{S}+0.75\right)$$ E42

The xy and z intensities can then be combined to determine the heat flux, as shown in Eq. (43),

$${\displaystyle \dot{Q}}\left(x,y,z\right)=-{\eta }_b{\eta }_e{P}_B\frac{{\varphi }_{\mathit{xy}}{\varphi }_z}{S}=-{\eta }_b{\eta }_e{P}_B\frac{4\mathit{ln}\left(0.1\right)}{\pi d_b^{2}S}{e}^{\frac{4\mathit{ln}\left(0.1\right)\left(x^2+y^2\right)}{d_b^2}}\left(-3{\left(\frac{z}{S}\right)}^2-2\frac{z}{S}+1\right)$$ E43

where ${\eta }_b$ is the beam control efficiency, ${\eta }_e$ is the energy conservation at the part surface, $P_B$ is the beam power, $S$ is the absolute penetration depth, and $d_b$ beam diameter [39].

5.3. Powder bed and blown powder models

The last element of the metal AM process that is typically necessary is the addition of material. In the wire feed DED process this is done by modeling the wire as a solid material, just as the substrate, and treating it in a similar fashion. When using either a powder bed or blown powder, this is not possible. Due to the stochastic nature of the powders, it can be challenging to model their behavior. When modeling the powder bed setup, there are two prevailing methods that have been used: discrete element method (DEM) and geometric methods.

The DEM technique tracks the powder particles on an individual basis to determine their final position in the build volume. This simulation technique typically begins by dropping particles (the blue particles in Figure 6a), or sets of particles, from a random x and y position but a designated height in the domain. From there, the particles position and velocity (the red vectors in Figure 6a) are tracked as the particle is subjected to the major forces of gravity, contact, and friction. In some cases, smaller interaction forces, such as Van der Waal forces [40] or JKR interaction forces [41], are added to increase the accuracy of the simulation. The particles are allowed settle to their resting point (the green particles in Figure 6a) and more particles are added as needed to the simulation. This technique is very appealing to the simulation powder spreading in powder bed AM due to the ease of adding a powder spreader to the simulation without much additional effort. This results in the ability to simulation the entire process from the layering of the powder to laser interaction to melting and solidification [42].

Figure 6.

Powder bed modeling techniques.

The geometric method is not as realistic for powder bed process but is computationally only a small fraction of the DEM technique. In this technique, the area to be filled is analyzed without regard to how the powder would flow. The first geometric method is referred to as the compression algorithm. In this technique, the particles are randomly spaced within the domain. The particles are then moved in the direction of compression, typically the direction of gravity, by the shortest distance that results in a collision of particles. This compressing is repeated until the potential energy is below a user-defined tolerance. The particles are then “shaken”, or moved laterally, and the process is repeated. The build volume is then refilled with more particles and the process is repeated until the volume is full [43]. Another geometric method, displayed in Figure 6b, focuses on a tetrahedral mesh. In this method, the build volume is meshed and particles are placed on the nodes and edges based on a set of rules. This allows for an efficient filling of the space that has a packing density approximately matching reality [44].

Each of these methods can be used to simulate the powder bed process, however, the selection can be based on a couple of main factors. If speed is required then a geometric approach should be used. This approach takes on the order of seconds for a full simulation whereas the DEM approach will take hours to days to simulation an identical setup. If physical accuracy is needed then it necessary to use the DEM approach. This is because it tracks the powder particles through time and space. They are subjected to the forces of nature that result in a realistic simulation, whereas the geometric approach is simply a packing problem where the particles are placed where they can fit. This can result in packing densities that are not representatives of natural occurrences.

Modeling of the blown powder DED has only been done with the DEM approach previously discussed. To apply the DEM to blown powder, the nozzle must be modeled as a boundary condition and the particles should be fed through the feed system and tracked to determine when and where they strike the melt pool.

6. Conclusion

In order to mathematically model the AM process, it is necessary to couple several distinct mathematical models. The necessary models are thermal, fluids modeling and energy input modeling. Other models that can be included, based on the user’s desires, are stress, microstructure, surface finish, and more attributes.

Some examples where this have been done are:

1. Fan and Liou [45] model heat transfer and fluid flow dynamics (VOF) for a blown powder DED AM method with a laser for Ti-64

2. Kumar et al. [46] model heat transfer and fluid flow dynamics (SPH) for a wire feed DED AM method with a laser for Ti-64

3. Lee and Zhang [47] model powder bed generation (DEM), heat transfer, fluid flow (VOF) and microstructure for powder bed AM with Nickel-based super-alloys

Acknowledgments

The support from Department of Energy Grant DE-SC0015207, and NSF grants CMMI 1625736 and EEC 1004839, and Product Innovation and Engineering, LLC., are appreciated. We also appreciate the financial support provided by the Center for Aerospace Manufacturing Technology (CAMT) at the Missouri University of Science and Technology.