Data-Driven Prediction of Mechanical Properties in Support of Rapid Certification of Additively Manufactured Alloys

Fuyao Yan , Yu-Chin Chan , Abhinav Saboo, Jiten Shah, Gregory B. Olson and Wei Chen

Abstract: Predicting the mechanical properties of additively manufactured parts is often a tedious process, requiring the integration of multiple stand-alone and expensive simulations. Furthermore, as properties are highly location-dependent due to repeated heating and cooling cycles, the properties prediction models must be run for multiple locations before the part-level performance can be analyzed for certification, compounding the computational expense. This work has proposed a rapid prediction framework that replaces the physics-based mechanistic models with Gaussian process metamodels, a type of machine learning model for statistical inference with limited data. The metamodels can predict the varying properties within an entire part in a fraction of the time while providing uncertainty quantification. The framework was demonstrated with the prediction of the tensile yield strength of Ferrium? PH48S maraging stainless steel fabricated by additive manufacturing. Impressive agreement was found between the metamodels and the mechanistic models, and the computation was dramatically decreased from hours of physics-based simulations to less than a second with metamodels. This method can be extended to predict various materials properties in different alloy systems whose processstructure-property-performance interrelationships are linked by mechanistic models. It is powerful for rapidly identifying the spatial properties of a part with compositional and processing parameter variations, and can support part certification by providing a fast interface between materials models and part-level thermal and performance simulations.

Keywords: Additive manufacturing, spatially-varying properties, Gaussian process modeling, statistical sensitivity analysis, maraging stainless steel, yield strength.

1 Introduction

For more than 20 years, additive manufacturing (AM) has been used to fabricate prototypes and models, but only in recent years has it become one of the advanced manufacturing techniques for mission-critical metallic parts. With computer-aided design tools and highprecision machines, the AM techniques, which build the parts layer-by-layer, provide unique capabilities for complex geometries that are difficult to create using other manufacturing processes.

Yet despite recent rapid advances in AM technology, its widespread adoption and commercialization are still hindered, as critical applications, such as those in aerospace,automobile, defense, etc. require a high degree of quality, performance, reliability, and repeatability [Bae, Diggs and Ramachandran (2018)]. The lack of standards and protocols for materials and equipment, as well as low testing efficiency, make complete quality assurance/quality control (QA/QC) one of the biggest constraints of product development for AM [Hessman (2015); Bae, Diggs and Ramachandran (2018)]. Due to compositional and processing variations, the conventional experiments required to certify an industrial product for commercialization are repetitive, cost- and time-consuming. Therefore, virtual certification methods are crucial to reducing the cost and amount of labor. In particular,robust numerical methods that rapidly predict the spatially-varying properties of the parts are highly sought to accelerate the evaluation of repeatability and consistency build-tobuild and machine-to-machine.

Load-bearing parts, for instance, are keys for critical applications. Assuring that the mechanical properties of the AM parts are equivalent or superior to their wrought counterparts is a topmost challenge for the adoption of AM technology. It urges the development of robust predicting tools to feedback the mechanical properties rapidly with respect to different AM thermal histories and post-AM heat treatment steps. Currently,stand-alone simulators that model AM thermal histories and microstructure evolutions only exist in research environments. Most simulators require highly-educated operators as well as expensive licenses and special running environments. More importantly, these simulations are extremely time-consuming, taking minutes to hours to ensure high accuracy, and not practical for certification with massive amounts of calculations. To ensure efficient QA/QC, it is necessary to develop a predicting tool that features: 1)accuracy and speed to perform large quantities of predictions in minutes or less; 2)seamless integration of the models that link the processing-structure-property-performance(PSPP) interrelations of the materials; 3) simple and economical operation that sets low barriers for the skills needed to access the tool.

Metamodels, also known as surrogate models, have grown in popularity as data-driven techniques to replace expensive computer simulations and experiments. Their ability to learn the relationship between sets of inputs and outputs have been well studied [Simpson,Poplinski, Koch et al. (2001); Jin, Chen and Simpson (2001)], and advantageously applied to metamodel-based design analysis and optimization of complex engineering systems[Zhu, Zhang and Chen (2009); Tao, Shintani, Bostanabad et al. (2017)]. Additionally,metamodels have been used to successfully bridge multiple length-scales to design materials systems [Bessaa, Bostanabad, Liu et al. (2017); Bostanabad, Liang, Gao et al.(2018); Hansoge, Huang, Sinko et al. (2018); Xu, Li, Brinson et al. (2016)].

