Numerical Estimation of Material Properties in the Electrohydraulic Forming Process Based on a Kriging Surrogate Model

High-speed forming processes, such as electrohydraulic forming, have recently attracted attention with the development of forming technology. However, because the high-speed operation (above 100 m/s) raises safety concerns, most experiments are conducted in a closed die, which hides the forming process. Therefore, the experimental process can only be observed in a numerical simulation with accurate material properties. The conventional quasistatic material properties are improper for high-speed forming simulations with high strain rates (>102 s−1). In this study, the material properties of Al 6061-T6, which reflect the deformation behavior in the high-strain-rate region, were investigated in a numerical approach based on a reduced order model and a surrogate model in which the numerical results resemble the experimental results. The strain rate effect on the material was determined by the Cowper–Symonds constitutive equation, and two strain rate parameters were predicted. The surrogate model takes two material parameters as inputs and outputs a weighting coefficient calculated by the reduced order model. The surrogate model is based on the Kriging method to reduce the simulation cost. Next, the optimal material parameters that minimize the error between the surrogate model and the experiments are estimated by nonlinear least-squares optimization using a genetic algorithm and the constructed surrogate model. The predicted optimal parameters were verified by comparing the results of the experiment, numerical simulation, and surrogate model.


Introduction
Forming processes under quasistatic conditions have developed alongside many technological advances in sheet metal forming at high speeds (above 100 m/s), such as electrohydraulic forming (EHF), electromagnetic forming (EMF), and explosive forming (EF) [1][2][3]. ese forming processes deform a sheet metal within 1 ms by a high-energy pulse, such as a high-voltage discharge or a high-energy explosion. When a sheet metal is deformed under highspeed conditions, the formability is improved by the diesheet interaction and the inertia effect of the sheet. Many researchers have demonstrated the improvement of the formability through various experiments and finite element analyses [3][4][5][6][7]. e EHF process has attracted more attention than other high-speed forming processes because it overcomes the bouncing effect of the sheet [8], which is a crucial problem in the EMF process [9], and is safer than EF, which directly uses explosives. EHF exploits electric energy as the forming resource and consists of a resistance-inductance-capacitance (RLC) circuit. Figure 1 presents a typical configuration of the EHF apparatus.
As shown in Figure 1(a), when the electrical energy (in kJ) in the typical EHF equipment is discharged from a capacitor bank, it increases the pressure in the water, and the water becomes a plasma. e quickly expanding plasma cavity creates shockwaves in the water, which propagate the high pressure and form a sheet into the die.
High-speed forming processes such as EHF cannot be observed during tests because the sheet deformation is finished within 1 ms and conducted in a closed room for safety purposes. erefore, we need a dependable finite element model to describe the forming process. To increase the reliability of the numerical model, the obtained material properties must minimize the error between the experimental and simulation results. e deformation behavior of the material depends on the strain rate [10], so it is important to understand the material properties according to specific strain rate condition. In general, material properties under high-strain-rate conditions (strain ranges 10 2 -10 4 s −1 ) are acquired in a Split-Hopkinson pressure bar test [11,12]. Besides, many experimental studies were conducted to obtain dynamic material properties [13][14][15][16][17]. Although the experimental approach directly obtains the material properties, it is both costly and time consuming.
To avoid this problem, we proposed a numerical approach that estimates the parameters in the material model by a surrogate model with nonlinear least-squares optimization. e numerical model estimates the parameters that minimize the error between the numerical and experimental results. e target material was Al 6061-T6, which is used in many industrial fields, such as vessels, airplanes, and automobiles. e material constitutive equation was the Cowper-Symonds model with two unknown parameters to be estimated. e parameters were estimated by a reduced order model (ROM) with weighting coefficients predicted by the surrogate model.
Previous study [18] described a parameter estimation method based on artificial neural network (ANN) method. However, ANN generates slightly different models according to internal weight factors and it constructs various models according to the number of nodes and hidden layers even at the same input and output, which means that ANN is suitable for the problem which should check the output tendency at given input, but it may not be appropriate when an accurate value should be obtained as in this paper. erefore, we adopted the original Kriging method to construct the surrogate model. Unlike ANN, the original Kriging model shows one deterministic trend value that can represent whole samples calculated by maximum likelihood estimation. As a result, from given input and output samples, one surrogate model is constructed, and we can estimate outputs at arbitrary inputs more clearly than when we use ANN surrogate model. e rest of this paper is structured as follows. Section 2 elucidates the theoretical background of the ROM and the surrogate model, which are based on principal component analysis (PCA) and a Kriging method, respectively. Section 3 describes the experiments and numerical simulations of the EHF free-bulging test. e basis vector and weighting coefficients based on the ROM simulation results are also delineated here. Section 4 verifies the constructed surrogate model on training and test samples. In Section 5, the optimized parameters that minimize the error between the results of the experiment and the surrogate model are predicted by a genetic algorithm. Section 6 concludes the paper and proposes future works.

