Crystal Plasticity Modeling of Al Alloy under Linear and Non-Linear Loading

The crystal plasticity (CP) model is widely used in many applications to link microstructure and mechanical properties. There are varying CP constitutive laws with phenomenological or physical-based formulation to cover a large range of loading conditions. In order to predict the deformation behavior of an Al alloy during the sheet metal forming process with either linear or non-linear strain path, both phenomenological and physical-based CP constitutive laws have been chosen, and the prediction performance of both models is compared. For the linear loading condition, the uniaxial tensile tests are performed on the smooth-dog-bone (SDB) specimens along rolling and transverse directions (RD/TD). The non-linear strain path is achieved by the Marciniak testing followed by uniaxial tension. In the first stage, the Marciniak testing is performed under the stress states of RD-uniaxial, plane strain, and biaxial tension. After being loaded to a certain strain level, mini-SDB specimens are cut along RD and TD from the uniform deformation region and reloaded under RD-uniaxial tension. The digital image correlation (DIC) technique is employed to measure the strain during testing. The electron backscatter diffraction (EBSD) technique is used to characterize the initial microstructure as well as the microstructure evolution of the specimens after the first stage loading in the non-linear strain path. A phenomenological power law and a dislocation-density-based hardening law have been employed in this study. The parameters are calibrated based on the flow curve of the RD uniaxial tension. The model performance is validated by stress–strain response under all the rest loading conditions including the non-linear loading path.


Introduction
In the last years, aluminum alloys have become widely used in the automotive industry. The main advantages of these materials are a high strength-to-weight ratio combined with good corrosion resistance. Aluminum is increasingly used in structural parts, body panels and various other components, and AA5754 alloy was chosen as the material for this study [1]. It is well known that the mechanical properties of metal materials are determined by the microstructure. Therefore, learning about the quantitative relationship between the microstructure and mechanical properties is very important [2]. Particularly, metal will be subjected to various linear or non-linear loading in the process of forming and application, which will cause different stress states to it and change its mechanical properties [3]. Giving materials different loading and studying the response behavior is of high interest for the sheet metal production and forming industry.
There has been a large number of studies focusing on bridging the microstructure and macroscopic mechanical properties by both experimental and numerical methods. The interested microstructural features [4] normally are the phase fraction, grain morphology, crystal orientation, secondary phase morphology, etc. While the focused mechanical properties can be various from microscopic level strain partition, macroscopic flow behavior, to fracture and fatigue properties, etc. For example, the effect of texture and dislocation structure on strain hardening anisotropy of aluminum alloy was studied by Lopes, et al [5]. J. W. Yoon, et al. [6] studied finite element simulations of the simple shear test were conducted for 1050-O and 6022-T4 aluminum alloy sheet samples. J. Wang, et al. [7] developed a crystal plasticity finite element model accounting for the microstructural features for simulating the fretting fatigue of AA7075-T651. However, for aluminum alloys that can be considered as a single phase, the challenges that hinder this study are mainly from the following two aspects: i) necessary to synthesize fine microstructure models, accurately characterize the microstructure of materials and accurately control the microstructure variables, and ii) two models are used to conduct efficient and reliable calibration of aluminum alloy under different stress conditions. Lian, et al. [8] developed a calibration strategy for synthetic microstructure and crystal plasticity parameters for producing fine grain structure aluminum alloys, which could systematically and quantitatively analyze the influence of microstructure characteristics.
This study aims to validate these complex pathways using the current formulation of either phenomenological constitutive models or the mechanism-based crystal plasticity models. We are intending to extend the experimental envelope to more general dimensions by combining the change of stress states and the loading angels, which are currently separated into two different research lines. This would create more general conditions that the data generated could be used to further develop, calibrate, and validate both the phenomenological models and mechanism-based models. Further analysis of the microstructure changes in this direction is also helpful to understand the plastic deformation mechanism from the microstructure and dislocation levels. A general research strategy is to decompose a complex path into several separate strain paths to simplify the problem [9]. Considering the mechanical properties of the material, the stress state and the loading direction are considered in the loading process. Generally, the loading design of the two stages is as follows. In the first stage, uniaxial, plane strain, and biaxial pre-stretching can be carried out on the specimen, and then small uniaxial tensile specimens are cut from the specimen after initial loading along different loading directions for the second step of the tensile test.

Material characterization
The material used in this study is a cold-rolled AA5754 H111 Aluminum alloy sheet of 1.5 mm thickness. H111 stands for the tempering type, and it recognizes that the alloy underwent some amount of cold strain hardening after annealing. The electron backscattered diffraction (EBSD) technique has been used to characterize the microstructure of the material, following the strategy used by Liu et al. [10]. The EBSD measurements were performed at 20 kV. The microstructure features like grain size, shape, and texture have been characterized. Samples have been obtained on two planes, i.e. the rolling direction (RD)-transverse direction (TD) plane and the RD-normal direction (ND) plane. In addition, varying investigation areas were chosen through thickness on the RD-ND plane to cover the possible thickness gradient. Grain reconstruction and analysis have been carried out by MATLAB/MTEX toolbox [11]. 15º misorientation threshold was used for the grain reconstruction. The investigated material has an average grain size of 16 μm and mean grain shape aspect ratio of 1:0.681:0.662. The typical fcc rolling texture is observed with a texture index of 1.24. Detailed information microstructure characterization methods and results is referred to Lian et al. [8].