In this work, a multi-response Gaussian process (GP) metamodeling method capable of predicting many locations simultaneously for nonlinear, high-dimensional datasets[Bostanabad, Kearney, Tao et al. (2017)] was chosen to replace the time-consuming microstructure evolution simulations of the AM process and subsequent heat-treating steps.In general, GP metamodels are machine learning models that provide effective statistical inference even when data is limited. Besides modeling highly nonlinear behavior, GP metamodels quantify the epistemic uncertainty of the predictions (due to lack of data or incomplete knowledge of the physics). The Gaussian assumption is a common one for spatial random processes, and used here to capture the prediction uncertainty based on the distance between the prediction site and the collected data, i.e. regions in the model input space with more training data points have less uncertainty. GP models can also be efficiently updated via Bayesian calibration and bias correction techniques when given experimental data. Although this work does not cover uncertainty quantification, robust design or calibration, such features are highly desirable in AM since the complex physics are not yet fully understood and considerable variations exist between builds and machines.In support of AM QA/QC, this work adopts AM Ferrium?PH48S maraging stainless steel and the prediction of its tensile yield strength in the service condition (heat-treated following AM processing) to illustrate the feasibility of GP metamodels as a replacement for complex physics-based mechanistic models and tools. This framework is proposed for AM techniques with a high-power heat source (such as laser and e-beam) such that the materials experience rapid melting and solidification during AM. The PSPP interrelationships are first identified with conventional mechanistic models using the system design approach. Next, a chained method involving intermediate oxide evolution,η precipitation strengthening and solid solution strengthening metamodels is presented to rapidly predict the total tensile yield strength. The metamodels are validated using both cross-validation and an independent dataset. Finally, the prediction accuracy and time consumption of mechanistic and metamodeling methods are compared and discussed.

2 Rapid prediction methodology

Before presenting the PH48S example, a general overview of the techniques used to fit,validate and analyze the GP metamodels is described in this section.

2.1 Optimal design of experiments

Prior to fitting the metamodel, a training dataset must be generated. To efficiently cover the model input space, Design of Experiments (DOE) can be utilized to select the points at which the mechanistic models are evaluated to obtain the desired output(s), or response(s). This work employs the optimal Latin hypercube sampling (OLHS) method, which uses a fast global search algorithm to maximize the uniformity of samples over the input space while minimizing the number of samples necessary [Jin, Chen and Sudjianto (2005)].

For a DOE with n samples, the generated set of d-dimensional inputs is denoted as X=[x1,x2,… ,xn]T, where xi=[xi1,xi2,… ,xid]T,i∈ 1,2,… ,n, and the outputs from the mechanistic models as y=[y1,y2,…,yn]T. The training dataset is then (X,y). In order to build a well-fitted metamodel, it is generally recommended to design a dataset with size n ≥ 10d. Here, n=20d samples are sufficient to ensure the metamodels’ accuracies.

2.2 Gaussian process modeling

In this work, the metamodels are built using the GP modeling method, which can smoothly capture highly nonlinear input-output relationships over a continuous space without needing as much data as other machine learning methods, especially when paired with OLHS. This stochastic technique assumes that the d-dimensional inputs are normally distributed random variables, and defines the GP as a collection of these variables such that any linear combination of them has a multivariate Gaussian distribution. In other words,the model response, y, also has a Gaussian distribution, implying that the mean and confidence interval of the mean prediction, i.e. prediction interval (PI), can be obtained.More specifically, the response can be characterized by its mean and covariance matrix as y(X)~GP(m(X),Σ). The ability to quantify the prediction uncertainty due to noise or lack of data is an advantage of GP models over many metamodeling techniques.

A GP model is formulated as the sum of a linear regression model and a stochastic process:

where fj(X) and βjare basis functions and their unknown coefficients [Rasmussen (2004)].The stochastic part, Z(X), can be taken as the error or “lack of fit” in the model from insufficient data. It has a zero-mean Gaussian distribution with a covariance function c(x,x′)= σ2(R(x,x′)), where σ2is an unknown variance and R(x,x′) is a spatial correlation function (SCF). As is common, the Gaussian SCF is used so that

where θ =[θ1,θ2,… ,θd] are “roughness” hyperparameters.

Therefore, m(X)=fT(X)β and Σij=cij= σ2Rij. Since the metamodels in this work are used for interpolation, not extrapolation, it is reasonable to assume a constant prior mean and set f(X)=1 (and p=1) without danger of reverting to the mean [Staum (2009)].

In order to fit the model, the optimal values of the hyperparameters, β,θ,σ2, can be found using maximum likelihood estimation (MLE), which maximizes the probability, or“likelihood”, of obtaining a prediction that is the same as the observed training response.This can be accomplished by recasting the optimization problem as the minimization of the negative log-likelihood function as follows:

