An Accurate Constitutive Model for AZ31B Magnesium Alloy during Superplastic Forming

: Accurate constitutive material models are essential for the realistic simulation of metal forming processes. However, for superplastic forming, mostly the material models found in the literature are based on fitting of the simple power law equation. In this study, an AZ31B constitutive model that takes into account microstructural evolution is introduced. This model takes into account grain growth and cavity formation in addition to strain and strain rate hardening. The model parameters were calibrated using the results of high temperature bulge forming tests and microstructural analysis. The Taguchi optimization method was used in the fitting process. In order to verify the model, simulations of the superplastic forming of two different geometries were carried out, and the results were compared with those obtained experimentally. Results show that the proposed model can accurately predict the formed geometry and thickness distribution.


Introduction
The ability of certain materials, under certain conditions, to undergo large uniform plastic strains is known as "superplasticity". This phenomenon had been applied in the manufacturing of the Damascus steels from 300 B.C. to the late nineteenth century [1][2][3]. However, it was not before the sixties of the twentieth century that it was commercially re-applied in modern industrial applications. Currently, superplastic forming is used for manufacturing lightweight complex components in niche applications; such as in the aerospace, luxury and sports cars, trains, and medical industries.
Magnesium (Mg) alloys have many advantages, such as high specific strength and specific stiffness, excellent thermal conductivity, outstanding shock absorption, and strong electromagnetic shielding. However, because of their typical hexagonal close packed crystal structure, they exhibit low ductility at Room Temperature (RT), making their sheet applications greatly limited. However, the Mg alloy AZ31B, which possess superplastic properties, has been demonstrated to exhibit improved formability at warm and high temperatures [4,5]. Therefore, it can be used in forming parts without the need for secondary joining operations. In fact, complex parts in a near-net shape can be manufactured with this alloy. Achieving its full potential will create a wide range of new applications in the electronics, transportation, aerospace and other industries. Therefore, there is a great practical significance in studying the superplastic formability of AZ31B. This is particularly important for the aerospace and space applications, in which lightweight structuring is a major requirement.
Despite the potential of this alloy, it is a known fact that, there are only a small number of available experimental and numerical studies on the forming of AZ31B. Most of the constitutive models predicting the deformation of this alloy are oversimplified, based on tensile experiments, or not validated with enough experimental results.
When incorporated in the simulation of superplastic forming (SPF), the simple material models were reported to produce unsatisfactory results [6][7][8][9][10]. The main reason many studies rely on simple material models is the difficulty in determining the material constants. Over the years, several procedures have been introduced in the literature. For constitutive models with only two material constants, the least square method is usually used in establishing the unknown parameters. However, for constitutive models involving more than two material parameters, more sophisticated tools should be used. For example, the Taguchi method is a recently new powerful tool to optimize multiple parameters in experimental design. The utilization of such a tool in constitutive modeling has been proven to provide significant benefits in material modeling.
Furthermore, uniaxial tensile data has been used in superplastic constitutive modeling in a number of studies, see for example [10]. However, forming experiments has shown that such an approach would produce inaccurate results. This is expected because the stress state in uniaxial testing does not represent the conditions during the forming of actual components. In fact, bulge forming data are better to be used when modeling the SPF process [11].
Whereas, many studies have been carried out related to the constitutive modeling of the superplastic aluminum alloy AA5083, only a few studies were concerned with the Mg alloy AZ31B. Therefore, there is a need for an accurate constitutive model for this material. In this study, a new model was proposed for AZ31B that take into account the microstructure evolution during deformation. Furthermore, a new numerical-experimental procedure was proposed and followed for establishing the material parameters. A successful AZ31B model can contribute to the literature enormously.