ROM Based on Principal Component Analysis.
To estimate the material parameters of Al 6061-T6, we adopted the ROM based on PCA and the surrogate model based on the Kriging method.
Reduced order model (ROM) refers to a model that has reduced dimension while maintaining the characteristics of the original high-dimensional model as much as possible for computational convenience. For example, in aerospace fields [19,20], fluid dynamics problems have degrees of freedom upward of 10 6 which requires hundreds of hours to complete the simulation. erefore, the method of reducing the dimension of the model while keeping the characteristics may be the best way to reduce the computational time.
is process can be performed by principal component analysis (PCA). PCA is a mathematical method for reducing highdimensional data to low-dimensional data by eigenvalue decomposition of covariance matrix of data. Correlated high-dimensional data are mapped to a lower-dimensional orthogonal space without linear correlation through orthogonal transformation [21,22].

Mathematical Problems in Engineering
Let y c (x; θ) ∈ R m be mean-centered simulation outputs at x ∈ R with respect to inputs θ ∈ R p . In the ROM, the vector y c can be spanned by bases v i (x) r i�1 ∈ R m×r , where r is the number of bases. erefore, we can present y c as a linear combination of v i ∈ R m as follows: where w i (θ) is the weighting coefficient of basis vector v i ∈ R m , which can be calculated by orthogonally projecting y c onto a basis vector v i . We can reduce the order of y c by dropping the r − q insignificant basis vectors. To this end, we first use the significant q basis vector v i (x) and approximate y c on a lowdimensional space as follows: e basis vectors for the ROM are computed by the PCA technique, which is widely adopted for dimensionality reduction and feature extraction.
Let Y � [y 1 , y 2 , . . . , y n ] ∈ R m×n be the total output and Y c be the mean-centered Y calculated by Y − Y mean , where Y mean � (1/n) n i�1 y i . e data Z, created by projecting the original meancentered data Y c into the lower-dimensional space, can be expressed as a linear combination of Y c , namely, as Z � v T Y c . To maximize the variance of Z, the following equation is solved by linear algebra: where Σ is the covariance matrix of Y c . Applying the Lagrange multiplier method under the constraint ‖v‖ � 1, we finally obtain In basic linear algebra [23], v and λ are defined as the eigenvector and eigenvalue of Σ, respectively. e eigenvectors are the principal components of v that maximize the variance of Y c . erefore, the basis vectors for the ROM are the eigenvectors of the covariance matrix of the original data Y c .

Surrogate Model Based on the Kriging Method.
e ROM process calculates the weighting coefficient w i (θ) in equation (1) by orthogonally projecting y c onto v i . e weighting coefficients are calculated by the vector projection formula as follows: where y c is the vector obtained by projecting y c onto the basis v i . For known input parameters θ j ∈ R p , such as the training and test samples in this study, the basis weighting coefficients w i (θ j ) are easily obtained using equation (5). However, when the inputs within a given domain Ω � [θ lb × θ ub ] ∈ R p are unknown, as in the present case, w i (θ j ) cannot be computed. erefore, they were estimated by the ordinary Kriging method.
e Kriging method, also known as Gaussian process regression, was developed by Georges Matheron as an estimation technique for geostatistical fields and has since been adopted in many industrial fields [24]. e Kriging method constructs the relationship between the input and  output training data and estimates the output for an unknown input as a weighted linear combination of the observed outputs. e Kriging model consists of two terms, a deterministic trend f(x) and a stochastic term Z(x): e Y(x) in equation (6) can be assigned as the weighting coefficients in the ROM.
In the Kriging formulation, a set of observations composed of the trend f(x) and a random variable Depending on their definition of f(x), Kriging formulations are classified as simple when f(x) � α (known), ordinary when f(x) � α (unknown), and universal when is study adopts the ordinary Kriging formulation, in which the deterministic trend is unknown and constant, that is, Because α is constant and deterministic, the covariance . e covariance between Z i and Z j is determined by correlating x i and x j as follows: where To calculate Ψ, the Kriging method uses the following Gaussian correlation function: Here, θ l ∈ R p is a hyperparameter that determines the correlation strength between x i and x j . A high θ l means that only the neighboring sample points are correlated, so the output Y(x) exhibits a highly nonlinear behavior.
In the Kriging method, Z(x) also varies locally and follows the multivariate follows the multivariate Gaussian distribution N(α1 n , σ 2 Ψ), where 1 n is a vector of n ones. e probability distribution of the output has three unknown parameters: α, σ 2 , and θ. In the original Kriging formulation, these parameters are estimated by the following log-likelihood function: Applying the maximum likelihood estimation (MLE), α and σ 2 are, respectively, determined as follows: Note that θ cannot be directly obtained from the MLE because Ψ is a function of θ. Instead, θ is estimated by nonlinear optimization using the concentrated log-likelihood function: Once all parameters α, σ 2 , and θ are estimated, we can predict y(x 0 ) at some unknown input x 0 . When predicting y(x 0 ), Kriging assumes that the augmented output y � [y; y(x 0 )] ∈ R n+1 has the same Gaussian distribution as the original output y.
erefore, y is distributed as Similarly to the previous parameter estimation process, y(x 0 ) at the unknown input x 0 is finally predicted by MLE using the following equation:  Mathematical Problems in Engineering In this expression, Ψ −1 (y − α1 n ) is generated by training the inputs, and ψ is calculated using the known training data and the unknown input x 0 .

