Crystal Plasticity Modeling for Non-ferrous Metals and its Engineering Applications

Abstract

Crystal plasticity models enable predictions of macroscopic deformation behavior as well as texture evolution of metallic materials based on mesoscopic deformation at the grain level. Owing to recent improvements in predictive accuracy, crystal plasticity models are expected to be used not only for academic purposes but also for industrial applications. There are several possible approaches for utilizing crystal plasticity models in industrial applications, including numerical material testing, in which the material parameters of phenomenological constitutive models are determined; alternative constitutive equations in simulations; and the development of innovative materials with improved formability. In this review paper, recent progress in crystal plasticity modeling, specifically in terms of engineering applications, is discussed. The focus is primarily on hexagonal close-packed (hcp) metals, including magnesium alloy and commercially pure titanium sheets, which exhibit strong anisotropic and asymmetric deformation behavior. On the basis of our recent progresses, the crystal plasticity modeling was first explained, followed by some application examples for a variety of loading conditions, including uniaxial tension and compression, reverse loading, and biaxial tension. The application to face-centered cubic (fcc) and body-centered cubic (bcc) metals and future prospects are also discussed.

1. Introduction

The demand for lightweight materials to reduce CO₂ emission from vehicles is increasing; a variety of metallic sheets are used for structural components of vehicles manufactured by press forming. In press forming, it is vital to set the forming conditions appropriately depending on the deformation characteristics of the sheet, but it is quite difficult to determine them experimentally. Therefore, in recent years, finite-element simulations are indispensable to understanding the press formability and designing process.

In finite-element simulations, elastoplastic deformation of sheet metals is generally described by using phenomenological constitutive models. The main purpose of such models is to mathematically represent the macroscopic deformation behavior of sheet metals obtained experimentally. Therefore, the following two technical problems are often faced. The first problem is that a huge number of material tests are required to utilize constitutive models effectively. Along with the growing need for accurate simulations, advanced constitutive models are often used¹⁾ in which several parameters are involved. In order to identify the parameters required in advanced constitutive models, it is essential to conduct a variety of material tests such as uniaxial tension, compression, reverse loading, and biaxial loading tests before finite-element simulations are conducted.^2,3) These experiments require significant effort and cost, which often hinders effective utilization of advanced constitutive models.

The second problem is that the deformation characteristics of target sheets can hardly be represented in simulations because appropriate constitutive models are still unavailable. For example, sheet metals with hexagonal close-packed (hcp) structure, including Mg alloy and α-Ti, have recently attracted attention as lightweight materials. In hexagonal metals, the deformation characteristics at the grain level, including critical resolved shear stress (CRSS) and work hardening, considerably differ depending on the slip system.^4,5) Moreover, deformation twinning with a polar character is also active.⁶⁾ These strong crystal anisotropies lead to notable asymmetry and anisotropy in the macroscopic deformation behavior.⁷⁾ However, material models suitable for these hexagonal metals have not been established; therefore, forming simulations of hexagonal metals are difficult.

Crystal plasticity models are expected to address the aforementioned technical problems and act as a bridge between mechanical properties of materials and forming simulation. The primary aim of crystal plasticity modeling is to predict macroscopic deformation behavior as well as texture evolution of metallic materials by modeling the activity of slip and twinning at the grain level in the framework of continuum mechanics.⁸⁾ The basis of crystal plasticity models was already established in the early 1980s.^9,10,11) Initially, crystal plasticity models were mainly utilized for fundamental studies, such as strain localization analyses and predictions of texture evolution. Thereafter, due to the improvements in computational power, crystal plasticity modeling, and simulation technique, the number of studies on quantitative prediction of macroscopic deformation behavior, including work hardening, contour of equal plastic work, and Lankford values, has recently increased especially for hcp metals;^{e.g.12,13,14,15,16,17,18,19,20,21,22)} as a result, crystal plasticity models are also drawing attention in industrial applications. Contrary to phenomenological constitutive models, one of the advantages of crystal plasticity models is that they can predict the macroscopic deformation characteristics even if they have not been observed experimentally. In other words, the aforementioned experiments necessary to determine material parameters could be replaced by crystal plasticity simulations, thereby determining the parameters quite efficiently. Furthermore, it is also expected that the use of crystal plasticity models instead of phenomenological models in forming simulations would further improve the predictive accuracy and deliver added value to finite-element simulations.

The present author and his co-workers have been developing crystal plasticity models that can be used to accurately predict plastic deformation behavior of different metallic sheets. Specifically, we focused our attention on metallic sheets with hcp structure, such as Mg alloy and commercially-pure Ti (CP–Ti) sheets. In this review paper, on the basis of our previous achievements,^{13,14,15,16,17,18,19,20,21,22,23,24)} crystal plasticity modeling and some application examples for prediction of deformation behavior under various loading paths are explained, especially focusing on Mg alloy and CP–Ti sheets. Thereafter, recent tendencies and future prospects for practical application of crystal plasticity models are given. The application to fcc and bcc metals is also discussed.

2. Crystal Plasticity Modeling

2.1. Constitutive Equations for Single Crystals

The crystal plasticity model employed in our study is briefly explained in this section. The reader is referred to literatures for detailed formulations.^9,10,11)

The velocity gradient tensor L is assumed to consist of the plastic part L^p and the elastic part L^* in the form   

$$\mathit{L}={\mathit{L}}^{\text{*}}+{\mathit{L}}^{\text{p}}.$$(1)

L^* is associated with rotation and distortion of crystal lattice and rigid body motions. The plastic velocity gradient tensor L^p is related to the activities of slip and twinning systems as follows   

$${\mathit{L}}^{\text{p}}=\sum _{\alpha =1}^N{\dot{\gamma }}^{\alpha }(s^{\alpha }\otimes {\mathit{m}}^{\alpha }),$$(2)

where ${\dot{\gamma }}^{\alpha }$ denotes the slip/shear-strain rate of the α slip/twinning system; s^α and m^α are the unit vectors that denote respectively the slip/twinning direction and the slip/twinning normal plane; and N is the number of slip and twinning systems in a grain.

Assuming that the Cauchy stress σ is governed by the elastic strain rate D^*, the elastic constitutive equations are given as   

$${\stackrel{\circ }{\sigma }}^{\text{*}}=\dot{\sigma }-{\mathit{W}}^{\text{*}}\cdot \sigma +\sigma \cdot {\mathit{W}}^{\text{*}}=\mathit{C}\text{:}{\mathit{D}}^{\text{*}},$$(3)

where C is the elastic constitutive moduli. The ° and · symbols over the elements denote, respectively, the Jaumann rate and material derivative. D and W are the symmetric and antisymmetric parts of L, respectively, and their additive decompositions are given in the following forms:   

$$\begin{array}{l}\mathit{D}={\mathit{D}}^{\text{*}}+{\mathit{D}}^{\text{p}},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ {\mathit{D}}^{\text{p}}=\sum _{\alpha =1}^N{\dot{\gamma }}^{\alpha }{\mathit{p}}^{\alpha },\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\mathit{p}}^{\alpha }=\frac{1}{2}({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha }+{\mathit{m}}^{(\alpha )}\otimes {\mathit{s}}^{(\alpha )})\end{array}$$(4)

  

$$\begin{array}{l}\mathit{W}={\mathit{W}}^{\text{*}}+{\mathit{W}}^{\text{p}},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ {\mathit{W}}^{\text{p}}=\sum _{\alpha =1}^N{\dot{\gamma }}^{\alpha }{\omega }^{\alpha },\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\omega }^{\alpha }=\frac{1}{2}({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha }-{\mathit{m}}^{\alpha }\otimes {\mathit{s}}^{\alpha })\end{array}$$(5)

Substituting Eqs. (4) and (5) into Eq. (3), we obtain the constitutive equations for single crystals:   

$$\begin{array}{l}\stackrel{\circ }{\sigma }=\dot{\sigma }-\mathit{W}\cdot \sigma +\sigma \cdot \mathit{W}\\ \hspace{1em}=\mathit{C}\text{:}\mathit{D}-\sum _{\alpha =1}^N{\dot{\gamma }}^{\alpha }(\mathit{C}:{\mathit{p}}^{\alpha }+{\omega }^{\alpha }\cdot \sigma -\sigma \cdot {\omega }^{\alpha }).\end{array}$$(6)

In the present study, the crystalline slip is assumed to follow Schmid’s law. Because most metallic materials exhibit some degree of rate sensitivity, it would be reasonable to assume that the slip rate ${\dot{\gamma }}^{\alpha }$ is given by a rate-dependent equation in the form   

$$\frac{{\dot{\gamma }}^{\alpha }}{{\dot{\gamma }}_0}={\left|\frac{{\tau }^{\alpha }}{{\tau }_{\mathrm{Y}}^{\alpha }}\right|}^{\frac{\text{1}}{m}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\tau }^{\alpha }),\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\tau }^{\alpha }=\sigma :({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha }),$$(7)

