This paper presents an analytical computation of temperature field evolved in a directed energy deposition process, using single-bead walls as illustrating examples. Essentially, the temperature field evolution during the deposition of a wall is computed by super-position of the temperature field generated by the laser source depositing the current bead and that induced from each of the past beads (layers). First, the transient solution to a point heat source in a semi-infinite body is applied to describe each individual temperature field. Then, to better describe temperature contribution from a past bead, a pair of virtual heat sources with positive and negative powers is assigned for each past bead to compute the temperature field under cooling. In addition, mirrored heat sources through a reflexion technique are introduced to define adiabatic boundaries of the part and to account for substrate thickness. In the end, three depositions of Ti-6AL-4V walls with different geometries and interlayer dwell times on an Optomec® laser engineering net shaping (LENS) system are used to validate the proposed analytical computation, where predicted temperatures at several locations of the substrate show reasonable agreement with the in situ temperature measurements with prediction error rate ranging from 12% to 27%. Furthermore, temperature distributions predicted by the proposed model are compared to finite element simulations. The proposed analytical computation for temperature field could be potentially used in model-based feedback control for thermal history in the deposition, which could affect microstructure evolution and other properties of the final part.

Introduction

Additive manufacturing (AM), also called three-dimensional (3D) printing, refers to a process that builds three-dimensional parts directly from complex computer aided design (CAD) files by depositing material layer by layer. Directed energy deposition (DED) is one type of AM processes, where focused energy (either an electron beam or laser beam) is used to melt material as the material is being deposited. Laser engineering net shape (LENS) and laser metal deposition-powder (LMD-p) are among the typical DED AM processes. In modeling and simulating thermo-mechanical behavior in a DED process, it often needs to compute the temperature field evolved in the deposition process since thermal history could affect part geometry as well as microstructure, material properties, residual stress, and distortion of the final part [13].

There has been a lot of research on developing finite element analysis (FEA) models for the DED process, especially the thermomechanical behavior in the deposition process [48]. Though traditionally heat loss due to convection was assumed to be negligible and thus not included in many FEA models, a few FEA models were able to capture convection in their DED simulations to further improve the agreement between the model predicted temperature fields and their measurements [8]. In addition, radiation loss could be included using the effective convective coefficient method (see Ref. [9] and references therein for detailed discussions on modeling approaches to account for radiation and convection losses). The above-mentioned FEA models, though providing high-fidelity simulations of the DED process, could be difficult to use for a model-based feedback control design that has become increasingly important for an AM process to improve precision in part geometry, process reliability, and quality of the final part. Existing model-based feedback control algorithms in general rely on lumped-parameter models rather than partial-differential (or partial-difference)-based models as in FEA. Hence, in this paper, we are motivated to develop an analytical computation of temperature field evolved in the deposition of a DED process, based on the solution to a moving point heat source. A thermal circuit network model was developed for the prediction of temperature history and part distortion [10]. Such thermal circuit network model is a part-scale model, whereas in this paper, we consider a moving source model. Essentially, the temperature field during the deposition of a part with multibead (multilayer) is computed by superposition of the temperature field generated by the laser source depositing the current bead and that induced from each of the past beads. Such idea of superposition of temperature fields induced by past beads was applied in estimating melt-pool geometry in our prior work [11], where Rosenthal's solution for the quasi steady-state temperature distribution subject to a point heat source was used. In this paper, we extend it to compute the temperature field evolution in the deposition process by adopting a transient solution to a point heat source. In addition, several modifications are made to improve the calculation of the temperature field, which include improved estimation of temperature cooling after depositing a past bead and applying a reflexion technique to account for adiabatic boundary conditions at certain surfaces of the substrate and the build.

First, the temperature contribution from each past bead is computed by assigning a pair of virtual heat sources with opposite laser powers to replace the original physical heat source depositing that bead but following the same speed and direction. The method of using a pair of heat sources with positive and negative powers was initially proposed by Rykalin [12] and Perret et al. [13] to simulate the temperature cooling after a welding process stops. Please be noted that the application of such technique was limited to depositing a single line in a welding process in the literature. To the best of our knowledge, this paper is the first one to assign a pair of virtual heat sources with positive/negative powers to simulate the temperature contribution from a past bead, i.e., to characterize thermal interaction across beads/layers in depositing a three-dimensional part in a DED process.

Second, we apply a reflexion technique to account for the thickness of the substrate as well as to define the adiabatic boundaries corresponding to the side surfaces of a wall being built. A reflexion technique was initially introduced by Rosenthal [14] to extend the solution of a point source on the surface of a semi-infinite plate to a finite plate of certain thickness. When two heat sources run in parallel, the heat flow in between is zero, leading to the adiabatic boundary condition. Then, the temperature distribution corresponding to a finite plate of thickness d can be computed by adding a fictitious moving heat source running in parallel with a distance of 2d from the real heat source, by which the adiabatic boundary condition on the lower side of the plate is achieved. Since the fictitious heat source and the real heat source are symmetric to the plane to be defined as adiabatic, the fictitious heat source is referred to as the mirrored heat source considering the adiabatic plane as the mirror. In addition to accounting for the substrate surface being adiabatic, this paper further extends the reflexion technique to introduce mirrored heat sources for the virtual heat-source pairs representing past beads to render the side surfaces of a wall adiabatic. Compared to the existing literature where a mirrored heat sources was defined for a real heat source in depositing a single line, to the best of our knowledge, this paper is the first one to extend the reflexion technique to account for adiabatic boundary conditions in depositing a multiple-bead (multilayer) part in a DED process.

To validate the proposed analytical computation of temperature field, three depositions of Ti-6AL-4V walls with different geometries and interlayer dwell times on an Optomec LENS system are used as illustrating examples. Specifically, predicted temperature histories at several locations on the substrate are compared to the in situ temperature measurements through thermocouples. Preliminary results on the experimental validation were submitted to the 2018 Dynamic Systems and Control Conference [15]. In this manuscript, we have added the comparison of the proposed model predicted temperature field with that simulated from the FEA software netfabblocalsimulation by Autodesk.

