Aluminum Alloy Sheet-Forming Limit Curve Prediction Based on Original Measured Stress–Strain Data and Its Application in Stretch-Forming Process

: A new method, by directly utilizing original measured data (OMD) of the stress–strain relation in the Marciniak–Kuczynski (M–K) model, was proposed to predict the forming limit curve (FLC) of an aluminum alloy sheet. In the groove zone of the M–K model, by establishing the relations of the equivalent strain increment, the ratio of shear stress to the ﬁrst principle stress and the ratio of the second principle stress to the ﬁrst principle stress, the iterative formula was established and solved. The equations of theoretical forming limits were derived in detail by using the OMD of the stress–strain relation. The stretching specimens of aluminum alloy 6016-T4 were tested and the true stress–strain curve of the material was obtained. Based on the numerical simulations of punch-stretch tests, the optimized specimens’ shape and test scheme were determined, and the tests for FLC were carried out. The FLC predicted by the proposed method was more consistent with the experimental results of FLC by comparing the theoretical FLCs based on OMD of the stress–strain relation and of that based on traditional power function. In addition, the inﬂuences of anisotropic parameter and groove angle on FLCs were analyzed. Finally, the FLC calculated by the proposed method was applied to analyze sheet formability in the stretch-forming process, and the predicted results of FLC were veriﬁed by numerical simulations and experiments. The fracture tendency of the formed parts can be visualized in the forming limit diagram (FLD), which has certain guiding signiﬁcance for fracture judgment in the sheet-forming process.


Introduction
Plastic forming of sheet metal keeps an important role in modern industrial production and is widely applied to the aerospace industry, automobile manufacturing, equipment manufacturing, and other fields [1][2][3][4]. Mechanical properties and plastic-forming abilities of alloy sheet were studied [5][6][7], and structural and textural evolutions in plastic forming were analyzed by the view of microstructure [8,9]. The forming limit is an important performance index and process parameter in sheet metal forming [10,11]. In the sheet formability evaluation method, the forming limit curve (FLC) is the most intuitive and widely used [12]. Therefore, predicting FLC accurately and effectively is worthy of further study and test [13].
To predict the FLC precisely, various related theoretical models have been developed. These theoretical models include Swift's theory of dispersion instability [14], Hill's concentrated instability theory [15], and Marciniak-Kuczynski's (M-K) theory based on thickness inhomogeneity assumptions, as well as the Storen-Rice theory [16] and the modified maximum force criterion (MMFC) theory [17,18]. achieve the complete fit between the sheet metal and the die, so the sheet metal was easily cracked [41]. In order to improve the forming quality, the integral clamps in a conventional stretch-forming process were replaced by discrete clamps in a stretch-forming process based on loading at multi-positions. Cai et al. proposed a new design method for a stretch-forming process based on a minimum deformation path, and which can achieve the manufacture of specified shape parts and the faultless product of large curvature [42]. In order to explore the optimal moving locus of discrete clamps, a series of numerical simulations and experiments were carried out by Yang et al. [43], and the results demonstrated that the strains and stresses of deformed parts presented different distributions when different loading trajectories were adopted. In this paper, the FLC was used to analyze the stretch-forming process of sheet metal, therefore, the related tensile test for stress-strain curve and the punch-stretch experiment for FLC were carried out under the temperature and strain rate similar to those of the stretch-forming process. The new method for FLC prediction directly uses original measured data (OMD) of the stress-strain relation to replace power approximations obtained by fitting the stress-strain data. In addition, the OMD of the stress-strain relation under room temperature was more realistic in the analysis on sheet formability in the stretch-forming process, because the stretch-forming process was generally performed under room temperature and the OMD of the stress-strain relation can strictly indicate the true stress-strain relation of the material. The forming limits by the proposed method were calculated by the detailed theoretical equations. By comparing the theoretical FLC calculated by different predicted methods with experiment results, the accuracy and superiority of the proposed method to predict the FLC of aluminum alloy sheet was verified. Finally, the predicted FLC was applied to analyze the fracture of the curved part in the stretch-forming process, which provided the theoretical basis for judging the fracture defects generated in the stretch-forming process.

Uniaxial Tensile Tests
To obtain the true stress-strain curves of 6016-T4 aluminum alloy sheet, uniaxial tensile tests on the specimens in 0 • , 45 • , and 90 • to the rolling direction were performed in a DNS100 microcomputer-controlled and universal material testing machine, and the direction of 0 • is the rolling direction. The specimen had the strain gauge length of 55 mm, the section width of 15 mm, and thickness of 1 mm. The true stress-strain curves of different orientations were plotted by analyzing experimental data in Figure 1. It was found that the true stress-strain curves had the same trend in three directions, and the difference was not particularly significant. So, the average value of the stress at a certain strain in each direction was the final measured data of stress-strain. The mechanical properties of 6016-T4 aluminum alloy sheet are shown in Table 1, and the mechanical properties of the material in different directions were not much different. The anisotropic parameter (r-values) can be calculated by measuring the gauge length and width of the specimen before and after deformation. The r-values in three orientations were similar and less than 1, which Metals 2019, 9, 1129 4 of 20 indicated the strain of the sheet metal was greater in the thickness direction than that in width direction. In engineering, the average value is generally used as the anisotropic parameter of sheet metal [20].

M-K Model Based on OMD of Stress-Strain
Marciniak and Kuczynski [13] proposed that necking is usually caused by geometric or structural imperfection and the original randomly distributed defects of sheet metal are assumed to be the weak area, i.e., the groove zone, and the thickness inhomogeneity is used to indicate the unevenness of the sheet metal. As shown in Figure 2, the sheet metal was divided into the uniform zone a and the weak zone b. The 1 and 2 axes indicate the principal stress directions in zone a, and the 3 axis indicates the thickness direction of the sheet metal. The normal and the tangential directions of the groove region are denoted by n and t, respectively. Moreover, the angle between the groove direction and the 2 direction is θ. For simplicity, the superscripts 'a' and 'b' are used to distinguish the corresponding variables in zones a and b.
Metals 2019, 9, x FOR PEER REVIEW 4 of 20 width direction. In engineering, the average value is generally used as the anisotropic parameter of sheet metal [20].

