Descriptor extraction on inherent creep strength of carbon steel by exhaustive search

ABSTRACT The alloying elements that control the creep rupture life of low carbon steel in the region of inherent creep strength were investigated by a data-driven model selection method. The experimental data published in NRIM Creep Data Sheet No. 7B was used. A model in which the Larson–Miller parameter is expressed as a linear sum of the logarithmic stress and the chemical composition of alloying elements was proposed on the basis of the analysis of the experimental data. There were 1023 ( ) candidate models containing one or more of the 10 alloying elements, and they were compared using two model selection methods: 10-fold cross validation (CV) and posterior probability based on Bayesian inference. The 10-fold CV suggested that Mo must be included in the model to improve the predictivity of the unknown data. A comparison among the Bayesian posterior probabilities of the candidate models showed that the Mo-only model has a predominantly high posterior probability of 79.76% compared with the other models. Furthermore, other element terms that co-occur with Mo evenly appeared in the top list after the Mo-only model, indicating that only Mo is crucial for increasing the posterior probability. From these results, it is concluded that the Mo-only model is the most appropriate as a data generation model. The selected Mo-only model was confirmed to have high predictivity of the creep life. Graphical abstract


Introduction
To develop heat-resistant steel for use in power plants and other applications, it is necessary to clarify the effects of chemical composition and process conditions on the long creep life. The chemical composition and process affect the initial microstructure, which influences the high-temperature strength. In addition, the microstructure changes upon exposure to for a long time, and hence the high-temperature strength also changes with time. Due to the involvement of these time-dependent factors, it is generally quite difficult to clarify the effects of chemical composition and process conditions on the long-term creep life. Kimura et al. proposed that there is a region where the creep strength is independent of time for ferritic steels, and they called this idea the inherent creep strength concept [1][2][3]. Specifically, in a class of ferritic steels, they state that precipitated carbides, which are effective in strengthening, become coarse and do not work as a strengthening factor owing to prolonged exposure to high temperatures, resulting in that the creep depends only on the inherent creep strength.
Recently, the number of power plants operated for a long period of time has been increasing, and it is important to clarify the behavior of materials, especially in the long-term side; thus, to clarify the controlling factors of the time-independent creep strength is becoming more important.
In the region where the inherent creep strength is reached, which is unaffected by time-dependent strengthening factors, it is expected that accurate estimates of the creep life can be made. Using low-carbon steel (JIS STB410) as an example, Kimura et al. have identified the range that is inferred to correspond to the inherent creep strength and investigated the factors that affect the creep life [2]. Possible factors include matrix strength, prior austenite grain size, and (thermodynamically stable) dispersion strengthening oxides, but since the prior austenite grain size is nearly constant in the material under analysis and does not include dispersion strengthening oxides, the matrix strength remains as a factor that should be considered [1]. Since the matrix strength is affected by solid solution strengthening, it can be assumed that only the solid solution elements affect the inherent creep strength as a material factor. On the basis of this assumption, the effects of trace elements on the lifetime were investigated, and a good linear relationship was found between the logarithm of the rupture time and the amount of Mo. Additionally, they mentioned that the amounts of Cr and Mn are positively correlated with the logarithm of the rupture time, although not as strongly as in the case of Mo [2].
It is expected that the effects of trace elements in the inherent creep strength region can be clarified for lowcarbon steels. However, it is not straightforward to derive conclusions from the plot of the rupture time vs. the fraction of each trace element. This is because the fractions of the other trace elements concurrently vary with the change of the trace element for which the rupture time is plotted. That is, the trend observed in the plot for a trace element can be merely an apparent one, which would be derived from the change in the fraction of the other elements. It would be possible if one could prepare additional data that allow to make a plot for each element under the condition that the fractions of the other elements are constant. However, this is not a practical solution because of the cost and time required for creep experiments. It is desirable to make full use of the already acquired valuable creep data to properly evaluate the effects of the elements.
In this study, we re-examine the effects of trace elements in the inherent creep strength region for STB410 on the basis of the NRIM Creep Data Sheet No. 7B [4] using data science methods. Specifically, the relationship between the mass fractions of trace elements and creep life is modeled, and a quantitative evaluation of which trace elements should be included in the model is carried out by an exhaustive search method [5]. In this study, two evaluation methods are compared: a model posterior probability based on Bayes' theorem to select an appropriate model for the data generation model and cross-validation to improve the predictivity of unknown data. These two evaluation methods are expected to give the same model selection results when the amount of data is sufficiently large, but they may not necessarily give the same results when the amount of data is limited, such as in the case of creep data. The overall results of these evaluations will be examined to identify the factors that govern the inherent creep strength of the material of interest.