Computation of Temperature Field in Deposition of a Single-Bead Wall

In this section, we illustrate the analytical calculation of the temperature field in the deposition process through the deposition of a single-bead wall. Section 2.1 introduces the Rosenthal's solution on the quasi steady-state temperature distribution for a point heat source moving on a semi-infinite plate as well as the transient temperature solution for a moving point source. Section 2.2 describes the process of computing the temperature field of a multibead/layer DED part, where the temperature at a given query point equals to the initial temperature plus the summation of temperature contributions from a set of heat sources including virtual heat sources assigned to all past beads. A single-bead wall is used as the illustrating example in this section. In Sec. 2.3, we improve the calculation of temperature contribution from a past bead by introducing a pair of virtual heat sources with opposite laser power to characterize the cooling effect occurred after the deposition of that past bead is finished. In Sec. 2.4, we further improve the calculation of the temperature field by (i) introducing mirrored heat sources below the substrate to define the bottom surface of the substrate to be adiabatic and (ii) introducing mirrored heat sources next to the front and back sides of the single-bead wall to define the front and back sides of the wall to be adiabatic.

Rosenthal's Solution.

We first introduce Rosenthal's solution [14]. Consider the coordinate systems defined in Fig. 1, where (x,y,z) denote the fixed coordinate. Assume that a point heat source travels with a constant speed v along the x-axis. Let (x1,y1,z1) denote the moving coordinate associated with the traveling laser heat source. Then Rosenthal's solution for the quasi steady-state temperature distribution produced by a steady-state point heat source traveling on the surface of a semi-infinite plate is given as follows [14]: 
Tss(x,y,z)=T0+q2πkRev(ξ+R)2a
(1)

where Tss(x,y,z) denotes the quasi-steady state temperature of a point of interest at coordinate (x,y,z), T0 denotes the initial temperature, q denotes the net power of the (laser) heat source (for a laser power Q with laser efficiency η, q=ηQ), k denotes the thermal conductivity, and a denotes the thermal diffusivity. Since the moving heat source is traveling at a speed of v, x1=vt. Let ξ denote the x coordinate of the point of interest in the moving coordinate (x1,y1,z1), then ξ=xx1. The parameter R denotes the distance from the point of interest to the laser heat source, i.e., R=[ξ2+(yy1)2+(zz1)2]1/2.

In the process of heat flow in arc welding, the stage before the quasi-stationary state, whose temperature field is given in Eq. (1), is referred to as the heat saturation process (transient state) during which the temperatures in the field, moving together with the heat source, continue to rise [12]. The temperature solution of a transient moving heat source in a semi-infinite body can be written as follows [12]: 
Ttr(x,y,z,t)T0=(Tss(x,y,z)T0)Γ(v,R,t)
(2)
where t denotes the time elapsed, Tss(x,y,z) denotes the quasi-steady state temperature at the same point given in Eq. (1), and Γ(v,R,t) denotes a transient transform function (or the coefficient of heat saturation), given as follows [13,16]: 
Γ(v,R,t)=12(Aexp(vRa)+B)
(3)

where A=(1erf(vt+R/2at)), and B=(1+erf(vtR/2at)), with erf() denoting the error function. It is straightforward to verify that when t = 0, Ttr(x,y,z,t)=T0, and limtTtr(x,y,z,t)=Tss(x,y,z) as the transient transform function approaches to unity in the stationary state.

It should be noted that either the Rosenthal's solution for the quasi-steady state temperature in Eq. (1) or the transient form of the temperature solution were derived using the linear heat equation, which was reduced from the original nonlinear heat equations by assuming temperature-independent material properties. For the linear heat equations, the superposition principle holds for the temperature solution and such superposition of temperature fields is applied in the rest of the paper.

Compute Temperature Field for a Multibead Deposition.

Consider a single-bead thin-wall structure, shown in Fig. 2, as an illustrative example. In Fig. 2, the fixed coordinate system is defined, with Y-axis omitted since the thin wall can be considered as two-dimensional with y = 0. Let j denote the layer index, which starts from j = 1 at the bottom layer and increases as the layer goes up. The deposition direction switches every layer as shown in Fig. 2. There could be interlayer dwell time at the end of each layer before depositing the next layer.

As illustrated in Fig. 2, consider a point of interest denoted by (x,y,z) at an arbitrary time instant t. Without loss of generality, assume that at time t, the jth layer is under deposition. Let t10=0,t20,,tj0 denote the time instants when the deposition of the first, second, , and jth layer starts. Note that in terms of the laser speed (assuming vj for depositing the jth layer), the length of the wall (L), and the interlayer dwell time, the time instants t20,,tj0 can be easily calculated. In addition, it is easy to determine the layer index j such that tj0t<tj+10. Then the temperature at the point of interest is affected by the laser heat source depositing the jth bead as well as the temperature contribution due to each of the past j − 1 beads (bead i=1,,j1). Corresponding to each past bead i, we assign a virtual heat source i to replace the original laser heat source as soon as it finishes depositing bead i. This virtual heat source i has the same power as the original laser heat source and continues to travel with the same speed along the same direction after the laser heat source ceases to deposit bead i. Then the temperature at the point of interest at time t can be computed by superposition of the temperature contribution from the laser source depositing the jth bead and temperature contributions from j − 1 virtual heat sources (corresponding to j − 1 past beads) as follows: 
T(x,y,z,t)=T0+ΔTjR+i=1j1ΔTiV
(4)

where ΔTjR denotes the temperature contribution from the real laser heat source depositing the jth bead and ΔTiV denotes the temperature contribution from the virtual heat source i. Here, we separate the temperature contributions from the real and virtual heat sources in Eq. (4) because in Sec. 2.3, the two types of temperature contributions will be modified in different ways.