Electrohydraulic Forming Test Using a Free-Bulging Die
3.1. Experiment. e experimental apparatus of the EHF test comprises two electrodes, a free-bulging die, a chamber, a deformable sheet, and a capacitor bank with a control system. Under the high-speed conditions of the EHF test, the formability of the sheet is influenced by contact with the die [3,7]. erefore, we applied a free-bulging die to avoid a contact between a sheet and a die. e experiment was configured in the RLC circuit presented in Figure 1(b). In this circuit, R t is the total resistance without the wire, R w is the wire resistance, L is the total inductance, and C is the capacitance of the capacitor bank. When the initial voltage is charged to V in the capacitor bank and then discharged, the resulting current destroys the wire located between the two electrodes set inside the chamber, creating high-pressure shockwaves.

Mathematical Problems in Engineering
When the shockwaves arrive at the sheet, they deform it into a die shape. e capacitance and maximum input voltage of the capacitor bank were 333 μF and 15 kV, respectively.
In the present study, the input voltage was charged to 8 kV, and the Al 6061-T6 sheet dimensions were 250 × 250 mm. e sheet thickness was 1 mm. e chemical composition of Al 6061-T6 is presented in Table 1.
When the 8 kV voltage was discharged from the capacitor bank, a high current flow was generated through the whole electric circuit of the EHF apparatus. A presented current curve in Figure 2 was measured from the EHF experiment by using Rogowski coil. Unlike EMF, which generates a sinusoidal current wave [9], the current curve of EHF shows a single peak because the electric wire placed between the two electrodes is broken by the high current after the voltage discharge. e reproducibility of the EHF experiment was confirmed in a previous study [25], so only one experimental result is presented here. e z-axis displacements of the deformed sheet were measured using a three-dimensional scanning machine. Because the deformed shape was symmetric about the z-axis, we present only half of the results, i.e., the results at x ≥ 0. Section 5 will compare our experimental results with those of other approaches.

