A Modiﬁed M-K Method for Accurate Prediction of FLC of Aluminum Alloy

: Forming limit curve (FLC) is an important failure criterion for sheet metals in sheet metal forming, while the M-K model is widely used for the prediction of the FLC. In the M-K model, such prediction is greatly inﬂuenced by the initial thickness imperfection factor and material properties, from which the original M-K model’s theoretical derivation is proposed as a solution to the above mentioned issue in this paper. Then the relationship between the M-K model and Keeler’s empirical formula is then studied, from which a new method to predict FLC is proposed that combines the M-K model with Keeler’s empirical formula according to the previous analyses. It turns out that this new method can simplify the calculation procedure. Finally, the experimental results of two kinds of aluminum alloys, AA6016 and AA5182, have veriﬁed the effectiveness of the proposed method.


Introduction
Lightweight materials are currently a hot issue in the automotive and aerospace industry. As light alloy material, aluminum alloy is of small destiny, has a high ratio of strength to weight and has high corrosion resistance ability [1][2][3]. Therefore, it is increasing in popularity. However, aluminum alloys are quite different from automotive steel sheets in ductility, anisotropy, yield strength and other indicators, so the study of its formability is of great significance. The forming limit curve (FLC) is the most intuitive and effective formability indicator for evaluating the forming performance of metal sheets. It is a curve determined by the instability limit strains of the sheet metal under different strain paths [4]. During the forming process, fracture occurs when the minor and major strain points of the sheet metal are above the FLC. As FLC is widely used to accurately predict the forming limit, great attention has been paid to it [5,6].
To quickly and easily obtain the FLC of metal sheets, the following theoretical models have been proposed, based on the plastic forming theory. The Swift model is based on a unique relationship between the root-mean-square values of shear stress and incremental strain [7]. The Hill model is derived from permissible discontinuities of stress, velocity, and surface slope [8]. The modified maximum force criterion (MMFC) predicts localized necking under non-proportional loading [9]. The M-K model takes the initial thickness imperfection into account [10]. It is widely used because of its simplicity and high precision, and its original model has already been improved significantly by scholars. For further researches on the M-K method, the angle between the groove and direction of the principal stress has been introduced in prediction [11][12][13][14]. Although scholars have done a lot of research on it, there are still some disadvantages. The predictions made by the M-K model are affected by many factors, such as the yield criterion of the materials and the termination conditions of the algorithm. The performance of the original M-K model on prediction of the FLC for aluminum alloy is in need of further research. So further study is needed to improve the M-K method.
Metals 2021, 11, 394 2 of 11 In this study, the M-K model's theoretical derivation is proposed. Besides, the effect of material properties and the termination condition of the algorithm are carried out, and then the relation between the M-K model and Keeler's empirical formula [15] is studied. Based on the previous analyses about such relation, a new method to predict FLC is proposed, which combines the M-K model with Keeler's empirical formula, and can greatly simplify the calculation. Finally, the experiments of two aluminum alloys, AA6016 and AA5182, are introduced to verify the effectiveness of this method.

Original M-K Model
The schematic diagram of groove model is shown in Figure 1. The plastic flow follows the normality rule, which states where f is yield function; σ 1 , σ 2 , and σ 3 are stress components that correspond to the principal axis of stress.
Metals 2021, 11, x FOR PEER REVIEW 2 of 12 prediction of the FLC for aluminum alloy is in need of further research. So further study is needed to improve the M-K method. In this study, the M-K model's theoretical derivation is proposed. Besides, the effect of material properties and the termination condition of the algorithm are carried out, and then the relation between the M-K model and Keeler's empirical formula [15] is studied. Based on the previous analyses about such relation, a new method to predict FLC is proposed, which combines the M-K model with Keeler's empirical formula, and can greatly simplify the calculation. Finally, the experiments of two aluminum alloys, AA6016 and AA5182, are introduced to verify the effectiveness of this method.

