AZ31 Sheet Forming by Clustering Ball Spinning-Analysis of Damage Evolution Using a Modiﬁed GTN Model

: Clustering ball spinning (CBS) forming is a novel approach to manufacturing a complex curved surface. In order to explore the forming limit of magnesium alloy in the CBS forming process, the modiﬁed GTN model was incorporated into the FE simulation to analyze the damage evolution. The theoretical analyses are conducted to investigate the deformation mechanism and to explore the stress state in the CBS forming process. The numerical results show that the modiﬁed GTN model can predict the result more accurately compared with the standard GTN model, and the damage parameters for GTN model are determined by the experiments. Besides, the forming limit of AZ31 magnesium alloy can be improved by CBS forming method. To explore the reason for the increased forming limit, the microstructure of curved surface was tested by electron backscattered scattering detection (EBSD). The results demonstrate that the deformation of magnesium alloy plate by the CBS method is dominated initially by the extension twinning, and non-basal slip systems are activated with the development of forming process.


Introduction
Clustering ball spinning (CBS) forming is a novel forming process used to produce complex curved surfaces. Compared with the traditional press forming process using male and female die, CBS forming exhibits its inherent advantages and flexibility such as simple tooling, low forming force, and improved mechanical properties. It can be seen from Figure 1 that the CBS method manufactures curved metal workpieces with discrete rigid balls instead of female die. Hu et al. [1,2] had demonstrated that CBS method has great forming performance for aluminum alloys, and proposed a multi-layer rigid ball (MLRB) model that can simulate the process of CBS accurately. Although CBS forming had high precision at room temperature, the theoretical study and the forming limit of magnesium alloy in CBS forming remained to be solved. There are highly non-linear tooling/sheet interactions and complex local stress states in the CBS forming process. For low-ductility materials such as magnesium alloy, cracking often took place during the CBS forming process due to plastic deformation surpassing the general limit of flow formability.
The FE simulation combined with ductile fracture criteria is an efficient way to analyze the ductile fracture of metal [3,4]. The ductile fracture criteria are divided into coupled and uncoupled models by defining whether the damage variable is included in the material constitutive model. The simulation results of uncoupled fracture criteria such as Rice-Tracey, McClintock and Cockcroft-Lathan are highly dependent on the size of mesh [5]. The other category is coupled ductile fracture criteria which couples the influence of plastic damage into the material constitutive model such as the Lemaitre model and Gurson-Tvergaard-Needleman (GTN) damage model [6]. The typical model is the GTN damage model which attributes the ductile failure process to micro voids nucleation, growth and coalescence [7]. Yan et al. [8] proved the GTN damage model is valid for silicon steel during the rolling process. Riedel et al. [9] supposed Gurson model had limitations arising from the fact that only spherical voids were considered and the model was generally too stable. In fact, the original GTN model is not appropriate for predicting damage evolution well under low stress triaxiality [10] and negative stress triaxiality. To overcome the limitation of stress triaxiality, the standard GTN constitutive model had been improved [11][12][13][14][15][16]. Nahshon et al. [17] modified the GTN model by introducing a phenomenological extension to the porosity evolution equation. SUN et al. [18] applied modified GTN model to analyze the ductile damage and failure behavior of steel sheet with edge defects under multi-pass cold rolling. In this paper, the GTN model modified by Nahshon and Hutchinson was incorporated into FE software (Abaqus) by user-defined material subroutine (VUMAT). In order to optimize the CBS forming process, it is necessary to analyze and predict the forming limit of sheet during the forming process. In the paper, a complex curved surface with the AZ31 sheet was fabricated by the CBS forming process and the forming limit of magnesium was predicted by FE model cooperated with modified GTN model.

Materials and Methods
In this experiment, the commercial magnesium alloy AZ31 with the size of 200 mm × 200 mm × 1.5 mm was employed, and the compositions were listed in Table 1. The forming device of CBS is shown in Figure 1, and fundamental parameters are listed in Table 2. As shown in Figure 2, the forming device included the following sections: rotatable gland stirrer, rigid balls, cylinder cavity, blank holder, and the bottom die. The forming process of AZ31 sheet was shown in Figure 2, and the detailed experimental steps were shown as follows: Prepared magnesium alloy sheet was placed between the blank holder and the cylinder cavity (applied with 2.5 MPa pressure), while allowing the rotatable gland stirrer to contact with rigid balls and the bottom die to contact with the magnesium alloy sheet. It can be seen from Figure 2 that the forming process is divided into two steps: Step 1 was that the bottom die rose while the rotatable gland stirrer was not moved and rotated.
Step 2 was the bottom die rose while the rotatable gland stirrer was moved down and rotated. It can be known that rigid balls have a complicated influence on the sheet forming, which needs further analysis by theoretical formula and FE simulation.

