Time-efficient Surrogate Models of Thermal Modeling in Laser Powder Bed Fusion

Two time-efficient surrogate models are proposed to emulate the nonlinear heat equation in the context of laser powder bed fusion, the performance of which is compared in accuracy and online execution time. Fast-computed numerical solvers are critical in implementing the digital twin framework in the additive manufacturing process addressing one of its main open problems: lack of quality assurance. The first surrogate model is the reduced Gaussian process emulator. It is a data-driven model equipped with a nonlinear dimension reduction scheme and manages to predict temperature profiles almost instantly (around 0.036s on average) with an accuracy of 95% for 99.38% of tests. Another surrogate model is the sketched emulator with local projection. It projects the accurate but high-dimensional finite element method solution on a low-dimensional basis and then bypasses the majority of costly computations for the temperature-dependent matrices in the projected model by randomized sketching. It has higher accuracy (97.78% of tests with relative errors below 1%) while spending comparably more time online (around 42.23s on average). Although compared with the finite element model both surrogates promote time efficiency with some minor controlled compromise in accuracy, the reduced Gaussian process emulator enables real-time implementation while the sketched emulator with local projection offers comparably higher levels of accuracy. A series of numerical experiments are carried out, which assumes a three-layer printing process with a fixed laser beam trajectory using a small number of printing control parameters as inputs, namely the laser power, scan speed, and time coordinates. Both surrogates are also principally feasible in other thermal-driven additive manufacturing to obtain better quality assurance with techniques like uncertainty management and closed-loop control.


Introduction
Laser powder bed fusion (LPBF), as a prevalent additive manufacturing (AM) technique in metal fabrication, successively spreads and selectively melts thin layers of powder to fabricate physical parts from 3D geometrical designs [1]. Though it is superior in realizing complex or personalized design resulting in increased application in fields like aerospace [2], automotive [3], and medical industry [4], quality assurance has been a persistent problem that makes digital twins of LPBF that help to approach the ''first-time-right'' goal receiving more and more attention [5]. As a thermal-driven AM, the thermal model of LPBF is a critical step in predicting the part properties such as microstructure and residual stress [6][7][8], according to which we can improve the part quality by designs like process optimization and closed-loop control [9,10]. Extensive researches with experimental validation have focused on the numerical solver with finite element method (FEM) such as the thermal model validated by an experiment with AlSi10Mg powder and argon atmosphere in [11,12] and the * Corresponding author.
nonlinear transient heat equation with phase changes and anisotropic thermal conductivity validated by the AM benchmark experimental set AMB2018-02 in [13]. Though the numerical simulators with FEM are precise, they are computationally expensive due to the high dimensionality caused by fine spatial-temporal discretization and the high nonlinearity caused by temperature-dependent thermal properties and boundary conditions. The high nonlinearity, however, means that we need to solve several high-dimensional linear equations to get the converged temperature result at a single time instant [14]. The expensive time cost detrimentally hurdles its application both before and during printing processes, thus a swift but accurate heat transfer modeling of LPBF becomes a valuable research area.
To address the challenge of time consumption, various approaches have been proposed. One type of approach is purely data-driven. In [15], a recurrent neural network was trained to predict thermal profiles in the directed energy deposition process from different geometries. However, it required 250000 data points taking extensive https://doi.org/10.1016/j.addma.2022.103122 Received 25 March 2022; Received in revised form 15 July 2022; Accepted 30 August 2022 time and storage in data preparation and training. Another type of approach focus on model order reduction (MOR) since the thermal model of LPBF is a large-scale full order model (FOM) with significant redundancy. MOR can be assorted in three types: structural, datadriven, and projection-based approaches. Structural approaches refer to adaptive meshes [12,16,17] and substructure coupling [18,19]. In [16], an adaptive mesh refinement strategy was suggested to refine the mesh around geometric components while having coarser mesh for other areas. Jensen et al. described the dynamic analysis with substructures coupling internal dynamic behaviors and fixedinterface normal modes [18]. The data-driven MOR uses machine learning algorithms to find and/or predict low-dimensional latent representations of high-dimensional training outputs. Nikolopoulos et al. integrated a convolution autoencoder and a feed-forward neural network as a surrogate of partial differential equations (PDEs) along the entire time history [20]. The projection-based MOR projects highdimensional PDEs into their reduced form where projection bases are obtained by methods such as moment matching [21], balanced truncation [22], and principal orthogonal decomposition (POD) [23]. Though the abovementioned attempts could expedite a thermal simulation, some obstacles remain. Data-driven approaches of exhaustive mapping generally rely on the representative quality of training data and normally need a large amount of training data which might be prohibitively expensive to generate. Projection-based approaches are challenging in finding proper and effective bases, and require additional approximations for nonlinear models such as discrete empirical interpolation method (DEIM) [24]. Nevertheless, a way out of the expensive time cost in thermal modeling of LPBF is surrogate models. The so-called surrogate model (SM) also known as emulator makes use of data generated by a high-fidelity physics-based model (the FEM in our case) and then swiftly predicts online with given design points, thereby inexpensively emulating the simulation [25].
Each SM individually balances model accuracy and online execution time. In this paper, we propose two SMs emulating the FEM with reasonable but different trade-offs. The first SM, the reduced Gaussian process (GP) emulator, learns a low-dimensional representation of the relative distance between temperature profiles which is later extrapolated to the prediction of high-dimensional temperature. The second SM, the sketched emulator with local projection, is a speedy numerical solver with FEM. It uses the relative distance predictor in the first SM to intelligently subsample the training temperatures and form effective local projection bases. The projected model is then sketched to suppress the time cost due to nonlinearity. Both SMs are time-efficient, but their run times and prediction accuracy are affected in different ways. The reduced GP emulator predicts temperatures very fast (0.036 s on average) with 99.38% of tests having relative errors below 5%, while the sketched emulator with local projection manages to predict 97.78% of tests with relative error less than 1% in an average execution time 42.23 s saving more than 82% of time required by FEM (4.23 min on average). Since both surrogates exhibit very good performance, we either choose one of them or indeed combine them to tackle problems requiring large and/or rapid thermal simulations in LPBF.