M-K Model Based on OMD of Stress-Strain
Marciniak and Kuczynski [13] proposed that necking is usually caused by geometric or structural imperfection and the original randomly distributed defects of sheet metal are assumed to be the weak area, i.e., the groove zone, and the thickness inhomogeneity is used to indicate the unevenness of the sheet metal. As shown in Figure 2, the sheet metal was divided into the uniform zone a and the weak zone b. The 1 and 2 axes indicate the principal stress directions in zone a, and the 3 axis indicates the thickness direction of the sheet metal. The normal and the tangential directions of the groove region are denoted by n and t, respectively. Moreover, the angle between the groove direction and the 2 direction is  . For simplicity, the superscripts 'a' and 'b' are used to distinguish the corresponding variables in zones a and b. In the M-K model, the force in the n direction is transmitted across zones a and b, therefore, the force equilibrium conditions can be obtained: where nn  and nt  are, respectively, normal stress and shear stress in the groove local coordinate system (n-t coordinate system), and f is the thickness inhomogeneity.
The equation can be obtained from Equation (1): In addition, the relation of stress components between the n-t coordinate system and principal axis coordinate system (1-2 coordinate system) can be expressed as follows:  In the M-K model, the force in the n direction is transmitted across zones a and b, therefore, the force equilibrium conditions can be obtained: where σ nn and σ nt are, respectively, normal stress and shear stress in the groove local coordinate system (n-t coordinate system), and f is the thickness inhomogeneity. The equation can be obtained from Equation (1): In addition, the relation of stress components between the n-t coordinate system and principal axis coordinate system (1-2 coordinate system) can be expressed as follows: σ nn = σ 11 cos 2 θ + σ 22 sin 2 θ + σ 12 sin 2θ σ nt = (σ 22 − σ 11 ) sin θ cos θ + σ 12 cos 2θ (3) There are not shear stress components in zone a, i.e., σ a 12 = 0. Based on Equation (3), the stress components in zones a and b satisfy: . By substituting Equation (4a) and (4b) into Equation (2), the equation can be found: where δ b is the function of α b . In the M-K model, the strain paralleling to the groove in zone b is constrained by the strain in uniform zone a, so the compatibility condition can be given by Moreover, the relation of strain increments between n-t and the 1-2 coordinates can be expressed as follows: dε tt = dε 11 sin 2 θ + dε 22 cos 2 θ − dε 12 sin 2θ (7) There are not shear strain increments in zone a, i.e., dε a 12 = 0. Based on Equation (7), the strain increments in zones a and b can satisfy: dε a tt = dε a 1 sin 2 θ + dε a 2 cos 2 θ dε b tt = dε b 11 sin 2 θ + dε b 22 cos 2 θ − dε b 12 sin 2θ (8) Combining with the plastic flow equation dε ij = dε(∂σ/∂σ ij ) () and the equivalent stress of Hill 1948 (σ = 3(1 + r)/[2(2 + r)] σ 2 11 − 2rσ 11 σ 22 /(1 + r) + σ 2 22 + 2(1 + 2r)σ 2 12 /(1 + r)), the major and minor strain increments in zones a and b can be expressed: ). By substituting Equations (8), (9a), and (9b) into Equation (6), the equation can be found: where, The measured stress-strain curve can be reflected by a series of OMD of stress-strain (ε i , σ i ), which can be recorded as When the OMD of stress-strain were recorded as n + 1 nodes and ε 0 < ε 1 < · · · < ε i < ε i+1 < · · · < ε n , as shown in Figure 3. For any equivalent strain ε ε i , firstly, determine the strain interval in which the strain is located, i.e., ε i < ε < ε i+1 , then, the equivalent stress at any ε can be calculated: where the equivalent stress at any equivalent strain can be calculated by OMD of stress-strain.
Metals 2019, 9, x FOR PEER REVIEW 6 of 20 The measured stress-strain curve can be reflected by a series of OMD of stress-strain ( , ) ii  , which can be recorded as When the OMD of stress-strain were recorded as n + 1 nodes and Figure 3. For any equivalent strain i   , firstly, determine the strain interval in which the strain is located, i.e., , then, the equivalent stress at any  can be calculated: where the equivalent stress at any equivalent strain can be calculated by OMD of stress-strain. can be transformed to where () a F  and () b F  can be calculated by Equation (12). Equation (13), the iterative solution equation can be found: The detailed steps of the Newton-Raphson method for solving the iterative solution equation (14) are described below: According to the M-K theory, zone a is loaded proportionally, and the plastic deformation conforms to the total plastic theory. Firstly, assume the equivalent strain increment a d  (  is a very small constant, such as 0.001

 
) in zone a. When solving the iterative equation for the kth time, the equivalent stress based on OMD of stress-strain relation can be expressed as: (14) where, 0  is the prestrain of the sheet metal.
In this case, the iteration equation is: In addition, based on the Equations (4a), (4b), and (12), σ a nn = f σ b nn can be transformed to where F(ε a ) and F(ε b ) can be calculated by Equation (12). By substituting f = t b /t a , f 0 = t b 0 /t a 0 , t = t 0 exp(ε 33 ), and Equations (5), (10), and (12) into Equation (13), the iterative solution equation can be found: The detailed steps of the Newton-Raphson method for solving the iterative solution equation (14) are described below: According to the M-K theory, zone a is loaded proportionally, and the plastic deformation conforms to the total plastic theory. Firstly, assume the equivalent strain increment dε a = ς (ς is a very small constant, such as ς = 0.001) in zone a. When solving the iterative equation for the kth time, the equivalent stress based on OMD of stress-strain relation can be expressed as: where, ε 0 is the prestrain of the sheet metal. In this case, the iteration equation is: where, .
Solving the iteration equation ψ k (α k b ) = 0, Newton's iteration formula is: in Equation (17), Iteratively solve Equation (17) to α k+1 b − α k b ≤ eps (eps is very small, such as eps = 10 −10 ), when the iterative equation converges, the unknown quantity α k b in zone b can be obtained. (5) and (10), and check dε k b /dε k a . If dε k b /dε k a ≥ λ (λ is an integer, such as λ = 10), the formulas for calculating the major and minor strains when the sheet metal is necked are: Otherwise, continue to solve equations ψ k+1 (α k+1 b ) = 0 by using Newton's iteration method until the condition of dε k+1 b /dε k+1 a ≥ λ can be satisfied. Change the groove angle θ, and the major and minor strains at different groove angles can be calculated separately. Finally, the minimum values of major and minor strains are taken as the ultimate strains.
Then, change α a (α a ∈ [0, 1]), and the theoretical prediction of the major and minor strains under different loading paths can be obtained, and the FLC can be drawn.
The detailed steps of the Newton-Raphson method for calculating the ultimate strain are described below in the corresponding flow chart, as shown in Figure 4.