Parameters Value
The radius of the rotatable gland stirer (mm) 75 The radius of rigid balls (mm) 4 The Rotational speed of rotatable gland stirere (r/min) 24 The radius of the bottom die 46

Contact Stress
In order to explore the stress state at the boundary of the die, membrane analysis of the contact stresses on the neutral surface was conducted to develop a theoretical model. Silva et al. [19] proposed a theoretical model for stress analysis based on membrane approach to address the material deformation in ISF process. Chang et al. [20] updated the membrane theory by taking into account the hardening of the material. In the present work, membrane analysis was applied to the CBS process. Figure 3a showed the two contact regions between rigid balls and the sheet in the forming regions and AF regions. Figure 3b showed the stress components on a small element in region A. As shown in Figure 3c,d, along the directions of thickness t, circumference θ and meridian ϕ, the corresponding stress components σ t , σ tb , σ θ , σ ϕ were considered in the membrane analysis. σ t represents the opposite stress of the rigid ball on the sheet, while σ tb defined as the stress of the bottom die acting on the sheet. t represents the thickness of the sheet, and r b represents the radius of the rigid ball. The assumptions of the theoretical analysis were summarized as, (1) Bending effect is neglected as the sheet thickness is much smaller than the tool radius.
(2) Shear stress τ tθ and τ tϕ are taken into account due to the shearing effect on the formability, while the other shear stress τ θ ϕ is neglected because it has no obvious effect on the plastic deformation.
(3) σ t is related to t and not related to θ: the through-thickness stress is linearly distributed along the thickness direction and evenly distributed along the circumferential direction.
Along the circumferential direction, the force equilibrium equation can be expressed as, In addition can be simplified by neglecting higher order terms and assuming rdθ ≈ t., Resolving the force equilibrium along the thickness direction, In addition can be simplified by, Considering the Levy-Mises constitutive equation in plane strain, the equivalent stress can be calculated by, By combining the above equations, the following equations can be derived as, The hydrostatic stress of the contact zone can be calculated as, The stress triaxiality of the sheet η t in region I during the CBS process can be obtained as: The stress component σ tb can be neglected in region II and the stress triaxiality of the sheet η AF can be obtained as, From Equations (9) and (10), it can be seen that the plate is in a lower stress triaxiality and the stress triaxiality decreased with the radius of rigid balls decreasing. In order to solve the problem that the GTN damage model cannot accurately describe material damage at low stress triaxiality, Nahshon and Hutchinson introduced a phenomenological extension to the porosity evolution equation.

Constitutive Model
Gurson proposed a continuous damage model by accounting the hydrostatic stress and the model was modified by Tvergaard and Needleman in order to encapsulate the nucleation, growth and coalescence of voids. The damage law can be expressed as follows: where the macroscopic Von-Mises stress is defined as, Cauchy stress tensor is represented by σ, σ h is the hydraulic stress, σ is the deviatoric stress, I is the identity tensor and σ s = σ s ε m q is a function of the equivalent plastic strain.
Nahshon and Huchinson modified the damage law to accommodate the damage contribution from shear since it could play a major role in the forming process.
The evolution of porosity and shear damage is given in Equations (12) and (13).
. f n is the void volume rate due to the nucleation of new voids. Once nucleated, porosity growth . f g is determined by the rate of the plastic strain, while . f s is determined by the shear deformation. The shear damage evolution is determined by damage parameter k w .
The effective volume fraction f * in the Gurson potential function is described by the following equation.
where f f denotes the shear damage at complete shear failure.
where J 3 = det(s) represents the third invariant of stress state, A is a function of the equivalent plastic strain, w(σ) is a function determined by the parameter L = − 27J 3 2σ 3 e . Yield stress material is considered as a dependent variable of plastic strain, strain rate and temperature as given in Equation (17).

Elastic-Plastic Behavior
A series of tests at different strain rates were carried out to determine the model parameters. The dominant forming temperature in the CBS forming process is 423 k. Figure 4 illustrated that the stress-strain relationship at different strain rates at 423 K. It can be seen from Figure 4 that the stress is in a relatively stable state when ε = 0.1, and the stress at this strain is called σ 0.1 .The equation can be obtained by taking the natural logarithm to Equation (17). ε . Figure 5b illustrated the relationship between ln σ 0.1 and temperature T. The stress-strain relationship at 0.01 s at 298 K was shown in Figure 5c. It can be concluded that b = 0.0026, m = 0.05544. The relationship between ln σ y and ln ε was plotted in Figure 5d, and it can be concluded that n = 0.1546. Parameter α can be obtained as follow: where k = 5.6647 was the intercept obtained from the curve of ln σ y and ln ε, and α = 1118.5629.