Notation and paper organization
We express matrices in capital letters, vectors and continuous functions in small case letter. For a matrix , is the ( , )th element, the th column is denoted as * , and * denotes its th row. For a vector , the th entry is denoted by .
In the next section, we describe the nonlinear transient heat equation model as well as its boundary conditions pertaining to LPBF and then outline its high-fidelity numerical solution with FEM. The section afterward illustrates the methods of two proposed SMs. The method of reduced GP emulator is sequentially depicted as the subsampling based on relative distance predictions and the approximation of highdimensional temperatures based on subsampled temperatures. The method of the sketched emulator with local projection is depicted by three subsections which contain local projection using relativedistance-based subsampling, the approximation of temperaturedependent matrices with randomized sketching, and the reconstruction of high-dimensional temperatures after solving the projected and sketched model. Then, the section of results and discussion follows, which details the numerical experiments comparing the performance of the two SMs. The prediction accuracy is compared in temperature profiles, thermal histories, and melt pool sizes, while the time efficiency is compared by average execution time. In the end, we have the section containing conclusions of designed models.

Numerical thermal modeling of LPBF
During the LPBF process, as shown in the schematic of Fig. 1, a high power-density laser beam selectively melts and fuses a thin layer of powder along a preset scan path. The printed object is then formed by bonding adjacent tracks and layers as new layers of powder are successively superimposed. Its thermal modeling entails a heat equation over a vertically extended domain, the boundary conditions of which consider a Gaussian heat flux, convection and radiation heat loss, and a temperature-controlled platform. The powder bed as a multi-layer spatial domain of heat equation is discretized by adaptive tetrahedron meshes, which yields a numerical solver with FEM.

Governing equations
The heat transfer in a 3D computational domain ∈ R 3 is governed by the nonlinear heat equation in Eq. (1). It is nonlinear due to the dependence of material properties on temperature. Let ∶= ( , ) be the temperature at spatial coordinate ∶= [ 1 , 2 , 3 ] in a 3D Cartesian frame, and temporal coordinate ∈ [0, ]. If , , and ⃗ are respectively the temperature-dependent density, specific heat capacity and thermal conductivity tensor of the material then the temperature in the domain is known to satisfy where the temperature dependence of the thermal properties is derived from material-specific experimental data. The thermal conductivity ⃗ is modeled as an anisotropic tensor field, and the latent heat effect during a solid-liquid phase change is introduced in terms of the model of specific heat capacity [13]. Further details of the models of nonlinear thermal properties are included in Appendix A. The heat is imparted to the domain by means of a laser beam, modeled as a moving Gaussian heat source [11] acting on the top surface of the domain . In Eq. (2), is the heat absorptivity, is the laser power, and is the effective laser beam radius, all three assumed invariant in space and time. In the same equation is the position of the center of the laser beam at time as it moves with a constant scan speed . The laser energy is mainly absorbed by the powder resulting in melt pools, while some of it escapes in the atmosphere as a heat loss satisfying wherêis the outward unit normal on the top surface and side boundary , while the respective convection and radiation constituents in = + are with ℎ > 0 is the heat convection coefficient, is the Stefan-Boltzmann constant, is the emissivity, and the ambient temperature [26,27]. To improve the final object's quality, such as mechanical properties,  microstructure and residual stress, a temperature is applied on the temperature-controlled bottom surface hence a Dirichlet boundary condition is imposed at all times. Finally we assume the initial condition so that Eqs. (1)-(6) admit a unique solution ( , ) ∈ × [0, ].

