## Abstract

Despite the advantages and emerging applications, broader adoption of powder bed fusion (PBF) additive manufacturing is challenged by insufficient reliability and in-process variations. Finite element modeling and control-oriented modeling have been shown to be effective for predicting and engineering part qualities in PBF. This paper first builds a finite element model (FEM) of the thermal fields to look into the convoluted thermal interactions during the PBF process. Using the FEM data, we identify a novel surrogate system model from the laser power to the melt pool width. Linking a linear model with a memoryless nonlinear submodel, we develop a physics-based Hammerstein model that captures the complex spatiotemporal thermomechanical dynamics. We verify the accuracy of the Hammerstein model using the FEM and prove that the linearized model is only a representation of the Hammerstein model around the equilibrium point. Along the way, we conduct the stability and robustness analyses and formalize the Hammerstein model to facilitate the subsequent control designs.

## 1 Introduction

Different from conventional subtractive machining, additive manufacturing (AM) builds a part directly from its digital model by joining materials layer by layer. In particular, by applying high-precision lasers or electron beams as the energy source, powder bed fusion (PBF) AM has enabled unprecedented fabrication of complex parts from polymeric and metallic powder materials. However, broader adoption of the technology remains challenged by insufficient reliability and in-process variations induced by, for example, uncertain laser-material interactions, environmental vibrations, powder recycling, imperfect interactions of mechanical components, and recursive thermal histories of materials [1–5].

A typical part in PBF is built from many thousands of thin layers.Within each layer, the energy beam is regulated to follow trajectories predefined by the part geometry in a slicing process. After one layer is finished printing, a new thin layer of powder will be spread on top, and then another cycle begins. Appropriate modeling of this sophisticated dynamic system plays a fundamental role in understanding and regulating the PBF and related techniques such as laser metal deposition (LMD). Current researches employ the finite element model (FEM) to explore energy deposition mechanisms and control-oriented modeling to build mathematical models for regulating in-process variations. Particularly, Masoomi et al. [7], Hussein et al. [10], and Foroozmehr et al. [11] adopt FEM to investigate the effects of various scan configurations on the thermal fields of the powder bed, the geometries of the melt pool, and the mechanical properties of the printed parts. In control-oriented modeling, Kruth et al. [2], Song and Mazumder [12], Cao and Ayalew [13], and Sammons et al. [14] apply low-order system models, and Cao and Ayalew [13], Sammons et al. [14], and Fathi et al. [15] further build nonlinear submodels of LMD from laser power and scan speed to layer height and melt pool temperature to cover more process dynamics. Based on the attained models, subsequent control algorithms such as proportional–integral–derivative (PID) control [2,16–18], sliding mode control [15], predictive control [12], and iterative learning control [19] have proved their efficiencies in improving the dimensional accuracy of the printed parts.

This paper establishes a new modeling and understanding of PBF by taking advantage of the FEM and control-oriented modeling. We first develop an FEM of the thermal fields to look into the convoluted thermal interactions during the PBF process. The developed FEM then serves as a simulation platform to provide data for verifying and identifying parameters of the proposed modeling schemes. In the control-oriented modeling of PBF, stepping beyond commonly used low-order system models, this paper develops a physics-based Hammerstein model that accommodates more of the complex spatiotemporal thermomechanical dynamics. The Hammerstein model is formulated by concatenating a memoryless nonlinear submodel derived from the Rosenthal equation to a linear model obtained from system identification techniques with the laser power as the input and the melt pool width as the output. We verify the accuracy of the Hammerstein model using the FEM and prove that the identified model is only a linear representation of the Hammerstein model around the equilibrium point. Along the way, we analyze the stability and robustness properties of the models and present a generic control scheme of the Hammerstein model.

The remainder of this paper is structured as follows. Section 2 builds the FEM of the thermal fields in PBF. Section 3 identifies the linear plant model from the FEM. Section 4 derives the closed-form expressions of the steady-state melt pool width and furthermore develops and analyzes the main Hammerstein model. Section 5 concludes the paper.