Materials and creep data
The data used are from the NRIM Creep Data Sheet No. 7B [4,6] (hereafter referred to as CDS No. 7B), consisting of nine heats of low-carbon steel STB410, which is similar to Kimura et al.'s study [2]. Their chemical compositions are shown in Table 1. Carbon concentration is specified in STB410 to have less than 0.32 mass%, while it is distributed within a narrow range of 0.20-0.24 mass% in the heats used in this study. For other contained elements, the maximum value and range of Si, Mn, P, and S are specified in JIS, and small amounts of Mo, Cr, Cu, Al, and N are contained, although they are not specified in the standard.
Out of all the 210 data points over a range of 673-773 K and 68.6-372.6 MPa, 80 data points whose stress is 150 MPa or less were considered to have reached inherent creep strength and were included in the analysis. We use the Larson-Miller parameter (LMP) [7] which is the logarithmic rupture time t r ½h� scaled by the test temperature T½K�, as follows: where we provide a constant term C ¼ 20 for all heats as in the previous study [2]. Figure 1 plots LMPs versus logarithmic stresses. The points belonging to the inherent creep strength region are represented by solid marks. In the inherent creep strength region, the LMPs are linearly related to the logarithmic stresses in each heat. Furthermore, the lines appear to be parallel between heats. In detail, the slope of the CAJ seems to be larger than that of the other heats, which may be related to the fact that the CAJ has larger amount of N and smaller amount of P (see Table 1). However, since the difference is small, we here assume that the slope is common among the heats and independent of the composition. Then, it is modeled that the difference between heats, i.e., the effect of compositions, appears on the intercept when the LMP is linearly regressed with logarithmic stress. In other words, this intercept is considered to have an LMP increment owing to solid solution strengthening. Furthermore, we approximate that the contribution from an element is proportional to the mass fraction of the element. Thus, we assume the following model in which the LMP is represented by a linear sum of terms proportional to the log stress and the mass fraction of each element:

Proposed models
where σ and w σ are the test stress [MPa] and the coefficient of logarithmic stress, respectively; the composition of the element k is represented by c k [mass%], and w k is its coefficient; and w 0 is the intercept. Given this model, we use a model selection method to identify which elements should be included as important factors. The number of elements is 10, and in order to consider all combinations containing more than one element, we prepare 1023 (=2 10 À 1) models, examples of which are listed in Table 2. For example, model 1 includes only Mo as an element term, and thus the other elements have no effect on the LMP. Model 2 includes Cr as an element term in addition to Mo, and only these two elements have an effect on the LMP. Model selection was carried out using Bayesian  posterior probability and by the cross-validation method.