where τ^α is the Schmid’s resolved shear stress, ${\tau }_{\mathrm{Y}}^{\alpha }$ is the slip resistance of the α slip system, ${\dot{\gamma }}_0$ is the reference slip rate, and m is the rate sensitivity exponent. This equation assumes that all slip systems are potentially active during deformation. To represent work hardening of the α slip system, the evolution of ${\tau }_{\mathrm{Y}}^{\alpha }$ is given as   

$${\dot{\tau }}_{\mathrm{Y}}^{\alpha }=\sum _{\beta }{q}_{\alpha \beta }h\left|{\dot{\gamma }}^{\beta }\right|,$$(8)

where q_αβ is the self (α = β) and latent (α ≠ β) hardening moduli and h is the rate of hardening. h is usually given in either phenomenological equations or dislocation-density-based models.^{e.g. 25,26)} In the present study, the following phenomenological laws²⁷⁾ are used   

$$h=h_0,$$(9)

  

$$h=h_0\left(1-\frac{{\tau }_0}{{\tau }_{\infty }}\right)\text{exp}\left(-\frac{h_0\overline{\gamma }}{{\tau }_{\infty }}\right),\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\overline{\gamma }=\sum _{\alpha }\int \left|{\dot{\gamma }}^{\alpha }\right|\mathrm{d}t,$$(10)

where τ₀ is the CRSS, and h₀ and τ_∞ are the material parameters. Modeling of twinning systems is described in the following section.

It should be noted that, in the case of bcc metals, sometimes the so-called non-Schmid effect is considered in crystal plasticity models. This effect includes some different phenomena,^28,29,30,31) such as the case where the slip/twinning activity is represented by not only τ^α but also other stress components, and the case where the slip resistance on the {112} slip plane is different between the twinning and anti-twinning directions. However, because the role of this effect in mesoscopic deformation is not well understood, its effect on macroscopic deformation behavior remains open to discussion.

2.2. Selection of Slip and Twinning Systems

In crystal plasticity modeling, appropriate slip and twinning systems selection depending on the material is significant. In the case of Mg alloys, Graff et al.²⁷⁾ concluded that at least some <c + a> slip should be considered in addition to the basal (0001)<1120> slip (basal slip) as the primary slip system. Agnew et al.⁵⁾ reiterated that some <c + a> slip should be considered to properly predict the evolution of Lankford value in rolled Mg alloy sheets. Following the previous studies, the basal slip, prismatic {1010}<1120> slip (prismatic slip), and pyramidal {1122}<1123> slip (pyramidal <c + a> −1 slip) systems were considered in our study.^{13,14,15,16,17,18)} The selection of slip systems is also similar in recent literature,^{32,33,34,35,36)} although sometimes the pyramidal {1101}<1 120> slip systems are considered in addition to the aforementioned combination.

The number of potentially active slip and twinning systems is much larger for CP–Ti than for Mg alloys. Following previous numerical and experimental studies,^{37,38,39,40,41,42,43)} the prismatic, basal, pyramidal {1011}<1120> (pyramidal <a>), pyramidal <c + a> −1, and pyramidal {1011}<1123> (pyramidal <c + a> −2) slip systems were considered for CP–Ti in our study.¹⁹⁾ Again, the choice of slip systems is similar in recent literatures despite of small variations.^44,45,46,47)

The difference in slip activities depending on the family was described by the hardening laws and material parameters. In our study, the saturation law (Eq. (10)) was assumed for the slip systems. However, for the Mg alloys, the linear hardening law (Eq. (9)) was assumed only for the basal slip systems because it was understood that the activity of these systems was easier compared to other slip systems.

Concerning twinning systems, it is established in Mg alloys that {1012} tensile twinning is active in the early stage of deformation, whereas {1011} compressive twinning is observed in the later stage, resulting in fracture.^48,49) Therefore, both twinning systems tend to be considered in recent literatures.^{32,33,34,35,36)} We focused our attention mainly on the deformation behavior within the uniform elongation range; thus, only {1012} tensile twinning was considered.^{13,14,15,16,17,18)}

In contrast, although activities of several types of twinning systems have been reported in CP–Ti, {1012} tensile twinning and {1122} compressive twinning are more active and could affect plastic deformation behavior even at small strains.^50,51) Therefore, in our study, these two families of twinning systems were considered.¹⁹⁾ This combination is almost equivalent to the used in the recent literature.^44,45,46,47)

The cases for cubic metals are rather simple. The {111}<110> slip systems are considered for face-centered cubic (fcc) metals. For body-centered cubic (bcc) metals, on the other hand, the {110}<111> and {112}<111> slip systems are usually considered, and sometimes also the {123}<111> slip systems. Normally, twinning is not considered for cubic metals.

2.3. Twinning and Detwinning Model

Observing twinning activities at the grain level, it is established that large shear strain and lattice rotation arise simultaneously in twined regions, and these regions expand with twinning activity. As described earlier, twinning is active only under a specific loading direction—i.e., in the case of Mg alloys and CP–Ti, {1012} twinning and {1122} twinning are active only under tension and compression of the c-axis, respectively. However, when {1012} twinning and {1122} twinning are subjected to reverse loading from tension to compression and from compression to tension of the c-axis, respectively, the so-called detwinning is active, i.e., the twined regions start shrinking and large shear strain and lattice rotation arise in the opposite direction of twinning processes. Because our purpose of crystal plasticity simulation was to model macroscopic deformation, this behavior was modelled rather coarsely at the mesoscopic level in our study,^{15,16,17,18,19)} as briefly explained in sequence.