Finite element method
To solve the nonlinear equations as Eqs.
(1)-(6) using the FEM, we respectively discretize the domain in space and time in linear tetrahedral elements and a number of time instants { ≐ } =0 for a given small . The imposed boundary conditions and the temperature −1 yields a model with degrees of freedom (DoF) for the temperature field in the form of the nonlinear system ( ) = ( ), = 1, 2, … where the temperature-dependent system matrix and right hand side vector at time have dimension which is typically very large. Details of the construction of these are given in Eqs. (B.1) and (B.2) in Appendix B.1, the matrices and vectors in which are defined as the integrals in Table B.1 and are numerically approximated with Gaussian quadrature rules [28,29] as in Table B.2 in Appendix B. The nonlinear heat equation as Eq. (7) is solved iteratively using the Picard algorithm [14] as we outline in Algorithm B.1 of Appendix B.4.
The three most time-consuming parts of the numerical solver with FEM are the high dimensionality caused by fine spatial and temporal discretization, the number of Picard iterations required to converge, and the computation of the FEM integrals involving high-order polynomials using Gaussian quadrature. To address these issues, two timeefficient SMs are proposed to significantly reduce the online execution time by exploiting offline data sets to minimize online computations.

Time-efficient surrogate models
The prevalent computational complexity associated with the highdimensionality and non-linearity of FEM renders it unsuitable for realtime simulation. To decrease this computational complexity, we propose two time-efficient SMs: the reduced GP emulator and the sketched emulator with local projection. The reduced GP emulator first builds a fast-computed data-driven model of relative distances between temperature profiles and then extrapolates the relative distance predictions as high-dimensional temperature predictions. In the sketched emulator with local projection, we propose low-dimensional local projections and randomized sketching to expedite FEM and form a SM. The local projection bases themselves are generated with the help of the relative distance predictor in the reduced GP emulator. Both SMs are significantly quicker than solving the FEM, and either of them may be used depending upon the desired trade-off between runtime and accuracy. Since a part of the reduced GP emulator is used to compute local projections for the sketched emulator, we begin this section by presenting the reduced GP emulator.

Reduced GP emulator
Data-driven approaches are frequently used in SMs where the result is obtained in near real-time as a direct prediction [20]. What hinders machine learning algorithms from emulating the thermal modeling in LPBF is the high dimensionality of thermal fields and the prohibitively expensive time in generating a large training data set. To address these challenges, we propose a reduced GP emulator where a subsampling scheme based on relative distance predictions is designed with Gaussian process regression (GPR) and nonlinear MOR, and a swift predictor is designed to extrapolate subsampled training temperatures and distance predictions to high-dimensional temperatures. The framework of the reduced GP heat emulator is as Fig. 2 where the module of subsampling and high-dimensional temperature predictor are respectively elaborated in Sections 3.1.1 and 3.1.2. In addition, the subsampling module is also used in the module of local projection demonstrated in the framework of the sketched emulator with local projection as Fig. 3 .

Subsampling
The high-dimensional temperature of a given test input can be predicted by a linear combination of some training temperatures deemed as the closest ones to this prediction, in which the weights of the linear combination are inversely proportional to relative distances between the temperature prediction and the selected training temperatures. We thus develop a subsampling scheme to select training temperature snapshots based on a relative distance predictor mapping input parameters to relative distances.
As a data-driven method, its prediction accuracy relies on finding representative training data set which is generated by the heat simulator with FEM. The training inputs are selected from the controllable parameters of a LPBF process including laser power, beam size, scan speed, time, preheating temperature, and so on. In this paper, the inputs matrix ∈ R × contains training points for a triplet ( = 3) of laser power , scan speed , and time . For each * ≐ [ , , ] , one can readily run the heat simulator with FEM as Eq. (7) by tuning the appropriate boundary conditions for times = 1 , … , until we obtain ( * ). Repeating the process for ≫ 3 sampling points produces a corresponding temperature snapshots matrix ∈ R × where * = ( * ) and is typically chosen to be smaller than . To make the learned parameters in GPR scale-invariant, we normalize the input quantities in the rows of to be in the [0, 1] rangẽ * = * − min( * ) max( * ) − min( * ) , As for the feature engineering of the training output, we approximate the pair-wise distances of temperatures in with nonlinear dimensionality reduction [30]. The relative distance between * and * for , = 1, … , , denoted as in the pair-wise distance matrix ∈ R × , is the shortest distance between the th and th node of a weighted -neighborhood graph. The construction of this weighted graph and the computation of its shortest path distance are further detailed in appendix C. We then find a reduced representation of by a rank truncation of its singular value decomposition (SVD) which is