Posterior probabilities of models [8]
The posterior probability of a model can be a measure that represents the closeness of this model to the actual model that generates a given data. Here, we describe how to obtain the posterior probability of a model. First, we describe the variables to be used. The data D consists of y and X, where N is the number of observations. The objective variable y is LMP, the explanatory variable X consists of the logarithm of the stress σ, the concentration c of the alloying elements, and the constant term, and w represents the corresponding coefficients. The logarithmic stress and each composition are normalized separately so that the sample mean and the uncorrected sample standard deviation are 0 and 1, respectively. The variables are illustrated below, but elemental indices may be present or absent for each model.
We cannot directly determine the posterior probability PðM i jDÞ of the i-th model M i given the data D, but we can use Bayes' theorem to transform it as follows: Here, the prior probability PðM i Þ of M i is assumed to be uniformly distributed. In other words, since the posterior probability of a model is proportional to the model evidence, it is possible to determine the posterior probability by obtaining the model evidence for all models. The model evidence PðDjM i Þ can be derived as the product of the likelihood function PðDjw i ; M i Þ and the parameter prior Pðw i jM i Þ marginalized with respect to w i , as shown below.
Since each observation is independent, PðDjw i ; M i Þ is assumed to be a normal distribution with an isotropic covariance matrix, � N;i ¼ β À 1 N;i E N;i . Also, assuming unimodality for each parameter, Pðw i jM i Þ is assumed to be a normal distribution with the isotropic covar- The product of the two normal distributions of Equation (5) is transformed as shown below ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi The negative logarithm of the model evidence is called Bayesian free energy FE i [5], which can be derived as The smaller the FE i , the better the pair of the training model and the prior distribution ð PðDjw i ; M i Þ; Pðw i jM i Þ Þ is likely to represent the data. To minimize the FE i , we optimized the pair of β N;i and β d;i by iterative calculations of the following self-consistent equations corresponding to After obtaining the free energy for all models, the posterior probability of the i-th model was determined by the following equation: The posterior distribution of the parameters Pðw i jD; M i Þ can be derived from Bayes' theorem as By taking the logarithm of both sides, one can see that the logarithm of the posterior distribution ln Pðw i jD; M i Þ is a function of w i that consists of the sum of the squared error term (the underlined part of the following equation) and the quadratic regularization term (the double underlined part of the following equation) as follows: This quadratic regularization term plays a role in suppressing overlearning in the estimation of w i . When the posterior probability is given, the predicted distribution of the objective variable y new over the value of an explanatory variable x new is expressed as

Cross-validation
Cross-validation is used to evaluate the ability to predict unknown data for a given model. Here, we performed 10-fold crossvalidation. That is, the data are randomly divided into 10 sets, the k-th test set S k is excluded from the training set to minimize FE i;k , and the resulting μ i;k is used to calculate the mean squared prediction error (MSPE). We compute the 10-fold cross-validation MSPE (CVE) as the sample mean of MSPE k;i with respect to k. The smaller this value is, the better the predictive performance of the model is.
Considering the dependence on sampling, the above procedure was performed 10 times to obtain the mean and standard deviation of CVE i .