Let (xjR,yjR,zjR) denote, at time t, the coordinates of the real heat source depositing the jth layer with laser power qj and speed vj, then by Eqs. (1) and (2) 
ΔTjR(x,y,z,t)=qj2πkRjevj(ξj+Rj)2aΓ(vj,Rj,ttj0)Γ(vj,Rj,ttj0)=12[Ajexp(vjRja)+Bj]
(5)
and 
Aj=1erf(vj(ttj0)+Rj2a(ttj0)),Bj=1+erf(vj(ttj0)Rj2a(ttj0))
(6)
with ξj=xxjR if j is odd or ξj=xjRx if j is even, and Rj=[ξj2+(yyjR)2+(zzjR)2].
Similarly, let (xiV,yiV,ziV) denote the coordinates of the virtual heat source i at time t in the fixed coordinate system (it is not difficult to compute such coordinates), then we have 
ΔTiV(x,y,z,t)=qi2πkRievi(ξi+Ri)2aΓ(vi,Ri,tti0)
(7)

where ξi=xxiV if i is odd or ξi=xiVx if i is even, and Ri=[ξi2+(yyiV)2+(zziV)2]1/2. Note that the determination of the coordinates of the ith virtual heat source, (xiV,yiV,ziV), should account for the interlayer dwell time during which the virtual heat source keeps traveling at the same speed as the original laser heat source.

Improve Temperature Calculation Due to Past Beads by Better Characterizing Cooling Effect.

In this section, we improve the calculation of ΔTiV, which denotes the temperature contribution from each past bead i, by better characterizing the cooling effect after the real laser heat source finishes depositing bead i and is turned off or moved to a new location. The simulation of the evolution of the temperature after a welding process is finished was first discussed by Rykalin [12] and later by Perret et al. [13]. As illustrated by Fig. 3, after the real heat source i with a constant power qi and speed vi is turned off at point B, the process of temperature cooling can be simulated by superposing two processes: the process for a virtual heat source i that represents the continuation of the original real heat source, and the process for a virtual heat sink starting at point B with equal but negative power value (−qi). In this paper, we refer such pair of virtual heat source/sink as positive virtual heat source i (denoted by qi+) and negative virtual heat source i (denoted by qi). In this section, such pair is used to replace each virtual heat source in Fig. 2 used in Sec. 2.2. Please be noted that the pair of virtual heat sources is located in the same position and move with the same speed and direction. They are plotted separately, rather than overlapped, in Fig. 3 for a better illustration.

Let ti0 denote the deposition time corresponding to point A in Fig. 3, and let L denote the track length (distance from point A to point B), then the time when the real heat source i with qi and speed vi stops acting at point B is ti0+L/vi, which is the time when the negative heat source (heat sink) qi starts. By Eqs. (1) and (2), for a point of interest (x,y,z), the temperature contribution due to the positive heat source qi+ can be computed as follows: 
ΔTiV(x,y,z,t|qi+)=qi2πkRievi(ξi+Ri)2aΓ(vi,Ri,tti0)
(8)

Note that the positive virtual heat source qi+ is considered as a continuation of the real heat source, and thus the time elapsed used in computing the temperature field due to qi+ in Eq. (2) is counted from ti0.

The temperature contribution due to the negative heat source qi, which starts at ti0+L/vi, can be computed as follows: 
ΔTiV(x,y,z,t|qi)=qi2πkRievi(ξi+Ri)2aΓ(vi,Ri,t(ti0+Lvi))
(9)

To this end, after bead i is finished depositing, the temperature at a point of interest contributed by the pair of positive and negative virtual heat sources i is given as follows:

For t>ti0+Lvi 
ΔTiV(x,y,z,t)=ΔTiV(x,y,z,t|qi+)+ΔTiV(x,y,x,t|qi)=qi2πkRievi(ξi+Ri)2aΓ(vi,Ri,tti0)+qi2πkRievi(ξi+Ri)2aΓ(vi,Ri,t(ti0+Lvi))
(10)

where ξi=xxiV if i is odd or ξi=xiVx if i is even, and Ri=[ξi2+(yyiV)2+(zziV)2]1/2.

Since qi+ starts earlier at ti0 whereas qi starts at ti0+L/vi, when the cooling starts, the temperature field due to qi+ is more close to the quasi-steady state than the temperature field due to qi, which is still in transient state and thus the absolute value of the corresponding temperature in the field continues to grow. Further, the absolute value of the temperature contribution due to qi is lower than that due to qi+, which leads the summation of the temperature contributions from the pair to be positive but decrease with the increase of the time, mimicking the temperature cooling. After a sufficient long time, the summation of the temperature contributions from the pair will approach zero since the temperatures due to qi+ and qi will reach to steady-state and they are the same in magnitude but opposite in sign. This corresponds to a complete cooled down where the temperature contribution from a past bead i equals zero.

Improve Temperature Calculation by Defining Adiabatic Boundaries at Surfaces of Substrate and Part.

Recall that in both the Rosenthal's quasi-steady state solution (1) and the transient solution (2), the plate on which the point heat source is moving is assumed semi-infinite. In practice, this assumption is not satisfied due to the finite geometry of the build part and the substrate. To simulate adiabatic boundary condition analytically, Rosenthal introduced the reflexion technique in Ref. [14]. As in Fig. 4, a pair of point heat sources is running in parallel with equal distance d from the plane AB. Then the heat flow across plane AB is equal to zero (i.e., adiabatic) if these two heat sources share the same power. These two heat sources are referred to as the original heat source and the mirrored heat source, with plane AB acting as a mirror. Consequently, the overall temperature field incorporating the adiabatic boundary condition can be approximated by initial temperature plus summation of temperature contributions from these two heat sources (the real one and its mirrored one), each of which can be calculated by Eq. (2).

In this paper, we define the bottom surface of the substrate underneath the built part, as well as the front and back side surfaces of the built part to be adiabatic. Figure 5 shows the side view of a single-bead wall with one-layer deposited. Corresponding to the real heat source, mirrored real heat sources (where bt, fr, and bk representing bottom, front and back mirrored heat sources) are defined to render the bottom surface of the substrate as well as the front and back side surfaces of the wall to be adiabatic (see the adiabatic boundaries marked in solid line segments in Fig. 5). For example, the real heat source and the “bt” mirrored heat source are running in parallel along the direction x (perpendicular to the yz plane shown in Fig. 5) with equal distance ziR+d from the bottom surface of the substrate, where ziR denotes the z-coordinate of the real heat source i and d denotes the thickness of the substrate.

