Numerical Material Testing Method for Hexagonal Close-Packed Metals Based on a Strain-Rate-Independent Finite Element Polycrystal Model

Abstract

The purpose of this study was to develop a numerical material testing method applicable to hexagonal close-packed (hcp) materials that can predict complex material behavior such as biaxial test results from relatively easy-to-perform uniaxial tests. The proposed numerical material testing method consists of a physical model that represents the macroscopic behavior of the material and a means of determining the included crystallographic parameters using macroscopic experimental data. First, as the physical model, the finite element polycrystal model (FEPM) previously applied by the authors for face-centered cubic (fcc) materials was applied and modified for hcp materials. A unique feature of the FEPM is that it avoids the use of strain-rate-dependent coefficients, whose physical meaning is ambiguous, because the deformation analysis can be performed while automatically determining the activity of all slip systems. The applicability of FEPM to numerical material testing methods was verified in hcp materials through this study. Then, a material parameter optimization process was developed using a genetic algorithm (GA). The proposed method was validated using literature values of magnesium alloy AZ31. First, the proposed optimization process was performed on cast AZ31 using uniaxial tensile and compressive stress—strain curves as teaching data to confirm that the stress—strain curves for the biaxial state could be predicted. Then, the proposed method was applied to rolled sheet AZ31, where the pseudo-anisotropic crystal orientations generated by numerical rolling were used as initial values. The prediction of unknown material data showed that, even in the case of sheets, the crystallographic parameters could be reasonably determined by the proposed optimization process.

1. Introduction

Due to their excellent properties such as their light weight, magnesium alloys are expected to be candidates for structural materials that can replace conventional materials used in engineering products, including automobiles and aircraft, where environmental friendliness is strongly required [1]. In particular, since significant ${\mathrm{C}\mathrm{O}}_2$ reduction regulations will be imposed on automobiles in the near future, the weight reduction of the car body, which has a direct impact, is considered to be a realistic policy from the viewpoint of cost effectiveness [2,3].

In general, magnesium alloys are often used in the manufacture of automotive components by casting, because they can be efficiently cast into complex shapes [4]. Although there are some features that need to be improved when they are used as sheet metal, such as their anisotropy and low ductility, there is the possibility of achieving a large weight reduction compared to steels and aluminum alloys by promoting their application in panel-like parts [5]. Therefore, it is of great engineering significance to understand the deformation characteristics of magnesium alloys and to model them appropriately. However, the complex deformation characteristics of magnesium alloys have not been sufficiently investigated, because they are subjected to plastic working less frequently than steels and aluminum alloys.

When conducting forming simulation, it is necessary to select an appropriate material model based on the yield function according to the application. Among such models, the most frequently used is the classical yield function proposed by von Mises [6]. It is applicable to many isotropic materials and, therefore, still plays an important role. For materials with strong in-plane anisotropy, such as aluminum alloy sheets, Hill’s yield function [7], which is a modification of the von Mises model, is often used. Various yield functions have subsequently been proposed for more accurate representation. For example, the Yld2000-2d model proposed by Barlat et al. [8] is often applied to anisotropic materials. Yld2000-2d has eight parameters, which are determined by solving simultaneous equations with multiple variables using eight experimental values. If the stress principal axis coincides with the anisotropic principal axis, the number of parameters decreases to six. In any case, yield stress in equi-biaxial tension as well as in-plane tensile tests are required. Although recent yield functions have shown high accuracy in practical use, the number of material parameters included in them has increased, and the burden of the material testing that they require has become non-negligible. For magnesium alloys, which are the focus of this paper, yield functions that take into account the strength differential (SD) effect are necessary, such as CPB06 [9] and that proposed by Verma et al. [10]. Many other yield function models have been proposed, but several advanced material models require not only uniaxial tensile and compression tests, but also biaxial material tests that necessitate specialized equipment and cannot easily be utilized in plastic forming simulations. To be more precise, the presence of friction and bulging in compression tests makes pure compressive material testing difficult. Biaxial testing is often performed using a cruciform testing machine or a hydraulic bulge tester, but these devices are not widespread. In addition, biaxial testing itself has problems to be considered in handling anisotropy and measuring plastic strain. Therefore, if at least some of the difficult-to-perform material tests, such as biaxial tests, could be replaced by numerical simulations, product design based on deformation analysis would be promoted, and the application of magnesium alloys would advance.

