Inverse parameter estimation to predict material parameters of the Cowper–Symonds constitutive equation in electrohydraulic forming process

The development of a reliable numerical simulation is essential for understanding high-speed forming processes such as electrohydraulic forming (EHF). This numerical model should be created based on the accurate material properties. However, dynamic material properties at strain rates exceeding 1000 s-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document} cannot be easily obtained through an experimental approach. Thus, this study predicted two material parameters in the Cowper–Symonds constitutive equation based on inverse parameter estimation, such that the parameters predicted using the numerical simulation corresponded well with those obtained from the experimental results. The target material was a 1-mm-thick Al 5052-H32 sheet. The comparison target included the final deformation shape of the sheet in the EHF-free bulging test at three input voltages of 6, 7, and 8 kV. For the inverse parameter estimation, the posterior distribution for the two parameters included a likelihood and a prior distribution. For the likelihood construction, a reduced-order surrogate model was developed in advance to substitute the numerical simulation based on ordinary Kriging and principal component analysis. Moreover, the error distribution of the bulge height between the experiment and reduced-order surrogate model was obtained. The prior distribution at 7 kV was defined as a uniform distribution, and the posterior distribution at 7 kV was employed as a prior distribution at 6, 7, and 8 kV. Furthermore, Markov chain Monte Carlo sampling was employed and the Metropolis–Hastings algorithm was adopted to obtain the samples following the posterior distribution. After the autocorrelation calculation for sample independence, the lag with an autocorrelation of ±0.02\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,0.02$$\end{document} interval was selected and every lagth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\mathrm {th}}$$\end{document} sample was obtained. The total number of acquired samples was 105\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{5}$$\end{document}, and the mean values were calculated from the obtained samples. Consequently, the numerical simulation with mean values displayed good agreement with the experimental results.


Introduction
Recently, several solutions have been developed in the automobile industry to improve fuel efficiency and reduce greenhouse gas emissions. To date, the most popular method is the weight reduction of the automotive parts, wherein the use of aluminum alloys and advanced high-strength steels and composites yields high-strength, lightweight products instead of the existing products obtained using mild steel. However, the application of these materials is limited owing to their low formability, which requires an advanced forming process in addition to the conventional quasistatic press forming (e.g., deep drawing). Owing to this increased demand, several high-speed forming processes have been developed [e.g., electrohydraulic forming (EHF), electromagnetic forming (EMF), and explosive forming (EF)]. In high-speed forming processes, the sheet undergoes deformation at speeds higher than 100 m/s and strain rates greater than 1000 s −1 . In particular, the formability of low-formability materials can be enhanced under this high-speed condition. Although the mechanisms of formability enhancement are not appropriately understood and have not been clearly verified, several prior studies [1][2][3][4][5][6][7][8][9][10][11][12][13] have attributed the primary causes of the formability enhancement to the inertial effect and void growth suppression in a sheet.
The EHF method uses a shockwave in the fluid as the forming pressure. The two electrodes in the fluid deliver a high current from the capacitor bank to the fluid (Fig. 1). Simultaneously, a high-pressure shock wave is generated between the two electrodes owing to the high electrical energy, and the shock wave propagates through the fluid toward the sheet. The EHF process poses several advantages as follows: (1) less capital investment as the EHF only requires a one-sided die; (2) unrestricted electrical conductivity of the sheet as the EHF process indirectly uses the electric energy through the fluid; (3) the sheet does not bounce back from the forming die owing to the fluid pressure. Thus, the EHF has garnered considerable attention in comparison to alternative high-speed forming processes (e.g., EMF and EF).
Most studies related to EHF have presented experimental and numerical simulation results [14][15][16][17][18][19][20]. The EHF process is conducted in a closed space for safety as the sheet is deformed at extremely high speeds > 100 m/s. Therefore, the experimental process cannot be directly observed, indicating that a numerical approach should be conducted to understand the EHF process more appropriately.
The most vital factor of the finite element analysis is regarding the reliable material properties for each numerical component to ensure adequate correspondence between the simulation and experimental results. In EHF, the material properties of a deformed sheet are essential for constructing a reliable numerical model. Generally, metallic materials such as aluminum alloys and steels exhibit strain-rate sensitivity, implying that the sheet characteristics vary with the strain rate. Therefore, the method for obtaining the material properties differs depending on the strain rate conditions. For instance, the tensile material properties of the sheet are simply obtained using dog-bone-shaped specimens on a tensile test machine at a high strain rate of approximately 10 2 -10 4 s −1 under quasistatic conditions, as specified in the ASTM standard [21]. Moreover, the split Hopkinson pressure bar test is a well-known experiment for obtaining dynamic material properties [13].  In the EHF process, the flow curves of the sheet should be obtained under as many conditions as possible prior to applying them in numerical simulations because the strain rate of the deformed sheet varies from quasistatic conditions to high-value conditions. However, the experimental method is inadequate in terms of time and cost because several experiments are required for each strain-rate condition.
Therefore, this study proposes a numerical approach for predicting the material properties of a deformed sheet in the EHF process based on the inverse parameter estimation technique, such that the final deformation shape of the sheet in the numerical simulation is almost identical to that in the experiment under the same conditions. A flowchart of the overall process is illustrated in Fig. 2.
The remainder of this paper is organized as follows. The equations required to perform the inverse parameter estimation are summarized in Sect. 2. The finite element analysis and reduced-order models required for developing the forward model of the EHF process are described in Sect. 3. Subsequently, the Markov chain Monte Carlo (MCMC) sampling technique applied to inversely estimate the parameters is described in Sect. 4. Lastly, the results of the MCMC-based inverse parameter estimation are discussed in Sect. 5, including the future scope of this research.