X. Li and N. Polydorides
where = is the × reduced representation after featuring. Here, ∈ R × , ∈ R × , and ∈ R × are respectively the left, right singular vectors, and the diagonal matrix of singular values.
∈ R × respectively contain the left, right singular vectors, and singular values corresponding to the largest singular values in . To guarantee the accuracy of reduction, the value of is chosen to ) > 99%. After normalization and featuring, the resulting training data set is {̃ * , * } =1 which represents a mapping from -dimensional normalized inputs to -dimensional reduced outputs with ≪ .
As the preprocessed training inputs̃and outputs are passed to the GPR module, GPs are respectively trained to fit {̃, * } =1 . GP is promising in providing accurate predictions with corresponding variance using less lengthy training comparing with some other machine learning algorithms like neural networks. The multiple inputs and single output regression betweeñ∈ R × and * = * − 1 * ∈ R 1× for = 1, … , is a discrete GP with zero mean and positive definite covariance matrix ( ) ∈ R × , which is * ≈ (0, ( ) ), where the covariance matrix ( ) ( * ) arises from the discretization of a covariance function (kernel) ( , ′ ; ( ) ) with hyper-parameters ( ) . More specifically, outlining the smoothness of GP as a squared exponential function of [31]. The strictly positive hyper-parameters * = [ ( ) 0 , ( ) 1 , … , ( ) , ] ∈ R +2 is optimally found by the maximum likelihood estimation. We hereafter respectively denote the th column of the strictly positive matrix̂∈ R ( +2)× and the covariance matrix̂( ) as the trained hyper-parameters and the trained covariance matrix of the th GP for = 1, … , . For a given test input ∈ R that is not in the training data set, we first normalized it as̃. Then, the prediction of its reduced representation̂∈ R and the corresponding variances Var(̂) arê where the matrix ∈ R × and the matrix ∈ R × are respectively It is worth noting that only the matrix is necessary to be computed online as it depends on the normalized test input̃, while the more computationally expensive matrix can be pre-computed offline. Therefore, the prediction of̂as Eq. (13) can be implemented instantaneously. We then reconstruct the full relative distancê∈ R aŝ where is the -dimensional orthonormal bases in Eq. (9), and̂is the relative distance prediction between the temperature ( ) and the th training temperature * . In other words, the smaller̂is the closer the training temperature * might be. We thereby could select a subset of training temperatures that are deemed to be more similar to the final temperature prediction ( ).

High-dimensional temperature prediction
According to the prediction of relative distanceŝ, we selectc losest training data in as̄∈ R ×̄w herēis normally selected as a comparably small value (like 10 or 20). To extrapolate the high-dimensional temperature we compute the weights ∈ R̄of̄as where the th column̄ * weighted by is selected from the th column in the original matrix of training temperatures , and eventually the temperature prediction of the test input iŝ

Sketched emulator with local projection
While the reduced GP emulator is entirely data-driven, the sketched emulator with local projection preserves the structure of the original heat simulator with FEM as Eq. (7) but expedites it by local projection and randomized sketching. Local projection targets at addressing the issue of high dimensionality with the help of subsampling described in Section 3.1.1. Randomized sketching, on the other hand, bypasses the majority of burdensome computation caused by high nonlinearity. The framework of the sketched emulator with local projection is as Fig. 3 where the module of local projection, randomized sketching and the high-dimensional temperature predictor are respectively explained in Sections 3.2.1, 3.2.2, and 3.2.3.

Local projection
Recall the numerical solver with FEM as in Section 2.2 where the temperature at = time instance is denoted as ∈ R on the discrete spatial domain with DoF. Instead of directly computing Eq. (7), we can assume ≈ for a low-dimensional vector ∈ Rã nd project the FEM equations to the subspace spanned by the columns of to yield where the projection basis ∈ R ×̃i s orthonormal with̃≪ . The basis should span a space containing the final temperature result with as less dimensioñas possible. In the reduced GP emulator, we select a subset of training temperatures and weight them according to the relative distance prediction. It indicates that the selected subset forms a basis of the final temperature prediction. Therefore, we develop a similar subsampling scheme to generate local projection bases. Given a test input = [ , , ] , we predict its relative distancef ollowing Eqs. (13)- (15). Then, temperature snapshots in the training temperature matrix corresponding to thẽsmallest values in̂are selected to form a subset̃∈ R ×̃. The local projection basis is the left singular vectors in the compact SVD of̃. We hereafter use boldface to represent matrices and vectors after local projection such as ∶= ∈ R̃×̃and ∶= ∈ R̃. The discretized heat equation after local projection is rewritten as One advantage of local projection is that it does not suffer from the curse of dimensionality as the number of inputs increases and/or each input has a larger range of interest. Both will cause a larger training data set (larger ) yielding more resource cost in training and storage, but once we get the relative distance predictor the dimension after pro-jectioñremains small so the online execution time does not increase. However, the construction of local bases needs to be implemented online since in subsampling we use the values of test inputs which are not known offline, while the construction of global basis with POD method [23] can be finished offline since it constructs the global basis as the left singular vectors of the compact SVD of all temperature snapshots available thus independent from test inputs. In other words, local projection sacrifices the online time cost required to generate local project bases to further reduce the projected dimension from tow ith̃≪ .
Nonetheless, the projected model remains time-consuming due to its high nonlinearity. After local projection, the computation of nonlinear temperature-dependent functions in both and are not reduced. The number of nonlinear computations still depends on the original fine spatial discretization, the order of nonlinearity, and Gaussian quadrature rules. The difficulty of nonlinearity and its solution with randomized sketching are further elaborated in the following section.