Results and discussion
Kimura et al. performed a linear regression of LMP on the logarithm of stress to estimate the rupture time at 773 K-100 MPa, and they found a linear relationship between the logarithm of the estimated rupture time and the amount of Mo, which was shown as a plot in the paper. On the other hand, a similar analysis is stated to have been carried out for other elements, although there were no plots for them [2]. Following the same method of analysis as that of Kimura et al., we investigated the correlation between element concentrations and the estimated rupture times at 773 K-100 MPa for all elements, as shown in Figure 2. We here calculate the Pearson correlation coefficients and show them in the upper part of each figure. The compositions of S and P were not included in the table in the report by Kimura et al. but are reported in the NIMS creep data sheet and are thus included in the present analysis.
As shown in Figure 2(f), there is a linear relationship between logarithmic rupture time and Mo composition. This is in agreement with the report of Kimura et al. The correlation coefficient for it is very high, 0.99. The correlation coefficients of Cr, Al, and Mn are also higher than 0.8, although not as high as that of Mo, indicating clearly that they have a high positive correlation. In the report by Kimura et al., a high correlation was pointed out for Cr and Mn, but there was no mention of the correlation coefficient for Al. Si and Cu were also not mentioned, but are shown in our present study to have relatively high negative correlations. The correlation coefficients of the other elements, C, N, S, and P, were less than 0.5 in absolute value. The correlation coefficients of the elements are summarized in Table 3.
The plots in Figure 2 only allow us to consider the relationship between logarithmic rupture time under a specific condition and one element separately. Furthermore, this plot does not take into account the influences of other elements, so it is not possible to exclude cases where the results happen to be highly correlated owing to the influence of another element. Here, as already mentioned, we introduce a model in which the LMP is expressed as the linear sum of logarithmic stress and concentration of elements, and treat it as a model selection problem. This will allow us to determine the effects of the elements in the entire region of the inherent creep strength. Figure 3 shows plots of FE versus CVE for all models in order to examine the differences in information obtainable from these two methods. For CVE, the uncorrected sample standard deviations obtained from 10-fold cross validation are represented as horizontal error bars. As mentioned earlier, the lower the values of FE and CVE, the better the model can be judged to be. Although the order is different in some places, there is a rough trend that a model with a high CVE also has a high FE. However, a closer look reveals that there is a difference between CVE and FE in their ways of ordering the top group of models. In the case of CVE, many of the models on the left side of the CVE axis are so close to each other that the error bars overlap, making it impossible to determine the order of the models. The top group in the range of CVE < 1:0e þ 4, i.e., the group by the left end of this plot, includes 608 of the 1023 models, and the order among them is unclear. In other words, these numerous models are similar in terms of predictivity for unknown data, so that it is difficult to decide which is better. In contrast, in the case of FE, these models are highly dispersed so that there is a possibility that the order can be clearly determined. This means that a more detailed model selection can be made by using FE, which is an indicator for choosing a model that is closer to the actual data generation model.
We divided this plot into Figure 4(a-j) for each compositional term, where models containing the   labeled compositional term are drawn with black dots. It is noteworthy that all models containing Mo are included in the top group in which they are near the smallest CVE. This observation indicates that the inclusion of Mo is important from the viewpoint of predictivity. In contrast, the models containing other elements are relatively more widely distributed. In particular, for interstitial alloying elements such as C, N, and P, the CVE of the models with and without these elements is spread over almost the entire area, suggesting that they are not essential for creep-life prediction in the inherent creep strength region. All of these elements are known to be important for steel materials because they affect their mechanical properties, but no conclusion could be drawn that they are important for the prediction in the inherent creep strength region, where the strengthening by microstructure has expired.
Considering the fact that the predictivity is higher when Mo is included and that elements other than Mo are not essential from the viewpoint of predictivity, it can be understood that the model with Mo only and the models with Mo and other elements cannot be distinguished from each other from the viewpoint of predictivity; therefore, the order of the CVE was unclear among the top models.
On the other hand, the models including Mo have different FE values, which may enable the ordering of the models. This is understood as follows. FE includes the regularization term derived from the prior distribution and the term that depends on the parameter number d (Eqaution (7)), so that the smaller the parameter number, the more advantageous. The addition of an ineffective element term does not change the likelihood, whereas the increment of the parameter number causes the increase of FE so that it makes the top models dispersive. Figure 5 shows a heat map of the 15 models from the lower FE with μ i , which consists of means of coefficients of the compositional terms. Red indicates positive, blue indicates negative, and blanks represent unused variables. Each column corresponds to each model, and the models are shown from left to right in the descending order of the posterior probability PðM i jDÞ shown at the bottom of the columns.
The fact that the models containing Mo occupy the top position shows that Mo is an important explanatory factor. The first-ranked model with a prominent posterior probability (79.76%) includes only Mo. The secondranked model and lower-ranked ones include elements other than Mo, but there is no consistency for the other elements that co-occur with Mo up to the twelfth position, where all the models with two elements that include Mo appear. For example, although the second-ranked model (14.66%) that consists of Mo and Cr has prominently higher probability than the lower-ranked models (less than 2%), other models containing Mo and Cr begin to appear after the 12th rank, after all the models with Mo and another element have appeared. In addition, it can be pointed out that the absolute values of the coefficients of other elements are smaller than those of Mo. In other words, the predominance of the model group containing Mo and Cr is lower than that of the model group containing Mo, and the results are insufficient to support the indispensability of Cr as an explanatory factor. From these results, it can be concluded that Mo alone is an important explanatory factor for the LMP in the inherent creep strength of carbon steel STB410. Mo was pointed out in a previous study [2] to be the stronger explanatory factor, but this time, we found that the influences of other elements are negligibly small.
Onodera et al. [9] studied the effect of Mo, Mn, and Cr on the creep strength in the inherent creep strength region in terms of the interaction with the interstitial element, i.e., C, and found that there is a linear relation between the logarithmic rupture time at 773 K-88 MPa and the concentration of the Mo-C and Mn-C atomic pairs. On the basis of the analysis, they proposed a mechanism that these atomic pairs reduce the climb velocity of dislocation due to their large interaction energies with dislocations. The proposed mechanism is partly supported by the present results at least for the Mo effect, that is, it is considered that the positive effect of Mo on the LMP can be explained by the effect of the Mo-C atomic pair on the dislocation climb motion. In order to examine the effect of the interaction with C, we additionally carried out the model selection among the candidate models including the cross terms of Mo, Mn, and Cr with C, as shown in Supplementary Fig. S1. We confirmed that the Moonly model is still selected and could not find clear evidence to support that the Mo-C and Mn-C crossterms are necessary. That is, the interaction with C and the effect of Mn were not clearly found. The difference from Onodera's analysis is presumably due to the limitation of the range of C concentration included in the present data, from 0.20 to 0.24 mass%, which is too narrow to find the effect of the interaction with Mo, Mn, and Cr.
The Mo-only model is shown at the bottom of Figure 3, where not only is FE minimal, but also CVE is small. That is, it exhibits a high predictivity for unknown data. On the basis of the above analysis, we propose an equation in which the LMP is expressed as a linear sum of logarithmic stress and Mo concentration as a creep life prediction model for the inherent creep strength region of STB410 carbon steel, as follows: where x new ¼ f 1 log σ c Mo ð Þ À 0 2:03 9:84e À 3 ð Þg (15b) � 1 ð1:01e À 1Þ À 1 ð4:52e À 3Þ À 1 À � ; x new μ ¼ ð2:65e þ 04Þ þ ðÀ 4:61e þ 3Þ log σ þ ð6:97e þ 4Þc Mo ; (15d) � 0 ¼ 1:03e þ 02 À 6:55e À 14 7:95e À 14 À 6:55e À 14 1:05e þ 02 À 1:34e þ 01 7:95e À 14 À 1:34e þ 01 1:05e þ 02 À ð4:11e þ 4Þ log σ þ ð5:00e þ 4Þ : In Equation (15a), the prediction y new for a new condition x new is given in the form of a normal distribution with the expectation x new μ and the variance s 2 , i.e., a predictive distribution. The variables used in Equation (15a) are written in Equation (15b)-(15f), which consist of the values obtained in the aforementioned training with 80 points and the new condition x new . The x new is prepared by scaling the new condition logσ; c Mo to be consistent with the normalization that was applied to the input variables when making the model, as shown in Equation (15b). μ in Equation (15c) is the expectation of the parameter posterior distribution, and the expectation of the prediction of LMP, x new μ, can be calculated using in Equation (15d). Equation (15f) represents the variance s 2 of the predictive distribution, and � 0 in it is the variance of the parameter posterior distribution whose value is represented in Equaiton (15e). The standard deviation s varies depending on x new , but is in the range of 89-97 for the same 80 conditions as those used in the training. When the range of x new μ � s that is one standard deviation away from the expectation is converted into that of creep life for 15 conditions at 723 K-137.2 MPa, it corresponds to a range of 0.74-1.34 times the expectation converted. For example, for CAB whose actual creep life is 10,120 h, the creep life converted from the expectation of LMP is 11,282 h, and the range converted from one standard deviation of LMP is 8549-15,426 h. Figure 2 suggests that other than Mo, strong positive correlations are also found for Cr, Al, and Mn, and negative correlations for Si and Cu. However, assuming the LMP linear model and applying the model selection method in terms of posterior probability based on the Bayesian theorem, the data in CDS No. 7B supported the Mo-only model. Figure 6 shows plots for Cr, Al, Mn, Si, and Cu similar to those in Figure 2, and the interpolation from the experiment is compared with the predictions of the Mo-only model proposed in the study. For any elements, the Mo-only model reproduces the experimental interpolation points well. That is, the observed correlations for these elements are merely apparent and can be explained by the changes of Mo behind the scenes. However, it should be pointed out that this conclusion is only based on the data of CDS No. 7B, and if further data are collected, it is possible that the effects of other elements could be considered essential. Moreover, the scope of application of the conclusions from the results of the present analysis based on the data in CDS No. 7B is presented in a probabilistic manner. For example, if a steel grade is designed with only Mo as a control factor, the creep life of this new steel grade is given in a predictive distribution rather than a single-point prediction. It will be important to examine the width of this predictive distribution for considering the applicability of the prediction.
Besides, the application range of our proposed model should be discussed in terms of Mo contents. Kimura et al. [2] concluded that the effect of Mo addition to increase the inherent creep strength saturates at about 0.03 mass% Mo from the detailed comparison of the LMP vs. logarithmic stress plots among a wide variety of ferritic steels with contents of Mo up to around 1 mass %. Onodera et al. [9] confirmed that the logarithmic rupture time drastically increases with Mo and exhibits a saturation around 0.03 mass%. These previous works provided a critical suggestion that the present Mo-only model is not applicable beyond the 0.03 mass % of Mo. More strictly, the applicable range should be less than 0.019 mass%, which is the highest end of the present data (see Table 1).
Finally, we would like to verify if the predictive distribution of the proposed model is certain. Due to the fact that it is impractical to carry out new creep rupture tests many times, here we compare the observed and predicted distributions in the following manner as a pseudo verification. Equation (15a) is the result of training using all 80 data points, whereas here we exclude the condition that we want to compare with the observed value from the training, i.e., the remaining 79 data points are used to obtain the predictive distribution for the condition. The black dots in Figure 7 show the observed values of LMP and the expectation x new μ of the predictive distributions, while each vertical bar shows an interval corresponding to one standard deviation of the predictive distribution.
The predicted values are in good agreement with the observations. There are 56 points of agreement within the standard deviation, which corresponds to 70% of the total 80 observed points. Considering that the range of one standard deviation corresponds to about 68%, the predictive distribution is a good representation of the actual predictivity. Thus, despite the simplicity of the model that is just a linear model including only logarithmic stress, Mo concentration, and an intercept, it gives excellent prediction performance for conditions with different compositions and stresses, and the prediction errors are properly estimated.