where R is an n×n matrix with elements calculated using Eq. (2). The problem can be simplified by formulating β and σ2as functions of θ so that the minimization is only with respect to the latter.

Recall that in GP models, the variables are assumed to be normally distributed. Thus, once the optimal hyperparameters are fitted, both the mean predictions and their associated uncertainties due to lack of data can be found. The mean prediction at a new point, x0, can be made using the closed-form equation for the Best Linear Unbiased Predictor of Eq. (1),

where r(x0) is a vector of size n×1 such that ri(x0)=R(x0,xi),i=1,…,n.

Additionally, the uncertainty of each prediction point can be easily quantified using the mean squared error (MSE):

where W=1-1n×1R-1r(x0).

Finally, the PI of the mean prediction can be obtained based on

where z is the z-score of the standard normal distribution. For a 95% PI, z=1.96.

The concept of mean predictors and of using PI to quantify uncertainty due to lack of data is shown in Fig. 1 with a toy 1-dimensional example. In Fig. 1(b), the 95% PI is smallnearly coincident with the mean prediction-whereas in Fig. 1(a), which is missing one data point at x=1.2, the PI is much wider.

Figure 1: Gaussian process models for a 1-dimensional toy example with training datasets(a) missing a sample at x=1.2, and (b) including the x=1.2 sample

Details of the mathematics behind GP models and hyperparameter learning methods can be readily found in several sources, including [Martin and Simpson (2005); Rasmussen (2006)].This work uses an enhanced GP modeling technique that can predict multiple responses(q-dimensional y) and is robust for large, high-dimensional and noisy datasets [Bostanabad,Kearney, Tao et al. (2018)]. An adaptive nugget parameter is utilized to smoothly model nonlinear behavior and represent the variance of noise. Like most GP models, the noise variance is assumed to be constant. This technique has been successfully applied to hierarchical woven fiber composites [Bostanabad, Liang, Gao et al. (2018)] and multidisciplinary vehicle suspension design [Tao, Shintani, Bostanabad et al. (2017)].

2.3 GP model validation

To validate the GP metamodels, Leave-One-Out Cross-Validation (LOO-CV) is adopted,where the error metric, here the Mean Absolute Error (MAE), is evaluated by leaving one training point out of the model fitting and repeating until all data points have been left out once. The LOO-CV error for each model is the average of these individual errors. While this would generally require the model to be fitted n times for a training dataset of size n,a closed-form equation for the MAE can be derived for GP models by using the leave-oneout predicted mean [Rasmussen (2006)]:

where ~i indicates that sample point i has been removed.

The definition of the LOO-CV MAE is then

Thus, using Eq. (8), the GP model need only be fitted (i.e. the hyperparameters optimized)once using the complete set of training data to calculate the LOO-CV error.

This work further validates the metamodels by comparing the predictions to those of the conventional mechanistic models. A separate OLHS validation dataset that is independent of the training set and not used for fitting the models is used.

2.4 Sobol statistical sensitivity analysis

With validated metamodels, Sobol statistical sensitivity analysis (SSA), also known as variance-based or global sensitivity analysis, can efficiently be used to study how uncertainties in the inputs contribute to variations in the output, i.e. how strongly each model input affects the response. In this method, the variance of the model output, y=g(x), is decomposed into functions of individual inputs and their interactions, then estimated via Monte Carlo integration, for which hundreds of thousands of model evaluations may be required [Sobol (2001); Saltelli, Annoni, Azzini et al. (2010)]. With fast predictions, the GP metamodels are therefore advantageous for SSA.

The decomposition is as follows:

where d is the number of inputs, Viis the variance when varying xionly, Vijis the variance when varying xiand xj, and so on. Note that the variances with respect to multiple variables can be seen as the variance due to the interaction of those variables.

The Sobol indices, which include the main, interaction and total indices, are found using the decomposed variances. The main sensitivity index (MSI) measures the effect of xialone and is calculated by

The total sensitivity index (TSI), or the effect of xias well as all of its interactions with any combination of the other variables, is

Data-Driven Prediction of Mechanical Properties 349

where ~i means all variables except xi, and E[·] is the expectation.

The sum of all interaction effects of xiwith the other variables can be found simply by the difference between its TSI and MSI. Finally, the input variables with larger indices can be interpreted as having larger influences on the output. The MSI can be utilized for this purpose, though it may be misleading if strong interactions exist.

3 Predictive models of PH48S tensile yield strength

To demonstrate the advantages of the data-driven rapid prediction methods above, the spatially-varying tensile yield strength of PH48S is taken as an example. Based on the PSPP paradigm, physics-based models and a framework to accelerate the prediction are presented in this section.

3.1 The material system