## 2 Finite Element Model of Thermal Fields in Powder Bed Fusion

*k*is the thermal conductivity,

*c*

_{p}the specific heat capacity,

*ρ*the effective density,

*t*the time,

*T*the temperature, and

*q*

_{s}the rate of local internal energy generated per unit volume [20]. When no confusion would arise in the context,

*T*(

*x*,

*y*,

*z*,

*t*) is abbreviated to

*T*in the remainder of this paper.

### 2.1 Nonlinear Phase Change and Temperature-Dependent Thermal Properties.

*L*

_{f}by introducing the effective heat capacity [21]:

*T*

_{0}is the ambient temperature,

*T*

_{sol}the solidus temperature,

*T*

_{m}the melting point,

*c*

_{p1}the heat capacity of the solid or powder, and

*c*

_{p2}the heat capacity of the liquid.

*k*,

*c*

_{p}, and

*ρ*in Refs. [6,7] for the solid and liquid materials. We generate the thermal properties of the powder material from those of the solid material by considering the porosity

*ϕ*[8,9]:

*ϕ*is expressed as

*ϕ*

_{0}denoting the initial porosity. Here, the heat capacity is assumed to be the same for the powder and solid materials [8].

### 2.2 Initial Condition, Boundary Conditions, and Laser Beam Profile.

*T*(

*x*,

*y*,

*z*, 0) =

*T*

_{0}. When the substrate (left plot of Fig. 2) is designed to be large enough compared to the heat affected zone, one boundary condition is established by assuming the bottom (

*z*=

*h*

_{sub}) of the substrate has no heat loss in the rapid laser scanning: $\u2212k(\u2202T/\u2202z)|z=hsub=0$. The other boundary condition considers surface conduction, convection, and radiation:

*Q*is the input heat flux,

*h*

_{c}the convection heat transfer coefficient,

*ɛ*the emissivity, and

*σ*

_{B}the Stefan–Boltzmann constant. Here, we assume

*Q*has a Gaussian laser beam profile: $Q\u2248(2q/\pi R2)e\u2212(2r2/R2)$, where

*q*is the laser power,

*R*the effective laser beam radius, and

*r*the radial distance from a certain point to the center of the laser spot. The Appendix has listed the process parameters used in this study.

### 2.3 Meshing and Scanning Schemes.

The left plot of Fig. 2 shows the built FEM with a substrate and a thin layer of powder bed. For melt pools in the scale of around 248 *μ*m in diameter, we use a selective meshing scheme to balance model accuracy with computation time: a fine quad-and-swept mesh with a maximum element size of 60 *μ*m is applied to the central powder bed region that directly interacts with the energy beam, whereas less finer tetrahedral mesh (1.1 mm) and triangular-and-swept mesh (0.83 mm) are applied to the substrate and the peripheral powder bed, respectively. We adopt the FEM practice to use coarser element sizes for the substrate and the peripheral powder bed that undergoes less significant heat transfer than the central melt pool. We consider both single and multiple scans of the laser beam on the powder bed. The left plot of Fig. 2 illustrates the bidirectional scan scheme used in this study.

The developed FEM has been validated experimentally and analytically in Ref. [22] and serves as a simulation platform in this paper. Later on the data generated from the FEM such as melt pool width will be used to identify and verify the accuracies of the proposed models. In FEM, melt pool widths are generated from the temperature distribution *T*(*x*, *y*, *z*, *t*) by searching around the position of the laser beam to find the maximum width of the melt pool geometry bounded by *T*_{m}.

## 3 Linear Model Around Quasi-Steady Equilibrium

