There is a need for the development of lumped-parameter models that can be used for real-time control design and optimization for laser-based additive manufacturing (AM) processes. Our prior work developed a physics-based multivariable model for melt–pool geometry and temperature dynamics in a single-bead deposition for a directed energy deposition process and then validated the model using experimental data from deposition of single-bead Ti–6AL–4V (or Inconel®718) tracks on an Optomec® Laser Engineering Net Shaping (LENS) system. In this paper, we extend such model for melt–pool geometry in a single-bead deposition to a multibead multilayer deposition and then use the extended model on melt–pool height dynamics to predict part height of a three-dimensional build. Specifically, the extended model incorporates temperature history during the build process, which is approximated by super-positioning the temperature fields generated from Rosenthal's solution of point heat sources, with one heat source corresponding to one bead built before. The proposed model for part height prediction is then validated using builds with a variety of shapes, including single-bead thin wall structures, a patch build, and L-shaped structures, all built with Ti–6AL–4V using an Optomec® LENSTM MR-7 system. The model predictions on average part height show reasonable agreement with the measured average part height, with error rate less than 15%.

Introduction

Additive manufacturing (AM) has the potential to produce three-dimensional (3D) parts directly from complex computer-aided design files. In the last two decades, various AM processes for metals have been developed, where parts can be built from metal wire [1] or powderized feedstock, either in a powder-feed process [24] or a powder-bed process [5,6]. Essentially, a 3D part represented by a computer-aided design file is sliced into layers, which encompass scan trajectories of the heat source. Then, a high power energy source (laser or electron beam) is used to heat and melt metal powder (or wire), which solidifies to form a dense layer. A number of such layers can then be repeated to form a complex shape. Despite of the huge potential of AM technology, it has not been widely applied in industry due to multiple challenges associated with the AM processes. One such challenge is on accuracy of part geometry, which motivates us to develop a dynamic model for melt–pool geometry that can be used for the design of model-based controller to achieve accuracy of part geometry.

There has been increasing research effort in developing analytical and numerical models for AM processes, especially in finite element (FE) modeling of AM processes [610]. Nevertheless, FE models use partial differential or difference equations, which are difficult to be used in real-time control design. The development of control-oriented models is in demand but still in its infancy. A set of studies aimed to develop empirical or statistical models using the experimental data, where linear transfer functions [1113], linear state space models [14,15], a semi-empirical nonlinear Hammerstein model [16,17], and statistical models such as analysis of variance [18] and Bayesian models [19] were developed to characterize how AM process parameters affect melt–pool temperature and height in a variety of AM processes. Another line on control-oriented models is to derive physics-based lumped-parameter models for melt–pool geometry and temperature. Based on the principles of mass, momentum, and energy conservation, Doumanidis and Kwak derived a lumped-parameter model for the melt pool geometry and temperature dynamics as a function of laser power, material feed rate, and laser scan speed, for a single-bead deposition [20]. A physics-based high-order time-varying nonlinear state space model of the melt pool thermal dynamics during laser cladding was derived by Devesse et al. [21], and then a linearization of such model was used in Ref. [22] to design a temperature feedback controller. In Ref. [23], the single-bead process model by Doumanidis and Kwak was extended to model a thin-wall structure by introducing a solidification rate in the equations of mass and momentum conservation, where the solidification rate was calculated using a simplified two-dimensional (2D) FE solution under a given laser travel speed and a given part height.

However, even for a single-bead deposition, the model of Doumanidis and Kwak did not provide an explicit mathematical expression on how laser power would affect the steady-state melt–pool height [20]; specifically, it can be shown that the steady-state melt–pool height derived from their model is independent of the laser power. In comparison, our past experimental results demonstrated that change of laser power only (by fixing other processing parameters) did affect melt–pool height and its steady-state value [24]. In Ref. [25], a similar work was performed as Doumanidis and Kwak, but assuming the powder catchment efficiency to be an exponential decay function with respect to the melt pool temperature. Though the resulting melt–pool height dynamics was related to the temperature and thus related to the laser power, the exponential function for the powder catchment efficiency was not validated by the experimental data. Experimental estimation of the powder catchment efficiency often had to rely on the measurement of melt–pool geometry or surface temperature and then reverse calculation of the amount of powder being absorbed by the melt pool [26]. Our prior work developed a melt–pool geometry model by approximating and parameterizing the entire material transfer (powder feeding) rate as a function of the process operating parameters [24]. Such derivation followed the modeling approach by Eagar and Tsai on dimensionless analysis of the size and shape of arc weld pools with respect to welding process variables and material parameters, considering a Gaussian heat distribution [27]. The resulting dynamic model of melt–pool geometry and temperature were validated using the experimental data or FE analysis for deposition of single-bead Ti–6AL–4V (or Inconel®718) tracks on a Laser Engineering Net Shaping (LENS) AM process.

In this paper, we extend our prior process model for melt–pool geometry in a single-bead deposition to a model for melt–pool geometry in a multibead multilayer deposition and then use this extended model to predict part height of a 3D build. Specifically, the extended model incorporates temperature history during the build process. Such temperature history is approximated by super-positioning the temperature fields generated from Rosenthal's solution of point heat sources, with one point heat source corresponding to one bead built before. Prediction on part height using the proposed model is then validated using measurements from builds with a variety of shapes, including single-bead thin wall structures, a patch build, and L-shaped structures. Preliminary results on the single-bead thin wall structures were presented at the 2017 American Control Conference [28]. In this paper, we have added the model prediction of part height for a patch build and L-shaped structures.

An Extended Model for Melt–Pool Height Dynamics in a Multibead Multilayer Deposition

In Sec. 2.1, we briefly review our prior lumped-parameter model for melt–pool geometry in a single-bead deposition on a LENS system [24]. Then, in Sec. 2.2, such model is extended to describe melt–pool height dynamics in a multibead multilayer deposition by accounting for temperature history during the build process.

Our Prior Lumped-Parameter Model for Melt–Pool Geometry.