A material can be regarded as a system of aggregates in a hierarchy of characteristic lengths and times. The hierarchical materials structures are closely related to the history of processing and directly influence the properties and performance [Olson (1997); Xiong and Olson (2016)]. To predict the mechanical properties of the materials with high confidence,it is critical to identify the contributing microstructural features that are of different length scales and to figure out the evolution of those features in response to the multi-step thermal processing, i.e. the PSPP interrelationships. These high-fidelity physics-based mechanistic models will lay the foundations for accurate metamodels.

Ferrium?PH48S maraging stainless steel (“PH48S” for short), whose nominal composition is listed in Tab. 1, possesses superior strength by precipitating finely-dispersed rod-shaped η-Ni3(Ti,Al) particles during the aging process. Due to its excellent solidification behavior,it is a desirable candidate for AM processing. (For more information about the alloy, please refer to https://www.questek.com/ferrium-ph48s.html.)

Table 1: Nominal composition of Ferrium? PH48S maraging stainless steel

The PSPP relationships in the AM PH48S system is outlined by the system design chart in Fig. 2. Key processing steps and their process parameters are linked to hierarchical microstructures, which in turn affect the properties that determine the overall performance.This work focuses on the strength subsystem of PH48S.

Figure 2: System design chart for AM PH48S. (Colored lines indicate the key PSPP relationships in PH48S when predicting the strength. Different colors represent different strength contributors that form after specific processing steps.)

AM-fabricated PH48S parts exhibit hierarchical solidification microstructures, such as segregated cellular/dendritic structures, columnar grain structure and fine oxide inclusions[Yan, Xiong, Faierson et al. (2018)], due to cycles of rapid melting/solidification. To achieve the best mechanical performance, a series of heat-treating steps are performed following the AM processing, including high-temperature homogenization and low-temperature aging, as indicated in Fig. 2. Homogenization aims to make the composition uniform throughout the whole part to ensure no heterogeneity in material performance. During homogenization,microsegregation and residual stress are eliminated, and grains become equiaxed through partial recrystallization pinned by sub-micron oxide inclusions (corundum oxide phase in PH48S). The following aging step enables the precipitation of fine η-Ni3(Ti,Al) particles,which can significantly strengthen the material. As indicated by the colored lines in Fig. 2,the strength of the material is primarily determined by the strengthening effects contributed by η-Ni3(Ti,Al) precipitates (blue), solute elements in the matrix (green), primary grain boundaries (pink) and oxide particles (red).

To obtain a good evaluation of the tensile yield strength, it is crucial to predict these key microstructures with regard to the related thermal histories during processing: 1)distribution (diameter and volume fraction) of η-Ni3(Ti,Al) precipitates that form during aging, as they contribute most to the strength through dispersion strengthening; 2)distribution of corundum oxide inclusions that form during the rapid solidification of AM,as they can slightly strengthen the material through dispersion strengthening, and can determine the recrystallized grain size that contributes to the strength by grain boundary strengthening; 3) effect of alloying elements on the matrix strength.

3.2 Physics-based microstructure and property models

A chart that illustrates the flow of predicting the tensile yield strength of PH48S with mechanistic processing-structure and structure-property models is shown in Fig. 3. With alloy composition, oxygen content and cooling rate during AM as the inputs, the oxide distribution and the matrix composition can be first obtained by a precipitation simulation.Oxide particles can directly strengthen the material, and can also determine the recrystallized grain size which then determines the grain boundary strengthening, a contributor to the total strength. With the composition of the homogenized matrix and the aging conditions, the distribution of η-Ni3(Ti,Al) precipitates and the composition of the BCC-Fe matrix can be calculated through an isothermal precipitation simulation, and the corresponding strengthening effects can be predicted. The total tensile yield strength is the superposition of all the strengthening effects. Details of the processing-structure and structure-property mechanistic models are introduced in the following sections.

Figure 3: Flow chart illustrating the prediction of tensile yield strength of PH48S with mechanistic processing-structure and structure-property models (numbers ① to ⑨)

3.2.1 Processing-structure models