Randomized sketching
The nonlinearity of the heat equation comes from temperaturedependent thermal properties and radiation heat loss, all of which are modeled as nonlinear functions. As stated in Appendices B.1 and B.2, the temperature-dependent parts of and , such as the stiffness matrix ∈ R × , the mass matrix ∈ R × , and the nonlinear part of radiation heat loss ∈ R × , need to integrate the corresponding nonlinear functions over the domain discretized into nodes and tetrahedral elements. For example, the mass matrix is defined as where and are respectively the basis function of the th and th node in FEM. The integration in Eq. (20) is approximated by taking integration points in each element based on Gaussian quadrature rules [28,29], which is where ( ) ∈ R × is the matrix of basis functions in FEM evaluated at the th integration point of each element, and ( ) ∈ R × is a diagonal matrix containing the product ( ) ( ) evaluated by the temperature at the th integration point of each element. Further details of integral definitions and Gaussian quadrature approximations are listed in Appendices B.2 and B.3. After local projection, the projected mass matrix ∈ R̃×̃becomes where ( ) = ( ) ∈ R ×̃f or = 1, … , can be pre-computed since they are not temperature-dependent. The diagonals of ( ) for = 1, … , , however, have a total of × values awaiting for nonlinear computation of density and specific heat capacity at integration points in each Picard iteration. It is time-consuming as finer spatial discretization will increase and higher order of nonlinearity will increase . Accordingly, the expensive time cost due to nonlinearity is not reduced by a local projection which is also true for other temperature-dependent matrices like and . It is a major obstacle in implementing the projected model in real-time. To address this issue, we introduce randomized sketching based on Bernoulli sampling. Take the projected mass matrix as an example, the sum in Eq. (22) can be rewritten as the sum of sparse matrices which is where the thin matrix  ∈ R ×̃a nd the diagonal matrix  ∈ R × are respectively We then sketch with only some of these sparse matrices weighted by a factor. Namely, where is a Bernoulli random variable with probability of success 0 ⩽ ⩽ 1. We only manage to dramatically reduce evaluating  without much accuracy compromise when many = 0 and ‖̂− ‖ is small with a high probability. To fulfill this, we set the probability according to the thin matrix  as Algorithm 1 [32]. The row selection scheme is accordingly generated as the row selection vector and the corresponding weight 1∕ for = 1, … , .

Algorithm 1
The row selection and corresponding probability with randomized sketching Input: a thin matrix  ∈ R ×̃, a constant . Output: a row selection vector ∈ R , a probability vector ∈ R . 1: ∈ R ×̃i s the left singular vectors in the compact SVD of  . 2: for = 1 to do 3:

4:
is generated as a Bernoulli random number with the success probability . 5: end for

High-dimensional temperature prediction
Similar to the randomized sketching of the projected mass matrix , we can sketch other temperature-dependent matrices in the locally projected model. The sketched emulator with local projection then becomeŝ̂=̂.
After solving it with Picard iteration, we obtain an accurate result of . The prediction of high-dimensional temperaturêis then reconstructed by the local projection basis followinĝ =̂.

Results and discussion
We develop a numerical experiment which is a printing process on a rectangular parallelepiped AlSi10Mg powder bed evolving from one to three layers in an argon atmosphere. The FEM-based numerical solver with fine spatial-temporal discretization, as a high-fidelity (groundtruth) heat simulation, is employed as a reference model. Its spatial domain is discretized with adaptive meshing as Fig. 4 to maintain the complexity to manageable levels as more layers are incorporated. While the printing area on the top layer keeps fine meshes being tetrahedrons with length 0.01 mm, the remaining meshes are gradually coarsened to tetrahedrons with length 0.1 mm. Under this adaptive mesh scheme, the DoFs in the one to three-layer domain are respectively 1 = 13860, 2 = 13650, and 3 = 14514. The scanning pattern plotted as red lines in Fig. 4 contains three straight lines scanned back and forth. The performance analysis is evaluated by model accuracy and online time cost compared between the reference simulator and the two SMs: the reduced GP emulator denoted as  1 and the sketched emulator with local projection denoted as  2 . All tests are run in Matlab R2020b on an Intel Core i5-9400F CPU at 2.90 GHz, 16 GB RAM computer.

Data generation
The thermal model with FEM following Eqs. (1)- (7) is validated with published experiments in [11] where LPBF processes with AlSi10Mg material in argon atmosphere were experimented and simulated. In [11], the computer-controlled selective laser melting process  Table 1 Model parameters of the data generator [11,12,33]. to construct the training data set of the -layer domain ( ) ∈ R 3×3630 and ( ) ∈ R ×3630 for = 1, 2, 3. The dimensionality of training temperatures , however, is the number of DoF in the -layer domain. Additionally, 9 simulations with different values of laser power (160, 220, 280) W and scan speed (115, 235, 355) mm/s are generated as the reference of tests, which yields 270 pairs of test inputs ( ) ∈ R 3×270 and the corresponding reference temperatures ( ) ∈ R ×270 in the -layer domain with DoF for = 1, 2, 3.

Dimension reduction
In the pre-processing of  1 , the number of GP is reduced by a nonlinear dimension reduction technique. There are two important values in X. Li and N. Polydorides Table 2 The dimension reduction in  1 .