The model was derived from the mass and energy conservations of the molten puddle and it differed from the lumped-parameter model by Doumanidis and Kwak [20] in the parameterization of the material transfer rate as a function of the process operating parameters [24]. First, mass conservation of the molten puddle requires that its mass change rate equals the material transfer rate by deposition minus material loss rate due to solidification, i.e., 
d(ρV(t))dt=ρA(t)v(t)+μf(t)
(1)

where ρ denotes the melt density, V(t) denotes the volume of the molten-puddle, A(t) denotes the maximum cross-sectional area of the molten puddle, v(t) denotes the laser scanning speed, μ denotes the powder catchment efficiency, and f(t) denotes the powder feeding rate. Assume that the molten puddle takes an ellipsoidal geometry, then V(t)=(π/6)w(t)h(t)l(t) and A(t)=(π/4)w(t)h(t), where h(t), w(t), and l(t) denote the puddle height, width, and length, respectively.

Note that the powder catchment efficiency often varies substantially, depending on the melt–pool temperature. If the laser power is not high enough to melt the powder, the catchment efficiency equals zero regardless of the powder feeding rate. For a given powder feeding rate, the higher the melt–pool temperature, the more powder the melt pool will catch and thus the higher the catchment efficiency becomes. In Ref. [24], for a given powder feeding rate, f(t) (but with varying laser power and scan speed), the material transfer rate μf(t) is parameterized by μf(t)β((η(QQc))/(πcl(TmTinitial))) if Q>Qc, and μf(t)=0 otherwise; where Q denotes the laser power, Qc denotes a critical value of the laser power, η denotes the laser transfer efficiency, β is a linear coefficient, cl is the molten material specific heat, Tm is the melting temperature of the material, and Tinitial denotes the initial temperature of the workpiece. Such parameterization was derived by extending the parameterization of dimensionless steady-state cross-sectional area of an arc weld pool with respect to a dimensionless operating parameter that captures the laser scan speed, the thermal diffusivity, and the net heat input to the workpiece [24].

Assuming that the melt pool's width equals to its length,2 i.e., w=l, and assuming a fixed ratio of width over height,3 i.e., w=rh with a constant r, the melt–pool height dynamics is then derived as follows: 
π2ρr2h2(t)dh(t)dt=π4ρrh2(t)v(t)+βη(Q(t)Qc)πcl(TmTinitial)
(2)

The parameters β and Qc can be calibrated using experimental data. In Ref. [24], for deposition of single tracks of Ti–6AL–4V on a LENS process, β=0.3026 and Qc=111.72W were first identified using measurements from a deposition with varying laser power but fixed scan speed and powder feeding rate. Then, such model parameter values were validated through prediction of melt–pool geometry in a deposition with varying scan speed but fixed laser power, where the powder feeding rate remains the same as in the deposition used for calibration.

The dynamic equation for average melt–pool temperature can be derived from the energy balance of the molten puddle (see Appendix  A). In simulating the melt–pool temperature dynamics in Ref. [24], the convection coefficient from the melt pool to solid material (substrate) was computed using FE analysis.

An Extended Model for Melt–Pool Geometry by Recomputing Initial Temperature Layer by Layer.

Note that in Eq. (2), the temperature Tinitial corresponds to the temperature of the substrate in the deposition of a single-bead track on the substrate. When more than one bead is built, the mass and energy conservations of the molten puddle still hold, but the temperature of the area before deposition is affected by the beads built before. For example, when the first bead is built on a substrate with room temperature, Tinitial is set to be the ambient temperature4 to simulate the melt–pool height in Eq. (2). During the deposition of the second bead (it can be deposited on top of the first bead or next to the first bead but in the same layer), the temperature at the predeposition point, Tinitial, is affected by the temperature distribution induced from deposition of the first bead.

In the remainder of this section, we describe how to compute the initial temperature Tinitial to account for the thermal history in the build process. Once suchTinitialis computed and updated bead by bead (or layer by layer), Eq.(2)serves as an extended bead-by-bead (or layer-by-layer) model for the height dynamics of a molten puddle. The total build height is then calculated by summing up the melt–pool height of each layer among all layers.

In calculating Tinitial in this paper, first we apply the Rosenthal's solution [30] to compute the temperature distribution induced from each bead built before. Specifically, a virtual point heat source is assigned to each past bead. This virtual heat source has the same power and velocity as the original laser heat source used to build the past bead, and it continues to travel in the same speed and direction even after the laser heat source is moved to build the next bead. Details will be given in Algorithm 1 later in this section. Second, we assume that such temperature distributions induced from past beads, which are calculated using the Rosenthal's solution under each virtual heat source, are additive, i.e., Tinitial at a point of interest can be computed by adding up the temperature contribution from each virtual source. Below, we summarize the two assumptions used in computing the initial temperature Tinitial before deposition:

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

Assumption 2. The superposition principle holds for the temperature distribution induced from the steady-state point heat source associated with deposition of each past bead (the temperature contributions of the virtual heat sources are additive).

Note that the Rosenthal's solution was derived using the linear heat equation assuming temperature-independent material properties. For the linear heat equation, the superposition principle holds for the temperature solution, hence, the temperature contributions of virtual heat sources are additive.

First, we introduce the Rosenthal's solution [30]. Let (x,y,z) denote the fixed coordinate in Fig. 1. Rosenthal's approach assumes that a laser heat source is a point source with all energy deposited into a single point and it 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, the Rosenthal's solution for the temperature distribution produced by a steady-state point heat source traveling on the surface of a semi-infinite plate is given as follows [30]:

 
T(x,y,z)=T0+q2πkRev(w+R)2a
(3)

where T(x,y,z) denotes the temperature of a point of interest at the 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 w denote the x coordinate of the point of interest in the moving coordinate (x1,y1,z1), then w=xx1. The parameter R denotes the distance from the point of interest to the laser heat source. For y1=0 and z1=0, R=(w2+y2+z2)1/2.