The twinning and detwinning model used in our study was based on the twinning model proposed by Van Houtte.⁵²⁾ One of the features of this model is that the occurrences of shear strain and lattice rotation due to twinning activity are modelled separately, as follows. As in the case of slip, the shear strain rate of the α twinning system, ${\dot{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }\hspace{0.17em}$, is assumed to be characterized by the resolved shear stress τ^α. Because the shear strain due to twinning within each grain may be considered to be mesoscopically similar to that of slip, the shear strain rate ${\dot{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }\hspace{0.17em}$ is given by Eq. (7). Moreover, slight work-hardening is numerically allowed by using the linear hardening law (Eq. (9)) to model the expansion of the twined region with increasing stresses. To represent the polar character of twinning activity, it is further assumed that {1012} twinning and {1122} twinning are active only when τ^α is positive and negative, respectively.

Using the cumulative shear strain ${\overline{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }=\int \left|{\dot{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }\right|\hspace{0.17em}\mathrm{d}t$, the volume fraction of the α twinning system $f^{\alpha }={\overline{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }/{\gamma }_{\mathrm{r}\mathrm{e}\mathrm{f}}$ is calculated during simulations, where γ_ref represents the reference shear strain due to twinning. When $f^{\alpha }=f_{\mathrm{t}\mathrm{h}}^{\alpha }$ is satisfied, where $f_{\mathrm{t}\mathrm{h}}^{\alpha }$ is the threshold value for the volume fraction, the orientation of this grain is rotated by using the rotation tensor in the form   

$${\mathit{R}}^{tw}=2{\mathit{m}}^{\alpha }\otimes {\mathit{m}}^{\alpha }-\mathit{I},$$(11)

where m^α is the unit vector normal to the α twin plane, and I is the identity tensor. The threshold value $f_{\mathrm{t}\mathrm{h}}^{\alpha }$(≤1) is determined randomly for each twinning system prior to simulations. This algorithm indicates that each grain can be configured as unrotated or rotated. It is also assumed that the shear strain rate ${\dot{\gamma }}_{\mathrm{t}\mathrm{w}}^{\alpha }\hspace{0.17em}$ cannot develop after the grain is rotated due to the α twinning system.

For grains that experienced the activity of the α twinning system, the present model assumes that detwinning can be active when the sign of τ^α is reversed. The shear strain due to detwinning, ${\dot{\gamma }}_{\mathrm{d}\mathrm{t}\mathrm{w}}^{\alpha }$, is also described by using Eq. (7). Detwinning can be active until ${\overline{\gamma }}_{\mathrm{d}\mathrm{t}\mathrm{w}}^{\alpha }={\gamma }_{\mathrm{r}\mathrm{e}\mathrm{f}}{f}_{\mathrm{t}\mathrm{h}}^{\alpha }$ is fulfilled, where ${\overline{\gamma }}_{\mathrm{d}\mathrm{t}\mathrm{w}}^{\alpha }=\int \left|{\dot{\gamma }}_{\mathrm{d}\mathrm{t}\mathrm{w}}^{\alpha }\right|\hspace{0.17em}\mathrm{d}t$, suggesting that detwinning can be active until the volume fraction of the α twinning system vanishes. This assumption is introduced based on our experimental observations.^53,54) If the grains already experienced lattice rotation due to twinning before the sign of τ^α was reversed, the orientations of the grains are re-rotated by using the rotation tensor ${({\mathit{R}}^{tw})}^{\mathrm{T}}$ when ${\overline{\gamma }}_{\mathrm{d}\mathrm{t}\mathrm{w}}^{\alpha }={\gamma }_{\mathrm{r}\mathrm{e}\mathrm{f}}{f}_{\mathrm{t}\mathrm{h}}^{\alpha }$ is fulfilled.

In recent years, several twinning models have been proposed because of the increasing attention to hcp metals.^{55,56,57,58,59,60)} Excluding the mentioned model, other well-known twinning models are introduced only briefly. Kalidindi⁵⁵⁾ proposed the following form of the plastic velocity gradient tensor L^p   

$$\begin{array}{l}{\mathit{L}}^{\mathrm{p}}=\left(1-\sum _{\beta }^{N_{O-tw}}{f}^{\beta }\right)\sum _{\alpha }^{N_{O-\mathrm{s}\mathrm{l}}}{\dot{\gamma }}^{\alpha }{({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha })}_{O-\mathrm{s}\mathrm{l}}\\ \hspace{1em}\hspace{0.5em}+\sum _{\alpha }^{N_{O-tw}}{\dot{\gamma }}^{\alpha }{\gamma }_{ref}{({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha })}_{O-tw}\\ \hspace{1em}\hspace{0.5em}+\sum _{\beta }^{N_{O-tw}}{f}^{\beta }\left(\sum _{\alpha }^{N_{tw-\mathrm{s}\mathrm{l}}}{\dot{\gamma }}^{\alpha }{({\mathit{s}}^{\alpha }\otimes {\mathit{m}}^{\alpha })}_{tw-\mathrm{s}\mathrm{l}}\right),\end{array}$$(12)

where the superscripts O-s1, O-tw, and tw-s1 denote the active slip systems in the matrix, active twinning systems in the matrix, and active slip systems in the twined regions, respectively. Unlike the Van Houtte model, one of the advantages of this model is that it can consider both matrix and twined regions within each grain. Because large lattice rotation due to twinning is not necessary in this model, simulations are more stable than in the Van Houtte model. On the other hand, twinning activities lead to increase in the number of slip systems in each grain; thus, the computational cost is sometimes high.

Proust et al.⁵⁶⁾ proposed a detwinning model for the first time based on their so-called composite-grain model. This model considers the compatibility conditions of stresses and strains between the matrix and the twined regions. Moreover, this model also considers a Hall-Petch like effect associated with the occurrence of twined regions.

2.4. Crystal Plasticity Finite-element Method

A variety of numerical techniques have been used to evaluate homogenized deformation behavior of polycrystalline metals, including the Taylor model,⁶¹⁾ self-consistent models,^62,63) crystal plasticity Fast Fourier Transform (FFT),^64,65) and crystal plasticity finite-element methods. In our study, the single crystal-plasticity model was introduced into the static finite-element method with explicit time integration.⁶⁶⁾ One of the advantages of finite-element methods is that inhomogeneous deformation fields at the mesoscopic level can be evaluated satisfying the mechanical compatibility conditions. Furthermore, it is easy to consider grain geometries in simulations. In contrast, mesh distortion often affects simulation results, especially at large strains. Therefore, when evaluating the deformation behavior to large strains, the effect of mesh distortion on the predictive accuracy should be carefully examined. High computational cost is also a notable problem in finite-element methods that should be solved for further applications.

For the stability of the explicit time integration, the rate tangent modulus method proposed by Peirce et al.⁶⁷⁾ was adopted. In addition, the size of the increment of each time step was adjusted by using the generalized r_min-strategy to prevent an excessive increase in nonequilibrium between external and internal forces. In fact, the generalized r_min-strategy was quite useful to prevent numerical instability due to lattice rotation of twinning.

A cubic model was utilized as a representative volume element that was uniformly divided in each direction using 8-node solid elements with selective reduced integration. Assuming that an element corresponds to a grain, the same initial orientation was assigned to all 8 integration points in an element. It should be noted that crystal orientations after plastic deformation may vary depending on the integration point even within an element because the crystal rotation was calculated independently for each integration point. Initial crystal orientations were prepared from results of electron backscattered diffraction (EBSD) measurements of the sample. More specifically, initial orientations were randomly extracted from the measured data and assigned to each integration point. To accurately represent macroscopic deformation behavior of polycrystalline metals using the crystal plasticity model, it is important to consider a sufficient number of initial orientations, i.e., the number of finite elements, in simulations. Our previous study confirmed that for Mg alloy¹³⁾ and CP–Ti¹⁹⁾ sheets sufficiently homogenized results could be obtained by dividing the cubic model into 7 to 10 elements in each direction, which corresponds to the number of initial orientations of 343 to 1000. This result is consistent with a statistical analysis conducted by Tadano.⁶⁸⁾

2.5. Identification of Material Parameters

It is important to appropriately determine the material parameters used in the crystal plasticity model, such as the parameters in Eqs. (7), (8), (9), (10) and the twinning model, to conduct reliable simulations. This problem is particularly notable for hcp metals because the activity of slip and twinning systems differs significantly depending on the deformation mechanism. In addition, it is also necessary to refer to a sufficient number of experimental data to systematically determine the parameters. For instance, in the case of Mg alloys, it was shown that several sets of material parameters could be determined to reproduce a stress-strain curve under uniaxial tension, but the evolution of Lankford value was completely different depending on the set, suggesting that a single stress-strain curve is not enough to determine the parameters.⁶⁹⁾

To address this problem, procedures to systematically determine the material parameters for rolled Mg alloy¹⁸⁾ and CP–Ti sheets¹⁹⁾ were proposed in our study. Utilizing the strong basal textures feature, as shown in Fig. 1, the proposed procedures enabled determining the parameters considering the role of each deformation mechanism by referring to experimental results of stress-strain curves under tension, compression, and reverse loading, and evolution of Lankford value along different directions.

Fig. 1.

Initial pole figures measured using EBSD. (a) Mg alloy sheet¹⁸⁾ and (b) CP–Ti sheet.⁵¹⁾ (Online version in color.)

For the twinning model, γ_ref and $f_{\mathrm{t}\mathrm{h}}^{\alpha }$ are also the fitting parameters. These parameters were determined with reference to texture evolution in addition to stress-strain curves. Rigorously speaking, γ_ref would correspond to the microscopic shear strain determined geometrically by the twinning activity. In contrast, the present twinning model focused rather on coarse-grained mesoscopic deformation behavior, as described earlier. Therefore, in our study, γ_ref was allowed to deviate somewhat from the geometrically determined shear strain.

One of the most difficult aspects of crystal plasticity modeling is determining the parameters of the latent hardening moduli q_αβ in Eq. (8). Because it is almost impossible to determine the parameters only by means of experiments, researchers utilized numerical simulations, including discrete dislocation dynamics, for fcc^70,71) and bcc metals.^72,73) However, these results are not conclusive and, moreover, the parameters for hcp metals are still hardly discussed. In our study, the parameters reported by Graff et al.²⁷⁾ were used for Mg alloy sheets, whereas for CP–Ti sheets the parameters were determined by trial and error¹⁹⁾ to fit the abovementioned macroscopic deformation behavior.

3. Prediction of Deformation Behavior by Crystal Plasticity Model

3.1. Application to Mg Alloy Sheets

3.1.1. Uniaxial Loadings¹⁸⁾

This chapter presents some simulation examples of deformation behavior of hcp metals by using the crystal plasticity model described earlier. In this section, experimental and simulation results of Mg alloy sheets are shown. In the experiments, rolled AZ31 Mg alloy sheets with a thickness of 1.0 mm were used unless otherwise noted. In the simulations, the material parameters determined using the procedure explained in section 2.5 were utilized. As an example, Fig. 2 shows the result of parameter identification, i.e., the results of the stress-strain curves under tension, compression, and reverse loading from compression to tension, and the evolution of Lankford value along the rolling direction (RD). The determined parameters are shown in Table 1. The simulation results accurately reproduced the characteristic deformation behavior of the Mg alloy sheet, including the tension-compression asymmetry and the sigmoidal curve after stress reversal, which shows that the parameters were properly determined for these conditions.

Fig. 2.

Deformation behavior of Mg alloy sheet under uniaxial loadings.¹⁸⁾ (a) Stress-strain curves and (b) evolution of Lankford value during tension. (Online version in color.)

Table 1. Material parameters determined for the Mg alloy sheet (MPa).¹⁸⁾

		Basal	Pris	Pyram <c + a>	Twin	Detwin
τ0	/MPa	25	95	100	40	35
τ∞	/MPa	–	230	190	–	–
h0	/MPa	10	1200	1200	5	5