When multiple layers are deposited, we extend the reflexion technique to further introduce mirrored heat sources for the virtual heat-source pairs representing past beads/layers to render the (extended) side surfaces of the multilayer wall adiabatic, i.e., following the moving of the pair of virtual heat sources replacing the real heat source in Fig. 5, the extended virtual side surfaces are defined to be adiabatic. In summary, the reflexion technique is applied to all heat sources, including the real heat source as well as the positive and negative virtual heat sources for all past beads.

In this paper, we omit the convective and radiative heat transfers during a DED process. In addition, we only account for adiabatic boundary conditions on two sets of planes in deposition of a wall: (1) bottom surface of the substrate, and (2) front and back sides of the wall (the ones along the deposition direction), assuming that adiabatic boundary conditions at other surfaces (e.g., the top surface of the substrate, and the two short edges of the wall) would have less dominant (or negligible) effect on the temperature field.

Consequently, ΔTjR in Eq. (4), originally only accounting for the temperature contribution from the real heat source, are now redefined in terms of contributions from both the real heat source and its mirrored heat sources as follows: 
ΔTjR=m=o,bt,ft,bkΔTjR(m)=m=o,bt,fr,bkqj2πkRj(m)evj(ξj+Rj(m))2aΓ(vj,Rj(m),ttj0)
(11)
where ΔTjR(o),ΔTjR(bt),ΔTjR(fr), and ΔTjR(bk) denote the temperature contribution from the original real heat source (not mirrored), mirrored heat source in terms of the bottom surface of the substrate, mirrored heat source in terms of the front side of the wall, and mirrored heat source in terms of the back side of the wall, respectively (see Fig. 5). Here, the front- and back-sides of the wall correspond to the wall sides along the x direction. Note that ξj=xxjR if j is odd or ξj=xjRx if j is even, and yjR=0 for a single-bead wall, then we have 
Rj(o)=[ξj2+y2+(zzjR)]1/2
(12)
 
Rj(bt)=[ξj2+y2+(z+2d+zjR)]1/2
(13)
 
Rj(fr)=[ξj2+(y+w)2+(zzjR)]1/2
(14)
 
Rj(bk)=[ξj2+(yw)2+(zzjR)]1/2
(15)

where d and w denote the thickness of the substrate and wall width, respectively.

Similarly, ΔTiV(qi+) and ΔTiV(qi) (which are shorthand notations for ΔTiV(x,y,z,t|qi+) and ΔTiV(x,y,z,t|qi)) in Eqs. (8) and (9) are redefined as follows to incorporate temperature contributions from their mirrored heat sources: 
ΔTiV(qi+)=m=o,bt,fr,bkΔTiV(qi+,m)=m=o,bt,fr,bkqi2πkRi(m)evi(ξi+Ri(m))2aΓ(vi,Ri(m),tti0)
(16)
 
ΔTiV(qi)=m=o,bt,fr,bkΔTiV(qi,m)=m=o,bt,fr,bkqi2πkRi(m)evi(ξi+Ri(m))2aΓ(vi,Ri(m),t(ti0+Lvi))
(17)

where ΔTiV(qi+,o),ΔTiV(qi+,bt),ΔTiV(qi+,fr) and ΔTiV(qi+,bk) denote the temperature contribution from the original positive virtual heat source i (not mirrored), mirrored positive virtual heat source i in terms of the bottom surface of the substrate, mirrored positive virtual heat source i in terms of the front side of the wall, and mirrored positive virtual heat source i in terms of the back side of the wall, respectively. Similarly, ΔTiV(qi,o),ΔTiV(qi,bt),ΔTiV(qi,fr) and ΔTiV(qi,bk) are defined to represent the temperature contributions from the negative virtual heat source i and its mirrored heat sources. The distance from the point of interest to each mirrored heat source Ri(m) can be defined in a similar fashion as in Eqs. (12)(15), simply replacing the coordinate (xjR,yjR=0,zjR) in Eqs. (12)(15) with the coordinate of the virtual heat source i(xiV,yiV=0,ziV).

In summarizing the modeling in Sec. 2.22.4, the temperature field evolved in a DED process is computed by superposition of temperature contributions from the real heat source depositing the current bead, the positive/negative virtual-heat-source pair representing each of the past beads, and their mirrored heat sources used to satisfy the adiabatic boundary conditions at the bottom surface of the substrate and sides of the wall. A list of the assumptions used in the modeling work in Sec. 2.22.4 is summarized as follows:

Assumption 1. Assumptions for Rosenthal's solution to hold [14].

Assumption 2. The superposition principle holds for the temperature distributions due to current heat source, virtual heat sources representing past beads, and their mirrored heat sources.

Assumption 3. Adiabatic boundary conditions at other surfaces than the substrate's bottom surface beneath the wall as well as the front- and back-side of the wall would have less dominant (or negligible) effect on the temperature field.

Assumption 4. Radiation and convection losses are omitted in the modeling.

Validate Model Predicted Temperature Histories Using In Situ Measurements From Thermal-Couples for Single-Bead Walls

Recall that Sec. 2.2 introduces the baseline calculation of the temperate field, Sec. 2.3 improves the temperature calculation by better characterizing the cooling effect after depositing each past bead, and Sec. 2.4 further improves the temperature calculation by adding adiabatic boundary conditions at the substrate bottom surface and at the front- and back-sides of the wall. In the following, we compare the model prediction from each afore-mentioned step (with increasing complexity) to the temperature measurements, to have a better understanding of the effect of increasing model complexity on the resulting accuracy in temperature prediction. Specifically, we consider the following four models with increasing complexity:

  • Model 1 (M1): Baseline model given by Eqs. (4)(7)

  • Model 2 (M2): Model 1 + updated cooling effect (replace Eq. (7) by Eq. (10))

  • Model 3 (M3): Model 2 + adiabatic boundary condition at bottom surface of the substrate

  • Model 4 (M4): Model 3 + adiabatic boundary conditions at front- and back-sides of the wall