Experiments' Determination on FLC
The FLC should include the loading path from uniaxial tension to plane strain to balanced biaxial tension, so that the complete FLC can be obtained. In previous studies, the numerical-experimental approach with the specimens of different shapes and widths was employed to calculate forming limits through uniaxial tensile, plane strain state, and biaxial tensile [44][45][46]. In this section, the punch-stretch tests for specimens of different shapes were simulated, then, the optimized specimens' shape and punch-stretch test scheme were determined. The punch-stretch tests were carried out, and the experimental results were compared with the theoretical FLCs based on different stress-strain data.

Experiments' Determination on FLC
The FLC should include the loading path from uniaxial tension to plane strain to balanced biaxial tension, so that the complete FLC can be obtained. In previous studies, the numerical-experimental approach with the specimens of different shapes and widths was employed to calculate forming limits through uniaxial tensile, plane strain state, and biaxial tensile [44][45][46]. In this section, the punch-stretch tests for specimens of different shapes were simulated, then, the optimized specimens' shape and punch-stretch test scheme were determined. The punch-stretch tests were carried out, and the experimental results were compared with the theoretical FLCs based on different stress-strain data.

Numerical Simulations and Experiment Scheme
The numerical simulations were performed by Dynaform. Specimens of different widths and shapes were used to obtain the forming limits under different strain paths. In the simulation process, it was considered that the concave die, the blank holder, and the spherical punch were defined as the rigid body and meshed with R3D4 elements; the sheet metal can be deformed and meshed with S4R elements. In addition, the mesh size of all parts was 2 mm, and the symmetric constraints were loaded at the symmetric boundary of parts. The contact state between the contact pairs is decided by a master-slave algorithm, and the surfaces of the rigid body were defined as the master surface, while the surface of the deformed sheet metal was defined as the slave surface. The blanking force of 20 kN was applied to the blank holder. After the process parameters were set, the numerical simulation task was submitted for analysis, and the simulation results were analyzed in the post-processing program.
The punch-stretch test was widely used to obtain forming limits. The effective forming area of the specimen was used to design a rectangular strip, however, the numerical simulation results showed that the specimen with the small width was easily broken at the blank holder, as shown in Figure 5a. Since the rectangular specimen was subjected to the large concentrated stress at the blank holder, the large major strain appeared there, which may eventually result in the reduction of forming limits. In order to prevent failure from occurring at the blank holder region, the specimen was designed as dumbbell-shaped in the effective forming area, that is, the intermediate area was designed as the arc, which can smoothly transit from the center to both sides. The maximum value of the major strain of the dumbbell-shaped specimen appeared in center region, as shown in Figure 5. The numerical simulations were performed by Dynaform. Specimens of different widths and shapes were used to obtain the forming limits under different strain paths. In the simulation process, it was considered that the concave die, the blank holder, and the spherical punch were defined as the rigid body and meshed with R3D4 elements; the sheet metal can be deformed and meshed with S4R elements. In addition, the mesh size of all parts was 2 mm, and the symmetric constraints were loaded at the symmetric boundary of parts. The contact state between the contact pairs is decided by a master-slave algorithm, and the surfaces of the rigid body were defined as the master surface, while the surface of the deformed sheet metal was defined as the slave surface. The blanking force of 20 kN was applied to the blank holder. After the process parameters were set, the numerical simulation task was submitted for analysis, and the simulation results were analyzed in the post-processing program.  The punch-stretch test was widely used to obtain forming limits. The effective forming area of the specimen was used to design a rectangular strip, however, the numerical simulation results showed that the specimen with the small width was easily broken at the blank holder, as shown in Figure 5(a). Since the rectangular specimen was subjected to the large concentrated stress at the blank holder, the large major strain appeared there, which may eventually result in the reduction of forming limits. In order to prevent failure from occurring at the blank holder region, the specimen was designed as dumbbell-shaped in the effective forming area, that is, the intermediate area was designed as the arc, which can smoothly transit from the center to both sides. The maximum value of the major strain of the dumbbell-shaped specimen appeared in center region, as shown in Figure 5.
Numerical simulation results of specimens with different sizes showed that, as the width of the specimen increased, the influence of the specimen shape on the simulation result was gradually weakened. As shown in Figure 5(b), when the width of the effective formed area was about 60 mm, the simulation results of the two kinds of specimens were not much different. Considering the convenience and feasibility of processing and testing, the specimens with a large width were designed as rectangle. Therefore, the shape and size of the specimens were finally determined, and the experiment scheme is as shown in Figure 6 and Table 2, respectively, where h = 90 mm. Numerical simulation results of specimens with different sizes showed that, as the width of the specimen increased, the influence of the specimen shape on the simulation result was gradually weakened. As shown in Figure 5b, when the width of the effective formed area was about 60 mm, the simulation results of the two kinds of specimens were not much different. Considering the convenience and feasibility of processing and testing, the specimens with a large width were designed as rectangle. Therefore, the shape and size of the specimens were finally determined, and the experiment scheme is as shown in Figure 6 and Table 2, respectively, where h = 90 mm.

Analysis of Experiments Results
The punch-stretch tests were carried out. As shown in Figure 7a, the equipment adopted a two-stage hydraulic stretching and crimping system, which applied uniform blanking force to the specimens and improved the stability of the punch-stretch process. The ultimate strains of the sheet were obtained on the ARGUS optical strain gauge produced by GOM. This equipment can obtain the high-precision 3D full field strain data of the deformed sheet, and the main compositions of the test equipment are shown in Figure 7b.

