Unambiguous Identification of Crystal Plasticity Parameters from Spherical Indentation

Abstract

Identification of elastic and plastic properties of materials from indentation tests received considerable attention in the open literature. However, unambiguous and automatic determination of parameters in the case of the crystal plasticity (CP) model is still an unsolved problem. In this paper, we investigate the possibility to unambiguously identify the CP parameters from spherical indentation tests using finite element method simulations combined with evolutionary algorithm (EA). To this aim, we check the efficiency and accuracy of EA while fitting either load–penetration curves, surface topographies, or both at the same time. By fitting the results against simulation data with known parameters, we can verify the accuracy of each parameter independently. We conclude that the best option is to fit both load–penetration curve and surface topography at the same time. To understand why a given fitting scheme leads to correct values for some parameters and incorrect values for others, a sensitivity analysis was performed.

1. Introduction

The optimal approach to establish the properties of metallic materials is to apply conventional mechanical testing such as tension, compression or torsion, possibly in multiple directions if the presence of anisotropy is expected; however, such an approach is often unavailable due to limited specimen dimensions. In the case when only a thin layer of the material is available, one commonly used solution is to use micro- or nanoindentation. Such a situation can occur, e.g., in the case of coatings [1,2] or when only a thin layer next to a surface of the bulk material was obtained by some surface treatment, e.g., surface heat treatment [3], carburizing or nitriding [4] or ion implantation [5,6]. This last case is often studied as a convenient way to mimic the microstructure and property changes introduced by neutron irradiation in fission and fusion nuclear reactors.

The strain and stress state in the indented material is complex and non-uniform; therefore, a convenient way to simulate this problem comes through the usage of the finite element method (FEM) [7,8], which makes it possible to apply complicated boundary conditions, including detailed treatment of contact [9], which is especially important in indentation simulations. In addition, the introduction of nonlinear material models including plasticity and damage is relatively straightforward, especially when automatic code generation software, such as AceGen [10,11] or MFRONT [12], are used.

In the case of metallic materials, at the micro- and nanoindentation scale, the appropriate microstructure level is the level of single crystals or grains of the polycrystal. This is why the coupling of the FEM with the crystal plasticity (CP) theory is an obvious approach. So far, such simulations have been reported in numerous contributions, see, e.g., [13,14,15] for an overview.