In this section, three Ti-6AL-4V single-bead walls in Ref. [8], which were deposited on an Optomec LENS MR-7 machine, are used to validate the analytical computation of temperature field evolution given by the above four models. Figure 6 shows the images of the three single-bead walls. Wall 1 and wall 3, each of which has 62 layers, were directly deposited on a substrate with the dimension of 76.2mm×25.4mm×6.4mm. There was no interlayer dwell time in deposition of wall 1 whereas a 20 s interlayer dwell time was used to deposit wall 3. Wall 2 was built by depositing another 62 layers on top of wall 1 after it was completely cooled down.

To collect in situ temperature measurements during the deposition process, three Omega GG-K-30 type K thermocouples (with response time ∼0.003 s) were deployed on the substrate or on the part (see Fig. 10 of Appendix  A for the locations of the thermocouples). Specifically,

  • TC1 was located at the bottom surface of the substrate, with coordinate (L/2,0,6.4mm)

  • TC2 was deployed on the top surface of the substrate and next to the single-bead wall, with coordinate (L/2,6.3mm,0)

  • TC3 was welded to the side of wall 1 after it was completely cooled down, with coordinate (L/2,w/2,10.5mm)

Where the above coordinates are given in term of the fixed coordinated system defined in Fig. 2, and L and w denote the length and width of the respective wall. Note that in Ref. [8], measurements from TC1 were plotted for wall 1 (see Fig. 8 of Ref. [8]); measurements from TC1 and TC3 were plotted for wall 2 (see Fig. 9 of Ref. [8]); and measurements from TC2 were plotted for wall 3 (see Fig. 10 of Ref. [8]). In this paper, these measurements on temperature history are used for comparison to predictions of the proposed models. More details regarding the experimental setup and thermocouples used in the experiment can be found in Ref. [8].

The goodness of fit for the predicted temperature history compared to the in situ temperature measurements from a thermocouple is characterized using the average prediction error rate e and the root-mean-square error RMSE defined as follows [8]: 
e=i=1n|(Tpred(i)Tmea(i))/Tmea(i)|n
(18)
 
RMSE=i=1n(Tpred(i)Tmea(i))2/n
(19)

where n denotes the total number of measurements with respect to time from the start to the end of the deposition. Tpred(i) (unit: °C) and Tmea(i)(unit: °C) denote the model predicted temperature and thermocouple measurement at each time instant i, respectively.

Table 1 lists the process parameter values used in the deposition and the measurements of the final parts. Table 2 provides simulation parameters and material property parameters for Ti-6AL-4V, where the ambient temperature T0 = 292 K, the melting temperature Tm = 1923 K, and the thermal diffusivity a = 2.48 × 10−6 are used; further, a constant thermal conductivity (k = 6.7 W/m K) is used in the simulation due to the assumption of temperature-independent material properties to allow the superposition of temperature distributions. The laser transmission efficiency η = 30% is used in the simulation. It should be noted that the laser transmission efficiency η = 45% was used in Ref. [8] for FEA simulations compared to measurements; this value of laser transmission efficiency (η = 45%) was determined using an inverse simulation method given in Ref. [17]. In this paper, a reduced value for laser transmission efficiency (lower transmission efficiency than η = 45%) is used to adjust for heat losses due to convection and radiation, noting that they are not considered in the proposed model. An optimization on η can be performed in the future work for a more systematic selection of η.

Figure 7 shows the comparison between the in situ temperature measurements from the thermocouples (TC1, TC2, and TC3) and the predicted temperature histories by model 1–4 at the respective location of the thermocouples. We can see that the predicted temperature drops after the cooling effect is better modeled by adding the negative virtual heat source, noting that the temperature history curve for model 2 is lower than that for model 1. Further, adding an adiabatic boundary condition in general drives up the temperature, noting that the temperature curve for model 3 is higher than that of model 2, and the temperature curve for model 4 is higher than that of model 3.

For both model predictions and measurements across almost all cases in Fig. 7, there is a visible high frequency oscillation associated with each temperature curve, where each cycle of the oscillation corresponds to the temperature variation within a layer of the wall when the laser travels from one end to the other end to build the layer. Specifically, corresponding to the laser scanning speed of 8.5 mm/s and wall length of 39.2 mm, it takes ∼4.6 s for laser to travel from one end of the wall to the other end. Due to the laser scan path, the temperature starts high at one end and then cools off when the laser moves to the other end. Such temperature change in the period of ∼4.6 s leads to the high-frequency oscillation. The larger temperature oscillations in Figs. 7(b) and 7(c) than in Fig. 7(a) are due to the locations of the thermocouples. Note that TC1 is located at bottom surface of the substrate, whereas TC2 is on the top surface and thus closer to the wall deposition, and TC3 is welded on wall 1 and thus close to the deposition of wall 2 as well. In addition, due to TC2 and TC3's closer proximity to deposition, their temperatures are affected by convection heat loss more than TC1 that is located further away from the deposition.

The magnitude of the high-frequency oscillation decreases as the layer goes up since the location of each thermocouple gets further away from the laser heat source as the layer goes up, which makes the temperature variation within a layer less significant. If removing the high-frequency oscillation from each temperature curve, we can see that the temperature increases to reach its peak after certain time (after building a number of layers) and then gradually decreases. This is consistent with the expectation that heat accumulates as the layer goes up, reaches to a steady-state, and then drops as the laser heat source moves further away.

Figure 7 shows that model 4 is able to capture the general shape of the measured temperature history in each case. Specifically, for wall 1 (see Fig. 7(a)), compared to the TC1 measurements, model 4 has an average prediction error rate e = 11.18% and root-mean-square error RMSE = 47.8 °C. For wall 2 (see Fig. 7(b)), compared to the TC1 measurements, model 4 has e = 27% and RMSE = 71.7 °C; and compared to the TC3 measurements, model 4 has e = 11.86% and RMSE = 76.04 °C. For wall 3 (see Fig. 7(c)), compared to TC2, model 4 has e = 13.37% and RMSE = 25.4 °C. In summary, the prediction error rate e by model 4 is less than 15% in each case except for TC1 in wall 2 for which the prediction error rate e has doubled. Note that wall 2 is built on top of wall 1, and thus wall 2 can be viewed as being built on a larger substrate with increased mass. In this case, mode 4's temperature prediction at the bottom surface of the substrate, though captures the general shape of TC1, is significantly lower than TC1 with e = 27.01%, possibly suggesting that η = 30% (which is set to adjust for not considering heat losses due to convection and radiation in the proposed models) may need to be retuned when there is a significant change in substrate.