To accurately understand the deformation behavior of polycrystalline metals, a simulation model that takes the crystal structure into account is necessary. The phenomenological yield function and work-hardening models are not suitable for this purpose because they ignore microscopic mechanisms. Therefore, in recent years, crystal plasticity models have been frequently used for material modeling and deformation analysis to understand the material properties. There are various methodologies for the homogenization process, such as the visco-plastic self-consistent (VPSC) model, elasto-plastic self-consistent (EPSC) model, and elasto-vicoplastic self-consistent (EVPSC) model based on the self-consistent approach, and the crystal plasticity finite element model (CPFEM) based on finite element discretization. A comprehensive review was provided by Nguyen et al. [11] concerning research into computational methods based on various crystal plasticity models, from basic theory to recent applications. In modeling based on the finite element method, one crystal grain is discretized by one or more elements, so that the behavior of each crystal and grain boundary can be treated. In addition, because the common finite element code is applicable, there are many application examples. Its applications to numerical material testing by Kraska et al. [12], Zhang et al. [13], and Sedighiani et al. [14] are especially pertinent to this study. They used material test data such as tensile test results to determine material parameters in a virtual material constructed based on the CPFEM and predicted material properties such as biaxial test results. Considering the increasing accuracy requirements for forming simulations and the expanding application range of magnesium alloys in the future, it is necessary to increase the reliability of such numerical material testing methods.

The majority of crystal plasticity models use formulations based on the strain-rate-dependent constitutive model [15]. In this model, the strain rate sensitivity is used as a parameter to calculate the strain rate, which assumes that all slip systems are activated at the same time. This situation does not occur in actual material deformation; therefore, the strain rate sensitivity parameter is used to compensate for this inconsistency. In other words, the strain rate sensitivity parameter does not necessarily represent the mechanism of a slip system in a physical sense. Thus, the strain rate sensitivity parameter is not determined by material tests but empirically. To construct a physically accurate alternative model, Takahashi et al. [16,17] developed a finite element polycrystal model (FEPM) based on a successive accumulation scheme that directly deals with Schmid’s law. In the FEPM, a criterion introduced by Takahashi et al. [16] for determining the activity is applied to all slip systems, the slip increment is calculated only for those slip systems that are determined to be active, and the plastic strain increment is calculated as the sum of these. In this method, it is not necessary to determine the activity or inactivity of a slip system in advance; it is just determined as a result of convergence. Therefore, there is no need to introduce physically ambiguous parameters such as the strain rate sensitivity parameter. In addition, an advantage is that the degree of activity of each slip system can be monitored throughout the deformation. Although the successive accumulation-based crystal plasticity model has some advantages, it is not as well validated as the conventional strain-rate-dependent model. The FEPM was recently revalidated by Oya and Araki [18], and its usefulness in numerical material testing for fcc (face-centered cubic) aluminum alloys was demonstrated. These authors succeeded in predicting unexecuted material test data from known material test data by optimizing the material parameters using a genetic algorithm. This study demonstrated the utility of a crystal plasticity scheme that does not use the common strain-rate-dependent model, but it has not been applied to bcc (body-centered cubic) and hcp materials, which exhibit more complex behavior.

Next, the modeling of magnesium alloys is explained. Magnesium is a material with a hexagonal close-packed (hcp) crystal structure. Different slip and twin modes are considered in the modeling of hcp metals, each of which is assumed to have different critical resolved shear stress (CRSS) and hardening parameters. In modeling plastic deformation behavior, the determination of these microscopic parameters for each deformation mode is usually carried out so that the simulated flow stresses fit the experimental stress–strain curve. The two main deformations at room temperature, when undergoing compression along the basal plane or tension along the basal plan normal, are basal slippage and extension twinning. In fact, these two deformations are characterized by the lowest CRSS. Pyramidal and prismatic slip systems are more likely to be activated during tensile deformation in the basal plane and compression along the basal plane normal, respectively [19]. To achieve a reliable determination of the parameters of the model, a loading path that would guarantee the activation of all the relevant deformation modes [20] should be used. It is clear that the parameters are not unique [21], since only a limited number of stress–strain curves can be used for fitting, and thus the prediction of deformation along other strain paths is inaccurate. The CRSS for each deformation system varies from study to study, because it strongly depends on not only the material properties such as the microstructure, but also the strain path, as well as the model chosen [22]. The hardening laws of the CRSS for $\left\{10\overline{1}2\right\}$ twinning is surmised to be one of the reasons for this, because it is considered constant in some studies [22,23] and decreasing [24] or increasing [25] in others. The properties of $\left\{10\overline{1}2\right\}$ extension twins, such as active twin variants, are highly dependent on the strain path. $\left\{10\overline{1}2\right\}$ twinning can be incorporated into the model in a variety of ways, because it introduces partial lattice reorientation and new interactions between twin boundaries and slip modes.

In this study, a new numerical material testing method for magnesium alloys with a hcp structure was developed by applying the FEPM as a crystal plasticity solver. This paper presents the first application of the FEPM to a numerical material testing method for hcp-structured materials and introduces a modeling method and a parameter determination method using a genetic algorithm (GA) for its construction. In Section 2, a brief description of the FEPM and its application to hcp materials is presented. In Section 3, we explain the process of using the GA to determine the material parameters included in the hardening law of the slip system in the crystal plasticity model, which are unknown parameters, and the proposed numerical material testing method is used to predict the unknown material properties, employing experimental data of cast and rolled sheets of magnesium alloy AZ31 obtained from the literature. The validation results of the proposed method are discussed in Section 4, and the conclusions are presented in Section 5.