Mechanical testing program
In this study, two categories of loading path experiments were designed, i.e. linear loading and non-linear loading respectively. Linear loading experiment for anisotropy investigation including the tensile test along different loading angles with respect to RD, i.e. 0º, 15º, 30º, 45º, 60º, 75º, and 90º, as shown in Fig. 1 (a). Quasi-static loading is applied with a constant strain rate of 10 -4 s -1 . For the non-linear strain paths loading experiment, the material has been firstly pre-strained under three stress states (i.e. RD-uniaxial, plane strain, and biaxial tension) using the hydraulic press and the Marciniak punch, as shown in Fig. 1 (b) and (c). After being loaded to a certain strain level, mini smooth dogbone (SDB) specimens were cut along RD and TD from the uniform deformation region of the Marciniak samples and reloaded under uniaxial tension. The dimensions of mini smooth dog-bone (SDB) tensile specimens are given in Fig. 1 (d). To prevent the possible edge crack formation at the edge of the cup before the targeted strain of the planar homogenous area has been reached, the stroke of the press was limited to just below the point of edge fracture. The obtained pre-strains were measured by the GOM Aramis DIC system. For each loading condition, 3 parallel tests were conducted to guarantee repeatability.

Artificial microstructure model generation
To simulate the microstructure of the chosen material, a representative volume element (RVE) has been employed with the statistically characterized microstructure features in terms of grain size, shape, and texture. This RVE is created by the software DREAM3D and can be adjusted in size and resolution. In our previous study [8], further optimization and filtering are needed to avoid exaggerated deviation of microstructure characteristics. Therefore, the microstructural representative evaluation criterion (MRAC), which can reduce the deviation to 0 [10]. It is defined in the equation below: Where is the characteristic values, and i represents the type of the microstructure features, and n is the total number of the characteristic values involved. The ref and RVE are the characterized values of the corresponding microstructural feature X from RVE input and output, respectively.
During the RVE generation process, the most representative and effective factors are element size and number, which are corresponding to the final RVE size. The finer element size and a large RVE will have a better representative but they are also will cause a higher computational cost. Considering the balance of RVE representativeness and the computational performance, the final element number is 64,000 elements with 40 elements per edge of the RVE, the element size is 3 μm and the whole RVE size is 120x120x120 μm 3 , respectively. The finally optimized RVE model is shown in Fig. 2, which contains 3165 Al grains. The detailed comparisons between optimal RVE and the reference material microstructure are referred to Lian, et al. [8]

Phenomenological crystal plasticity model
The crystal plasticity model and the governing equations are based on the dislocation slip mechanism theory [12]. In the phenomenological crystal plasticity model, the slip rate ̇ is calculated by the kinetic law of the slip plane system (Eq. 2).
0 is the reference shear rate and the rate sensitivity. α is the resolved shear stress on the slip system and c α the critical resolved shear stress. α is defined by Eq. 3. (4) ℎ is the hardening matrix and is given by Eq. 5.
ℎ 0 , a and c s are the slip hardening parameters. αβ descrives the effect of self-hardening ( = ) and latent hardening ( ≠ ). Finally, the hardening law of the slip system is determined by Eq. 6.
0 is the initial critical shear stress. Considering the quasi-static deformation, the calibrated crystal plasticity parameters of Eq. 6 are 0 , ℎ 0 , c s and . As the strain rate effect is not considered in this study, the strain rate sensitive parameters, 0 and , have been taken from the literature [10].

Dislocation-based crystal plasticity model
In the dislocation-based crystal plasticity model, the dislocation densities ( e , edge dislocation density and d , dipole dislocation density) and their evolution are involved. The Orowan equation is used to define the shear rate of the slip system α.
In Eq. 7, b s , Burger's vector length for slip, represents the distance and direction of a dislocation in the crystal plane. Dislocation glide velocity, 0 , represents the velocity of the dislocation, which is dependent on the applied shear stress and temperature. , the activation energy for dislocation slip, means the minimum amount of energy needed for a dislocation slip to occur. Previously mentioned (9) where, G and ´ are material constants, defined as the shear modulus and interaction matrix between α and α' slip planes, respectively.
The evolution of edge and dipole dislocation densities on α' slip system is defined as, where, the mean free path for slip is defined as and climb dislocation velocity given by, where, climbing activation volume is the fitting parameter to control the volume required by the material in the process of surpassing the energy barrier. The activation energy for climb, the FCC aluminum self-diffusivity coefficients ( , and 0 ) are physical constants that define the motion of a single molecule in the same molecule layer. ̂α , is the maximum slip plane distance required for two dislocations to form a dipole and is the minimum distance required for the elimination of two edge dislocations, which is a fitting parameter. ̂α is defined in the following equation: The dislocation mean free path will describe hardening during the straining process. Considering the slip as the deformation mechanism, it can be described as: From here two new parameters are obtained. d is the average grain size of material, and slip is another fitting parameter.