Next, we describe the computation of Tinitial by applying the Rosenthal's solution (Eq. (3)), using a single-bead thin-wall structure as an illustrative example. The fixed coordinate system is defined as shown in Fig. 2, with Y axis omitted since the thin-wall can be considered as 2D 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. The length of the wall is denoted by L. There could be an interlayer dwell time at the end of each layer before building the next layer. Here, for the simplicity of explanation, we assume the same dwell time between building consecutive layers, and let Δt denote such interlayer dwell time.

The initial temperature before deposition, Tinitial, is then computed as follows:

  1. (1)
    j = 1, where the first layer with track length L is built with a laser power q1 and a scan speed v1 on the substrate with room temperature T0 (see Fig. 3(a)). Let (x1,z1=0) denote the moving coordinate associated with the traveling point heat source of laser. Then, the initial temperature Tinitial at (x1,z1) before deposition is set as 
    Tinitial(x1,z1)=T0
    (4)
    If a dwell time of Δt is applied before the start of the second layer, during which the laser is turned off. The temperature of a point of interest during the dwell time is computed using Eq. (3) by defining a virtual heat source 1 (to replace the original laser heat source), which has power q1 and continues to travel with the speed v1 along the x direction5 (see Fig. 3(b)). For example, the temperature of a point of interest at the end of the dwell time can be computed by setting x1=L+v1Δt in Eq. (3).
  2. (2)
    j = 2, where the second layer is built on top of the first layer, and the laser position is moved in the z-direction to z=h¯1. The laser travels along the −x direction, with the laser power q2 and scan speed v2 (see Fig. 3(c)). Let (x2,z2=h¯1) denote the moving coordinate associated with the traveling laser heat source. The initial temperature Tinitial at (x2,z2=h¯1) before deposition is computed using the temperature distribution induced from the virtual heat source 1. As illustrated in Fig. 3(c), when the laser has traveled for time t2 on the second layer, the virtual heat source 1 has traveled v1(Δt+t2) in the x-direction from the end of the track of layer 1. As a result, by applying Eq. (3), the initial temperature Tinitial at (x2,z2=h¯1) before deposition is given as follows: 
    Tinitial(x2,z2)=T0+q12πkR12ev1(w12+R12)2a
    (5)
     
    w12=v2t2v1(Δt+t2),R12=w122+h¯12
    (6)
    where w12 denotes the x -coordinate of the point (x2,z2) in the moving coordinate (x1,z1), and R12 denotes the distance from (x2,z2) to (x1,z1).
  3. (3)
    j = 3, where the third layer is built on top of the second layer, and the laser position is moved in the z-direction to z=h¯1+h¯2. The laser travels along the x direction, with laser power q3 and scan speed v3 (see Fig. 3(d)). Let (x3,z3=h¯1+h¯2) denote the moving coordinate associated with the traveling laser heat source. Define a virtual heat source 2 to replace the original laser heat source associated with layer 2 and let it continue to travel along the −x direction with speed v2 and power q2. Then, the initial temperature Tinitial at (x3,z3=h¯1+h¯2) before deposition is now affected by both virtual heat source 1 (associated with layer 1) and virtual heat source 2 (associated with layer 2). As illustrated in Fig. 3(d), when the laser has traveled for time t3 on the third layer, Tinitial at (x3,z3) can be computed as the summation of T0, temperature rise contributed from virtual heat source 1, and temperature rise contributed from virtual heat source 2 
    Tinitial(x3,z3)=T0+q12πkR13ev1(w13+R13)2a+q22πkR23ev2(w23+R23)2a
    (7)
     
    w13=L+v3t3v1(2Δt+t2+t3),R13=w132+(h¯1+h¯2)2
    (8)
     
    w23=v3t3v2(Δt+t3),R23=w232+h¯22
    (9)
    where w13 and w23 denote the x -coordinate of the point (x3,z3) in the moving coordinate (x1,z1) and coordinate (x2,z2), respectively; R13 and R23 denote the distance from (x3,z3) to (x1,z1) and the distance from (x3,z3) to (x2,z2), respectively.

In general, let (xj,zj) denote the moving coordinate associated with the traveling laser heat source in building the jth layer of the single-bead thin-wall structure with laser power qj and scan speed vj. The initial temperature Tinitial at (xj,zj) before deposition can be calculated in a similar fashion by super-positioning of temperature contributions of virtual heat sources using the Rosenthal's solution, with one virtual heat source corresponding to one bead built before. A pseudo-code of the algorithm for computing Tinitial(xj,zj) iteratively for a single-bead thin-wall structure is given in Fig. 4.

Results

In this section, the proposed model on melt–pool height dynamics (Eq. (2)) together with the computation of initial temperature Tinitial will be used to predict the melt–pool height layer-by-layer, and then to predict part height for several geometries, including single-bead walls (see Sec. 3.1), a patch build (see Sec. 3.2), and L-shaped builds (see Sec. 3.3). Algorithm 1 in Fig. 4 for the single-bead wall will be extended to the patch build in Sec. 3.2 and to the L-shaped build in Sec. 3.3, and the implementation of such algorithms will also depend on the build plan (hatch pattern) of each build. All these parts were built of Ti–6Al–4V using an Optomec® LENS MR-7 system with a 500 W IPG photonics fiber laser. Material property parameters used in the model simulations are given in Appendix  B. The model prediction for part height of each build is then compared to its measured part height as well as to the scanned image of the build. Image scanning in this study was performed through a FARO Edge (2.7 m, seven axis) coupled with a FARO Laser Line Probe ES, which has ±0.041  mm volumetric accuracy [31]. Since in general the scanning directions by which the parts go through the FARO Edge may not align with the deposition directions used to build the parts, the raw scanning data from the FARO have to be processed with a rotation matrix. Then, the resulting image data (after rotation and shifting) are plotted using the matlab command “scatter” with the color of the “point” marker defined as a function of height, to be compared to the model prediction.

Single-Bead Thin Wall.