The simulation run time used to compute each temperature history in Fig. 7 is given in Table 3. All simulations were conducted using matlab® R2016a on a laptop with Intel® Core™ i5-3210M CPU at 2.50 GHz and 8.00 GB RAM. The sampling time of the simulation was 50 ms, which was selected to match the sampling time of the thermocouple measurement. The run time for wall 2 was counted from the start of depositing the top 62 layers. The run time for wall 3 accounted for the interlayer-dwell time.

Comparison to Finite Element Analysis Predicted Temperature Distribution

Note that the experimental validation through thermocouple measurements is performed at a few locations. In this section, taking wall 1 as an illustrative example, we compare the predicted temperature distributions by the proposed models versus that simulated from the FEA software netfabblocalsimulation by Autodesk. Netfabb Local Simulation, whose predecessor was CUBES by Pan Computing LLC (State College, PA), is a nonlinear Newton–Raphson based code and solves the discretized transient thermal equilibrium of a component subject to a body heat source. Modeling details used in the Netfabb Local Simulation (and its predecessor CUBES) including the activation strategy and boundary condition application method can be found in Ref. [7].

In Ref. [8], FEA analysis using CUBES was performed to compute the temperature history in each of the three single-bead wall depositions discussed in Sec. 3, and the predicted temperature history at TC1–TC3 was compared to the measurements from the respective thermocouple. Specifically, FEA with two different convection models were run in Ref. [8]: (1) A free convection model assuming a uniform coefficient of convection on all surfaces equal to 10 W/m2/°C and (2) a forced convection model that incorporates convection measurements. Readers are referred to Figs. 810 of Ref. [8] for the comparison among the measurements, temperature history predicted by the FEA with the free convection model, and that by the FEA with the forced convection model. In this paper, we compare our model predicted temperature distributions with the FEA simulation with no convection (since convection is not considered in the proposed models in this paper) and the same laser transmission efficiency (η = 30%). Readers are referred to Fig. 5 and Sec. 4 of Ref. [8] for the finite element mesh and other FEA simulation details used for the single-bead wall depositions. Due to the limit of the space, we do not repeat them in this paper.

For the sake of easier illustrating the point-to-point comparison of temperature distribution predicted by the proposed models versus that of the FEA, we consider the initial/start temperature at the deposition point before laser hits it and assume that the temperature at a distance ahead of the moving laser source has yet been heated up by the laser heat source. Such initial/start temperature at the deposition point essentially represents the thermal contributions due to depositions of all past beads, excluding the effect of depositing the current bead. We choose such initial/start temperature as a proxy to compare the capability of the proposed models versus FEA in predicting temperature field during the deposition process due to the following reasons: (i) such initial temperature at the deposition point, when following the laser heat source to move from layer to layer, contains both temporal and spatial information, and (ii) the calculation of such initial temperature due to the past beads follows the same procedures as illustrated in Sec. 2 excluding the temperature term due to heat source depositing the current bead. For example, consider the initial temperature at the deposition point in Fig. 2. In terms of the proposed model 1, the initial temperature at the deposition point is equal to the ambient temperature plus the summation of temperature contributions from the virtual heat sources 1–4. In terms of the proposed model 4, the calculation of temperature contributions from past beads 1–4 is updated by replacing each of the virtual heat sources 1–4 with a pair of positive and negative virtual heat sources and their mirrored heat sources.

Next, we describe how to pick the temperature value from the FEA simulation and use it to compare to the afore-mentioned model 1's or model 4' predicted initial temperature at the deposition point. Figure 8 shows a snapshot of the FEA simulation, with the temperature contour given in Fig. 8(a) and the temperature along line A given in Fig. 8(b). When the laser travels from left to right to deposit a bead, the temperature pulse along line A will travel from left to right accordingly. The temperature pulse along line A indicates that the temperature at a point to the right side of the pulse is not much affected by the laser. Let r denote the laser spot radius, r = 1.5 mm is used in this study. The circle point in the temperature curve in Fig. 8(b), with a distance s=4r=6mm ahead of the laser, is one example point whose temperature is picked to represent the initial temperature at the deposition point due to past beads.

Figure 9 shows the comparison of the predicted initial temperature at the deposition point by the proposed model 1 (represented by dashed lines), prediction by model 4 (represented by solid lines), and prediction from the FEA simulation (represented by circles). For the predicted temperatures given by the FEA, temperature of the point with s=4r=6mm along line A is picked to represent the initial temperature at the deposition point. We have explored using temperature value at a point with a different distance s to represent the initial temperature at the laser deposition point. In Appendix  B, Figs. 11 and 12 show the comparison of initial temperatures predicted by model 1 and model 4, as well as initial temperature profiles generated by the FEA using s=2r=3mm and s=6r=9mm, respectively.

For the clarity of the illustration, only a number of example layers are shown in Fig. 9. Figure 9(a) shows the initial temperature at the deposition point along the wall length for representative even layers, where the laser scans from right to left. Figure 9(b) shows the initial temperature at the deposition point at representative odd layers, where the laser scans from left to right. First, we can see that compared to the FEA predictions, the proposed model predictions capture the right trend and shape of the temperature variation following the laser scan direction. For example, when laser moves from the right end of the wall toward left in depositing an even layer of the wall, the initial temperature at the deposition point decreases from right to left since the laser has just heated up the right end during the previous pass. Similarly, when laser moves from the left end of the wall toward right in depositing an odd layer of the wall, the initial temperature at the deposition point decreases from left to right since the laser has just heated up the left end during the previous pass. The curved right ends in model 4's predictions for even layers, as well as the curved left ends for odd layers are artifacts caused by applying a pair of positive and negative virtual sources simulating the cooling effect, where the transient temperature due to the negative heat source decreases in a much faster rate than that due to the positive heat source near the end of the wall.