Analysis of Experiments Results
The punch-stretch tests were carried out. As shown in Figure 7(a), the equipment adopted a two-stage hydraulic stretching and crimping system, which applied uniform blanking force to the specimens and improved the stability of the punch-stretch process. The ultimate strains of the sheet were obtained on the ARGUS optical strain gauge produced by GOM. This equipment can obtain the high-precision 3D full field strain data of the deformed sheet, and the main compositions of the test equipment are shown in Figure 7  In this experiment, the meshing method used was the circular grid division. The circular grid is the most widely used, and it is convenient to determine the major and minor strains. At the end of the deformation of the aluminum alloy sheet, the shape of the circular grid became the elliptical grid, and the major and minor strains of the sheet metal were calculated by the following formula: 1 1 0 where d0 is the original diameter of the circular grid, and d1 and d2 are the long and short axes of the ellipse grid after deformation, respectively. In this experiment, the meshing method used was the circular grid division. The circular grid is the most widely used, and it is convenient to determine the major and minor strains. At the end of the deformation of the aluminum alloy sheet, the shape of the circular grid became the elliptical grid, and the major and minor strains of the sheet metal were calculated by the following formula: where d 0 is the original diameter of the circular grid, and d 1 and d 2 are the long and short axes of the ellipse grid after deformation, respectively. According to the designed test scheme, punch-stretch tests of different specimens were carried out. The deformed specimens are shown in Figure 8. According to the fracture position of the specimen and the deformation degree of the grid, the deformed regions of the sheet metal were divided into three categories: Fracture region, critical region, and safety region. It can be seen that the cracking degree of the specimen and the deformation degree of the mesh were different, so a mass of data needs to be extracted to predict the FLC of the sheet metal reasonably. Metals 2019, 9,   According to the designed test scheme, punch-stretch tests of different specimens were carried out. The deformed specimens are shown in Figure 8. According to the fracture position of the specimen and the deformation degree of the grid, the deformed regions of the sheet metal were divided into three categories: Fracture region, critical region, and safety region. It can be seen that the cracking degree of the specimen and the deformation degree of the mesh were different, so a mass of data needs to be extracted to predict the FLC of the sheet metal reasonably. The ARGUS optical strain gauge can obtain the deformed mesh distribution on the surface of the deformed specimen. Then, the major and minor strains' distributions of the deformed specimens can be further analyzed by the software. The major and minor strain of the specimens with the width of 30 mm and 60 mm are shown in Figure 9. The cracking degree of specimens with different sizes is different. Therefore, in order to predict the FLC of the 6016-T4 aluminum alloy, the strain data in the The ARGUS optical strain gauge can obtain the deformed mesh distribution on the surface of the deformed specimen. Then, the major and minor strains' distributions of the deformed specimens can be further analyzed by the software. The major and minor strain of the specimens with the width of 30 mm and 60 mm are shown in Figure 9. The cracking degree of specimens with different sizes is different. Therefore, in order to predict the FLC of the 6016-T4 aluminum alloy, the strain data in the safety and fracture zones of specimens were extracted. It can be found that the major strains of the fracture zone were significantly larger than those of the safety zone, and the smaller the width of the specimen, the more likely the minor strain was negative value. The extraction locations of major and minor strains in the safety and fracture zones are also marked in Figure 9.  According to the designed test scheme, punch-stretch tests of different specimens were carried out. The deformed specimens are shown in Figure 8. According to the fracture position of the specimen and the deformation degree of the grid, the deformed regions of the sheet metal were divided into three categories: Fracture region, critical region, and safety region. It can be seen that the cracking degree of the specimen and the deformation degree of the mesh were different, so a mass of data needs to be extracted to predict the FLC of the sheet metal reasonably. The ARGUS optical strain gauge can obtain the deformed mesh distribution on the surface of the deformed specimen. Then, the major and minor strains' distributions of the deformed specimens can be further analyzed by the software. The major and minor strain of the specimens with the width of 30 mm and 60 mm are shown in Figure 9. The cracking degree of specimens with different sizes is different. Therefore, in order to predict the FLC of the 6016-T4 aluminum alloy, the strain data in the For the punch-stretch test, the major and minor strains in the safety and fracture zones are expressed in Figure 10. The strain data in the fracture and safety zones were extracted, and the strain data distributions were scattered, so it was difficult to present a clear FLC. The majority of strain data distributed in the fracture and safety zones, and very few strain data were located in the intermediate critical region. Instead, the FLC band was presented with a certain width. The two curves can be plotted to limit the dispersed strain points to a certain zone, that is, the critical zone is between the upper limit curve of the safety zone and the lower limit curve of the fracture zone. The FLC band can be regarded as the critical zone of fracture. safety and fracture zones of specimens were extracted. It can be found that the major strains of the fracture zone were significantly larger than those of the safety zone, and the smaller the width of the specimen, the more likely the minor strain was negative value. The extraction locations of major and minor strains in the safety and fracture zones are also marked in Figure 9. For the punch-stretch test, the major and minor strains in the safety and fracture zones are expressed in Figure 10. The strain data in the fracture and safety zones were extracted, and the strain data distributions were scattered, so it was difficult to present a clear FLC. The majority of strain data distributed in the fracture and safety zones, and very few strain data were located in the intermediate critical region. Instead, the FLC band was presented with a certain width. The two curves can be plotted to limit the dispersed strain points to a certain zone, that is, the critical zone is between the upper limit curve of the safety zone and the lower limit curve of the fracture zone. The FLC band can be regarded as the critical zone of fracture.