In Ref. [10], two single-bead walls of Ti–6AL–4V with 62 layers were deposited using an Optomec® LENS MR-7 system with a 500 W IPG photonics fiber laser. The images of the walls are shown in Fig. 5. The deposition direction switched every single layer in building the walls. There was no interlayer dwell time in building wall 1 (Fig. 5(a)), and there was a 20 s interlayer dwell time in building wall 2 (Fig. 5(b)). The process parameters used in the experiments and the resulting wall geometric measurements are given in Table 1.

We then simulate the melt–pool height dynamics using Eq. (2), where the initial temperature distribution at each layer is computed from Algorithm 1 in Fig. 4. In addition, the model parameter values for β and Qc used in the simulation are set to be β=0.3026 and Qc=111.72W, as explained in Sec. 2.1. An estimated melt–pool width over height ratio is computed using the wall width and average layer height (wall height divided by 62 layers) of wall 2 listed in Table 1, i.e., r=(w/h)=(2.2/(10.7/62))=12.75. This ratio for r is used in simulation for both walls. The simulated initial temperature distribution Tinitial with respect to wall distance and layer number is plotted in Fig. 6 for both walls. Tinitial is set to the ambient temperature for the first layer. The deposition direction for the first layer is from left to right, and then from right to left for the second layer, hence Tinitial at the second layer starts from a higher value at the right end of the wall and then gradually decreases when the deposition moves to the left. Similarly, the deposition direction for the third layer is from left to right, thus Tinitial at the third layer starts at a higher value at the left end of the wall and then decreases when the deposition moves to the right. As expected, the initial temperature of wall 2 at a given point of a layer is lower than that of wall 1 due to the 20 s interlayer dwell time.

Figure 7 shows the model predicted wall height after 62 layers along wall distance versus the measured wall height at the middle of the wall, which was measured using a caliper and used as a reference value for comparison. The length of wall 1 is longer than that of wall 2. For wall 1, two humps are predicted at the two ends of the wall. Specifically, the predicted build height along wall distance less than 3 mm or larger than 37 mm is lower than the height at the middle section of the wall due to the transient response of the melt–pool height dynamics; the over-builds along wall distance 3–10 mm and 30–37 mm are due to a higher initial temperature Tinitial at the start of every layer than the rest of the layer (from the build path). The prediction of humps occurring at both ends of wall 1 is in agreement with the image of wall 1 (see Fig. 5(a)). When 20 s dwell time is applied, the predicted top of wall 2 is relatively flat, and no hump is predicted at the two ends of wall 2. Along the wall distance 6.5–33 mm for wall 1, the predicted wall height is 0.8–2.0% higher than the measured midwall height; for wall 2, the predicted wall height is ∼3% lower than the measured midwall height along the wall distance 10–27 mm. Figure 7 also shows the model prediction for wall height after every ten layers of deposition. As expected, the build height of wall 2 at a given layer is lower than that of wall 1 at the same layer due to the 20 s dwell time.

Figure 8 shows comparison of the predicted height contour of wall 2 versus its scanned image. Wall 1 was cut for other analyses and thus not available for comparison. Figure 8 shows that the predicted height contour (represented by the solid line) agrees well with the scanned image (represented by a scatter plot of the image data) except missing prediction of overbuilds at the two ends of wall 2. Miss-prediction of overbuilds at the two ends of wall 2 is possibly attributed by the way by which the temperature distribution is computed during the cool-down period after finishing each bead, including the dwell time. Future work will be conducted to improve the thermal modeling during cool down.

Sensitivity Analysis ofTinitialwith respect to Laser Power and Scan Speed: Note that the initial temperature Tinitial is affected by both laser power and scan speed, as shown by Eq. (3). For wall 1, under a fixed scan speed of 8.5 mm/s, the initial temperature distribution Tinitial(xj,zj) is recalculated for each laser power from the set Q=[369W,389.5W,410W,430.5W,451W] (corresponding to 90–110% of the nominal laser power 410 W) to show how varying laser power would affect Tinitial(xj,zj). Then, the analysis is repeated for a fixed laser power of 410 W but with the laser scan speed from a set v=[7.65,8.075,8.5,8.925,9.35] (mm/s) (corresponding to 90–110% of the nominal laser scan speed 8.5 mm/s). Figure 9(a) shows Tinitial of a few layers at x = half of the wall length, as a function of the laser power, and Fig. 9(b) shows the resulting wall height profile under different laser power. Figure 10(a) shows Tinitial of a few layers at x = half of the wall length, as a function of the laser scan speed, and Fig. 10(b) shows the resulting wall height profile under different laser scan speed. The experimental results for building the single-bead wall under different laser power and scan speed are not available to validate the sensitivity analysis here. In the future work, FE analysis will be conducted to compute the initial temperature distribution as a function of different laser power and scan speed, and the results will be compared to the model prediction and sensitivity analysis here.

Patch Build.

A schematic plot of a patch deposit of Ti–6AL–4V using an Optomec® LENS MR-7 [32] is shown in Fig. 11. The square patch has a width of about 25.4 mm (1 in), deposited on a substrate with 76.2 mm in length, 50.8 mm in width, and 2.54 mm in depth (see Figs. 11(a) and 11(b)). A parallel hatch pattern (build plan) was used where the deposition direction switched every single bead as shown in Fig. 11(c). A total of 36 beads were deposited in each layer, and three layers were built. The measured patch height is around 1.64 mm. Table 2 lists the process parameters.