In general, obtaining the correct set of material model parameters in CP is more challenging than in macroscopic plasticity (e.g., J2 plasticity), even when conventional mechanical testing is applied. Therefore, a number of automatic optimization strategies were developed so far. Gradient optimization was used in [16] to determine the parameters of Mg alloy. The Newton–Raphson algorithm was used to identify the parameters of AISI 316LN austenitic stainless steel [17]. The Levenberg–Marquardt method was used to optimize the parameters of Inconel 718 in [18]. Bayesian optimization was applied to obtain parameters of high strength steel 50CrMo4 in [19]. Particle swarm optimization was used to determine the CP parameters of polycrystalline Cu in [20]. Evolutionary algorithms (EAs) were used to optimize the parameters of polycrystalline materials possessing face centered cubic (FCC) lattice (copper: [21,22]), body-centered cubic (BCC) (interstitial-free steel [21], martensitic steel [23], Nb and dual-phase steel [24]), hexagonal close-packed (HCP) lattice (Zr alloy [25], Ti (pure and alloy) [24,26,27], Mg alloy [21,28,29], and Zn (pure and alloy) [30,31]). As can be seen, the EA method is the most commonly used one for obtaining the CP parameters of polycrystalline metals and alloys. Concerning the data used for fitting, the stress–strain curves were typically used [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. However, crystallographic texture [24,26,31], digital image correlation (DIC) data [17] and other experimental data [24,25] were also used for fitting in some research.

Obtaining the optimized macroscopic plasticity parameters through the usage of indentation results has also been subjected to considerable attention, both in the case of pyramidal and spherical indenters. For example, in [3], the macroscopic plasticity model parameters of Z38CDV5 steel were identified based on load–displacement curves coming from Vickers microindentation tests. A similar procedure for spherical indentation was reported in [32]. The procedure was used to determine the yield stress and strain hardening exponent of 100C6 steel. Fitting one-parameter analytical expression for strain hardening based on spherical indentation results for carburized and nitraded steel was reported in [4]. The Kalman filter method was used in [33] to determine two parameters relating compositional variation and stress–strain transfer in elasto-plastic functionally graded layer. The procedure used only load–displacement data as an input and using two values for the tip’s radius were necessary. Several methods to directly infer the parameters of power-law plasticity models were developed by Kucharski and Mróz, cf. e.g., [34,35] and references therein. A method using spherical indentation to determine the yield strength and hardening parameters for various models was presented in [36]. On the other hand, the methodology to obtain elastic and plastic parameters from sharp indentation was presented in [37].

So far, the determination of crystal plasticity model parameters from indentation was also a subject of considerable attention. In [38], the CP parameters were obtained by fitting the experimental load–displacement curves using a trial and error approach. In [39], the analytic expressions to determine the hardening exponent from spherical indentation based on three indicators (slope of indentation curve in logarithmic scale, contact area, and ring-based pile-up/sink-in volume) were developed.

Chakraborty and Eisenlohr [40] studied in detail the influence of the choice of the objective function on the possibility to obtain the correct set of crystal plasticity parameters in the CPFEM study of conospherical indentation. The Nelder–Mead (NM) simplex algorithm was chosen as an optimization algorithm. A rate-dependent power-law CP model was used. The paper reports using fitting of either the (1) load–displacement (LD) curve, (2) surface topography (ST), or (3) both as the objective function. Three parameters were optimized. It was concluded that the best results can be obtained using both load–penetration and surface topography, while using only surface topography yields the worst results. Fitting two crystallographic orientations at the same time did not increase the accuracy of the optimized parameters and lead to increasing the computational cost. Engels, Vajragupta, and Hartmaier [41] also used the NM algorithm to optimize the CP parameters from indentation; however, they fitted the topography of the deformed surface only along the lines going through pile-up maxima rather than considering the entire surface as in [40]. The optimization task was also more challenging due to the more complicated material microstructure (lath martensite).

The trust region reflective algorithm was applied in [42] in order to optimize the parameters of ARMCO iron subjected to nanoindentation. The experimental load–penetration and surface topography data were used to find four parameters of the nonlocal crystal plasticity model; however, only one surface profile (one cross-section) was used and, due to large differences between experimental and measured response, it is hard to say whether the pile-up or sink-in pattern reflecting the crystal symmetry was correctly reproduced. It is also questionable to calibrate the parameters accounting for size-independent and size-dependent effects at the same time. Logically, the more sound approach would be to find the parameters of the conventional crystal plasticity first (using the experimental data in the size-independent regime) and then calibrate the parameter accounting for size-dependence by choosing the smaller tip’s radius or shallower depth. Finally, it was shown that the initial set of parameters affects the values of the calibrated parameters considerably. The parameter ambiguity study (checking if the same set of parameters is obtained when fitting the simulation to simulation data with known parameters) was also not performed.

In this paper, we performed an analysis of the possibility to unambiguously determine the hardening parameters of the standard crystal plasticity model by using evolutionary algorithm coupled with finite element simulations. For this, both load–displacement data and surface topographies are taken into account. In Section 2, the CPFEM formulation and EA details are described. The Section 3 presents results of the performed optimizations and discussion. The Section 4 contains conclusions.

2. Methods

2.1. Crystal Plasticity Finite Element Method

The classical crystal plasticity model is applied. The particular implementation of the rate-dependent crystal plasticity in the total Lagrangian framework was already described in Section 2.1 of [43] and in Section 3.1 of [44]. Therefore, here we only recapitulate the hardening model as this is crucial for understanding which parameters are subjected to calibration.

The metallic material possessing face centered cubic (FCC) structure and deforming on 12 $\left\{111\right\}〈110〉$ slip systems is considered. The rate of shearing ${\dot{\gamma }}^r$ on a given slip system r depends on the ratio of the resolved shear stress and the critical resolved shear stress (CRSS). The CRSS value evolves with the amount of accumulated plastic shear $\Gamma =\int {\sum }_r|{\dot{\gamma }}^r|dt$ according to the Voce-type hardening law:

$${\dot{{\tau }_c}}^r=H(\Gamma )\sum _{s=1}^M{h}_{rs}\left|{\dot{\gamma }}^s\right|,$$

(1)

where:

$$H(\Gamma )=\frac{d\tau (\Gamma )}{d\Gamma },$$

(2)

and

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

(3)

The parameters that are subjected to optimization in the present paper are the initial CRSS ${\tau }_0$, the initial hardening rate ${\theta }_0$, asymptotic hardening rate ${\theta }_1$, and the back extrapolated CRSS (minus initial CRSS) ${\tau }_1$. The latent hardening parameters $h_{rs}$ are fixed and equal to 1 for self-hardening and latent hardening on coplanar systems, while for latent hardening on non-coplanar systems they are equal to 1.4.

The code of the finite element containing the CP model formulation was implemented using the AceGen software. The AceGen package [10,11] takes advantage of automatic code generation (from symbolic input in Wolfram Mathematica), automatic differentiation and expression optimization. The FEM simulations were performed using the software AceFEM. The mesh used for EA optimization had to be sufficiently refined so as to smoothly reproduce the load–displacement curves and surface topography. At the same time, the number of elements should not be too large so that multiple simulations can be conducted in a reasonable time. The mesh used in the simulations of indentation is shown in Figure 1. It consists of 7200 linear hexahedral elements, 576 linear quadrilateral contact elements, and 192 quadratic quadrilateral elements accounting for the hanging nodes formulation.

The following boundary conditions were prescribed. The horizontal displacements of the nodes belonging to the lateral surfaces of the domain were forced to be zero. The vertical displacements of the nodes belonging to the bottom surface were prescribed to increase until reaching the value of $h_{max}$ = 200 μm while the position of the indenter was fixed. Obviously, this was fully equivalent to fixing the position of the bottom plane and moving the indenter down; however, the former approach as applied here was easier to implement. Similarly as in [44] the formulation of frictionless contact with spherical tip is used. The tip is treated as a rigid body and described analytically. The contact is enforced using contact elements. Moreover, a 4-noded quadrilateral contact element with Lobatto quadrature and one additional degree of freedom per node was implemented using augmented Lagrangian multipliers as described in [45]. The only difference with respect to [44] is that here we are using tip radius equal to 200 μm whereas it was 5 nm in op. cit. However, this does not affect the results as we used the same ratio of element size to radius, and the crystal plasticity model is size-independent. The hanging nodes formulation makes it possible to join elements of different size by allowing to place a node of the smaller element in places where the larger element does not have any node. The details of the hanging nodes implementation can be found in Section 2.4 of [15].

2.2. Optimization

In order to find the parameters of the CP model, the evolutionary algorithm [25] was applied. The implementation of EA was written in Wolfram Mathematica but it is structured in the same way as the previous Python implementation of EA that was applied to calibrate parameters of VPSC [28,29,31] and SEVPSC [22] mean-field models. Nevertheless, the EA is described here again for the sake of paper’s self-containment.

The principle of EA’s operation is as follows, cf. Figure 2. First, the hardening parameters of every individual (within the population of $N_{ind}$ individuals) are randomly generated within the specified ranges. This way, the first generation can be run. The fitness evaluation step contains the most computationally intensive part, which is running the CPFEM simulations of spherical indentation for each individual. Having the results of each individual, it is possible to asses the fitness of its parameters. This is going to be explained in detail later. After the fitness evaluation is completed, the highest-ranked individuals (the ones having lowest difference with respect to the goal) are selected. Next step is the stop criterion. Various choices can be made at this stage, e.g., mean standard deviation of the fitness among the individuals being lower than some prescribed value, cf. [31] or the saturation of the increase in fitness in subsequent generations. However, in this paper we have applied the simplest possible stop criterion, which is reaching the prescribed number of generations $N_{gen}$. If the stop criterion is met, the operation of EA is terminated. Otherwise, the parameters of the next generation are produced. Two steps are used to this aim. The first one is the crossover. It consists in randomly selecting the pairs of individuals from the set of highest-ranked individuals and averaging their parameters. The second is mutation and consists of changing some of the parameters of the next generation individuals. Specifically, the loop goes over all parameters and for each generates the random real in the range from 0 to 1. If this number is lower than the mutation probability $P_{mut}$, the given parameter is again randomly generated from the range specified in the beginning. To sum up, the hyper-parameters of EA are the number of generations $N_{gen}$, number of individuals $N_{ind}$, number of highest-ranked individuals (allowed to produce offspring) $N_{hri}$, and the probability of mutation $P_{mut}$. Here, $N_{gen}=10$, $N_{ind}=40$, $N_{hri}=4$, and $P_{mut}=0.15$. The stop criterion used here was only attaining the required $N_{gen}$, but the algorithm’s performance was interactively monitored by the user.

What is important for the present contribution is the fitness evaluation step, cf. Figure 2. To the best of the authors’ knowledge, the fitness evaluation accounting for load–penetration curves and surface topographies from indentation was not used in any EA applications (although similar fitness evaluations were already used together with other algorithms, as mentioned in the introduction). The LD curves are compared against the reference one by calculating the difference between load value in subsequent penetration increments. On the other hand, the surface topography difference is calculated by creating a regular mesh in x and y directions (Figure 3a) and finding the surface position (z coordinate) by interpolation. Then, the difference between the topography of a given individual and the reference topography is just the sum of differences of z values in every point of the mesh (Figure 3b). Note that mesh here does not refer to the FEM mesh. Note also that Figure 3 is only illustrative—the true regular grid used $401\times 401$ points. The highest-ranked individuals selected for reproduction in the subsequent generation are those with lowest differences with the reference.

3. Results and Discussion

3.1. Parameter Optimization

In order to demonstrate that the evolutionary algorithm is able to find a correct set of crystal plasticity parameters and in order to check if the determined set is unambiguous, the EA was used to fit the load–penetration curve and surface topography against data coming from simulation. Obviously, the same material model was used in the reference simulation and the hardening model parameters were known, see

Table 1. The parameters of the reference simulation (pure Nickel, cf. [44,46]).

${\mathit{\tau }}_0$ [MPa]	${\mathit{\theta }}_0$ [MPa]	${\mathit{\tau }}_1$ [MPa]	${\mathit{\theta }}_1$ [MPa]
8	240	142	7.5

. The values of the parameters are intended to reproduce the behavior of pure Ni and were adopted from the literature [44,46]. The simulations were run for 001 orientation only. The ranges of the possible values for parameters are presented in

Table 2. The ranges for the parameters of the EA optimization.

${\mathit{\tau }}_0$ [MPa]	${\mathit{\theta }}_0$ [MPa]	${\mathit{\tau }}_1$ [MPa]	${\mathit{\theta }}_1$ [MPa]
5–15	100–400	50–250	5–10

.

The following EA runs were performed:

• Case 1: fitting only LD curve,
• Case 2: fitting only surface topography,
• Case 3: fitting LD curve and surface topography.

In order to check the reproducibility of the results, two different initial populations (Population 1 and Population 2) were examined. Thus, each case was run for two initial populations, thus resulting in six EA runs in total. Figure 4 shows the normalized difference (difference divided by the difference in the first generation) for each EA run. All the simulations converged maximally after generation 6, with the exception of Case 3 in Population 1, where significant difference drop appeared in the 10th generation. However, it was checked that for that particular case additional five generations (results not shown) did not result in any further improvement.

Figure 5 shows the LD curves obtained in each EA run compared against the reference simulation. In the cases where the LD curve was used for parameter fitting (Case 1 and Case 3), the LD curve was fitted very well. On the other hand, using only surface topography information (Case 2) did not result in correct fit: the level of force was either too low (Population 1) or too high (Population 2). This shows that surface information alone is not sufficient to correctly reproduce LD curve. This highlights the fact that various choices of initial parameters can lead to very similar surface topographies.

Figure 6 shows the surface topography in the reference solution and in each EA run. At first glance, it may seem that every map is very similar. Closer examination reveals that Case 1 (i.e., fitting only LD curves) results in somewhat different surface topography. In order to see it better, another figure was prepared. Figure 7 shows the maps of point-wise differences between the solution obtained in each EA run vs. the reference simulation. In addition,

Table 3. The minimum, maximum, and mean values (in μm) of the differences between surface topographies obtained in each of the EA run and the reference solution.

	Minimum			Maximum			Mean
	Case 1	Case 2	Case 3	Case 1	Case 2	Case 3	Case 1	Case 2	Case 3
Population 1	−7.275	−0.101	−0.049	0.671	0.073	0.074	−0.700	−0.044	0.008
Population 2	−4.466	−0.154	−0.023	1.042	0.755	0.185	−0.126	0.154	0.017

shows the minimum, maximum, and mean values of the differences. The differences are now clearly visible. The largest differences are present in the Case 1 (fitting only LD). This shows that fitting LD curve is not sufficient to achieve excellent surface topography fit, which is a counterbalance to the previously observed fact that fitting only surface topography is not sufficient to obtain the correct LD curve.

The particular advantage of calibrating the parameters against the simulation data is that the parameters of the reference simulation are known. Therefore, one can check whether it is possible to obtain the correct set of parameters, cf. Figure 8. Regarding the initial CRSS (${\tau }_0$), initial hardening rate ${\theta }_0$ and saturated CRSS ${\tau }_1$, the value closest to the reference was given by Case 3 regardless of population. Second best values were found by Case 1, while Case 2 gave the worst answer. Concerning ${\theta }_1$, a good approximation was obtained with Case 2 and Case 3 when starting from Population 1, but when starting from Population 2, none of the EA runs resulted in correct value. Moreover, in the latter case, the values filled the investigated range almost completely. In order to understand the efficiency of various EA simulations in determination of correct parameters, a sensitivity analysis was performed.

3.2. Sensitivity Analysis

When performing a sensitivity analysis by changing a given parameter in the range specified for EA optimization, other hardening parameters are kept fixed and equal to the reference values (see Table 1). Figure 9 shows the sensitivity analysis for LD curves. The highest variability can be seen for ${\theta }_0$ and ${\tau }_1$. ${\tau }_0$ shows moderate variability, while ${\theta }_1$ shows almost no variability in the specified range. Figure 10 presents the sensitivity analysis for surface topographies in terms of differences with respect to the reference simulation. The raw surface topographies are shown in Supplementary Figure S1. The highest variation is present for ${\tau }_1$, moderate variations are shown for ${\tau }_0$ and ${\theta }_0$, and ${\theta }_1$ presents the smallest variation.

It is well accepted in the literature (cf., e.g., [47]), that strain-hardened materials (with little potential for further strain hardening) tend to produce pile-up, while materials undergoing considerable strain-hardening develop sink-in. Generalizing, it can be stated that plastic zone size is correlated with the reserves for strain-hardening and the pile-up height is inversely correlated with it. The origin of such a behavior is that in strain-hardened materials the movement of dislocations is restricted due to high density of obstacles so the plastic zone size is small and pile-ups develop. On the other hand, in the material that was not strain-hardened the dislocations can move deeper into the material. Thus, the size of the plastic zone increases and pile-ups tend to be small or absent. Figure S1 shows that the results of CP simulations are consistent with this general understanding. Figure S1a presents the variability of ST to ${\tau }_0$. With low value of ${\tau }_0$ (high strain-hardening potential), pile-ups and sink-ins are extended in plane. With high ${\tau }_0$ (strain-hardened case), pile-ups are sharp and close to the indent. ${\theta }_0$ (Figure S1b) has a similar influence as ${\tau }_0$. High value of saturation stress ${\tau }_1$ increases the strain-hardening reserves so the influence of this parameter is opposite to the influence of ${\tau }_0$ (Figure S1c). Finally, ${\theta }_1$ has very little influence on the results in the investigated range (Figure 10d and Figure S1d).

Let us now try to understand the efficiency of each EA in determining the proper parameters. Case 1 that considered only LD curves, was able to find a close guess for ${\tau }_0$, ${\theta }_0$ and ${\tau }_1$. Its prediction of ${\theta }_1$ was however the worst. The correct determination of ${\tau }_0$, ${\theta }_0$, and ${\tau }_1$ is understandable in view of Figure 9a–c. On the other hand, it is not surprising that the value of ${\theta }_1$ was not determined correctly (remembering that LD curve was insensitive to this parameter alone). It was already outlined that Case 2 (considering only surface topography) provided the worst parameter calibration. Its prediction of ${\theta }_1$ is however better than the one provided by Case 1. Although the surface topography was relatively insensitive to the value of this parameter (cf. Figure 10d), still the variability of surface topography with respect to ${\theta }_1$ was greater than the variability of LD curve in that case. Thus, it seems that this variability enabled Case 2, to perform better in finding ${\theta }_1$, in particular when starting from Population 1. Finally, Case 3 was able to find reasonably good values of ${\tau }_0$, ${\theta }_0$, and ${\tau }_1$ regardless of the initial population, which is understandable considering that both LD curves and surface topographies were sensitive to these parameters. As LD curves are insensitive to ${\theta }_1$ and surface topographies are almost insensitive to this parameter, it seems also reasonable that determination of its value when using both information for fitting heavily depended on the initial population.

In order to check if choosing a wider range of asymptotic hardening rate leads to higher sensitivity, we conducted another sensitivity analysis with ${\theta }_1$ ranging from 1 to 20 MPa. The results are shown in Figure 11. It appears that the LD curves present a higher variability but it is still small as compared to ${\theta }_0$ and ${\tau }_1$. On the other hand, the variability of surface topography does not increase significantly even with the wider parameter range. It can be thus concluded that finding precise value of ${\theta }_1$ using instrumented indentation data is a challenging task and one should rather focus on the determination of the other three parameters while keeping the value of the asymptotic hardening rate fixed.

3.3. Discussion

When analyzing the results obtained in this study, it is easy to see that the results of the EA optimization greatly depend on the (here random) choice of the initial population. This should give the hint when seeking for the parameters based on experimental data. Namely, one should ideally perform several EA runs with different initial populations and pick the set of parameters that gives both the LD curve and surface topography in the closest agreement with experimental data. Another option would be to enlarge the initial population, however, giving a hint on the number of required individuals in the initial population is beyond the scope of the present paper.

All the performed optimization were performed for one crystallographic orientation (001) only. One may wonder if changing the orientation could affect the effectiveness or accuracy of the developed EA approach. To our understanding, this is not the case, although we did not run any additional optimizations to prove this. In the slips-system-based CP model, the symmetry of the pile-ups and sink-ins in plane is dictated by the crystallographic orientation. On the other hand, the corresponding heights and depths are affected by hardening parameters. Therefore, it seems that changing the orientation or even adding additional orientations (the latter case impacts the computational cost proportionally to the number of orientations) will not change the results considerably. In fact, our understanding of this issue is consistent with the observations of [40]. The authors studied the influence of using two orientations simultaneously in their optimization. They concluded that although the computational cost was doubled, there was no noticeable improvement in the performance of their optimization algorithm.

All the optimizations presented in the article were performed using the frictionless contact with rigid indenter. One may argue that friction is important and should not be neglected. Although modeling of friction is a broad and interesting subject by itself, here we restrict our attention to the simplest case of Coulomb friction that was considered in majority of the related literature. In [38,48], the authors concluded that in spherical indentation friction has a negligible effect on load–displacement curves but noticeably affects the pile-up heights. We have performed our own study of the effect of friction. To this aim, the implementation of the contact element with additional degrees of freedom responsible for tangential contact was used. Figure 12 shows the effect of friction on LD curves. Indeed, consistently with op. cit., the curves are almost the same until about $\frac{h}{R}=0.4$. However, they differ considerably at a later stage. As in [38], the ratio $\frac{h_{max}}{R}$ was about 0.1, it is logical that their results did not show the difference that can be observed at extremely large penetration.

Figure 13 shows the changes of deformed surface maps with increasing friction coefficient. The higher is the friction, the broader is the zone affected by plastic deformation, which is consistent with a simple mechanistic understanding of the problem and the analysis presented in [38]. Although the code of the contact element that includes friction is more complicated, including friction does not pose a major computational problem (it was verified that a time of a single simulation does not increase by more than 10%). Friction could be thus included in two different ways. In the first one, friction could be considered as an additional parameter that is being optimized. Such an approach would make the studied problem more ambiguous. In the second, one could use a fixed friction coefficient; however, as the value of the friction coefficient is not known, one would have to take some arbitrary value (e.g., 0.3 was used in [40]). Finally, in the studied example of fitting against the reference simulation with known parameters taking arbitrary value of friction coefficient would lead to exactly the same answer and all the conclusions would stay unchanged; thus, we claim that restricting our attention to the frictionless case is justified in the present paper; however, in the future studies devoted to obtaining parameters based on real experimental data, the effect of friction will have to be carefully considered.

Obviously, the CP model considered in this paper is one of the simplest. In order to take into account some additional effects, typically models also contain additional parameters. In the recent years, some attention was paid to modeling of indentation size effect (ISE). In the case of spherical indentation ISE manifests itself by an increase in hardness with a decrease in indenter’s radius (with a fixed depth-to-radius ratio). In general, two types of approaches were presented in order to take ISE into account. In the first one, the size effect is assumed to result from the incompatibility of accumulated plastic deformation. Such an approach was used to simulate indentation in the presence of ISE in [49,50,51,52,53]. In the second approach, the ISE stems not only from the incompatibility of accumulated plastic deformation, but also from the history of incremental plastic incompatibility. Simulations of indentation following this point of view were presented in [46,54,55]. In this paper, we performed simulations with indenter’s radius equal to 200 μm, so we do not expect any significant ISE in this regime. There are however no fundamental restrictions in the case one would like to use the optimization method proposed here to obtain the parameters of a size-sensitive model. One should only consider two things. First is the increased number of parameters, which is always challenging in optimization; however, it seems that using information from indenters with different radii should help circumvent this problem. The other thing that should be carefully taken into account is the computational time. Note that here, the evolutionary algorithm was conducted on a standalone workstation. These particular simulations could be also performed on any modern PC with reasonable performance. The time spent on the simulations would probably then increase several times but would also stay in the same order of magnitude. As the simulation time in the case of the strain gradient CP model is typically much larger than in standard CP (due to additional degrees of freedom), the cost of performing an optimization may become prohibitive not only on standard PCs, but also using an efficient workstation. Thus, one should consider switching to a high-performance computer when dealing with a size-sensitive CP model.

As mentioned in the beginning, nanoindentation is a useful way to study the behavior of ion-implanted surface layers. Wang et al. [56] simulated indentation of irradiated pure zirconium. Nie et al. [57] applied the dislocation density-based CP model to simulate the behavior of irradiated A508-3 steel. Xiao et al. [52,58] applied the strain gradient CP model to irradiated copper and tungsten. In most of the CP models accounting for irradiation, the hardening stemming from irradiation defects is dependent on their number and size. It can be also easily omitted, e.g., by setting the number of defects to zero; therefore, the suggested solution for this application is to perform two subsequent optimizations. In the first one, the parameters of a virgin material should be established. In the second, only the additional parameters accounting for irradiation effects should be established.

The other interesting effect, which in the context of continuum modeling, was examined in one paper only [44], is the effect of chemical disorder in concentrated solid-solution alloys (CSAs). In op. cit., this effect was taken into account by merging the isotropic plasticity and CP models using a single parameter. The simulation time increases to some extent, but the increase is moderate and performing the optimization in a reasonable time on a standard machine should be possible. The more important challenge is the number of parameters: the model includes the hardening parameters of both isotropic and crystal plasticity models together with the additional parameter accounting for their proportion. The problem of numerous parameters is also inherent to hexagonal close packed (HCP) materials (such as e.g., Ti, Mg, Zn, and Zr) where there are many active slip and twinning systems and each needs its own set of hardening parameters. As we have seen, in [40] only three parameters were subject to optimization while four parameters were considered here. In order to unambiguously determine more parameters one should also increase the amount of information for fitness evaluation. Suggesting which additional information could be added is however outside the scope of the present paper.

4. Conclusions

We have performed a numerical analysis of the possibility to unambiguously determine the hardening parameters present in the crystal plasticity model using evolutionary algorithms with two types of data coming from spherical indentation. Namely, the load–displacement curves or surface topographies were fitted. We concluded, in agreement with [40] (but implementing a different approach), that using both information results in best parameter values and using only LD curves is better than using only ST.

Based on the knowledge gained in the course of this study, we plan to perform a similar analysis using experimental data. Aside from calibrating the parameters of standard CP model, the parameters responsible for some specific effects, such as irradiation or chemical disorder, will be also calibrated; this is however outside the scope of this contribution.