Layer
Reduced dimension Reduction percentage (%) is the singular value diagonal matrix in Eq. (9). With the training temperature snapshots ( ) for = 1, 2, 3, the reduced dimension, as in Table 2, are respectively 91, 152, and 234. In comparison with the number of training data = 3630, it is a dramatic reduction respectively cutting down 97.49%, 95.81%, and 93.55% of GPs required thus less effort in both training and prediction. It is also observed that the percentage of dimension reduction is decreasing as more layers are superimposed. One reasonable explanation is that the temperature and phase information of previous layers are kept, so the information is less redundant in the domain with more layers. Though this decreasing trend exists, it is not significant enough to have a detrimental effect on online time cost.

Prediction accuracy
Given test inputs ( ) in the -layer domain for = 1, 2, 3, we predict the corresponding temperatures aŝ( ) =  1 ( ( ) ). The prediction accuracy hereafter is evaluated by relative errors. For example, the relative error of the th test in the -layer domain is ‖ ( ) * −̂( ) * ‖∕‖ ( ) * ‖. To evaluate the accuracy of all predictions, we graphically plot the spread of the relative errors of temperature profiles as the box plots in Fig. 5 where relative errors of all tests in the -layer domain are plotted as the th box. In a box, the 25th percentile ( 1 ), median ( 2 ), and 75th percentile ( 3 ) are respectively drawn as its bottom, red and top line. Two whiskers are additionally extended with a distance of 1.5 × ( 3 − 1 ) where 3 − 1 is also known as the interquartile range (IQR). Red crosses beyond whiskers represent outliers. We can tell from Fig. 5 that the majority of temperatures can be accurately predicted. Specifically, there are respectively around 17.78% and 99.38% of tests having relative errors below 1% and 5%, while the maximum relative error is around 6.54%.
Despite the relative error of temperature profiles, we also evaluate the prediction accuracy of thermal histories and melt pools. Take the printing process with laser power 280 W and scan speed 235 mm/s as an example, three positions selected from different scanned lines and different layers are used to compare the simulated and emulated results. The coordinates of these 3 positions are specified as  Fig. 6. Thermal histories of ( ) for = 1, 2, 3 throughout the three-layer printing are compared with the reference simulator in Fig. 7, the relative errors of which are respectively 4.49%, 4.51%, and 4.41%. The thermal history of one position starts from the layer it locates to the end of printing. In each layer of printing, three peaks are reached as the laser beam following the scan path of three straight lines approaches. The highest peak, however, appears when the laser beam passes right through the selected position yielding melt pools. It is also where obvious prediction discrepancy occurs while the general trend of curves is accurately approximated. The temperature distribution around melt pools at position ( ) for = 1, 2, 3 are compared as Fig. 8 where the spatial distribution of relative error are also displayed. Higher relative errors concentrate around the areas with higher temperature, in which the maximums of relative errors are respectively 8.82%, 8.04%, and 8.38%. The highest temperature of melt pools at ( ) for = 1, 2, 3 are respectively 2252.0 • C, 2308.2 • C, and 2462.8 • C in the FEM simulator, while with  1 they are respectively emulated as 2390.7 • C, 2446.4 • C, and 2590.3 • C with relative error 6.16%, 5.99%, and 5.18%. Based on the comparison of the temperature colormaps around melt pools, it is indicated that the shape of melt pools is generally emulated. In specific, we compare the size of melt pools as Table 3 where the length, width, and depth of melt pool are respectively the maximum span of the melted area in 1 , 2 , and 3 directions. At positions ( ) for = 1, 2, 3, the melt pool sizes are respectively predicted with relative error below 4.56%, 6.44%, and 3.10% but none of the emulated sizes is predicted with high precision (< 1%).

Time cost reduction
With the level of accuracy compromise described above, a dramatic reduction of online time cost is accomplished since only swift computations are left for online execution. In specific, the comparison of time cost is listed as Table 4. Though the preparation time of data generation and training are required offline, the online execution time is respectively reduced from 3.88 min, 4.21 min, and 4.61 min to 0.024 s, 0.035 s, and 0.049 s on average in the one to three-layer domain. The dramatic reduction of time cost benefits from three advantages of the reduced GP emulator. First, the purely data-driven approach skips the burdensome computation caused by nonlinearity like Picard iteration and Gaussian quadratures. Second, only a comparably small number of GPs are required after nonlinear dimension reduction, which cuts the time cost of offline training and online prediction. Third, the prediction as Eq. (13) manages to simultaneously predict all GPs by executing the dot product between matrices instead of computing the univariate GP individually. While the performance in reducing online execution time is satisfactory, the compromise of accuracy is evident. Therefore,  1 is suitable for applications that ask for real-time implementations but are endurable to some accuracy compromise.

Table 3
The melt pool size comparison between FEM and  1 .

The performance of sketched emulator with local projection
For the performance analysis of the SM  2 , we develop temperature predictions with the same test inputs ( ) for the convenience of comparison. For each pair of test input, specifically ( ) * for = 1, … , 270, a local projection basis is formed by selecting 30 temperature snapshots from ( ) that are predicted as the closest ones to ( ) * . After projecting and sketching, the temperature profiles are predicted aŝ( ) =  2 ( ( ) ) for = 1, 2, 3.

