Advancements in and Applications of Crystal Plasticity Modelling of Metallic Materials

Abstract

Integrated Computational Materials Engineering (ICME) is a set of methodologies utilized by researchers and engineers assisting the study of material behaviour during production processes and/or service. ICME aligns with societal efforts for the twin green and digital transitions while improving the sustainability and cost efficiency of relevant products/processes. A significant link of the ICME chain, especially for metallic materials, is the crystal plasticity (CP) formulation. This review examines firstly the progress CP has made since its conceptualization and secondly the relevant thematic areas of its utilization and portraits them in a concise and condensed manner. CP is a proven tool able to capture complex phenomena and to provide realistic results, while elucidating on the material behaviour under complex loading conditions. To this end, a significant number of formulations falling under CP, each with their unique strengths and weaknesses, is offered. It is a developing field and there are still efforts to improve the models in various terms. One of the biggest struggles in setting up a CP simulation, especially a physics-based one, is the definition of the proper values for the relevant parameters. This review provides valuable data tables with indicative values.

1. Introduction

Nowadays, the twin green and digital transitions force society to change products, processes, and “lifestyle”, targeting circularity and taking action towards raw material preservation. Integrated Computational Materials Engineering (ICME) is aligned with these goals as it can be applied in the development of new alloys while addressing the improvement in sustainability and cost efficiency of alloys/products, taking into consideration all requirements from their design and processing to their service.

Gradually, with the development of computer processing power, the complexity of systems and geometries has significantly increased to the point where ICME approaches today have become more realistic even for materials/systems of great complexity.

This dynamic environment creates a need for developing materials and processes guided by a deep understanding of the mechanisms and phenomena occurring on a microscopic level. This can be achieved through the coupling of different computational approaches. Such approaches, involving among other things the thermodynamic behaviour of the system and the deformation behaviour of the material, need to occur in fast, cost-effective, and sustainable development cycles. The structural design process needs to expand now to also include the material instead of narrowly focusing on the geometry of the part, as has been the case [1].

ICME provides improved understanding of the mechanics of the material, through which two main aspects of a material’s lifecycle are addressed, namely the production process and the operational behaviour. The potential benefits of utilizing ICME in the process optimization are as follows: (i) development of products with higher added value, (ii) process optimization resulting in reduced energy consumption, and lastly, (iii) improved product properties, which also target reduced raw material consumption. On the other hand, through the understanding of the behaviour during service, the requirements that the microstructure needs to meet become more clear, potentially leading to light-weight design and increased safety without compromising on the part’s life. At the same time, the operating conditions for materials constantly gets more demanding. Additionally, to tailor the properties of materials to their respective applications, microstructures are becoming more complex. Another emerging challenge derives from circularity. The increasing level of impurities from scrap, following the recycling rates of materials, may negatively affect the microstructure through the build-up of residual elements, with some of them being critical for the integrity of the microstructure.

These potential benefits highlight the alignment between the adaptation of ICME approaches and the Sustainable Development Goals set by the United Nations, namely goals 4—Quality Education; 8—Decent work and economic growth; 9—Industry, innovation and infrastructure; 12—Responsible consumption and production; and 17—partnerships for the goals.

There are several methodologies under the umbrella of ICME. These can either focus on the “process–microstructure” (e.g., phase field) or the “microstructure–properties” (e.g., crystal plasticity) relations. The coupling of the two approaches allows the consecutive study of the effect of a process on the microstructure (e.g., thermal treatment) and then the use of the affected microstructure as input and prediction of its mechanical properties, enabling the tailoring of the properties to each application. Indeed, these approaches have been a subject of study for decades, with examples such as the Hall–Petch effect [2,3] and the work of Rhines [4], where critical microstructural parameters, such as the volume fraction of phases, have been connected with the mechanical properties. CP modelling falls under the microstructure–properties category [5], allowing the impact of factors such as the β-precipitation of AA6061 on mechanical properties to be obtained [6], yet it can also be used to examine the effect of mechanical loading on grains.

The CP approach is based on two fundamental axioms. Firstly, each grain is considered a continuum, and secondly, the total deformation is controlled by the interaction of the active slip systems of each grain. The latter is identified as the most important differentiator between CP and the other computational approaches [7]. Yu et al. [8] attribute the significance of the CP methodology to (i) its strong mathematical foundations, (ii) the correlation between slip systems and geometry, and lastly (iii) the high-fidelity algorithms solving such problems. These factors, in combination with the knowledge accumulated around continuum solid mechanics [9], and CP’s lower costs when compared to physical experiments [10], have enabled the prediction of mechanical properties in relation to microstructural characteristics under complex loading conditions [11], making CP a widely used tool to study the mechanical responses of alloys [8], such as the deformation mechanics, e.g., during metal forming processes [12,13,14,15] or Forming Limit Diagrams, FLDs [10,16]. Additionally, when CP is coupled with other methodologies, the frameworks can expand to capture more phenomena such as the behaviour under cyclic loading [17]. This enables an enhanced understanding of how the material behaves during service and how to effectively tailor the process. Through this understanding, the development of complex, tailored materials with desired/advanced mechanical properties can be realized [16].

As mentioned, the deformation of each grain is influenced by the orientation of the slip systems in relationship with the applied load. Polycrystalline materials have grains with different orientations. This combination of the different orientations present in the material is known as crystallographic texture. This affects how the applied loading influences the material’s total deformation due to the synergetic slip of all the grains deformed. This explains why the quantification and study of anisotropy, especially in relation to crystal texture, is of importance [18,19].

Texture is measured in reference to a coordination system (c.s.). Usually, the selection is between the three commonly used c.s., which are as follows [20]:

i.

