Proper Orthogonal Decomposition–Radial Basis Function Surrogate Model-Based Inverse Analysis for Identifying Nonlinear Burgers Model Parameters From Nanoindentation Data

Conventional mechanical testing methods, such as uniaxial, flexural, or bend tests, are widely used to characterize the mechanical behavior of different materials [1–3]. However, these tests require special specimen preparation due to extensive size or shape restrictions, and are unable to characterize localized changes in material behavior. Nanoindentation, which grew in popularity after the initial work of Doerner and Nix and Oliver and Pharr, can characterize the nanoscale mechanical behavior of a material if a relatively flat test surface with small surface roughness is provided [4,5].

Even though nanoindentation is a useful technique for materials that are a challenge to characterize using conventional tests, the highly nonlinear nature of nanoindentation load–displacement data requires analytical or numerical analysis before the data can be interpreted. For simpler material models, such as elastic or linear viscoelastic, analytical solutions can be used to identify the model parameters [5–7]. However, for complex material behavior, such as nonlinear viscoelasticity, analytical solutions are not available [8]. In these cases, inverse analysis coupled with finite element (FE) analysis can be adopted to calibrate the constitutive relationship for observed mechanical behavior [9–15]. Inverse analysis finds the set of model parameters through optimization for which the value of the objective function (difference between simulated and experimental data) is the minimum.

For material models for which no a priori information is available, a global optimization technique such as the genetic algorithm (GA) or particle swarm optimization is usually recommended. These techniques are very useful at finding the global minima in a parametric space, albeit the computational expense associated with the process is rather high. This is a challenge in FE-based inverse analysis since computationally expensive FE simulations are required at every iteration. Moreover, the computation expense of the process increases tremendously if the number of model parameters required to capture material behavior is high, which is often the case for complex material behavior such as nonlinear viscoelasticity.

One way to minimize the computational expense is to use a numerical model that approximates the input–output relationship of the FE simulation. In this process, called surrogate modeling or meta-modeling, the numerical model is trained with FE data and used inside an optimization routine instead of the actual FE simulation. A wide range of mathematical tools, such as artificial neural network [16], support vector regression [17,18], surrogate model accelerated random search [19], Kringing [20], radial basis function's (RBFs) [21], and nonuniform rational B-spline [22], have been used in the past to create surrogate models. Studies concerning the overall performance of interpolative surrogate models show that RBF is highly capable as compared with other meta-modeling techniques [23–25]. RBF-based surrogate models are capable of providing good approximations for a wide variety of functions, making them ideal for problems where definitive knowledge about the system function is absent.

Proper orthogonal decomposition's (POD) ability to reduce model size combined with RBF interpolative power has been previously used to identify material model parameters using nanoindentation data [26–28]. Although these studies have demonstrated the usefulness of the POD–RBF surrogate modeling approach in solving nanoindentation inverse problem, its applicability has been limited to simple time-independent material models. Time-dependent complex material behaviors, such as nonlinear viscoelasticity, require higher number of model parameters to describe the material response. In addition, in the case of time dependence, a mixed-mode interpolation problem arises due to the involvement of both spatial and temporal variables. This is especially challenging for RBF-based surrogate models due to the scaling difference between spatial and temporal variables.

In published literature, a lack of systematic approach to nanoindentation-based identification of nonlinear viscoelastic model parameters can be observed. The few studies that have tackled the problem considering a large number of model parameters have not reported or emphasized the parameter identification procedure [29,30]. Since the nanoindentation-related inverse problem is ill-posed in nature, using unnecessarily high number of model parameters can often lead to a local minima. Furthermore, to the best of our knowledge, none of the previous studies have demonstrated the applicability of surrogate model-based parameter calibration for a nonlinear viscoelastic constitutive relationship. Therefore, this study is designed to identify nonlinear viscoelastic model parameters using nanoindentation data through a systematic POD–RBF surrogate modeling approach, and to overcome the limitation of the RBF-based surrogate model for spatial–temporal mixed-mode interpolation problem.

To demonstrate the applicability of the POD–RBF technique for determining nonlinear viscoelastic model parameters, nanoindentation was carried out on a candidate material. The FE model was constructed using a commercial finite element package abaqus (Dassault Systémes, Providence, RI). The constitutive model was implemented in a user-defined subroutine (UMAT) via a fortran script.

Using a combination of Taguchi and analysis of variance (ANOVA) methods, the sensitivity of FE simulation output to the material model parameters was identified. This information was utilized to reduce the number of levels in which individual model parameters were varied. FE simulations were carried out for a finite number of model parameter sets where at least one parameter was varied within the parametric space. The simulation data were utilized to create the POD–RBF surrogate model. The model parameters were then identified through a GA-based optimization algorithm that minimized the difference between the experimental data and the surrogate model's prediction.