The precipitation of oxide particles during the rapid solidification of AM and of the η particles during aging are keys to tensile yield strength prediction. The evolution of precipitates represented by the dependency of particle size and number density on time can be simulated by the precipitation simulator, PrecipiCalc software. PrecipiCalc adopts the numerical LS (Langer-Schwartz)-KW (Kampmann-Wagner) framework in modeling precipitation in the multicomponent system. It builds on CALPHAD-based databases(thermodynamic and mobility) to capture fundamental mechanistic features of multicomponent alloys. The LS model provides an integrated approach to nucleation,growth and coarsening, and KW model provides a numerical solution to the LS model [Jou,Voorhees and Olson (2004); Olson, Jou, Jung et al. (2008)]. Key inputs to the precipitation simulation are a temperature-time profile and alloy composition. This simulation typically takes the most time out of the whole modeling process, since each time step is calculated based on several previously calculated points, and the thermodynamic calculations have to be invoked every few time steps to feedback the equilibrium phase relations. Corundum oxide phase precipitating in the liquid-Fe phase during solidification is simulated by PrcipiCalc with a provided AM cooling rate (model ①). Meanwhile, the evolution of η phase in BCC-Fe phase during aging (isothermal) is also simulated by PrecipiCalc (model④), which not only outputs the distribution of particles but the matrix composition that is used to evaluate the solid solution strengthening by solute elements (Cr, Ni, Co, Mo).

Due to the pinning effect of oxide particles on recrystallized grain boundaries [Yan, Xiong,Faierson et al. (2018)], the recrystallized grain size D can be predicted using the size distribution of oxide particles (mean particle radius R and volume fraction F) through the Zener-Gladman boundary pinning model [Olson (1987)] (model ②), as in Eq. (12),

where the grain size distribution parameter Z is 1.53.

3.2.2 Structure-property models

The prediction of the tensile yield strength of alloys has been a prevailing goal of materials science, and several mechanistic analytical models have been established to describe different strengthening mechanisms with microstructural features.

The primary contributor to total tensile yield strength is η precipitation strengthening, which is achieved by particle-dislocation interactions. When the η particles are fine, dislocations can shear through the particles. Coherency strengthening and modulus strengthening are major strengthening mechanisms, which can be linearly superpositioned as described in Eq.(13) with the η particle distribution (average particle radius r and volume fraction f) simulated by PrecipiCalc,

where the BCC-Fe shear modulus G is 71 GPa, the misfit strain ε is 0.0148, and the Burgers vector b is 0.219 nm.

When the η particles grow larger and become impenetrable, the passing dislocation then loops around the particle, whose strengthening effect can be described by the Orowan equation with rod-to-spherical shape correction, as in Eq. (14),

where

Competing between the dislocation-shearing and dislocation-looping mechanisms, the absolute η strengthening (model ⑥) is given by

Since oxide particles are typically large (～0.1 μm), strengthening by oxide particles is dominated by dislocation looping (model ③), which can be described by the Orowan equation with the PrecipiCalc-simulated oxide distribution (mean particle radius R and volume fraction F) as the input, as in Eq. (17),

where the Taylor factor M is 2.8.

The total precipitation strengthening from η and corundum oxide precipitates (model ⑦)is given by

where the exponent p is 1.71 [Wang, Mulholland, Olson et al. (2013)].

The grain boundary (specifically, martensite packet boundary) strengthening effect (model⑤) can be represented by the Hall-Petch equation [Wang, Mulholland, Olson et al. (2013)]as,

where the Hall-Petch constant for martensite packets K is 0.2 MPa m1/2, and Dpis the martensite packet size, which is linearly related to the recrystallized austenite grain size[Morito, Yoshida, Maki et al. (2006)] as in Eq. (20),

Another strengthening contribution is solid solution strengthening, determined by the solute concentrations in the matrix (model ⑧) which can be expressed by Eq. (21) [Qian (2007)],

where ciis the atomic fraction of the element i in the matrix and kiis the strengthening coefficient listed in Tab. 2. The matrix composition after aging is obtained from the PrecipiCalc simulation of η precipitation.

Table 2: Solid solution strengthening coefficient in BCC-Fe (MPa)

In maraging steels, the dislocation density is ~1014m-2in the as-quenched condition[Galindo-Nava, Rainforth and Rivera-Díaz-del-Castillo (2016); Zhu, Yin and Faulkner(2011)]. Dislocations recover during aging, and the dislocation density is estimated to drop to 0.5×1014m-2[Zhu, Yin and Faulkner (2011)]. Dislocation strengthening is then predicted using the dislocation density ρ by Eq. (22),

Finally, assuming it is the linear superposition of all strengthening contributions, the total yield strength (model ⑨) is given by Eq. (23),

where the Peierls stress of pure BCC-Fe, τα-Fe, is 64 MPa [Wang, Mulholland, Olson et al.(2013)].

3.3 Rapid prediction framework using GP metamodels