The c.s. defined by the deformation axes (e.g., for the case of rolling the RD (rolling direction), TD (transverse direction), and the ND (normal direction).

ii.

The c.s. based on the crystallography, which rotates along the grains during deformations. This is especially useful for single-crystal systems.

iii.

An external c.s. free of deformation, which is used for simplicity purposes (e.g., x, y, z).

Modern CP models, utilizing either a Finite Element Solver (FEM) or a Spectral Element Solver (SEM), can successfully predict, among others, the TRIP and TWIP effects of HSLA steels and the fatigue behaviour of dual-phase Ti-alloys [21] and Co-Cr alloys [22]. Other applications include adiabatic shear banding of FCC crystals [23] and the effect of grain size on the flow stress [24], with more models being tested even today, indicatively [25,26], and even the effect of activated slip systems on fatigue behaviour [27]. Additionally, it is often the case that the research is focused, e.g., [28], on micro-mechanical testing. This is a result of three main factors. Firstly, models need validation tests on single crystals, secondly, due to the lesser complexity of the system, they have lower computational cost, and, lastly, conducting physical experiments, under controlled conditions, on this scale is complex. Mechanical testing on single crystals is rarely preferred and most models rely on macroscopical quantities, such as stress–strain curves, for validation [29].

The importance of the CP methodology today is evidence by the breadth of research topics and alloying systems. Indicatively, some recent subjects of study, categorized by material, are as follows:

i.

Aluminium (Al) alloys: The effect of grain topology on crack propagation for aerospace-grade materials [30], modelling of the Continuous Dynamic Recrystallization (CDRX) [31], and the behaviour of such alloys under cyclic compression loading [32], as well as their behaviour under high-cycle fatigue following 3D printing [33]. Moreover, they have been successfully utilized to study warm forming conditions [34] and even the Bauschinger effect on single crystals [35].

ii.

Iron (Fe) and steel alloys: The effect of non-metallic inclusion on the steel failure, e.g., of a 16MnCrS5 steel, has been studied [36]. The effect of martensite on forming has been studied [10,14]. Deep drawing of dual-phase (DP) steels [37], sintering of transformation-induced plasticity (TRIP) [38], and twinning-induced plasticity (TWIP) [39] steels have also been subjects of study. Recently, a formulation modification was proposed for studying the effect of the grain size and the slip system interaction in austenitic stainless steels [40]. Another example is the study of the behaviour of an austenitic steel under multiaxial loading [41].

iii.

Magnesium (Mg) alloys: The evolution of twin formation [42,43], the effect of heat treatments on the mechanical behaviour of the twins [44], the evolution [45] and effect [46] of twins in low-cycle fatigue conditions, as well as the effect of grain size on texture evolution under mechanical loading have been studied [47]. Cheng et al. expanded a model for it to be able to capture the effect of hotspots in the local deformation of twin bands [48]. A recent review detailing the applications for Mg alloys can be found in [49].

iv.

Titanium (Ti) alloys: The anisotropic plastic deformation of commercial alloys has been studied [50]. The behaviour of dual-phase Ti-alloys has also been studied [21]. Additionally, the results regarding the micro-mechanic response of a TNM alloy have been validated through nano-indentation [51]. Results from CP modelling have also been compared to High-Resolution EBSD (Electron Back-Scatter Diffraction) and High-Resolution DIC (Digital Image Correlation) during micro-slip [52] and stress/strain localization [53]. Another difficulty in simulating the response of highly anisotropic materials is the effect of kink bands during forming, for which an approach for Ti alloys has been proposed [54]. Lastly, the texture evolution during cold forming has also been studied, e.g., [55].

v.

Nickel (Ni) super-alloys: The effect of nano-indentation on the distribution of dislocation density of a single-crystal [56], mechanical testing on specimens produced by direct casting [57], applying the weak link methodology in order to evaluate the fatigue behaviour of the material, while taking into account the part size [58], multiaxial fatigue behaviour [59], and the effect of grain size on fatigue behaviour have been studied [60]. Apart from fatigue behaviour, the creep behaviour has also been studied, e.g., [61,62].

vi.

Other alloys: The behaviour of copper (Cu) alloy oligo-crystals under mechanical shock [63], the effect of additive manufacturing on the evolutions of residual stresses of tungsten (W) alloys [64], the low-cycle fatigue behaviour of cobalt–chromium (CoCr) alloys [22], and, lastly, the evolution of texture during severe plastic deformation (SPD) in high-entropy alloys through an expansion of the Taylor model have also been studied.

However, despite the constant progress of CP and the utility of its results, determining the input parameters is still a matter of study [9,26].

This review is structured as follows: After the introduction to the subject (Section 1), the motivation for this work and our contribution to the field is presented (Section 2). This contribution does not only reflect on how the crystal plasticity methodology has evolved over the years (Section 3), but in addition provides a discussion regarding the proper model selection including the constitutive equations (phenomenological or physics-based), the solver (Finite Elements or Fast Fourier Transformations), and the most used packages (Section 4). Then, carefully selected parameters crucial for the execution of a CP simulation, which are absolutely necessary for every researcher working in the field, are gathered, assessed, and coded in very handy table form (Section 5). Furthermore, the required critical steps (i.e., geometry, material, and loading conditions) for a successful simulation execution also considering important parameters affecting the process are presented and discussed as well. Finally, this work concludes with our closing remarks (Section 6).

2. Motivation and Contribution

Crystal plasticity (CP) has significantly evolved since its conceptualization. What began as a tool for studying the deformation of single crystals under uniaxial deformation has become capable of capturing complex phenomena in multiphase polycrystals, even in materials with a gradient crystal structure under complex loading conditions, while at the same time improving computational efficiency and accuracy. There are several papers that perform a review on CP and relevant methodologies, indicatively, Beyerlein et al. [65] reviewed extensively the various models and the difference in their mathematical formulation, Lucarini et al. [66] followed a very detailed, textbook-like approach in explaining Fast Fourier Transformations in relation to micro-mechanics, Yaghoobi et al. [49] details the challenges of magnesium alloy modelling and highlights the relevant literature, whereas, Roters et al. [67] showcased the extent of the possible applications. These reviews provide the reader with a detailed description of the mathematical formulation and mostly focus on literature directly relevant to the specific subject. This review differentiates from those mentioned above (and similar approaches) in the following aspects:

i.

It systematically presents the evolution and recent advancements of the CP method through selected milestones and applications with an emphasis on metallic crystalline materials while providing explanation of the basic formulations; in-depth mathematical analysis is omitted.

ii.

It is designed to be as short as possible while showcasing most of the critical aspects of running a CP simulation and providing relevant resources.

3. Evolvement of Crystal Plasticity Approaches

Taylor [68] described the slip of single-crystal materials and he was the first to present crystal plasticity in 1982 [69]. By 1993, Beaudoin et al. [70] were already using three-dimensional (3D) models. Figure 1 [71,72,73,74,75,76,77,78] highlights some, non-exhaustive, milestones in the evolution of the mathematical models describing the plastic behaviour of materials from the first mention of dislocations up to today. The formulations connecting plastic and hardening behaviours to dislocations only became accepted in 1986 by Estren and Kubin [79]. The dislocations were introduced in the study of polycrystalline materials by Arsenlis et al. [80] in the early 21st century and quickly started developing [81]. These early adaptations restricted dislocation movement across multiple grains [81], with this limitation to be overcome by Ma, Roters, and Raabe [82,83].

The first attempts to utilize CP were restricted to simple geometries and significant assumptions were made to simplify the problem. Gradually, the complexity of the systems and geometries has significantly increased to the point where they have become realistic even for more complex materials. The evolvement of CP models is partly dependent on atomic-scale models, since their results can help improve the CP formulation. This can be seen by the expansion of the phenomenological models to include BCC materials [84,85], from FCC where they were originally applied. Additionally, through atomic-scale studies, it was proven that the work-hardening of FCC metals was mainly affected by the collinear dislocation interaction [86], i.e., the result of attraction of dislocations moving on the same slip plane. This created a need to incorporate the dislocation interaction into CP models, resulting in what is now known as physics-based CP. Parameters from various scales, such as dislocations, are used in a statistical manner to calculate the final macroscopic deformation and stress distribution. Due to such parameters being statistically integrated, instead of being explicitly studied, physics-based CP still falls under continuum solid mechanics (mesoscale modelling). The interaction of mobile dislocations with dislocation forests was studied via a multiscale approach, combining the dislocation density with the resolved shear stress on each slip system [87]. This was followed by efforts to capture the effect of deformation-induced heat on the material behaviour for small deformations [88]. CP is also a capable tool for the study of gradient crystal materials [89,90].

There are currently three multiscale modelling approaches to study alloys [91]. Firstly, the embedded approach, according to which the lower-scale models are embedded in the higher-order ones, which has high accuracy and is often used in physics-based CP. According to the hierarchical approach, a separate model is used at each scale and the outputs of lower scales are iteratively used as input for the higher scales, and lastly, the Adaptive Sampling Method (ASM) according to which, the lower-scale models executed only once, and their results are used in many different scenarios.

Heterogeneities within the material have been successfully studied, firstly by Hashin et al. [92] using a framework for homogenized and fully elastic (linear) microstructural models [93]. A shift from the iso-strain homogenization scheme to the Relaxed Grain Cluster (RGC) model [94,95] can allow the study of the interaction between the grains within the representative volume element (RVE). RGC is a homogenization technique used to reduce the computational cost while retaining sufficient information [95] at an improved speed [94]. This is achieved through clustering grain behaviour under mechanical loading yet allowing a degree of freedom to minimize the energy of the cluster. The impact of the homogenization scheme on texture evolution was studied by Eisenlohr et al. [96].

3.1. General Taxonomy of Crystal Plasticity Models

There are many ways to categorize crystal plasticity models. A branch of the categories described by Trusov et al. [97] are depicted in Figure 2. They can be categorized as follows: Statistical, which is better suited for uniformly deformed volumes [98,99,100], Self-Consistent, which can study steady-state deformation [65], and Direct models, which are mostly used for nonuniform deformations [97]. The Self-Consistent approach was considered separate to crystal plasticity [101] until Lebensohn et al. [102] integrated them. Each approach can follow a local or a non-local approach while calculating the critical shear stress. The local approach takes into account the stress applied in each calculation point, whereas the non-local approach utilizes either [103] the deformation gradient (e.g., [24,81,104]) or continuum dislocation dynamics (e.g., [56,105,106,107,108,109,110,111,112]).

Taking the strain gradient into consideration, studying the effect of the RVE size on the results is allowed. Furthermore, non-local strain-based models can be further categorized as lower or higher order, depending on the formulation of the boundary conditions [53,54]. Lastly, they are described as mean or full field (e.g., [113]) by the selection of the solver, as discussed in Section 4.2. Utilization of a full-field formulation translates to the ability to study local heterogeneities in the macroscopical mechanical behaviour of the materials [114] and a prediction of localized stresses [18].

3.2. Approaches for Further Improvement of CP

The CP framework is constantly being upgraded and optimized in its ability to utilize and predict complex phenomena, e.g., [25,26]. Other approaches for improvement can be categorized as follows:

i.

Improving execution time and accuracy. As mentioned, CP problems can be solved with either a FEM approach or a FFT approach. The CP-FFT methodology is significantly faster, mostly due to the low mesh sensitivity and the mathematical formulation. However, efforts have been made to further improve their computational efficiency [115] and accuracy [116]. Ling et al. [117] decreased the mesh sensitivity of the FEM solvers, allowing the application of a coarser grid while retaining result accuracy and indirectly accelerating the calculation process. Other efforts are focused on the optimization of resource allocation [118] via the modification of CP to run on a GPU and creating a hybrid CPU-GUP architecture to increase computational efficiency, while others focus on simplifying the FFT controls [114]. Zecevic et al. [119] introduced an LS-EVP-FFT (large-strain elasto-viscoplastic) model, which allows for increased accuracy in quantifying the effect of anisotropy on the creation of hotspots, through the modification of Green’s operator and introduction of more grids per material point. Admal et al. [120] depicted grain boundaries as a special category of geometrically necessary dislocations (GNDs) to better depict their impact. Recently, Romanov et al. [11] presented a statistical model, thus reducing execution time, capable of studying the ECAP test. Another topic of interest is the evaluation of the stability of the developed models [121,122].

ii.

Combining the CPFEM and CPFFT approaches. Yu et al. [8] combined the two solvers (FEM and FFT), allowing the expansion of possible applications. Alharbi et al. [123] used the results of the FFT approach as input for the CPFEM, achieving faster results.

iii.

Optimizing input geometry. Apart from the geometry and the number of grains, an important parameter of the RVE file has proven to be the total volume of the synthetic microstructure [124,125]. This topic will be discussed in detail in Section 5.1.

iv.

Streamlining the parameter identification process. This can be achieved through in situ tensile tests [126] or genetic algorithms [127]. Furthermore, some CP models are being designed to be material-independent [26] to increase the spectrum of potential applications.

v.

Creation of tailored models. To achieve higher accuracy, some models are being created to be tailored to specific materials, e.g., [26], where a material-invariant approach of mesoscale parameters has been proposed, or processes, e.g., for the effect of the hardening parameters on the crystallographic texture evolution during the rolling of aluminium alloys [128].

vi.

Other efforts focus on expanding the field of potential applications, e.g., to other crystal systems by accounting for the dislocation slip and twinning of Hexagonal Closed Packed (HCP) systems [129] or even capturing the effect of sample size to the yield stress through an embedded, sub-routine for Discrete Dislocation Dynamics in a CPFEM model [28].

3.3. Coupling CP to Machine Learning Algorithms

A different approach is the utilization of Artificial Intelligence (AI) and machine learning (ML) methodologies in conjunction with crystal plasticity modelling. A description of such approaches can be found in [130], but they will not be discussed in detail here. The usage of these approaches in the field of materials modelling is expected to grow [131]. One of the most important issues preventing their widespread application [132] is the volume of data needed for training and validation [127]. The coupling of these tools with CP usually has one of the following goals:

i.

Model calibration. This approach can solve one of the most important issues of applying CP models, the determination of the correct values for the input parameters. Among others, Galan Lopez et al. [18], Chakraborty et al. [19], Plowman et al. [10], and Sedighiani et al. [127] have all performed a calibration of the input parameters needed for CP simulations through an ML algorithm, applied to data from tensile testing. A similar approach was followed by Sahoo et al. [29] for the phenomenological model. There are examples of calibration being performed under cyclic loading conditions [133].

ii.

Improved accuracy while maintaining the computation cost. Due to the increasing complexity of the Partial Differential Equations (PDEs), resulting from raising the accuracy and realism of the models, the implementation of such an approach is also expected to show an increase in the following years [134]. Adopting methodologies such as Convolutional Neural Networks (CNNs) can result in a reduction in the computational cost of macroscopic properties by 99% when compared to a two-dimensional (2D) CP simulation [135]. Applying Deep Neural Networks (DNNs) for the local response of an FFT solver, carried out by Mianroodi et al. [134], lead to an acceleration in the time needed for obtaining the results by ×8300 times for heterogeneous materials. Another example of achieving high accuracy, while maintaining a low computational cost, is the prediction of the local response in industrial aluminium alloys using a CNN-CPFEM methodology [136].

iii.

Surrogate modelling. Some indicative works are the one of Saidi et al. [137], in which they utilized the Taylor model to train an ML algorithm for the rolling of aluminium alloys. Khorrami et al. [138] calibrated a CNN model for tensile tests with results from VSPC (Viscoplastic Crystal plasticity), and accurately predicted the von Mises stresses.

4. Methodologies

This section discusses the main difference in the formulation of phenomenological CP models and physics-based CP models, includes a comparison between the Fast Fourier Transformation (FFT) and the Finite Element (FEM) solvers, and, lastly, refers to some indicative commonly used CP packages.

4.1. Phenomenological and Physics-Based CP

Regardless of the chosen model category, the shear rate of each slip system is correlated with the deformation of the whole volume under consideration. Equation (1) [20] provides the formulation for the speed of the plastic deformation gradient (L_p), where m is the vector normal to the slip system, a represents the slip system, n is the total number of independent slip systems (i.e., 12 for FCC, 5 for BCC, and 3 for HCP), and ${\dot{\mathsf{\gamma }}}^{\mathsf{\alpha }}$ is the shear rate of the ath slip system.

$${\mathrm{L}}_{\mathrm{p}}=\sum _{\mathrm{a}=1}^{\mathrm{n}}{\dot{\mathsf{\gamma }}}^{\mathsf{\alpha }}·{\mathrm{m}}^{\mathsf{\alpha }}\otimes {\mathrm{n}}^{\mathsf{\alpha }}$$

(1)

The expression of the shear rate is one of the most fundamental differences between the phenomenological and the physics-based formulations. For the phenomenological models, the shear rate is given by Equation (2). In Equation (2), the variable m differs from the one in Equation (1) as it a case-specific parameter with no physical meaning, modified to achieve the best fit with experimental observations.

$${\dot{\mathsf{\gamma }}}^{\mathsf{\alpha }}={\dot{\mathsf{\gamma }}}_0·{\left|{\displaystyle \frac{{\mathsf{\tau }}^{\mathsf{\alpha }}}{{\mathsf{\tau }}_{\mathrm{c}}^{\mathsf{\alpha }}}}\right|}^{{\displaystyle \frac{1}{\mathrm{m}}}}·\sinh \left({\mathsf{\tau }}^{\mathsf{\alpha }}\right)$$

(2)

${\mathsf{\tau }}_{\mathrm{C}}^{\mathsf{\alpha }}$ is the critical shear stress on the slip plane, while τ^α is the saturation shear stress [139] and a is the number of slip systems in the crystal.

On the other hand, the realism provided by the formulation incorporating physics laws leads to a more a more complex expression, with a significant number of parameters [127]. An example of such a formulation is the work of Kords [140], shown in Equation (3). It can be seen that, additionally to stress-related terms, it includes quantities such as the dislocation density (ρ) and speed (${\mathrm{v}}_0$), the burger vector (b), solid solution activation energy (Q_s), the Bolman’s constant (k_B), and temperature (T), plus two parameters (p, q) controlling dislocation movement behaviour. Most of these terms have their own mathematical formulation, resulting in more parameters being involved.

$$\dot{\mathsf{\gamma }}=\mathsf{\rho }·{\mathrm{b}}_{\mathrm{s}}·{\mathrm{v}}_0·{\mathrm{e}}^{\left[-{\displaystyle \frac{{\mathrm{Q}}_{\mathrm{s}}}{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}}·{\left\{1-{\left({\displaystyle \frac{\left|{\mathsf{\tau }}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right|}{{\mathsf{\tau }}_{\mathrm{s}\mathrm{o}\mathrm{l}}}}\right)}^{{\mathrm{p}}_{\mathrm{s}}}\right\}}^{{\mathrm{q}}_{\mathrm{s}}}\right]·\mathrm{s}\mathrm{g}\mathrm{n}\left(\mathsf{\tau }\right)}$$

(3)

Considering complex interactions, even in a statistical manner, is a computationally expensive process. The mathematical foundations for reducing the cost of capturing the dislocation movement in CP modelling were laid by Arsenlis et al. [105] and Ganghoffer et al. [141], who developed differential equations for the calculation of the dislocations, with higher degrees of freedom. Later, Ma, Roters, and Raabe [142] suggested an alternative, which is more efficient under complex loads. Subsequent developments further enhanced the accuracy of dislocation movement [56].

Physics-based CP models, although they can receive input parameters (e.g., atomic volume and dislocation density) from lower scales (i.e., atomic scale), still examine the grain as a continuum body. Phenomena such as dislocation interaction are considered in a statistical manner [143] instead of being individually tracked. Thus, just like phenomenological models, physics-based ones also fall under the mesoscale modelling category, since they consider each grain a continuum.

4.2. FEM vs. Spectral Methods

As mentioned, each of the CP models can be solved either with the FEM or with a spectral approach (FFT), which was proposed as an alternative [144], resulting in CPFFT and CPFEM categories. Even though FFT was suggested relatively early, its use on polycrystalline materials is rather new, firstly proposed by Lebensohn [145], from where it was expanded to anisotropic materials with significant texture [91]. Spectral solvers (FFT) have two main basic advantages over the FEM [146], which are as follows [25,115,147,148,149]:

i.

Computational efficiency: They are faster, by around an order of magnitude, and can therefore achieve higher resolution/grid density on the RVE at the same time.

ii.

Due to their low mesh sensitivity, they enable the use of a coarser grid without sacrificing result accuracy [93,146]. An example of this is [38], where they used one calculation point per grain.

On the other hand, these benefits involve some trade-offs that need to be considered for each case. Firstly, due to the periodic boundary conditions applied it is difficult to study the materials near the free surfaces. This is a result of the fact that the volume under study is virtually extended infinitely around itself in three-dimensional space [150].

To demonstrate the periodic boundary conditions, Figure 3 shows an RVE created based on the Voronoi methodology containing 30 grains with a grid of 500³ voxels visualized through ParaView [151,152] (Figure 3a). Grain morphology is also shown (Figure 3b–d) whereas in Figure 3e,f, a face of the original RVE is repeated in space. Since the material is constantly in contact with other identical volumes, RVEs studied under the periodic conditions’ assumptions have no free surfaces.

To overcome the limitation introduced by the lack of free surfaces, Maiti and Eisenlohr [153] proposed to enclose the RVE of interest in another material, simulating the interaction with the free surface, resulting, however, in a significantly more complex RVE. Additionally, it can only work with the use of a normal/regular (x = z = y) grid [154,155] and the remeshing processes have proven to be rather complex [156]. It should be noted that the issue regarding the grid sensitivity of the FEM solver can in some cases be overcome by utilizing the micro-morphic crystal plasticity technique [23]. A more detailed comparison of the approaches can be found in [157].

4.3. Indicative Available CP Packages

A significant number of CP-related software/codes exist, with research groups trying to develop in-house/custom packages, with some of them being openly available. The most widely used, to the best of the authors knowledge, are as follows:

i.

DAMASK, developed by the Max Plank Institute for Sustainable Materials GmbH—R&D centre, released in 2011 [67,158].

ii.

PRISMS, developed by the University of Michigan [159,160].

iii.

CAPSUL and FFTMAD, both developed by the IMDEA Research centre, released in 2012 [161,162].

iv.

ρ-CP, co-developed by the Indian Institute of Technology and Georgia Institute of Technology, released in 2023 [163].

v.

AMITEX_FFTP, developed by the coalition Maison de la Simulation, consisting of the French National Centre for Scientific Research, the French Alternative Energies and Atomic Energy Commission (CEA), Université Paris-Saclay, and Université Versailles Saint-Quentin, also released in 2023 [164].

vi.

Other packages available in the literature are the ones found in [66,104,165] as well as in [166], implementing CPFFT and CPFEM, respectively [167].

Table 1. Comparison between the selected packages.

	Solver		Formulation		Access
	FEM	FFT	Phenomenological	Physics-Based	Open Access (OA)/ Commercial Licence (CL)
PRISMS	X		X		OA
CAPSUL	X		X	X	CL
FFTMAD		X			OA/CL
ρ-CP	X			X	OA
AMITEX		X			OA
DAMASK	X	X	X	X	OA

provides a concise overview of the packages in terms of (i) the solver used, FEM vs FFT as previously discussed, (ii) the formulation used, phenomenological or physics based, and, lastly, (iii) if they are open access or require a commercial licence.

5. Critical Aspects for CP Simulations

This section is focused on the most critical aspects of setting up a CP simulation. These aspects include the RVE dimensions, the parameters that describe the material behaviour, and the loading conditions.

5.1. RVE Modelling

5.1.1. RVE Creation

Proper selection of geometry-related parameters is the first, and crucial [18], step towards the design of the synthetic microstructure. The most common ways to create an RVE are either through (i) a 1:1 representation based on data provided by techniques such as serial sectioning [124], Atom Probe Tomography (APT), 3D EBSD, assisted by software, such as NEPER [168], DREAM 3D [169], and Voro++ [170], through the cellular automata approach [124] or (ii) creating the RVE based on typical data obtained by standard metallography procedures, via the aforementioned software, most commonly utilizing the Voronoi methodology [171]. These two approaches are not interchangeable in their results. Indicatively, Tu et al. [22] studied the effect of the grain morphology of the RVE under fatigue and found that the RVEs with an accurate morphology had a 10% better fit with the experiments in comparison to Voronoi-based microstructures. As explained in [93], the experimentally produced RVEs have some deviation in capturing grain-boundary-related phenomena, whereas the computationally calculated RVEs mostly consider grain boundaries to be planes. Moreover, during serial sectioning, the information between the sections is lost [124].

There are studies regarding measurement accuracy, e.g., [172,173,174] through focused ion beam-scanning electron microscopy (FIB-SEM) and automation for larger scales and APT [175] for smaller scales, which, however, greatly increases the cost of the RVE extraction. An extended review on advanced Microstructure Characterization and Reconstruction (MCR) techniques was performed by Bostanabad et al. [176].

5.1.2. RVE Size

A critical parameter when designing the RVE is the size. Specifically, small volumes, depending on the formulation used, may behave like oligo-crystals, losing information about the interaction between grain and phase boundaries, while large volumes have a higher computational cost [177]. For most cases focusing on deformation, a volume of a few hundreds of grains can suffice [24] or equivalently an edge of the RVE equal to a few hundred micrometres [128]. Studies involving statistical quantities (e.g., crystallographic texture) require a significantly larger number of grains [155].

Size also incorporates the thickness of the RVE. Macroscopically, there is evidence [178] suggesting that the fit with the validation points may be sufficient regardless of the thickness. However, the thickness can influence the accuracy on a local level, especially when the comparison is between a 2D and a 3D volume [179,180,181]. Mirhosseini et al. [182] schematically demonstrate the effect of the size for 2D and 3D RVEs with grains in the range between 9 and 225 and from 8 to 125, respectively. Specifically, it was shown that 3D RVEs require less grains in total for a similar accuracy. An approach for transferring the results of 2D experiments to the 3D space is described in the recent work of Tseng et al. [179]. The effect of the selected size on the calculated mechanical properties, as well as the differences between RVEs and statistical RVEs (S-RVEs or SVEs), is discussed in [183].

For each case, the proper size needs to be determined a posteriori through a dependency analysis [124], especially in the cases of fatigue analysis, where size has been found to have a significant impact on the accuracy of the fatigue prediction [58]. This means using progressively larger volumes until the results become stable.

5.1.3. Grid Density

Grid, or mesh, density refers to the number of calculation points within the volume of a given RVE. The required grid density mostly depends on the scale of interest. Specifically, if focus of the study is on macro-mechanical properties (e.g., stress–strain curves) a few calculation points per grain are considered enough, whereas at the mesoscale 10²–10³ [156] calculation points per grain would be needed, significantly increasing calculation time.

Texture as a statistical quantity is considered macroscopical and, thus, it is considered to be mesh-independent if there are calculation points in all grains. Yet, Yu et al. [184] observed that when utilizing only some of the original voxels as an RVE (an approach called sub-modelling), effectively increasing the resolution in this area, the prediction of micro-texture and in grain variations was enhanced.

A grid sensitivity analysis is required to determine the proper grid size of each system. Moreover, due to the heterogenous nature of the local deformation and strain localization, the grid is subjected to degradation [156]. To counter this, an adaptive remeshing approach, i.e., changing the mesh along with RVE size, has been proposed.

5.2. Material Related Input Parameters

Defining the input parameters (calibration) for the simulations remains a challenging [9,26,29] and time-consuming task [185]. Some parameters, such as self-diffusion, can be easily found in the relevant literature, e.g., [186], whereas others have proven to be more difficult. This is due to their dependency on loading conditions and material properties [26,187], chemistry, e.g., [21], as well as on the methodology/formulation utilized.

The efforts made to surpass these challenges [127,188] and their quantification can be considered to be a separate field of study [19]. Zeng et al. [189] quantified the effect of each parameter of their model. Moreover, for some models, e.g., [127], it has been proven that each one of the various parameters for the physics-based model takes a single value for each case. As was discussed in Section 4.1, the use of physics-based models entails the quantification of a significant number of parameters. The parameters used in such models can be categorized as follows:

i.

Dislocation-Related Parameters. Here, included are parameters such as the dislocation density per dislocation type, the dislocation glide and transmissivity, and various other parameters that assist with the statistical representation of the dislocations.

ii.

Thermodynamical Parameters. These can be related to the solid solution, the diffusion, etc.

iii.

General Parameters. These can be, for example, the mechanical parameters of the material.

iv.

Case-specific parameters need to be calibrated for each case (fitted to experimental data). Parameters that do not necessarily have a physical meaning need to be fitted in each specific case, accounting both for the material and the process.

One of the most difficult-to-measure parameters is the dislocation density, usually measured in m/m³ [190]. It is the carrier of plastic deformation [20] and thus heavily influences the mechanical response of the material. Some formulations, e.g., [140], differentiate between edge and screw dislocations and further divide them into dipolar (immobile) and monopolar (mobile) dislocation density as well as positive and negative. In one study [105], monopolar ρ is also found as GNDs (Geometrically Necessary Dislocations), while the dipoles are found as SSDs (statistically stored dislocations).

Their importance is not restricted to mechanics, since phenomena such as grain growth are heavily dependent on them. Recently, their effect on static recrystallization has been studied through a coupling of crystal plasticity with phase field [191]. Phase field modelling for nickel super alloys revealed that under high values of shear stresses, their movement speed is constant, whereas for lower values, it was found to be periodic [192].

Wang et al. [190] calculated the dislocation density for cold-rolled aluminium alloys through the Taylor model and compared their results to experimental measurements via Transmission electron microscopy (TEM), with a deviation of 30%. Significant change, around 500%, in their population as the forming progresses is also observed.

Table 2. Indicative required equipment for the measurement of the values of parameters used by physics-based CP models.

	Tensile Test	Nanoindentation	XRD	Elec. Cond.	EDS	Diffusion Couple	TEM/High Resolution	Atom Probe Tom.	Atomic Force Mi.	Fitting	Atomic-Level Simulation
Elastic constants		X									X
Dislocation core radius							X				X
Atomic volume			X								X
Attempt frequency				X							X
Mean free P path							X			X	X
Width (double kink)			X				X		X		X
Drag coefficient											X
Solid solution activation energy						X					X
Atomic concentration					X
Transmissivity parameters										X
Energy for dislocation climb	X										X
Peierls stress		X									X
Minimum dislocation density of the material							X				X
Starting dislocation density							X				X
Fitting parameters										X
Energy required for a solute atom to move								X			X
Contribution of edge dislocations to the multiplication of dislocations										X	X
Dislocation interaction			X				X		X		X
Poisson’s ratio	X										X

summarizes some of the experimental methods and techniques used for calculating specific parameters, based on the earlier categorization. Parameters that lack a physical meaning, usually non-dimensional, such as the energy barrier profile constants, can be determined through macroscopic (e.g., comparing the texture after the applied deformation and stress—strain curves) or microscopic (e.g., using Digital Image Correlation, DIC, to compare the stresses developed) tests until the result agrees with experimental results. Parameters referring to lower scales, such as kink width, are more difficult to measure experimentally and require access to equipment such as APT. Most of the required parameters can also be calculated from atomistic simulations such as ab initio or molecular dynamics. An example of parameter calculations through the ab initio approach is [193].

Experimentally determined parameters necessary for physics-based CP simulation execution are difficult to achieve. Alternative possible approaches to calibrate the models are as follows: (i) an iterative one, either via hand or via a dedicated algorithm until the results demonstrate good fit with experimental validation points such as the stress–strain curve, (ii) use machine learning algorithms to approximate the values of the required parameters, and (iii) finding relevant data in the available literature. The manual–iterative approach is extremely time-consuming, especially for formulations that have a significant number of required inputs. The automated–iterative approach seems to provide a good balance of cost–benefit, yet the, limited, available codes may be designed for a different formulation and, thus, require some basic programming skills. Arguably, the faster approach is the machine learning, albiet it requires significant amounts of data. Lastly, a dedicated literature review requires a significant investment of time in order to properly retrieve valuable data.

Table 3. Indicative references per alloy category.

Alloy Category	Indicative References
Fe and Steel	[13,36,139, 179194]
Ti	[21,114]
Ni	[24,56]
Cu	[24,25,63,87,195,196,197,198,199]
Al	[6,24,30,114,128,140,187,200,201,202,203,204,205]

shows indicative articles providing the values used for the various parameters organized by material. The parameters included in these articles are also compared to each other to highlight the variance they may exhibit for the same material.

Table 4. Indicative crystal plasticity parameters for Fe (BCC—BCT) and β-Ti (BCC).

Parameter		α-Fe	α′-Fe		β-Ti
First elastic stiffness constant with normal strain	C11	233.3	417.4	[13,195]	160	[114]	GPa
				[36]	135	[21]
Second elastic stiffness constant with normal strain	C12	135.5	242.4	[13,179]	87	[114]	GPa
		235.5		[36]	113	[21]
First elastic stiffness constant with shear strain	C44	128.0	211.1	[13,179]	54	[114]	GPa
				[36]	54.9	[21]
Shear strain rate	${\dot{\mathsf{\gamma }}}_0$	0.001	0.001	[13]	0.001	[114]	s−1
		0.56		[36]
Initial shear resistance on [111]	S0 [111]	95	406	[13,179]			MPa
				[36]
Saturation shear resistance on [111]	S∞ [111]	222	873	[13,179]			MPa
				[36]
Initial shear resistance on [112]	S0 [112]	96	457	[13,179]			MPa
				[36]
Saturation shear resistance on [112]	S∞ [112]	412	971	[13,179]			MPa
				[36]
Slip hardening parameter/self-hardening coefficient	h0	1.0	563	[13,179]	0.1		GPa
				[36]
Interaction hardening parameter/hardening matrix	hα,β	1.0	1.0	[13,179]			-
				[36]
Stress exponent	N	20	20	[13,179]			-
		3		[36]
m = strain rate sensitivity component	1/m				20	[114]	-
Curve fitting parameter	W	2.0	2.0	[13,179]			-
				[36]
Critical resolved shear stress	${\mathsf{\tau }}_{\mathrm{c}}^{\mathsf{\beta }}$
Initial slip hardness	τ0α				60	[114]	MPa
					120 × 103	[21]
Saturation value of slip resistance	τs				450	[114]	MPa
Hardening exponent	A				2.25	[114]	-
(Self and coplanar slip systems)	qab						-
(Non-coplanar slip systems)	qab				1.4	[114]	-
Burgers vector	B				2.86 × 10−10	[21]	m
Reference dislocation velocity	v0				10−3	[21]	m s−1
Exponent	P				0.71	[21]
	Q				1.1	[21]
Total initial dislocation density	ρ0				10+13	[21]	m−2
Critical radius for edge annihilation	Re				11.5	[21]	nm
Critical radius for screw annihilation	Rs				58	[21]	nm
Strength interaction coefficients	g0				0.5	[21]	-
	g1				0.5	[21]	-
	g2				0.8	[21]	-
	g3				0.8	[21]	-
	g4				0.8	[21]	-
	g5				0.8	[21]	-
	g6				0.8	[21]	-
Segment length interaction coefficients	h0				0.0	[21]	-
	h1				0.0	[21]	-
	h2				0.2	[21]	-
	h3				0.02	[21]	-
	h4				0.01	[21]	-
	h5				0.18	[21]	-
	h6				0.02	[21]	-
Free energy of activation	F0				3.1 × 10−23	[21]	J K−1

provides indicative literature for the BCC- Fe, BCT-Fe (α and α′ respectively) and the BCC-Ti (β), whereas the

Table 5. Indicative parameters for aluminium alloys.

Description	Symbol	Values	Reference	Units
First elastic stiffness constant with normal strain	C11	106.75	[114,128,140]	GPa
		100	[127]
		108	[24,203]
		108.2	[204]
Second elastic stiffness constant with normal strain	C12	60.41	[114,128,140]	GPa
		60	[127]
		61.3	[24,203,204]
First elastic stiffness constant with shear strain	C44	28.34	[114,128,140]	GPa
		30	[127]
		28.0	[24]
		28.5	[203,204]
Isotropic shear modulus	Μ	26.27	[140]	GPa
		25.0	[24,203]
Poisson ratio	ν	0.345	[140]	-
Burger vector	b	0.286	[24,140,200,205]	nm
Atomic volume	Ω	0.017	[140]	nm3
		1.7 × 10−29	[200]
Shear strain rate	${\dot{\mathsf{\gamma }}}_0$	0.001	[19,114,127,128,204]	/s
		400	[205]
Slip hardening parameter/self-hardening coefficient	h0	75	[114]	ΜPa
		400	[19]
		190	[204]
		80	[127]
Stress exponent	N			-
m = strain rate sensitivity component	1/m	20	[114,127]	-
		60	[204]
		25	[19]
		20	[187]
		6.66
		4
		333.33	[187]
		5.88
Curve fitting parameter	W			-
Slip resistance	τ0	30	[127]	MPa
		58.5	[205]
		47	[204]
		31	[19]
Critical resolved shear stress	${\mathsf{\tau }}_{\mathrm{c}}^{\mathsf{\beta }}$	31	[128]	MPa
Saturation value of slip resistance	τs	63	[19,114,128]	MPa
		95	[204]
		60	[127]
Hardening exponent	a	2.25	[19,114,128]	-
		2	[127]
(Self and coplanar slip systems)	qab	1	[128]	-
(Non-coplanar slip systems)	qab	1.4	[128]	-
Latent/self-hardening ratio		1.4	[204]
Minimum edge dipole separation	${\check{\mathrm{d}}}_{\mathrm{e}}$	1.6	[140]	nm
		1 × 10−9	[200]	m
Minimum screw dipole separation	${\check{\mathrm{d}}}_{\mathrm{s}}$	10	[140]	nm
		1 × 10−9	[200]	m
Dislocation multiplication constant	λ0	60	[140]
		100	[200]
Edge contribution to multiplication	k1	0.1	[140]
Initial overall dislocation density	ρ0	6 × 1010	[140]	m−2
Self-diffusivity (at T = 300 K)	DSD	7 × 10−29	[140]	m2 s−1
Solid solution activation energy	QSol	1.25	[140]	eV
Activation energy for dislocation climb	Qcl	3 × 10−19	[200]
Solid solution concentration	cat	1.5 × 10−6	[140]
Solid solution size	dobst	0.572	[140]	nm
Peierls stress	τP	0.1	[140]	MPa
Double kink width	wk	2.86	[140]	nm
Energy barrier profile constants	p	1	[140]
		0.233	[205]
	q	2	[205]
		1	[140]
Attack frequency	να	50	[140]	GHz
Dislocation viscosity	η	0.01	[140]	Pa s
Edge jog formation factor	k3	1	[140]

, provides similar information for FCC-Al (α).

5.3. Macroscopic Loading Conditions

Information about the loading conditions, the deformation (F), the applied stresses (P), and their rates ($\dot{\mathrm{F}}$, $\dot{\mathrm{P}}$) is also required. An example of such input is the following, where the macroscopic load is applied steadily for the duration of the experiment:

$$\mathrm{F}=\left[\begin{array}{ccc}\mathrm{x}& \mathbb{R}& \mathbb{R}\\ \mathbb{R}& \mathrm{x}& \mathbb{R}\\ \mathbb{R}& \mathbb{R}& \mathrm{x}\end{array}\right] \mathrm{P}=\left[\begin{array}{ccc}\mathbb{R}& \mathrm{x}& \mathrm{x}\\ \mathrm{x}& \mathbb{R}& \mathrm{x}\\ \mathrm{x}& \mathrm{x}& \mathbb{R}\end{array}\right]$$

In the tables above, $\mathbb{R}$ stands for a known real number and x means that the value is undefined and to be determined by the algorithm. It should be noted that these matrices refer to the macroscopic values of strain and pressure applied to the material. For example, a tensile test with a known enforced total deformation (ε_tot) along one axis can be formulated as follows [154]:

$$\mathrm{F}=\left[\begin{array}{ccc}{\mathsf{\epsilon }}_{tot}& 0& 0\\ 0& \mathrm{x}& 0\\ 0& 0& \mathrm{x}\end{array}\right] \mathrm{P}=\left[\begin{array}{ccc}\mathrm{x}& \mathrm{x}& \mathrm{x}\\ \mathrm{x}& 0& \mathrm{x}\\ \mathrm{x}& \mathrm{x}& 0\end{array}\right]$$

Alternatively, the process can be described by the applied strain rate $(\dot{\mathsf{\epsilon }})$ through the rate of the strain gradient $(\dot{\mathrm{F}})$.

During flat rolling, the material deforms in two directions. The rolling direction (RD) and in the Normal to the rolling direction (ND). Macroscopically, it is not subjected to shearing, nor is the width of the plate affected. Thus, the proper formulation would be as follows:

$$\mathrm{F}=\left[\begin{array}{ccc}{\mathsf{\epsilon }}_{\mathrm{RD}}& 0& 0\\ 0& 0& 0\\ 0& 0& {\mathsf{\epsilon }}_{\mathrm{ND}}\end{array}\right]$$

By defining all the nine elements of the deformation matrix, all the elements of the stress matrix should be left undefined. For more complex metal forming processes, e.g., extrusion, these matrices can be calculated through FEM simulations.

6. Closing Remarks

Since its introduction, crystal plasticity has widened its spectrum of applications from FCC single crystals under tensile testing to realistic, multiphase polycrystalline materials under multiaxial loading conditions. The range of the possible applications have piqued the interest of the scientific community, as has become evident through the volume of relevant published work, as well as by the variety of models that have been developed. Yet, there is a need for further improvement of the models to reduce their computational costs, to overcome remaining limitations, to further increase their accuracy, to streamline their set up, to reduce the difficulty of their calibration, and to link them with other models.

It is evident that CP simulations can provide valuable insight into the material behaviour and the impact of micromechanics during forming processes or for fatigue damage. For each case, the selection of the proper type of CP and solver are the first steps that need to be guided by a deep understanding of the process and the assumptions each choice inevitably entails. Each formulation comes with its own set of parameters that need to be determined either experimentally or via an extensive literature review. Regarding the solver (FFT or FEM), the choice is guided by the area of interest. For complex geometries or free surfaces, the FEM solver is usually best suited, compromising, however, the execution speed. For achieving reliable results, proper dependency analysis for the RVE size and mesh density also need to be performed. The selection of the constitutive model,, i.e., phenomenological or physics-based, relies on (i) the area of interest, since for macroscopical phenomena the phenomenological approach is significantly faster and delivers reliable results, whereas for local phenomena the physics-based models seems to be better suited, and (ii) the overall modelling approach, since the models incorporating parameters with physical meaning can be better integrated with various models covering all length scales.

One of the biggest struggles of implementing ICME solutions is the calibration of the models utilized. The first adopters of these methodologies did highlight the importance of creating repositories for the values of the various parameters used. This is also the case for CP. Despite the fact that increasingly more authors publish their parameters, there is a lack of a widespread repository or review articles providing this information per material and process. This should be weighed against the competitive advantage of having in-house databases for the material/processes, which as discussed are not as easy to produce. Utilizing machine learning raises the issue of proper data storage, labelling where needed, and of course it usually requires a significant amount of data. On the other hand, when possible, its use can significantly accelerate the R&D cycle.

Finally, regarding Sustainable Development, CP has a lot to offer. As a challenging field of study, researchers utilizing CP need to be highly trained (SDG 4). Its application can help optimize industrial processes (SDG 9) creating new job positions in the industry (SDG 8) and reducing the environmental footprint of the product (SDG 12 and SDG 13). CP is often provided as a service from research centers to the industry (SDG 17) and its exact impact, qualitative or quantitative, must be calculated for each case individually.