Theoretical Analysis
The schematic diagram of groove model is shown in Figure 1. The plastic flow follows the normality rule, which states where f is yield function; 1 σ , 2 σ , and 3 σ are stress components that correspond to the principal axis of stress. According to the hypothesis of incompressibility, we have In order to take sheet metal's normal anisotropy into account, the Hill'48 plastic yield [16] criterion is applied. The following formulas can be obtained. According to the hypothesis of incompressibility, we have In order to take sheet metal's normal anisotropy into account, the Hill'48 plastic yield [16] criterion is applied. The following formulas can be obtained.
where r is normal anisotropic coefficient. α and β are defined as follows: Considering stresses are parallel to principal axis of stress, the relation between dε i and dε 1 can be derived by the equivalent plastic work which is expressed as From Equation (1), β is derived as a function of α, while β = β(α). It is supposed that the strains of the zone A and the zone B are equal in the 2-direction: force balance condition in 1-direction is expressed as: And the error is set as: Based on the above analysis, the stresses and strains of zone A can be obtained when the equivalent strain increment dε A i and the ratio α A of zone A are given. The equivalent strain increment dε B i of zone B can be determined by Newton iterative algorithm with the objective function (9). The strain of zone B is considered as the fracture strain when the main strain of the zone B is far greater than that of zone A. Setting then the termination of the algorithm can be written as ϕ > ϕ T . The recommended value of ϕ T is stated in this paper. The flow chart of each path is shown in Figure 2.

Relationship between M-K Model and Keeler's Empirical Formula
The hardening model is expressed as

Relationship between M-K Model and Keeler's Empirical Formula
The hardening model is expressed as During each iteration, σ A 1 t A is a constant, while the stress component σ B 1 and the thickness t B are determined by the strain increment of zone B dε B i . That is to say, e is a function of dε B i . From Equations (2) and (6), the strain increment dε 3 is a function of dε B i . Accordingly, we have Divide both sides of Equation (10) by K, we have Therefore, the calculated FLCs are independent of strength coefficient K, which is similar to Keeler's empirical formula [15]. However, there is a slight difference between Keeler's empirical formula and M-k method. The former is dependent on the thickness of the sheet metal, which is expressed as Equations (15) and (16), while the latter introduces the thickness imperfection.

The Modified M-K Model
The contrastive study points out the predicted forming limit curves, obtained by the original M-K model, do not match the experiment results for either AA6016 or AA5182. However, the value of FLC 0 accords with the experiments. In order to predict the FLC more precisely, the M-K method is combined with Keeler's empirical formula, where FLC 0 is calculated by the M-K method, and the accurate FLC is calculated according to Equation (15).
In the following equations, the calculation of FLC 0 will be simplified. When the minor strain is 0, from Equations (4) and (5), we have According to Equations (10), (12) and (18), Equation (14) can be written as where t B 0 and t A 0 are the initial thickness of zone B and A, let f 0 = t B 0 /t A 0 denote the initial thickness imperfection factor. In each iteration, e is dependent only on dε B i . Equations (14) and (19) have similarities. Meanwhile, Equation (19) expresses that the FLC 0 only depends on the hardening exponent and the initial thickness imperfection factor.

Experiments to Determine the Material Properties
Aluminum alloy 6016-T4 and 5182-O (sheet thickness is 1 mm) are applied in the experiments of this paper. The geometrical dimension of the dog bones is shown in Figure 3, designed according to the standard GB/T 228.1-2010 [17]. The stretch testing machine is SANS (50 KN, MTS Industrial Systems (China) CO., LTD., Shenzhen, China), and the strains are measured by the extensometer with 25 mm-gauge-length. During the test, the chuck moves at a constant speed of 3 mm/s. The stress-strain response of AA6016 and AA5182 is expressed in Figure 4. The stressstrain response are fitted by power function and the adjusted R-squared is greater than 0.995. The fitting results are expressed in Table 1 and Figure 4.  Figure 4. The stress-strain response of AA6016 and AA5182.
The forming limit tests are carried out in accordance with GB/T15825.631995 [18]. The drawing of the test tooling is shown in Figure 5. The strain path (ratio of major to minor strain) is changed according to the change of specimen width and the lubrication condition. The test pieces are rectangular sheets shown in Figure 6. The width 1 w is in a range from 20 mm to 160 mm and the interval is 20 mm. The value of w is 20 mm greater than 1 w . Only one experiment is carried out for each specimen. Each point in the forming limit diagram is only one sample. The stress-strain response of AA6016 and AA5182 is expressed in Figure 4. The stressstrain response are fitted by power function and the adjusted R-squared is greater than 0.995. The fitting results are expressed in Table 1 and Figure 4. The stress-strain response of AA6016 and AA5182 is expressed in Figure 4. The stressstrain response are fitted by power function and the adjusted R-squared is greater than 0.995. The fitting results are expressed in Table 1 and Figure 4.  Figure 4. The stress-strain response of AA6016 and AA5182.
The forming limit tests are carried out in accordance with GB/T15825.631995 [18]. The drawing of the test tooling is shown in Figure 5. The strain path (ratio of major to minor strain) is changed according to the change of specimen width and the lubrication condition. The test pieces are rectangular sheets shown in Figure 6. The width 1 w is in a range from 20 mm to 160 mm and the interval is 20 mm. The value of w is 20 mm greater than 1 w . Only one experiment is carried out for each specimen. Each point in the forming limit diagram is only one sample.  The forming limit tests are carried out in accordance with GB/T15825.631995 [18]. The drawing of the test tooling is shown in Figure 5. The strain path (ratio of major to minor strain) is changed according to the change of specimen width and the lubrication condition. The test pieces are rectangular sheets shown in Figure 6. The width w 1 is in a range from 20 mm to 160 mm and the interval is 20 mm. The value of w is 20 mm greater