Nanoindentation experiments were conducted on an MTS Nanoindenter XP (Agilent Technologies, Santa Clara, CA) using a load-controlled scheme with a Berkovich tip. The maximum load was set to 0.5, 0.75, and 1.0 mN for the experiments. A triangular loading profile was chosen with 30, 45, 60, or 240 s total duration. The duration was kept constant for the loading and unloading segments. The acceptable thermal drift rate was chosen to be 0.15 nm/s and the experiment was started only after the system was stable and within the acceptable thermal drift rate. The raw load–displacement data were corrected for thermal drift before being used in the optimization algorithm.

An epoxy polymer, fabricated using EPON 862 and Epikure 3274 (Miller-Stephenson Chemical Company, Inc., Danbury, CT) mixed at 100:40 weight ratio, was selected for the nanoindentation experiments. Surface preparation of the samples was following standard metallographic techniques.

The commercial finite element software, abaqus, was used for modeling the nanoindentation experiment and for solving the FE problem. The Berkovich tip was modeled as a 3D discrete rigid body while the sample was modeled as a 3D deformable body. A finer mesh was provided to the sample near the contact region to ensure good convergence and also to improve the quality of the finite element solution.

Contact between the indenter and the sample was defined as surface-to-surface contact. The indenter and the sample were assigned as the master and the slave surfaces, respectively. The element type was chosen from the eight-node brick element family (C3D8). The FE problem consisted a total of 1323 elements and 1817 nodes. Figure 1 shows a schematic of the abaqus finite element model in this investigation.

A spring–dashpot model developed by Marin and Pao was used as the constitutive material model [31]. For the linear case, this model is generally called the four-parameter Burgers model and it is formed by a serial connection of a Maxwell element to a Voigt element [32]. Figure 2 shows the schematic of the nonlinear Burgers model.

where E and ν are the Young's modulus and Poisson ratio, respectively; J₂ is the second invariant of the deviatoric stress tensor s; C_s, C_t, m_s, m_t, and t_ε are the nonlinear material parameters; σ is the Cauchy stress tensor; i and j are the indices ranging across 1, 2, and 3; and δ_ij is the Kronecker delta.

A stress increment $Δσ′=CΔε′$ is then added into the system equation to account for the creep strain in Eq. (4), with C being the Jacobian stiffness matrix calculated from Eq. (3).

An important factor to consider in nanoindentation experiments is the friction generated between the indenter tip and the sample surface. This can be taken into consideration in the analysis by defining a sliding contact with a finite friction coefficient between the surfaces.

To simplify the FE model, many researchers have opted to assume frictionless contact between the tip and the sample surface [34,35]. This assumption was based on the fact that the nanoindentation response was insensitive to friction. However, a few researchers using FE simulation have shown that the influence of friction depends on other factors, such as the material model used and the geometry of the tip [36,37]. If these factors change from one study to another, studying the effect of friction becomes a necessity.

Previous studies that used the nonlinear Burgers model in a nanoindentation-based FE study did not investigate the effect of friction on the load–displacement data [29,30]. Hence, in order to improve the understanding, this study performed a parametric study, where the value of the coefficient of friction was varied in four steps ranging from 0 to 0.5.

A sensitivity study of the parameters was conducted, before generating FE simulation data, by varying the model parameters. This information helps to reduce the number of FE simulations used for training [38].

The nonlinear Burgers model has seven independent parameters, i.e., E, ν, C_s, m_s, C_t, m_t, and t_ε. In general, the nanoindentation load–displacement output is not influenced by the Poisson's ratio, ν [39–41]. Therefore, in order to keep the number of independent parameters to a minimum, ν was given a constant value of 0.34 and was not included in the sensitivity analysis scheme.

Sensitivity analysis was carried out using the ANOVA technique. The data required for ANOVA were generated using the Taguchi design of experiments method. In this study, the six nonlinear model parameters were varied in three equidistant levels. A statistical software, minitab (Minitab Inc., State College, PA), was used to design the experiments. For six parameters, where each parameter was varied in three levels, Taguchi L₂₇ orthogonal array design is appropriate. Table 1 shows the levels of the six individual parameters of the nonlinear Burgers model.

In Eq. (5), i = 1, 2, 3, …, n, and n is the number of data points in a single nanoindentation simulation or experiment.