Since the physics-based simulations are intractably expensive when performing large amounts of prediction, this framework substitutes mechanistic models at each processing step with GP metamodels, which can then be chained together to predict the end quantity of interest. As outlined in blue in Fig. 4, for the PH48S example, the metamodels replace the oxide and η precipitation simulations, which are the most time-consuming. In total,four metamodels are developed and labelled in Fig. 4: A) oxide precipitation, B) η precipitation strengthening, C) solid solution strengthening, and D) tensile yield strength.The inputs include the alloy composition as well as the critical processing parameters that determine the property-controlling microstructures. By selecting their ranges to reflect composition tolerance and processing variations, the metamodels can be harnessed for robust optimization under uncertainty.

For the oxide precipitation metamodel, the inputs are alloy composition, oxygen content and the AM thermal history, i.e. the cooling rate during solidification. The alloy composition varies spatially in powder AM, but the concentration of each alloying element should be within the chemical tolerance listed in Tab. 1. In conventional manufacturing of steels, oxygen is a trace element with low concentration (~10 ppm). However, the powder feedstock for AM typically contains oxygen with an order of magnitude higher (100~300 ppm), which enables the formation of great amounts of oxide particles in the AM counterparts. Therefore, oxygen content is a key input, ranging within 0.01~0.03 weight %(wt.%). Another factor that influences oxide precipitation is the cooling rate during solidification. Cooling rate differs from location to location within a part and varies between parts when different AM processing parameters are chosen, but it is typically ranged between 103~106K/s.

For the η precipitation and solid solution strengthening metamodels, besides the alloying elements, aging temperature and aging time are chosen as inputs to reflect the thermal condition. They are varied around the standard thermal condition (520oC for 2 h) to represent practical furnace temperature variations and holding time inconsistency. A detailed specification of metamodel inputs and ranges is listed in Tab. 3.

Figure 4: Data-driven framework to replace mechanistic models with GP metamodels

Table 3: Ranges of the inputs used to train metamodels

After selecting the inputs and ranges, the sample points of the training dataset are optimized using OLHS. For each training sample, the physics-based models are used to calculate the oxide distribution, as well as the η precipitation and solid solution strengthening contributions. With these outputs, individual GP metamodels are fitted for oxide precipitation, η precipitation strengthening and solid solution strengthening (models A, B and C in Fig. 4). Using these intermediate metamodels in lieu of the mechanistic models,the oxide distribution and strength contributions for each sample of the training dataset are predicted, and the total tensile yield strength is calculated using Eq. (23). Thus, the contribution metamodels are “chained” together and the final metamodel (model D) is fitted (chained method).

Although it is more straightforward to fit the tensile yield strength metamodel directly from mechanistic predictions (direct method), including intermediate contribution metamodels in the framework allows separate Bayesian calibration or bias correction, if only specific models need to be improved. Calibration against experimental data can also be used to combat any uncertainty in the simulations or analytical equations, since those are deterministic and assumed here to be free from error. Furthermore, detailed SSA can be performed to identify which inputs effect each strengthening contribution rather than the total strength. A comparison of the metamodels built with the same input DOE but different ways of obtaining the output (direct vs. chained methods) is discussed in Section 4.

4 Results and discussions

4.1 Metamodels

A training dataset with ten inputs that combine the inputs to oxide precipitation, η precipitation strengthening and solid solution strengthening (Tab. 3) and 200 samples was generated using OLHS. For each sample, the oxide distribution descriptors, η precipitation strengthening and solid solution strengthening were obtained from PrecipiCalc simulations and analytical models. As mentioned previously, the tensile yield strength was calculated both analytically (direct method) and through the nonlinear superposition of intermediate metamodels (chained method).

While training the metamodels, the LOO-CV MAE for each response was calculated using Eq. (8) and is shown in Tab. 4. Overall, the errors are low; the highest (0.0235) belongs to the mean radius of the oxide precipitation metamodel. The chained tensile yield strength metamodel achieves a better, i.e. lower, MAE than the direct method metamodel.

Table 4: Summary of metamodels and their accuracies

Once the metamodels were fitted, response surfaces were plotted to visualize the relationship between the inputs and outputs (Figs. 5-7). Since only two inputs can be plotted on a single surface, all other inputs were held at the median of their ranges (Tab.3). The mean predictions of the two most influential inputs according to SSA (Figs. 11-15)were plotted for each metamodel response as gradient-colored (blue-green-yellow) surface and 2D contours, with the corresponding colorbar shown to the side. In addition, the upper and lower bounds of the 95% PI were plotted as red surfaces. The response surface for the direct method tensile yield strength metamodel is not shown, as it is nearly identical to that of the chained method metamodel (Fig. 7).

In all plots, the prediction intervals are close to the mean predictions, indicating very small uncertainty and causing the red surfaces of the PI to overlap the gradient-colored surface of the mean. The exception is the oxide metamodel, which shows relatively higher uncertainty at portions of the model for both responses. The accuracies of the metamodels are further studied in the next section.