Parameters' Effect on the Prediction of the FLC
From Equations (4) and (5), the yield criterion directly determines the stress and strain increment tensor of each iteration. Accordingly, the yield criterion has a direct impact on the prediction results. In part 2.1.1 (Theoretical Analyses), the Hill'48 yield criterion was applied with the normal anisotropy. When the normal anisotropic coefficient was equal to one, the Von Mises yield criterion was applied. For both of the above mentioned yield criteria, the initial thickness imperfection factor was 0.999.
The predicted FLCs of AA6016 and AA5182, based on the M-K method, are shown in Figure 7. The lowest point of the FLC, predicted by the M-K method, is independent of the normal anisotropic coefficient. This outcome is consistent with the previous theoretical derivations in part 2.1.3. When the normal anisotropic coefficient changes, there is a subtle change in the shape of the prediction curve at the left side (negative minor strain) while significant change is seen on that at the right side (positive minor strain). The curve rises with the decrease of the normal anisotropic coefficient.
From the analysis in part 2.1.1, the predicted results, based on the M-K method, are dependent on the hardening exponent, but independent on the strength coefficient. According to Table 1, the hardening exponent of 6016 is less than that of 5182. Therefore, when the hardening exponent increases, the predicted FLC based on M-K method moves up, as shown in Figure 7.
By comparing the experimental data of AA6016 and AA5182, the predicted FLC is independent on the strength coefficient under the M-K method. The FLC tends to move up when the strength coefficient increases. In addition, the 0 FLC relies on the hardening exponent.

Parameters' Effect on the Prediction of the FLC
From Equations (4) and (5), the yield criterion directly determines the stress and strain increment tensor of each iteration. Accordingly, the yield criterion has a direct impact on the prediction results. In part 2.1.1 (Theoretical Analyses), the Hill'48 yield criterion was applied with the normal anisotropy. When the normal anisotropic coefficient was equal to one, the Von Mises yield criterion was applied. For both of the above mentioned yield criteria, the initial thickness imperfection factor was 0.999.
The predicted FLCs of AA6016 and AA5182, based on the M-K method, are shown in Figure 7. The lowest point of the FLC, predicted by the M-K method, is independent of the normal anisotropic coefficient. This outcome is consistent with the previous theoretical derivations in part 2.1.3. When the normal anisotropic coefficient changes, there is a subtle change in the shape of the prediction curve at the left side (negative minor strain) while significant change is seen on that at the right side (positive minor strain). The curve rises with the decrease of the normal anisotropic coefficient.
From the analysis in part 2.1.1, the predicted results, based on the M-K method, are dependent on the hardening exponent, but independent on the strength coefficient. According to Table 1, the hardening exponent of 6016 is less than that of 5182. Therefore, when the hardening exponent increases, the predicted FLC based on M-K method moves up, as shown in Figure 7.
By comparing the experimental data of AA6016 and AA5182, the predicted FLC is independent on the strength coefficient under the M-K method. The FLC tends to move up when the strength coefficient increases. In addition, the 0 FLC relies on the hardening exponent.  (4) and (5), the yield criterion directly determines the stress and strain increment tensor of each iteration. Accordingly, the yield criterion has a direct impact on the prediction results. In Section 2.1.1 (Theoretical Analyses), the Hill'48 yield criterion was applied with the normal anisotropy. When the normal anisotropic coefficient was equal to one, the Von Mises yield criterion was applied. For both of the above mentioned yield criteria, the initial thickness imperfection factor was 0.999.

