2.3. Constitutive Modeling
Accurate constitutive material models are essential for the realistic simulation of the SPF process. Many established models found in the literature are identified as phenomenological material models [
6,
7,
8,
9,
10,
16,
17,
18,
19,
20,
21]. These are viscoplastic models that represent the macroscopic behavior of deformation, usually at a constant forming temperature. They reflect mainly the relationship between the effective flow stress and the effective strain and strain rate. The constants in these models are fitted using uniaxial tensile [
10] or bulge forming experimental data [
9,
20,
21].
The superplastic behavior in many cases is represented by the simple power law equation
where (
) is the effective flow stress, (
) is the material strength coefficient, (
) is the effective strain rate, and (
) is the strain rate sensitivity index which should be greater than 0.3 for superplastic behavior to exist, and thus, the high sensitivity to strain rate.
As a first stage in establishing a material model for AZ31B, calibration of the constitutive in Equation (1) was carried out using data from the bulge forming experiments and following the procedure by Enikeev and Kruglov [
22]. For two different pressures (
) and (
), the forming times (
) and (
) are defined as the corresponding forming times that lead to the same pole height (
). According to Reference [
22] the value of (
) can be determined using the following equation:
A schematic representation of the free bulging process is shown in
Figure 6. The angle (
) is defined as half of the angle subtended by the dome surface at the center of curvature.
In this study, instead of relying only on two sets of data points (
,
), an averaged strain rate sensitivity index (
) was determined using six different values of (
), see
Table 5 and
Table 6. These correspond to six different values of the pole height (
) as calculated by the geometric relationship
where (
) is the die radius.
It is noteworthy that, in
Table 5 and
Table 6, the first experimental data corresponds to a value of (
). This is due to the high rate of change in (
) at the initial stages of the test, which limits the practicality of collecting data for (
).
From the data in
Table 2,
Table 3 and
Table 4 and using the linear interpolation method, the forming times corresponding to the six values of the pole height were calculated for the constant pressures (
) and (
). Then Equation (2) was used to generate six different values of the strain rate sensitivity index (
). At last, an average value of 0.4 from the six distinct (
) values was obtained. A summary of these calculations is found in
Table 5.
According to Reference [
22], the parameter (
) is the average of (
) and (
) as calculated from the expressions:
where (
) is the initial thickness of the sheet, which is equal to 1.2 mm. Actually, the integral on the right side of Equation (5) is not expressed through elementary functions; however, it can be easily calculated numerically.
Table 6 shows the values of (
) and (
) calculated from Equations (4) and (5). An average value of the material strength coefficient (
) was determined to be
. Notice that the model constants (
) and (
) were calculated directly with no need for any assumptions or trial and error.
Following is a relatively more sophisticated constitutive model that has been used in a number of studies [
23,
24,
25,
26] for superplastic materials. This model involves the effects of strain hardening, strain rate sensitivity, grain growth and cavitation:
where (
) is the effective strain, (
) is the initial grain size, (
) is the initial area fraction of voids, (
) is the void growth index, (
) is strain hardening exponent, (
) is a grain growth constant, and (
) is a material constant. However, in the aforementioned studies [
23,
24,
25,
26], the values of the parameters were determined by trial and error.
In order to calibrate this model for the AZ31B material, the following procedure was carried out. The value of the strain rate sensitivity index (
) was taken to be the same as that calculated for the simple model Equation (1). The values of (
) and (
) were based on the relationship between the grain size and the strain established from the microstructure analysis. The strain hardening exponent (
) was assigned a value of 0.2, similar to its value in a previous study by Li et al. [
27]. Since no experimental results were available on the volume fraction of voids, (
) and (
), the values of (
A,
and
ψ) were determined using the Taguchi design method.
Taguchi method is a powerful tool to optimize multiple parameters in experimental design that produces the minimum experimental units [
28,
29,
30]. This method was implemented in the commercially available optimization software modeFRONTIER version 2017R2 (ESTECO, AREA Science Park, Trieste, Italy).
Figure 7 shows the workflow chart for the optimization of parameters (
,
, and
). At first, a Design of Experiment (DOE) was set up, and the acceptable range for the three parameters was established.
An EasyDriver function was created with the input file, and the material Fortran-subroutine file for ABAQUS added to this EasyDriver. Parameter () was linked with the Fortran file, and EasyDriver performed the loop operations of ABAQUS version 6.13-1 (Dassault Systemes, Velizy-Villacoublay Cedex, France) simulations changing the value of () automatically.
When the loop operations were finished, the simulated output (ODB) files were imported into ABAQUS through the Transfile function, then the simulated time (Time-S), through thickness strain of the pole elements (LE22-1 and LE22-2) and vertical displacement of the same element (U2) were obtained from ABAQUS.
Time-S, LE22, and U2 were imported into the Calculator box, in which, the relationship between the simulation pole thickness (Thickness-S) and LE22 were defined, as well as for the relationship between the dome height (Height-S) and U2. Then the simulation-generated pole thickness and dome height corresponding to the experimental time (Time-E) were obtained through interpolation.
Finally, the mean square error (MSE) of the pole thickness was calculated by comparing the Thickness-S with experimental pole thickness (Thickness-E) through a (CurveFitting-1) function. The MSE of the dome height (Mean-square-error-2) was obtained by comparing the Thickness-H with the experimental dome height (Height-E) through CurveFitting-2.
In the first round of optimization, the acceptable ranges for the parameters (
), (
) and (
) were set as 100–500, 0.01–0.1, and 0–1.8, respectively. An arbitrary combination of the values of (
), (
) and (
) were set in 30 different simulated experiments.
Figure 8 shows the MSE of both the pole thickness and dome height corresponding to different values of the parameter (
) at
P = 0.15 MPa. From the figure, it can be concluded that the optimum value of (
) is roughly in the range of 230 to 320. The corresponding figures for (
) and (
) were not conclusive.
The second round of optimization was carried out with the range of parameter () reduced to 230-320. The ranges of parameters () and () were still 0.01–0.1 and 0–1.8, respectively. The step sizes for the parameters (), () and () were 10, 0.01 and 0.2, respectively. For a conventional Design of Experiment (DOE), 1000 simulations need to be performed to achieve the optimum values for the three parameters, because parameter (), () and () have 10 different values, respectively. However, the number of experiment sets was decreased to 100 through the use of the Taguchi method. The three parameters were arrayed orthogonally using the Taguchi method. A total of 100 simulations were conducted for the three different pressures (0.15 MPa, 0.2 MPa and 0.3 MPa).
Figure 9 demonstrates the influence of (
) on the MSE of the pole thickness and dome height for
P = 0.15 MPa. Both the pole thickness and dome height almost have the same trend with a minimum value of MSE at (
= 0.02). Corresponding to the results considering the pressure values of 0.2 MPa and 0.3 MPa, not shown here for briefness, a final average value of 0.01 was set for (
).
Figure 10 displays the influence of (
) on the MSE of the pole thickness and the dome height for a forming pressure of 0.15 MPa. An optimum value of (
) of 0.6 can be easily identified from the figure. Almost the same result was obtained for
P = 0.2 MPa and
P = 0.3 MPa.
Figure 11 shows the effect of (
) on the MSE of the pole thickness and dome height for
P = 0.15 MPa. For this pressure value, from the figure, the optimum value of (
) is around 280. The optimum of (
) was about 270 for
P = 0.2 MPa and in the range between 260 and 270 for
P = 0.3 MPa. The range of parameter (
) was further narrowed down with the third round of optimization.
For the last round of optimization, the values of (
) and (
) were fixed to 0.01 and 0.6, respectively. The range of parameter (
) was chosen to be between 245 and 280 for the three pressures. After running the analysis, the optimum value of parameter (
) was found to be 256. Therefore, the final obtained form of the AZ31B material model is given by