The mechanisms that yield the pronounced tension-compression asymmetry and the sigmoidal curve can be explained using the evolution of relative activities. The relative activity is the relative contribution of each family of slip/twinning systems to the plastic deformation. The relative activity of the family of the slip/twinning system i, r_i, is given in the form   

$$r_i=\frac{{{\displaystyle \sum }}_n^{n_s}{{\displaystyle \sum }}_k\left|\mathrm{\Delta }{\gamma }^{(n,k)}\right|}{{{\displaystyle \sum }}_n^{n_s}{{\displaystyle \sum }}_j\left|\mathrm{\Delta }{\gamma }^{(n,j)}\right|},$$(13)

where k is the number of slip and twinning systems of the family i, j is the number of all slip and twinning systems, and n_s is the number of grains. Figure 3 shows the evolution of relative activities during tension and reverse loading depicted in Fig. 2. Basal and prismatic slips were dominant during tension (Fig. 3(a)), while basal slip and twinning were active during compression (Fig. 3(b)). The large activity of twinning occurred during compression because the strong basal texture was developed in the initial sheet; thus, a lot of grains were subjected to tensile stresses in the c-axis direction during compression. This result suggests that the stress level was much lower under compression than under tension because the CRSS of twinning is lower than that of prismatic slip (Table 1).

Fig. 3.

Evolution of relative activities.¹⁸⁾ (a) Tension and (b) reverse loading from compression to tension. (Online version in color.)

After the stress reversal in Fig. 3(b), detwinning was dominant from the beginning of reverse loading, but it rapidly decreased at a cumulative strain of approximately 0.07 and, alternatively, prismatic slip became active. Interestingly, the subsequent activities were similar to those of tension (Fig. 3(a)). This is due to the fact that detwinning completed at a cumulative strain of approximately 0.07 and the deformation mode shifted to that of monotonic tension. This result indicates that the sigmoidal curve was exhibited, i.e., the stress level increased rapidly in the later stage, because the CRSS of prismatic slip was much larger than that of detwinning. It is emphasized that, using the crystal plasticity model, these complex work-hardening behaviors could be simply reproduced as a result of mesoscopic deformation behavior at the grain level.

3.1.2. Two-step Loadings¹⁸⁾

Since sheets often undergo complex strain-path changes during sheet forming processes, material models need to reproduce deformation behaviors with strain-path changes. To verify the availability of the present crystal plasticity model, the following two-step loading tests were conducted.

Figure 4 shows the schematic diagram of the two-step loading tests. In these tests, 6% tensile or compressive strain was first applied to a large sample prepared along the RD. After small samples were prepared from the pre-strained large sample along an angle θ, 10% tensile strain was applied to the small samples as the second step. θ was set to either 0, 30, 45, 60, and 90°. The purpose of this test was to examine how the work-hardening behavior and the texture evolution differ depending on the angle. Figure 5 shows the stress-strain curves during the second loading. When pre-tension was applied, the stress level was the largest for θ = 0° and tended to decrease as the angle θ increased. The work-hardening rate abruptly changed after yielding for θ = 0°, while it gradually changed for the other conditions. A certain degree of in-plane anisotropy appeared, but it was not significant. On the contrary, a remarkable in-plane anisotropy appeared for the pre-compression condition. An apparent sigmoidal curve was exhibited for θ = 0° because this condition corresponded to the simple reverse loading from compression to tension, shown in Fig. 2(a). On the other hand, the sigmoidal tendency became less pronounced with increasing θ and was not visible for θ = 60° and 90°. The aforementioned work-hardening behavior was well reproduced in the simulation results, verifying the model.

Fig. 4.

Schematic diagram of the two-step loading test. (Online version in color.)

Fig. 5.

Stress-strain curves of Mg alloy sheet during the second loading.¹⁸⁾ (a) Pre-stretched and (b) pre-compressed sheet. (Online version in color.)

The in-plane anisotropy of work-hardening behavior observed in the pre-compressed sheet can be explained as follows. Figure 6 shows the evolution of relative activities of each slip and twinning systems during the second loading for θ = 0°, 45° and 90°. It should be noted that only the activities in the matrix region are shown. As described earlier (Fig. 3(b)), the condition of θ = 0° showed a rapid transition from detwinning-dominated to slip-dominated deformation. The detwinning activity became less pronounced with increasing θ. For θ = 90°, the detwinning activity was quite small from the initial stage of the second loading, whereas the activities of basal and prismatic slips were dominant; thus, rapid transition of deformation mode did not arise as happened with θ = 0°. These results suggest that the sigmoidal tendency became less pronounced with increasing θ because of the decrease of detwinning activity.

Fig. 6.

Relative activities in matrix region during the second loading for θ = (a) 0°, (b) 45° and (c) 90°.¹⁸⁾ (Online version in color.)

It is well known that the twinning activity yields rapid texture evolution. Figure 7 shows the (0001) pole figures before and after the second loading for the pre-compression condition. Strong peaks appeared in the vicinity of the RD before the second loading. This is because of the large activity of twinning during the first loading (compression), as shown in Fig. 3(b). After the second loading, the peaks in the RD disappeared for θ = 0° and resulted in the texture quite similar to the initial one (Fig. 1(a)). This notable change occurred because of detwinning activity. On the other hand, for θ = 45° and 90°, the peaks remained and their intensities tended to become larger as θ increased because of the decrease of detwinning activity. The simulation results were in good agreements with the experimental results. Moreover, these results were consistent with the detwinning activities observed in Fig. 6, which supported the abovementioned deformation mechanism.

Fig. 7.

(0001) pole figures measured before and after the second loading in the pre-compressed Mg alloy sheet.¹⁸⁾ Results (a) before the second loading and after the second loading for (b) θ = 0°, (c) θ = 45°, and (d) θ = 90°. (Online version in color.)

3.1.3. Unloading¹³⁾

Figure 8 shows the stress-strain curves during cyclic tension-unloading. It is well known that Mg alloy sheets show a remarkable nonlinear behavior during unloading, as shown in Fig. 8. Because the unloading behavior is considered important in springback, which is one of the major defects in sheet metal forming, it is crucial to properly model the unloading behavior in simulations. In classical elastoplasticity theories, the unloading process is often modeled as elastic deformation, i.e., the stress-strain relationship during unloading is assumed to be linear with a gradient of Young’s modulus. In other words, if the nonlinear behavior during unloading needs to be modeled using classical elastoplasticity theories, it is essential to assume in advance that nonlinear behavior occurs. In contrast, the nonlinear behavior can be well reproduced without giving any special assumptions if the crystal plasticity model is used, as shown in Fig. 8. Although detailed descriptions are not provided here, the simulation result showed that the nonlinear behavior during unloading could be explained by the following mechanism. Prismatic slip was active during tension (Fig. 3(a)), suggesting that the stress level during tension was governed primarily by the prismatic slip resistance. When the sheet was macroscopically unloaded thereafter, basal slip was active. This is because of the fact that basal slip had a remarkably low slip resistance compared to prismatic slip (Table 1); thus, resolved shear stresses of basal slip could be sufficiently large to activate even during macroscopic unloading. Eventually, the activity of basal slip yielded nonlinear behavior during unloading. This result indicates that the nonlinear deformation behavior during unloading of the Mg alloy sheet is due to plastic deformation. It should be noted that the nonlinearity is further pronounced for the compression-unloading condition because of the pronounced activity of detwinning during unloading.^16,74)

Fig. 8.

Stress-strain curves of Mg alloy sheet under cyclic tension-unloading.¹³⁾ (Online version in color.)

3.1.4. Application to Cast Sheet¹⁷⁾

One of the major advantages of crystal plasticity models is that the effect of the initial crystallographic orientations on the macroscopic mechanical behavior can be examined numerically. In this section, an application example of the present crystal plasticity model to a cast AZ31 Mg alloy sheet with random crystallographic orientations is shown. In the experiment, a sheet with a thickness of 1.0 mm was prepared from an ingot of AZ31Mg alloy.

The simulation model was prepared as follows. Initial orientations were randomly assigned to each element. Because the parameter identification procedure explained in Section 2.5 is applicable only for rolled sheets with strong basal texture, it was difficult to determine precisely the parameters of the cast sheet. Therefore, we first assumed that the parameters identified for a rolled sheet could be applied to the cast sheet. Moreover, assuming that the Hall-Petch law for Mg alloy sheets⁷⁵⁾ was also applicable for CRSS, τ₀ was adjusted to account for the difference in grain size between the rolled and cast sheets.

Figure 9(a) shows the stress-strain curves during tension and compression. It should be noted that both axes are in absolute values to directly compare the results of tension and compression. The stress level was higher during tension than during compression, exhibiting the tension-compression asymmetry. The asymmetry was predicted qualitatively well in the simulation results. Figure 9(b) shows the stress-strain curves during reverse loadings from tension to compression and from compression to tension. A small sigmoidal tendency occurred after stress reversal as in the case of the rolled sheet (Fig. 2), when the loading direction was reversed from compression to tension, whereas such a trend was not observed when reversed from tension to compression. The strain-path dependency was predicted qualitatively well in the simulations. These results display that although the parameters were determined rather roughly, the crystal plasticity model could capture the effect of initial orientations on the work-hardening behavior.