Sensitivity of the nanoindentation output was also determined in a different way, where the difference in output between lowest and highest limit of the individual parameter levels was determined. Unlike the Taguchi–ANOVA procedure, here only the indenter depth at maximum load and the depth after unloading were studied.

This type of parametric sensitivity analysis has been used previously in understanding nanoindentation experiments. In this study, this sensitivity analysis was performed in order to complement the Taguchi–ANOVA procedure, and to get an understanding of how individual parameters contributed to the variance of an indentation load–displacement plot's two key features.

The POD theory was developed to approximate a function over some domain of interest based on the known input–output relationships [42–44]. This study followed the POD–RBF surrogate training procedure outlined by Buljak [45] and Rogers et al. [46]. The POD–RBF method requires creating snapshots (input–output relationships of the system) from which the surrogate model can be established. The more snapshots or training data points that can be utilized to generate the surrogate model the better the approximation.

Nonetheless, the computational burden associated with generating a large number of snapshots becomes the limiting factor in obtaining very high-fidelity predictions from the surrogate model. Sensitivity analysis could be utilized to reduce the number of snapshots without the loss of approximation error [38]. Hence, in this study, a similar approach was adopted to reduce the computation burden of training the surrogate model for the nonlinear Burgers model.

In this study, FE training data were generated by simulating four different experimental conditions. In these experimental conditions, the maximum load was kept constant at 1.0 mN, while the strain rate was varied from 1/30 s⁻¹ to 1/240 s⁻¹. Separate surrogate models were trained using FE data to represent each of these experimental conditions. Creating separate surrogate models ensured that each of these numerical models contained only spatial variables. This was done to overcome RBF's limitation in dealing with mixed-mode scaling differences. The approximations from each surrogate model were compared against their own experimental indenter displacement data to form the objective or error function. Based on a recent investigation by Hamim and Singh, a multiquadratic RBF with the shape parameter value of c_j = 0.5 was chosen for this study [38].

A multi-objective GA-based optimization procedure was used to identify the parameters of the nonlinear Burgers model. The procedure was implemented using matlab's (Mathworks Inc., Natick, MA) global optimization toolbox. Scores of the first and all subsequent generations were determined by evaluating the fitness function that was submitted to the program via a matlab script.

A double vector initial population of 200 was randomly created with a uniform distribution. Selection of the next generation parents was carried out via a tournament of size 2. Eighty percentage of the next generation population was produced via intermediate crossover, while the remainder was created through Gaussian mutation. The values of scale and shrink parameters were set to 1.

In case of intermediate crossover, the creation of children from two parents is controlled by a single parameter ratio. The value of this parameter was selected to be 1 for this study. Next generation children were created through a random weighted average of the parents.

The forward migration parameters, fraction and interval, were selected to be 0.2 and 20, respectively. A total number of generations for the optimization algorithm were chosen to 100× number of parameters, i.e., 100 × 6 = 600 for this study. The fitness (error) function tolerance was chosen to be 10⁻⁴.

The effect of the friction coefficient on maximum and residual depths attained during nanoindentation was analyzed. Figure 3 shows the effect for conditions represented by a fixed maximum load of 1.0 mN with different loading–unloading times. For any given value of the friction coefficient, the values of both maximum and residual depths decreased as compared to the corresponding frictionless case of indentation. The plots represent the reduction of depths between the simulations of a frictionless condition and a particular friction coefficient (e.g., f = 0.125, 0.25, or 0.5). All other conditions, e.g., boundary conditions, maximum load, loading–unloading time, and material model parameters, were kept constant. Both Figs. 3(a) and 3(b) are plotted at the same vertical scale for ease of comparison.

Figure 3 shows that, for any value of the friction coefficient within the studied range, the value of maximum and residual depths was reduced in comparison to the frictionless condition. This behavior was found to be true for all studied conditions with varying loading–unloading times. When friction is considered in a nanoindentation study, part of the energy that could be utilized to displace the material is dissipated as frictional energy. This loss of energy leads to a reduced displacement of the indenter. Similar behavior has been observed for the simulation of elastic–plastic indentation. DiCarlo et al. observed that the introduction of friction in the model increased the calculated hardness by lowering the indenter displacement at maximum load [47].

Figure 3(a) also illustrates that for a given friction coefficient, the reduction in maximum depth varied as a function of loading–unloading time. The greater the loading–unloading time, the lower the maximum depth observed in comparison to the frictionless condition.