Numerical Simulation of Electrohydraulic Forming.
A numerical model of the EHF free-bulging test was built in LS-DYNA. EHF analysis requires a different approach from general forming analysis because both the structural and fluid parts move simultaneously. Here, the fluid-structure interaction (FSI) was implemented by an arbitrary Lagrangian-Eulerian (ALE) approach. In LS-DYNA, the FSI is simply defined by a constraint keyword named "Constraint Lagrange in Solid" [26]. e fluid parts (plasma and water) were modeled as ALE materials, and the structural parts (sheet, chamber, and die) were modeled as Lagrangian elements. e whole numerical model of EHF is presented in Figure 3. As described in the previous section, the deformation shape of the sheet is symmetric about the z-axis. erefore, we established a ¼ model and applied symmetric conditions in the x-z and y-z planes.
Electrical energy was applied to the fluid parts of the plasma part (1/4 of a sphere with a radius of 1 mm) to replicate the electrical circuit in the experiment. e current and voltage obtained from the experiment were inputted into the plasma part by an equation-of-state keyword, including the electric energy term [26]. e free-bulging die and the chamber were assumed as rigid bodies, and the sheet was modeled as a deformable four-node shell model to reduce the computational time.
e high-speed deformation of the sheet under the highstrain-rate effect was modeled by the following Cowper-Symonds constitutive equation: where σ eff and ε eff denote the effective stress and strain, respectively, and _ ε eff is the effective strain rate. A, B, and n are the quasistatic material parameters, and C and p are the high-strain-rate material parameters. e material parameters of Al 6061-T6 were determined in a quasistatic tensile test conducted on an Instron universal testing machine at 0.0007 s −1 , yielding A � 291.8 MPa, B � 451.5 MPa, and n � 0.666 [18]. Both C and p were unknown parameters θ ∈ R 2 to be found by inverse estimation in the Kriging surrogate model. Before constructing the surrogate model,    Figure 5. e legend in Figure 5 represents the number of the training and test samples, respectively. From the z-axis displacements of the training samples, we determined the basis vectors of the ROM process. Twenty eigenvalues and eigenvectors were obtained from the collection of z-axis displacements Y ∈ R 39×20 by the snapshot method. e cumulative sum of the normalized eigenvalues is presented in Figure 6. e first eigenvector has an eigenvalue of 0.9886, meaning that the variance of the original sample data is preserved by more than 98% (Figure 6 (inset)). As presented in Figure 7, the shape of the first eigenvector closely resembles the deformation of the sheet in Figure 5, indicating that the first eigenvector is the most   erefore, by applying the ROM process and choosing only one eigenvector, the high-dimensional numerical simulation of the EHF was reduced to a one-dimensional space. In terms of one basis vector, equation (1) becomes e weighting coefficient w 1 (θ) was calculated from the orthogonal projections of the 20 training samples. e Kriging surrogate model was then constructed using the 20 training samples and 20 weighting coefficients as the inputs and outputs, respectively. Finally, the constructed surrogate model was verified on the 20 test samples. To improve the prediction accuracy, we normalized the inputs and outputs.
Because the Kriging surrogate model provides only the scalar output, it should be constructed as many times as the dimensions of the desired outputs. However, as the output of the present surrogate model is only one weighting coefficient, we constructed a single Kriging surrogate model. Figure 8 shows a flowchart of the Kriging-and ROMbased surrogate model.

Surrogate Model Verification
e Kriging surrogate model for the EHF simulation was constructed in MATLAB 2019b, and the estimated hyperparameters are presented in Table 2.
As the hyperparameter θ is more influential at larger values, the material parameter C affects the weighting coefficient more than p [28]. e predicted weighting coefficients are plotted against the actual values in Figure 9. e actual values were calculated by orthogonally projecting the EHF simulation results onto the basis vector, and the predicted values were obtained from the surrogate model. All weighting coefficients determined by the model lie on the actual-predicted line of slope 1.0 for the training samples, but those of the test samples are scattered around the actual-predicted line. e coefficients of determination (R 2 ) were 1.0 for the training samples and 0.77012 for the test samples. R 2 is typically between zero and one, with high values indicating good predictability of the surrogate model. Although the test samples achieved a lower R 2 value than the training samples, all results were clustered around the actual-predicted line (Figure 9(b)) and explained more than 77% of the estimation rate. We concluded that the weighting coefficient prediction capacity of the surrogate model was quite acceptable.
Using the weighting coefficients obtained from the surrogate model, we verified the z-axis displacement y i of the deformed sheet calculated using equation (15). e accuracy of y i was assessed by two error values: R 2 and the root mean squared error (RMSE): When the RMSE value is close to zero, a mathematical model makes accurate predictions. Table 3 presents the numerical verification results on the training and test datasets, and Figure 10 shows a graphical representation of these results. All R 2 values exceeded 0.998 (maximum 0.99997) on the training samples and 0.992 (maximum 0. 99982) on the test samples.
e RMSE values of all 40 samples were below 1.0 mm. Given that the z-axis displacement of the sheet is approximately 30 mm, RMSE values below 1 mm confirm the high predictability of the surrogate model.
Observing the predicted weighting coefficients in Figure 9(b), the R 2 value was lower for the test samples than for the training samples. However, as mentioned above, the test samples follow the trend line, so the surrogate model based on the ROM technique accurately predicted the z-axis displacement, as presented in Table 3.
e performance degradation on the test samples was caused by overfitting due to the high nonlinearity of the predicted weighting coefficients. Accordingly, the weighting coefficients of the first input C (see Figure 11) exhibited no obvious trend.

