1 Introduction

Most of time, the creep fracture of weld joints of ferritic heat-resistant steel is caused by Type IV cracking after creep damage accumulation in the fine-grain heat-affected zone (FGHAZ, hereafter referred to as HAZ) during fusion welding [1,2,3]. The creep strength of HAZ is low, and the creep rupture life of a simulated HAZ material, in which the HAZ microstructure is simulated by heat treatment that reproduces the thermal history during welding, can be as short as 1/10 of that of the base metal [1, 4]. For this reason, the creep rupture life of an actual weld joint is between the creep rupture lives of the base metal and the simulated HAZ material [5, 6].

We have investigated a framework for accurately predicting the creep rupture life of steel weld joints with relatively simple HAZ geometries based on the creep properties of the base metal and the simulated HAZ material within the framework of the time exhaustion rule [7, 8]. As a result, it was found that the creep rupture life of a weld joint can be predicted well by evaluating the creep damage accumulation on the basis of the time exhaustion rule using the Huddleston stress. Detailed analysis revealed that the rate at which creep damage accumulates in the weld joint depends on the degree to which the HAZ creep deformation is constrained by the base and weld metals [9]. The closer the interface between HAZ and the base or weld metal is, the greater the deformation restraint is, and the more the interface is perpendicular to the loading axis, the more the damage accumulation is suppressed. Conversely, in the HAZ region where the deformation constraint is small, damage tends to accumulate as the equivalent stress increases. Thus, the deformation constraint changes depending on the cross-sectional shape of HAZ, which eventually changes the creep rupture life. In other words, the creep performance of a weld joint can be improved if HAZ has a shape that can restrain deformation effectively. From this viewpoint, we came up with an idea to use the cross-sectional shape of HAZ as a clue for tuning welding conditions that will result in a long creep rupture life of the weld joint.

Table 1 shows the previous studies on the creep rupture of heat-resistant ferritic Cr steel weld joints featuring Type IV cracking, focusing on the HAZ geometry. The previous studies can be mainly classified into two on the basis of the approaches used. One is to investigate the effect of the HAZ geometry by experiment and calculation for simple HAZ geometries, and the other is to perform finite element analysis by simulating the complex HAZ geometries that appear in actual weld joints. In the former approach, Francis et al. [10] reported that the HAZ angle has a significant effect on the creep rupture life, on the basis of the results of creep tests conducted on specimens with different HAZ angles obtained by changing the sampling angle from the weld joint. The effect of the HAZ angle was examined by Parker [11], Koiwa et al. [12], and Tanner et al. [13] by finite element analysis, the result of which supports the trend that the more inclined HAZ is from the load axis, the shorter the creep rupture life becomes. On the other hand, Divya et al. [14] studied the effect of the HAZ width by comparing the result of experiments and finite element analysis, and they reported that the effect of the HAZ width is more pronounced at high stresses. In the latter approach, Watanabe et al. [15] developed a finite element model that precisely copied the HAZ geometry in an actual weld joint creep specimen, performed finite element analysis, and compared the results with real creep test results. Hongo et al. [16] similarly performed finite element analysis to simulate experiments and compared the results with microstructural observations from a creep interruption test, noting that the areas where the creep strain and stress triaxiality are concentrated accumulate much creep damage. In addition, Kawashima et al. [17] studied changing material properties in the vicinity of HAZ in a stepwise manner.

Table 1 List of papers investigating creep property of steel welded joints focusing on HAZ shape 
Full size table

Thus, it can be seen that the creep rupture life caused by Type IV cracking varies depending on the HAZ geometry. For simple geometries of HAZ, it is suggested that the rupture life can be described by a relatively prospective relation depending on the HAZ angle width. For general complex HAZ geometries, similarly, it is possible to describe the rupture life by an equation, e.g., a linear regression model, which clearly shows the relationships between HAZ shape and the rupture life if the geometric features (HAZ shape factors) are appropriately designed. We have already established a computational workflow that simulates multilayer welding and estimates the rupture life for gas tungsten arc welding (GTAW), as described in the next section and elsewhere [18]. Using this workflow as a solver of the forward problem that calculates a rupture time from a welding condition, we can consider an inverse problem to search for welding conditions that increase the rupture life in combination with a sequential optimization method.

In this study, we construct a tandem type Bayesian inference-based model which predicts the rupture life from the welding conditions by utilizing the data generated by the computational workflow; here, the rupture life is given by a linear regression from the HAZ shape factors. Then, an optimization method based on this model is developed. The target material is 2 1/4Cr–1Mo steel, which is heat-resistant ferritic Cr steel, and creep tests of a weld joint with three-layer welded plate I-bevel were simulated by GTAW.

2 Proposed concept

2.1 Development of a tandem Bayesian model

We have developed the computational workflow to determine the HAZ region by heat transfer analysis for any welding conditions and then calculate the creep damage of weld joints including HAZ. This is reported separately [18]. The workflow is implemented on the MInt system [19, 20] developed in the SIP project [20,21,22,23,24,25,26,27,28,29,30,31,32] and can be executed via an application programming interface, and the calculation results can be retrieved. In this study, we substituted each step of this workflow using two Bayesian statistical models via the HAZ geometry and constructed a prediction model (tandem Bayesian model) that links the two models in a Bayesian manner.

The schematic structure of the workflow is shown in the lower part of Fig. 1. The workflow specifies the HAZ shape in the steady-state cross-section of the weld on the basis of the input welding conditions by a two-dimensional heat transfer analysis that simulates the cross-section of the specimen perpendicular to the weld line. Next, a 2.5-dimensional finite element model is constructed by extruding this shape in the depth (weld line) direction, and the creep rupture life is calculated by creep damage calculation. The workflow is divided into two parts with the HAZ shape as the connecting point. Let \(wf1\) be the workflow that returns the HAZ shape from the welding conditions and \(wf2\) be the workflow that returns the rupture life from the HAZ shape.

Fig. 1
figure 1

Schematic diagram of a tandem Bayesian model in which the forward computation method is transferred to the Bayesian model 

Full size image