Fractional Factor of GTN Model
The GTN constitutive parameters q 1 = 1.5, q 2 = 1 and q 3 = q 2 1 were determined by Tverggard and were applied in much research. The damage parameters of AZ31 alloy sheet according to the previous research were listed in Table 3. For most structural alloys, k w should lie in the range 1 < k w < 3. The value of k w is related to the stress triaxiality and the Lode parameter. The formula for the Lode parameter can be expressed as follow: Combining Equation above, the Lode parameter in the CBS forming process can be obtained as follow: The stress component σ tb can be neglected in region II and the Lode parameter of the sheet η AF can be obtained as, It can be found from Equation (22) that the stress triaxiality is related to the sheet thickness and the radius of rigid balls. It can be calculated that in region II the stress triaxiality is about 0.26 and the Lode parameter is about 0.68. In addition, the parameter k w is set to 3 according to the value of the stress triaxiality and the Lode parameter.

The FE Model of Shear-Tensile
In order to verify the accuracy of the constitutive model, the shear-tensile test was established and simulated by ABAQUS/Explicit. The GTN damage model was integrated with ABAQUS/Explicit through user-defined material parameters. It can be seen from Figure 6 that three models with different numbers of meshes were built separately, in order to compare the influence of the number of meshes on the results of numerical simulations. Theoretically, the constitutive model in the paper requires a comparatively low number of meshes, and it is important to choose the suitable amount of mesh. Figure 7 showed the Load-Displacement curves results and crack shapes of models with different numbers of mesh. When the model was meshed with 6852 elements and 3 layers, it can be seen that the simulation results are completely different from the other two and do not match the actual situation in terms of crack shape. The most suitable number of meshes is 14,720 with 4 layers, and the increasing number of meshes had a little influence on the calculation results when compared with the model of 37,420 elements both in crack shape and the Load-Displacement curves. The simulation of the CBS forming process is more complex and reducing the number of meshes can greatly improve the efficiency of the simulation without influencing the accuracy of the simulation. In order to explore the effect of k w on the simulation results, the values of parameter k w are set to 0, 1, 2, 3, respectively. It can be seen that k w = 2 is a more accurate damage parameter in the shear-tension model. It can be seen from Figure 8 that when k w increased, the elongation of the material decreased.

Step 1
The result of modified GTN model during CBS forming process was discussed here in order to optimize the process. It was proposed in Section 2 that the forming process was divided into two steps.
Step 1 was that the bottom die raised up while the rotatable gland stirrer was not moved and rotated, and the temperature of magnesium alloy sheet was 293 K.
As the sheet cracked in the CBS forming mainly depended on void volume fraction and damage parameter k w in the modified GTN model, the damage parameter k w was set as 0, 1, 2, 3, respectively, to compare the effect of k w on material damage. The distribution of effective void volume faction f * at different k w was shown in Figure 9. The maximum void volume fraction f * predicted by GTN model increased with the value of k w . In the modified GTN model with the damage parameter k w = 3, the sheet cracked as the bottom die reached 13.9 mm. In the GTN model with k w = 0, the sheet cracked when the displacement of the bottom die reached 15 mm. Besides, it can be noted that the maximum void volume fraction was insufficient at a low deformation and increased dramatically at a high deformation. The crack occurred when the bottom die displacement reached 14 mm in the CBS forming experiment. Compared the results of simulation with different values of k w in Step 1, the modified model with k w = 3 was similar with the experiment well.
The purpose of Step 1 was to make the sheet contact with the bottom die and increase the size of the forming region. The effective forming radius used to define the size of the forming region and the variation of maximum void volume fraction of the sheet in Step 1 was shown in Figure 10. It should be noticed the maximum void volume fraction increased rapidly when the displacement of the bottom die increased to 10 mm and the growth rate of effective forming radius decreased. It can be concluded the suitable displacement of the bottom die was 5 mm and the growth rate of the effective forming radius reached the maximum without increasing void volume fraction.