From the developed FEM, we identify the linear plant model as *P*(*s*) = 0.001671/(*s* + 1055) from laser power changes *δq* to melt pool width changes *δw* around the equilibrium point at (*q*_{0} = 60 W, *w*_{0} = 248.41 *μ*m). Here, *q* = *q*_{0} + *δq* and *w* = *w*_{0} + *δw* are the actual laser power fed to the FEM and the melt pool width generated from the FEM, respectively. The input signals fed to the FEM include a pseudorandom binary sequence signal and multiple sinusoidal signals (10–300 Hz), with a magnitude of 20 W and an add-on DC component of *q*_{0} = 60 W. As shown in Fig. 3, the frequency responses of the measured and identified systems match well with each other. From a physics viewpoint, the low-pass dynamics is attributed to the high-density energy deposition of the laser and the first-order temporal dynamics of the temperature evolution in Eq. (1). Under the sampling time *t*_{s} of 0.5 ms, the zero-order-hold equivalent of the plant model is *P*_{d}(*z*) = 6.493 × 10^{−7}/(*z* − 0.5901). For a more general and unified analysis (Fig. 4), *P*_{d}(*z*) can be further normalized to have a unit DC gain by *P*(*z*) = *P*_{d}(*z*)/*c* = 0.4099/(*z* − 0.5901), where *c* is the DC gain of *P*_{d}(*z*).

## 4 Hammerstein Model in Powder Bed Fusion

In this section, we show the limit of the linear model subject to the complicated nonlinear thermomechanical dynamics of PBF and build a new physics-based Hammerstein model to address the limitations. After that, we analyze the stability and robustness of the models. The Hammerstein model is conventionally employed in system identification for nonlinear systems, consisting of a nonlinear static element followed by a linear dynamic element. Recent studies of the Hammerstein model aim at parameter estimation and neural network-based solutions [23–25]. Here, we repurpose the method for identifying a nonlinear model for PBF.

### 4.1 Core Physics of the Melt Pool at Quasi-Steady-State.

*ξ*,

*y*,

*z*) is a coordinate system attached to the moving source,

*u*

_{x}is the laser scanning speed, $r=\xi 2+y2+z2$, and

*κ*=

*k*/(

*ρc*

_{p}).

Some assumptions and simplifications in deriving the Rosenthal equation are:

*The material’s physical coefficients such as k, ρ, and c*_{p}are independent of temperature. Using an average value provides a reasonable approximation and enables a closed-form solution to be obtained.The internal heat generation is neglected, i.e.,

*q*_{s}= 0.The workpiece material is homogeneous and isotropic.

When the powder bed is processed long enough, a quasistationary state is reached, that is, the temperature undergoes no change with time with respect to the coordinate system attached to the heat source, i.e., (

*ξ*,*y*,*z*).A point heat source rather than Gaussian distribution is used.

The effect of latent heat of fusion is negligible since the absorbed latent heat evolves later on.

Assumptions in deriving Eq. (5) are