Constitutive equation for a sheet
In general, most metallic material sheets (e.g., aluminum alloys and steels) exhibit various characteristics depending on the deformation environment, especially at various temperatures and strain rates. Unlike in quasistatic forming processes, material sheets in high-speed forming processes (e.g., EHF) undergo deformation at extremely high strain rates greater than 1000 s −1 under various strain-rate conditions (e.g., deep drawing). Therefore, the constitutive equation considering strain-rate sensitivity should be applied for the sheet material. In context, constitutive equations have been established for high-speed forming simulations such as the Cowper-Symonds model [22], Zerilli-Armstrong model [23], Steinberg model [24], and Johnson-Cook model [25]. In this study, the Cowper-Symonds constitutive equation was adopted for the Al 5052-H32 sheet. The Cowper-Symonds model displays appropriate agreement with the experimental results at various strain-rate conditions, especially for aluminum and titanium materials [26][27][28]. The equation can be expressed as follows: where σ denotes the true stress; σ 0 represents the quasistatic true stress; ε indicates the true strain; A, B, and n denote the material constants;ε represents the strain rate; and C and p denote the strain-rate parameters. The value of σ 0 for Al 5052-H32 was obtained from a quasistatic tensile test as follows: A = 188 MPa, B = 327 MPa, and n = 0.58. The strain-rate sensitivity parameters C and p are unknowns that should be estimated using inverse parameter estimation.

Bayesian statistics based on conditional probability
The inverse parameter estimation was based on Bayesian statistics, which represents the probability of an event as a conditional probability. For a given measurement condition y ∈ R m in Bayesian statistics, the posterior probability distribution π (θ | y) of a parameter set θ ∈ R n is proportional to the product of a prior distribution p(θ ) and a likelihood l(y|θ), as described in Eq. (2).
Bayesian statistical inference involves parameter uncertainties; hence, a prior distribution p(θ ) is determined by the subjective belief of the analyst. However, it lacks objectivity, yet it can be a useful analytical tool if the population does not exist or is uncertain [29]. A likelihood l(y|θ) is the probability density of a parameter set θ for fixed measurement data y.
The inverse parameter estimation based on Eq. (2) proceeds through the following three steps: (1) Determination of the likelihood l(y|θ).
The posterior distribution of parameter θ is not expressed as a specific value but represented as a probability, including uncertainties. In this study, the desired values of θ constitute the material parameters in the Cowper-Symonds constitutive equation, and the measurement data y denote the deformed shapes of the sheet at 6, 7, and 8 kV in the experiment.

Measurement error in EHF experiment
Prior to the likelihood construction, the EHF experiment was performed under various input voltage conditions. The EHF experimental apparatus comprising the two electrodes is illustrated in Fig. 1: a 1-mm-thick sheet, chamber, free bulging die, and a capacitor bank with 32 kJ at the maximum input voltage of 15 kV.
Subsequently, input voltages of 6, 7, and 8 kV were applied, and the deformation profile at each input voltage was measured using three-dimensional (3D) scanning equipment (Fig. 3). For comparison with the 1/4 numerical simulation model, only half of the shapes presented in Fig. 3 were used and represented as a single vector (Fig. 4), which was used as a target in the inverse parameter estimation. The vector index in Fig. 3 represents the number of points for the experimental results. Overall, each experimental result had 39 points at various voltage conditions, totaling to 117 points.
The experimental equipment and optimization of the apparatus have been detailed in a previous paper [30].

Formulation for likelihood
According to Calvetti and Somersalo [29], the measurement data y are expressed as a forward model f (θ ) by mapping f :R n → R m . In this study, n denotes the number of material parameters C and p, and m represents the number of points on the deformed sheet at 6, 7, and 8 kV. In particular, n = 2 and m = 117. Thus, y can be rewritten with the measurement noise e and forward model error e f (θ ) as follows: The measurement noise e is independent of θ and acts a random variable in a probability space Ω. As θ and y are unknown, they are treated as random variables [29]. Moreover, a forward model is dependent on θ ; thus, e f (θ ) is dependent on θ as well. Based on Eq. (3), Eq. (2) can be rewritten as follows: Equation (4) signifies that the likelihood in Eq. (2) is a convolution of the measurement noise e and forward model error e f (θ ). The measurement noise in the EHF experiment contains two errors, namely, the 3D scanning noise and manual error when using the Plot Digitizer program. However, the 3D scanning noise was neglected because the precision of the 3D scanning equipment was less than 10 mm. In addition, the Plot Digitizer program was employed to manually extract the results of the bulge height at a desired point from the 3D scan result of the deformed sheet. Moreover, the error at each × point was assumed to follow a normal distribution with a mean value of 0 and a 95% confidence interval of ± 1 mm. Therefore, the bulge height variance at each point was inversely calculated as 0.5102.
Consequently, each experimental result (6, 7, and 8 kV) presented in Fig. 4 was combined into a single vector form, and the 117 points were plotted to obtain an error distribution that followed a normal distribution with the bulge height as the mean value and 0.5102 as the variance.

Prior distribution
The posterior distribution estimated herein depicts the parameter distribution, which is in good agreement with the experimental results at input voltages of 6, 7, and 8 kV. Accordingly, we divided this process into two stages: (1) estimation of the posterior distribution, which demonstrated the consistency of the experiment and the numerical simulation at a median input voltage of 7 kV, and (2) prediction of the posterior distribution based on the measurement data at 6, 7, and 8 kV.
A non-informative prior distribution was applied in the absence of prior knowledge of the parameters, which signified a uniform distribution. Therefore, the posterior distribution at 7 kV, π 7 (θ | y), can be estimated as where C lb = 10, C ub = 5000, p lb = 2, and p ub = 50.

Background of a surrogate model
A forward model is a numerical model with a function that yields the results based on the concerned parameters. Therefore, the forward model used herein is a finite element analysis model based on the strain-rate parameter θ in the Cowper-Symonds constitutive model. A forward model that can represent the deformed shape of the sheet for any θ within the defined range is required to create the likelihood, as described in Sect. 2. However, a single numerical simulation of the EHF process is time-consuming (i.e., more than 13 h). Therefore, the simulation could not be possibly conducted for a number of arbitrary θ within a predetermined range. Moreover, this study used the analysis results under three input voltage conditions (i.e., 6, 7, and 8 kV) for a single θ , which required three times as long for execution. This problem can be solved by constructing a surrogate model for the numerical simulation.
In addition, the surrogate model for the EHF numerical simulation was constructed based on 40 samples for θ , where 20 samples were extracted by Latin hypercube sampling [31] for use as the training sets. The remaining samples were selected by random sampling for use as test sets. The selected samples are presented in Fig. 5.
In this study, the deformation height of the sheet only varies in the order of a few millimeters according to the two parameters, as presented in Fig. 7. Therefore, instead of applying the classical 80/20 rule for samples, an equal number of test and training samples was used to validate the prediction accuracy of the surrogate model over the entire parameter range.

Numerical simulation in LS-DYNA
The numerical simulation was conducted using 40 samples at 6, 7, and 8 kV. The numerical model for the EHF process comprised the structural (i.e., sheet, chamber, and die) and fluid (i.e., plasma and water) components, as depicted in Fig. 6. Moreover, the EHF model dimensions were the same dimensions as the equipment used in the experimental results discussed in Sect. 2. Upon considering the apparatus characteristics, the 3D full-model is not required because it is axisymmetric in shape for the z-axis. However, a 1/4 model was applied because the 3D electrical energy required to form the sheet cannot be converted to a two-dimensional (2D) model. Subsequently, a symmetric boundary condition was applied at the x-z and y-z planes.
Although Eulerian elements are conventionally used to model the fluid parts, they cannot appropriately define the boundary of the component, which is a crucial limitation in defining the conditions of contact with the structural components [32]. Moreover, the Lagrangian element cannot be used owing to the large deformation occurring in the EHF simulation. Thus, this problem was solved using an arbitrary Lagrangian-Eulerian (ALE) element.

Numerical model for fluid components using ALE
The ALE scheme was developed to overcome the limitations of pure Lagrangian and Eulerian elements [32][33][34]. The ALE approach comprises two steps: a Lagrangian step and an advection step. In the Lagrangian step, the nodes created by the ALE formulation progressed as a Lagrangian element, and the distorted element solutions were mapped into the spatially remapped reference ALE meshes that exhibited arbitrary motion, unlike the Eulerian element that prevented the distortion of the ALE element. Moreover, ALE can manage problems with a moving boundary or a free surface, such as the fluid-structural coupling problems.
The electrical energy generated in the RLC circuit must act as an input to describe the shockwave in the water part. Therefore, a spherical plasma component was modeled with the same density as water, and the equation of state (EOS) was employed, as depicted in Fig. 6. The plasma component was assumed as an adiabatic ideal gas with a radius of 1 mm. The internal energy of the plasma can be expressed as follows: where E denotes the internal energy, p denotes the density, γ represents the adiabatic index, and V signifies the relative volume calculated by V = ρ 0 /ρ (ρ 0 denotes the initial density). The keyword "LINEAR POLYNOMIAL WITH ENERGY LEAK EO" was adopted to express Eq. (4) in LS-DYNA [35]. In Eq. (6), the internal energy E for the plasma component is described using the electrical energy and defined as the product of the current and the voltage in the RLC circuit over time. Although the current can be conveniently obtained with a Rogowski coil during the EHF experiment, the high-voltage measurement at the kV unit is dangerous and requires expensive equipment. Therefore, the voltage was assumed to drop linearly from the discharge and become zero corresponding to the maximum value of the current [36].

Structural parts for EHF simulation
All the structural parts in the EHF model were constructed with shell elements, as depicted in Fig. 6. The die and chamber were assumed as rigid bodies for time efficiency in the numerical simulation. The material properties of the sheet parts of Al 5052-H32 with 1 mm thickness were assumed to follow the Cowper-Symonds model, as described in Sect. 2.1.

Numerical simulation results for training samples
The z-axis displacements of the sheet for the 20 training samples are presented in Fig. 7. For a given training sample, the simulations were conducted at three input voltages (6, 7, and 8 kV). Moreover, each deformed shape of the sheet was combined into a single vector to describe the deformed shapes for comparison with the experimental results in Fig. 4.

Kriging formulations
As mentioned in Sect. 3.1, a surrogate model was developed to predict the free bulging deformation shape of the sheet in the EHF process with arbitrary strain-rate sensitivity parameters based on the finite element analysis model and 20 training samples. Recently, supervised learning models based on artificial neural networks (ANNs) have been applied to construct a surrogate model. However, each time a model is created with ANNs, the results of each model display a small variation. On the contrary, the Kriging model yields a deterministic result that can represent all the samples calculated by maximum likelihood estimation, so a unique surrogate model can be constructed from the given input and output samples. Therefore, we can use the Kriging approach for more appropriately estimating the outputs for arbitrary inputs as compared to the ANN surrogate model. Kriging is a Gaussian process that constructs a relationship between the input and output of the training data and can estimate the outputs at unknown inputs as a weighted linear combination of the trained outputs [37][38][39].
The Kriging model comprises two terms: a deterministic term f (x) and a stochastic term Z (x).
including the trend f (x) and a random variable . This study employed the ordinary Kriging formulation, where f (x) = α ∈ R (unknown). Moreover, as f is constant and deterministic, the covariance of Y is similar to that of Z .
The covariance between Z i and Z j can be calculated as follows: where S 2 represents the variance of Z and Ψ = Corr(x i , x j ) ∈ R n×n . In particular, the Kriging method uses the following correlation function to calculate Ψ : where θ l ∈ R n describes the hyperparameters that determine Ψ . As a high value of θ l yields correlated neighboring sample points, the output Y (x) exhibits a highly non-linear behavior. Furthermore, p l denotes a "smoothness" parameter that impacts the correlation. In this study, p l was set as 2. Equation (9) is called the Gaussian correlation function.
In particular, Z (x) follows the multivariate Gaussian distribution Z (x) ∼ N (0 n , S 2 Ψ ), where 0 n denotes a vector of n zeros. In addition, Y (x) follows the multivariate Gaussian distribution N (α1 n , S 2 Ψ ), where 1 n is a vector of n ones. Consequently, the probability distribution of the output included three unknown parameters: α, S 2 , and θ . In the original Kriging formulation, these parameters were estimated using the following log-likelihood function: Moreover, α and S 2 were evaluated based on the maximum likelihood estimation (MLE).
However, θ cannot be directly obtained from the MLE because Ψ is a function of θ . Instead, θ can be estimated using non-linear optimization using the concentrated log-likelihood function as follows: Upon determining all the parameters (i.e., α, S 2 , and θ), y(x 0 ) can be predicted at an unknown input x 0 . In the case of predicting y(x 0 ), Kriging assumes that the augmented outputỹ = [y; y(x 0 )] ∈ R n+1 exhibits the same Gaussian distribution as the original output y. Therefore,ỹ is distributed as N (α1 n+1 ,Ŝ 2Ψ ), whereΨ ∈ R (n+1)×(n+1) is the augmented correlation matrix obtained with {x 1 , x 2 , ..., x n , x 0 ∈ R n+1 }. This matrix is defined asΨ = [Ψ, ψ; ψ T , 1], where ψ is the correlation matrix between {x 1 , x 2 , ..., x n } and x 0 , Similar to the previous parameter estimation process, y(x 0 ) can be predicted at the unknown input x 0 using the MLE based on Eq. (13): In Eq. (13), Ψ −1 (y −α1 n ) was generated by training the inputs, whereas ψ was evaluated using the known training data and unknown input x 0 .
In conclusion, the surrogate model for the EHF can be generated using Eqs. (8)-(13) based on the 20 training data points.

Reduced-order model based on principal component analysis
The strain-rate sensitivity parameters C and p are the input values in the Kriging surrogate model, whereas the deformed shapes of the sheet under three distinct voltage conditions constitute the output value. However, a critical drawback of the Kriging model is that it provides only the scalar output; thus, we had to construct the Kriging model for every output. As depicted in Fig. 7, 117 deformation results were obtained for a particular parameter of the sheet. Consequently, the Kriging models must be created 117 times, requiring considerable time for model construction.
Therefore, a reduced-order model (ROM) was employed such that the order of the output could be reduced to obtain prediction results that are similar to the original results. The ROM is a mathematical model that can reduce the output dimensions while preserving the characteristics of the high-dimensional original output as much as possible. The ROM is usually employed in aerospace or fluid dynamics problems with a degree of freedom greater than 10 6 . The fundamental equation for ROM can be expressed as follows: where Z c denotes the mean-centered original data Z ∈ R m×n expressed as Z minus the sample mean Z m = 1 n n j=1 Z x;θ j ∈ R m ; w i θ j denotes the weighting factor corresponding to the basis vector v i x ∈ R m of the original data; and n denotes the number of samples. Thus, the dimension of Z can be reduced from m to s by discarding the m − s non-major basis vectors. Therefore, Eq. (14) can be rewritten as As expressed in Eq. (15), the values of the mean-centered data with small but important basis vectors could be similar to those of the original data. Thus, the basis vectors must be extracted from the original data to determine the important basis vectors, which can be performed using principal component analysis (PCA).
PCA is a mathematical technique based on linear algebra, which projects the original data in a high-dimensional space into a low-dimensional space and preserves the variance of the original data as much as possible [40][41][42]. Thus, PCA determines orthogonal basis vectors from the original data and reduces the dimension of the data into certain basis vectors based on an orthogonal projection of the original data to preserve the variance of the original data.
Let the mean-centered original deformation data be Y c (x;θ ) = [y 1 , y 2 , ..., y n ] ∈ R m×n at the x-axis coordinate x ∈ R and the material parameters θ ∈ R n . In particular, Y c was calculated by subtracting Y m = (1/n) n i=1 y i from the original data Y. Based on linear algebra, the basis vectors of Y c are the eigenvectors of the covariance matrix Σ of Y c , whereas the eigenvectors are the column vectors of Σ. The order of importance of eigenvectors is determined by arranging the eigenvalues in descending order, corresponding to each eigenvector. When the basis vectors from 2D data such as Y c (x;θ ) ∈ R m×n are extracted, a snapshot method is generally used instead of m basis vectors to obtain smaller n basis vectors for an efficient calculation [43][44][45]. Thus, the basis vectors must be extracted from the covariance matrix of Y T c , Σ snap , instead of the covariance matrix Σ of Y c . The eigenvalues and cumulative values of the 20 training samples are presented in Fig. 8. For convenience, the eigenvalues were normalized and listed in descending order. The first and maximum eigenvalues were greater than 0.991, indicating that the variance of the original data could be retained by more than 99% when the original data were orthogonally projected to the first basis vector. Moreover, the data characteristics in the new space were almost identical to those in the existing space. The absolute value for the first eigenvector in Fig. 9 exhibited almost the same geometry as the deformation shape of the numerical simulation shown in Fig. 7. Therefore, Fig. 9 implies that the EHF simulation results were represented based on the first eigenvector.
Furthermore, Eq. (15) can be rewritten as The weighting factor w 1 was calculated based on the orthogonal projection as follows: The weighting factors for the known parameter θ j can be obtained using Eq. (17). However, the weighting factor cannot be computed within the given domain of θ j , Ω = [θ lb × θ ub ] ∈ R 2 for an unknown parameter θ . Therefore, the Kriging method was employed to predict the weighting factor instead of directly estimating the sheet deformation for unknown parameters.
Based on the 20 training samples, the Kriging surrogate model was developed using the MATLAB software to estimate the sheet deformation at arbitrary material parameters θ . The actual-predicted plot for the weighting factors is illustrated in Fig. 10. The actual values were calculated using Eq. (17), whereas the predicted values were obtained using the Kriging surrogate model. Moreover, the predicted-actual line was plotted for comparison. All Furthermore, we evaluated the validity of this model for arbitrary parameters in the given range by conducting a validation using 20 test samples. The actual-predicted plot for the 20 test samples is shown in Fig. 11. Although certain test sample points were situated farther from the comparison line than the training samples, most of the test sample points were in the vicinity of the comparison line. In addition, the coefficient of determination (R 2 ) and the root mean square error (RMSE) were compared to further evaluate the prediction results in detail (Table 1).
where y i denotes the observed data,ȳ i represents the mean of y i ,ŷ i denotes the predicted value in the surrogate model, and n denotes the number of samples. A prediction model with high reliability exhibits an R 2 proximate to 1 and an RMSE almost equal to 0. According to Forrester et al. [46], in case the number of samples is more than 10, a reasonable surrogate model is required to have an RMSE of less than 0.1. Both the validation results for the training and test samples displayed appropriate prediction results. In addition, the variation in the height of the deformation at a particular point in this surrogate model is in the order of several millimeters according to the parameters, so an RMSE of less than 1, which was calculated using the error for every point, was considered as an acceptable result. In conclusion, the constructed Kriging surrogate model can appropriately predict the weighting factors within a defined range. The actual-predicted plots for the sheet deformation are presented in Fig. 12 using Eq. 8 and the predicted weighting factors for 40 samples. The numerical validation was conducted for R 2 and RMSE, and the results are listed in Table 2, which reveal that the results were reasonable. In conclusion, the reduced-order Kriging surrogate model was constructed with adequate reliability for application to the likelihood as a forward model for inverse parameter estimation.

Probability distribution for forward model error
The forward model is required to be generated and the error distribution of the model ought to be represented using equations to completely express Eq. (4). The error can be conveniently calculated from the actual and predicted deformation results of the sheet for both the training and test samples at each data point. However, the surrogate model was constructed based on the training samples; thus, the error distribution could not be reasonably determined using the training samples. Only 20 test samples were used to determine the error distribution of the forward model.
The error was calculated at each x point for all the test samples. Moreover, the normality was verified using error results at each point. In addition, a chi-square test was performed to validate the normality of the error distribution. The chi-square goodness-of-fit test is a statistical analysis based on the chi-square distribution that compares the frequency predicted from a specific distribution with the actually observed frequency. Therefore, it is a useful fit test for discrete data rather than continuous data, as applicable in this study.  Consequently, the chi-square normality test verified that all the 117 points followed a normal distribution. The mean and variance of each point are presented in Fig. 13. A number of bulge heights from the surrogate model varied from the actual values owing to the application of the ROM. Therefore, the mean value was not zero at most points.
Both errors (i.e., measurement and surrogate model errors) should be convoluted based on Eqs. (4) to complete the likelihood construction. Both error distributions follow a normal distribution. Therefore, the likelihood follows a normal distribution.
As discussed in Sect. 2.4, the posterior distribution for 7 kV was used as the prior distribution of the posterior for 6, 7, and 8 kV. The posterior distribution for 7 kV is expressed as follows: where the numbers 40 and 78 denote the 40th to 78th points displaying the deformation height at 7 kV among all the 117 points representing the sheet deformation height. Based on Eq. (20), the final posterior distribution for all input voltage conditions can be expressed as π θ | y ∝ π y | θ · π θ = π y | θ · π θ | y 7kV 4 Inverse parameter estimation using MCMC

Background of MCMC
A Markov chain is a stochastic process in which a future event is explained only by the present state and not any of its past states. In addition, the Monte Carlo method is an algorithm that stochastically predicts the parameters using random variables. The combined random-sampling process, called Markov chain Monte Carlo (MCMC), can extract random samples from a desired distribution based on a Markov chain. The samples obtained from MCMC sampling cannot be completely independent because of the Markov chain characteristics. Therefore, the samples are associated with at least one previous point in the chain [29]. Accordingly, several researchers have suggested using every nth sample from the original samples to minimize sample dependency [47]. The Metropolis-Hastings (MH) [48] and Gibbs [49] algorithms are widely used for numerous statistical or stochastic problems. The Gibbs algorithm requires a conditional distribution for each parameter. Thus, the complexity of the algorithm increases as the number of parameters to be predicted increases, and the algorithm cannot be used if representing the conditional distribution as an equation is challenging. In this study, the Gibbs algorithm could not be used because the conditional distribution of the two parameters is unknown. Therefore, the MH algorithm was employed in this study. If the distribution is f (x) at the location of the desired extraction, the MH algorithm contains two simple steps that are repeated.
(1) Given x i , create new sample y using a candidate generating kernel q(y|x i ).
(2) Select x i+1 using following this procedure: where the acceptance probability ρ (x i , y) = min f (y)q(x i |y) f (x i )q(y|x i ) , 1 . (3) Repeat steps 1 to 2 for n times.
Moreover, the white-noise random-walk function was employed as the candidate generating the kernel q.
The white-noise random-walk proposal was symmetric; thus, ρ was simplified as follows: In Eq. (22), r denotes a random-walk step size that significantly correlates with the overall acceptance ratio of the MH algorithm. If r is excessively small, the MH process exhibits a high acceptance ratio, and vice versa. Based on several prior studies, the optimal acceptance ratio in MCMC sampling is known as 23.4-44% [50][51][52]. Therefore, the random-walk step size of each parameter was adjusted to satisfy the optimum ratio. The posterior distribution of the material parameters that displayed an appropriate agreement between the numerical simulation and the experimental results (only for 7 kV) was used as the prior distribution of the posterior distribution at 6, 7, and 8 kV. Therefore, the samples following the posterior distribution at 7 kV were acquired prior to using the MH algorithm. The autocorrelation was calculated, as depicted in Fig. 14, to validate the independence of the samples. The samples were assumed as independent when the autocorrelation was within 0 ± 0.02. As can be observed from Fig. 14, both the autocorrelations were almost zero at lag = 20 for the two parameters and were within the defined range. Therefore, we obtained 10 5 samples by extracting every 20th sample and setting the burn-in period as 10 4 . The histograms of the obtained samples are presented in Fig. 15. 4.3 MCMC sampling at 6, 7, and 8 kV