2. Material Modeling Based on the FEPM

2.1. Overview of the Finite Element Polycrystal Model

In this research, as presented originally by Takahashi et al. [16,17] and used in our previous work [18,26,27], the finite element polycrystal model (FEPM), which is a strain-rate-independent crystal plasticity model, was applied. A crystal is modeled as a cube consisting of five tetrahedral element, and N crystals with different orientations form the polycrystalline aggregate. Unlike the strain-rate-dependent crystal plastic constitutive equation, the FEPM calculation proceeds by determining the activity of all slip systems. The plastic strain component is calculated by integrating the slip increments of the determined active slip systems, and the finite element analysis proceeds. The finite element formulation based on the initial strain method is a general, as also used in [16,17], and is not explained in detail in this paper. A brief description of the numerical calculation procedures is provided in Appendix A.

2.2. Expression of the Plastic Strain Increment

For the finite element formulation, it is necessary to calculate the microscopic plastic strain increment of an element. To model deformation based on slip and twinning systems, the following formulation [28] was used in this study:

$${\dot{\mathit{\epsilon }}}^p=(1-\sum _{\beta =1}^{N^{\mathrm{tw}}}{g}^{\beta })\sum _{\alpha =1}^{N^{\mathrm{sl}}}{\dot{\gamma }}^{\alpha }{\mathit{L}}_{\mathrm{sl}}^{\alpha }+\sum _{\beta =1}^{N^{\mathrm{tw}}}{g}^{\beta }{\dot{\gamma }}_{\mathrm{tw}}^{\beta }{\mathit{L}}_{\mathrm{tw}}^{\beta },$$

(1)

where $\mathit{L}$ is the plastic velocity gradient, which is additively decomposed into the slip velocity gradient ${\mathit{L}}_{\mathrm{sl}}$ and twinning velocity gradient ${\mathit{L}}_{\mathrm{tw}}$; numbers $N^{\mathrm{s}}$ and $N^{\mathrm{tw}}$ denote the number of active slip and twinning systems; ${\dot{\gamma }}^{\alpha }$ and ${\dot{\gamma }}_{\mathrm{tw}}^{\alpha }$ denote the shear strain increment on slip system $\alpha$ and twinning system $\beta$, respectively; and ${\sum }_{\beta =1}^{N^{\mathrm{tw}}}{g}^{\beta }$ represents the total volume fraction of the deformation twinned region in the grain.

For the calculation of the plastic strain increment, the slip increment for each slip system must be obtained explicitly. However, this is not easy due to the lack of constitutive laws for single crystals, and various methods have been investigated. Currently, the well-known strain-rate-dependent model [15] is often used:

$$\begin{array}{c}\hfill {\dot{\gamma }}^{\alpha }={\dot{\gamma }}_0{\left(\frac{{\tau }^{\alpha }}{{\tau }_c^{\alpha }}\right)}^{\frac{1}{m}}\mathrm{sign}\left({\tau }^{\alpha }\right),\end{array}$$

(2)

where ${\dot{\gamma }}_0$ is the reference shear strain, ${\tau }^{\alpha }$ is the working resolved shear stress, ${\tau }_c^{\alpha }$ is the CRSS, and m is the strain rate sensitivity parameter. The model assumes that all slip systems are active, and the parameter m adjusts for discrepancies with the actual degree of activity. While this model is efficient, it is not physically accurate, and the results depend on the value of the parameter m, which is determined empirically.

2.3. Determination of Shear Strain Rate by the Successive Accumulation Method

Once the plastic strain increment of each element is determined, virtual external forces are obtained, and the stiffness equation of the FEM can be solved. The plastic strain increment ${\dot{\epsilon }}_{ij}^p$ can be calculated by

$${\dot{\epsilon }}_{ij}^p=\sum _r{L}_{ij}^{\alpha }{\dot{\gamma }}^{\alpha },L_{ij}^{\alpha }=(a_i^{\alpha }{b}_j^{\alpha }+a_j^{\alpha }{b}_i^{\alpha })/2,$$

(3)

where $\dot{\gamma }$ is the shear strain increment and $L_{ij}$ is the Schmid tensor component, which is expressed by $a_i$, the unit vector normal to a slip plane, and $b_i$, the unit vector in the slip direction.

The resolved shear stress on the $\alpha$th slip system can be determined by

$${\tau }^{\alpha }=L_{ij}{\sigma }_{ij}.$$

(4)

The yield condition of a slip system is expressed by

$$\begin{array}{ccc}\hfill {\zeta }^{\alpha }{\dot{\gamma }}^{\alpha }\ge 0& \mathrm{if}& \hfill {\zeta }^{\alpha }{\tau }^{\alpha }\ge {\tau }_c^{\alpha },\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\zeta }^{\alpha }{\dot{\gamma }}^{\alpha }=0& \mathrm{if}& \hfill {\zeta }^{\alpha }{\tau }^r<{\tau }_c^{\alpha },\hfill \end{array}$$