Figure 12 shows the patch image. In Fig. 12(a), A, B, C, D denote four cross-sectional line segments along the y-direction, where x = 18 mm for A, x = 9 mm for B, x = 4.5 mm for C, and x = 0.4 mm for D. Two quarter patches were cut from the half patch in Fig. 12(a) for other analyses before being placed next to each other for scanning (see Fig. 12(b)). Corresponding to line segments A, B, C, and D in Fig. 12(a), four segments A1/A2, B1/B2, C1/C2, and D1/D2 are identified from the scan image of the two quarter patches. To ensure enough image data to be extracted, rather than a line segment, each of A1/A2, B1/B2, and C1/C2 is defined to be a thin-band segment with a width of 0.8 mm, and D1/D2 is a thin band with a width of 0.2 mm. Then, A1, B1, C1, D1 are stitched to A2, B2, C2, D2, respectively, to form the scan images of segments A, B, C, D used for comparison to model prediction, where the material loss in cutting the two quarter patches is considered negligible to allow the image stitching.

Initial Temperature Calculation for Patch Build: The algorithm for computing initial temperature Tinitial for a 2D structure such as a single-bead wall in Sec. 2.2 can be easily extended to compute the initial temperature for the 3D patch with a parallel build plan, where the deposition direction aligns with one side of the patch. Consider a fixed coordinate (x,y,z) (in right-handed coordinates) with the x direction aligns with the deposition direction of the first bead and the origin at the corner where the first bead starts to deposit. Similar to the algorithm used for the single-bead thin wall, a virtual heat source is defined to be associated with each bead built before (either below or next in the same layer). Then, the initial temperature Tinitial(xj,yj,zj) at the deposition point (xj,yj,zj) can be computed by recalculating w and R in Eq. (3) in terms of each virtual moving heat source.