Figure 3(b) shows the reduction in residual depth values between a frictionless simulation and a finite friction coefficient simulation. Here, for any given value of the friction coefficient, the difference diminished with increasing loading–unloading time. When lower loading–unloading time is used in an indentation experiment, the viscoelastic creep response is subdued. Hence, the elastic response has relatively higher dominance on the overall deformation behavior. The observed behavior may mean that friction has more effect on the residual depth when viscoelastic behavior has lower dominance over the nanoindentation data.

For both maximum and residual depths, the reduction in indentation depths was observed to be very small for all loading–unloading times. For instance, in the case of t = 240 s, the condition which showed the highest deviation for maximum depth, the reduction was found to be ≈ 0.5%. On the other hand, the highest reduction in residual depth was found to be ≈1.8%.

Figure 4 shows the effect of the friction coefficient on maximum and residual depths when the loading–unloading time was kept constant and the maximum load was varied from 0.5 to 1.0 mN. Similar to the results discussed earlier, where loading–unloading time was varied, the values of maximum and residual depths were found to have decreased from the values obtained for the frictionless condition.

Figure 4(a) shows that increase in the friction coefficient resulted in a higher reduction in maximum depth in comparison to the frictionless counterpart. This observation was common for all three maximum load conditions. Similar behavior was observed for residual depth reductions, as illustrated in Fig. 4(b). Since a higher friction coefficient would lead to greater frictionally dissipated energy, greater reductions in maximum depth, as compared to the frictionless conditions, would be expected. Nevertheless, the overall differences were very small. For f = 0.5, which provided the maximum differences, reductions in maximum and residual depths were found to be ≈0.5% and 1%, respectively.

Another observation common for both Figs. 4(a) and 4(b) was that between different maximum load conditions, there were little differences for a given coefficient of friction. This could mean that within the given range of loads (0.5–1.0 mN), the maximum load has no effect over the friction behavior. However, determining whether the maximum load insensitivity is a universal fact requires further investigation.

This parametric study shows that the inclusion of friction in the finite element model leads to changes in the indentation load–displacement response. Nonetheless, the variations are small for the conditions of interest. An actual nanoindentation experiment can never be entirely frictionless. Therefore, this study included the effect of friction in the model by using a coefficient f = 0.25 for all sensitivity analysis and surrogate model development purposes.

Table 2 shows the result of sensitivity analysis carried out using Taguchi-based design of experiments and ANOVA.

The “% Contribution” data, a measure of variation contributed by individual parameters, show that, except for t_ε, all other parameters contributed toward the overall variation, although majority of the contribution was originating from only two parameters, C_s and m_s.

Figure 5 shows the sensitivity to the model parameters for the indentation depths at maximum load and at the end of unloading, i.e., maximum and residual depths. Similar to the Taguchi–ANOVA-based sensitivity results, it is seen that t_ε had little to no impact on maximum or residual depths. Meanwhile, C_s and m_s had the most significant impact for both. Furthermore, the level of significance was much more pronounced for residual depth as compared to maximum depth.

Another observation that could be drawn from these results was that the elastic modulus has a positive correlation with the strain rate. In other words, when the strain rate was higher, both the maximum and residual depths were comparatively more strongly influenced by variations in the elastic modulus. One explanation of this observation is that when the strain rate is higher, viscoelastic response gets subdued due to the inherent time lag between the elastic and viscoelastic responses. As the strain rate becomes smaller and smaller, the viscoelastic (or creep) response catches up with the elastic response. Hence, the elastic part of the displacement becomes less dominant in the overall displacement pattern of the material.

Figure 5 shows another important observation, which is contrary to ANOVA results. The two parameters, C_t and m_t, show opposite trends in these two sensitivity tests. In the case of ANOVA, C_t showed substantial influence toward the output, while it was fairly insignificant in Fig. 5. Contrary to the effect of C_t, the m_t parameter exhibited a high level of sensitivity in Fig. 5 but insignificant influence in the ANOVA results. This is because Fig. 5 represents the sensitivity of individual parameters toward two points in the nanoindentation load–displacement plot, namely the maximum and residual depths. Although these two points are very important in understanding a material's response, they do not represent the entire load–displacement plot during indentation. It is possible that two plots, distinct in every other way, can have the same maximum and residual depths. That is why having multiple complementary means of determining sensitivity can provide a broader view of the problem.

The findings from the sensitivity analysis were taken into account to revise the number of levels for each nonlinear Burgers model parameter. As discussed, t_ε showed no influence over the output of the nanoindentation simulations. This implies that either t_ε cannot be accurately determined from a Berkovich nanoindentation experiment or that it is a redundant parameter in describing the material response. Hence, t_ε was given a constant value to reduce computational expense.