We used the data calculated from these workflows as training data and constructed two Bayesian statistical models as surrogate models for these workflows, as shown in the upper part of Fig. 1. First, we formulated an appropriate reduced representation of the HAZ shape that was expected to represent the creep rupture life, i.e., HAZ shape factor. Then, we constructed two types of Bayesian statistical models to output the HAZ shape factors from the welding conditions and a Bayesian statistical model to predict the rupture life from the HAZ shape factors. Model selection based on Bayesian statistics was carried out for the latter model to identify the HAZ shape factors with high explanatory capability for rupture life.

2.2 Predicted distribution for tandem Bayesian models

In this section, we describe the Bayesian statistical models used in the tandem Bayesian model and the method used to calculate the prediction distribution. The first half of the Bayesian statistical model, which predicts each of the HAZ shape factors from the welding conditions, employs a Gaussian process (\(GP\)) regression owing to the highly nonlinear nature of the welding process. For the latter, the Bayesian linear regression model (\(BL\)) was adopted to easily understand the relationship between the two, as described in the “Introduction” section. Moreover, HAZ shape factors with high explanatory capability for rupture life were selected using the model selection method based on the Bayesian free energy [33, 34], as detailed below.

Figure 2 shows the graphical representation of causality for the tandem Bayesian model. The set of welding conditions is \(x\), the set of HAZ shape factors is \({y}_{{\text{shape}}}\), and the rupture life is \({y}_{Tr}\). The encircled variables are random variables. Figure 2a, b, c shows the causality representations corresponding to \(GP\) and \(BL\) and the tandem Bayesian model combining \(GP\) and \(BL\), respectively. Here, \(GP\) and \(BL\) are the probability models for \({y}_{{\text{shape}}}\) and \({y}_{Tr}\), respectively.

Fig. 2
figure 2

Graphical representation of causality for a surrogate model 1, b surrogate model 2, and c the present tandem Bayesian model

Full size image

The simultaneous distribution \(p\left({y}_{Tr},{y}_{{\text{shape}}},BL,GP;x\right)\) of the tandem Bayesian model can be written as 

$$p\left({y}_{Tr},{y}_{{\text{shape}}},BL,GP;x\right)=p\left({y}_{Tr}|{BL,y}_{{\text{shape}}}\right)p\left({y}_{{\text{shape}}}|GP;x\right)p\left(BL\right)p\left(GP\right)$$
(1)

on the basis of Fig. 2c. Using the training data \(\{x,{y}_{{\text{shape}}},{y}_{Tr}\}\), let us consider the predicted HAZ shape factors \({y}_{{\text{shape}}}^{{\text{new}}}\) and predicted creep rupture life \({y}_{Tr}^{{\text{new}}}\) for a new welding condition \({x}^{{\text{new}}}\). The predicted distribution of \({y}_{{\text{shape}}}^{{\text{new}}}\) is the prior distribution of \({y}_{Tr}^{{\text{new}}}\). The predictive distribution \(p\left({y}_{Tr}^{{\text{new}}}|{y}_{Tr},{y}_{{\text{shape}}} ;x,{x}^{{\text{new}}}\right)\) of \({y}_{Tr}^{{\text{new}}}\) is given by the marginalization of the joint distribution \(p({y}_{Tr}^{{\text{new}}},{y}_{{\text{shape}}}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}},GP,BL;x,{x}^{{\text{new}}})\) by \(GP\), \(BL\), \({y}_{{\text{shape}}}^{{\text{new}}}\) and \({y}_{Tr}^{{\text{new}}}\) as

$$\begin{aligned} p\left({y}_{Tr}^{{\text{new}}}|{y}_{Tr},{y}_{{\text{shape}}} ;x,{x}^{{\text{new}}}\right)=\frac{p({y}_{Tr}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}};x,{x}^{{\text{new}}})}{p({y}_{Tr},{y}_{{\text{shape}}};x,{x}^{{\text{new}}})}\\ \propto p\left({y}_{Tr}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}};x,{x}^{{\text{new}}}\right)\\ =\int {{\text{dy}}}_{{\text{shape}}}^{{\text{new}}} {\text{d}}BL\, {\text{d}}GP\,\boldsymbol{ }p\left({y}_{Tr}^{{\text{new}}},{y}_{{\text{shape}}}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}},GP,BL;x,{x}^{{\text{new}}}\right)\end{aligned}$$
(2)

(\(p\left({y}_{Tr},{y}_{{\text{shape}}};x,{x}^{{\text{new}}}\right)\) appearing in the denominator is a constant for a fixed training data \(\{x,{y}_{{\text{shape}}},{y}_{Tr}\}\) and is able to be omitted hereafter.) Following the graphical representation of causality in Fig. 2c, we take apart Eq. (2) into terms and partition out to integrals by each integral variable. The contents of each integral can be transformed as follows:

$$p\left({y}_{Tr}^{{\text{new}}},{y}_{{\text{shape}}}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}},BL,GP;x,{x}^{{\text{new}}}\right)=p\left({y}_{Tr}^{{\text{new}}}|BL,{y}_{{\text{shape}}}^{{\text{new}}}\right)p\left({y}_{Tr}|BL,{y}_{{\text{shape}}}\right)p\left(BL\right)p\left({y}_{{\text{shape}}}^{{\text{new}}}|GP\right)p\left({y}_{{\text{shape}}}|GP\right)p\left(GP\right)$$
(3)

(Henceforth, the definite variables \(x\) and \({x}^{{\text{new}}}\) are omitted.) Thus, we finally obtain the following:

$$\begin{aligned} p\left({y}_{Tr}^{{\text{new}}}|{y}_{{\text{shape}}},{y}_{Tr};x,{x}_{{\text{new}}}\right)\propto \int {\text{d}}{y}_{{\text{shape}}}^{{\text{new}}} {\text{d}}GP\, {\text{d}}BL\, p\left({y}_{Tr}^{{\text{new}}},{y}_{{\text{shape}}}^{{\text{new}}},{y}_{Tr},{y}_{{\text{shape}}},BL,GP;x,{x}^{{\text{new}}}\right)\\=\int{{\text{d}}y}_{{\text{shape}}}^{{\text{new}}}\left[\int{\text{d}}GP\,p\left({y}_{{\text{shape}}}^{{\text{new}}}|GP\right)p\left({y}_{{\text{shape}}}|GP\right)p\left(GP\right)\times \int {\text{d}}BL\,p\left({y}_{Tr}^{{\text{new}}}|BL,{y}_{{\text{shape}}}^{{\text{new}}}\right)p\left({y}_{Tr}|BL,{y}_{{\text{shape}}}\right)p\left(BL\right)\right]\end{aligned}$$
(4)