Optimizing the Material Parameters
In the constructed surrogate model, we found the optimum material parameters C and p that minimize the error between the experiment and the numerical simulation. With a misfit formula based on the L2 norm, we define the objective function as follows: where y i,est is the z-axis displacement determined by the surrogate model and y i,exp is the displacement measured in the experiment. Here, common error values, such as the mean squared error and RMSE, were unsuitable because we aimed at observing a large variation in the error values as the two parameters were varied. Because equation (16) contains a nonlinear term y i,est determined by the material parameters C and p, our selected optimization algorithm must handle nonlinear problems. erefore, we solved equation (16) by a GA, a widely used global optimization algorithm based on natural selection [29]. Starting from initial random candidates, the GA optimally solves the objective function through the genetic operations of selection, crossover, and mutation. e global nature of the GA avoids trapping the solutions in local minima; moreover, the initial values are randomly generated and need not be preset. Figure 12 presents the two-and three-dimensional plots of the error as a function of C and p, calculated using equation (16). After the GA optimization, the error was  e acquired optimal parameters were validated in three cases: experiment, LS-DYNA, and the surrogate model. e results in the three cases are compared in Figure 13, and the error values, R 2 , and RMSE, are listed in Table 4.
As presented in Table 4, all errors were acceptable. In the LS-DYNA versus surrogate model, the R 2 was 0.99828, confirming that the surrogate model based on the Kriging method can obtain reliable prediction results without expensive numerical simulations.
In conclusion, the material properties for Al 6061-T6 in the EHF process were determined as C � 8904.6 and p � 15.832. Using the obtained parameters, we can predict the deformation behavior of the sheet in the EHF process at arbitrary strain rates.
When the experimental result is compared with the numerical simulation with only quasistatic material property as shown in Figure 14, the comparison in the peak z-axis displacement shows almost 37 % error, which can support that Al 6061-T6 has strain rate dependency and the numerical simulation requires material properties that take into account the strain rates. e flow curves for Al 6061-T6 are described in Figure 15 according to different strain rates. Lee et al. [30] showed that the yield stress for Al 6061-T6 has a value between 450 and 650 MPa in the strain rate of 10 3 s −1 or more. Because the yield stress at a strain rate of 1000 s −1 in Figure 15 is about 545 MPa, it is concluded that Figure 15 shows reasonable results compared to the previous study [30].

Conclusions
In this study, the material parameters for Al 6061-T6 under the electrohydraulic forming process were estimated by a surrogate model based on a ROM and the Kriging method. Works and findings from this study are as follows: (1) Using 20 training samples extracted from a specific domain of the two material parameters C and p, we performed numerical simulations in LS-DYNA and extracted the eigenvectors and eigenvalues by principle component analysis. e first eigenvector represented the whole numerical model by preserving over 98% variance of the original simulation data, which was taken as the basis vector for the numerical simulation. (2) We constructed a Kriging surrogate model that inputs the material parameters and outputs the weighting coefficients of the basis vector, calculated for the training samples. e predictability of the constructed model was evaluated by the R 2 and RMSE of 20 test samples. e R 2 of the weighting coefficients was approximately 0.77 for the test samples, which is lower than for the training samples. However, the R 2 value of the z-axis displacement was very high (>0.9 for both training and test samples), endorsing the acceptable reliability of the surrogate model."   purpose, we employed a genetic algorithm that minimizes the error between the experiment and the numerical simulation, finally obtaining C � 8904.6 and p � 15.832.
However, this study has two shortcomings as follows: (1) e Kriging model gave poorer prediction results for the test sample than for the training sample, because the high nonlinearity of the predicted weighting coefficients led to the overfitting problem. In addition, because the original Kriging formulation estimates the trend line using a constant, it cannot properly predict highly nonlinear data. (2) e obtained material parameters may not be accurate under other conditions because they satisfy the objective function under a specific condition, 8 kV. In particular, the error contours differ within 1% in the upper white valley of Figure 12, indicating that several solutions satisfy the objective function.
erefore, it is necessary to adopt another prediction model that can accurately predict nonlinear data, such as a universal Kriging formulation or an artificial neural network.
In addition, we presented a methodology that acquires the material parameters in high-speed forming processes. To improve the accuracy of the parameter estimation, we suggest two improvements. First, the results under various conditions (such as varying voltage) should be obtained to find the parameters that satisfy all conditions. Second, a probabilistic approach, which considers the uncertainties in the experiment and the surrogate model, is required.

Data Availability
e data used to support the findings of the study are confidential and are not publically available.

Conflicts of Interest
e authors declare that they have no conflicts of interest.