Parameter Calibration
In both phenomenological and dislocation-based models, Hooke's law is used for explaining the elastic behavior and the same elastic parameters from the literature [13] have been used. In addition, these two models also share some same settings for the calculation of plastic deformation. For example, the solid solution strength, sol or τ0. For the rest of the fitting parameters, Fig. 3 shows a schematic drawing of the calibration process. For both models, the parameters controlling the yield point shall be calibrated first. After the experimental yield strength has been picked in the simulation, the rest hardening parameters can be calibrated according to the flow curve obtained from RD uniaxial tension. Table 1 lists the parameter set used for the phenomenological model.  -In addition to the fitting parameters, there are other microstructural parameters and material constants in the dislocation-based model that need to be fixed first. The microstructural parameters, like grain size (d) and dipole dislocation density ( d ), can be analyzed based on the EBSD measurement or obtained from the literature. The rest of the material constants such as Burger's vector, dislocation glide velocity, activation energies for slip and climb, and self-diffusion coefficient are obtained from the literature. All parameters and constants used for the dislocation-based model are summarized in Table 2.
To calibrate all six parameters in the dislocation-based mode, a parametric study has been carried out to investigate the effects of each parameter. The flow curve from tensile testing along RD is used for the parameter calibration. It is indicated that p influences the yielding point of the material and while the parameter value increases, the yield strength also increases. While the parameter q has opposite effects with the p parameter, that is, with the q value decreasing, the yield strength increases. sol is the third parameter controlling the yield strength, a higher value for solid solution strength is the equivalence of a higher-yielding point. The islip parameter controls the plastic hardening region and leads to the different hardening coefficients. With a lower value of islip, the strength of the material will be higher. The effects of and α are similar and just control the end part of the flow curve. It can be seen that the calibration process of the dislocation-based model is more difficult than that of the phenomenological model in this section. It involves more parameters and a lot of effort. Therefore a quantitative evaluation criterion is necessary to assess the parameter calibration performance.

Prediction on linear loading case
The experimental result in Fig. 4 shows an inhomogeneous plastic deformation, which results in a serrated stress-strain curve. This phenomenon is attributed to the Portevin-Le Chatelier (PLC) effect. This effect, called also dynamic strain aging, is a plastic instability observed in many alloys (such as Al-Mg alloys) when deformed at certain ranges of strain rates and temperatures [17]. The origin of the PLC effect lies in dynamic interactions between mobile solid solution atoms and dislocations.
For analyzing the anisotropic behavior of investigated material, the stress-strain responses of different loading angles have been used for testing the prediction of the two models and comparing their performance (See Fig. 4, as similar results have been obtained in all loading angles, here only 0º/45º/90º results are presented.). It can be seen in the case of the dislocation-based model, early strain rates seem quite good under different loading angles. However, the prediction results of this method are not very good at medium and large strain levels, which may be related to the low quality of parameter calibration in this range. While the phenomenological model shows good adaptability with the loading angle change. It makes a good prediction of large strain values at almost all loading angles. But for the early strain rates, the dislocation-based model can match the experimental data better than the phenomenological model. Particularly for 90º, the stress prediction results at large plastic strain are notably lower than the experimental data. For 45º, both of the phenomenological and dislocationbased models show a good performance for the prediction. (c) 90° Fig. 4 Anisotropic prediction of phenomenological and dislocation-based CP models.

Prediction on non-linear loading case
In this part, the experimental values were compared with the predicted values of the CP model, as shown in Fig. 6. The experimental data come from two parts, one is the values of the linear tensile test mentioned in the previous section, the other is the values obtained by the tensile test of the SDB that were cut from the material after three kinds of pre-strain of the material. There are three kinds for pre-strain the biaxial, plane strain, and uniaxial pre-strain. The simulation results of the CP model also come from two parts, one is the phenomenological model and the other is the dislocation-based model. In addition, all experiments and simulations have been carried out from both RD and TD directions. From the experimental results, it is clear that the stress response has changed after a certain prestrain loading. Particularly for TD , the stress values of the materials after the three kinds of pre-strain are much higher than the experimental values of linear loading at the same strain level. However, the stress values of the uniaxial strain path are only slightly higher than the initial specimen, especially for the RD uniaxial mini-SDB specimen (Fig. 5 c), as the stress states are identical in two stages.
Overall, both models performed slightly better along the RD. This is because the calibration of parameters has been carried out based on the initial RD linear loading result. For the biaxial loading path (Fig. 5 a), the simulation results of the dislocation-based model have better prediction than the phenomenological model at the smaller strain values. However, it shows a large deviation under large strain conditions. On the contrary, phenomenology shows higher quality at larger strain values. For