The marginalization for \(GP\) is the predictive distribution for a Gaussian process, and the marginalization for \(BL\) is the predictive distribution for a Bayesian linear regression model. The remaining marginalization for \({y}_{{\text{shape}}}^{{\text{new}}}\) was obtained by averaging over the samples by sampling for a Gaussian process described below.

2.3 Details of the Bayesian statistical model construction

\(GP\) Was constructed using GaussianProcessRegressor from scikit-learn [35] and a Gaussian process with a RationalQuadratic kernel. The training data were normalized so that both explanatory and objective variables had a mean of 0 and a standard deviation of 1, and hyperparameters were optimized for each time when training data were added. Hyperparameters were lengthscale, alpha, and mix. A grid search for these hyperparameters was conducted within the grid shown in Table 2, and the maximum value of the pseudo-log likelihood [36] defined by the CVE of tenfold CV for the training data was adopted. Considering the possibility that the optimal value of lengthscale may vary greatly depending on the data range, coarse tuning was performed according to the initial values shown in Table 2. If the maximum value of the pseudo-log likelihood was at the edge of the search range, the grid was shifted in parallel so that the maximum value was at the center of the grid, and when the maximum value was within the search range, fine-tuning was performed around the maximum value so that the grid precision was set to two decimal places.

Table 2 Grids range of hyperparameters of Gaussian processes
Full size table

\(BL\) was constructed using a Bayesian linear regression model by the model selection method [34]. As described below, seven HAZ shape factors were formulated and used as candidate explanatory variables. The Bayesian free energy was calculated from the negative logarithm of the marginal likelihood for the creep rupture life, and the Bayesian linear regression model with the lowest Bayesian free energy was selected. The training data were normalized so that both explanatory and objective variables had a mean of 0 and a standard deviation of 1, as was the case with \(GP\). As a result, the mean of the objective variable was zero, so the linear regression model was constructed with no intercept.

2.4 Proposed optimization algorithm

We proposed an efficient optimization algorithm to maximize the rupture life as follows. First, the probability \(p\left({y}_{Tr}^{{\text{new}}}>T{r}_{{\text{target}}}|{y}_{{\text{shape}}}^{{\text{new}}};{x}^{{\text{new}}}\right)\) that the predicted rupture life \({y}_{Tr}^{{\text{new}}}\) exceeds a target rupture life \(T{r}_{{\text{target}}}\) for a given welding condition \({x}^{{\text{new}}}\) was defined as a figure of merit \((FoM)\), and candidates with high \(FoM\) were selected. These candidates were calculated in the workflow, the results were added as new training data, and the model was updated. The \(FoM\) was then recalculated according to the updated model, and the next candidate was selected. In this way, optimization was performed sequentially. The distinctive feature of the present optimization algorithm is that it took advantage of the fact that the prediction model was a tandem one. Moreover, the algorithm consisted of double loops: the first loop updating only the workflow \(wf1\) and the second loop updating the workflow \(wf2\). After the first loop was repeated as long as new candidates continued to be selected, the second loop was activated. This double-loop algorithm aimed to update the surrogate model for \(wf1\) efficiently, which is computationally relatively less expensive and requires more training data owing to its strong nonlinearity to surrogates.

The pseudocode of the optimization algorithm proposed in this study is shown in Fig. 3. The training data for \(GP\) is \({D}_{{\text{train}}}^{GP}=\{{x}^{GP},{y}_{{\text{shape}}}^{GP}\}\), consisting of the welding condition \({x}^{GP}\) and the HAZ shape factors \({y}_{{\text{shape}}}\). The training data for \(BL\) is \({D}_{{\text{train}}}^{BL}=\{{y}_{{\text{shape}}}^{BL},{y}_{Tr}\}\), consisting of the HAZ shape factors \({y}_{{\text{shape}}}^{BL}\) and the rupture life \({y}_{Tr}\). The space of welding conditions to be explored is denoted as \({x}_{{\text{all}}}\). The welding conditions for the initial training data are extracted from \({x}_{{\text{all}}}\). HAZ shape factors and creep rupture lives are obtained to generate the initial training data \({D}_{{\text{train}}}^{GP}\) and \({D}_{{\text{train}}}^{BL}\) from the welding conditions, using \(wf1\) and \(wf2\). The validity of the calculation results is verified by analyzing the damage initiation position at which rupture occurs, as described below. The maximum creep rupture life in \({D}_{{\text{train}}}^{BL}\) is set as \(T{r}_{{\text{target}}}\).

Fig. 3
figure 3

Pseudocode of the search method of this research, who consists of the optimizing double loop of step and round. The Bayesian linear regression model BL is optimized at the outer loop(round) and the Gaussian process GP is optimized at the inner loop(step). In this study, n is set as 10

Full size image