Materials and Methods
The superplastic Mg alloy AZ31B is known to exhibit high ductility at elevated temperatures and low strain rates. As an example of such ductility, Figure 1 shows a comparison between uniaxial AZ31B specimens tested at various temperatures. The tests were conducted at a constant cross head speed of 0.01 mm/s. The original specimen of the as-received (AR) material is shown at the top of the figure. The specimen tested at RT showed a percent elongation of only 15.3%. The specimens tested at 100 °C and 450 °C showed a percent elongation of 24.1% and 113.7%, respectively. The experiments in this study were conducted using sheets of AZ31B. The chemical constituents of the as-received AZ31B alloy used in this study are shown in Table 1. Since it resembles the initial stages in SPF of real components, high-temperature bulge forming is usually used for judging the formability of superplastic materials. In this study, an in-house developed SPF experimental setup, shown in Figure 2, was used for carrying out the bulge forming of AZ31B. The setup uses high pressure argon gas to blow the sheet specimen downwards into an open cylindrical die. The applied gas pressure is precisely controlled using the pressure control system, shown in Figures 2 and 3. It contains two lines: A low pressure gas line and a high pressure one. For the low pressure line, a manual pressure reducing regulator (Pilot Regulator) is installed to limit the maximum pressure to 0.83 MPa. Gas from this line is used as input to an electronic regulator. The pressure output of the electronic regulator is tuned via a computer signal and is used to actuate an air loaded regulator. In the high pressure line, the pressure gets controlled using the air loaded regulator. In addition, a feedback transducer was added to improve the control precision of the system. Finally, the output from the high pressure line is used to perform the actual pneumatic forming.
Before applying the gas pressure, the die and the blank holder were both heated to a temperature of 450 °C by a ceramic band heater. According to recent studies by Sorgente et al. [12] and Kaya et al. [13], this temperature represents a valuable compromise between forming times, microstructure and formability. The 1.2 mm thick AZ31B sheet was then clamped between the die and the blank holder with force enough to seal the chamber between the blank holder and the sheet. It is assumed that the sheet will reach the forming temperature after clamping it for ten minutes.  Bulge forming experiments were carried out, according to the ASTM E2712-15 standard [14], using constant gas pressures of 0.15, 0.2, and 0.3 MPa which are similar in magnitude to those considered by Franchitti et al. [15] for the same alloy. Furthermore, for each pressure, the rupture time of the AZ31B sheet was first determined by taking the average of three trials. Figure 4 shows an AZ31B sheet bulge formed up to rupture at a temperature of 450 °C and using a constant gas pressure of 0.15 MPa. Afterwards, the rupture time for each pressure level was divided into six intervals in order to determine the times at which to terminate the tests. At the end of each test, the specimen was removed and allowed to air cool to room temperature. The dome height was then recorded, and the specimen was cut into two halves for measuring the pole thickness. It is noteworthy that for each constant pressure and at each forming time three different trials were conducted. The complete results of the bulge forming experiments are shown in Tables 2-4.

Microstructure Analysis
Samples were cut from the bulge formed specimens around the pole area. The samples were then prepared for optical microscopy using a 5 megapixel digital camera microscope (OLYMPUS DP27, MEA FZ-LLC, Dubai, United Arab Emirates). This was accomplished by a sequence of grinding and polishing steps followed by the application of acetic-picral as an etchant. Considering samples from forming under different pressures up to different times, a relationship was established between the grain size and the strain, as shown in Figure 5. When extrapolated, the linear relationship crosses the y-axis at a grain size of 10.3 micrometers. This value agrees with the typical required initial grain size for superplastic materials.

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 Tables 5-6. These correspond to six different values of the pole height (ℎ) as calculated by the geometric relationship where ( = 51 mm) is the die radius. It is noteworthy that, in Tables 5-6, the first experimental data corresponds to a value of ( ≅ 53.1°). 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 ( < 53.1°).
From the data in Tables 2-4 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 100 MPa · sec . 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.       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