From Equations
The predicted FLCs of AA6016 and AA5182, based on the M-K method, are shown in Figure 7. The lowest point of the FLC, predicted by the M-K method, is independent of the normal anisotropic coefficient. This outcome is consistent with the previous theoretical derivations in Section 2.1.3. When the normal anisotropic coefficient changes, there is a subtle change in the shape of the prediction curve at the left side (negative minor strain) while significant change is seen on that at the right side (positive minor strain). The curve rises with the decrease of the normal anisotropic coefficient.  Table 2. Therefore, the iteration termination condition can be taken as 3 ϕ > . That is to say, 1  From the analysis in Section 2.1.1, the predicted results, based on the M-K method, are dependent on the hardening exponent, but independent on the strength coefficient. According to Table 1, the hardening exponent of 6016 is less than that of 5182. Therefore, when the hardening exponent increases, the predicted FLC based on M-K method moves up, as shown in Figure 7.
By comparing the experimental data of AA6016 and AA5182, the predicted FLC is independent on the strength coefficient under the M-K method. The FLC tends to move up when the strength coefficient increases. In addition, the FLC 0 relies on the hardening exponent.
Based on previous relevant studies, the value of ϕ T in the termination condition of the algorithm is in a range from 7 to 10. However, the principle of determining this value is unclear. This paper discusses the influence of the value on predicted FLC. Figure 8 reveals the ratios of equivalent strains in region B to those in region A. The values of α are chosen for the purpose that the strain paths are evenly distributed throughout the forming limit diagram. The value ϕ changes slowly at the initial stage and then increases rapidly when ε, the equivalent strain, reaches a certain value. To research the influence of ϕ T on predicted result, let ε ϕ T denote the equivalent strain when ϕ reaches ϕ T . If setting the values of ∆ε are very small under different stress paths. The ratios of ∆ε to ε 10 are less than 2%, which are listed in Table 2. Therefore, the iteration termination condition can be taken as ϕ > 3. That is to say, dε B 1 /dε A 1 > 3 is regarded as the termination condition of the algorithm under the M-K method. Based on previous relevant studies, the value of T ϕ in the termination condition of the algorithm is in a range from 7 to 10. However, the principle of determining this value is unclear. This paper discusses the influence of the value on predicted FLC. Figure 8 reveals the ratios of equivalent strains in region B to those in region A. The values of α are chosen for the purpose that the strain paths are evenly distributed throughout the forming limit diagram. The value ϕ changes slowly at the initial stage and then increases rapidly when ε , the equivalent strain, reaches a certain value. To research the influence of T ϕ on predicted result, let T ϕ ε denote the equivalent strain when ϕ reaches T ϕ . If setting the values of ε Δ are very small under different stress paths. The ratios of ε Δ to 10 ε are less than 2%, which are listed in Table 2. Therefore, the iteration termination condition can be taken as 3 ϕ > . That is to say, 1    In addition, the predicted FLCs for AA6016 under different initial thickness imperfection factors are expressed in Figure 9. The FLC curve calculated from the M-K theory is greatly influenced by initial thickness imperfection factor. The FLC tends to move up when the initial thickness imperfection increases. In addition, the  In addition, the predicted FLCs for AA6016 under different initial thickness imperfection factors are expressed in Figure 9. The FLC curve calculated from the M-K theory is greatly influenced by initial thickness imperfection factor. The FLC tends to move up when the initial thickness imperfection increases. In addition, the 0 FLC relies on the initial thickness imperfection factor.

Performance of the Improved M-K Model
According to the analysis in part 2.2, both Keeler's empirical formula and the M-K method depend on the hardening exponent. However, the M-K method depends on initial thickness imperfection factor, while Keeler's empirical formula depends on the thickness of the sheet metal. FLC experiments are carried out to explore the relationship between these parameters and the predicted FLC's accuracy. The experiments' results are stated in Figures 10 and 11. Meanwhile, the predicted FLCs, based on M-K method and Keeler's empirical formula, are presented.
For AA6016 and AA5182, the prediction results, based on the M-K method and Keeler's empirical formula, respectively, are very different. Under Keeler's empirical formula, the shape of the FLC never changes, while only the lowest point  FLC 0 relies on the initial thickness imperfection factor.