Fig. 9.

Stress-strain curves of a cast Mg alloy sheet.¹⁷⁾ (a) Tension and compression and (b) reverse loading. Solid and broken lines are respectively experimental and simulation results. (Online version in color.)

It is interesting that both tension-compression asymmetry and strain-path dependency also occurred in the cast sheet with random orientations. Although detailed mechanism is not explained here, the simulation results suggested that the tension-compression asymmetry and the strain-path dependency arose due to the difference in twinning activities between tension and compression, as in the case of rolled sheets. This result indicated that the tension-compression asymmetry and the strain-path dependency are inevitable in Mg alloy sheets regardless of the initial orientations.

3.2. Application to a CP–Ti Sheet¹⁹⁾

3.2.1. Uniaxial Loadings

In sequence, application examples to CP–Ti sheets are described. In the experiment, a rolled CP–Ti grade 1 sheet with a thickness of 1.0 mm was used. Figures 10(a) and 10(b) show the stress-strain curves under uniaxial loadings and the evolution of Lankford value under tension, respectively. The two axes of Fig. 10(a) are in absolute values. The initial yield stress and subsequent work hardening between tension and compression and between the RD and the transverse direction (TD) considerably differed. The Lankford value in the RD was approximately 1.5 and remained almost unchanged during deformation, while in the TD, it was remarkably high at the beginning of deformation and then decreased rapidly. In the case of the CP–Ti sheet, these experimental results were used to determine the parameters of the crystal plasticity model shown in Table 2. As a result of the parameter identification, the simulation results accurately captured the anisotropic deformation behavior.

Fig. 10.

Deformation behavior under uniaxial loadings of the CP–Ti sheet.¹⁹⁾ (a) Absolute true stress-absolute true strain curves. Solid and broken lines are simulation and experimental results, respectively. (b) Evolution of Lankford value under uniaxial tension. Closed and open circles are experimental and simulation results, respectively. (Online version in color.)

Table 2. Material parameters determined for the Cp–Ti sheet (MPa).¹⁹⁾

	Basal	Pris	Pyram <a>	Pyram <a + c>−1	Pyram<a + c>−2	Ten. Twin	Comp. Twin
τ0	133	62	81	145	145	90	140
τ∞	180	160	210	270	270
h0	1950	1050	580	2050	2050	350	350

Figure 11 shows the (0001) pole figures measured after the uniaxial loadings. When the sheet was loaded in the RD, peaks occurred in the vicinity of the RD for both tension and compression, but the intensity was stronger in compression. In contrast, when loaded in the TD, weak peaks occurred in the RD in tension, while the pole figure hardly changed from the initial one (Fig. 1(b)) in compression. These results depicted that the texture evolution strongly depended on the loading direction. These tendencies are well predicted in the simulation results. Although detailed results are not provided here, the simulation results suggested that the difference in texture evolution depending on the loading condition would be because of the differences in the twinning activities and in the types of active twinning systems.

Fig. 11.

(0001) pole figures measured on samples subjected to uniaxial loadings to an absolute strain of 10% of the CP–Ti sheet.¹⁹⁾ (a) Tension along RD, (b) compression along RD, (c) tension along TD, and (d) compression along TD. (Online version in color.)

Figure 12 shows the stress-strain curves under reverse loading in the RD. When the sheet was subjected to reverse loading from compression to tension, a small sigmoidal tendency occurred after reversal. In contrast, in the case of reverse loading from tension to compression, such tendency did not occur, showing that the CP–Ti sheet also exhibited a strain-path dependency as in the case of Mg alloy sheets. These tendencies were well predicted in the simulation results. Although detailed results are omitted, the simulation results showed that the sigmoidal tendency observed under reverse loading from compression to tension was due to detwinning activity, suggesting that the mechanism that yielded the strain-path dependency was quite similar between Mg alloy and CP–Ti sheets.

Fig. 12.

True stress-true strain curves under reverse loading along RD of the CP–Ti sheet.¹⁹⁾ Solid and broken lines are simulation and experimental results, respectively. Loading directions were inverted (a) from compression to tension and (b) from tension to compression. (Online version in color.)

3.2.2. Contours of Equal Plastic Work

Because sheet metals are often subjected to multiaxial loadings in practical forming processes, it is also important to understand and to model deformation behavior under multiaxial stress conditions. Figure 13 shows the simulation results of the contour of equal plastic work of the CP–Ti sheet. Both axes are non-dimensionalized using the stress under uniaxial tension in the RD, σ₀. Noteworthy asymmetry occurred between the RD and the TD in the initial contour. Specifically, the stresses were globally higher in the TD than in the RD. However, the non-dimensionalized contour tended to shrink with the progress of plastic deformation. Interestingly, the degree of shrinkage was more noticeable in the TD than in the RD. Thus, the contour resulted in a nearly symmetrical shape at a uniaxial tensile strain of approximately 0.085. The initial asymmetry and the anisotropic hardening were in qualitatively good agreements with the experimental results reported by Ishiki et al.,⁷⁶⁾ displaying that the crystal plasticity model is also useful to predict deformation behavior under multiaxial stress conditions.

Fig. 13.

Simulation results of contours of equal plastic work of the CP–Ti sheet.¹⁹⁾ ${\epsilon }_0^p$ denotes plastic strain under uniaxial tension along RD.

From the simulation results, the mechanism that yielded the strong initial asymmetry can be explained as follows. The most active slip systems were different depending on the stress ratio because of the strong rolling texture. Specifically, as illustrated in Fig. 13, the activities of prismatic slip, pyramidal<a> slip, and either basal slip or pyramidal<a> slip were dominant for the conditions σ_RD>σ_TD, σ_{RD =}σ_TD, and σ_RD<σ_TD, respectively. The CRSSs of basal slip and pyramidal<a> slip were larger than the CRSS of prismatic slip; thus, the stresses were higher in the TD than in the RD. Similarly, the contour was more expanded in the σ_RD<σ_TD region than in the equibiaxial condition because the CRSS of basal slip was larger than that of pyramidal <a> slip. These results suggested that the strong initial asymmetry occurred because of the strong rolling texture and the notable difference in the CRSSs between the slip systems.

It should be noted that the contour of the plastic work in Mg alloy sheets could also be well predicted by using the same crystal plasticity model,¹⁴⁾ indicating that deformation behavior of different materials can be predicted by using a same simulation program; it can be achieved if the active slip and twinning systems are selected according to the material and the material parameters are properly determined. This advantage shows the applicability of crystal plasticity models for a wide range of sheet metals with less effort.

3.3. Applications to FCC and BCC Metals

As described in the previous sections, anisotropic and asymmetric deformation behavior of hcp metals could be approximately reproduced by using the aforementioned crystal plasticity model. One of the reasons of this achievement is that the noticeable crystal anisotropy of hcp structure, such as the large differences in CRSS and twinning activity, is important in the macroscopic anisotropy and asymmetry, and moreover, the resulting macroscopic anisotropy and asymmetry are also quite pronounced. On the contrary, the crystal anisotropy is much less pronounced in fcc and bcc metals, and therefore, the resultant macroscopic anisotropy is also much less pronounced than that of hcp metals; thus, more precise modeling, including parameter identification, is required for fcc and bcc metals. Some recent studies on crystal plasticity analyses of Al alloy sheets are introduced in sequence.

Yoshida et al.⁷⁷⁾ studied the predictive accuracy for a 3000-series Al alloy sheet and reported that anisotropic work-hardening behavior could be qualitatively predicted if a dislocation-density based hardening model was utilized. However, the degree of anisotropy was much less pronounced in the simulation. Zhang et al.⁷⁸⁾ used several crystal plasticity models to reproduce the in-plane anisotropy of AA3103-H18 and AA3103-O sheets and found that the Alamel-type models and crystal plasticity FEM yielded better predictive accuracy than the other models used. Brahme et al.⁷⁹⁾ predicted the work-hardening behavior of an AA5754 Al sheet using different crystal-plasticity models and found that the predicted work-hardening for large strains was different depending on the model. This result suggested that the parameters determined by using a stress-strain curve until ultimate strength did not always accurately reproduce the work-hardening behavior under various loading paths. Hu et al.⁸⁰⁾ used different crystal plasticity models to reproduce stress–strain curves and texture evolution of an AA5754 Al sheet. They reported that the predictive accuracy of the mechanical responses was dependent on the loading path.