Specifically, for the jth bead,6 the initial temperature Tinitial(xj,yj,zj) at the deposition point (xj,yj,zj) is computed as 
Tinitial(xj,yj,zj)={T0,j=1T0+i=1j1qi2πkRijevi(wij+Rij)2a,j2
(10)

where wij=xjxi denotes the x coordinate of the deposition point (xj,yj,zj) in the moving coordinate (xi,yi,zi) (associated with the virtual heat source i defined for bead i, i < j), and Rij denotes the distance from (xj,yj,zj) to (xi,yi,zi).

Simulation Results: In simulating the height dynamics, the model parameters β and Qc were set to be the same as in Sec. 3.1, and the melt–pool width over height ratio r was recalibrated using an estimated melt–pool width and height which led to r = 3.57. A 3D surface plot of the simulated three-layer patch build is displayed in Fig. 13, where the top surface after each layer being deposited is plotted. Figure 13 shows that the four corners of the patch taper off a bit due to the transient response of the melt–pool height dynamics. It also shows that the top surface of the patch is not flat, with slight overbuilds on ends.

In Fig. 14, we compare the model predicted height of each line segment A, B, C, and D defined in Fig. 12(a) to the scanned image of the corresponding thin-band segment A1/A2, B1/B2, C1/C2, and D1/D2 defined in Fig. 12(b), respectively. Figure 14 shows that the average height predicted by the model is very close to the average height of the actual build (less than 15%), but the model always over predicts the height at the two edges of each segment.

L-Shaped Builds.

In this section, two L-shaped structures of Ti–6AL–4V are considered, which were built, with different interlayer dwell time, on substrates of 76.2mm×76.2mm×6.35mm using an Optomec® LENS MR-7 system with a 500 W IPG Yb-doped fiber laser. As shown in Fig. 15, each L-shaped structure consists of a one-bead-wide leg and a three-bead-wide leg, with a target length of 25.4 mm (1 in) and a target height of 50.8 mm (2 in) for each leg. The scanning images in Figs. 15(c)15(d) are rotated, isometric view of Figs. 15(a) and 15(b). A laser power of 450 W was applied to build the one-bead leg and, in order to achieve roughly equal height for the two legs, 350 W was applied to build the three-bead leg. Table 3 summarizes the process parameters used in the experiment. The experimental environment setup was the same as Ref. [29], and readers are referred to Ref. [29] for more details. A total of 286 layers were deposited, with the eight-hatch build plan shown in Fig. 16, where hatches 1–4 and 5–8 denote the deposition path in odd and even layers, respectively. The first L-shaped structure (L-1) was built with no interlayer dwell time, and the second structure (L-2) was built with 4 s interlayer dwell time. It should be noted that in Fig. 16, the x–y coordinates are defined in terms of the deposition directions of the legs, and they are not aligned with the scanning x–y directions shown in Fig. 15. The latter indexed the scanning data obtained from the FARO Edge. A rotation matrix was then calculated to map the scanning data in the scanning x–y directions into the x–y coordinates defined in Fig. 16.

Initial Temperature Calculation for L-shaped Structure: At a point of interest (xj,yj,zj), the initial temperature Tinitial(xj,yj,zj) can be decomposed into two parts: (i) the temperature induced by past beads deposited along the x direction (any past hatches 2–7), and it is denoted by TinitialX(xj,yj,zj) that can be computed by adding up the temperature contributions of virtual heat sources defined for past hatches 2–7; and (ii) the temperature induced by past beads deposited along the y direction (any past hatches 1 and 8), denoted by TinitialY(xj,yj,zj) that can be computed by adding up the temperature contributions of virtual heat sources defined for past hatches 1 and 9. Then, due to the assumption of superposition on temperature contributions of virtual heat sources 
Tinitial(xj,yj,zj)=T0+TinitialX(xj,yj,zj)+TinitialY(xj,yj,zj)
(11)
where TinitialX(xj,yj,zj) can be calculated similarly as Eq. (10) 
TinitialX(xj,yj,zj)=iHXqi2πkRijev(wijx+Rij)2a
(12)

HX denotes the set of any past hatches 2–7 before the current hatch, wijx=xjxi denotes the x coordinate of the deposition point (xj,yj,zj) in the moving coordinate (xi,yi,zi) (associated with the virtual heat source i), and Rij denotes the distance from (xj,yj,zj) to (xi,yi,zi).

Note that the virtual heat source associated with each past hatch 1 moves in the y direction, and the virtual heat source associated with each past hatch 8 moves in the −y direction, the initial temperature affected by the set of past hatches 1 and 8 can be computed as 
TinitialY(xj,yj,zj)=iHYqi2πkRijev(wijy+Rij)2a
(13)

where HY denotes the set of any past hatches 1 and 8 before the current hatch, wijy=yjyi denotes the y coordinate of the deposition point (xj,yj,zj) in the moving coordinate (xi,yi,zi) (associated with the virtual heat source i defined for each past hatches 1 and 8), and Rij denotes the distance from (xj,yj,zj) to (xi,yi,zi).

Simulation Results: In simulating the melt–pool height dynamics, the same model parameters, β=0.3026, Qc=111.72W, and r=(w/h)=12.75, as used for simulating the thin walls in Sec. 3.1 are applied here. Figure 17 shows the 3D plot of the simulated L-shaped builds. Similar to the thin walls and the patch build, Fig. 17 shows that the outside edges of the two legs of the L-shaped structures taper off a bit as a result of the transient response. Further, there is slight overbuild at the intersection of the two legs due to the repeated laser heating according to the build path.

Figure 18 shows comparison of the model predicted height for each leg versus its scan image for both L-1 and L-2. For the three-bead leg, the vertical component built by the set of hatches 2 and 7 is referred to as “wall 2,” the vertical component built by the set of hatches 3 and 6 is referred to as “wall 3,” and the vertical component built by the set of hatches 4 and 7 is referred to as “wall 4.” It should be noted that in Figs. 18(b) and 18(d), walls 2, 3 and 4 cannot be differentiated in the scan image, but the model predicted height for each individual walls 2, 3, and 4 is plotted. For L-1 (no dwell time), the average predicted height for the one-bead leg is 52.9 mm (2% higher than the measured average height of 52.17 mm), and the average predicted height for the three-bead leg (considering pointwise maximum of the predicted height of walls 2, 3, and 4) is 48.9 mm (6% lower than the measured average height of 51.83 mm). For L-2 (4 s dwell time), the average predicted height for the one-bead leg is 50 mm (3% lower than the measured average height of 51.3 mm), and the average predicted height for the three-bead leg (considering pointwise maximum of the predicted height of walls 2, 3, and 4) is 45.8 mm (7% lower than the measured average height of 49.1 mm).

Discussion

This paper extends our prior model for melt–pool height dynamics in a single-bead deposition to a model for melt–pool geometry in a multibead multilayer deposition, by taking into account the thermal effect on the initial temperature Tinitial due to the deposition of past beads. The calculation of the temperature distribution Tinitial before deposition was approximated via applying super-position to Rosenthal's solution on moving point heat source associated with each past-built bead. Hence, such approximation inherits all assumptions and limitations associated with the classic Rosenthal's solution. For example, the model assumes the point laser source which allows to apply the closed-form Rosenthal's solution but ignores the effect of laser diameter on melt pool geometry (one study on how laser diameter affected melt pool geometry given fixed laser scan speed and power can be found in Ref. [33]). Furthermore, Rosenthal's solution was derived from the linear heat equation assuming temperature-independent material properties. However, material properties, especially thermal conductivity, vary substantially when temperature changes. This is a limitation of the model.

It should also be noted that the computation of Tinitial in this paper did not take into account many other factors, e.g., the thickness of the substrate would affect Tinitial; in addition, when the laser moves slowly, a small area in front of the laser will be heated up (in addition to the thermal effect due to deposition of past beads). Last, the proof of concept study in this paper was limited to several simple builds such as single-bead walls, patch builds, and L-shaped structures. Future work will extend the proposed approach to more realistic complex shapes. We believe that it is possible to design a suite of virtual point heat sources to account for the thermal history during the build process for a certain class of shapes. In the future work, a FE analysis will be conducted to compute the initial temperature distribution before deposition, as well as perform sensitivity analysis of Tinitial with respect to laser power and laser scan speed, for builds with more complex shapes. The FE results will be used as comparison basis to guide the improvement of the model proposed in this paper.

The simulations in this paper were run on a Lenovo G780 Laptop with Intel® Core i5-3210M processor 2.50 GHz. The simulation for each of the single-bead walls took 1.5 s, the simulation for the patch build took 15 s, and the simulation for each of the L-shaped structures took 229 s.

Though not discussed in this paper, we expect that the average melt–pool temperature of each bead in a 3D part can be simulated by extending the temperature dynamics in Appendix  A to account for thermal history due to beads built before. Such modeling and simulation will be investigated in our future research as well.

Conclusion

In this paper, our prior developed lumped-parameter model for melt–pool geometry of single-bead deposition is extended to a model for melt–pool height dynamics in multibead multilayer depositions by incorporating the temperature history during the build process. Such temperature history is computed using the temperature field generated from super-positioning of Rosenthal's solution of point heat sources, with one heat source corresponding to one bead built before. Such model for melt–pool height dynamics is then used to predict part height for several multilayer builds, including single-bead thin walls, a patch build, and L-shaped structures, which were all built of Ti–6AL–4V using an Optomec® LENS MR-7 system. The model predictions for part height match the experimental measurements very well with low error rates for all three types of structures, and in addition, the model could predict the edge tapering-off due to transient response as well as overbuild at corners most of the time. We acknowledge that the proof-of-concept work in this paper is based on many assumptions and approximations, and a comparison to the temperature history computed by a finite element analysis is deemed for the future work. Nevertheless, developing such a lumped-parameter model for melt–pool geometry in a multibead multilayer deposition not only allows quick prediction of part height without resorting to lengthy FE computation, but also allows the design of model-based control systems to improve geometric accuracy of the final build in a directed energy deposition process.

Funding Data

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

Appendix A: Dynamic Equation for Melt–Pool
Average Temperature

Doumanidis and Kwak derived the lumped-parameter model governing the dynamics of melt–pool average temperature in terms of energy conservation in a molten puddle [20], where the internal energy rate of accumulation in the puddle is equal to the difference of power inflow of material addition and power outflow loss to solidified material, plus total external heat transfer rate from the puddle surface. By assuming no preheating of the powder (i.e., the power inflow of material addition equals zero), the energy balance of the puddle leads to 
d(ρV(t)e)dt=ρA(t)v(t)eb+Ps
(A1)
where e denotes the specific internal energy of the melt pool 
e=cs(TmT0)+hSL+cl(T(t)Tm)
(A2)
where cs denotes the solid material specific heat, T0 denotes the ambient temperature, hSL denotes the specific latent heat of solidification, and T denotes the average melt–pool temperature. The specific energy of the solidified material, denoted by eb in Eq. (A1) satisfies 
eb=cs(TmT0)
(A3)
The total thermal transfer at the puddle surface, denoted by Ps in Eq. (A1), can be computed as follows: 
Ps=ηQ(t)Asαs(T(t)Tm)AG×[αG(T(t)T0)+εσ(T4(t)T04)]
(A4)

where η denotes the laser absorption efficiency for the laser power Q, αs denotes the convection coefficient, αG denotes the heat transfer coefficient, ε denotes the surface emissivity, and σ is the Stefan–Boltzmann constant. The parameters As and AG denote the areas of melt–pool interface to the substrate and to the free surface, respectively. Assuming ellipsoidal geometry of the melt–pool, As=(π/4)w(t)l(t) and AG=(π/23)[w(t)h(t)l(t)]2/3.

When such dynamic equation of melt–pool temperature was applied in our prior work [24], the convection coefficient to solid material (substrate) αs was estimated using finite element analysis, and the heat transfer coefficient αG was calibrated using experimental measurements.

Appendix B: Material Properties of Ti–6AL–4V
Used in Model Prediction

2

In Ref. [29] (Table 4), measurements from coaxial ThermaViz® images of melt-pool for L-shaped builds showed that the average values of melt-pool length-to-width ratio at different build locations varied within [0.873, 1.085].

3

This is an assumption used to simplify the model derivation. It should be noted that such ratio in general varies, and we will work on dropping this assumption in the future.

4

Assuming that the area in front of the laser has not been heated up by the laser yet; we do acknowledge that if the laser moves slowly, the immediate area in front of the laser should be heated up by the laser.

5

This is an approximation, and we acknowledge that it does not truly reflect the temperature distribution during the cool-down period.

6

This index corresponds to the ordered index of the entire set of beads across all layers, e.g., first bead of layer 2 in Fig. 11(c) is counted as the 37th bead in Eq. (10).

7

Note that the build plan (x, y, z) here are left-handed coordinates.

References

References
1.
Heralić
,
A.
,
Christiansson
,
A.-K.
, and
Lennartson
,
B.
,
2012
, “
Height Control of Laser Metal-Wire Deposition Based on Iterative Learning Control and 3D Scanning
,”
Opt. Lasers Eng.
,
50
(
9
), pp.
1230
1241
.
2.
Mazumder
,
J.
,
Koch
,
J.
,
Nagarathnam
,
K.
, and
Choi
,
J.
,
1996
, “
Rapid Manufacturing by Laser Aided Direct Deposition of Metals
,”
Adv. Powder Metall. Part. Mater.
,
4
(
1
), pp.
15
107
.
3.
Griffith
,
M. L.
,
Schlienger
,
M. E.
,
Harwell
,
L. D.
,
Oliver
,
M. S.
,
Baldwin
,
M. D.
,
Ensz
,
M. T.
,
Essien
,
M.
,
Brooks
,
J.
,
Robino
,
C. V.
, and
Smugeresky
,
J. E.
,
1999
, “
Understanding Thermal Behavior in the Lens Process
,”
Mater. Des.
,
20
(
2
), pp.
107
113
.
4.
Beuth
,
J.
, and
Klingbeil
,
N.
,
2001
, “
The Role of Process Variables in Laser-Based Direct Metal Solid Freeform Fabrication
,”
JOM
,
53
(
9
), pp.
36
39
.
5.
Das
,
S.
,
Wohlert
,
M.
,
Beaman
,
J. J.
, and
Bourell
,
D. L.
,
1998
, “
Producing Metal Parts With Selective Laser Sintering/Hot Isostatic Pressing
,”
JOM
,
50
(
12
), pp.
17
20
.
6.
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
.
7.
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
.
8.
Marimuthu
,
S.
,
Clark
,
D.
,
Allen
,
J.
,
Kamara
,
A. M.
,
Mativenga
,
P.
,
Li
,
L.
, and
Scudamore
,
R.
,
2013
, “
Finite Element Modelling 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
.
9.
Michaleris
,
P.
,
2014
, “
Modeling Metal Deposition in Heat Transfer Analyses of Additive Manufacturing Processes
,”
Finite Elem. Anal. Des.
,
86
(
1
), pp.
51
60
.
10.
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
(
1
), pp.
9
19
.
11.
Toyserkani
,
E.
, and
Khajepour
,
A.
,
2006
, “
A Mechatronics Approach to Laser Powder Deposition Process
,”
Mechatronics
,
16
(
10
), pp.
631
641
.
12.
Kruth
,
J.-P.
,
Mercelis
,
P.
,
Van Vaerenbergh
,
J.
, and
Craeghs
,
T.
,
2007
, “
Feedback Control of Selective Laser Melting
,”
Third International Conference on Advanced Research in Virtual and Rapid Prototyping
, Leiria, Portugal, Sept. 24–29, pp.
521
527
.https://lirias.kuleuven.be/bitstream/123456789/185342/1/krufcs.pdf
13.
Tang
,
L.
, and
Landers
,
R. G.
,
2010
, “
Melt Pool Temperature Control for Laser Metal Deposition Processes—Part I: Online Temperature Control
,”
ASME J. Manuf. Sci. Eng.
,
132
(
1
), p.
011010
.
14.
Song
,
L.
,
Bagavath-Singh
,
V.
,
Dutta
,
B.
, and
Mazumder
,
J.
,
2012
, “
Control of Melt Pool Temperature and Deposition Height During Direct Metal Deposition Process
,”
Int. J. Adv. Manuf. Technol.
,
58
(
1–4
), pp.
247
256
.
15.
Song
,
L.
, and
Mazumder
,
J.
,
2011
, “
Feedback Control of Melt Pool Temperature During Laser Cladding Process
,”
IEEE Trans. Control Syst. Technol.
,
19
(
6
), pp.
1349
1356
.
16.
Fathi
,
A.
,
Khajepour
,
A.
,
Toyserkani
,
E.
, and
Durali
,
M.
,
2007
, “
Clad Height Control in Laser Solid Freeform Fabrication Using a Feedforward PID Controller
,”
Int. J. Adv. Manuf. Technol.
,
35
(
3–4
), pp.
280
292
.
17.
Fathi
,
A.
,
Khajepour
,
A.
,
Durali
,
M.
, and
Toyserkani
,
E.
,
2008
, “
Geometry Control of the Deposited Layer in a Nonplanar Laser Cladding Process Using a Variable Structure Controller
,”
ASME J. Manuf. Sci. Eng.
,
130
(
3
), p.
031003
.
18.
Urbanic
,
R. J.
,
Saqib
,
S. M.
, and
Aggarwal
,
K.
,
2016
, “
Using Predictive Modeling and Classification Methods for Single and Overlapping Bead Laser Cladding to Understand Bead Geometry to Process Parameter Relationship
,”
ASME J. Manuf. Sci. Eng.
,
138
(
5
), p.
051012
.
19.
Jin
,
Y.
,
Qin
,
S.
, and
Huang
,
Q.
,
2016
, “
Offline Predictive Control of Out-of-Plane Shape Deformation for Additive Manufacturing
,”
ASME J. Manuf. Sci. Eng.
,
138
(
12
), p.
121005
.
20.
Doumanidis
,
C.
, and
Kwak
,
Y.-M.
,
2001
, “
Geometry Modeling and Control by Infrared and Laser Sensing in Thermal Manufacturing With Material Deposition
,”
ASME J. Manuf. Sci. Eng.
,
123
(
1
), pp.
45
52
.
21.
Devesse
,
W.
,
De Baere
,
D.
, and
Guillaume
,
P.
,
2014
, “
The Isotherm Migration Method in Spherical Coordinates With a Moving Heat Source
,”
Int. J. Heat Mass Transfer
,
75
, pp.
726
735
.
22.
Devesse
,
W.
,
De Baere
,
D.
, and
Guillaume
,
P.
,
2014
, “
Design of a Model-Based Controller With Temperature Feedback for Laser Cladding
,”
Phys. Procedia
,
56
, pp.
211
219
.
23.
Sammons
,
P. M.
,
Bristow
,
D. A.
, and
Landers
,
R. G.
,
2013
, “
Height Dependent Laser Metal Deposition Process Modeling
,”
ASME J. Manuf. Sci. Eng.
,
135
(
5
), p.
054501
.
24.
Wang
,
Q.
,
Li
,
J.
,
Gouge
,
M.
,
Nassar
,
A. R.
,
Michaleris
,
P.
, and
Reutzel
,
E. W.
,
2016
, “
Physics-Based Multivariable Modeling and Feedback Linearization Control of Melt-Pool Geometry and Temperature in Directed Energy Deposition
,”
ASME J. Manuf. Sci. Eng.
,
139
(
2
), p.
021013
.
25.
Cao
,
X.
, and
Ayalew
,
B.
,
2015
, “
Control-Oriented MIMO Modeling of Laser-Aided Powder Deposition Processes
,”
American Control Conference
(
ACC
), Chicago, IL, July 1–3, pp.
3637
3642
.
26.
Lin
,
J.
, and
Steen
,
W. M.
,
1998
, “
An In-Process Method for the Inverse Estimation of the Powder Catchment Efficiency During Laser Cladding
,”
Opt. Laser Technol.
,
30
(
2
), pp.
77
84
.
27.
Eagar
,
T. W.
, and
Tsai
,
N. S.
,
1983
, “
Temperature Fields Produced by Traveling Distributed Heat Sources
,”
Weld. J.
,
62
(
12
), pp.
346
355
. http://eagar.mit.edu/Publications/Eagar036.pdf
28.
Li
,
J.
,
Wang
,
Q.
,
Michaleris
,
P.
, and
Reutzel
,
E. W.
,
2017
, “
Model Prediction for Deposition Height During a Direct Metal Deposition Process
,”
American Control Conference
(
ACC
), Seattle, WA, May 24–26, pp.
2188
2194
.
29.
Kriczky
,
D. A.
,
Irwin
,
J.
,
Reutzel
,
E. W.
,
Michaleris
,
P.
,
Nassar
,
A. R.
, and
Craig
,
J.
,
2015
, “
3D Spatial Reconstruction of Thermal Characteristics in Directed Energy Deposition Through Optical Thermal Imaging
,”
J. Mater. Process. Technol.
,
221
, pp.
172
186
.
30.
Rosenthal
,
D.
,
1941
, “
Mathematical Theory of Heat Distribution During Welding and Cutting
,”
Weld. J.
,
20
(
5
), pp.
220
234
.
31.
FARO, 2015, “Laser Trackers: Defining Accuracy,” FARO, Lake Mary, FL, accessed Apr. 9, 2017, http://www.faro.com/edge
32.
Applied Research Laboratory at the Pennsylvania State University
,
2015
, “
America Makes—Reborn in the USA: Laser/Powder Fed Directed Energy Deposition, Task 2.1 Process Development
,” Pennsylvania State University, State College, PA, Technical Report No. 3-11-15.
33.
Cheng
,
B.
, and
Chou
,
K.
,
2013
, “
Melt Pool Geometry Simulations for Powder-Based Electron Beam Additive Manufacturing
,”
Annual International Solid Freeform Fabrication Symposium
(
SFF
), Austin, TX, Aug. 12–14, pp.
644
654
.https://sffsymposium.engr.utexas.edu/Manuscripts/2013/2013-51-Cheng.pdf