The loop that updates \(GP\) by repeatedly adding data to \({D}_{{\text{train}}}^{GP}\) is called the step, and the loop that updates \(BL\) by adding data to \({D}_{{\text{train}}}^{BL}\) is called the round. First, the \(FoM\) is calculated for \({x}_{{\text{all}}}\) from \(BL\) and \(GP\) obtained in the first round. The welding conditions \({x}_{to{p}_{n}}\) from the one with the largest \(FoM\) to the \(n\) th are selected as candidates, and the HAZ shape factors are obtained by driving \(wf1\) for them and added to \({D}_{{\text{train}}}^{GP}\), and \(GP\) is updated. Then, the \(FoM\) is calculated again, and a candidate \({x}_{to{p}_{n}}\) is selected similarly. This is a step, and the steps are repeated as long as a new candidate that is not included in \({D}_{{\text{train}}}^{GP}\) is selected as \({x}_{to{p}_{n}}\). When no more data to be added appear, the loop of a step exits and \({D}_{{\text{train}}}^{BL}\) is rebuilt for the next round, where \(wf2\) is driven to find rupture lives for the welding conditions in \({x}_{to{p}_{n}}\) that are not included in \({D}_{{\text{train}}}^{BL}\) and rupture life data are added to \({D}_{{\text{train}}}^{BL}\). \(T{r}_{{\text{target}}}\) is also updated to the new value if it actually updates the previous maximum creep rupture life. In the next round, \(BL\) is updated with the new \({D}_{{\text{train}}}^{BL}\), and steps are repeated with the updated \(BL\) and \(T{r}_{{\text{target}}}\). When no new candidates are proposed and no more existing data are added to \({D}_{{\text{train}}}^{BL}\) in the course of repeating the above steps and rounds, the calculation is completed with the \(T{r}_{{\text{target}}}\) at that time as the longest rupture life \(T{r}_{{\text{best}}}\) and the welding condition that gives \(T{r}_{{\text{best}}}\) as the best welding condition \({x}_{{\text{best}}}\). The training data \({D}_{{\text{train}}}^{GP}\) and \({D}_{{\text{train}}}^{BL}\) at the end of the calculation become \({D}_{{\text{final}}}^{GP}\) and \({D}_{{\text{final}}}^{BL}\), respectively.

Sampling for \({y}_{{\text{shape}}}^{{\text{new}}}\) is carried out to calculate the \(FoM\). \(s\) sample is obtained for the HAZ shape factors \({y}_{{\text{shape}}}^{{\text{new}},(i)}\) \((i=1,\dots ,s)\) predicted by \(GP\) from \({x}^{{\text{new}}}\) with the sampling function of GaussianProcessRegressor. The predicted distribution of a predicted value \({y}_{Tr}^{{\text{new}},(i)}\) from a sample \({y}_{{\text{shape}}}^{{\text{new}},(i)}\) is obtained by solving Eq. (4). When this distribution is written as \(\mathcal{N}({y}_{Tr}^{{\text{new}},\left(i\right)};\mu \left({y}_{{\text{shape}}}^{{\text{new}},\left(i\right)}\right),\Sigma \left({y}_{{\text{shape}}}^{{\text{new}},\left(i\right)}\right))\), the \(FoM\) is defined as the sample mean of the values integrated from \({Tr}_{{\text{target}}}\) to \(\infty\), as follows:

$$FoM\left({x}^{{\text{new}}}\right)=\frac{1}{n}\sum_{i=1}^{s}{\int }_{T{r}_{{\text{target}}}}^{\infty }{{\text{dy}}}_{Tr}^{{\text{new}},\left(i\right)} \mathcal{N}({y}_{Tr}^{{\text{new}},\left(i\right)};\mu \left({y}_{{\text{shape}}}^{{\text{new}},\left(i\right)}\right),\Sigma \left({y}_{{\text{shape}}}^{{\text{new}},\left(i\right)}\right))$$
(5)

In this study, we set the number of samples \(s=100\).

3 Computation procedures

3.1 Weld calculation

The groove geometry of the weld, the Goldak-type heat source model, and the finite element model outline for the rupture life calculation assumed in this study are the same as those in another report [18]. A single-pass three-layer weld was simulated for an I-shaped groove with a plate thickness of 8 mm and a groove width of 4 mm, and welding simulation was performed using Code-Aster [37] to obtain a two-dimensional distribution of the maximum temperature reached during welding in the plane perpendicular to the weld line via heat transfer analysis. The HAZ region was defined as the range where the maximum temperature attained during welding is between \(Ac1=775^\circ{\rm C}\) and \(Ac3=875^\circ{\rm C}\), on the basis of these results and with the observed result of temperature tests during welding and post-weld microstructure observation.

3.2 Creep damage calculation

The two-dimensional cross-sectional HAZ shape was extruded in the direction of the weld line to generate a 2.5-dimensional finite element model assuming a 4-mm-deep quadrangular prismatic specimen, and the finite element model was constructed as a 1/4 model with symmetry planes set on the \(x\)- and \(z\)-axes. HAZ was assigned the material properties of the simulated HAZ, whereas the base and weld metals were both assigned the material properties of the base metal. The parameters for material and creep properties are the same as those in a previous report [8].