Comparison of FLCs Calculated by OMD of Stress-Strain and Power Function
Uniaxial tensile tests were carried out in the directions of 0°, 45°, and 90° with the rolling direction, and three sets of true stress-strain curves with little differences were obtained. In order to improve the accuracy of the experiment, the average values were taken as OMD of stress-strain, as shown in Figure 11. The equivalent stress at any equivalent strain based on OMD of stress-strain can be calculated by Equation (12). Then, the FLC based on OMD of stress-strain can be calculated by the theoretical derivation in the third section.
The power function stress-strain curve (  Figure 11. It can be seen that the power function of the stress-strain curve was restricted by the equation form, which deviated from the OMD of stress-strain.

Comparison of FLCs Calculated by OMD of Stress-Strain and Power Function
Uniaxial tensile tests were carried out in the directions of 0 • , 45 • , and 90 • with the rolling direction, and three sets of true stress-strain curves with little differences were obtained. In order to improve the accuracy of the experiment, the average values were taken as OMD of stress-strain, as shown in Figure 11. The equivalent stress at any equivalent strain based on OMD of stress-strain can be calculated by Equation (12). Then, the FLC based on OMD of stress-strain can be calculated by the theoretical derivation in the third section. The FLCs based on different stress-strain relations were calculated, and shown in Figure 12. It was found that the theoretical FLCs based on power function of stress-strain curves were higher than the experimental FLC under the negative-strain path. The left side of FLCs based on power functions of stress-strain curves was above the critical zone of punch-stretch tests, but the right side of FLC was in the critical zone measured by punch-stretch tests. Moreover, the theoretical FLC calculated based on the power function ( . Finally, the theoretical FLC calculated by OMD of stress-strain was between Figure 11. OMD of stress-strain and power function fitting curves. The power function stress-strain curve (σ = Kε n ) was obtained by fitting experimental stress-strain data. The parameters K and n in power function based on fitting the complete stress-strain curve (including elastic and plastic deformation) were 501.99 MPa and 0.227, respectively. In addition, since the fracture occurred in the plastic deformation stage, the power function based on the stress-strain relation of plastic deformation stage was σ = 412.45ε 0.142 , as shown in Figure 11. It can be seen that the power function of the stress-strain curve was restricted by the equation form, which deviated from the OMD of stress-strain.
The FLCs based on different stress-strain relations were calculated, and shown in Figure 12. It was found that the theoretical FLCs based on power function of stress-strain curves were higher than the experimental FLC under the negative-strain path. The left side of FLCs based on power functions of stress-strain curves was above the critical zone of punch-stretch tests, but the right side of FLC was in the critical zone measured by punch-stretch tests. Moreover, the theoretical FLC calculated based on the power function (σ = 412.45ε 0.142 ) was closer to the test value than that based on σ = 501.99ε 0.227 . Finally, the theoretical FLC calculated by OMD of stress-strain was between the safety and fracture zones of punch-stretch tests, that is, the predicted major and minor strains were in the critical zone of the necking. By comparison, it was found that the stress-strain curves of sheet metal had significant impact on FLC under the negative-strain path, and the theoretical FLC based on OMD of stress-strain was more consistent with the experimental results, which verified the accuracy and superiority of the proposed method in this paper. The FLCs based on different stress-strain relations were calculated, and shown in Figure 12. It was found that the theoretical FLCs based on power function of stress-strain curves were higher than the experimental FLC under the negative-strain path. The left side of FLCs based on power functions of stress-strain curves was above the critical zone of punch-stretch tests, but the right side of FLC was in the critical zone measured by punch-stretch tests. Moreover, the theoretical FLC calculated based on the power function ( . Finally, the theoretical FLC calculated by OMD of stress-strain was between the safety and fracture zones of punch-stretch tests, that is, the predicted major and minor strains were in the critical zone of the necking. By comparison, it was found that the stress-strain curves of sheet metal had significant impact on FLC under the negative-strain path, and the theoretical FLC based on OMD of stress-strain was more consistent with the experimental results, which verified the accuracy and superiority of the proposed method in this paper.

Influence of Anisotropic Parameter and Groove Angle on FLCs
It is well known that the von Mises yield function considers the anisotropic parameter (r) equals 1, and the Hill 1948 yield function can include a non-zero anisotropic parameter. In Figure 13, PF1 and PF2 represent power functions by fitting the complete (

 
), respectively. It was found that the anisotropic parameter had significant impact on FLC, and the ultimate strains were greater when r = 1 than when r < 1 in three prediction methods, especially, in the negative-strain path. It can be found from Table 1 that the

Influence of Anisotropic Parameter and Groove Angle on FLCs
It is well known that the von Mises yield function considers the anisotropic parameter (r) equals 1, and the Hill 1948 yield function can include a non-zero anisotropic parameter. In Figure 13, PF1 and PF2 represent power functions by fitting the complete (σ = 501.99ε 0.227 ) and plastic deformation stress-strain curve (σ = 412.45ε 0.142 ), respectively. It was found that the anisotropic parameter had significant impact on FLC, and the ultimate strains were greater when r = 1 than when r < 1 in three prediction methods, especially, in the negative-strain path. It can be found from Table 1 that the aluminum alloy 6016-T4 was anisotropy. It was more reasonable to predict the FLCs by using r < 1. Moreover, the FLCs calculated by the M-K model combining OMD of stress-strain were more consistent with the punch-stretch test results than those based on power functions. aluminum alloy 6016-T4 was anisotropy. It was more reasonable to predict the FLCs by using 1 r  .
Moreover, the FLCs calculated by the M-K model combining OMD of stress-strain were more consistent with the punch-stretch test results than those based on power functions. From the above analysis, it can be seen that, when the FLCs of the sheet metal were predicted by the M-K model combining with OMD of stress-strain and power functions, the FLCs under the negative-strain path were observably influenced by stress-strain relation. The FLCs under the negative-strain paths were minimized under non-zero groove angles. However, the FLCs under the positive-strain path were minimized when the groove angle was zero. In order to analyze the effect of the groove angle on the FLCs, the left side of the FLC was calculated with zero and non-zero groove angles, as shown in Figure 14(a). It can be seen that the FLCs predicted by three methods were higher in zero angle than in non-zero angle, and the predicted FLCs were closer to the experimental data in non-zero groove angle. Therefore, in order to accurately predict the FLCs, the FLCs in the negative-strain paths should be calculated at non-zero groove angles. In addition, when the ultimate strains take the minimum values, the variations of groove angle with different strain paths are shown in Figure 14(b). It was found that, when the ultimate strain took the minimum value, the groove angle increased firstly and then decreased when sheet deformations changed from the uniaxial tension to the plane strain state. The groove angle at the plane strain state was zero with the minimum strain, which was consistent with the minimum value of the forming limit in the right side when the groove angle was zero. Therefore, the groove angle had a certain influence on the left side of the FLC.  From the above analysis, it can be seen that, when the FLCs of the sheet metal were predicted by the M-K model combining with OMD of stress-strain and power functions, the FLCs under the negative-strain path were observably influenced by stress-strain relation. The FLCs under the negative-strain paths were minimized under non-zero groove angles. However, the FLCs under the positive-strain path were minimized when the groove angle was zero. In order to analyze the effect of the groove angle on the FLCs, the left side of the FLC was calculated with zero and non-zero groove angles, as shown in Figure 14a. It can be seen that the FLCs predicted by three methods were higher in zero angle than in non-zero angle, and the predicted FLCs were closer to the experimental data in non-zero groove angle. Therefore, in order to accurately predict the FLCs, the FLCs in the negative-strain paths should be calculated at non-zero groove angles. In addition, when the ultimate strains take the minimum values, the variations of groove angle with different strain paths are shown in Figure 14b. It was found that, when the ultimate strain took the minimum value, the groove angle increased firstly and then decreased when sheet deformations changed from the uniaxial tension to the plane strain state. The groove angle at the plane strain state was zero with the minimum strain, which was consistent with the minimum value of the forming limit in the right side when the groove angle was zero. Therefore, the groove angle had a certain influence on the left side of the FLC. aluminum alloy 6016-T4 was anisotropy. It was more reasonable to predict the FLCs by using 1 r  .
Moreover, the FLCs calculated by the M-K model combining OMD of stress-strain were more consistent with the punch-stretch test results than those based on power functions. From the above analysis, it can be seen that, when the FLCs of the sheet metal were predicted by the M-K model combining with OMD of stress-strain and power functions, the FLCs under the negative-strain path were observably influenced by stress-strain relation. The FLCs under the negative-strain paths were minimized under non-zero groove angles. However, the FLCs under the positive-strain path were minimized when the groove angle was zero. In order to analyze the effect of the groove angle on the FLCs, the left side of the FLC was calculated with zero and non-zero groove angles, as shown in Figure 14(a). It can be seen that the FLCs predicted by three methods were higher in zero angle than in non-zero angle, and the predicted FLCs were closer to the experimental data in non-zero groove angle. Therefore, in order to accurately predict the FLCs, the FLCs in the negative-strain paths should be calculated at non-zero groove angles. In addition, when the ultimate strains take the minimum values, the variations of groove angle with different strain paths are shown in Figure 14(b). It was found that, when the ultimate strain took the minimum value, the groove angle increased firstly and then decreased when sheet deformations changed from the uniaxial tension to the plane strain state. The groove angle at the plane strain state was zero with the minimum strain, which was consistent with the minimum value of the forming limit in the right side when the groove angle was zero. Therefore, the groove angle had a certain influence on the left side of the FLC.

Application of FLC in the Stretch-Forming Process Based on Loading at Multi-Positions
The FLC was applied to evaluate whether the three-dimensional curved parts appear in the fracture defect during the stretch-forming process, and the fracture defect can be predicted by FLC.
The stretch-forming process is an important method for forming three-dimensional curved parts. It is widely used to form the skin parts of aircraft and spacecraft, the front of high-speed trains, the outer plates of ships, and the large-curvature curved parts of modern buildings. The integral clamps in a conventional stretch-forming process were replaced by discrete clamps based on loading at multi-positions, as shown in Figure 15. The displacement load was applied at a series of discrete clamps on the two ends of the sheet metal, and the load of each loading point was controlled in real time, so that the deformed sheet could actively deform along with the die profile. The discrete clamps made the sheet metal bend and stretch according to the die surface, and ultimately realized the complete fit between the sheet metal and the die, forming high-precision surface parts.

Application of FLC in the Stretch-Forming Process Based on Loading at Multi-Positions
The FLC was applied to evaluate whether the three-dimensional curved parts appear in the fracture defect during the stretch-forming process, and the fracture defect can be predicted by FLC. The stretch-forming process is an important method for forming three-dimensional curved parts. It is widely used to form the skin parts of aircraft and spacecraft, the front of high-speed trains, the outer plates of ships, and the large-curvature curved parts of modern buildings. The integral clamps in a conventional stretch-forming process were replaced by discrete clamps based on loading at multi-positions, as shown in Figure 15. The displacement load was applied at a series of discrete clamps on the two ends of the sheet metal, and the load of each loading point was controlled in real time, so that the deformed sheet could actively deform along with the die profile. The discrete clamps made the sheet metal bend and stretch according to the die surface, and ultimately realized the complete fit between the sheet metal and the die, forming high-precision surface parts. The loading trajectory had a great influence on forming results of a curved-surface part, and stress and strain distributions of curved parts were different under different loading trajectories. In this section, numerical simulation and experimental verification were carried out for spherical parts under different loading paths. Combining with the prediction of sheet metal-forming limits, the fracture defects of curved-surface parts were predicted by FLC obtained by the proposed method.

Loading Trajectories for the Stretch-Forming Process Based on Loading at Multi-Positions
It is well known that the loading trajectories of the discrete clamps have great influence on the forming precision and quality of a three-dimensional curved part. From a geometric point of view, there are countless loading trajectories from flat sheet to three-dimensional curved parts. In order to compare and analyze the influence of different loading paths on the forming results of curved parts, the stretch-forming process based on equal and variable elongation was designed to form curved parts.  Since the influence of the sheet width on the longitudinal strain was relatively small, it was The loading trajectory had a great influence on forming results of a curved-surface part, and stress and strain distributions of curved parts were different under different loading trajectories. In this section, numerical simulation and experimental verification were carried out for spherical parts under different loading paths. Combining with the prediction of sheet metal-forming limits, the fracture defects of curved-surface parts were predicted by FLC obtained by the proposed method.

Loading Trajectories for the Stretch-Forming Process Based on Loading at Multi-Positions
It is well known that the loading trajectories of the discrete clamps have great influence on the forming precision and quality of a three-dimensional curved part. From a geometric point of view, there are countless loading trajectories from flat sheet to three-dimensional curved parts. In order to compare and analyze the influence of different loading paths on the forming results of curved parts, the stretch-forming process based on equal and variable elongation was designed to form curved parts.
Since the influence of the sheet width on the longitudinal strain was relatively small, it was considered that the sheet metal was composed of several longitudinal strips with the same width as the discrete clamps, and the corresponding strips were deformed by controlling the displacements at the loading points. The deformation process of ith longitudinal strip is shown in Figure 16. The loading trajectory was the moving locus of the end of the material line during the process. It was seen that the material line changed from its initial length to an intermediate length and finally to the length corresponding to the arc length of the stretching die. In the stretch-forming process based on equal elongation (SF-e), it was assumed that the elongation of the longitudinal strip of the sheet metal was equal, and the displacement of each loading point was calculated separately. In the stretch-forming process based on variable elongation (SF-v), the elongations of the longitudinal strips of the sheet metal varied and were determined by cubic functions in the forming process. The details for deriving the process of equal elongation loading mode and variable elongation loading mode can be found in reference [43].
It is well known that the loading trajectories of the discrete clamps have great influence on the forming precision and quality of a three-dimensional curved part. From a geometric point of view, there are countless loading trajectories from flat sheet to three-dimensional curved parts. In order to compare and analyze the influence of different loading paths on the forming results of curved parts, the stretch-forming process based on equal and variable elongation was designed to form curved parts. Since the influence of the sheet width on the longitudinal strain was relatively small, it was When the sheet length was 400 mm, the sheet width was 300 mm and the spherical part with the radius of 500 mm was formed in numerical simulations. The loading trajectories of different loading points in different loading methods (SF-e and SF-v) are shown in Figure 17. considered that the sheet metal was composed of several longitudinal strips with the same width as the discrete clamps, and the corresponding strips were deformed by controlling the displacements at the loading points. The deformation process of ith longitudinal strip is shown in Figure 16. The loading trajectory was the moving locus of the end of the material line during the process. It was seen that the material line changed from its initial length to an intermediate length and finally to the length corresponding to the arc length of the stretching die. In the stretch-forming process based on equal elongation (SF-e), it was assumed that the elongation of the longitudinal strip of the sheet metal was equal, and the displacement of each loading point was calculated separately. In the stretch-forming process based on variable elongation (SF-v), the elongations of the longitudinal strips of the sheet metal varied and were determined by cubic functions in the forming process. The details for deriving the process of equal elongation loading mode and variable elongation loading mode can be found in reference [43]. When the sheet length was 400 mm, the sheet width was 300 mm and the spherical part with the radius of 500 mm was formed in numerical simulations. The loading trajectories of different loading points in different loading methods (SF-e and SF-v) are shown in Figure 17.

Numerical Simulation Results and Experimental Validation
The spherical parts with the radius of 500 mm were formed by two loading methods, and the major strain distributions of the numerical simulations are shown in Figure 18. It can be seen from Figure 18 that when different loading trajectories were adopted to form spherical parts of the same curvature, the major strain distributions in the effective formed region had the same tendency, but the spherical part in SF-v had more uniform strain distribution than in SF-e. The maximum major strain of the spherical part in SF-e was 0.4993, which was much larger than 0.1502 in SF-v. Moreover, the spherical part formed by SF-e had the larger strains in the middle region that was adjacent to the discrete clamps, and the larger strains may have caused the formation of fracture defect in the region. The appearance of the fracture defect reduced the accuracy of the curved part and even produced unacceptable curved parts.

Numerical Simulation Results and Experimental Validation
The spherical parts with the radius of 500 mm were formed by two loading methods, and the major strain distributions of the numerical simulations are shown in Figure 18. It can be seen from Figure 18 that when different loading trajectories were adopted to form spherical parts of the same curvature, the major strain distributions in the effective formed region had the same tendency, but the spherical part in SF-v had more uniform strain distribution than in SF-e. The maximum major strain of the spherical part in SF-e was 0.4993, which was much larger than 0.1502 in SF-v. Moreover, the spherical part formed by SF-e had the larger strains in the middle region that was adjacent to the discrete clamps, and the larger strains may have caused the formation of fracture defect in the region. The appearance of the fracture defect reduced the accuracy of the curved part and even produced unacceptable curved parts. In order to study the fracture defect under different loading modes, the major and minor strains of the spherical parts in a critical section of fracture were extracted. It can be seen from Figure 19 that, in SF-e, the strains in a critical section of fracture exceeded the theoretically calculated ultimate strains and were located in the fracture region measured by experiment. The regions of strain points above the FLC appeared in fracture defect, and the spherical part in SF-e was prone to fracture in plane strain state. The fracture defects reduced the accuracy of the formed parts, and the defective parts were manufactured in SF-e. However, in SF-v, all strain points in the critical section of fracture were below the FLC and were located in the safe region measured by experiment, which indicated that the possibility of fracture defect in SF-v was especially small. Meanwhile, it was found that the FLC predicted by combining the M-K theory and the OMD of stress-strain can quantitatively analyze the fracture defect, which has certain guiding significance for the sheet-forming process. In order to verify the numerical simulation results of the stretch-forming process in different loading trajectories, the stretch-forming experiments for spherical parts of aluminum alloy 6016-T4 were carried out on a homemade apparatus by Jilin University, as shown in Figure 15. The critical sections of fracture of experimental parts are shown in Figure 20. It can be seen that the large deformation was produced in SF-e, and the severe fracture defect appeared in the center region that was adjacent to the discrete clamps. However, in SF-v, as indicated by experimental result, the target surface part was perfectly manufactured without fracture defect. In conclusion, the experiment results verified the accuracy of the numerical simulation and FLC prediction, and the fracture defect appeared in the region where the strain points exceeded the theoretically calculated ultimate strains. In order to study the fracture defect under different loading modes, the major and minor strains of the spherical parts in a critical section of fracture were extracted. It can be seen from Figure 19 that, in SF-e, the strains in a critical section of fracture exceeded the theoretically calculated ultimate strains and were located in the fracture region measured by experiment. The regions of strain points above the FLC appeared in fracture defect, and the spherical part in SF-e was prone to fracture in plane strain state. The fracture defects reduced the accuracy of the formed parts, and the defective parts were manufactured in SF-e. However, in SF-v, all strain points in the critical section of fracture were below the FLC and were located in the safe region measured by experiment, which indicated that the possibility of fracture defect in SF-v was especially small. Meanwhile, it was found that the FLC predicted by combining the M-K theory and the OMD of stress-strain can quantitatively analyze the fracture defect, which has certain guiding significance for the sheet-forming process.  In order to study the fracture defect under different loading modes, the major and minor strains of the spherical parts in a critical section of fracture were extracted. It can be seen from Figure 19 that, in SF-e, the strains in a critical section of fracture exceeded the theoretically calculated ultimate strains and were located in the fracture region measured by experiment. The regions of strain points above the FLC appeared in fracture defect, and the spherical part in SF-e was prone to fracture in plane strain state. The fracture defects reduced the accuracy of the formed parts, and the defective parts were manufactured in SF-e. However, in SF-v, all strain points in the critical section of fracture were below the FLC and were located in the safe region measured by experiment, which indicated that the possibility of fracture defect in SF-v was especially small. Meanwhile, it was found that the FLC predicted by combining the M-K theory and the OMD of stress-strain can quantitatively analyze the fracture defect, which has certain guiding significance for the sheet-forming process. In order to verify the numerical simulation results of the stretch-forming process in different loading trajectories, the stretch-forming experiments for spherical parts of aluminum alloy 6016-T4 were carried out on a homemade apparatus by Jilin University, as shown in Figure 15. The critical sections of fracture of experimental parts are shown in Figure 20. It can be seen that the large deformation was produced in SF-e, and the severe fracture defect appeared in the center region that was adjacent to the discrete clamps. However, in SF-v, as indicated by experimental result, the target surface part was perfectly manufactured without fracture defect. In conclusion, the experiment results verified the accuracy of the numerical simulation and FLC prediction, and the fracture defect appeared in the region where the strain points exceeded the theoretically calculated ultimate strains. In order to verify the numerical simulation results of the stretch-forming process in different loading trajectories, the stretch-forming experiments for spherical parts of aluminum alloy 6016-T4 were carried out on a homemade apparatus by Jilin University, as shown in Figure 15. The critical sections of fracture of experimental parts are shown in Figure 20. It can be seen that the large deformation was produced in SF-e, and the severe fracture defect appeared in the center region that was adjacent to the discrete clamps. However, in SF-v, as indicated by experimental result, the target surface part was perfectly manufactured without fracture defect. In conclusion, the experiment results verified the accuracy of the numerical simulation and FLC prediction, and the fracture defect appeared in the region where the strain points exceeded the theoretically calculated ultimate strains.

Conclusions
In this paper, the new method by directly utilizing OMD of stress-strain relation in the M-K model was proposed to predict the FLC of an aluminum alloy sheet, and the predicted FLC in this method was compared with that of the M-K model based on power functions and punch-stretch tests. In addition, the predicted FLC in the proposed method was used to analyze the sheet formability in the stretch-forming process. The conclusions of this paper can be summarized as follows: (1) In this paper, the new method by directly utilizing OMD of stress-strain relation in the M-K model was proposed to predict FLC of an aluminum alloy sheet, and the detailed theoretical derivations for calculating forming limits were presented.
(2) The stress-strain curve had significant effect on the FLC prediction of sheet metal. The power functions were fitted by the OMD of stress-strain, but the true stress-strain relation of the material was not strictly reflected by power function, which may have affected the prediction results of the sheet-forming limits. However, the OMD of stress-strain can strictly reflect the true stress-strain relation of the material, so the FLC prediction results based on OMD of stress-strain can be more accurate.
(3) By comparing the FLCs based on OMD of stress-strain and power functions with experiment result, it was found that the theoretical FLCs based on power functions were higher on the left side than the experimental FLC, and the right side of FLC was in the critical zone measured by experiment. However, the theoretical FLC based on OMD of stress-strain were in the critical zone measured by experiment. The predicted FLC in the proposed method were more consistent with the experimental results, and the accuracy and superiority of the proposed method were verified in this paper.
(4) Anisotropic parameter affected the prediction accuracy of FLC, and the groove angle had significant effect on the FLC under the negative-strain path. The predicted forming limits with zero groove angle were higher than that when the groove angle was not zero. In this paper, the forming limits calculated in the conditions of r<1 and 0   were more consistent with the experimental FLC.
(5) The FLC calculated by the proposed method was applied to analyze sheet formability in the stretch-forming process, and the predicted results of FLC were verified by numerical simulations and experiment of the stretch-forming process. The fracture tendency of the weakened region of the formed parts could be intuitively visualized and quantitatively analyzed in FLD, which had certain guiding significance for fracture judgment in the sheet-forming process.

Conclusions
In this paper, the new method by directly utilizing OMD of stress-strain relation in the M-K model was proposed to predict the FLC of an aluminum alloy sheet, and the predicted FLC in this method was compared with that of the M-K model based on power functions and punch-stretch tests. In addition, the predicted FLC in the proposed method was used to analyze the sheet formability in the stretch-forming process. The conclusions of this paper can be summarized as follows: (1) In this paper, the new method by directly utilizing OMD of stress-strain relation in the M-K model was proposed to predict FLC of an aluminum alloy sheet, and the detailed theoretical derivations for calculating forming limits were presented.
(2) The stress-strain curve had significant effect on the FLC prediction of sheet metal. The power functions were fitted by the OMD of stress-strain, but the true stress-strain relation of the material was not strictly reflected by power function, which may have affected the prediction results of the sheet-forming limits. However, the OMD of stress-strain can strictly reflect the true stress-strain relation of the material, so the FLC prediction results based on OMD of stress-strain can be more accurate.
(3) By comparing the FLCs based on OMD of stress-strain and power functions with experiment result, it was found that the theoretical FLCs based on power functions were higher on the left side than the experimental FLC, and the right side of FLC was in the critical zone measured by experiment. However, the theoretical FLC based on OMD of stress-strain were in the critical zone measured by experiment. The predicted FLC in the proposed method were more consistent with the experimental results, and the accuracy and superiority of the proposed method were verified in this paper.
(4) Anisotropic parameter affected the prediction accuracy of FLC, and the groove angle had significant effect on the FLC under the negative-strain path. The predicted forming limits with zero groove angle were higher than that when the groove angle was not zero. In this paper, the forming limits calculated in the conditions of r < 1 and were more consistent with the experimental FLC.
(5) The FLC calculated by the proposed method was applied to analyze sheet formability in the stretch-forming process, and the predicted results of FLC were verified by numerical simulations and experiment of the stretch-forming process. The fracture tendency of the weakened region of the formed parts could be intuitively visualized and quantitatively analyzed in FLD, which had certain guiding significance for fracture judgment in the sheet-forming process.