(6)

where ${\zeta }^{\alpha }$ is defined as ${\zeta }^{\alpha }=\mathrm{sign}\left({\tau }^{\alpha }\right)$, and ${\tau }_c^{\alpha }$ is the CRSS on the $\alpha$th slip system that is expressed by a certain hardening law.

An unknown quantity, the shear strain increment ${\dot{\gamma }}^{\alpha }$, needs to be determined, but this is made difficult by the fact that the combination of active slip systems is not predetermined. Takahashi solved this problem by introducing the successive accumulation method in order to reproduce the actual slip phenomenon. Takahashi assumed the following: When a load increment is applied, the material first responds elastically. Subsequently, the shear stress in the active slip system exceeds the yield stress, which acts as a driving force to initiate dislocation movement. When even one grain is plastically deformed, the balance of forces in the entire aggregate changes, and the stress in the plastically deformed grain relaxes. On the other hand, the yield shear stress increases due to work-hardening. When the decrease in applied stress and increase in deformation resistance are balanced, the slippage within the crystal stops. The successive accumulation method is an algorithm for numerically reproducing this process.

For all grains and slip systems, the following calculation is conducted:

$$\begin{array}{c}\hfill {\chi }_i^{\alpha }={\chi }_{i-1}^{\alpha }+\left({\zeta }^{\alpha }{\tau }^{\alpha }-{\tau }_c^{\alpha }\right)\frac{\Delta \rho }{2G},\end{array}$$

(7)

where i is a repetition counter, and $\Delta \rho$ determines the repetition step. This evaluation begins with an initial value ${\chi }_0^{\alpha }=0$. Then, the equivalent external force ${\mathit{F}}^p$ is updated at each iteration step, and the stiffness equation is solved using this to obtain the plastic strain. During the iterative process, the following condition is imposed: If ${\chi }^{\alpha }<0$ occurs during the iterations, set ${\chi }^{\alpha }=0$ for the slip system and continue the entire iteration. When the value of ${\chi }^{\alpha }$ no longer changes for all slip systems and all grains, the solution is fixed; ${\chi }^{\alpha }$ converges to a positive value or 0, corresponding to active and inactive slip systems, respectively. Therefore, this method eliminates the need to determine in advance whether each slip system is active or inactive, which is automatically determined by the convergence results during the calculation.

2.4. Hardening Model for Slip System

In order to proceed with the calculations using the methods described so far, a hardening law is needed to update the CRSS in Equation (7). The work-hardening expression can be represented by the sum of the self- and latent hardening that is caused by the interaction among different slip/twin systems. The increment of the threshold shear stress ${\dot{\tau }}_c^{\alpha }$ is expressed by

$${\dot{\tau }}_c^{\alpha }=\frac{\mathrm{d}{\tau }_c^{\alpha }}{\mathrm{d}\Gamma }\sum _{\beta }{h}^{\alpha \beta }\left|{\dot{\gamma }}^{\beta }\right|,$$

(8)

where $\Gamma ={\sum }_{\beta }\left|{\dot{\gamma }}^{\beta }\right|$ represents the accumulated shear strain in the grain, and $h^{\alpha \beta }$ provides empirical coefficients for the latent hardening to express the obstacles on system $\alpha$ associated with system $\beta$. The threshold stress ${\tau }_c^{\alpha }$ is characterized by a Voce hardening law [22,29]:

$${\tau }_c^{\alpha }={\tau }_0^{\alpha }+({\tau }_1^{\alpha }+{\theta }_1\Gamma )\left(1-\exp \left(-\frac{{\theta }_0}{{\tau }_1^{\alpha }}\Gamma \right)\right)$$

(9)

where ${\tau }_0$, ${\tau }_1$, ${\theta }_0$, and ${\theta }_1$ are the initial and back-extrapolated increase in the threshold stress and the initial and final slope of the hardening curve, respectively.

2.5. Formulation of Twinning

Magnesium alloys do not have a sufficient number of independent slip systems to accumulate arbitrary deformations. Therefore, especially at room temperature, twinning is an important deformation mechanism that requires special treatment in the simulation of magnesium alloy deformation. In this study, the predominant twin reorientation (PTR) scheme [28] was applied to formulate the plastic strain increment. The model based on the PTR scheme requires the evaluation of the volume fraction of the twin deformation. The following formulation was used, based on Tang et al. [30].

In each deformation step, the incremental twin volume fraction of the grain n associated with each twinning system $\beta$ is

$$\Delta g^{n,\beta }=\frac{\Delta {\gamma }^{n,\beta }}{{\gamma }^{\beta }}.$$

(10)

The fraction of twinning is then accumulated after each deformation step:

$$g^{n,\beta }=\sum _{steps}\Delta g^{n,\beta },$$

(11)