To the author’s understanding, the predictive accuracy for bcc metals is worse than that of fcc metals and needs to be notably improved. Jeong et al.⁸¹⁾ used a viscoplastic self-consistent crystal plasticity model to predict the anisotropic hardening of plastic contours in an interstitial-free steel sheet. However, the simulation results were not in good agreement with the experimental results. Eyckens et al.,⁸²⁾ Coppieters et al.⁸³⁾ and Hama et al.⁸⁴⁾ reported that for steel sheets the anisotropic hardening of the plastic contour, especially for the equibiaxial tension condition, could not be reproduced regardless of the crystal plasticity model they utilized. One of the reasons for this difficulty in the simulation of bcc metals is that the deformation characteristics at the crystalline level, including the difference in activities of the {110} and {112} slip systems,^31,84) the effect of non-Schmid law,^28,29,30) and the latent-hardening/interaction matrix,^70,71) are not understood and modelled well, as described earlier. For instance, a same CRSS value is usually assigned to both {110} and {112} slip systems in crystal plasticity simulations, but its validity is not yet confirmed. Hama et al.⁸⁴⁾ reported that the anisotropic hardening of plastic contours could be qualitatively reproduced when different CRSS values were assigned to the {110} and {112} slip systems, as shown in Fig. 14. In addition, they also showed that the latent hardening parameters remarkably affected the anisotropic hardening. To verify these results and to further improve the predictive accuracy, it is indispensable to understand these microscopic deformation characteristics, especially at the grain level.

Fig. 14.

Contours of equal plastic work of a steel sheet.⁸⁴⁾ (a) Experimental results and (b) simulation results. In the simulation, different CRSS values were assigned to {110} and {112} slip systems. (Online version in color.)

4. Concluding Remarks: Recent Tendencies and Future Prospects

As concluding remarks, recent tendencies and future prospects of crystal plasticity models are discussed with a focus on practical applications, particularly for sheet metal forming processes.

The author considers that one of the most realistic industrial applications at the present stage is numerical material testing in which the material parameters of phenomenological constitutive models are determined. As mentioned in the introduction, several experiments are typically required to select an appropriate constitutive model and determine its material parameters. If crystal plasticity analyses could replace the experiments, time, cost, and effort required to conduct a variety of material tests could be considerably reduced; thus, it is expected that the modeling accuracy would further improve. Some attempts have already been reported. Inal et al.,⁸⁵⁾ Zhang et al.,⁸⁶⁾ and Hashimoto et al.⁸⁷⁾ determined the parameters of phenomenological constitutive models for Al alloy sheets by using crystal plasticity models. Yamanaka et al.⁸⁸⁾ conducted a simulation of hydraulic bulging of a 5000 series Al alloy sheet using the material parameters determined by using a crystal plasticity model and showed that the simulation results were as accurate as those obtained using the parameters determined experimentally. A similar study on a cold-rolled steel sheet was also conducted by Coppieters et al.,⁸³⁾ verifying their ALAMEL model. Roters et al.⁸⁹⁾ developed the open source software “DAMASK,” which contained a variety of crystal plasticity models to describe plasticity of metals and constitutive models. Han et al.⁹⁰⁾ evaluated the material parameters of phenomenological models and their evolution with anisotropic hardening of a 2090-T3 Al alloy sheet by using DAMASK and simulated a deep drawing process using the determined parameters. They exhibited that the prediction of earing could be improved when the evolution of anisotropy was considered. Zhang et al.⁹¹⁾ also employed DAMASK to identify phenomenological models for AA3104 Al alloy sheets. It is expected that this approach would be a solution to the first technical problem of present finite-element simulations described in the introduction. Moreover, various advanced material tests can be performed quite easily by using crystal plasticity models, leading to a deeper understanding of plastic deformation behavior of sheet metals and eventually to an improvement of phenomenological constitutive modeling. This approach is especially useful for materials for which phenomenological constitutive models have already been well established.

Utilizing crystal plasticity models directly in forming simulations as alternative constitutive equations would also be an effective approach, although significant effort is still required for practical use from the viewpoint of computational costs. For instance, this approach would be quite effective for hcp metals for which phenomenological model has not been established yet, but the complex deformation behavior can be accurately predicted by crystal plasticity models. Hama et al.⁹²⁾ recently performed FEM simulations of cylindrical cup deep drawing processes of a CP–Ti sheet by using the aforementioned crystal plasticity model and showed that the simulation results were in good agreements with experimental results in terms of earing formation, thickness strain distribution, and texture evolution. Because the texture evolution could be predicted simultaneously, its effect on the evolution of mechanical anisotropy could also be considered in the simulations; thus, the improvement of the predictive accuracy can be expected. Moreover, Hama et al.⁹³⁾ conducted FEM simulations of time-dependent springback of a CP–Ti sheet by using the crystal plasticity model and showed that it occurred primarily because of slip and twinning activities driven by the residual stresses. As understood from these examples, one of the advantages of this approach is that the deformation mechanism at the crystalline level during forming processes can be examined numerically.

There are also several application examples for cubic metals, including stamping simulation of Al alloy and steel sheets by Nakamachi et al.,⁹⁴⁾ simulations of deep drawing processes of a 17% Cr stainless steel sheet by Tikhovskiy et al.,⁹⁵⁾ simulations of deep drawing processes of an Al alloy sheet by Van Houtte et al.,⁹⁶⁾ Tikhovskiy et al.,⁹⁷⁾ and Inal et al.,⁹⁸⁾ and simulations of ECAP processes of pure copper by Zecevic et al.⁹⁹⁾ It is expected that further improvement of this approach will lead to a solution to the second technical problem of the present finite-element simulations described previously.

In addition to the aforementioned approaches, it may also be useful to utilize crystal plasticity models for the development of innovative materials with improved formability. Using the advantage of crystal plasticity models that macroscopic and mesoscopic deformation can be correlated seamlessly, textures with rich formability can be numerically explored. Yoshida et al.¹⁰⁰⁾ investigated the forming limit of Al alloy sheets by using a crystal plasticity model and explored textures with excellent forming limit. Nakamachi et al.¹⁰¹⁾ studied textures with rich bendability in Al alloy sheets. Ikawa et al.¹⁰²⁾ and Muhammad et al.¹⁰³⁾ also studied the bendability of Al alloy sheets using crystal plasticity models and discussed the effects of microstructure and enhancement of the bendability from the microstructural viewpoints. These examples show that crystal plasticity models are beneficial not only for forming processes of metallic materials but also for exploring advanced materials, if the simulation results could be effectively fed back to the development of materials.

However, further efforts are still essential to solve several remaining technical issues and eventually to widely utilize crystal plasticity models to practical applications. One of the remarkable challenges is to establish procedures of parameter identification for a variety of materials. There were some reports that studied the identification of CRSS or work-hardening parameters in addition to those proposed in our study for rolled Mg alloy¹⁸⁾ and CP–Ti sheets.¹⁹⁾ Steglich et al.¹⁰⁴⁾ determined the work-hardening parameters for Mg alloy sheets by referring to anisotropic hardening in contours of plastic work. Recently, Baudoin et al.²⁴⁾ studied appropriate ratios in CRSS values of CP–Ti by referring not only to macroscopic stress-strain curves but also to in-plane strain distributions in an oligocrystal. For bcc metals, the difference of the parameters between the {110} and {112} slip systems were studied both experimentally and numerically,^31,84) as described earlier.

Furthermore, determining latent-hardening/interaction moduli parameters is more problematic. Kubin et al.⁷⁰⁾ and Madec et al.^71,72,73) used discrete dislocation dynamics to determine the parameters for fcc and bcc metals. Graff et al.²⁷⁾ evaluated the parameters by trial and error for pure Mg. Hiura et al.¹⁰⁵⁾ experimentally evaluated latent hardening under basal self and coplanar dislocation interactions in pure Mg and reported that hardening of the coplanar slip systems would increase proportionally to the dislocation density accumulated in the other coplanar slip system. Zecevic and Knezevic¹⁰⁶⁾ extended their dislocation-based crystal plasticity model to incorporate the latent-hardening effects for prediction of anisotropic deformation behavior of AA6022-T4 sheets. However, such results are not conclusive and are still open to discussion. Because the parameters of latent-hardening/interaction moduli would be important especially in the prediction of anisotropic hardening under biaxial loadings,⁸⁴⁾ further studies are essential from academic as well as practical viewpoints.

Computational cost is also an important issue for practical applications. Specifically, finite-element methods are quite time consuming although heterogeneous strain fields can be evaluated reasonably. For this problem, effective calculation algorithms^{e.g., 107)} and effective homogenization techniques have been proposed, including the self-consistent methods^62,63) and the FFT method.^64,65) It is expected that these studies can accelerate the development of practical applications for crystal plasticity models.

Acknowledgements

The author appreciates the help of the graduate students of Kyoto University with conducting our studies reported in the paper. This project was partially supported by JSPS KAKENHI grant numbers 26289271, 17H03428 and 17K06858, and the Amada Foundation grant numbers AF-2010021 and AF-2015020.