Prediction accuracy
Relative errors between ( ) and̂( ) spread as the boxplot in Fig. 9, from which we can tell that  2 is very accurate. 97.78% of relative errors are less than 1%, and the maximum relative error is 2.09%. We then use the same example (positions (1) , (2) , and (3) in the thermal simulation with laser power = 280 W and scan speed = 235 mm/s) to show the accuracy of  2 in emulating thermal histories and melt pools, which are respectively drawn in Figs. 10 and 11. At positions ( ) for = 1, 2, 3, the thermal histories are respectively predicted with relative errors 0.24%, 0.27%, and 0.25%. Throughout the entire printing process, the temperature curves are finely recovered indicating precise predictions of the peak temperature, heating and cooling rate. The spatial distribution of temperatures around melt pools at ( ) for = 1, 2, 3 are delineated in Fig. 11 where relative errors of temperature prediction are respectively below 1.26%, 0.90%, and 1.12%. In comparison with FEM, the highest temperatures estimated as 2256.2 • C, 2310.0 • C, and 2462.9 • C with  2 respectively have relative errors 0.19%, 0.08%, and 0.006%. Listed in Table 5, with  2 the width and depth of (1) , the width of (2) , and the width of (3) are exactly emulated while the maximum relative error of other melt pool size predictions is 0.82%.

Time cost reduction
Since the accuracy compromise is small, it is expected that  2 needs more running time. While the preparation time of data generation and training is the same as  1 , the total online execution time, listed as 41.50 s, 41.67 s, and 43.51 s on average for the one to three layers domain in Table 6, is composed of three parts: projection, sketching, and prediction. It is noted that in the -layer domain for = 1, 2, 3 the average execution time of projection (2.17 s, 2.66 s, and 2.78 s) and sketching (6.19 s, 6.52 s, and 6.88 s) is much less than prediction (33.14 s, 32.49 s, and 33.85 s) since projection and sketching only needs to be computed once with one set of test inputs while the prediction part includes the Picard iteration taking several runs to converge. As the number of layers increases, the execution time of projection, sketching, and prediction slightly increases but meanwhile the percentage of time cost reduction also increases. It means that compared with FEM the online running time of  2 is less sensitive to more elements in the discretized domain. However, the online execution time is still observed with a significant reduction respectively saving 82.17%, 83.50%, and 84.27% of time cost in the one to three layers domain. The promising advantage of  2 is that it expedites the high-fidelity simulator with a small compromise of accuracy.  2 manages to make most predictions with relative error less than 1% using no more than 18% of online running time required by FEM, thus its performance is encouraging for applications that require accurate thermal simulations but do not excessively demand real-time implementations.

Performance comparison
Two time-efficient SMs: the reduced GP emulator  1 and the sketched emulator with local projection  2 are compared with the same training data set { ( ) , ( ) } and test data set { ( ) , ( ) } in the -layer domain for = 1, 2, 3. Both SMs aim to reduce the online computation time while maintaining accuracy. However, they are, by design, different trade-offs between accuracy and time cost.  1 manages to predict temperatures very fast (0.024 s, 0.035 s, and 0.049 s on average for one to three layers) with relatively higher prediction errors (17.78% and 99.38% of tests have relative errors less than 1% and 5%), while  2 produces very accurate results (97.78% of tests have relative errors less than 1%) in a comparably longer execution time (41.50 s, 41.67 s, and 43.51 s on average for one to three layers). In other words,  1 is biased on time cost reduction while  2 is biased on prediction accuracy. As a data-driven method with nonlinear dimension reduction,  1 highly relies on finding representative training data set to guarantee accuracy. This reliance also exists in  2 since it uses the same relative distance predictor in  1 to construct local projection bases. Nonetheless, the generator of local projection bases has less accuracy requirement in relative distance prediction as all subsampled temperature snapshots equally contribute to local projection bases, while in  1 the exact value and order of relative distance prediction matters since they are used in computing proper weights to extrapolate the high-dimensional temperature predictions. This difference makes  2 more robust in a different sampling of training data.
We can tell from the analysis above that both SMs have good performance in the reduction of time cost and the maintenance of accuracy. The less accurate  1 still suppresses most prediction error less than 5% while the slower  2 still manages to save more than 82% of time cost. We could select  1 or  2 based on the requirement of applications. If only a general description of temperature profiles is enough,  1 is a more effective choice. Wherever precise analysis is necessary,  2 is preferable when the execution time in dozens of seconds is acceptable.