Conclusion
We assumed the framework where LMP can be expressed as a linear sum of logarithmic stress and mass fraction of alloying elements in the inherent creep strength region for nine heats of carbon steel STB410. On the basis of this assumption, we carried out model selection by 10-fold CV and using Bayesian free energy (posterior probability) for possible 2 10 À 1 candidates with one or more of 10 alloying elements, and obtained the following conclusions: • The 10-fold CV revealed that the models containing Mo had better prediction performance than other models. It is concluded that Mo concentration is an essential factor in predicting the unknown data. Each point represents the expectation of predictive distribution, while each vertical error bar represents one standard deviation. Each experimental data corresponding to a condition to predict is excluded from the training set.
• Based on the Bayesian free energy criterion, the Moonly model showed the highest posterior probability of 79.76%, and models containing Mo were at the top of the list. Furthermore, no element had a dominant position as a co-occurring element with Mo. From these results, it is concluded that the Mo-only model is the most appropriate as the data generation model. • The proposed Mo-only model has high creep life predictivity. When the range of 1 standard deviation is converted into creep life, it corresponds to 0.74-1.34 times the expectation. • The creep life at 773 K-100 MPa that was interpolated from the experimental data showed a strong correlation with Mo, and the correlation coefficient was 0.99. There were also strong positive correlations with Cr, Mn, and Al, for which the correlation coefficients were higher than 0.8, and relatively strong negative correlations with Si and Cu. However, it is concluded that the correlations for elements other than Mo are merely apparent, which can be explained by the Mo-only model.