The temperature inside the part in powder bed fuison

Powder Bed Fusion (PBF) is a type of Additive Manufacturing (AM) technology that builds parts in a layer-by-layer fashion out of a bed of metal powder via the selective melting action of a laser or electron beam heat source. The technology has become widespread, however the demand is growing for closed loop process monitoring and control in PBF systems to replace the open loop architectures that exist today. This paper demonstrates the simulated efﬁcacy of applying closed-loop state estimation to the problem of monitoring temperature ﬁelds within parts during the PBF build process. A simpliﬁed LTI model of PBF thermal physics with the properties of stability, controllability and observability is presented. An Ensemble Kalman Filter is applied to the model. The accuracy of this ﬁlters’ predictions are assessed in simulation studies of the temperature evolution of various test parts when subjected to simulated laser heat input. The signiﬁcant result of this study is that the ﬁlter supplied predictions that were about 2.5x more accurate than the open loop model in these simulation studies.


Introduction
Powder Bed Fusion (PBF) belongs to a class of technologies known as additive manufacturing (AM). Commonly referred to as "3D printing," these technologies have rapidly grown in popularity and market size due to their ability to produce near net-shape parts of complex geometry, with engineering properties meeting or exceeding those produced by conventional techniques, while removing the majority of the overhead costs normally associated with production [1][2][3].
The PBF process iteratively builds three-dimensional parts out of layers of metal powder, using a build cycle consisting of three stages: 1) sweeping a thin layer of powder over a base of metal feedstock or previously-applied powder, 2) selectively melting a pattern of desired geometry into the powder by application of a high-powered laser or electron beam, and 3) lowering the build − platform in the z direction to accommodate a fresh layer of powder. Schematics of PBF are given in Fig.  1 along with typical input and output channels that are available to the controls engineer.
The PBF process is not without flaws. It is well-documented in the literature that components manufactured with PBF display high levels of residual stresses ( [6,7]), porosity ( [8][9][10]) , and anisotropy in material properties ( [2,[10][11][12][13][14]), and that these defects are a direct consequence of the heat transfer throughout the part, which is manifested by cooling rates within the part interior. The application of thermal model-based process monitoring and control techniques could help detect and mitigate these defects. Process monitoring suites offered by commercial PBF systems, as reviewed in [15], typically assess the presence of defects based on exhaustive calibrations that "train" the model to accept certain measured values as defect-free. Predictive model-based process monitoring with a minimum of necessary calibration remains elusive. This paper advances the goal of predictive model-based process monitoring for PBF. We con-tinue the work shown in [16] and present a linear time-invariant (LTI) state space model of PBF conductive heat transfer physics with established stability, controllability, and observability. This model is based on first principles and thus requires a minimum of training/calibration to perform. We reduce the model and express it in discrete time, and then demonstrate the application of a state estimator known as an Ensemble Kalman Filter to the reduced-order model. We conduct simula-tion studies of this state estimator by constructing models for simulated test parts when subjected to simulated laser heat input. We assess the performance of the state estimator by comparing the estimator accuracy with that of the open loop model, with respect to a reference simulation repre-senting a "true" evolution of the part thermal history. We show that the state estimator provided predictions that were approximately 2.5x as accurate as the open loop model, which we believe constitutes justification for further research into this topic.

PBF model assumptions and LTI model construction
In this work we build upon the spatiotemporal model constructed in [16]. What follows here is a statement of the method results, readers interested in a complete description on the model construction and its properties should consult [16]. We examine a simplified model of PBF thermal physics in which only Fourier conduction within the fused material is considered. The Fourier conduction BVP is stated in (1).
Here, K represents the material thermal conductivity, c represents the material specific heat, and ρ represents the material density. V represents the domain spanned by the (possibly unfinished) build, Λ collects all faces of the build in contact with the base plate, Ω collects all surfaces of the build exposed to the environment, and Γ collects all remaining faces. T0 is the assumed-isothermal temperature of the base plate.
The BVP (1) was transformed into the linear state-space thermal model shown in (2). We use the transformation described in [16], which converted the PDE into a set of coupled ODEs via the Finite Element Method (FEM). Each ODE governs the evolving temperature at a single node in the FEM mesh, and the system input was quantized to hold a constant value over each element surface belonging to Ω. In Thermal Model (2), x collects the temperature signals at all nodes in the mesh, A maps the degree of conductive heat flow between nodes (0 for nonadjacent nodes), B maps the degree of distribution of laser energy input onto nodes located on Ω, and C selects the nodes belonging to Ω as system output in keeping with our assumption that only exposed faces of the build are available for measurement.
It was shown in [16] that Thermal Model (2) is unconditionally stable, stabilizable, and de-tectable, and that it is both structurally controllable and observable provided that at least one node in the FEM mesh exists on the exposed build surface.

Model order reduction via balanced realization
Thermal model (2) was developed to mitigate the problem of model scale when attempting to represent PBF physics on the macroscale, however the quantization of the system heat input can produce cumbersome node counts and therefore system sizes. Model order reduction (MOR) is necessary to reduce the impact of these issues. We have chosen to perform MOR via residualiza-tion [17,18].
The residualization algorithm requires stability, controllability and observability of the system, which we have shown in [16]. The algorithm begins by performing the linear state transformation ✓ ∈ ∈ z2 z = Tx, which puts the system in its' balanced realization. The user then selects the first r (largest) Hankel Singular Values (HSVs) of the system, which are demonstrated in (3) [19,20]. Each HSV is λi (WcWo); λi denotes the i th eigenvalue of a matrix, and Wc and Wo are the controllability and observability grammians, respectively.
Partitioning the HSVs in this manner also partitions z into two groups: z1 R r , which consti-tute the "dominant" modes in the system input/output dynamics, and z2 R n−r , which constitute the negligible modes. The partitioned system takes the form shown in (4) The residualization algorithm assumes that the "weak" modes stored in z2 operate at quasi-steady state. In other words, it assumes that on the time scales of interest, ż2 = 0. This assumption allows for the algebraic solution of z2, which reduces (4) to the form shown in (5).
With residualization, one only needs to solve r coupled differential equations instead of n, and the original state x may be reconstituted during postprocessing by algebraically calculating z2 and performing the inverse transformation x = T −1 [z1, z2] t . By [18,21], the (structural) controllabil-ity/observability of (2) implies that (5) is also (structurally) controllable and observable.

Discretization of continuous-time model
This study utilized a Kalman filter to estimate the state of Thermal Model (2), which required it to be implemented in discrete time. The continuous-to-discrete time conversion of model (2) for a discrete time step Δt is shown in (6): The discretization scheme provided in (6) is based on the complete discretization method given in [22] for a system with constant parameters. The approximation of the matrix exponential e A r Δt in (6) is based on the Bilinear Transform as given in [23].

Ensemble Kalman filter (EnKF) implementation
Effective implementation of a Kalman filter requires knowledge of the covariances of the pro-cess and measurement noise of the system under consideration. The process noise covariance is denoted as Q, and the measurement noise covariance is denoted as R. We intend to for model (6) to be implementable for arbitrary build layer geometry V and under arbitrary external loading q(v, t). These constraints make the prediction of Q intractable, while also limiting the utility of experimentally determining Q due to the near-infinite variety of allowable geometric and loading conditions.
To overcome this limitation, we seek a means of approximating Q and R in-situ from the measured data y and reconstructed state estimate x. The Ensemble Kalman filter detailed in [24] provides the optimal means of accomplishing this goal. Appendix A gives a brief overview of the EnKF algorithm, readers interested in a complete description of EnKF theory should consult [24].
Our process model as described in (6) differs from the structure typically assumed for the oper-ation of conventional Kalman filters (13) due to the presence of direct feedthrough. Accounting for independent process noise wk N (0, Q) and measurement noise vk N (0, R), our stochastic process takes the form shown in in (7). The subscript 1 attached to z1,k indicates that the filter is estimating z1 of the reduced order system as defined in (5), (6).
We follow the modified Kalman filter architecture supplied in [25], which acts to provide simultaneous unbiased minimum-variance estimates for the system input uk and state z1,k in the presence of systems with direct feedthrough. To this end, we define an ensemble of inputs corresponding to Z1,k, Uk = u 1 , u 2 ,..., u N . [25] assumes that the values of all u i are unknown. We simulate this condition by making the process input uk stochastic via direct injection of the system process noise into the input for each ensemble member: In (8), each ensemble member of Uk is defined as 1,k−1 represents the filter estimates for the member at time step k 1. Each member represents a sample from the randomly-distributed system inputs and measurements, respectively. This stochastic treatment of the process input reflects the uncertain nature of PBF processing conditions discussed in [16]. (8) can be expressed in more conventional form by rearrangement: As shown in [26], The filter architecture specified by [25] has three steps. When combined with the EnKF archi-tecture specified by [24], these steps take the following form: 1. Predict. Each ensemble member updates its ensemble values z i according to (8), using u i as inputs, while collecting measurements y i . The ensemble is used to calculate sample k k ¯ f ¯ ¯ f estimates P k and Rk according to (15) and (16), respectively. Here, P k represents the sample estimate for the covariance associated with the estimated state error.
3. Estimate state. Minimum-variance estimates for ensemble member state estimates ẑ i calculated according to an expression similar to (19). The procedure is outlined below: In (9) and (10), Mk and Kk are calculated with pseudoinverses because R k becomes singular if the number of measurements p is greater than the number of ensemble members N, as noted in [24].

Filter algorithm summary
Construct a reduced order linearized model from the governing FEM-discretized heat trans-fer equation as done in [16].
• Express the model in discrete time as done in (6).
• Define initial temperature distribution throughout the part, express in terms of z1,k=0.
• Define an ensemble of N parallel instances of model (8), each having ensemble members z i , y i , and u i . ble members to generate all z 1,k . Ensemble member measurements y k are constructed by corrupting the system measurement yk with independent instances of white noise: Collect into ensembles Z1,k, Yk, and Uk, respectively.
• Compute P f and R k from (17) and (18).
• Construct estimated ẑ1,k by taking the sample average of Ẑ k .

Case Studies Description
Two case studies were conducted to assess the performance of the state estimator constructed in Sections 2 and 3. These tests constitute a continuation of the work performed by the authors in [27], in which data from assumed-accurate simulations were used in place of physical test parts. We adopt the same basic "virtual test" procedure here.

Test procedures
The simulated tests utilized the following procedure: 1. Construct a linearized state space model corresponding to test parts according to the proce-dures outlined in Section 2 and Fig. 2. The temperature along isothermal boundary Λ was set to 0. All test parts used material properties corresponding to Aluminum ore, tabulated in Tb. 1. Initial temperature was uniformly 0.

2.
Construct a heat conduction simulation in ANSYS utilizing the same test part geometry, mesh, and isothermal boundary, but with the full nonlinear treatment of the heat source. Initial temperature was uniformly 0. These simulations were used as surrogate "true" data, to test the amount of error incurred by the linearization process, tk.
(a) Tb. 1 shows that the material properties used to construct the LTI model and those used to construct the reference ANSYS simulation differed substantially. This was done to assess the effectiveness of the filter in purging modeling errors from the predicted ˆt in a worst-case scenario. The material properties used to construct the reference ANSYS data were those corresponding to Aluminum at its' melting point, while those used to construct the linearized LTI model were those corresponding to Aluminum at room temperature, and therefore represent the maximum possible modeling error. (b) Process noise, visible in Fig. 3, was added to the ANSYS laser data, to reflect the inevitability of uncertainty in the time-varying process inputs during PBF. This contin-uous time noise had power equal to 1 10 5 mW for Part (a) and 2 10 5 mW for Part (b).
Construct uk for all test parts as Gaussian laser beams of a given power P and variance σ 2 moving over the top surface Ω of the parts in a raster scanning path with average speed v, incorporated into the LTI model according to [16].

Run the Ensemble Kalman
Filter algorithm with the procedure described in Section 3.3, using the constructed uk as an input and using ANSYS data along surface Ω of the test parts from Step 2 as system measurements yk. Recover temperature estimations xk.
(a) Ensembles with N = 100 members were used for both parts.
(b) The continuous-time process noise power injected into the input ensemble members uk was 2 10 9 mW for Part (a) and 2 10 12 mW for Part (b). These noise powers were high because testing showed that EnKF performance improved as the injected noise power tended to infinity, since doing so increased the difference z i − zk for all ensemble k ¯ f members in (17), therefore increasing the accuracy of computing P k for large N. The continuous-time measurement noise power injected into yk for both parts was 1 K.

5.
Define the EnKF state estimation error (the "closed loop" estimation error) as Error(t)= xk − xk. Plot and animate this error.

Test parts
The two test parts utilized in this study are depicted in Fig. 2. Information pertaining to the mesh and system size for these test parts is shown in Tb. 2. The slight discrepancy between node count and n is due to temperature-constrained nodes along Λ being removed from the system by ANSYS. Table 1: Comparison of material properties used to construct LTI model data and reference ANSYS data for all test parts [28] 2.7 × 10 −9 2.5 × 10 −9  Fig. 2 illustrates the basic path of the laser beam across surface Ω for both test parts, as well as its nominal power P. The laser power for both parts was distributed across Ω according to theprocedure outlined in [16]. The time-varying laser centerpoints were calculated according to the expressions below.

Test loading conditions
The laser centerpoint x(t) for Part (a) was calculated according to the expression: x max −x min 2 2 Where v = 954 mm/s, and xmax and xmin were taken from the part geometry as shown in Fig.  2. The end time of the simulation was set to be tfinal = 4 ms. The variance of the beam was set to be σ 2 = 0.01.
The laser centerpoint {x(t), y(t)} for Part (b) was calculated according to the expression: Where v = 954 mm/s, and xmax, xmin ymax, and ymin were taken from the part geometry as shown in Fig. 2.
The end time of the simulation was set to be tfinal = 20 ms. The variance of the beam was set to be σ 2 = 0.0225. Fig. 3 displays these loads for sample time steps in the simulation.

Simulation Results
Fig. 4 plots the OL model error (ErrorOL(t)) and CL estimation error (Error(t)) for both test parts. It reveals two important aspects of the filters' performance. The addition of the EnKF reduced the model error in both test parts by approximately a factor of 2.5. Additionally, as Fig. 4 shows, ErrorOL(t) for Part (a) was unbounded. This phenomena is due to the geometry of Part (a). The Neumann boundary condition that models surrounding powder insulating the build [29] trapped heat inside the "arm" of Part (a), therefore causing a monotonic temperature increase in that region. The rate of increase was a function of the thermal diffusivity of the material being modeled, which was constructed to contain the maximum possible error. This disparity in rates of temperature increase produced an unbounded ErrorOL(t). As Fig. 4 shows, the filter stabilized this error such that it was bounded. Fig. 4 shows that both ErrorOL(t) and Error(t) appeared to be strongly periodic for Part (b). This phenomena is explained by Fig. 5, which illustrates the EnKF state estimation error for both test parts for sample time steps in the simulation. One clearly observes from the figure that the region of maximum error closely "followed" the laser centerpoint. This result was expected, given that the region nearest to the heat source produced the most extreme thermal response and thus deviated the furthest from the operating point of the linearization. What appears to be periodic behavior in Fig. 4 is actually the superimposition of select few state components in Part (b)ie nodal temperatures corresponding to physical locations in the geometry -experiencing a "pulse" of substantial error as the laser beam passed over them at staggered moments in time, before the heat diffused away and the regions decayed back to approximate ambient temperature with little error.
Fig. 5 also demonstrates that Error(t) oscillated about the laser point center. Error(t) in the "neck" of Part (a) oscillated between positive and negative as the laser beam transitioned from close to far away from the x = 0 plane. Fig. 5 shows that this tendency was also present in Error(t) in Part (b), with regions of positive and negative estimation error alternating radially outward from the laser centerpoint. It is currently unknown what produces this wave-like pattern in Error(t).

Conclusions
This paper showed the feasibility of applying state estimation to the problem of acquiring internal temperature field predictions for parts being manufactured via PBF. It demonstrated the application of an Ensemble Kalman Filter to enforce discrete-time reduced LTI model accuracy in the presence of uncertain system parameters and uncertain model error/noise statistics. It demon-strated that the implementation of the Ensemble Kalman Filter in simulation studies improved the accuracy of these model predictions by a factor of 2.5, even in the presence of worst-case model-ing error. These results show that pursuing a controls-based approach to improving the accuracy of simplified predictive models of PBF thermal physics holds great promise and warrants further research.
We intend to pursue several avenues of research in light of the results of this test. First, the performance of the filter will be tested experimentally against the temperature evolution of test coupons that are subjected to varying load conditions in open source PBF machines. We intend to explore the theoretical performance limitations of our reduced model, Ensemble Kalman Fil-ter approach. We also intend to research model reduction techniques available for more complex, time-varying models of the system. We anticipate that this research will present a marked contri-bution toward the goal of realizing closed-loop, model-based monitoring of the PBF process.

Acknowledgements
Financial support was provided by the member organizations of the Smart Vehicle Concepts Center, a Phase III National Science Foundation Industry-University Cooperative Research Center (www.SmartVehicleCenter.org) under grant NSF IIP 1738723. The authors acknowledge technical support from ANSYS.
unknown true value z t . z t would equal the value of zk obtained from (13) if no noise were present.

k k
The random variable P f : zk is distributed with the state error covariance Similarly, the measurement yk is treated as a random variable with some corresponding un-known, noise-free true value y t . It is clear that yk is distributed with measurement error covariance equal to R: As shown in [24], the EnKF defines an ensemble of N parallel instances of (14), denoted as Zk = z 1 , z 2 ,..., z N , with corresponding measurement ensemble Yk = y 1 , y 2 ,..., y N . Each ensemble member z k , y k is generated by running (14) with independent instances of wk and vk. Therefore, Zk, Yk collect N samples of the random variables zk, yk. By defining sample averages zk = 1 ∑ z i and ȳk = 1 ∑ y i , one may construct the sample estimations of P f and R as defined in N k N k k (15) and (16), respectively: [24] note that ensembles of size N = 100 or greater estimate the true values of P f and Rk to consistently acceptable accuracy. Model order reduction is absolutely essential to avoid unreason-able computational burden when running several dozen concurrent models.
Having defined P f and R k , [24] runs the standard Kalman filter update for every ensemble member: This process may be represented compactly by operating on the ensembles Zk and Yk: