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

Authors: Mina Woo; Kyunghoon Lee; Woojin Song; Beomsoo Kang; Jeong Kim
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–3]. These 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 high-speed conditions, the formability is improved by the die-sheet 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–7].The 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.Figure 1
Typical configuration of the electrohydraulic forming process. (a) Schematic of electrohydraulic forming. (b) Electrical circuit of electrohydraulic forming.
(a)(b)As shown in Figure1(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. The 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. Therefore, 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. The 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 102–104 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–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. The numerical model estimates the parameters that minimize the error between the numerical and experimental results.The target material was Al 6061-T6, which is used in many industrial fields, such as vessels, airplanes, and automobiles. The material constitutive equation was the Cowper–Symonds model with two unknown parameters to be estimated. The 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.Therefore, 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.The rest of this paper is structured as follows. Section2 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. The 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.


## 2. Formulation

### 2.1. 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 106 which requires hundreds of hours to complete the simulation. Therefore, the method of reducing the dimension of the model while keeping the characteristics may be the best way to reduce the computational time. This process can be performed by principal component analysis (PCA). PCA is a mathematical method for reducing high-dimensional 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].Letycx;θ∈ℝm be mean-centered simulation outputs at x∈ℝ with respect to inputs θ∈ℝp. In the ROM, the vector yc can be spanned by bases vixi=1r∈ℝm×r, where r is the number of bases. Therefore, we can present yc as a linear combination of vi∈ℝm as follows:(1)ycx;θ=∑i=1rwiθvix,where wiθ is the weighting coefficient of basis vector vi∈ℝm, which can be calculated by orthogonally projecting yc onto a basis vector vi.We can reduce the order ofyc by dropping the r−q insignificant basis vectors. To this end, we first use the significant q basis vector vix and approximate yc on a low-dimensional space as follows:(2)ycx;θ≈∑i=1qwiθvix.The basis vectors for the ROM are computed by the PCA technique, which is widely adopted for dimensionality reduction and feature extraction.LetY=y1,y2,…,yn∈ℝm×n be the total output and Yc be the mean-centered Y calculated by Y−Ymean, where Ymean=1/n∑i=1nyi.The dataZ, created by projecting the original mean-centered data Yc into the lower-dimensional space, can be expressed as a linear combination of Yc, namely, as Z=vTYc. To maximize the variance of Z, the following equation is solved by linear algebra:(3)maxvarZ=maxvarvTYc=maxvarvT∑v,where Σ is the covariance matrix of Yc. Applying the Lagrange multiplier method under the constraint v=1, we finally obtain(4)Σ−λv=0.In basic linear algebra [23], v and λ are defined as the eigenvector and eigenvalue of Σ, respectively. The eigenvectors are the principal components of v that maximize the variance of Yc. Therefore, the basis vectors for the ROM are the eigenvectors of the covariance matrix of the original data Yc.


### 2.2. Surrogate Model Based on the Kriging Method

The ROM process calculates the weighting coefficientwiθ in equation (1) by orthogonally projecting yc onto vi. The weighting coefficients are calculated by the vector projection formula as follows:(5)y¯c=yc⋅vivivivi=wiθjvi⟶wiθj=yc⋅vi,where y¯c is the vector obtained by projecting yc onto the basis vi.For known input parametersθj∈ℝp, such as the training and test samples in this study, the basis weighting coefficients wiθj are easily obtained using equation (5). However, when the inputs within a given domain Ω=θlb×θub∈ℝp are unknown, as in the present case, wiθj cannot be computed. Therefore, they were estimated by the ordinary Kriging method.The 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]. The 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.The Kriging model consists of two terms, a deterministic trendfx and a stochastic term Zx:(6)Yx=fx+Zx.TheYx in equation (6) can be assigned as the weighting coefficients in the ROM.In the Kriging formulation, a set of observationsyxii=1n is represented by a random variable Yxii=1n composed of the trend fx and a random variable zxii=1n.Depending on their definition offx, Kriging formulations are classified as simple when fx=α (known), ordinary when fx=α (unknown), and universal when fx=xTβ. This study adopts the ordinary Kriging formulation, in which the deterministic trend is unknown and constant, that is, fx=α∈ℝ.Becauseα is constant and deterministic, the covariance of Yxii=1n is that of Zxii=1n.The covariance betweenZi and Zj is determined by correlating xi and xj as follows:(7)CovZi,Zj=σ2Ψ,where Ψ  =  Corrxi,xj∈ℝn×n.To calculateΨ, the Kriging method uses the following Gaussian correlation function:(8)Ψ=exp−∑l=1kθlxi.l−xj.l2.Here,θl∈ℝp is a hyperparameter that determines the correlation strength between xi and xj. A high θl means that only the neighboring sample points are correlated, so the output Yx exhibits a highly nonlinear behavior.In the Kriging method,Zx also varies locally and follows the multivariate Gaussian distribution Zx∼N0n,σ2Ψ, where 0n is a vector of n zeros. Consequently, Yxii=1n follows the multivariate Gaussian distribution Nα1n,σ2Ψ, where 1n is a vector of n ones. The 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:(9)lnLα,σ2,θ=−n2ln2π+lnσ2+lnΨ−12σ2y−α1nTΨ−1y−α1n.Applying the maximum likelihood estimation (MLE),α and σ2 are, respectively, determined as follows:(10)α^=1nTΨ−1y1nTΨ−11n,σ^2=y−α^1nTΨ−1y−α^1nn.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:(11)lnLcθ=−n2lnσ^2−12lnΨ.Once all parametersα, σ2, and θ are estimated, we can predict yx0 at some unknown input x0. When predicting yx0, Kriging assumes that the augmented output y˜=y;yx0∈ℝn+1 has the same Gaussian distribution as the original output y. Therefore, y˜ is distributed as Nα^1n+1,σ^2Ψ˜, where Ψ˜∈ℝn+1×n+1 is the augmented correlation matrix obtained with x1,x2,…,xn,x0∈ℝn+1. This matrix is defined as Ψ˜=Ψ,ψ;ψT,1, where ψ=ψx1,x0,…ψxn,x0.Similarly to the previous parameter estimation process,yx0 at the unknown input x0 is finally predicted by MLE using the following equation:(12)y^0=α^+ψTΨ−1y−α^1n.In this expression,Ψ−1y−α^1n is generated by training the inputs, and ψ is calculated using the known training data and the unknown input x0.


## 2.1. 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 106 which requires hundreds of hours to complete the simulation. Therefore, the method of reducing the dimension of the model while keeping the characteristics may be the best way to reduce the computational time. This process can be performed by principal component analysis (PCA). PCA is a mathematical method for reducing high-dimensional 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].Letycx;θ∈ℝm be mean-centered simulation outputs at x∈ℝ with respect to inputs θ∈ℝp. In the ROM, the vector yc can be spanned by bases vixi=1r∈ℝm×r, where r is the number of bases. Therefore, we can present yc as a linear combination of vi∈ℝm as follows:(1)ycx;θ=∑i=1rwiθvix,where wiθ is the weighting coefficient of basis vector vi∈ℝm, which can be calculated by orthogonally projecting yc onto a basis vector vi.We can reduce the order ofyc by dropping the r−q insignificant basis vectors. To this end, we first use the significant q basis vector vix and approximate yc on a low-dimensional space as follows:(2)ycx;θ≈∑i=1qwiθvix.The basis vectors for the ROM are computed by the PCA technique, which is widely adopted for dimensionality reduction and feature extraction.LetY=y1,y2,…,yn∈ℝm×n be the total output and Yc be the mean-centered Y calculated by Y−Ymean, where Ymean=1/n∑i=1nyi.The dataZ, created by projecting the original mean-centered data Yc into the lower-dimensional space, can be expressed as a linear combination of Yc, namely, as Z=vTYc. To maximize the variance of Z, the following equation is solved by linear algebra:(3)maxvarZ=maxvarvTYc=maxvarvT∑v,where Σ is the covariance matrix of Yc. Applying the Lagrange multiplier method under the constraint v=1, we finally obtain(4)Σ−λv=0.In basic linear algebra [23], v and λ are defined as the eigenvector and eigenvalue of Σ, respectively. The eigenvectors are the principal components of v that maximize the variance of Yc. Therefore, the basis vectors for the ROM are the eigenvectors of the covariance matrix of the original data Yc.


## 2.2. Surrogate Model Based on the Kriging Method

The ROM process calculates the weighting coefficientwiθ in equation (1) by orthogonally projecting yc onto vi. The weighting coefficients are calculated by the vector projection formula as follows:(5)y¯c=yc⋅vivivivi=wiθjvi⟶wiθj=yc⋅vi,where y¯c is the vector obtained by projecting yc onto the basis vi.For known input parametersθj∈ℝp, such as the training and test samples in this study, the basis weighting coefficients wiθj are easily obtained using equation (5). However, when the inputs within a given domain Ω=θlb×θub∈ℝp are unknown, as in the present case, wiθj cannot be computed. Therefore, they were estimated by the ordinary Kriging method.The 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]. The 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.The Kriging model consists of two terms, a deterministic trendfx and a stochastic term Zx:(6)Yx=fx+Zx.TheYx in equation (6) can be assigned as the weighting coefficients in the ROM.In the Kriging formulation, a set of observationsyxii=1n is represented by a random variable Yxii=1n composed of the trend fx and a random variable zxii=1n.Depending on their definition offx, Kriging formulations are classified as simple when fx=α (known), ordinary when fx=α (unknown), and universal when fx=xTβ. This study adopts the ordinary Kriging formulation, in which the deterministic trend is unknown and constant, that is, fx=α∈ℝ.Becauseα is constant and deterministic, the covariance of Yxii=1n is that of Zxii=1n.The covariance betweenZi and Zj is determined by correlating xi and xj as follows:(7)CovZi,Zj=σ2Ψ,where Ψ  =  Corrxi,xj∈ℝn×n.To calculateΨ, the Kriging method uses the following Gaussian correlation function:(8)Ψ=exp−∑l=1kθlxi.l−xj.l2.Here,θl∈ℝp is a hyperparameter that determines the correlation strength between xi and xj. A high θl means that only the neighboring sample points are correlated, so the output Yx exhibits a highly nonlinear behavior.In the Kriging method,Zx also varies locally and follows the multivariate Gaussian distribution Zx∼N0n,σ2Ψ, where 0n is a vector of n zeros. Consequently, Yxii=1n follows the multivariate Gaussian distribution Nα1n,σ2Ψ, where 1n is a vector of n ones. The 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:(9)lnLα,σ2,θ=−n2ln2π+lnσ2+lnΨ−12σ2y−α1nTΨ−1y−α1n.Applying the maximum likelihood estimation (MLE),α and σ2 are, respectively, determined as follows:(10)α^=1nTΨ−1y1nTΨ−11n,σ^2=y−α^1nTΨ−1y−α^1nn.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:(11)lnLcθ=−n2lnσ^2−12lnΨ.Once all parametersα, σ2, and θ are estimated, we can predict yx0 at some unknown input x0. When predicting yx0, Kriging assumes that the augmented output y˜=y;yx0∈ℝn+1 has the same Gaussian distribution as the original output y. Therefore, y˜ is distributed as Nα^1n+1,σ^2Ψ˜, where Ψ˜∈ℝn+1×n+1 is the augmented correlation matrix obtained with x1,x2,…,xn,x0∈ℝn+1. This matrix is defined as Ψ˜=Ψ,ψ;ψT,1, where ψ=ψx1,x0,…ψxn,x0.Similarly to the previous parameter estimation process,yx0 at the unknown input x0 is finally predicted by MLE using the following equation:(12)y^0=α^+ψTΨ−1y−α^1n.In this expression,Ψ−1y−α^1n is generated by training the inputs, and ψ is calculated using the known training data and the unknown input x0.


## 3. Electrohydraulic Forming Test Using a Free-Bulging Die

### 3.1. Experiment

The 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]. Therefore, we applied a free-bulging die to avoid a contact between a sheet and a die.The experiment was configured in the RLC circuit presented in Figure1(b). In this circuit, Rt is the total resistance without the wire, Rw 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. When the shockwaves arrive at the sheet, they deform it into a die shape. The 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 were250×250 mm. The sheet thickness was 1 mm. The chemical composition of Al 6061-T6 is presented in Table 1.Table 1
Composition of the Al 6061-T6 sheet used in the present experiments (unit: %).
MaterialAlSiFeCuMnMgCrZnTiRemainderValue97.370.680.410.220.130.910.160.050.040.03When 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 Figure2 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.Figure 2
Current curve in the electrohydraulic forming test at 8 kV.The reproducibility of the EHF experiment was confirmed in a previous study [25], so only one experimental result is presented here. The 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.


### 3.2. 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]. The fluid parts (plasma and water) were modeled as ALE materials, and the structural parts (sheet, chamber, and die) were modeled as Lagrangian elements. The 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. Therefore, we established a ¼ model and applied symmetric conditions in the x-z and y-z planes.Figure 3
One-quarter numerical model of electrohydraulic forming.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. The 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].The 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.The high-speed deformation of the sheet under the high-strain-rate effect was modeled by the following Cowper–Symonds constitutive equation:(13)σeff=A+Bεeffn1+ε˙effC1/p,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. The 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 θ∈ℝ2 to be found by inverse estimation in the Kriging surrogate model. Before constructing the surrogate model, we created the training and test samples for the reduced order EHF simulation by Latin hypercube and random sampling, respectively (see Figure 4). As the training samples characterize the surrogate model in a defined domain, Latin hypercube sampling is the best choice because it distributes the samples evenly within the domain [27]. Setting the specific domain as ΩθC,p=1000,2×12000,20∈ℝ2, we selected 20 training and 20 test samples. Training samples were used to construct the surrogate model based on the Kriging method, and test samples were employed to determine whether the surrogate model based on the training data can predict well for any point in a given domain. After numerical simulations on the 40 samples, the final z-axis displacement of the sheet was extracted as presented in Figure 5. The legend in Figure 5 represents the number of the training and test samples, respectively.Figure 4
Samples for constructing the surrogate model. (a) Training samples. (b) Test samples.
(a)(b)Figure 5
Results of the (a) training samples and (b) test samples in the numerical simulation.
(a)(b)From thez-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∈ℝ39×20 by the snapshot method. The cumulative sum of the normalized eigenvalues is presented in Figure 6. The 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 significant vector representation of the training samples. Therefore, 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(14)ycx;θ≈w1θv1x.Figure 6
Normalized eigenvalues and cumulative sum of eigenvalues.Figure 7
Profile of the basis vector of the numerical simulation output.The weighting coefficientw1θ was calculated from the orthogonal projections of the 20 training samples. The 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.Figure8 shows a flowchart of the Kriging- and ROM-based surrogate model.Figure 8
Flowchart for constructing a surrogate model based on ROM and Kriging.


## 3.1. Experiment

The 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]. Therefore, we applied a free-bulging die to avoid a contact between a sheet and a die.The experiment was configured in the RLC circuit presented in Figure1(b). In this circuit, Rt is the total resistance without the wire, Rw 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. When the shockwaves arrive at the sheet, they deform it into a die shape. The 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 were250×250 mm. The sheet thickness was 1 mm. The chemical composition of Al 6061-T6 is presented in Table 1.Table 1
Composition of the Al 6061-T6 sheet used in the present experiments (unit: %).
MaterialAlSiFeCuMnMgCrZnTiRemainderValue97.370.680.410.220.130.910.160.050.040.03When 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 Figure2 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.Figure 2
Current curve in the electrohydraulic forming test at 8 kV.The reproducibility of the EHF experiment was confirmed in a previous study [25], so only one experimental result is presented here. The 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.


## 3.2. 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]. The fluid parts (plasma and water) were modeled as ALE materials, and the structural parts (sheet, chamber, and die) were modeled as Lagrangian elements. The 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. Therefore, we established a ¼ model and applied symmetric conditions in the x-z and y-z planes.Figure 3
One-quarter numerical model of electrohydraulic forming.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. The 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].The 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.The high-speed deformation of the sheet under the high-strain-rate effect was modeled by the following Cowper–Symonds constitutive equation:(13)σeff=A+Bεeffn1+ε˙effC1/p,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. The 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 θ∈ℝ2 to be found by inverse estimation in the Kriging surrogate model. Before constructing the surrogate model, we created the training and test samples for the reduced order EHF simulation by Latin hypercube and random sampling, respectively (see Figure 4). As the training samples characterize the surrogate model in a defined domain, Latin hypercube sampling is the best choice because it distributes the samples evenly within the domain [27]. Setting the specific domain as ΩθC,p=1000,2×12000,20∈ℝ2, we selected 20 training and 20 test samples. Training samples were used to construct the surrogate model based on the Kriging method, and test samples were employed to determine whether the surrogate model based on the training data can predict well for any point in a given domain. After numerical simulations on the 40 samples, the final z-axis displacement of the sheet was extracted as presented in Figure 5. The legend in Figure 5 represents the number of the training and test samples, respectively.Figure 4
Samples for constructing the surrogate model. (a) Training samples. (b) Test samples.
(a)(b)Figure 5
Results of the (a) training samples and (b) test samples in the numerical simulation.
(a)(b)From thez-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∈ℝ39×20 by the snapshot method. The cumulative sum of the normalized eigenvalues is presented in Figure 6. The 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 significant vector representation of the training samples. Therefore, 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(14)ycx;θ≈w1θv1x.Figure 6
Normalized eigenvalues and cumulative sum of eigenvalues.Figure 7
Profile of the basis vector of the numerical simulation output.The weighting coefficientw1θ was calculated from the orthogonal projections of the 20 training samples. The 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.Figure8 shows a flowchart of the Kriging- and ROM-based surrogate model.Figure 8
Flowchart for constructing a surrogate model based on ROM and Kriging.


## 4. Surrogate Model Verification

The Kriging surrogate model for the EHF simulation was constructed in MATLAB 2019b, and the estimated hyperparameters are presented in Table2.Table 2
Estimated hyperparameters of the Kriging model.
Hyperparametersθ^Cθ^pα^σ^2Values100.8182100.04010.43340.1342As the hyperparameterθ^ is more influential at larger values, the material parameter C affects the weighting coefficient more than p [28].The predicted weighting coefficients are plotted against the actual values in Figure9. The 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. The coefficients of determination (R2) were 1.0 for the training samples and 0.77012 for the test samples. R2 is typically between zero and one, with high values indicating good predictability of the surrogate model. Although the test samples achieved a lower R2 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.Figure 9
Actual-predicted plots of weighting coefficients for the (a) training samples and (b) test samples.
(a)(b)Using the weighting coefficients obtained from the surrogate model, we verified thez-axis displacement yi of the deformed sheet calculated using equation (15). The accuracy of yi was assessed by two error values: R2 and the root mean squared error (RMSE):(15)yi=wiv1+yc.When the RMSE value is close to zero, a mathematical model makes accurate predictions. Table3 presents the numerical verification results on the training and test datasets, and Figure 10 shows a graphical representation of these results. All R2 values exceeded 0.998 (maximum 0.99997) on the training samples and 0.992 (maximum 0. 99982) on the test samples. The 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.Table 3
Numerical verification results ofz-axis displacements.
Training samplesTest samplesR2RMSE (mm)R2RMSE (mm)19.9948E − 012.3659E − 019.9806E − 014.2400E − 0129.9924E − 012.6144E − 019.9965E − 011.9028E − 0139.9985E − 011.3088E − 019.9664E − 015.8167E − 0149.9983E − 011.2709E − 019.9858E − 013.6856E − 0159.9986E − 011.1556E − 019.9811E − 014.1098E − 0169.9985E − 011.1969E − 019.9601E − 016.0066E − 0179.9997E − 015.2793E − 029.9856E − 013.7306E − 0189.9989E − 011.0223E − 019.9260E − 019.2354E − 0199.9996E − 016.4213E − 029.9598E − 016.7832E − 01109.9990E − 011.1268E − 019.9883E − 013.5514E − 01119.9997E − 016.0162E − 029.9930E − 012.5394E − 01129.9994E − 017.6831E − 029.9969E − 011.8252E − 01139.9976E − 011.5338E − 019.9634E − 015.8759E − 01149.9962E − 011.9055E − 019.9952E − 012.2083E − 01159.9978E − 011.4466E − 019.9539E − 016.4169E − 01169.9990E − 011.0184E − 019.9975E − 011.5894E − 01179.9994E − 017.9496E − 029.9480E − 017.7122E − 01189.9993E − 019.5267E − 029.9840E − 014.1125E − 01199.9919E − 012.7595E − 019.9982E − 011.3214E − 01209.9994E − 017.2983E − 029.9797E − 014.4194E − 01Min9.9919E − 015.2793E − 029.9260E − 011.3214E − 01Max9.9997E − 012.7595E − 019.9982E − 019.2354E − 01Avg9.9979E − 011.2871E − 019.9770E − 014.3541E − 01Figure 10
Graphical results of the numerical verification. (a) Coefficient of determination (R2). (b) Root mean squared error (RMSE).
(a)(b)Observing the predicted weighting coefficients in Figure9(b), the R2 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. The 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.Figure 11
Normalized weighting coefficient as a function of normalized parameterC.