The finite element solver for creep damage calculations is FrontISTR [https://www.frontistr.com/index.php], which incorporates damage calculation functions using the time exhaustion rule framework. The Norton law [38] was used as the creep constitutive equation. The Huddleston stress [9] was used as the scalar stress in the time exhaustion rule, and the rupture life was defined as the time when the damage of a certain element in the finite element model reaches 1. Calculation time steps were set to 1 h up to 100 h and every 10 h thereafter. The creep test conditions simulated in the calculation were a test temperature of 873 K and a tensile stress of 100 MPa.

3.3 Search space, initial training data space

The range of search space \({x}_{{\text{all}}}\) is shown in Table 3. The welding parameters for the search are layer thickness (mm), heat input (J/mm), and virtual heat source width (mm) for each of the three welding layers. The dimension of the search space is 8, excluding the third layer thickness (\(=8 mm-(\mathrm{first layer thickness} + \mathrm{second layer thickness})\)). The 78 = 5,764,801 pairs of grid points were used as the search space \({x}_{{\text{all}}}\). The L49 orthogonal table of design of experiments was used to extract 49 conditions so that all levels were uniformly selected. For these conditions, \(wf1\) and \(wf2\) were driven, and the initial training data were extracted. It was confirmed that 27 of the results for which \(wf2\) was driven resulted in the free-edge fracture [18], where the fracture initiation element was located at the interface between the HAZ and the weld metal or base metal on the HAZ free-edge surface (free-edge fracture), we judged that the rupture life had not converged. Many of the free-edge fractures had short rupture lives, so we decided to exclude them from the initial training data of \({D}_{{\text{train}}}^{BL}\) because they do not contribute to the purpose of this study. As a result, the number of initial training data for \({D}_{{\text{train}}}^{BL}\) was 22. Since it is not possible to automatically determine whether an additional training data is a free-edge fracture or not, we decided not to apply such exclusion to the additional training data.

Table 3 Range of search welding conditions
Full size table

3.4 Computational environment

The Python 3.6 program controlled by bash scripts built on Ubuntu 18.04 was used to drive the tandem Bayesian model. Bayesian free energy and \(FoM\) calculations were performed using functions in NumPy [39]. Forward calculation for \(wf1\) and \(wf2\) were conducted by the MInt system [19, 20] via an application programming interface.

4 Results

4.1 Initial training data and proposed HAZ shape factors

Table 4 shows the welding conditions, rupture life, and HAZ shape factors for the initial training data. The 22 conditions in bold are the data for \({D}_{{\text{train}}}^{BL}\). The longest and shortest rupture lives within \({D}_{{\text{train}}}^{BL}\) are 580 h and 290 h, respectively.

Table 4 Initial training data, obtained from L49 diagonal table of experimental design method
Full size table

Figure 4 shows the cross-sectional HAZ shape included in \({D}_{{\text{train}}}^{BL}\), arranged in order of rupture life. The HAZ shape with a longer creep rupture life is smaller in width or the area is small in width, and is generally upright perpendicular to the tensile direction. The HAZ shape with a short creep rupture life is generally wide and the entire HAZ is inclined.

Fig. 4
figure 4

Cross-sectional shapes of the HAZ obtained by heat conduction calculation for 22 welding conditions in which it was determined that the rupture life is calculated properly at the creep damage calculation of the 49 initial data, by the rupture initiation position analysis. The calculation was performed at the testing condition of 873 K, 100 MPa

Full size image

From these observations, the seven HAZ shape factors are shown in Fig. 5 as they were developed as expected to govern the creep rupture life. To mathematically define the HAZ shape factors, the following variables are introduced. The HAZ width \({w}_{i}\) is defined as the difference between the loading direction coordinate values of the weld and base metal sides in the horizontal direction at a certain thickness height position \({y}_{i}\). Similarly, let \({l}_{Ri}\) and \({l}_{Li}\) be the lengths of the interfaces of the base metal and HAZ, and weld metal and HAZ between \({y}_{i}\) and \({y}_{i+1}\), respectively, and \({\theta }_{i}\) be the angle that the interface makes with the shear direction (45°). The seven HAZ shape factors are defined as follows:

  1. 1.

    MAX angle of the HAZ boundary: For all \({\theta }_{i}\), the maximum value of \(\left|{\text{sin}}(2(45^\circ -{\theta }_{i}))\right|\); the value is observed to be 0 when \({\theta }_{i}\) coincides with the shear direction (45°) and 1 when \({\theta }_{i}\) is perfectly perpendicular to the load axis, when shearing does not occur at the interface.

  2. 2.

    HAZ global inclination: Absolute value of the difference between the center positions of the HAZ interface at the top and bottom free edges.

  3. 3.

    HAZ width (standard deviation): Standard deviation of \(\{{w}_{i};i=1\dots m\}\).

  4. 4.

    MAX HAZ width at the free surface: Larger HAZ length in the x-axis direction at the upper and lower free edges.

  5. 5.

    HAZ length: Root mean square of the interface length \(\{{l}_{R}=\sum_{i}{l}_{Ri}, {l}_{L}=\sum_{i}{l}_{Li}\}\) between the base metal and weld metal sides.

  6. 6.

    HAZ width (maximum): Maximum value of \(\{{w}_{i};i=1\dots m\}\).

  7. 7.

    HAZ width (average): Average value of \(\{{w}_{i};i=1\dots m\}\).

Fig. 5
figure 5

Schematic figure of seven proposed HAZ shape factors for a HAZ shape

Full size image

The minimum HAZ length is 8.0 mm, which is equal to the specimen height. The maximum of MAX HAZ width at the free surface and the HAZ width (average) do not exceed the HAZ width (maximum). In this study, we set \(m=100\).

4.2 Tandem Bayesian model search results

Using the HAZ shape factors described above, we search for the welding conditions with longer rupture life according to the algorithm described in Fig. 3 based on the tandem Bayesian model. The number \(n\) of selected candidates \({x}_{{\text{all}}}\) for a single step is \(10\).

Table 5 shows the results of the search for each \(BL\) obtained in each round: the coefficient for the normalized HAZ shape factors (0 at the coefficient column means that the HAZ shape factor was not selected), the number of steps in each round, the number of additional training data for both \({D}_{{\text{train}}}^{GP}\) and \({D}_{{\text{train}}}^{BL}\), the posterior probability, and the maximum rupture life obtained until that round \(T{r}_{{\text{target}}}\), which is the longest rupture life obtained up to that round.

Table 5 Coefficients for BL for normalized HAZ shape factor, number of steps, number of added \(GP\) and \(BL\) model, and number of added training data whose \(Tr(FEM)>T{r}_{target}\) and \(T{r}_{target}\) at each round
Full size table

The total numbers of data added in the search are 137 in total for \({D}_{{\text{final}}}^{GP}\) and 43 for \({D}_{{\text{final}}}^{BL}\), which is about 0.02% of the size for \({x}_{{\text{all}}}\), the search space. While adding this small number of forward-calculated data, \(T{r}_{{\text{target}}}\) increased by 70 h from 580 to 650 h, and a weld shape with a 12% longer rupture life than the highest value of the initial training data could be found. The longest creep life was found in round 2, and 82 forward calculations from the weld conditions to HAZ shape and 20 calculations to rupture life were added up to that point.

Figure 6 shows the distribution of the rupture life for the data added in each round. The \(FoM\) of the data at the last step of each round is shown in a color with reference to the color bar on the right side of the figure. \(T{r}_{{\text{target}}}\) of each round is indicated by a black arrow. The maximum value of the training data is indicated by larger cross marks. Each time \(T{r}_{{\text{target}}}\) is updated, the explored rupture life moves toward the longer side, especially in round 2, where the leap is particularly large. The \(FoM\) is gradually decreasing, indicating that it is becoming difficult to update the maximum creep rupture life in further searches. Ten candidates were proposed per step, but the number of additional data was smaller than that of steps because of overlap with the data already included in the training data. In round 2, 63 data were added in 9 steps, indicating that the tandem Bayesian model changed and new candidates were proposed as a result of the \(BL\) update in round 1. In round 3, the number of steps and the number of additional data decreased significantly. In round 4, no more data were added in the first step, and neither \(GP\) nor \(BL\) was updated, and the search was terminated. The posterior probabilities and coefficients of each explanatory variable of the \(BL\) model selected in rounds 3 and 4 remained almost unchanged. This indicates that the search using this algorithm was heading toward early convergence with a small number of updates.

Fig. 6
figure 6

The rupture life distribution of the search and the evolution of the \(T{r}_{target}\) at each round. The calculation was performed at the testing condition of 873 K, 100 MPa. The color bar indicates the value of \(FoM\) for each added data

Full size image

4.3 HAZ shape factors explaining rupture life and predictability of rupture life

The results of model selection for \(BL\) are analyzed as follows. The HAZ width (maximum) and HAZ length were always chosen as explanatory variables for all rounds. Their coefficients are negative for each model selected for every round, which means that the smaller the selected HAZ shape factor, the longer the rupture life. This is consistent with the qualitative observation in the initial training data that the HAZ shape with smaller HAZ widths and interface bends has a longer rupture life. The coefficient of HAZ width (maximum) on the normalized rupture life does not change at around \(-0.5\), whereas the coefficient of the HAZ length gradually decreases in their absolute values from \(-0.94\) to \(-0.61\) as the round advances. Instead, the absolute value of the coefficient of the MAX HAZ width at the free surface, which is selected after round 2, increases. However, the posterior probability of the model finally selected in round 4 is only 21.6%, so it is difficult to examine the superiority of these three explanatory variables over the others only by the examination of this selected model.

Figure 7 shows the coefficient heat maps for the top 20 by the Bayesian free energy, calculated in each round in decreasing order. The posterior probabilities for each model are shown at the bottom. In all rounds, the Bayesian linear regression model including HAZ length is always selected, and HAZ length is an essential explanatory variable. HAZ maximum width is also an essential explanatory variable, as it is adopted in most of the top models, although it is in a competitive relationship with HAZ average thickness. MAX HAZ width at the free surface is rarely selected in Fig. 7a, but in Fig. 7c, d. In most cases, the coefficient of MAX HAZ width at the free surface is estimated to be positive initially, whereas it is estimated to be negative in round 1 and later. This may be due to the addition of data with relatively long rupture lives in round 1. After round 1, there is no significant change in the trend of model selection or positive/negative coefficients, and finally, in round 3, the MAX HAZ width at the free surface is selected for all the top models, which is an essential explanatory variable. The posterior probabilities of the selected \(BL\) models, which consist of the three HAZ shape factors mentioned above, are significantly higher than those of the other models, especially in rounds 2 and 3. These results indicate that these three explanatory variables are selected as the factors necessary for explaining rupture life compared with the other explanatory variables.

Fig. 7
figure 7

Coefficient heat map matrices of top 20 (20 smallest FE) Bayesian linear regression (BL) model selection on initial (a), round 1 (b), round 2 (c), and round 3 (d). Color indicates the magnitude of the coefficients versus the normalized explanatory variables. The posterior probability for each model is shown at the bottom of each column

Full size image

Figure 8 shows the relationship between the two HAZ shape factors, HAZ width (maximum) and HAZ length, for the explored data. The color of each plot shows the rupture life with the reference to the color bar on the right side of the figure. The data added in the search are indicated by cross marks which are concentrated in the range of 0.4 to 0.8 mm for the maximum HAZ width and in the range of 8.0 to 8.3 mm for the HAZ length. The range is shown in the inset. If the HAZ shape can be freely designed, the shape with the smaller values for both HAZ shape factors, i.e., in the direction of the arrow shown in the inset, should be targeted in order to extend the rupture life. In practice, it is difficult to reduce both HAZ shape factors simultaneously, and the proposed candidate points of the search are concentrated near the gray dashed line. The line appears to be a Pareto front that can be realized by weld conditions. As will be discussed in the next section, reducing the heat input of welding can reduce the HAZ width, but the HAZ length tends to increase owing to bending at the HAZ interface. It is considered that there is probably no feasible HAZ shape within \({x}_{{\text{all}}}\) below this region to the left.

Fig. 8
figure 8

Plot of two explanatory variables for Bayesian linear regression, Max HAZ width vs HAZ length in the training data for each round. The color bar indicates the rupture time for each data by FEM, at the testing condition of 873 K, 100 MPa. The inset figure is the enlargement of areas where additional data from round 1 to 3 exists. The thick dashed line within the inset shows the Pareto front of the search, and the arrow indicates the direction to longer Tr

Full size image

Figure 9 shows the rupture life prediction of the tandem model for the initial data and for each round up to the third, versus the finite element analysis results of the data contained in \({D}_{{\text{train}}}^{BL}\), at the creep test conditions of 873 K and 100 MPa. Training data sets are as follows: the initial data set (a), the data set of a and the data obtained from the first round (b), the data set of b and the data obtained from the second round (c), and the data set of c and the data obtained from the third round (d). The initial data in Fig. 9a and the training data obtained from each round in Fig. 9b, c, and d are shown by cross marks. For the data added in round 2 (Fig. 9c) that have a rupture life than 600 h (shown within a red circle), the predicted rupture life is 600 h, which is an underestimation. This trend remains unchanged in round 3 (Fig. 9d). This suggests that it is difficult to predict from the short rupture life side to the long rupture life side with a single \(BL\) model. It suggests that a separate \(BL\) model is needed for the long rupture life side.

Fig. 9
figure 9

Tandem Bayesian model prediction results versus the FEM result on initial (a), round 1 (b), round 2 (c), and round 3 (d), at the testing condition of 873 K, 100 MPa. Training data sets are as follows: initial data set (a), the data set of a and the data obtained from first round (b), the data set of b and data obtained from the second round (c), the data set of c and data obtained from the third round (d). The cross marks show the plot of added training data for each round

Full size image

4.4 Weld conditions and rupture life predictability of additional training data

Table 6 shows the welding conditions, HAZ shape factors, \(FoM\) at the last step of each round, and rupture life of the data added to \({D}_{{\text{train}}}^{BL}\) for each round. Data where the rupture life exceeds the \(T{r}_{{\text{target}}}\) for the round are shown in bold.

Table 6 Welding parameters, HAZ shape factors, and rupture time of predicted and calculated for added training data at each round
Full size table

For the ten welding conditions added in round 1, the first layer thickness DL1 ranged from 2.8 to 3, and the second layer thickness DL2 was 3.2 (upper limit). The heat inputs HV1 and HV3 for the first and third layers were both 1000 J/mm (lower limit), and the heat input HV2 for the second layer ranged from 1400 to 1800 J/mm. In other words, the welding conditions were chosen to keep the heat input per unit volume at the weld low in each layer while it slightly increases it in the second layer. By increasing the heat input in the second layer, the search tends to try to realize HAZ that uprights to the tensile direction while avoiding bending, and the HAZ length is actually suppressed to the lower limit value of 8 mm [18]. However, the results of calculations in round 1 using \(wf2\) for these candidates showed that only two of the ten conditions achieved a long rupture life that exceeded \(T{r}_{{\text{target}}}\), and were given that the probability of exceeding \(T{r}_{{\text{target}}}\) (i.e., \(FoM\)) was 0.6 or greater for all ten conditions, the accuracy of the proposal was low. This may be because the initial training data were not sufficient for capturing the response phase over a wide area.

In contrast, almost all the welding conditions after round 2 are in the upper limits for both DL1 and DL2, whereas HV1, HV2, and HV3 are all in the lower limits. As a result, the HAZ width (maximum) is greatly reduced. In detail, only two round 2 candidates exceeded 0.5 mm, compared with 0.6 to 0.7 mm for round 1. On the other hand, owing to the effect of decreasing W2, which is the second layer heat source width, bending became more prominent and the HAZ length increased slightly; in fact, the HAZ length in round 1 was around 8.1 mm, and it increased to around 8.25 mm after round 2. Thus, the trend of the proposed candidates changed significantly from round 2. Note that in round 2, all of the 10 proposals have longer rupture lives than \(T{r}_{{\text{target}}}\), including the longest rupture life obtained in this study. The above results show that it is easier to obtain a longer creep life by prioritizing the HAZ width, even if the interface bends slightly. In addition, note that the tandem Bayesian model learning from the additional data in round 1 was able to quickly capture this trend.

Figure 10 shows the cross-sectional HAZ shape for the additional training data for \({D}_{{\text{train}}}^{BL}\), arranged in order of increasing rupture life. All HAZ shapes are narrower than in the initial training data in width, and the overall shape appears to be more upright perpendicular to the tensile direction. The HAZ shapes with rupture lives of 610 h or longer are even narrower than those with rupture lives of less than 610 h, which is considered to be achieved owing to the low heat input. Two bends protruding toward the weld metal exist at the interface between the base and weld metals, which are considered to be caused by the boundary of the weld layer.

Fig. 10
figure 10

Cross-sectional HAZ shapes of all models added to Bayesian linear regression training data by the end of the calculation, arranged by rupture time at 873 K, 100 MPa

Full size image

5 Discussion

The proposed optimization algorithm for rupture life based on the tandem Bayesian model in this study works well and succeeds in finding the welding condition which gives the longest possible rupture life up to 650 h with a small number of trials from the initial data. To discuss the validity of the tandem Bayesian model in terms of generalizability, we compared the Bayesian free energy between the tandem Bayesian model and two other methods that predict the rupture life directly from the welding conditions.

Bayesian free energy calculations were performed on \({D}_{{\text{final}}}^{BL}\) to compare these three models, the final prediction model of the tandem Bayesian model, the model obtained by a Gaussian process (indicated as \(DirectGP\)), and the model obtained by a Bayesian linear regression model (indicated as \(DirectBL\)) as a model to directly predict the rupture life \({y}_{Tr}\) from the welding conditions.

The framework for calculating the Bayesian free energy of the tandem Bayesian model is shown in Fig. 11. Although \({D}_{{\text{final}}}^{GP}\) is used in the construction of the tandem Bayesian model to construct \(GP\), note that only the 43 points included in \({D}_{{\text{final}}}^{BL}\) appear in the Bayesian free energy evaluation, and data on \({D}_{{\text{final}}}^{GP}\) do not appear explicitly in the overall prediction of rupture life with the welding condition \(x\) as input. Note that the prior distribution of \(BL\) is the posterior distribution \(p({y}_{{\text{shape}}}|x)\) of the \(GP\) in the first part, and the likelihood \(p({M}_{{\text{tandem}}}|{y}_{Tr})\) of the tandem Bayesian model \({M}_{{\text{tandem}}}=\{BL,GP\}\) is as follows:

Fig. 11
figure 11

Calculation scheme of Bayesian free energy of tandem model

Full size image
$$p({M}_{{\text{tandem}}}\left|{y}_{Tr}\right. )\propto p({y}_{Tr}\left|{M}_{{\text{tandem}}} \right.)p({M}_{{\text{tandem}}} )\propto \int {\text{d}}{y}_{{\text{shape}}}\mathrm{ d}BL {\text{d}}GP p({y}_{{\text{shape}}},y\_Tr,BL,GP;x)p(BL)p(GP)$$
(6)

We can rewrite the integrals by distributing \({y}_{{\text{shape}}}\), \(BL\), and \(GP\) according to the graphical representation causality of the tandem Bayesian model in Eq. (6) as:

$$\int {\text{d}}{y}_{{\text{shape}}} [\int {\text{d}}GP p\left({y}_{{\text{shape}}}|GP;x\right)p\left(GP\right)\int {\text{d}}BL\, p\left({y}_{Tr}|BL,{y}_{{\text{shape}}};x\right)p\left(BL\right)]$$
(7)

\(\int {\text{d}}GP p\left({y}_{{\text{shape}}}|GP;x\right)p\left(GP\right)\) is a predict distribution for \(GP\) and \(\int {\text{d}}BL p({y}_{Tr}|BL,{y}_{{\text{shape}}};x)p(BL)\) is a likelihood function of \(BL\). Equation (7) follows an integration of \({y}_{{\text{shape}}}\) as:

$$\int {\text{d}}{y}_{{\text{shape}}}p({y}_{Tr}\left|{y}_{{\text{shape}}},x\right)p\left({y}_{{\text{shape}}}|x\right)$$
(8)

But Eq. (8) is nothing but the expected value \({{\mathbb{E}}(p\left({y}_{Tr}|{y}_{{\text{shape}}},x\right))=\langle p\left({y}_{Tr}|{y}_{{\text{shape}}},x\right)\rangle }_{p({y}_{{\text{shape}}}|x)}\) for the probability of \({y}_{Tr}\), \(p\left({y}_{Tr}|{y}_{{\text{shape}}},x\right)\), in the sense of \(p\left({y}_{{\text{shape}}}|x\right)\). Let \(\{{x}_{k},{y}_{{\text{shape}},k},{y}_{Tr,k}\}\) be a training data contained in \({D}_{{\text{final}}}^{BL}\). The expected value of \(p\left({y}_{Tr,k}|{y}_{{\text{shape}},{\text{k}}},{x}_{k}\right)\)(i.e., Eq. (8)) for this data was calculated as the m-point sampling average \(\frac{1}{m}{\sum }_{i=1}^{m}p({y}_{Tr,k}^{i}=T{r}_{{\text{FEM}},k};{x}_{k})\), which was then multiplied for all 43 data to obtain the neighborhood likelihood \(p({M}_{{\text{tandem}}}|{y}_{Tr})=\prod_{k=1}^{43}{\mathbb{E}}(p\left({y}_{Tr,k}|{y}_{{\text{shape}},k},{x}_{k}\right))\) to determine the Bayesian free energy as shown below:

$$FE\left({M}_{{\text{tandem}}}\right)=-{\text{log}}\,p({M}_{{\text{tandem}}}|{y}_{Tr})$$
(9)

Here, we set \(m=100000\); and furthermore, the final value was calculated as the average of 10 calculations of Eq. (9) with different seed values for random sampling for \(GP\) corresponding to each HAZ shape factor.

For \(DirectGP\) construction, hyperparameter tuning was performed for the kernel function using the same method (described in 2.3) as that used for the \(GP\) in the first stage of the tandem Bayesian model, and the Bayesian free energy was calculated using log_marginal_likelihood function which is implemented in scikit-learn. \(DirectBL\) was constructed in the same way as \(BL\) for the tandem Bayesian model, with model selection performed on the welding conditions as explanatory variables to construct a model with the smallest Bayesian free energy.

Table 7 shows the RMSE for the Bayesian free energy and the posterior probability obtained from the normalized rupture life for each model. The tandem Bayesian model has the lowest Bayesian free energy, followed by \(DirectGP\) and \(DirectBL\). The difference between the Bayesian free energy of the tandem Bayesian model and that of \(DirectGP\) is 17, which clearly indicates that the tandem Bayesian model is valid feasible for predicting rupture life from welding conditions. Within these three models, the posterior probability of the tandem Bayesian model was 1 and those of the other two models were both 0. The tandem Bayesian model explains the data very well, and it can be concluded that the method is highly generalizable. On the other hand, the RMSE, which indicates the prediction error for the training data, is the smallest for \(DirectGP\), followed by tandem Bayesian model and \(DirectBL\). At first glance, the fitting performance of the tandem Bayesian model appears to be inferior to that of \(DirectGP\), but \(DirectGP\) has a higher degree of freedom in proportion to the number of training data, which may cause overlearning. The tandem Bayesian model proposed in this study can incorporate a linear model by passing an appropriate HAZ shape factors in between, rather than directly predicting the rupture life from the welding conditions with a nonlinear model. It is inferred that this suppresses overlearning and improves generalization performance.

Table 7 Bayesian free energy, marginal likelihood, root-mean-square error of normalized rupture life, LOOCV error and ratio of LOOCV error versus RMSE for the tandem Bayesian model, DirectGP, and DirectBL
Full size table

6 Conclusion

The inverse problem was tackled to find the optimum welding conditions that give a longer creep rupture life in Type IV cracking for the heat-resistant steel weld joint. The target was 2 1/4Cr–1Mo weld joints by GTAW with three-layer I-shaped grooves in a plate. The search space was huge, with more than 5.76 million conditions, and it was practically impossible to compute it exhaustively. Therefore, a method of sequential optimization by constructing a surrogate model was adopted. As a surrogate model, the tandem Bayesian model consisting of the nonlinear Gaussian process \((GP)\) for the HAZ shape factors from the welding conditions and the Bayesian linear regression model \((BL)\) for the rupture life from the HAZ shape factors was developed. The original double-loop optimization algorithm combining \(GP\) updating (step) and \(BL\) updating (round), which corresponded to the tandem Bayesian model, was proposed and used to search for welding conditions that maximize the creep rupture life of the weld joint starting from the initial training data of 49 conditions selected by the design of experiment.

  • •In the second round (\(BL\) update), we were able to find a welding condition that improved the rupture life by 12% over the initial training data. The number of forward calculations required to obtain this result was 63 from the welding condition to the HAZ geometry, which is relatively light in calculations, and 20 up to the rupture life, which is heavy in calculations.

  • •The search converged in round 4, where 137 forward calculations from welding conditions to HAZ geometry and 43 calculations to rupture life were added before convergence. This is only about 0.02% of the total search space.

  • •The Bayesian linear regression model selected for this method included the HAZ width (maximum), HAZ length, and MAX HAZ width at the free surface, and their coefficients were negative. It was found that a HAZ shape with a small overall width and a straight interface perpendicular to the tensile direction is advantageous in extending the rupture life.

  • •The HAZ shape that exhibited the longest rupture life was characterized by a small HAZ width, an overall upright toward the tensile direction, and slight bending. These findings are consistent with the conventional qualitative understanding of HAZ shapes that gives a long creep life.

  • •The distributions of the HAZ width (maximum) and HAZ length obtained from the search suggest that there was a trade-off between these two HAZ shape factors and that it was difficult to minimize both at the same time.

  • •The Bayesian free energy of the tandem Bayesian model was lower than that of the model predicting the rupture life directly from the welding conditions by the Gaussian process, indicating that the present tandem model was more reasonable as a surrogate model.