−(ln(

*r***N*)/*r***M*) ≈ 0, that is,*r** ≈ 1/e*N*, where*r** is the value of*r*at the width of the melt pool,*M*=*u*_{x}/2*κ*, and*N*= (2*πk*(*T*_{m}−*T*_{0}))/*q*.*r***M*≫ 1.The approximation of

*q*is found to be improved by accounting for the zero-speed power in Eq. (4), that is, the first term on the right-hand side of Eq. (5).

Under typical PBF configurations, the first two assumptions are reasonably valid for all alloys except for the alloy AlSi10Mg.

### 4.2 Structure of the Hammerstein Model.

*P*(

*z*) that has unit DC gain (see Sec. 3). As shown in Fig. 4, the Hammerstein model upgrades

*P*

_{d}(

*z*) by replacing the constant

*c*with the nonlinear closed-form expression of the steady-state melt pool width

*f*(·). In Eq. (5), the values of parameters

*k*,

*ρ*, and

*c*

_{p}are to be determined. Substituting the equilibrium point (

*q*

_{0},

*w*

_{0}) to Eq. (5) gives $q0=\pi k(Tm\u2212T0)w0+e\pi \rho cp$$(Tm\u2212T0)uxw02/8$, that is,

*B*= 1/[

*π*(

*T*

_{m}−

*T*

_{0})

*w*

_{0}] is a constant. In Eqs. (5) and (6),

*ρ*and

*c*

_{p}are multiplied together and related to

*k*. Based on the first assumption in Sec. 4.1, we choose

*k*= 40 W/(m k) from Fig. 1. Substituting Eq. (6) in Eq. (5) yields

*P*(

*z*) and letting

*δw*

_{ssi}=

*w*

_{ssi}−

*w*

_{0}and

*δq*=

*q*−

*q*

_{0}. Certainly, due to simplifications in deriving Eq. (5), the direct solution Eq. (8) only works at specific input conditions. Under the input signal of

*δq*

_{10}= 20 W sin (2

*πft*

_{s}

*n*), where

*f*= 10 Hz,

*n*is the discrete-time integer time index, and

*t*

_{s}is the sampling time, we can tell from Fig. 5 that the output of the Hammerstein model (dashed line) deviates from the FEM result (solid line).

*w*

_{ssi}in Eq. (8) with a compensation factor

*α*(

*q*):

*α*(

*q*) is a quadratic function that passes through three points (60 W, 1) (i.e., no compensation at the equilibrium point), (80 W,

*α*

_{1}) (i.e., the maximum laser power), and (40 W,

*α*

_{2}) (i.e., the minimum laser power). We identify the parameters

*α*

_{1}and

*α*

_{2}, respectively, as 0.8507 and 1.1973 using the parameter estimation tool in matlab. The nonlinear least-square regression is used to minimize the sum of squared errors between the FEM data and the output of the updated Hammerstein model with compensation (solid and dashed-dotted lines in Fig. 5).

The compensated Hammerstein model is achieved by using Eq. (9) instead of Eq. (8) and letting *δw*_{ss} = *w*_{ss} − *w*_{0}. From Figs. 5 and 6, we can tell that this Hammerstein model (dashed-dotted line) gives a better approximation (41% more accurate in terms of root-mean-square (RMS) errors) of the system dynamics than the identified linear model *P*_{d}(*z*) at 10 Hz.More generally, as shown in Fig. 6, under different input frequencies, the compensated Hammerstein model consistently yields smaller RMS errors with respect to the FEM result than the linear model. Figure 7 illustrates the results of the compensated Hammerstein model and the FEM under the input frequency of 60 Hz. Following the procedures from Eq. (6) to Eq. (9), we can adaptively build the Hammerstein model for each specific equilibrium point and structurally draw the complete model map for the entire task space of PBF.

We compare in Fig. 8 how *δw*_{ss} changes with respect to *δq* under different modeling schemes. In the identified linear model *P*_{d}(*z*), the gradient of the dashed-dotted line that links *δw*_{ss} with *δq* is the constant *c* (Fig. 4). From Fig. 8, we can tell that *P*_{d}(*z*) is only a linear representation of the nonlinear Hammerstein model (solid line) near the equilibrium point. It is remarkable how *P*_{d}(*z*) identified from the FEM data coincides tangentially with the Hammerstein model derived from the governing equation. Next we will conduct the stability and robustness analyses to investigate when *P*_{d}(*z*) would fail in representing the Hammerstein model.

### 4.3 Stability and Robustness.

*P*

_{d}(

*z*) that are commonly used in practice. Let

*H*=

*P*

_{d}(

*z*)(1 + Δ), where

*H*is the Hammerstein model and Δ is the bounded model uncertainty. From Fig. 4, we have

*δw*

_{ss}·

*P*(

*z*) =

*δq*·

*cP*(

*z*)(1 + Δ), which gives

Standard robust-stability analysis indicates that a closed-loop system containing *P*_{d}(*z*) is stable if and only if both of the following conditions hold:

Nominal stability condition is satisfied, that is, the closed loop is stable when Δ = 0.

Robust-stability requirement is met by applying the small gain theorem [27]: for any frequency Ω in Hz, $|\Delta \u22c5T(ej\Omega Ts)|<1$, that is, $|\Delta |<1/|T(ej\Omega Ts)|$, where

*T*(*z*) is the complementary sensitivity function [1].

Note that |Δ| in Eq. (10) is positively correlated to the control signal *δq*, that is, more laser power deviation from the equilibrium point yields a larger |Δ|. Under a certain frequency Ω, we need to make sure the maximum |Δ| is less than $1/|T(ej\Omega Ts)|$. When the condition is violated, *P*_{d}(*z*) will no longer be a valid representation of the Hammerstein model.

### 4.4 Control Implementation.

Although the focus of this paper is on the modeling of the complex physics in PBF, we have additionally validated the proposed model in closed-loop controls. Limited in space, we will briefly discuss the key concept. Figure 9 presents a typical feedback loop when applying the Hammerstein model. An *f*^{−1}(·) block is added and connected with the block of the Hammerstein model. Combining these two blocks together yields the linear model *P*(*z*). It is thereafter a standard practice to design the control algorithms for *P*(*z*). To find the inverse of the nonlinear element, one approach is to use a high-order polynomial to approximate *f* (solid line in Fig. 8). Besides, Ref. [23] proposes an approximate method for the cases when *f* is not invertible.

## 5 Conclusion

In this paper, we first build a FEM to simulate the thermal fields during the PBF process. Using the FEM data, we identify around the equilibrium point the linear system model from the laser power to the melt pool width. In addition, deriving from the Rosenthal equation, we reach a nonlinear closed-form expression of the steady-state melt pool width. Concatenating the nonlinear expression to the identified linear model, we develop the main Hammerstein model that captures more of the convoluted thermomechanical dynamics of PBF. We prove that the Hammerstein model gives a better approximation (e.g., 41% increasing at 10 Hz) of the FEM result than the linear model. From there, we analyze the stability and robustness properties of the models and present a generic control scheme for the Hammerstein model.

## Acknowledgment

This material is based upon work supported in part by the National Science Foundation under Award No. 1953155.

## Conflict of interest

There are no conflicts of interest.

### Appendix

Defined parameters of the Finite Element Model

Parameters . | Value . | Parameters . | Value . |
---|---|---|---|

Powder bed size | 15 mm × 15 mm × 50 μm | Material | Ti6Al4V |

Substrate size | 15 mm × 15 mm × 5 mm | Track length L | 5 mm |

R | 220 μm/2 | Time-step T_{s} | 0.5 ms |

Absorptance | 0.25 | Emissivity | 0.35 |

T_{sol} | 1873 K | Scan speed u_{x} | 100 mm/s |

L_{f} | 295 kJ/kg | Laser power P | 60 W |

T_{0}/T_{m} | 293.15 K/1923.15 K | ϕ_{0} | 0.4 |

h_{c} | 12.7 W/(m^{2} K) | k, c_{p}, and ρ | See Fig. 1 |

Parameters . | Value . | Parameters . | Value . |
---|---|---|---|

Powder bed size | 15 mm × 15 mm × 50 μm | Material | Ti6Al4V |

Substrate size | 15 mm × 15 mm × 5 mm | Track length L | 5 mm |

R | 220 μm/2 | Time-step T_{s} | 0.5 ms |

Absorptance | 0.25 | Emissivity | 0.35 |

T_{sol} | 1873 K | Scan speed u_{x} | 100 mm/s |

L_{f} | 295 kJ/kg | Laser power P | 60 W |

T_{0}/T_{m} | 293.15 K/1923.15 K | ϕ_{0} | 0.4 |

h_{c} | 12.7 W/(m^{2} K) | k, c_{p}, and ρ | See Fig. 1 |