Performance of the Improved M-K Model
According to the analysis in Section 2.2, both Keeler's empirical formula and the M-K method depend on the hardening exponent. However, the M-K method depends on initial thickness imperfection factor, while Keeler's empirical formula depends on the thickness of the sheet metal. FLC experiments are carried out to explore the relationship between these parameters and the predicted FLC's accuracy. The experiments' results are stated in Figures 10 and 11. Meanwhile, the predicted FLCs, based on M-K method and Keeler's empirical formula, are presented. empirical formula is more consistent with the experimental results. Therefore, these two methods are suggested to be combined. Calculating 0 FLC under the M-K method, and then the FLC curve according to Equation (15), the FLCs are presented in Figures 10 and  11. By comparison of the experimental data and the predicted FLC under the proposed method, the predicted FLC under the proposed new method fits the experimental data very well for both AA6016 and AA5182 when the initial thickness imperfection factor is suitable.

Conclusions
In this study, the theoretical derivation of the M-K method is proposed. Besides, the relationship between the M-K method and Keeler's empirical formula is presented. Under this relation, a new method to predict FLC, which combines both of these two aspects, is proposed. Finally, the experiments about the forming limit curve of AA6016 and AA5182 validate the proposed new method. Consequently, it is concluded that: 1. By comparison with the experimental data of AA6016 and AA5182, the predicted FLC is independent on the strength coefficient under the M-K method. The FLC tends to move up when the strength coefficient or the initial thickness imperfection increases. In addition, the 0 FLC relies on the hardening exponent and the initial thickness imperfection factor. empirical formula is more consistent with the experimental results. Therefore, these two methods are suggested to be combined. Calculating 0 FLC under the M-K method, and then the FLC curve according to Equation (15), the FLCs are presented in Figures 10 and  11. By comparison of the experimental data and the predicted FLC under the proposed method, the predicted FLC under the proposed new method fits the experimental data very well for both AA6016 and AA5182 when the initial thickness imperfection factor is suitable.

Conclusions
In this study, the theoretical derivation of the M-K method is proposed. Besides, the relationship between the M-K method and Keeler's empirical formula is presented. Under this relation, a new method to predict FLC, which combines both of these two aspects, is proposed. Finally, the experiments about the forming limit curve of AA6016 and AA5182 validate the proposed new method. Consequently, it is concluded that: 1. By comparison with the experimental data of AA6016 and AA5182, the predicted FLC is independent on the strength coefficient under the M-K method. The FLC tends to move up when the strength coefficient or the initial thickness imperfection increases. In addition, the 0 FLC relies on the hardening exponent and the initial thickness imperfection factor. For AA6016 and AA5182, the prediction results, based on the M-K method and Keeler's empirical formula, respectively, are very different. Under Keeler's empirical formula, the shape of the FLC never changes, while only the lowest point FLC 0 does with the corresponding parameters. On the contrary, both FLC 0 and the shape of the FLC under the M-K method change with the parameters. Therefore, although there are some similarities between Keeler's empirical formula and the M-K method, their predicted results are quite different in both the lowest point of the FLC, FLC 0 , and the shape of the FLC.
Compared with experimental results, the value of FLC 0 under the M-K method, tends to be more accurate when the initial thickness imperfection factor is appropriate (for example, the thickness imperfection factor of AA6016 is between 0.9 and 0.999, while the factor of AA5182 is between 0.7 and 0.9). In contrast, the shape of the FLC under Keeler's empirical formula is more consistent with the experimental results. Therefore, these two methods are suggested to be combined. Calculating FLC 0 under the M-K method, and then the FLC curve according to Equation (15), the FLCs are presented in Figures 10  and 11. By comparison of the experimental data and the predicted FLC under the proposed method, the predicted FLC under the proposed new method fits the experimental data very well for both AA6016 and AA5182 when the initial thickness imperfection factor is suitable.

Conclusions
In this study, the theoretical derivation of the M-K method is proposed. Besides, the relationship between the M-K method and Keeler's empirical formula is presented. Under this relation, a new method to predict FLC, which combines both of these two