## 5. Optimizing the Material Parameters

In the constructed surrogate model, we found the optimum material parametersC 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:(16)eC,p=min∑i=1myi,estC,p−yi,exp2,where yi,est is the z-axis displacement determined by the surrogate model and yi,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 yi,est determined by the material parameters C and p, our selected optimization algorithm must handle nonlinear problems. Therefore, 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. The global nature of the GA avoids trapping the solutions in local minima; moreover, the initial values are randomly generated and need not be preset.Figure12 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 minimized at 2.9405 for C = 8904.6 and p = 15.832. The acquired optimal parameters were validated in three cases: experiment, LS-DYNA, and the surrogate model. The results in the three cases are compared in Figure 13, and the error values, R2, and RMSE, are listed in Table 4.Figure 12
Two- and three-dimensional error contours in the experiment and the surrogate model. (a) 2D contour and (b) 3D contour.
(a)(b)Figure 13
Comparison ofz-axis displacements in the experiment, LS-DYNA, and the surrogate model.Table 4
Numerical verification: comparisons of the experimental, LS-DYNA, and surrogate model results.
Experiment-surrogate modelExperiment-LS-DYNALS-DYNA-surrogate modelR29.9763E − 019.9639E − 019.9828E − 01RMSE4.7086E − 015.7931E − 014.0122E − 01As presented in Table4, all errors were acceptable. In the LS-DYNA versus surrogate model, the R2 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 asC = 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 Figure14, 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.Figure 14
Comparison ofz-axis displacements in the experiment and numerical simulation with quasistatic material property.The flow curves for Al 6061-T6 are described in Figure15 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 103 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].Figure 15
Stress-strain curves for Al 6061-T6 based on Cowper–Symonds model.


## 6. 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 parametersC and p, we performed numerical simulations in LS-DYNA and extracted the eigenvectors and eigenvalues by principle component analysis. The 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. The predictability of the constructed model was evaluated by theR2 and RMSE of 20 test samples. The R2 of the weighting coefficients was approximately 0.77 for the test samples, which is lower than for the training samples. However, the R2 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.”(3)
In the constructed surrogate model, we refined the material parameters in the Cowper–Symonds equation by nonlinear optimization. For this purpose, we employed a genetic algorithm that minimizes the error between the experiment and the numerical simulation, finally obtainingC = 8904.6 and p = 15.832.However, this study has two shortcomings as follows:(1)
The 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)
The 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 Figure12, indicating that several solutions satisfy the objective function.Therefore, 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.



---

*Source: 3219829-2020-06-10.xml*