The two parameters that were the most influential, C_s and m_s, were varied at four levels. Meanwhile, moderately influential parameters, E and C_t, were varied at three levels. According to ANOVA, m_t did not significantly influence the output. However, the parametric study showed that m_t had some influence over the maximum and residual depths. For this reason, instead of assigning a constant value to m_t, it was varied at two levels.

Table 3 shows the corresponding levels for each parameter that was selected based on the sensitivity analysis. In a full factorial basis, a total of 3 × 4 × 4 × 3 × 2 × 1 = 288 finite element simulations were carried out in order to generate the surrogate model for every single experimental condition. For each of these simulations, 100 load–displacement data points were obtained. Since there were four individual experimental conditions to represent, a total of four surrogate models were developed. The snapshot matrix used to generate each of these surrogate model had dimensions of 100 × 288.

After the POD model reduction process was carried out and the RBF coefficients were calculated, the POD–RBF surrogate models were used to approximate FE outputs within the specified parametric space. The objective functions were written in matlab where each surrogate model's output was compared against the corresponding experimental data. These objective functions were the used within the matlab global optimization toolbox to run multi-objective GA-based global optimization. Table 4 shows the result from the global optimization algorithm.

The optimized set of parameters found using the optimization algorithm is listed in Table 4. In addition, Fig. 6 shows the comparison of predicted and experimental data for all four experimental cases. These were the experimental conditions that were closely followed in creating finite element models and were used to train the predictive or surrogate model. From Fig. 6, it can be seen that all four surrogate model outputs were very close to the corresponding experimental data. This demonstrated that the multi-objective genetic algorithm-based optimization procedure was successful in finding a common minima taking the constraints into consideration.

Although, surrogate model predictions mostly agreed with the experimental data, a few inconsistencies were observed. For example, the final unloading portion data for the loading–unloading time t = 30 s did not match very well with the experimental data. Similar behavior was observed for t = 45 s, even though qualitatively the difference between prediction and experiment had diminished. For higher loading–unloading times, i.e., t = 60 s and 240 s, the differences were noticeably much smaller.

Figure 7 compares the finite element model outputs for the identified Burgers model parameters with the corresponding nanoindentation experimental data. The surrogate model developed with finite element simulation data was not trained for these experimental conditions. Therefore, these conditions were used to validate the optimized set of Burgers model parameters.

In training the surrogate model for approximating nonlinear Burgers model output, experimental conditions with varying loading–unloading times were used. On the other hand, these validation experiments used maximum loads for which the model was not trained. This decision was deliberately made in consideration of the fact that nonlinear viscoelastic behavior depends not only on the strain rate but also on the strain levels associated with the experiment. It was assumed that changing the load levels would result in a change of strain levels.

From the validation plots for untrained experimental conditions, it could be observed that the finite element simulation output closely matched with the experimental data. In two of the cases (Figs. 7(b) and 7(d)), however, a portion of the unloading curve showed some discrepancies in a qualitative sense.

The variation between different plots for the optimized set of model parameters could stem from different factors. One such factor could be the friction coefficient used in the finite element model. As it can be seen, the effect of the friction coefficient varied depending on the experimental conditions that were being replicated. Therefore, the error associated with using a particular friction coefficient also varied from one experimental condition to another. Since the whole process of inverse analysis depended on numerical manipulations, different errors in the finite element simulation data could have skewed the parameter optimization in one way or the other.

The inverse problem of identifying nonlinear viscoelastic model parameters was solved using an FE simulation trained POD–RBF surrogate model and a multi-objective GA optimization algorithm. A constitutive material model, known as the nonlinear Burgers model, was used to represent the nonlinear viscoelastic behavior of a material under indentation loading. RBF's limitation of properly scaling temporal and spatial variables was circumvented by using multiple surrogate models representing experimental conditions with different loading–unloading times. The effect of the friction coefficient between the sample and the tip surface was investigated for the particular constitutive model. The individual model parameter's impact toward simulation output was identified using Taguchi and ANOVA method, and this information was used to reduce the number of training points.

From the sensitivity analysis, it was observed that various parameters had different impact on the load–displacement output. The nanoindentation response was most sensitive to the change in parameters, C_s and m_s, and least sensitive to t_ε. The effect of the friction coefficient on the simulation study was found to be small; the maximum difference for the studied conditions was approximately 2.0%. It was also found that multi-objective GA coupled with separate surrogate models representing different experimental conditions was capable of determining the material model parameters within the parametric space. The identified parameter set provided good agreement between numerical or simulation and experimental data.

• Office of Experimental Program to Stimulate Competitive Research (NNX14AN41A)

• Oklahoma State University (C.F. Colcord Endowed Chair)