Figure 9 also shows that when only a few layers being built (e.g., at the second layer and third layer), predicted temperatures from model 1 are close to the FEA predictions, indicating that Rosenthal's solution (or superposition of Rosenthal's solutions of past beads) can be a pretty good prediction on the temperature field. However, as the layer goes up, model 1 starts to underestimate the temperatures compared to the FEA predictions and such underestimation grows significantly with the increase of the layer (model 1's predictions are ∼50% lower than the FEA predictions at the 61st or 62nd layer). By considering the adiabatic boundary conditions, model 4 significantly overpredicts temperatures at the first several layers compared to the FEA predictions (∼45% over-prediction at the second or third layer), but the overprediction decreases as the layer goes up, with less than 15% prediction error at wall distance from 5 mm to 30 mm for the 61st layer and at wall distance from 35 mm to 10 mm for the 62nd layer.

Discussions

The main objective of this paper is to conduct a proof-of-concept study of a simplified analytical computation of temperature field in a DED process, with the ultimate goal to use such simplified model of temperature field for real-time optimization and feedback control. Noting that the proposed model has an analytical expression, it is more amenable for designing model-based optimization and feedback control than FEA models. We understand that such analytical model is based on many assumptions and simplifications (especially the assumptions to allow the Rosenthal's solution to hold) and hence try to evaluate how the predictions of the proposed analytical computation are compared to FEA simulations and experimental measurements. In comparison to experimental measurements, the model prediction error rate could vary between 12% and 27%. As far as control systems design, though higher modeling accuracy often leads to better control performance, complexity of the model itself can often affect greatly the easiness and robustness of the implementation of the control system in practice. Very often, reduced-order or linearized models are adopted for control system design and implementation in practice as long as the modeling uncertainties can be quantified, with the expectation of applying robust control algorithms to provide robustness with respect to modeling uncertainties. Note that as discussed in Sec. 2, heat losses due to convection and radiation are not accounted in the models in this paper, future work will investigate improving the modeling accuracy by incorporating such heat losses.

Very simple part geometry such as single-bead walls is used in this paper as illustrating examples. In our prior work [11], computation of initial temperature at the deposition point with a steady-state version of model 1 of this paper was applied in the prediction of build height for an L-shaped structure that has a one-bead leg and a three-bead leg as well as for a patch build. This indicates that the proposed modeling work in this paper can be scaled to builds with multiwall beads, though experimental validation for parts with more complex geometries needs to be investigated in the future. In addition, since Rosenthal's solution is used, the models in this paper are limited to constant laser scanning speed. Future work will also investigate the possibility of generalizing the approach of placing virtual heat sources with different speeds to relax such limitation.

Conclusion

In this paper, we present an analytical computation of temperature field in a DED process, by superposition of the temperature field generated by the laser source depositing the current bead and the temperature fields individually induced from each of the past beads. A pair of positive and negative virtual heat sources is assigned to each past bead to characterize the temperature contribution from each past bead, and mirrored heat sources are introduced to account for adiabatic boundary conditions at certain surfaces of the substrate and the part. Three single-bead walls deposited using an LENS machine with Ti-6AL-4V are used to validate the proposed analytical models. The model predicted temperature histories at several locations of the substrate show reasonable agreement with the in situ temperature measurements, with the average prediction error rate ranging from 12% to 27%. In the end, the temperature distributions predicted by the proposed models are compared to FEA simulations, where the initial temperature at the deposition point (representing temperature contributions from past beads) is used as a proxy for comparison of temperature field as laser moves from layer to layer. The temperature distributions predicted by the proposed models are able to capture the shape and trend of the FEA predictions, where model 1' predictions have a better agreement with FEA predictions at lower layers, whereas model 4's predictions tend to have a better agreement with FEA predictions at higher layers.

Note that the proposed models for computing temperature field have not accounted for heat loss due to convection as considered in Ref. [8], neither accounted for radiation loss, but instead, a reduced value for laser transmission efficiency is used to compensate such model inadequacy. The future modeling work should incorporate heat losses of convection and radiation. Though the analytical computation of temperature field presented in this paper is based on many assumptions, we believe that it provides the initial proof of concept and it could potentially be used for the development of model-based feedback control for thermal history in deposition, which could largely affect microstructural evolution and mechanical properties of the final part.

Acknowledgment

The authors would like to thank Mr. Jeff Irwin from Autodesk to help run the Netfabb Local Simulation for this paper.

Funding Data

  • National Science Foundation (Grant No. CMMI 1563271).

Appendix A: In Situ Temperature Measurements by Thermocouples

As described in Ref. [8], in situ temperature was measured at several locations of the substrate using Omega GG-K-30 type K thermocouples. See Fig. 10 for the schematic plot of the locations of the measurements. TC1 was located at the bottom surface of the substrate, and TC2 was located on the top surface near the wall. TC3 was welded to the face of wall 1 before the second wall was built.

Appendix B: Finite Element Analysis Prediction on Initial Temperature at Deposition Point

Figures 11 and 12 show the comparison of initial temperatures predicted by model 1, model 4, and FEA, which uses temperature at s=2r=3mm or s=6r=9mm ahead of the laser to represent the initial temperature at the laser deposition point. From Fig. 8(b), we can see that the point 3 mm ahead of the laser is inside the temperature pulse generated by the moving laser and thus its temperature is not a good representative of the initial/start temperature at the laser deposition point. As expected, Fig. 11 confirms that by considering the 3 mm offset, the initial temperature profiles predicted by FEA at low layers (e.g., second or third layers) are ∼30% higher than the predicted initial temperatures by model 1. In addition, Fig. 11 shows that by considering the 3 mm offset, the FEA predicted initial temperatures along the wall distance at each example layer, though follow a similar shape as the FEA predicted initial temperature profiles using the 6 mm offset (shown in Fig. 9), demonstrate small ripples affected by the temperature pulse of the moving laser.

In comparison, the point 9 mm ahead of the laser is outside the temperature pulse generated by the moving laser (as shown in Fig. 8(b)) and thus its temperature is much less affected by the moving laser. At each example layer, compared to the FEA predicted initial temperature profiles using s=4r=6mm in Fig. 9, the FEA predicted initial temperature profiles using s=6r=9mm show a similar shape but with a much smoother curve. In addition, similar as demonstrated in Fig. 9, at low layers (second or third layer), the FEA predicted initial temperatures along the wall distance have good agreement with the predicted initial temperatures by model 1. However, with the increase of the layer, model 1 starts to underestimate the temperatures and such underestimation grows with the increase of the layer. By considering adiabatic boundary conditions, model 4 significantly overpredicts temperatures compared to the FEA predictions at the second or third layer, but such overprediction decreases as the layer goes up, with a much reduced prediction error at the 61st or 62nd layers.

References

References
1.
Kelly
,
S.
, and
Kampe
,
S.
,
2004
, “
Microstructural Evolution in Laser-Deposited Multilayer Ti-6AL-4V Builds—Part I: Microstructural Characterization
,”
Met. Mater. Trans. A
,
35
(
6
), pp.
1861
1867
.
2.
Kelly
,
S.
, and
Kampe
,
S.
,
2004
, “
Microstructural Evolution in Laser-Deposited Multilayer Ti-6Al-4V Builds—Part II: Thermal Modeling
,”
Met. Mater. Trans. A
,
35
(
6
), pp.
1869
1879
.
3.
Irwin
,
J.
,
Reutzel
,
E. W.
,
Michaleris
,
P.
,
Keist
,
J.
, and
Nassar
,
A. R.
,
2016
, “
Predicting Microstructure From Thermal History During Additive Manufacturing for Ti-6AL-4V
,”
ASME J. Manuf. Sci. Eng.
,
138
(
11
), p.
111007
.
4.
Hodge
,
N. E.
,
Ferencz
,
R. M.
, and
Solberg
,
J. M.
,
2014
, “
Implementation of a Thermomechanical Model for the Simulation of Selective Laser Melting
,”
Comput. Mech.
,
54
(
1
), pp.
33
51
.
5.
Yang
,
Y.-P.
, and
Babu
,
S. S.
,
2010
, “
An Integrated Model to Simulate Laser Cladding Manufacturing Process for Engine Repair Applications
,”
Weld. World
,
54
(
9–10
), pp.
R298
R307
.
6.
Marimuthu
,
S.
,
Clark
,
D.
,
Allen
,
J.
,
Kamara
,
A. M.
,
Mativenga
,
P.
,
Li
,
L.
, and
Scudamore
,
R.
,
2013
, “
Finite Element Modeling of Substrate Thermal Distortion in Direct Laser Additive Manufacture of an Aero-Engine Component
,”
Proc. Inst. Mech. Eng., Part C
,
227
(
9
), pp.
1987
1999
.
7.
Michaleris
,
P.
,
2014
, “
Modeling Metal Deposition in Heat Transfer Analysis of Additive Manufacturing Processes
,”
Finite Elem. Anal. Des.
,
86
(
1
), pp.
51
60
.
8.
Heigel
,
J. C.
,
Michaleris
,
P.
, and
Reutzel
,
E. W.
,
2015
, “
Thermo-Mechanical Model Development and Validation of Directed Energy Deposition Additive Manufacturing of Ti–6Al–4V
,”
Addit. Manuf.
,
5
, pp.
9
19
.
9.
Gouge
,
M.
, and
Michaleris
,
P.
,
2017
,
Thermal-Mechanical Modeling of Additive Manufacturing
, Science & Technology, Elsevier, Amsterdam, The Netherlands.
10.
Peng
,
H.
,
Go
,
D. B.
,
Gong
,
S.
,
Shankar
,
M. R.
,
Gatrell
,
B. A.
,
Budzinski
,
J.
,
Ostiguy
,
P.
,
Attardo
,
R.
,
Tomonto
,
C.
,
Neidig
,
J.
, and
Hoelzle
,
D.
, 2016, “
Part-Scale Model for Fast Prediction of Thermal Distortion in DMLS Additive Manufacturing—Part 1: A Thermal Circuit Network Model
,”
Solid Freeform Fabrication Symposium
, Austin, Texas, Aug. 8–10, pp. 297–382.
11.
Li
,
J.
,
Wang
,
Q.
,
Michaleris
,
P. (Pan)
,
Reutzel
,
E. W.
, and
Nassar
,
A. R.
,
2017
, “
An Extended Lumped-Parameter Model of Melt–Pool Geometry to Predict Part Height for Directed Energy Deposition
,”
ASME J. Manuf. Sci. Eng.
,
139
(
9
), p.
091016
.
12.
Rykalin
,
N. N.
,
1960
,
Calculation of Heat Processes in Welding
,
Moscow, Russia
.
13.
Perret
,
W.
,
Schwenk
,
C.
, and
Rethmeier
,
M.
,
2010
, “
Comparison of Analytical and Numerical Welding Temperature Field Calculation
,”
Comput. Mater. Sci.
,
47
(
4
), pp.
1005
1015
.
14.
Rosenthal
,
D.
,
1946
, “
The Theory of Moving Sources of Heat and Its Application to Metal Treatments
,”
Trans. ASME
,
68
(
8
), pp.
849
866
.
15.
Li
,
J.
,
Wang
,
Q.
, and
Michaleris
,
P.
,
2018
, “
Towards Computational Modeling of Temperature Field Evolution in Directed Energy Deposition Processes
,” ASME Paper No. DSCC2018-8973.
16.
Cao
,
Z.
,
Dong
,
P.
, and
Brust
,
F.
,
2000
, “
Fast Thermal Solution Procedure for Analyzing 3D Multi-Pass Welded Structures
,”
Weld. Res. Counc. Bull.
,
455
(1), pp.
12
21
.
17.
Denlinger
,
E.
,
Heigel
,
J.
, and
Michaleris
,
P.
,
2015
, “
Residual Stress and Distortion Modeling of Electron Beam Direct Manufacturing Ti-6AL-4V
,”
Proc. Inst. Mech. Eng., Part B
,
229
(
10
), pp.
1803
1813
.