SPF a Cup with a Cylindrical Boss
The goal here is to form an AZ31B part with a primary and a secondary feature. Figure 12 shows a three-dimensional representation of this relatively complicated part. It consists of a cup (a primary feature) with a cylindrical boss (a secondary feature) at its base. It is needless to say that this component is difficult to form using conventional forming techniques. To produce this part, a closed die with an inner cavity was used. The surface of the inner cavity has the same geometry as the desired shape of the finished part. Figure 13 shows a finished part after 40 min of SPF using a constant gas pressure of 0.15 MPa.  The constitutive model developed in the previous section was incorporated in a finite element analysis carried out using ABAQUS. An axisymmetric model was used for the simulation of the SPF of the part that is shown in Figure 14. The 102-mm diameter circular die with a die entry radius ER of 3 mm was modeled as a rigid analytical surface. The 1.2-mm thick sheet was modeled using four layers of axisymmetric elements. Each layer contained 300 solid elements. The sheet elements between the die and the blank holder were assumed to be fixed. This was necessary in order to simulate the pure stretching condition in the SPF process. Since no lubricant was used in the actual forming, a coefficient of coulomb friction of 0.8 at the die-sheet interface was considered in the simulations [24]. The height of the boss (h) is 6 mm. The diameters of the boss are 22.8 mm and 40.6 mm for (d1) and (d2), respectively. (R1) and (R2) is 3.6 mm and 6 mm, respectively. A comparison between the part shape predicted by the simulation and the experimentally measured shape is shown in Figures 15-17 for forming using a pressure of 0.15 MPa considering the cup height (H) values of 25, 30, and 35 mm, respectively. The results demonstrate excellent agreement between the experimental and the simulation results for the three considered cases. However, the results of one of the forming experiments, indicated by the blue triangles in Figure 16, differ from those predicted by simulation at the die bottom radius. This occurs at the last stage of the forming process and can be attributed to variations in the material properties of the sheet. When the experiments were repeated, as indicated by the red circle in the same figure, matching results were obtained.      In Figures 18-20, arrows were drawn pointing at the regions corresponding to the cup entry radius, cup bottom radius and the boss entry radius. From the experimental results, abrupt changes in the thickness distribution are obvious around the cup entry radius and the boss entry radius. This effect was clearly predicted by the simulation results. The behavior of thinning in the vicinity of the die entry radius has been reported by many researchers [24,[31][32][33]. It is attributed to the combined compression and bending effects caused by the small die entry radius, which usually takes place at the beginning of the forming process. However, due to the high value of the coefficient of friction, since no lubricant was used, a near-sticking condition will develop once the sheet gets in contact with the die. Therefore, less material will flow towards the remaining regions of the sheet and more stretching will develop there. Consequently, localized thinning will occur at the die bottom corner, since that region of the sheet is the last to get in contact with the die surface before the forming process ends. This phenomenon was also demonstrated in both the experimental and simulation results in Figures 18-20.

SPF of a Cup with a Car-shaped Boss
Another complicated geometry was considered for the verification of the constitutive material model. It consisted of a cup with a car-shaped boss. The actual formed part is shown in Figure 21. The forming time was 46.7 min, and the temperature was kept at 450 °C. The part profile and the thickness along two directions (direction 1 and direction 2) were measured for each experiment. A three-dimensional model in ABAQUS was used to simulate the SPF of the cup with a car-shaped boss. The circular sheet was meshed using 2500 linear quadrilateral shell elements. The die, on the other hand, was modeled using three-dimensional rigid elements.
A comparison between the part shape predicted by the simulation and the experimentally measured shape is shown in Figures 22 and 23 for forming using a pressure of 0.3 MPa considering a cup height (H) value of 20 mm. The results demonstrate an excellent agreement between the experimental and the simulation results.  The predicted thickness distribution of the formed car-shape part is shown in Figures 24-25. It is interesting to see that the thickness profile and the geometry of the cross-section along each direction show some resemblance. This is expected, since, at the higher parts of the geometry, the sheet will get in contact with the die surface early during the forming process before any considerable stretching occurs. Due to the high coefficient of friction, once in contact no more material flow will happen there, and therefore, the thickness will be relatively higher. Exactly the opposite can be said about the lower parts of the geometry, as they are the last to be in contact with the sheet, a low value of sheet thickness will take place at those regions.
Most of the region along the formed part surface displays a good correspondence between the predicted results and the experimental ones. However, there is a noticeable gap at the bottom radius for direction 1. The prediction of the thickness distribution for direction 2 is slightly better, with only a small difference between the simulation and experimental results at the top of the car-shape profile.

Conclusions
The nonlinearities incorporated in the superplastic forming process alongside with the complexity in the material behavior, makes constitutive modeling of such materials a hard challenge. In this study, a constitutive model for AZ31B was established based on a comprehensive numerical-experimental approach. The model takes into account grain growth and cavity formation. The model was calibrated by fitting to data from free bulge forming and microstructure analysis. The fitting process was based on analytical and optimization techniques. Taguchi method was applied in the optimization process. This procedure can be implemented in constitutive modeling of other superplastic materials.
In order to verify the established model, two relatively complicated geometries were simulated, and the results were compared with those obtained from actual SPF experiments. Results show that the material model is capable of predicting the thickness distribution and the shape of the formed part accurately. In fact, the close accord between the experimentally obtained results and the simulation predictions demonstrate the significance of including the microstructure evolution in the AZ31B material model. It is needless to say: This model can be easily implemented in the simulation of the production of complex parts, usually through a material subroutine.