Kernel density estimation
The distribution can be expressed as a function to employ the posterior distribution at 7 kV as a prior distribution for 6, 7, and 8 kV. The samples shown in Fig. 15 were similar to the normal distribution, but the chi-square goodness- of-fit test demonstrated that the samples did not follow a normal distribution. Therefore, kernel density estimation (KDE) was conducted to represent the sample histograms as functions of the probability distribution. The kernel function is a non-negative symmetrical function with an integral value of 1. Typically, kernel functions are Gaussian and uniform distributions. For each observed data, we generated a kernel function centered on the data value, added all the created kernel functions, and divided the result into the total number of data to estimate the probability density. Discrete data (e.g., histograms) can be represented as a continuous distribution in the form of a probability density function by employing KDE. Figure 16 shows the estimated probability density for parameters C and p.

MCMC sampling at 6, 7, and 8 kV
The final posterior distribution of the material parameters satisfying 6, 7, and 8 kV was obtained based on the probability density at 7 kV. As discussed in Sect. 4.2, the autocorrelation was evaluated from the first sampling (Fig. 17), and the lag was determined to be 30 with the autocorrelation existing within a given range. The acceptance ratio was approximately 29%. The sample history and histogram for each parameter are presented in Figs. 18 and 19, respectively. Figure 18 illustrates that both the sample histories appear as a white noise signal, which is completely uncorrelated. In conclusion, the obtained samples were almost independent. The normalized mean values for each parameter were selected to select a parameter set from the sample distribution for application in the numerical simulation, as follows:C = 9.2675E−01 andp = 4.8475E−01. In case these values were converted to the original values,C = 4.6345E+03 andp = 2.3060E+01.

Validation of estimated parameters
The mean values from the MCMC sampling were employed as the optimal material parameters in the Cowper-Symonds model for Al 5052-H32. Moreover, a numerical simulation was conducted using these parameters. The bulge heights for the experiment, numerical simulation, and reduced-order surrogate model with optimal parameters are comparatively presented in Fig. 20, wherein all the three cases displayed a similar deformation geometry under the three input voltage conditions. In addition, a numerical validation was performed for a more detailed investigation, as presented in Table 3.
The comparison between the numerical simulation and the surrogate model displayed the best agreement for both the error values, implying that the surrogate model based on the reduced-order and ordinary Kriging models can estimate the bulge height of the sheet. Furthermore, the other cases demonstrated reasonable error values, indicating   that the material parameters for the Cowper-Symonds model obtained from the inverse parameter estimation were reliable for predicting the bulge height in the EHF process at 6, 7, and 8 kV.

Conclusions
In this study, the Al 5052-H32 sheet material parameters in the Cowper-Symonds equation were characterized based on inverse parameter estimation, such that the numerical simulation for EHF displayed adequate correspondence with the experimental results. The comparison of the experiment and the analysis for only one experimental condition (i.e., one voltage condition) may not be valid under other conditions. Therefore, the material parameters that minimized the errors between the experiment and numerical simulation were predicted under three voltage conditions. The experiment was performed using an EHF apparatus at 6, 7, and 8 kV. The bulge height of the deformed sheet was measured and compared with that of the finite element analysis.
A finite element model with the same geometry as the experimental apparatus was constructed using the LS-DYNA program. The fluid components, plasma, air, and water were modeled based on an ALE element. In addition, the EOS with an energy term was employed to create a shockwave in the fluid components.
A reduced-order surrogate model was constructed based on an ordinary Kriging method and the PCA technique.
In developing the prior distribution at the 7 kV input voltage, a uniform distribution was selected for each parameter considering the absence of information. Accordingly, the posterior distribution at 7 kV was defined as the prior distribution for 6, 7, and 8 kV. As the equation was complex, we could not directly identify the probability distribution of the parameters. Moreover, the MCMC sampling was applied such that the distribution characteristics could be predicted by inversely acquiring the samples that follow the posterior distribution. Furthermore, the MH algorithm was employed, and the samples were obtained based on the random-walk proposal distribution. The material parameters of the Cowper-Symonds model were predicted using MCMC sampling asC = 4.6345E+03 andp = 2.3060E+01.
The cases (i.e., experiment, numerical simulation, and reduced-order surrogate model) were compared to validate the reliability of the estimated parameters. Consequently, the error values were found to be reasonable. This study has some limitations. In this study, only the final deformation shapes were employed as the objective, and the strain rate range of the parameter C was set to as low as 5000, considering the characteristics of the EHF equipment used. The parameter distributions satisfying the condition may exist outside the range specified. Therefore, the predicted material parameters may not be effective at higher strain-rate conditions. In future studies, more reliable material parameters must be predicted using a deformation shape over time instead of a final deformation shape for inverse parameter estimation.