where $\Delta {\gamma }^{n,\beta }$ is the shear strain increment obtained by the contribution of the twinning system $\beta$ in the grain n, and ${\gamma }^{\beta }$ is the characteristic shear strain associated with twinning, whose value is a constant of 0.13 [31]. As a simplified version of the procedure presented in the literature [30], the twinned volume fraction $f^n$ of a grain n can be defined by summation over the whole twinning system:

$$f^n=\sum _{\beta }^{N_{tw}}{g}^{n,\beta }.$$

(12)

Although it can be a variable dependent on the twinning process, in this study a threshold value for reorientation was defined as an empirical constant [27]:

$$F_T=0.9.$$

(13)

The constant was set at this value to better characterize the sudden hardening of the model due to the complete reorientation of the lattice via twinning. The twinning fraction $f^n$ of each grain is compared to the threshold $F^T$ at each deformation step. If $f^n$ is greater than $F_T$, the grain is completely reoriented toward the twin, and the twin system stops growing in the next step.

$${\dot{\gamma }}^{n,\beta }=\left\{\begin{array}{ccc}{\dot{\gamma }}^{n,\beta }& \mathrm{if}& f^n<F_T\\ 0& \mathrm{if}& f^n\ge F_T\end{array}.\right.$$

(14)

3. Numerical Material Testing Method and Its Demonstration

3.1. Preliminary Analysis Using Cast AZ31

To confirm the applicability of the proposed model, it was necessary to compare the simulation results with experimental data. In particular, to evaluate the effectiveness of the model, the model was fit and optimized to the results of uniaxial tension and compression experiments, and then the curve for biaxial compression was predicted as a demonstration of the proposed numerical material testing method. Experimental data were obtained from the literature [32]. In the referenced paper, cast AZ31 magnesium alloy (2.90 wt% Al, 0.89 wt% Zn) was supplied as 89 mm diameter cylindrical rods. The specimens were annealed at 523 K for 2 h, and a random microstructure was observed. The average grain size of the cast alloy was about 180 μm. Tensile specimens were machined from the rods in $7\times 7{\mathrm{mm}}^2$ square cross-sections. Cubic specimens with edge lengths of 7 mm were also machined for uniaxial and equi-biaxial compression tests. All surfaces of the specimens were abrasive-finished before mechanical testing. Uniaxial tensile tests were performed at a crosshead speed of 5 mm/min. Uniaxial and biaxial compression tests were performed at an average strain rate of $5\times {10}^{-4}$/s. In the compression tests, silicone grease mixed with boron nitride was applied to the die/specimen interfaces to reduce friction. All tests were conducted at room temperature.

A polycrystal model consisting of 1000 elements (grains) was used to simulate mechanical deformation. Figure 1 shows the polar figures of the texture of the cast material used in the simulation. Since the cast material was expected to exhibit isotropic behavior, the crystal orientation was randomly assigned to simulate the properties of the cast material. The random orientation generated was the same as when we previously validated another model [27].

Considering the lack of studies on constitutive models of cast magnesium alloys, the initial values to start fitting the curves were initially chosen by trial and error to achieve the best fit with the experimental data, as shown in

Table 1. Hardening parameters used to produce curves shown in Figure 2. These values were chosen as the initial values of the investigation to reproduce the behavior of the cast AZ31.

Mode	${\mathit{\tau }}_0$ MPa	${\mathit{\tau }}_1$ MPa	${\mathit{\theta }}_0$ MPa	${\mathit{\theta }}_1$ MPa
Basal	10.0	3.0	10.0	0
Prismatic	25.0	5.0	25.0	0
Pyramidal	50.0	9.0	35.0	0
Twin	20.0	8.0	15.0	0

. On the other hand, the latent hardening parameters were considered constant during the optimization, as shown in

Table 2. Empirically chosen self- and latent hardening parameters. These values were constant during the optimization process.

${\mathit{h}}^{\mathit{\alpha }\mathit{\beta }}$	Basal	Prismatic	Pyramidal	Twin
Basal	1.0	0.2	0.2	1.0
Prismatic	0.2	1.0	0.2	1.0
Pyramidal	0.2	0.2	1.0	1.0
Twin	0.3	0.3	0.3	1.0

. The initial results using these parameters are shown in Figure 2. As can be seen from Figure 2, the simulated stress–strain curves were close to the experimental data, but the accuracy was not sufficient. Therefore, it was necessary to adapt the microscopic parameters through an optimization process.

3.2. Material Parameter Optimization Procedure Based on a Genetic Algorithm for Cast AZ31

To determine the parameters of the crystal plasticity model, we constructed an optimization scheme using a genetic algorithm. The following strategies were considered in the construction of the optimization scheme: Since optimizing 12 parameters at once would be time-consuming and iterative, and the risk of obtaining local minima would be high, it was necessary to find a way to optimize fewer parameters at each step. Also, by dividing the optimization range, the error during the computation could be reduced in a non-monotonic manner, so that the local minimum was not reached. Therefore, we examined the influence of each parameter independently and investigated the regions in the flow curve that were highly related to each parameter. In fact, since ${\tau }_0$ had the greatest influence in the early stage of deformation, the values were optimized to just fit the simulation and the experimental results in the range of $0\le \epsilon \le 0.05$. ${\theta }_0$ had the greatest influence after yielding occurred during deformation, and ${\tau }_1$ was found to have a greater impact at the end of the deformation. Figure 3 summarizes these strategies. Following this scheme had another advantage. It did not require all simulations to be performed until the end, but only up to the point of interest for parameter-based optimization, thus reducing the computation time. Figure 4 shows a more detailed flowchart of the chromosome representation and the optimization process. As shown in this figure, the parameters corresponding to each slip or twin system were grouped by their type and optimized separately in three stages. In one optimization process of the proposed method, 48 chromosomes were used, and the total number of generations was set to 80. The reason for this setting was that a larger number of generations and chromosomes increased the computation time and did not necessarily improve the accuracy of the simulations. The fitness function in the GA was the sum of the errors in the simulated and experimental curves for each chromosome, calculated as the root mean square error. Since optimization by GA does not always converge to a physically correct solution no matter how many steps are taken, it is necessary to select appropriate initial values. When searching for parameters that represent material properties, as in numerical material testing, starting optimization from random initial values almost always results in failure, so appropriate parameters were determined by trial and error, such that the shapes of the two curves were close.

For the proposed method, the results of the tensile and compression uniaxial tests were used as teaching data, and the results of the equi-biaxial compression tests were used to validate the prediction results. The obtained hardening parameters are shown in

Table 3. Hardening parameters for the cast AZ31 obtained with the optimization process.

Mode	${\mathit{\tau }}_0$ MPa	${\mathit{\tau }}_1$ MPa	${\mathit{\theta }}_0$ MPa	${\mathit{\theta }}_1$ MPa
Basal	10.9	3.6	10.2	0
Prismatic	25.0	5.8	19.0	0
Pyramidal	48.7	8.7	33.4	0
Twin	25.2	8.1	11.9	0

. The stress–strain curves calculated using the determined parameters are shown in Figure 5. Note that simulations were performed for both compression and tensile tests in all three axes $(x,y,z)$, but due to the isotropy of the target material, only the results of the tensile test along the z axis are shown. Comparing the results of Figure 2 and Figure 5, it is clear that both the reproduction of the curves used as teaching data and the prediction results of the stress–strain curve for biaxial compression improved considerably, indicating that the optimization scheme constructed was effective.

3.3. Analyses Using Rolled Sheet of AZ31

For the modeling of the rolled sheet of magnesium AZ31, we used the data from [33] and compared them with the simulation results. The initial material was magnesium AZ31 (3% Al, 1% Zn, 0.6% Mn) in the form of a 60 mm thick rolled sheet, annealed at 400 °C for 2 h. Compressive blocks with dimensions of $15\times 10\times 12$ mm and dog-bone-type tensile specimens with a gauge length of 20 mm, cross-section of $7\times 1$ mm, and total length of 45 mm were spark-cut in the following directions: rolling direction (RD, corresponding to the x axis); transverse direction (TD, corresponding to the y axis); normal direction (ND, corresponding to the z axis); and 45° diagonal direction between the RD and ND. Tensile and compression tests were performed at an initial strain rate of 0.0005 s⁻¹, and each type of specimen was tested three times at room temperature. Initially, the crystalline texture was numerically obtained randomly, as shown in Figure 6. Then, an anisotropic texture was obtained by simulating in-plane tension along the RD until reaching Euler angles comparable to pole figures of the texture measured by the EBSD procedure. This was obtained because the texture of a rolled sheet has a strong basal orientation. Moreover, during the compression of the polycrystals of magnesium alloys, the c axes of the crystals composing the material tend to align closer to the direction of the compression (in this case, the compression was along the c axes, corresponding to the ND), as seen in [34]. This provided a composition of a similar texture after rolling. The pole figures of texture used in the following simulations are shown in Figure 7.

The self- and latent hardening matrix coefficients used are presented in

Table 4. Empirically chosen self- and latent hardening matrix coefficients for rolled sheet AZ31.

${\mathit{h}}^{\mathit{\alpha }\mathit{\beta }}$	Basal	Prismatic	Pyramidal	Twin
Basal	1.0	1.0	1.0	0
Prismatic	1.0	1.0	1.0	0
Pyramidal	1.0	1.0	1.0	0
Twin	0	0	0	1.0

. Due to the lack of a method to measure latent hardening between slips, all other slip system hardening coefficients $h^{\alpha \alpha }$ were set to 1, the same value used for self-hardening. A study by Bertin et al. [35], who applied discrete dislocation dynamics simulations to pure magnesium, reported that different latent hardening factors must be used. In particular, a basal slip is strongly hardened by pyramidal dislocations, whereas pyramidal dislocations are barely hardened by other slip modes. However, these results seem to be incompatible with the parameters determined in this study (explained in more detail later), in which ${\tau }_{\mathrm{basal}}$ slowly increased, whereas ${\tau }_{\mathrm{pyramidal}}$ rapidly increased. According to the literature, the latent hardening parameters between the twin and the slip systems are usually given as $h_{\alpha \beta }=1.2$, and twins induce the hardening of the slip systems compared to the base material. This hardening is considered to come from twin boundaries acting as barriers to dislocation propagation [36] and dislocation transmutation [24]. Most simulations carried out using crystal plasticity models have also shown the effect of higher latent hardening coefficients for slips caused by $\left\{10\overline{1}2\right\}$ twinning. The latent hardening induced by the twin to slip system is an adjustment to the hardening indirectly resulting from the crystal rotation induced by Hall–Petch effect hardening and the Basinski hardening mechanism. Here, the Hall–Petch effect refers to the effective grain size refinement and the reduction of the dislocation mean free path, and the Basinski hardening mechanism predicts an increment in strength that develops as a result of dislocation transmutations. The utilization of these mechanisms for the hcp crystal plasticity model was mentioned in the work of Patel et al. [37], who introduced the dislocation density model into the VPSC model. In this research, however, the crystal rotation of the bulk caused by the twinning was considered sufficient to simulate these hardening effects, keeping the latent hardening coefficient for twinning at 0.

3.4. Optimization of Microscopic Parameters Using GA for Rolled Sheet AZ31

Our goal was to predict the deformation behavior along the ND through thickness compression corresponding to the biaxial stress condition [38] using compression and tensile tests along the RD and TD. A similar optimization scheme to that designed for the cast material was developed and implemented for the rolled sheet material. The initial values were determined so that the simulation results were close to the experimental results, using the values in Table 3 as a reference. The determined values are shown in

Table 5. Hardening parameters as initial values for the optimization of the rolled sheet AZ31.

Mode	${\mathit{\tau }}_0$ MPa	${\mathit{\tau }}_1$ MPa	${\mathit{\theta }}_0$ MPa	${\mathit{\theta }}_1$ MPa
Basal	15.0	1.0	9.0	7.0
Prismatic	65.0	1.9	36.7	0
Pyramidal	80.0	5.2	57.0	0
Twin	27.0	4.0	10.0	0

, and the obtained curves are shown in Figure 8. Although the simulated and experimental values agreed in their trends, the results presented in Figure 8 were not quantitative enough. Therefore, as in the case of the cast material, further optimization was performed using the values in Table 5 as initial values. Again, a GA scheme was designed for this optimization process, as shown in Figure 9, in which 180 generations and 48 chromosomes were used. The strategy of this scheme is explained below.

As can be seen in Figure 8, the overall fit of the curve was achieved, so it was not appropriate to change the CRSS, which had a significant effect on the overall shape of the curve, and a constant value was assumed for ${\tau }_0$. Then, the optimization process was split into different slip/twin systems. The optimization of the parameters of each slip system was conducted separately with only the curves where the considered slip/twin was active as the teaching curve. This was because optimizing 12 parameters at once would take a considerable amount of computing time and would not necessarily converge to a physically significant solution, since reaching a local minimum is easier than reaching a global minimum. In addition to dividing, the order of the slips to be optimized was also important. That is, as shown in Figure 9, the slip systems that were active on fewer curves were optimized first, taking into account their relative activity, and then the slip systems that affected all curves, especially the basal slips, were optimized last. By separating the optimization process for each slip system in this way, it was physically possible to converge to a reasonable solution rather than optimizing all variables at once.

In order to perform proper optimization using a GA, the fitness criterion must be well defined. In this study, the fitness function was defined as a linear combination of the least-squares errors between experimental data and simulation results for the four stress–strain curves, as shown in the following equation:

$${\mathrm{Err}}_{\mathrm{total}}={\alpha }_1{\mathrm{Err}}_{\mathrm{RD}-\mathrm{T}}+{\alpha }_2{\mathrm{Err}}_{\mathrm{RD}-\mathrm{C}}+{\alpha }_3{\mathrm{Err}}_{\mathrm{TD}-\mathrm{T}}+{\alpha }_4{\mathrm{Err}}_{\mathrm{TD}-\mathrm{C}},$$

(15)

where Err represents the least-squares error, and its subscript indicates the type of curve under consideration. For example, RD-T refers to the tension in the RD. ${\alpha }_i$ are the coupling coefficients. Although these could be adjusted appropriately, they were found not to be highly sensitive in the present study, so they were simply all set equal at ${\alpha }_1={\alpha }_2={\alpha }_3={\alpha }_4=0.25$.

The parameters obtained by the proposed optimization method are shown in

Table 6. Hardening parameters of the rolled sheet AZ31 obtained by the proposed optimization procedure.

Mode	${\mathit{\tau }}_0$ MPa	${\mathit{\tau }}_1$ MPa	${\mathit{\theta }}_0$ MPa	${\mathit{\theta }}_1$ MPa
Basal	15.0	0.8	22.9	6.7
Prismatic	65.0	2.3	25.9	0
Pyramidal	80.0	3.5	132.0	0
Twin	27.0	0.8	7.5	0

, and the stress–strain curves drawn using them are shown in Figure 10. As can be seen from Figure 10, there was a clear improvement in accuracy compared to Figure 8, where no optimization was performed.

4. Discussion

The following is a discussion of the investigations that were conducted on cast and rolled AZ31 specimens. The difference between the two materials was that the cast material was isotropic, i.e., had a random texture, while the rolled material was anisotropic and had a strongly oriented texture. In the case of either material, since it is not easy to determine the multiple microscopic parameters of a polycrystalline model that are not directly measurable, some learning or optimization process was required. Even with a genetic algorithm, there was no certainty that a perfect optimization could be achieved, and the reliability of the results was highly dependent on the initial values and optimization scheme. In addition, a sufficient number of experimental data must be available to perform reliable optimization, but material tests other than uniaxial tensile tests are not common. This problem is particularly acute for hcp materials, where several different types of slip systems are active. In addition to the usual uniaxial tensile tests, biaxial tensile and compression tests may also be required, making material testing even more difficult. Therefore, if the numerical material testing method proposed in this paper can predict experimental data for material tests that are difficult to perform, the benefits would be especially great for hcp materials.

Although this study dealt with two types of materials, cast and rolled, from the viewpoint of the importance of deformation analysis, the main target was the prediction of the material properties of rolled sheet materials, which are subjected to various press forming operations. The crystal plasticity model used in this study contained 12 hardening parameters. It would be very difficult to determine all of them to fit the behavior of rolled AZ31. Therefore, we first estimated the parameters of cast AZ31, whose material properties are considered isotropic. Although the initial optimization variables were obtained by trial and error, judging from the stress–strain curves, the determined parameters were reasonable. Therefore, cast AZ31, although an hcp material, had relatively simple properties and was an appropriate first optimization target.

Next, using the parameters obtained for cast AZ31 and applying them directly to rolled AZ31, errors were naturally observed. Therefore, we designed another optimization scheme and performed optimization using the initial parameters separately prepared. In this scheme, pseudo-anisotropy data obtained by numerical rolling were used for the crystallographic orientation. In addition, the treatment of twinning systems in this study should also be explained. As was already described, latent hardening coefficients were considered necessary to explain indirect hardening due to crystal rotation by twins. However, in doing so, the stresses were overestimated. This resulted in significantly lower CRSS values for pyramidal and prismatic slips and the underestimation of stresses during TTC. From these results, it could be said that the hardening effect of twinning on the other slips in the proposed model could be fully explained only by the lattice rotation during deformation. These efforts improved the simulation results of all the stress–strain curves so that they were close to the experimental results in the range of up to $\epsilon =0.1$. These results indicate that the polycrystalline model and its microscopic parameter optimization scheme have physical validity and are reliable as a numerical material testing method.

It should also be noted that the strain-rate-dependent coefficient was not used in the verification process described above. This meant that in the deformation analysis conducted in this study and in the optimization process using it, the activity of the slip systems was evaluated without any discrepancy between the analytical and experimental results, which is an example that proves the usefulness of the FEPM with less physical ambiguity.

However, errors were not entirely absent, and the largest gap was in the ND-T prediction curve. Pyramidal slip behavior is not very relevant in uniaxial stress, but it becomes a major slip system under biaxial stress conditions, such as thickness direction compression (NT-D). Therefore, one of the reasons for this error may be that there is no teaching curve that can optimize only the parameters related to it. In addition, the treatment of latent hardening requires further validation. Since the direct experimental verification of latent hardening is very difficult, the validity of the proposed method could be enhanced by applying it to more macroscopic material test data.

5. Conclusions

The FEPM applied in this study was based on a successive accumulation algorithm, which conducted the deformation analysis while determining the activity of each slip system. This eliminated the need to use the physically ambiguous parameter of strain rate sensitivity and reduced the empirical aspect of crystal plasticity analysis. In this study, the FEPM was extended to be applicable to hcp materials, which was then applied to a numerical material testing method. The proposed method was used to represent and analyze the mechanical behavior of cast and rolled sheets of magnesium alloy AZ31.

The validation of the proposed method using literature values demonstrated the usefulness of the FEPM model and the numerical material testing method based on the proposed parameter determination method for hcp materials. These results suggest that the macroscopic parameters required for FEM-based forming simulation can be obtained by numerical simulation on a computer instead of conducting actual material tests. It should be emphasized that EBSD observations were not performed to obtain these results, and it can be assumed that the proposed optimization scheme successfully estimated the appropriate crystal orientation. Future work is needed to investigate different loading path cases with additional de-twinning modes and to further investigate other hcp material data to verify the applicability of the FEPM.