Conclusions
Together with a nonlinear thermal model of LPBF using FEM, two time-efficient SMs: the reduced GP emulator and the sketched emulator with local projection are proposed. To display the performance of two SMs, a numerical experiment is displayed, which is a three-layer printing process scanning three straight lines back and forth with AlSi10Mg material and argon atmosphere. One SM is a reduced GP emulator which first predicts relative distance between temperature profiles with GPR and nonlinear dimension reduction and then extrapolates high-dimensional temperature predictions based on this prediction. It is implemented in real-time (the average execution time is 0.036 s)   . 11. The melt pool comparison between FEM and  2 . while respectively keeping 17.78% and 99.38% of tests with relative errors below 1% and 5%. Another SM, a sketched emulator with local projection, is an expedited FEM solver. It first locally projects the FEM with bases formed by subsampled training data based on the relative distance prediction and then bypasses most nonlinear computation via randomized sketching. It is a very accurate SM (97.78% of tests with relative errors below 1%) while taking a comparably longer time online (the average execution time is 42.23 s). However, it still significantly reduces more than 82% of execution time compared with FEM. Despite the accuracy and time cost of temperature profile predictions, we also compare thermal histories and melt pools at selected positions. Both are finely emulated by the sketched emulator with local projection while being roughly recovered by the reduced GP emulator with visible deviation. The performance of the two SMs indicates that they are different trade-offs between accuracy and time cost: the reduced GP emulator outperforms in online execution time while the sketched emulator with local projection outperforms in model accuracy. In addition, though the data-driven relative distance predictor is used in both SMs, the subsampled training data based on relative distance predictions equally contribute in local bases of the sketched emulator while the exact values of relative distance predictions are used to weight the subsampled training data in the reduced GP emulator. This makes the sketched emulator with local projection less reliant on finding representative training data since it is more robust to the prediction error of relative distances. Overall, it is reasonable to select one SM based on the demand in precision and time cost for different applications. Both provide fast computations of temperature profiles, thermal histories, and melt pool sizes, the values of which are necessary for the models of microstructure and residual stress yielding the prediction of the final part's mechanical properties, deformation, and fatigue life. They, therefore, take one step further in facilitating the quality assurance of final parts. Moreover, though these two SMs are proposed for LPBF, they are essentially applicable for other thermal-driven AM by altering domain shape, heat source model, and so on. The future work will focus on how to properly apply these two SMs in the thermal analysis of LPBF such as defects prediction, uncertainty management, parameter optimization, and closed-loop controls.

Declaration of competing interest
One or more of the authors of this paper have disclosed potential or pertinent conflicts of interest, which may include receipt of payment, either direct or indirect, institutional support, or association with an entity in the biomedical field which may be perceived to have potential conflict of interest with this work. For full disclosure statements refer to https://doi.org/10.1016/j.addma.2022.103122. Xiaohan Li reports financial support was provided by The University of Edinburgh School of Engineering.

Data availability
Data will be made available on request.

Acknowledgments
We are grateful to the Principal's Career Development and Edinburgh Global Research Scholarship from the School of Engineering at the University of Edinburgh, UK for sponsoring this research. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

A.1. Latent heat
The latent heat effect during the phase change between solid and liquid is characterized in the model of specific heat capacity , which is where is the specific latent heat of materials. = 1∕(1 + exp (− ( − ( + )∕2))) is a temperature-dependent sigmoid function, in which > 0 is the logistic growth rate controlling smoothness. and are respectively the lower and upper bound of mushy area [13].

A.3. Polynomial fitting
The temperature-dependent thermal properties including thermal conductivitȳ, density , and specific heat capacitȳare modeled by fitting experimental data. These thermal properties are modeled as piecewise polynomial functions in three temperature ranges: ⩽ , < < , and ⩾ , which respectively indicates the solid phase, the mushy area, and the liquid phase [11]. The fitted polynomials of AlSi10Mg and argon are detailed in Table A.1.

A.4. Thermal properties of powder
For the temperature range ⩽ , the materials have two phases in LPBF: powder and solid. Different from thermal properties of solid which are fitted polynomials based on experimental data, thermal properties of powder denoted as̄, and̄in Eq. (A.3) are modeled according to the porosity of powder bed , the materials, and the inert gas atmosphere [34,35]. where the density and specific heat capacity of the surrounded inert gas is similarly modeled as polynomials fitted by experimental data.

B.4. The algorithm for Picard iteration
As stated in Algorithm B.1, the nonlinear heat equation is solved via Picard iterations which takes several iterations of solving -dimensional linear equations to converge to the solution of the original nonlinear equation with error ‖ − ‖ below a small tolerance, e.g. 10 −5 [14].

Appendix C. The computation of the relative distance matrix
To get the relative distance matrix ∈ R × , we first construct a weighted -neighborhood graph with vertices representing the  The comparison between the thermal model with FEM in this paper and experiment results in [11] where laser power 150 W and scan speed 200 mm/s.

Type of results
Thermal model with FEM Results in [11] Highest temperature ( training temperature profiles and with edges weighted as the adjacency matrix ∈ R × which is given by = { ‖ * − * ‖ if * is a neighbor of * ∞ otherwise, for , = 1, … , , where * is a neighbor of * if their Euclidean distance belongs to the smallest values in ‖ * − * ‖ for = 1, … , . The pair-wise distance matrix is computed as the shortest path matrix of this graph via the Floyd-Warshall algorithm [40].