Figure 5: Response surfaces of the oxide precipitation metamodel for two responses, (a)mean radius and (b) volume fraction

Figure 6: Response surfaces of the (a) η precipitation strengthening and (b) solid solution strengthening metamodels

Figure 7: Response surface of the chained tensile yield strength metamodel

4.2 Validation of metamodels

To further test the metamodels and compare them to the mechanistic models, an OLHS validation dataset was created with a size that is 20% of the training set (40 samples) and consisting of the same ten input variables. These samples were not used to fit the metamodels. The mechanistic outputs y was plotted against the metamodel predictions y^in Figs. 8-10, and the R2, a common measure for the “goodness of fit”, was calculated for each using

where μyis the mean of the mechanistic outputs, y. A value of 1 for R2means a perfect fit.

Figure 8: Comparison of mechanistic models and oxide metamodel for (a) mean radius and (b) volume fraction

Figure 9: Comparison of mechanistic models and (a) η precipitation strengthening and (b)solid solution strengthening metamodels

Figure 10: Comparison of mechanistic models and total tensile yield strength metamodels that were trained using the (a) direct method and (b) chained method

The results show remarkable agreement, with R2values above 0.98 for all metamodels.Notably, the R2of the tensile yield strength metamodel fitted using the chained method is slightly better than when fitted via the direct method (Fig. 10), corroborating with their respective LOO-CV errors. An inspection of Eq. (12) and Eq. (19) sheds light on why the chained method did not magnify the uncertainty of the chained metamodels. While the uncertainty of the sum of independent Gaussian processes is additive, the responses of the oxide precipitation metamodel are not linearly related to the final tensile yield strength.Thus, the uncertainty was not propagated linearly and in fact decreased. This may also indicate that chaining intermediate metamodels captured more of the underlying physics than directly learning the input-output relationship.

4.3 Sobol sensitivity analysis

Using the validated metamodels, SSA was performed using 106Monte Carlo samples on each metamodel, including the intermediate metamodels. The resulting Sobol indices are shown in Figs. 11-15, with the inputs of the metamodels ordered from highest to lowest TSI values, i.e. from most to least influential. In these figures, cooling rate has been shortened to CR in the interest of space. The results for the direct method tensile yield strength metamodel are not shown because they are nearly identical to that of the chained method metamodel (Fig. 15).

For all metamodels, there are only small differences between the MSI and TSI, indicating that there is little interaction between input variables. Also, there is one dominant input variable in each analysis, explaining the almost linear behavior in Figs. 5-7. This may be partly due to the small ranges of the alloy composition variables (Tab. 3), which lead to those inputs having much weaker influence.

Fig. 11 and Fig. 12 imply that the size of oxide particles precipitating during solidification is primarily determined by the cooling rate, while the volume fraction is controlled by the oxygen level in the melt. This means that the oxide distribution is not related to the PH48S composition variations (Tab. 3), but is primarily controlled by processing, such as the cooling rate and oxygen uptake in the initial powder feedstock and the melt pool. Fig. 16 demonstrates the impacts of cooling rate on the size of oxide inclusions by experimentally comparing PH48S fabricated by the laser engineered net shaping (LENS?) AM process to casting. It is obvious that AM PH48S, which experiences a typical cooling rate within 103~106K/sec, exhibits much finer and more dispersed corundum oxide particles than cast PH48S with a cooling rate of only 0.05 K/sec, which agrees well with the SSA results.

Figure 11: SSA results for log10(oxide radius (m)) using the oxide precipitation metamodel

Figure 12: SSA results for oxide volume fraction using the oxide precipitation metamodel

The SSA results for η precipitation strengthening metamodel (Fig. 13) show that aging temperature/time and Ni/Al/Ti contents determine the strengthening effect by the η phase.This is expected as Ni/Al/Ti are major constituents of the η-Ni3(Ti,Al) phase and aging process parameters are key to precipitation evolution. Since aging temperature plays the most important role of all, it is critical to control the aging temperature precisely to achieve uniform strengthening. Similarly, solid solution strengthening is also affected by aging parameters,since the BCC-Fe matrix composition is what remains after the η phase precipitates. As seen in Tab. 2 and Eq. (21), Mo has the highest solid solution strengthening coefficient in BCCFe because its perturbation in concentration affects solid solution strengthening contribution the most, which matches with the SSA results in Fig. 14.

When combining all the strengthening contributions as in Fig. 15, the aging processing parameters, especially aging temperature, are the most significant factors to the total yield strength, since η precipitation strengthening contributes most to the total yield strength.Because cooling rate is another important variable, location-to-location properties within an AM PH48S part are well worth investigating as cooling rate may differ greatly at different locations within the same AM part. However, due to much less strengthening contributed by oxide strengthening and grain boundary strengthening in PH48S compared to η precipitation strengthening, the variations and uncertainty in oxide precipitation are less significant.