Step 2
When the displacement of the bottom die reached 5 mm, it was time to move and rotate the rotatable gland stirrer in order to increase the size of the effective forming radius. The temperature of the magnesium sheet raised to 423 K due to the spinning of rigid balls, and the formability of magnesium sheet had been improved. The logarithmic strain distribution and void volume fraction under different damage parameter k w are shown in Figure 11. It can be found that the maximum displacement of the bottom die decreased from 30.8 mm to 27.3 mm with the increasing value of k w . It can be concluded that the forming limit of the AZ31 sheet was improved by Step 2, and the damage parameter k w = 3 are similar with the result of experiment. According to the damage evolution Equations (15) and (17) during plastic forming, the equivalent plastic strain was an important factor that led to material failure. It can be seen from Figure 12b the equivalent plastic strain of all regions increased with the displacement of the bottom die at different regions and the maximum strain took place at Region B, which is the boundary of forming region and AF region. Besides, the equivalent plastic strain of Region C increased due to the pressure of rigid balls. In fact, the effective forming radius increased by the extra pressure. When the bottom die reached the displacement of 27.3 mm, the void of fraction reached the maximum at Region B and crack occurred. Besides, the effective radius reached 25 mm which is similar with the result of the experiment (24.7 mm). It can be concluded that CBS forming process can improve the forming limit of the sheet and the effective forming radius. Figure 12d shows the evolution of stress triaxiality in different regions and the stress triaxiality decreased when the rigid ball pressed on the surface of the sheet. It can be found that the stress triaxiality influenced by the rigid balls and it is necessary to apply modified GTN model to the FE simulation. Figure 13 shows the photographs of cracked AZ31 magnesium alloy in CBS forming process. It can be seen from Figure 13 that the displacement of the bottom die increased from 14 mm to 27 mm and the maximum effective increased from 12.5 mm to 24.7 mm. It can be concluded that the rotation of rigid balls can improve the forming limit of AZ31 magnesium alloy and increase the effective forming radius. In order to explore the mechanism of deformation, the microstructure of curved surface was tested by EBSD.  The EBSD inverse pole figure map (IPF) of as-received AZ31 sheet about the distribution of the grain orientations was shown in Figure 14a. Figure 14b showed the (0002) pole figure plots of the as-received AZ31 sheet. It can be seen from Figure 14a,b that the c-axes in the grains were slightly tilted away from the normal of the sheet and lie on the ND-RD plane. The EBSD Inverse pole figure map (IPF) of the deformed AZ31 sheet was shown in Figure 14c, and the (0002) pole figure plots of the deformed AZ31 sheet was shown in Figure 14d. It can be seen from Figure 14c,d that the c-axes in the grains were rotated to 86 • and lay nearly parallel to the RD. It was suggested that the extension twinning participate fully in the deformation process. Meanwhile, it can be seen from Figure 14b that the c-axes in the grains had a small angle rotation which was caused by the shear deformation during CBS forming process. During CBS forming process, a large amount of friction heat improved the temperature of the sheet up to 423 K, which would further provide the energy to active the pyramidal slip system <c+a> and induce the shear deformation of the sheet. In all, the deformation of magnesium alloy plate using the CBS method is dominated initially by the extension twinning while non-basal slip system and shear deformation take place with the spinning of rigid balls.

Conclusions
In the paper, the CBS forming process was optimized by FE simulation and simulations. A modified GTN model extended to shear deformation was implemented in ABAQUS/Explicit to analyze the damage evolution. The main conclusions are summarized as follows: 1.
The material parameters of the GTN model for AZ31 magnesium alloy were determined by experiment, and the reliability of modified GTN model was verified through the comparison of numerical simulation and experimental observation on the sheartensile process with the AZ31 magnesium alloy. Theoretical formula is conducted to contact stress analysis in the CBS forming process, and it can be found that the stress triaxiality is related to the radius of rigid balls and the thickness of the sheet. Besides, the stress triaxiality was about 0.26 and the Lode parameter was about 0.68 in the AF region, and the forming region was under the negative stress triaxiality in this simulations and experiments. 2.
The numerical results are compared under different damage parameter k w , it can be found that the results of damage parameterk w = 3 are similar with that of experiments both in Step 1 and Step 2. The maximum displacement of the rotatable gland stirrer reached 27.3 mm and the effective forming radius reached 24.7 mm in Step 2, which was much higher than that in Step 1 (The maximum displacement was 13.9 mm, and the effective forming radius was 12 mm in Step 1). With spinning of rigid balls, forming limit was improved and the effective forming radius increased. 3.
To explore the increased forming limit in the CBS forming process, the microstructure of curved surface was tested by electron backscattered scattering detection (EBSD). It is indicated from the microstructure that the deformation of AZ31 sheet using the CBS method is dominated initially by the extension twinning and non-basal slip system is activated with development of forming and the improvement of the temperature. Meanwhile, it can be found that the c-axes in the grains had a small angle rotation which was caused by the shear deformation during CBS process.