References

• 1)   F.  Barlat and  T.  Kuwabara: J. Jpn. Soc. Technol. Plast., 57 (2016), 230.
• 2)   T.  Kuwabara,  A.  Ishiwatari,  M.  Kuroda,  T.  Koizumi,  T.  Wu,  H.  Takizawa,  T.  Hama,  H.  Hamasaki,  K.  Yoshida,  J.  Yanagimoto and  A.  Yamanaka: Bull. Iron Steel Inst. Jpn., 24 (2019), 586 (in Japanese).
• 3)   T.  Kuwabara: Int. J. Plast., 23 (2007), 385.
• 4)   H.  Yoshinaga and  R.  Horiuchi: Trans. Jpn. Inst. Met., 5 (1964), 14.
• 5)   S. R.  Agnew and  O.  Duygulu: Int. J. Plast, 21 (2005), 1161.
• 6)   J. W.  Christian and  S.  Mahajan: Prog. Mater. Sci., 39 (1995), 1.
• 7)   X. Y.  Lou,  M.  Li,  R. K.  Boger,  S. R.  Agnew and  R. H.  Wagoner: Int. J. Plast., 23 (2007), 44.
• 8)   M.  Kuroda and  K.  Shizawa: J. Jpn. Soc. Technol. Plast., 43 (2002), 299 (in Japanese).
• 9)  D. Peirce R. J. Asaro and A. Needleman: Acta Metall., 31 (1983), 1951.
• 10)   R. J.  Asaro and  A.  Needleman: Acta Metall., 33 (1985), 923.
• 11)   R. J.  Asaro: Adv. Appl. Mech., 23 (1983), 1.
• 12)   A.  Staroselsky and  L.  Anand: Int. J. Plast., 19 (2003), 1843.
• 13)   T.  Hama and  H.  Takuda: Int. J. Plast., 27 (2011), 1072.
• 14)   T.  Hama and  H.  Takuda: Comput. Mater. Sci., 51 (2012), 156.
• 15)   T.  Hama and  H.  Takuda: Steel Res. Int., Special Edition, (2012), 1115.
• 16)   T.  Hama,  N.  Kitamura and  H.  Takuda: Mater. Sci. Eng. A, 583 (2013), 232.
• 17)   T.  Hama,  T.  Mayama and  H.  Takuda: Rom. J. Tech. Sci., Appl. Mech., 60 (2015), 24.
• 18)   T.  Hama,  Y.  Tanaka,  M.  Uratani and  H.  Takuda: Int. J. Plast., 82 (2016), 283.
• 19)   T.  Hama,  A.  Kobuki and  H.  Takuda: Int. J. Plast., 91 (2017), 77.
• 20)   T.  Hama,  R.  Matsudai,  Y.  Kuchinomachi,  H.  Fujimoto and  H.  Takuda: ISIJ Int., 55 (2015), 1067.
• 21)   T.  Hama,  K.  Kojima,  M.  Kubo,  H.  Fujimoto and  H.  Takuda: ISIJ Int., 57 (2017), 866.
• 22)   T.  Hama,  T.  Suzuki,  S.  Hatakeyama,  H.  Fujimoto and  H.  Takuda: Mater. Sci. Eng. A, 725 (2018), 8.
• 23)   M.  Kubo,  T.  Hama,  Y.  Tsunemi,  Y.  Nakazawa and  H.  Takuda: ISIJ Int., 58 (2018), 704.
• 24)   P.  Baudoin,  T.  Hama and  H.  Takuda: Int. J. Plast., 115 (2019), 111.
• 25)   C.  Teodosiu,  J. L.  Raphanel and  L.  Tabourot: Proc. Int. Seminar MECAMAT’91, A.A. BALKEMA, Rotterdam, (1993), 153.
• 26)   B.  Peeters,  M.  Seefeldt,  C.  Teodosiu,  S. R.  Kalidindi,  P.  Van Houtte and  E.  Aernoudt: Acta Mater., 49 (2001), 1607.
• 27)   S.  Graff,  W.  Brocks and  D.  Steglich: Int. J. Plast., 23 (2007), 1957.
• 28)   C. R.  Weinberger,  C. C.  Battaile,  T. E.  Buchheit and  E. A.  Holm: Int. J. Plast., 37 (2012), 16.
• 29)   A.  Patra,  T.  Zhu and  D. L.  McDowell: Int. J. Plast., 59 (2014), 1.
• 30)   T.  Yalcinkaya,  W. A. M.  Brekelmans and  M. G. D.  Geers: Model. Simul. Mater. Sci. Eng., 16 (2008), 085007.
• 31)   P.  Franciosi,  L. T.  Le,  G.  Monnet,  C.  Kahloun and  M.-H.  Chavanne: Int. J. Plast., 65 (2015), 226.
• 32)   C.  Paramatmuni and  A. K.  Kanjarla: Int. J. Plast., 113 (2019), 269.
• 33)   R.  Zhou,  A.  Roy and  V. V.  Silberschmidt: Comput. Mater. Sci., 170 (2019), 109140.
• 34)   W. G.  Feather,  S.  Ghorbanpour,  D. J.  Savage,  M.  Ardeljan,  M.  Jahedi,  B. A.  McWilliams,  N.  Gupta,  C.  Xiang,  S. C.  Vogel and  M.  Knezevic: Int. J. Plast., 120 (2019), 180.
• 35)   H.  Zhang,  A.  Jérusalem,  E.  Salvati,  C.  Papadaki,  K. S.  Fong,  X.  Song and  A. M.  Korsunsky: Int. J. Plast., 119 (2019), 43.
• 36)   J.  Singh,  M.-S.  Kim and  S.-H.  Choi: Int. J. Plast., 117 (2019), 33.
• 37)   N. E.  Paton and  W. A.  Backofen: Metall. Trans., 1 (1970), 2839.
• 38)   H.  Conrad: Prog. Mater. Sci., 26 (1981), 123.
• 39)   M. J.  Philippe,  M.  Serghat,  P.  Van Houtte and  C.  Esling: Acta Metall. Mater., 43 (1995), 1619.
• 40)   X.  Wu,  S. R.  Kalidindi,  C.  Necker and  A. A.  Salem: Acta Mater., 55 (2007), 423.
• 41)   J. L. W.  Warwick,  N. G.  Jones,  K. M.  Rahman and  D.  Dye: Acta Mater., 60 (2012), 6720.
• 42)   N.  Benmhenni,  S.  Bouvier,  R.  Brenner,  T.  Chauveau and  B.  Bacroix: Mater. Sci. Eng. A, 573 (2013), 222.
• 43)   A.  Roth,  M. A.  Lebyodkin,  T. A.  Lebedkina,  J.-S.  Lecomte,  T.  Richeton and  K. E. K.  Amouzou: Mater. Sci. Eng. A, 596 (2014), 236.
• 44)   A.  Marchenko,  M.  Mazière,  S.  Forest and  J.-L.  Strudel: Int. J. Plast., 85 (2016), 1.
• 45)   M.  Wronski,  M.  Arul Kumar,  L.  Capolungo,  R. J.  McCabe,  K.  Wierzbanowski and  C. N.  Tomé: Mater. Sci. Eng. A, 724 (2018), 289.
• 46)   K.  Frydrych and  K.  Kowalczyk-Gajewska: Metall. Mater. Trans. A, 49 (2018), 3610.
• 47)   J.  Wang,  M.  Zecevic,  M.  Knezevic and  I. J.  Beyerlein: Int. J. Plast., 125 (2020), 294.
• 48)   J.  Koike: Metal. Mater. Trans. A, 36 (2005), 1689.
• 49)   M. R.  Barnett: Mater. Sci. Eng. A, 464 (2007), 8.
• 50)   T.  Hama,  H.  Nagao,  A.  Kobuki,  H.  Fujimoto and  H.  Takuda: Mater. Sci. Eng. A, 620 (2015), 390.
• 51)   N.  Yi,  T.  Hama,  A.  Kobuki,  H.  Fujimoto and  H.  Takuda: Mater. Sci. Eng. A, 655 (2016), 70.
• 52)   P.  Van Houtte: Acta Metall., 26 (1978), 591.
• 53)   T.  Hama,  Y.  Kariyazaki,  N.  Hosokawa,  H.  Fujimoto and  H.  Takuda: Mater. Sci. Eng. A, 551 (2012), 209.
• 54)   T.  Hama,  H.  Nagao,  Y.  Kuchinomachi and  H.  Takuda: Mater. Sci. Eng. A, 591 (2014), 69.
• 55)   S. R.  Kalidindi: J. Mech. Phys. Solids, 46 (1998), 267.
• 56)   G.  Proust,  C. N.  Tomé,  A.  Jain and  S. R.  Agnew: Int. J. Plast., 25 (2009), 861.
• 57)   H.  El Kadiri and  A. L.  Oppedal: J. Mech. Phys. Solids, 58 (2010), 613.
• 58)   R.  Glüge,  A.  Bertram,  T.  Böhlke and  E.  Specht: J. Appl. Math. Mech., 90 (2010), 565.
• 59)   M.  Homayonifar and  J.  Mosler: Int. J. Plast., 27 (2011), 983.
• 60)   H.  Abdolvand and  M. R.  Daymond: Acta Mater., 60 (2012), 2240.
• 61)   G. I.  Taylor: J. Inst. Met., 62 (1938), 307.
• 62)   R. A.  Lebensohn and  C. N.  Tomé: Acta Metall. Mater., 41 (1993), 2611.
• 63)   H.  Wang,  P. D.  Wu,  J.  Wang and  C. N.  Tomé: Int. J. Plast., 49 (2013), 36.
• 64)   R. A.  Lebensohn: Acta Mater., 49 (2001), 2723.
• 65)   P.  Eisenlohr,  M.  Diehl,  R. A.  Lebensohn and  F.  Roters: Int. J. Plast., 46 (2013), 37.
• 66)   T.  Hama,  T.  Nagata,  C.  Teodosiu,  A.  Makinouchi and  H.  Takuda: Int. J. Mech. Sci., 50 (2008), 175.
• 67)   D.  Peirce,  C. F.  Shih and  A.  Needleman: Comput. Struct., 18 (1984), 875.
• 68)   Y.  Tadano: Int. J. Mech. Sci., 52 (2010), 257.
• 69)   T.  Hama,  N.  Hosokawa and  H.  Takuda: Proc. NUMISHEET 2014, AIP Publishing, Melville, New York, (2014), 692.
• 70)   L.  Kubin,  B.  Devincre and  T.  Hoc: Acta Mater., 56 (2008), 6040.
• 71)   R.  Madec and  L.  Kubin: Acta Mater., 126 (2017), 166.
• 72)   R.  Madec and  L. P.  Kubin: Proc. IUTAM Symp. on Mesoscopic Dynamics of Fracture Process and Materials Strength, Kluwer Academic Publishers, Dordrecht, (2004), 69.
• 73)   R.  Madec and  L. P.  Kubin: Proc. 3rd Int. Conf. on Computational Modeling and Simulation of Materials Part A, Techna Group, Faenza, (2004), 671.
• 74)   T.  Hama,  K.  Ochi,  N.  Kitamura,  H.  Fujimoto and  H.  Takuda: Steel Res. Int., Special Edition, (2011), 1054.
• 75)   Y.  Nakaura,  A.  Watanabe and  K.  Ohori: J. Jpn. Inst. Light Met., 58 (2008), 22 (in Japanese).
• 76)   M.  Ishiki,  T.  Kuwabara,  M.  Yamaguchi,  Y.  Maeda,  Y.  Hayashida and  Y.  Itsumi: Trans. Jpn. Soc. Mech. Eng. Ser. A, 75 (2009), 491 (in Japanese).
• 77)   K.  Yoshida,  A.  Ishii and  Y.  Tadano: Int. J. Plast., 53 (2014), 17.
• 78)   K.  Zhang,  B.  Holmedal,  O. S.  Hopperstad and  S.  Dumoulin: Model. Simul. Mater. Sci. Eng., 22 (2014), 075015.
• 79)   A. P.  Brahme,  K.  Inal,  R. K.  Mishra and  S.  Saimoto: Comput. Mater. Sci., 50 (2011), 2898.
• 80)   L.  Hu,  A. D.  Rollett,  M.  Iadicola,  T.  Foecke and  S.  Banovic: Metall. Mater. Trans. A, 43 (2012), 854.
• 81)   Y.  Jeong,  M. A.  Iadicola,  T.  Gnäupel-Herold and  A.  Creuziger: Acta Mater, 112 (2016), 84.
• 82)   P.  Eyckens,  H.  Mulder,  J.  Gawad,  H.  Vegter,  D.  Roose,  T.  van den Boogaard,  A. V.  Bael and  P. V.  Houtte: Int. J. Plast., 73 (2015), 119.
• 83)   S.  Coppieters,  T.  Hakoyama,  P.  Eyckens,  H.  Nakano,  A.  Van Bael,  D.  Debruyne and  T.  Kuwabara: Int. J. Mater. Form., 12 (2019), 741.
• 84)   T.  Hama,  H.  Fujimoto and  H.  Takuda: Procedia Manuf., 15 (2018), 1808.
• 85)   K.  Inal,  R. K.  Mishra and  O.  Cazacu: Int. J. Solids Struct., 47 (2010), 2223.
• 86)   K.  Zhang,  B.  Holmedal,  O. S.  Hopperstad,  S.  Dumoulin,  J.  Gawad,  A.  Van Bael and  P.  Van Houtte: Int. J. Plast., 66 (2015), 3.
• 87)   K.  Hashimoto,  A.  Yamanaka,  J.  Kawaguchi,  T.  Sakurai and  T.  Kuwabara: J. Jpn. Inst. Light Met., 65 (2015), 196 (in Japanese).
• 88)   A.  Yamanaka,  K.  Hashimoto,  J.  Kawaguchi,  T.  Sakurai and  T.  Kuwabara: J. Jpn. Inst. Light Met., 65 (2015), 561 (in Japanese).
• 89)   F.  Roters,  M.  Diehl,  P.  Shanthraj,  P.  Eisenlohr,  C.  Reuber,  S. L.  Wong,  T.  Maiti,  A.  Ebrahimi,  T.  Hochrainer,  H.-O.  Fabritius,  S.  Nikolov,  M.  Friak,  N.  Fujita,  N.  Grilli,  K. G. F.  Janssens,  N.  Jia,  P. J. J.  Kok,  D.  Ma,  F.  Meier,  E.  Werner,  M.  Stricker,  D.  Weygand and  D.  Raabe: Comput. Mater. Sci., 158 (2019), 420.
• 90)   F.  Han,  M.  Diehl,  F.  Roters and  D.  Raabe: J. Mater. Process. Technol., 277 (2020), 116449.
• 91)   H.  Zhang,  M.  Diehl,  F.  Roters and  D.  Raabe: Int. J. Plast., 80 (2016), 111.
• 92)   T.  Hama,  K.  Hirano,  R.  Matsuura,  Y. P.  Korkolis and  H.  Takuda: Proc. NUMIFORM2019, TMS, Pittsburgh, PA, (2019), 601, CD-ROM.
• 93)   T.  Hama,  T.  Sakai,  Y. P.  Korkolis and  H.  Takuda: J. Phys.: Conf. Ser., 1063 (2018), 012122.
• 94)   E.  Nakamachi,  N. N.  Tam and  H.  Morimoto: Int. J. Plast., 23 (2007), 450.
• 95)   I.  Tikhovskiy,  D.  Raabe and  F.  Roters: Mater. Sci. Eng. A, 488 (2008), 482.
• 96)   P.  Van Houtte,  J.  Gawad,  P.  Eyckens,  A.  Van Bael,  G.  Samaey and  D.  Roose: JOM, 63 (2011), 37.
• 97)   I.  Tikhovskiy,  D.  Raabe and  F.  Roters: J. Mater. Process. Technol., 183 (2007), 169.
• 98)   K.  Inal,  P. D.  Wu and  K. W.  Neale: Int. J. Plast., 16 (2000), 635.
• 99)   M.  Zecevic,  R. J.  McCabe and  M.  Knezevic: Mech. Mater., 84 (2015), 114.
• 100)   K.  Yoshida,  T.  Ishizaka,  M.  Kuroda and  S.  Ikawa: Acta Mater., 55 (2007), 4499.
• 101)   E.  Nakamachi,  T.  Honda,  H.  Kuramae,  Y.  Morita and  H.  Morimoto: Int. J. Mech. Sci., 87 (2014), 89.
• 102)   S.  Ikawa,  M.  Asano,  M.  Kuroda and  K.  Yoshida: Mater. Sci. Eng. A, 528 (2011), 4050.
• 103)   W.  Muhammad,  A. P.  Brahme,  U.  Ali,  J.  Hirsch,  O.  Engler,  H.  Aretz,  J.  Kang,  R. K.  Mishra and  K.  Inal: Mater. Sci. Eng. A, 754 (2019), 161.
• 104)   D.  Steglich,  Y.  Jeong,  M. O.  Andar and  T.  Kuwabara: Int. J. Solids Struct., 49 (2012), 3551.
• 105)   F.  Hiura and  M.  Niewczas: Scr. Mater., 106 (2015), 17.
• 106)   M.  Zecevic and  M.  Knezevic: Int. J. Plast., 105 (2018), 141.
• 107)   H.  Zhang,  X.  Dong,  Q.  Wang and  Z.  Zeng: Comput. Mater. Sci., 54 (2012), 208.