Figure 13: SSA results for η precipitation strengthening metamodel

Figure 14: SSA results for solid solution strengthening metamodel

Figure 15: SSA results for total tensile yield strength using the chained metamodel

Figure 16: Oxide inclusions in (a) AM PH48S and (b) cast PH48S with a cooling rate of 0.05 K/sec

4.4 Comparison of prediction time

Since the stand-alone mechanistic models can be integrated and replaced by a single GP metamodel, the prediction of tensile yield strength becomes not only easier and more accessible, but also considerably faster. The time to individually predict the tensile yield strength of the 40 validation samples using the integrated mechanistic models (Fig. 2)versus the chained tensile yield strength GP metamodel (Fig. 4) is compared in Fig. 17. All calculations were performed on the same computer (Lenovo TD350 server with 32 GB memory and two 2.40 GHz Intel? Xeon? processors) to exclude the influences of the computer configurations. Although slight variations in computation time exist within each group, Fig. 17 clearly shows that the average time consumed for one integrated mechanistic simulation (300 s) is dramatically cut down to fraction of a second (0.0025 seconds) by the GP model when predicting one location at a time.

However, a compelling advantage of the GP metamodel is its ability to predict all locations simultaneously in one calculation. With just one run on the same computer, the metamodel obtained the tensile yield strength of all 40 validation points in 0.006 s, whereas it took over 6 h of repetitively invoking the mechanistic models for each point.

Figure 17: Comparison of time to predict the tensile yield strength of one location using physics-based mechanistic models and GP metamodel (chained method)

5 Conclusions

To support the development of virtual QA/QC methods for AM processing, this work has proposed a data-driven rapid prediction framework that replaces location-specific physicsbased mechanistic models with GP metamodels that can predict the spatially-varying properties within an entire part in a fraction of the time. With the system design approach,this work can rapidly identify the PSPP interrelationships within a material system. The prediction of the tensile yield strength of AM PH48S was taken as an example mechanical property to evaluate the feasibility of constructing a seamless and integrated robust prediction tool using GP metamodels. The strength subsystem has been demonstrated to be related with hierarchical microstructures (i.e. oxide particles, η precipitates, grain structure and solutes in the matrix) that form during different processing steps. GP metamodels for oxide precipitation evolution, η precipitation strengthening and solid solution strengthening have been developed, validated and utilized for statistical analysis. The result of the accuracy comparison between metamodels and mechanistic models was impressive,while the computation of tensile yield strength was sped up dramatically from hours of physics-based simulations to less than a second with GP metamodels.

It must be noted that a linear cooling is assumed in this framework in terms of oxide precipitation, while in reality AM materials may experience multiple cycles of heating/cooling. Although a nonlinear thermal profile could be used as a functional metamodel input [Shi and Choi (2011)], it was chosen to simplify the GP modeling by using a scalar value. Nonetheless, this is a valid assumption since the oxide particles precipitate within a high temperature range that has close to a linear profile during cooling due to the high melting point of the oxide ceramic.

In addition, the almost linear behavior and high accuracy of each metamodel can be attributed to the strong dominance of a single variable. This may not be the case for other material systems with stronger, more nonlinear interaction between inputs. However, the nugget-leveraging GP metamodeling method used in this work has been shown to be accurate and efficient even for high dimensions and very nonlinear behavior [Bostanabad,Kearney, Tao et al. (2018)]. In addition, the chaining of intermediate GP metamodels improved the accuracy of tensile yield strength predictions as a result of nonlinear uncertainty propagation. While this again may not be the case for all materials, it is an appealing result, especially as the intermediate metamodels can be individually studied and improved via SSA and Bayesian techniques, respectively.

This methodology is suitable for the rapid identification of spatially-variant properties within a part with compositional and processing parameter variations. The speed and builtin uncertainty quantification of the GP metamodels is advantageous for guiding the sampling of new points and optimization under uncertainty, as well as efficient Bayesian calibration and bias correction against experimental data. Furthermore, this method can predict properties at all part locations simultaneously in one calculation within seconds,allowing a fast interface between part-level thermal and performance simulations. Such features are invaluable for AM, a technology that is becoming ever more prevalent but in critical need of robust and rapid certification methods.

Acknowledgment:This work was supported by the Digital Manufacturing and Design Innovation Institute (DMDII) through award number 15-07-07. This material is also based upon the work of Ms. Yu-Chin Chan supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1842165. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.