Parameter Estimation of a Nonlinear Burgers Model Using Nanoindentation and Finite Element-Based Inverse Analysis

List of Figures

Chapter 1 Introduction

The advent of 21st century has led to the development of several new materials for different engineering applications. Metals or metallic alloys which have been predominantly used for structural applications are being replaced by lighter and stronger composite counterparts. Most often polymers are used as the matrix materials in these composites. Polymers are also being increasingly used in thin film applications. Thin films have already found considerable industrial applications, e.g. in the production of plane and automobile, as well as in electronic, optical, medical and chemical devices. Unlike metals or ceramics, which exhibit simple elastic–plastic behavior, mechanical behavior of polymeric materials are very complex. Understanding the mechanical behavior of polymers, which exhibit time–dependent responses under applied load, is an important issue in predicting the performance of these materials while in use.

Biomaterials are another important area where knowing the mechanical behavior would greatly improve the quality of understanding of these material systems. The onset of various diseases such as breast cancer [6, 7], atherosclerosis [8, 9], fibrosis [10], and glaucoma [11] has been found to be related with change in tissue compliance. The application of fast and reliable characterization of biomaterials would not only be beneficial for disease progression, but also in designing improved artificial organs, building virtual surgical simulators and automated robotic surgeon [12, 13].

Therefore, one of the most important question in today’s materials science is understanding mechanical behavior of materials under different loading conditions at different length scales. Traditional mechanical testing methods such as tension, compression, flexural, and bend tests can only provide the macroscale mechanical behavior [14, 15, 16, 17]. Macroscale test data often differs from the nanoscale test data if the material system is non-homogeneous [1]. Moreover, these testing methods require the specimen to be of certain size or shape, which is often difficult to obtain. Complex fixture design and gripping issues are a few more challenges to overcome in case of traditional testing methods.

Nanoindentation is one of the most promising material characterization technique that has the potential to overcome the complexities of conventional testing methods. Nanoindentation involves probing a material with a very small, hard diamond tip of known geometry, while the load and the displacement experienced by the tip is recorded continuously. This load and displacement data is a direct function of material’s inherent mechanical properties, and thus makes it theoretically possible to attain mechanical properties from nanoindentation data. The biggest advantage of nanoindentation, which is driving the use of this technique, is that it removes the size or shape restriction placed by the macro or bulk testing techniques.

However, indenting a material and recording loads and displacements is just the preliminary step in obtaining mechanical properties from the nanoindentation process. As loads and displacements are the only experimentally measurable variables, in order to extract mechanical properties, suitable analytical or numerical methods that relate indentation loads and displacements to material properties are required [18, 19]. This is a challenging task because unlike traditional uniaxial testing methods, nanoindentation load–displacement data comes from complex multi-axial loading, thus making it much more difficult to analyze and subsequently interpret in terms of mechanical properties.

Past developments in this area has reached to the point where nanoindentation measurements could be related to mechanical properties for materials exhibiting simple elastic or elastoplastic material behavior. However, for materials exhibiting complex material behavior, such as time–dependent material behavior, suitable analysis technique is not available. Understanding mechanical behavior of materials is a root problem, and it carries forward to severely limit applications. For example, accurate mathematical descriptions of the mechanical behavior of soft tissues remain the limiting factor in the advancement of realistic medical simulations and non-invasive diagnostic tools as soft tissues exhibit nonlinear stress-strain behavior at large deformations.

Developing an analysis technique for nanoindentation of soft materials, such as, polymers, gels, metals at high temperature, and biomaterials, is especially challenging due to the inherent time–dependent mechanical behavior [20]. Time-dependent mechanical behavior, which is known as viscoelasticity or viscoplasticity needs to be taken into consideration in order to accurately predict material behavior under service [21]. In case of a viscoelastic or viscoplastic material, the stress state not only depends on the strain, but also the strain rate.

In chapter 2, a comprehensive review of the existing nanoindentation-based analysis techniques is presented. This includes both analytical and inverse approach for analysis of load–displacement nanoindentation data. Based on the state-of-the-art review few questions are raised, the answers to which if known could significantly improve the applicability of nanoindentation technique for material property characterization.

In chapter 3, the development of a technique that can be used to characterize nonlinear viscoelastic behavior of soft materials is described. The theories and challenges of the specific techniques is also provided to improve the understanding of the effectiveness of each constituent of the overall technique.

In chapter 4, an application of the developed technique is presented for an elastoplastic material behavior. This case study is used to understand the overall numerical technique in the context of the nanoindentation experiment. By extending this understanding of the numerical technique, the problem of determining nonlinear viscoelastic constitutive model parameters is solved. Later part of chapter 5 is utilized to draw a conclusion of the study, as well as to report about the possible future works that could improve the robustness and the general applicability of this technique.

Chapter 2 Literature Survey

Nanoindentation, also known as depth-sensing indentation (DSI), is a very popular technique for determination of mechanical properties such as elastic modulus and hardness. It has been extensively used to study the behavior of metallic or ceramic materials in the past couple of decades. Local mechanical properties at the micro- and nanoscale can be effectively characterized by nanoindentation, which is the major advantage of using this technique [22, 23]. This also makes it ideal to study materials that are otherwise not characterizable by conventional testing methods e.g. thin films, coatings, and localized surface modification of materials [24, 25, 26]. Nanoindentation has also attracted interest for biological material characterization, since it may be used to assess mechanical properties on the cellular scale [27].

Two different approaches have been primarily used for mechanical characterization of materials by nanoindentation [28]. The first approach is based on analytical or semi-analytical solutions arising from mathematical contact theories. The second approach, which is popularly known as ‘inverse analysis’, utilizes a combination of finite element methods and numerical optimization algorithms. In inverse analysis the difference between experimental and numerical nanoindentation data, called the objective or error function, is minimized with respect to the material model parameters using numerical optimization. Subsequently, the parameters of the constitutive models are identified as the optimized material properties. Inverse analysis has been found to be applicable in tackling a wide range of problems by the research community [29]. In the next few subsections, a brief review of the nanoindentation based studies is presented for both methods.

2.1 Nanoindentation Analysis: Oliver–Pharr Method

Theoretical studies to characterize the material properties by indentation were first conducted by Hertz. He developed a relationship between the load and indentation depth for spherical elastic bodies. Later Sneddon extended Hertz’s work to derive expressions for load, displacement, and contact depth for elastic contacts between a rigid, axisymmetric punch with an arbitrary smooth profile and an elastic half-space [30]. The first study to use Sneddon’s analytical solution and measure the mechanical properties from nanoindentation experiment was conducted by Doerner and Nix [31]. Their study demonstrated that hardness and Young’s modulus could be calculated based on the information provided by nanoindentation load–displacement plot. They also pointed out that with the help of suitable analytical procedure plastic properties of a material can also be obtained from nanoindentation.

In subsequent years, Oliver and Pharr modified the method proposed by Doerner and Nix to find elastic properties of materials [2]. This method has since been cited for more than 13000 times and became more of an unofficial standard for nanoindentation testing. The underlying assumption of this method is that unloading curve of a nanoindentation plot is purely dominated by the elastic properties of the material. Using this method for time-dependent materials would provide inaccurate results since the original assumption does not remain valid. To provide a better understanding of nanoindentation technique a brief overview of this method is followed in next paragraphs.

Figure 2.1: a) Typical nanoindentation load-displacement plot, b) schematic of the material surface before and after loading [2]

Figure 2.1 shows a typical nanoindentation load–displacement plot. In order to extract mechanical properties, such as Young’s modulus and hardness, values of contact stiffness, contact depth, and area of contact are required from the nanoindentation plot. The contact stiffness S is the slope of the unloading curve, while the contact depth (the depth of actual contact between the indenter and the material) is calculated by Eq. 2.1.


Here, is the maximum depth of penetration including elastic deformation of the surface under load, is the maximum force, and is a geometrical constant associated with the geometry of the indenter [2]. This method of determining the contact depth is commonly referred to as the Oliver–Pharr method; a schematic in Fig. 2.1 shows and .

Once is determined, the projected area A of actual contact can be calculated using the cross-sectional shape of the indenter along its length. Determining accurate contact area is found to be crucial for elastic analysis of nanoindentation data [32]. This area function could be determined by direct measurement of the imprint geometry under a scanning microscope [33], but in practice is normally determined by indenting a reference sample and iteratively fitting the results. The relationship between contact area, A and the contact depth, for a Berkovich tip is generally expressed by the following equation–


here, the coefficients can be determined by iterative fitting to indentation measurements conducted on reference material such as fused silica. The first term in Eq. 2.1 represents the area–depth relationship for a perfectly sharp Berkovich indenter, while the other terms account for tip imperfections e.g. tip roundness. Once the area of contact A is determined, hardness is found using a simple equation-


It is important to note that this hardness is defined using the projected area of contact under load, while macroscopic definition of hardness is force divided by the area of the residual imprint left by the indenter. For most materials the two definitions yields very similar values. However, in case of a material showing little to no plastic flow, the hardness calculated by Eq. 2.1 tends to be lower than the macroscopic definition.

Once contact area, A and contact stiffness, S is known, Sneddon’s solution can be adapted independent of the geometry of the punch, and Young’s modulus can be calculated using the following equation for reduced modulus:


here, is a small correction for the non-axisymmetric indenter shape (e.g. = 1.034 for a Berkovich tip). For a perfectly elastic–plastic material with no other form of deformation present, the unloading curve is purely dominated by the elastic recovery of the material. As a result, Young’s modulus determination from an unloading curve of a nanoindentation experiment becomes possible. Later Field and Swain developed means of extracting both Young’s modulus and yield strength from load–displacement curves of a spherical indentation [34].

Approximation of Sneddon’s solution is that the indenter is rigid, and therefore, deformation of the indenter is small and insignificant compared to the material being tested. As long as this approximation is valid Sneddon’s solution can yield good results for the reduced modulus. However, in case of testing very hard materials, such as diamond-like carbon the deformation experienced by the indenter is substantial, thus violating the approximation Sneddon’s solution is based upon.

Oliver–Pharr method and subsequent developments provided a means for extracting few key material parameters from a nanoindentation plot, namely Young’s modulus, hardness, and yield strength. These were groundbreaking developments in terms of characterizing elastic–plastic material behavior. However, these methods are unusable for viscoelastic materials as the underlying assumption of unloading curve purely dominated by elastic recovery no longer holds for viscoelastic materials. In addition to that, hardness, modulus, and yield strength properties are inadequate to represent the full spectrum of behavior for these materials.

2.2 Adaptation of Oliver–Pharr Method for Time–Dependent Behavior

The most widely used indenter load or displacement profiles are the triangular, where the load or displacement is ramped at a certain rate to the maximum value and then unloaded back to zero, as shown in Fig. 2.2. For elastoplastic material systems (e.g. most metals and ceramics) exhibiting little to no time–dependent behavior, the load–displacement nanoindentation curve is insensitive to loading or unloading rates; thus, triangular profiles can be effectively used to characterize these materials.

However, this is not the case for testing of viscoelastic materials such as polymers and biomaterials due to the fact that viscous behavior of these materials dramatically affect the load–displacement curve. The inherent time–dependency in mechanical response of these materials make the unloading curve of the nanoindentation experiment noticeably different by producing a“nose” [33, 35].

Figure 2.3 shows a typical nanoindentation load–displacement plot for a viscoelastic material. The nose results from excessive creep of a material under the indenter, which dominates over the elastic recovery of the material as the tip retracts from the surface. Applying Oliver–Pharr method on a nanoindentation plot exhibiting nose often provides a negative value for contact stiffness, S, and prevents extracting elastic modulus altogether. Even without the appearance of the nose, the presence of viscoelasticity often leads to overestimation of Young’s modulus.

Figure 2.2: Load-displacement profiles a) triangular, b) trapezoidal [3]
Figure 2.3: Typical nanoindentation plot for a viscoelastic material [4]

A trapezoidal load or displacement profile that implements a long enough hold before the unloading has been found to suppress the creep behavior near the initial unloading part [35, 36, 37, 38, 39]. The holding period ensures complete relaxation of the material, and minimizes the viscoelastic recovery during the unloading.

An useful modification to Oliver–Pharr method was proposed by Ngan et al. so that true value of Young’s modulus for viscoelastic materials can be extracted from nanoindentation experiment [35]. According to their study, for a load-controlled indentation test with hold period prior to unloading, contact stiffness can be corrected using the following equation—


where dh/dt is the indenter displacement rate at the end of the load hold just prior to unloading, S is the contact stiffness found via original Oliver–Pharr method, and = is the initial unloading rate.

Although, Ngan et al.’s method provided an useful way to use Oliver–Pharr method for characterizing viscoelastic solids, it simply cannot address the various other important properties of a viscoelastic material [40]. As a result developing dedicated analysis techniques for viscoelastic material characterization via nanoindentation has been one of the most popular research area of the past decade.

2.3 Analytical Approaches for Viscoelastic Materials

Analytical solutions capable of characterizing viscoelastic material behavior from nanoindentation load–displacement plot originated from the early works of Radok [41]. He was the first to tackle the linear viscoelastic contact problem using the method of functional equations or hereditary integrals, which was later completed by Lee and Radok [42]. This method of functional equations solved the viscoelastic problem by replacing the elastic constant with their corresponding viscoelastic operators. This is why this method is also known as ‘Correspondence Principal’. Radok extended the ‘Laplace transform method’ formulated by Lee to eliminate the explicit time dependence of the viscoelastic contact problem and solved it in the Laplace domain [43]. Before Lee and Radok, Laplace method was only applicable to problems where displacement and stress boundary conditions are unchanged e.g. flat punch indentation problem.

The method of functional equations proved to be very successful in formulating analytical solutions for viscoelastic bodies; however, the solutions were only valid as long the penetration depth in a viscoelastic indentation monotonically increased [42]. Hence, this method is only valid for the loading portion of the nanoindentation plot. Many researchers attempted to remove this restriction. Hunter was able to remove it for spherical indentation, while Ting’s implicit equations were able to remove it altogether for any linear viscoelastic material tested under any axisymmetric indenter shape [44, 45]. However, except for few specific cases applying Ting’s formulation is a challenge. As a result, closed form solutions for linear viscoelastic problems are still being formulated using Radok’s method of functional equations.

In 1985, Johnson summarized the correspondence analysis of spherical indentation replacing the elastic constants by the Boltzmann viscoelastic hereditary integral operators. Based on these approaches, in several contributions [46, 47, 48, 49, 50, 51, 52], the viscoelastic analytical solutions of nanoindentation with different indenter tips were presented. However, the ‘Correspondence Principal’ method is restricted by yielding accurate identification only for specific linear viscoelastic models under fixed experimental processes.

2.3.1 Spring–Dashpot Based Viscoelastic Models

A key aspect of viscoelasticity is that mechanical behavior of a material can be successfully modeled using a combination of springs and dashpots. Viscoelastic materials demonstrate both elastic and viscous behavior in the same material. If spring represents the elastic behavior and dashpot represents the viscous behavior, a combination of two could be able to model the behavior of a viscoelastic material. The biggest advantage of using spring–dashpot based model is that the viscoelastic models can be tailored to suit specific observations [53].

By putting this idea to use several spring–dashpot model has been proposed in the literature. For a viscoelastic material, stress level is related to both strain level and strain rate in the following general form


where, and are the strain and stress levels, respectively, and t is the time. and are the coefficients that determine the linear or even non-linear stress–strain behavior.

In the most simplest of forms, where one spring element and one dashpot element is used to create a model, this technique leads to two well known models, namely Kelvin–Voigt and Maxwell model. These models assume linear stress–strain relationship. Figures 2.5 and 2.4 show Kevin–Voigt and Maxwell model, respectively.

Figure 2.4: Kelvin-Voigt model
Figure 2.5: Maxwell model

Using the corresponding constitutive equations for spring and dashpot, the following equation can be developed for Kelvin–Voigt model.


here, =E is the creep part and is the recovery part of the model; = is the strain rate; E is the rigidity modulus, and is the coefficient of viscosity.

Under constant stress conditions, the strain response of the material can then be captured as an exponential decay function


However, under constant strain rate conditions (stress relaxation part), the Kelvin–Voigt provides unrealistic linear elastic behavior for the viscoelastic material. The Maxwell model, however, povides better approximation for constant stress relaxation. For Maxwell model, the constitutive equation comes in the following form-


In case of stress relaxation (=0), an exponential decay of stress is found,


while in a recovery experiment (=0), the model predicts the basic equation of pure Newtonian flow.

Figure 2.6: Schematic of Generalized Maxwell model

Real world materials are much more complex in their behavior, and these simplistic models are not sufficient in representing that. Previous studies show that modeling creep and relaxation behavior of complex viscoelastic behavior requires an assembly of multiple spring and dashpot in the model [54]. One such model is generalized Maxwell model. Various studies have reported that generalized Maxwell model worked well in terms of modeling the viscoelastic behavior [55, 51, 56, 57, 58]. For this model relaxation can be written in the general form


where, Y(t) is the relaxation function.

The relaxation function can be represented using Prony series having the following expression:


where, is the Prony constant, is the Prony retardation constant, is the instantaneous modulus. Prony coefficients are usually found by nonlinear regression, which allow adjusting the model with respect to the observed behavior.

2.4 Inverse Approach Based Material Behavior Modeling

Inverse problems are defined as the problems where the output is known and the input or source of output remains to be determined. They are contrary to the direct problems, in which output or response are determined using information from input [59]. In order to analyze inverse problems, experimental data obtained under known boundary conditions are compared with the calculated ones. The combination of nanoindentation and FEA has proved to be a powerful analysis tool for soft polymers such as gels, and coatings, and for soft tissues [60, 61, 62, 63].

Inverse analysis requires an optimization algorithm to extract the set of parameters for which the objective function (difference between simulation and experimental load-displacement data) attains the minimum value. The choice of the optimization algorithm for minimizing an objective function is a topic of interest. Whenever possible it is better to employ global optimization techniques. There are many variants of Simulated Annealing or Genetic Algorithm based global optimization scheme, such as evolutionary algorithms, or deterministic algorithms like the Simplex method. These algorithms have proven to be very useful in case of optimization problems where user has no prior information about the location of the solution in the parameter space, thus incapable of making a priori choices about the initial estimates.

However, in case of finite element analysis, where time required to run one single analysis can range from few minutes to even days, the success of global optimization methods come at a price of astronomical computational cost. In these cases, local optimization algorithms could prove to be useful given that the quality of initial estimates are good. However, these algorithms are gradient-based and involves computationally costly calculations of second order partial derivatives of the objective function. For an objective function f, Hessian matrix, H is defined by the Eq. 2.4.


Determining H often requires high computational cost. The cost involving calculation of Hessian matrix can be reduced if it is approximated using the Jacobian matrix, J (Eq. 2.4).


The biggest disadvantage of gradient-based scheme is that algorithms can sometime get trapped inside a local minima. In addition to that, due to the ill-conditioned nature of inverse analysis, identifying the correct minima from a set of local minima is troublesome. Singularity in approximated Hessian matrix and non-covergence are few other problems that often trouble the local optimization techniques.

Figure 2.7: An illustration of trade-off between fidelity and computational cost [5]

2.5 Review of Existing Literature

2.5.1 Analytical Approach in Nanoindentation

A large number of studies have been conducted in an effort to determine viscoelastic behavior using nanoindentation. Cheng et al. derived the analytical solutions for linearly viscoelastic deformation and provided a method to measure viscoelastic properties described by a three-element standard linear solid (SLS) model using a flat-punch indenter [46]. Hernandez-Jimenez et al. studied the viscoelastic behavior of PMMA and PTFE using Maxwell model [64]. Lu et al. developed methods to measure the creep compliance of PMMA and PC polymers using spherical or Berkovich indenter by deducing closed form analytical solutions using Generalized Maxwell model [51]. Prony series parameters for the stress relaxation or creep compliance was found by curve fitting only the loading portion of the nanoindentation plot. Measurement of creep compliance from conventional tension and shear tests were compared with the nanoindentation technique, where reasonable agreement between the values from different techniques was observed. Fisher-Cripps developed creep compliance analytical closed form solutions for three-element Voigt model, four-element Maxwell–Voigt model [23]. Gonda et al. spherical and conical indentations on a thin polymer film on a substrate, where viscoelastic properties found by analytical equations through correspondence principal and the results were verified using finite element modeling [65]. Vanlandingham et al. investigated linear viscoelastic material analytical solutions for epoxy, PMMA and PDMS polymers, and compared the values obtained from nanoindentation with values from rheometry measurements [52]. Cheng and Cheng derived closed form solutions relating the initial unloading slope, instantaneous relaxation modulus, and contact depth for linearly viscoelastic material under a rigid and arbitrary axisymmetric shape [66]. In another study, they also derived the relationship between unloading slope, contact depth and instantaneous modulus for conical indentation [47]. Three parameter Standard Linear Solid (SLS) model has also been used to derive equations for spherical nanoindentation of viscoelastic materials [67]. Cheng et al. addressed the linear viscoelastic material indentation using three-parameter Maxwell solid [68]. Finite element simulations were conducted using these relationships to verify the solutions. Zhou and Lu developed methods to measure creep compliance and relaxation modulus in both time and frequency domains under constant and ramped loading conditions using a spherical indenter [69]. Vandamme et al. modeled using 3-parameter Maxwell model, the 4-parameter Kelvin-Voigt model and the 5-parameter combined Kelvin-Voigt-Maxwell model and derived the closed form analytical solution for conical indenter [70]. Wei et al. studied the viscous behavior of PMMA and PU materials using a combination of Kelvin–Voigt model and a dashpot [71]. Their model accounted for the irreversible delayed plastic (viscoplastic) deformation, irreversible viscous deformation, and reversible delayed elastic (viscoelastic) deformations. Oyen and Cook examined the creep displacements as a function of time for PMMA and a few other polymers using constant loading and unloading rates [72]. They also examined the effect of triangular and trapezoidal loading profiles. For pyramidal indentation tests, a new method for estimation of time-constant was proposed. Liu et al. developed a model based on Burgers model and applied to understand the viscoelastic behavior of soft polymers like PMMA [54]. According to them, Burgers model provided the best agreement with the experimental data in comparison to simple Maxwell or Kelvin-Voigt model. In addition to that, they also indicated that the nose formation at the beginning of the unloading stemmed from the decrease of the viscosity parameter. Jager et al. characterized viscoelastic properties of bitumen using different spring-dashpot models for real tip geometry of the indenter [50]. Linear viscoelastic analysis based on spherical indentation experiment has also been carried out on human tympanic membrane [58]. Lin et al. studied viscoelastic behavior of PDMS micro pillars using uniaxial, DMA, nanoindentation tests, where generalized Maxwell model was used to describe the viscoelastic behavior of the material [57].

Mencik et al. analyzed the viscoelastic–viscoplastic behavior of material under indentation for different indenter profiles [73]. They found that materials under sharp indenter undergoes high stresses and exhibits viscoplastic effects. Chen et al. used dimensional analysis and finite element analysis to understand the effects of residual stress, substrate, and creep behavior on the load-displacement data [74]. Through development of some analytical solutions they showed that it is possible to obtain not only Young’s modulus and hardness, but also viscoelastic properties and residual stress.

2.5.2 Inverse Approach in Nanoindentation

The very first instance of applying inverse method for an indentation-based study was probably by Knapp et al., where they studied the elastic–plastic behavior of Al under nanoindentation [75]. Their study showed that it was possible to extract modulus, yield strength, and hardening coefficient from the nanoindentation data of thin films using FEA based inverse analysis independent of the effect of substrate. Later, Huber et al. employed Artificial Neural Network (ANN)-based inverse analysis to extract material parameters from an indentation experiment of metals [76]. From that point onwards, inverse FEA-based analysis has been used to extract material properties for different classes of materials, such as, isotropic and anisotropic elastic–plastic materials [77, 78, 79, 80], linear viscoelastic materials [81, 82, 83, 84], hyper-elastic materials [85, 86, 87, 88], nonlinear viscoelastic materials [89, 90], etc.

After being introduced by Knapp et al., inverse FEA technique has been reported in numerous publications dealing with material property extraction for elasto-plastic materials. On the contrary, the number of studies that tackled viscoelastic nanoindentation using inverse analysis is found to be very low. The probable reasons could include the lack of understanding about the viscoelastic constitutive relationships, the high number of model parameters needing to be optimized, etc. While elasto-plastic behavior in materials has been studied for a long time, viscoelasticity is being studied only recently fueled by the recent interests in understanding polymers and biomaterials. As we are interested in the viscoelastic materials, this part of the literature will only review the related studies in the area of viscoelasticity.

Ovaert et al. studied on viscoelastic properties of thin polymer coatings using three- and four-parameter viscoelastic models using indentation and inverse analysis [61]. Kim and Srinivasan et al. extracted Fung’s QLV model parameters for soft tissues using two step parameter optimization process [91]. Hartmann et al. used uniaxial test data for viscoplastic parameter identification and validated those using indentation test data [92]. Samur et al. studied the viscoelastic behavior of pig liver tissues using inverse analysis [93]. Resapu et al. extracted Prony series parameters for the relaxation behavior of PVC and PE in indentation tests [94]. Guessasma et al. determined viscoelastic properties of biopolymer composite materials [81]. Liu et al. characterized viscoelastic behavior of soft gels using Kelvin–Voigt model [82]. Rauchs identified viscous hyper-elastic and elasto-viscoplastic material parameters from indentation tests [95, 88]. Abyaneh et al. characterized porcine cornia using linear viscoelastic model [84]. Viscoelastic Arruda–Boyce constitutive model has also been studied with AFM indentation and inverse FE analysis for porcine zone pellucida [96]. Rayleigh dissipative function has been used by Abetkovskaia et al. to develop AFM based viscoelastic characterization of soft materials [97]. Valdez-Jasso et al. used inverse analysis to characterize viscoelastic behaviors of ovine aorta, where the viscoelastic behavior was modeled using arctangent and sigmoid viscoelastic models [83]. Recently, Kucuk et al. used nonlinear Burgers model to characterize the nonlinear viscoelastic behavior of PMMA and PVAc [89, 90].

Inverse analysis of nanoindentation data is challenging due to various reasons. One of this big challenge is to find out unique solution. In case of non-unique solutions, two approaches were found to be effective. In the first approach, additional information from the nanoindentation experiment is gathered, and used in the objective function. These information can include imprint geometry [77, 98, 99] or pile-up/sink-in [80] information. The other approach is to use multiple indenters with different geometry. This method has also been found to beneficial in providing unique solutions from the inverse analysis [100, 78, 28].

2.6 Motivation and Objective of the Study

Nanoindentation has the potential to become a very effective material characterization tool given that appropriate analysis technique with proper constitutive relationship is used. In the last two or three decades this technique has come a long way in terms of applicability for metallic or ceramic material characterization.

However, suitable analysis technique for materials such as polymer or soft tissues is still lacking. Table 2.1 summarizes the results found from various studies that used nanoindentation technique for characterizing soft tissues. It can be seen that, for almost all the studied tissues the value of Young’s modulus varied by few orders of magnitude. Part of the variability comes from the difference in experimental design and sample preparation, while most of it stemmed from the fact that these materials exhibited time–dependent deformation behavior [1]. If consistent strain rate were to be used in the experiment a more consistent Young’s modulus could probably be found.

Even if we consider that the Young’s modulus could be extracted reliably independent of viscoelastic influence, it would only serve as a partial knowledge about the material system. Young’s modulus only quantifies the intrinsic elastic behavior of a material, which limits its usefulness only to metals or crystalline solids.

Tissue Range
(kPa) Average
modulus (kPa) Reference
Liver and Kidney 0.6-760   190  [101, 102, 103]
Artery and Vein 6.5-560   125  [104, 105, 106]
Skin 6-222   85  [101, 107, 108]
Cornea anterior base 7.5-50  29  [109]
Breast tissue 0.167-29   8  [6, 7, 110]
Muscle 2-12   7  [106, 111]
Spinal cord and gray matter 0.2-7   3  [112, 113]
Table 2.1: Young’s modulus of soft tissues, measured by indentation [1]

In the attempts to understand or characterize the time–dependent properties inherent to soft polymers and biomaterials, most researchers simplified the behavior of these materials as linear viscoelastic. In fact most of the studies that used nanoindentation to characterize viscoelastic materials used linear viscoelastic theory developed through ‘Correspondence Principal’. In addition to the fact that material behavior is simplified as linear viscoelastic, correspondence principal based analytical solutions has further limitations i.e. useful only till contact area increases monotonically (loading portion of the nanoindentation plot). As a result, this method fails to address how viscous behavior affects the unloading curve of the nanoindentation experiment, although substantial amount of information about the material behavior is present in the unloading portion of the curve.

To best of our knowledge, no analytical or closed-form solutions (in either differential or integral form) exist for indentation of quasi-linear or nonlinear viscoelastic material. However, soft tissues and polymers are generally nonlinearly viscoelastic [114, 56], where the creep compliance or relaxation modulus are a nonlinear function of both time and applied stress or strain. In these cases, an appropriate constitutive law should be used to describe the distinct behaviors of these materials [115]. Due to the fact that, no closed form solution can be obtained for nonlinear viscoelastic behavior, many researchers tried modeling the behavior of the material using Fung’s Quasilinear Viscoelastic (QLV) model [114]. Fung’s model however considers the material to be nonlinear only with respect to strain [116], and fails to represent the full spectrum of nonlinearity of the material.

The closest work that tackled nonlinear viscoelastic behavior of the material was by Kucuk et al. [89, 90]. In these studies, a nonlinear viscoelastic model based on modified Burgers model was used. The unknown model parameters were then obtained using inverse analysis. However, the authors did not provided any information about the inverse analysis procedure that was followed. Without such key information obtaining the values of the model parameter for other material system is difficult. In addition to that, their study utilized quite a high number of parameters in the nonlinear model without providing any information about whether all the parameters were required to capture the behavior or not.

To understand the full spectrum of mechanical behavior in soft biomaterials and polymers, an study is thus required which would improve on the limitations of previous studies. Because without understanding the mechanical behavior, it would be impossible to predict the behavior of these materials under complex loading scenarios.

Figure 2.8: Viscoelasticity a) linear and b) nonlinear

In order to develop a nonlinear viscoelastic model for soft materials, we propose to implement finite element analysis with inverse analysis. The process is called the inverse analysis because it is the opposite of an ordinary simulation (i.e. solving for forces or displacements given material parameters and boundary conditions). The inverse method permit us to treat any material models with nonlinear properties and to include further affecting factors in the numerical model. The rate-dependent properties of materials can be more accurately identified using the inverse method.

The whole process of developing nonlinear viscoelastic model can be subdivided in few steps. In the first step, a finite element model of the nanoindentation experiment is required which can effectively simulate the experiment. For this work, we have chosen commercially available finite element analysis software–ABAQUS. Confidence was established on the ABAQUS representation of the nanoindentation experiment by comparing the simulation results with the well established analytical solution from contact theory.

In the next step, an appropriate spring–dashpot system for describing these kind of materials needs to be developed. The associated mathematical model for the spring–dashpot system has to be incorporated in the ABAQUS simulation of the nanoindentation experiment via user-defined subroutine called UMAT.

In the final step, an optimization based algorithm needs to be established, which will be able to minimize the difference between the simulated and experimentally found load–displacement data. This study will use both the loading and unloading portion of the nanoindentation experiment in the model development process; because unloading curve would provide additional constraints which a successful model must satisfy. In addition, this will provide additional validation of the viscoelastic model parameters extracted from nanoindentation data.

One of the main issue in inverse analysis based model development is the high computation cost associated with the optimization process. In case of a nonlinear model, the number of model parameters that needs to be optimized is usually high. In addition to that, nonlinear FEA analysis requires considerably higher computational cost due to the continuous updating procedure of global stiffness matrix. This updating process is a serious drawback for FEA-based realtime optimization applications [12]. That is why, improving the computational efficiency in the inverse analysis process is another important objective of this study.

Chapter 3 Development of the Technique

3.1 Nonlinear Viscoelastic Mathematical Model

There are a few nonlinear viscoelastic models in the literature but to date it appears that none of these models can describe the nonlinearly viscoelastic behavior of a polymer under all loading and environmental conditions [117]. Under a given set of loading conditions, however, an appropriate nonlinearly constitutive model can be used to model the viscoelastic response.

In this study, a spring–dashpot model suggested by Marin and Pao [118] was used. In linear case this model is generally called four-parameter Burgers model [119] and it is formed by a serial connection of a Maxwell element to a Voigt element as seen in Fig. 3.1. For an increased relaxation spectrum, the viscoelastic response can be modeled by increasing the number of Voigt elements.

The nonlinear characteristic is introduced when the dashpot constants ( and ) take values other than unity. In the three-dimensional model, the total strains are calculated as the summation of the elastic (e), transient creep (t), and steady creep strains (s[120]. In this study, the nonlinear creep deformation is assumed to be incompressible. Under these assumptions, the three-dimensional nonlinearly viscoelastic law can be expressed as:

Figure 3.1: Schematic of nonlinear Burger’s model

where , are the Young’s modulus and Poisson ratio, respectively; is the second invariant of the deviator stress tensor s; , , , , are the nonlinear material parameters. is the Cauchy stress tensor; i, j are the indices ranging among 1, 2 and 3. is the Kronecker delta which used in the context of summation convention with the well-known property when i = j and . Small deformations are assumed in the formulation. When more than one Voigt element is included in the model, the total strain components can be given as the sum of elastic, steady creep, and transient creep components for all Voigt elements,


where n is the number of Voigt elements shown in Fig. Equations 3.1 and 3.1 can also be written in integral form:


ABAQUS/Standard finite element code is used as the implementation platform. Although ABAQUS has a rich material library for various applications, a nonlinearly viscoelastic model suitable for this work was not available. In this study, a UMAT was developed in order to implement the nonlinear Burgers model. UMAT requires the tangent stiffness matrix of the material model for finite element calculations. For implementation of the nonlinear Burgers viscoelastic model, the UMAT involves mainly temporal discretization. This was done following the procedure implemented by Kucuk et al. [89, 90].

A simple, stable integration operator for these equations is the central difference operator:


where is a function, is its value at the beginning of the increment, is the change in the function over the increment, and is the time increment.

Tangent stiffness matrix of the constitutive model, with being the stress increments and the strain increments, can be derived by applying central difference operator to the rate-dependent constitutive equations (Eq. 3.13.1).

Applying the central difference method to the elastic strain component as depicted in Eq. 3.1, yields


If the elastic Hooke’s law is defined by Eq. 3.1, the elastic compliance matrix, C is defined by Eq. 3.1.


Similar procedure as applied to Eq. 3.1 for steady creep component gives


Assuming , we have


Since , we have


The compliance matrix of steady creep then can be written as


Finally for the transient creep component as defined in Eq. 3.1, we have


The compliance matrix of transient creep can then be written as


From Eq. 3.1, the total compliance is now


By investigating the total compliance matrix, system tangent stiffness matrix (Jacobian matrix) can be obtained from Eq. 3.1. It should be noted that the Jacobian matrix in Eq. 3.1 accounts only for the elastic deformation and creep deformation caused by load or stress increment. It is seen from Eq. 3.1 and 3.1 that the aforementioned creep strain is just a small part of the total steady and transient creep strain. The rest of the creep strain is developed over the time period during the time increment and controlled by the applied stress. An artificial stress increment is introduced to include this creep strain in the system equation. This part of creep strain can be extracted from Eq. 3.1 and 3.1 as


A stress increment is then added into the system equation to account for the creep strain in Eq. 3.1, with C being the Jacobian stiffness matrix calculated from Eq. 3.1.

3.2 Finite Element Modeling

The 3D finite element model of nanoindentation experiment was constructed using commercial finite element package ABAQUS (Dassault Systémes, Providence, RI). The indenter in a nanoindentation experiment is made with diamond and possess very high Young’s modulus. So, it is possible to model the indenter as analytical rigid body. Finite element solver does not require calculating stress and strains in an analytical rigid body, hence reduces the computational time.

Berkovich indenter can also be modeled as a 2D axisymmetric conical indenter with an effective cone angle [23]. The effective cone angle is calculated in a manner so that it provides the same area to depth relationship as the actual Berkovich indenter. The benefit of using a 2D model is that it requires less computation time compared to a 3D model. Nonetheless, a 3D model was implemented in this study to achieve higher accuracy in simulating the nanoindentation experiment.

Even after adopting few simplifications, modeling nanoindentation experiment is still very challenging due to the several nonlinearities associated with the experiment (boundary, geometrical and material nonlinearity). Studies showed that, in case of modeling complex geometries, it is beneficial to model rigid elements as discrete rigid body rather than analytical rigid body. So, the Berkovich indenter was modeled as discrete rigid body, while the sample was modeled as deformable body.

To ensure accuracy of the simulation results, the sample was modeled with finer mesh near the contact area where the stress and strain generated was much higher due to the singularity dominated zone. The contact between the indenter and the sample was defined as surface-to-surface contact, where the indenter was designated as master surface and the sample was as designated as slave surface. The element types for the sample was chosen from the eight-node brick element family (C3D8). Material behavior of the sample was defined in the model using a subroutine called (UMAT). The mathematical development of a nonlinear viscoelastic constitutive relationship was required to code the UMAT, which development was covered in the previous section.

3.3 Inverse Analysis

In order to facilitate the identification of global solution in the parameter space, our study implemented surrogate modeling approach. Surrogate models, also known as metamodels are particularly useful in case of finite element based inverse analysis.

Figure 3.2 shows the typical workflow of an inverse analysis for nanoindentation based model parameter extraction. Due to the fact that in every iteration of the inverse analysis one finite element analysis is required, the high computational cost involving the inverse analysis becomes the limiting factor in determining the correct solution. If finite element analysis can be replaced with a surrogate model, which is a numerical approximation of the input–output relationship, the total computational cost can be dramatically reduced. In a nutshell, use of surrogate model can effectively reduce the computational cost while still keeping the fidelity of the solution adequately high.

Figure 3.2: Typical inverse method flowchart

In this study, surrogate model is built by utilizing two numerical techniques named as Proper Orthogonal Decomposition (POD) and Radial Basis Function (RBF). Proper Orthogonal Decomposition (POD), also known as Principal Component Analysis (PCA) technique can be used with either experimental or simulated field data to derive a reduced-order set of basis functions capable of being used in a numerical representation of the system [121]. POD reduced-order approximation has been shown to provide accurate numerical representations for complex systems with minimal computational cost [122, 123, 124]. In addition, POD has been applied to several inverse problem methodologies, such as optimal control  [125, 126, 127, 128], and nondestructive testing and system identification  [129, 130, 131, 132, 133]. However, work has yet to be shown (to our knowledge) for using POD reduced-order modeling for inverse viscoelastic material characterization from quasi-static indentation testing.

As a means of correlating unclear data using only spatial lines and planes, the concept of proper orthogonal decomposition (POD) was first developed over a century ago as a statistical tool by Pearson [134]. Since then POD has been redeveloped under various names and has been used in numerous different applications from signal processing and control theory, human face recognition, data compression, parameter estimation and many others [135]. POD is also known as Karhunen–Loeve Decomposition (KLD), Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) [136, 135, 137]. In the recent past, POD has been increasingly used in many engineering applications ranging from computational fluid dynamics (CFD) to modeling of heat transfer problems due to its ability to reduce computational burden while maintaining adequately high accuracy.

For simple understanding of the POD technique, one should imagine a collection of vectors inside a Cartesian coordinate system. If these vectors are parallel to one another it could be assumed that these are correlated. On the other hand, uncorrelated would mean that these vectors are orthogonal (or perpendicular) to one another.

POD’s major objective is to rotate the coordinate system in such a manner so that the least amount of coordinates can be used to define the system. As an example, we know that a vector in cartesian coordinate system requires two projections (x- and y-axis projection) to be effectively defined. However, if the coordinate system is rotated only one projection can define the same vector. In case of complicated data sets, the number of rotated coordinates would be higher for effective representation of the data. In such cases, POD captures the maximum projection of the vectors in the first rotated coordinate, which is commonly referred to as the first principal component. The second axis in the POD frame, called the second principal component, captures the next orthogonal direction with the largest projection and so on.

POD is completely data dependent and does not assume any prior knowledge of the process that generates the data. This property is advantageous in situations where a priori knowledge of the underlying process is insufficient. POD does not neglect the nonlinearities of the original vector field. If the original system is nonlinear, then the resulting POD reduced order model will also be nonlinear.

3.3.1 Proper Orthogonal Decomposition (POD) Theory

If a function is needed to be approximated over some domain of interest, it can be written as the following equation through a linear combination of few basis functions .


Here represents the unknown coefficients. Once basis functions are known, the coefficient values are obtained in a least square means.


For any function , number of choices can be made regarding the selection of basis function. Based on one’s expertise and prior knowledge about the system being represented, one can often opt for a basis constructed from polynomial, trigonometric, or exponential functions. Proper Orthogonal Decomposition (POD) is one such technique that can be used to construct the optimal basis for a function under investigation in a least square sense.

The derivation of POD presented in subsequent paragraphs refers to arbitrary case of vectors with dimensionality N>2. The notations presented in this section is congruent with Buljak [138].

POD starts with the idea of snapshots. Snapshots can be defined as an one-to-one relationship between the input and output of a system. In a typical scenario of an inverse finite element nanoindentation simulation, snapshots are the relation of material model parameters and output load–displacement data. In more concrete definition, a snapshot will be a collection of N discrete values of a certain state variable resulting from a simulation (which represents a system) collected in vector , corresponding to some input parameters (collected in vector ) on which the solution depends. A system can also be represented by an experiment, where snapshot will store the measurements taken from an experiment.

Figure 3.3: Input output relationship in a typical system

Further, a set of M different snapshots, corresponding to different input parameters, can be collected in a rectangular NM matrix U, called the snapshot matrix.

Therefore, a snapshot matrix represents a collection of responses of one system, under given conditions, corresponding to different values of parameters on which the solution depends. This snapshot matrix can be interpreted as a set of M, N-dimensional vectors. Each vector corresponds to one parameter combination. In the context of inverse finite element analysis for model parameter extraction, it can be said that the inputs to the system that are changing from one snapshot to another are some parameters entering into the constitutive model of material, while the boundary conditions and initial conditions are the same for all of the snapshots. So, u contains N number of individual displacement data for N number of corresponding load increments, while the whole snapshot matrix, U contains M number of individual finite element simulations.

It is reasonable to expect that there will be a strong correlation between these snapshot vectors since they represent the outputs of the same system where just some material model parameters are changed. The POD theory can be effectively applied on the snapshot matrix, allowing to construct a new basis in which the dimensionality can be drastically cut-down to K N. POD finds the most accurate representation in some subspace W with the dimension of KN. If we denote , ,…, as the orthonormal basis of the subspace W, then each vector from the original set can be written as


where are amplitudes corresponding to vector in new subspace W, and is matrix that collects all the orthonormal basis of the subspace . In a least square sense the error of approximation then becomes


Eq. 3.0 provides the error for only the vector. For all the vectors in the snapshot, total error is expressed by the Eq. 3.0:


The orthonormal basis has to be chosen in such a manner so that total error is minimized. To do that, first order derivate of total error with respect to all the unknowns (namely ) are needed. Taking partial derivative of total error:


By substitution of in Eq. 3.0:


Few more mathematical manipulation provides:


where C is called the covariance matrix defined as . The first part of Eq. 3.0 is a scalar constant which depends on the original set of snapshots. So, in order to reduce the error of approximation one has to maximize under the constraint of orthonormality of the new basis . By using Lagrange multipliers method the constrained problem can be converted into


In order to maximize Eq. 3.0 first order derivatives with respect to is required. By doing that we find


Eq. 3.0 is only satisfied if is eigenvector and is the corresponding eigenvalue of matrix C. Now from taking Eq. 3.0 and Eq. 3.0 into consideration the total error equation can be changed to


Recalling that the first term is a constant it results that the error of approximation is minimized if the new basis is constructed of K eigenvectors that are corresponding to the first K largest eigenvalues of covariance matrix C


If the subspace W is constructed with all the eigenvectors of matrix C, there is no error of approximation because in that case all the vectors are just expressed in a different coordinate basis. Approximation in any other subspace that uses smaller number of eigenvectors the error of approximation is found using the following equation


which is the ratio between the summation of kept eigenvalues and summation of all the eigenvalues.

In this study, POD is used to determine the displacement of the indenter tip inside the material, by finding the correction from results of FE simulations of the nanoindentation experiment with different material model parameter sets. This process is called the method of snapshots [138]. The snapshot matrix, U then consists of the resulting indenter displacement that are expected to be somewhat correlated.

3.3.2 Radial Basis Function Theory

Radial Basis Functions (RBF) are very effective in providing an output approximation of a multivariable function for an unknown input point in the parameter space through interpolation of information from the known points [139]. In this section a very brief description of RBF is provided. The procedure through which RBFs can be combined with the information from POD to solve the inverse problem is also be discussed in the following paragraphs.

As mentioned earlier RBF is a very effective interpolation technique. To illustrate the idea of RBF, let us assume a function f(x) for which we only know N number of input–output relations. Let us also assume, x is a point in the parameter space for which we want to approximate the function’s value, where x is a M-dimensional vector. Classical interpolation methods use only the information around the point x to provide the approximation. The biggest difference that RBF provides in a similar scenario is that it uses all the N number of input–output relationship to build one continuous function over the whole domain. Therefore, the actual function f(x) is approximated as a linear combination of some function


where are coefficients of this combination. This equation is complete when the basis functions and the coefficients are known. Various Radial Basis Functions can be chosen as basis function . Most notable few are given below—


For an unknown point x in parameter space, the linear spline RBF will provide the basis using the following manner


For determining the coefficients , known N values of the function in the nodes are used in such a manner so that the RBFs approximate exact values of the function at the known points. This is solved using the following equation


where are the known values of the function. In compact manner Eq. 3.0 can be written as



Eq. 3.0 can be solved for unknown interpolation coefficients , which can then be used to obtain approximated values of the function in any given points in the parameter space. For a particular sampling set N, is only need to be determined once. In matrix notation Eq. 3.0 can be written as


As RBF takes into account the whole set of input–output relationship of a system, it can provide much more informed approximation compared to the classical local interpolation schemes. Another important advantage of using RBF is that the sampling of N in the parameter space need not to follow any particular distribution (in other words, can be scattered). However, particular distribution of sampling points help to keep the error of approximation under control.

3.3.3 Combining POD–RBF for Approximation

The ability of POD is to create a reduced order model by truncating orthogonal basis or dimensions. In a manner, POD works as a data compression tool where the loss of data is negligible. On the other hand, RBF provides the ability to approximate a function with high fidelity in between the known values in a multivariable parameter space. If both techniques are combined we can get a tool that can essentially provide high quality output approximation without incurring the huge computational cost associated with finite element analysis during an inverse analysis.

In context of nanoindentation study, let us assume vector p collects the material model parameters and u collects the output of the simulation (load or displacement data). Our goal is to find a function such that f(p) = u. This function needs to approximate the output of the simulation over some domain in parameter space. Following the theory related to POD, a reduced dimension model of snapshot matrix, U can be developed where represents the reduced order of amplitudes. In reduced dimension, the aforementioned equation can be written as


where, the relation between the reduced model and full model is given by the following equation


If RBF is applied Eq. 3.0 can be expressed in following manner


Once the basis function is known, interpolation coefficients collected in matrix B is solved in the reduced space using the following equation


Then the final equation that will provide the approximation of the system response for any arbitrary set of parameter in the subspace is given by


Eq. 3.0 involves simple matrix multiplication, and thus can provide much faster turnaround time when compared to finite element simulation. This is particularly useful for inverse analysis where a large part of computational effort is directed towards running simulation inside the optimization loop. It is also a simple approach, where the training of the POD–RBF (obtaining the matrices and B) is done only once. Moreover, once trained this technique can provide high enough computation accuracy, which can even be improved with a larger sampling points.

3.4 Taguchi Design of Experiments for Sensitivity Analysis

In any process where the output is influenced by multiple number of parameters, there is a need for the information that how individual parameters affect the overall output. In other words, it is useful to know the sensitivity of an output to an input parameter change. This need gave rise to an independent area of research inside statistics called Design of Experiments (DOE).

Traditionally, researchers used to carry out experiments where only one of the parameter was changed within a certain range while keeping the other parameters constant. Then the same process was to be replicated for the other parameters. This method is called full factorial experimental design, where the number of experiments required to perform the sensitivity analysis is astronomical. On the contrary, Taguchi applied the concept of orthogonal arrays, where all factors are changed simultaneously. For an experiment involving three parameters changing in four levels, the number of experiments required by Taguchi method is only 16, while full factorial design requires 64 independent experiments.

To perform a systematic sensitivity analysis first an experimental design is required. It is done by choosing an orthogonal array depending on the degrees of freedom (Eq. 3.0):


If is the number of levels for factor A, then = – 1. The experiments are conducted based on the chosen orthogonal array. By employing suitable analysis technique, such as Analysis of Variance (ANOVA) one can determine the contribution of individual parameters contribution towards the output. ANOVA is an useful statistical tool for quantitative determination of influence of any given input parameter and it can be used to interpret experimental data.

Chapter 4 Application of the Proposed Technique

The main objective of this study was to use an inverse analysis technique to extract material model parameters for a nonlinear viscoelastic model from a nanoindentation experiment. There were three big challenges to this problem —

  1. Modeling nanoindentation experiment using finite element analysis

  2. Incorporating nonlinear viscoelastic model in the finite element simulation

  3. Developing the optimization routine to extract the model parameters

In last chapter modeling of the nanoindentation experiment for a Berkovich tip using commercial software package ABAQUS was described. To verify that the ABAQUS model was in fact able to simulate the nanoindentation experiment, a simple elastic indentation simulation was performed. From the simulation corresponding load and displacement data were obtained, which was compared against Hertzian analytical solution provided by Sneddon. Figure 4.1 shows the comparison of analytical and simulated load–displacement data. It can be seen that, although it was not a perfect match, simulated data closely followed the analytical data. Attaining perfect match between simulation and analytical data is not very practical as it means going for very fine meshing in the simulation model thus increasing the computational expense exponentially.

The second challenge which was to incorporate nonlinear viscoelastic model in the finite element simulation has also been solved. This required discretizing the constitutive mathematical model for time step, . A detailed description of the numerical approach that was involved in developing the equations required for finite element approach has been presented in the previous chapter.

Figure 4.1: Load–displacement behavior of a Berkovich nanoindentation

The third challenge, which was the most critical of the three, has been solved by combining two separate technique, i.e. Proper Orthogonal Decomposition and Radial Basis Functions. The nonlinear viscoelastic model of our choice has seven unknown parameters of interest that needs to be extracted using the POD–RBF technique. As discussed in the previous chapter, POD–RBF technique needs snapshots of the system to become trained in approximating the system.

The process of training a surrogate model is often referred as sampling. A simple way of sampling the parameter space can be the grid system, where the distance between the sampling points for a parameter is kept constant over the domain of interest. If every unknown parameter is sampled m times over its domain, for n number of unknown parameters a total of simulations will be required. This would be a large number of finite element simulations to handle. Hence, in order to verify the POD–RBF technique’s ability in solving an inverse problem, first a rather simple problem where the number of unknown parameters are less were chosen. This was done in order to investigate the key properties of a POD–RBF based surrogate model, which could later be utilized to increase the confidence in POD–RBF technique for the ultimate application i.e. parameter identification of nonlinear Burger’s model.

The performance or the ability of the POD–RBF based surrogate model to precisely approximate the FE simulation depends on couple of parameters, namely number of training points used and the choice of RBF. Higher number of training points relating input parameters to system output improves the quality of surrogate model’s approximation at the expense of higher number of FE simulations. Although being offline or outside optimization loop, optimizing the number of training points is crucial since it directly relates to the computational cost of the overall inverse analysis.

There are only handful of articles in the literature that have tackled the inverse problem of nanoindentation-based material model calibration using POD–RBF based surrogate models  [140, 77, 141, 142, 143, 144, 145, 146]. Furthermore, to our best of our knowledge, none of the previous studies reported if the performance of the surrogate model could be optimized with respect of number of training points. In addition, the choice of an RBF, which affects the performance of the surrogate model, has also not been investigated at depth. Prior studies have typically employed only one kind of RBF in an ad hoc manner without providing much analysis into comparative benefits of using different types of RBFs to solve a given problem.

Since a well-trained surrogate model is at the root of solving the nanoindentation-based inverse problem, this study was designed to facilitate the understanding of a POD–RBF based surrogate model’s performance with respect to the number of training points and the choice of an RBF. It was expected that the findings of this study would provide a general framework for solving nanoindentation-based material modeling inverse problem using POD–RBF technique.

4.1 Case Study

In this study, nanoindentation was conducted on a standard metallic material. The nanoindentation experiment was then modeled with a finite element analysis software, where a custom elastic–plastic material behavior was incorporated. A range was selected for each parameter within which the values of the parameter would be altered. A Taguchi orthogonal array-based experimental design was formulated by varying each parameter within the range at four equidistant levels. The analysis of variance (ANOVA) technique was employed to recognize the influence of the parameters over the output. The number of levels for the unknown parameters within the specified range were optimized based on the ANOVA results. Training data were generated in a full factorial basis by varying each parameter of the custom material model for the initial and optimized model. A random noise of 1% and 5% was appended to the training data to investigate the stability of each surrogate model.

4.1.1 Nanoindentation

The nanoindentation experiment was conducted using an MTS Nanoindenter XP (Agilent Technologies, Santa Clara, CA, USA) using a load-controlled scheme with a Berkovich tip. The maximum load was set to be 4.9 mN for the experiments. A triangular loading profile was chosen with a 15 s duration for both the loading and unloading segments. Before conducting the actual experiments the Berkovich tip was calibrated using a fused silica reference material. Also, the acceptable thermal drift rate was chosen to be 0.15 nm/s.

The nanoindentation experiment was conducted on a reference material i.e. single crystal aluminum. This sample is commonly used to check the performance of a nanoindenter. The single crystal aluminum sample has Young’s modulus of 70.4 GPa and Possion’s ratio of 0.345 as provided by the supplier (Agilent Technologies, Santa Clara, CA, USA).

4.1.2 Finite Element Simulation

A commercial finite element software (ABAQUS, Dassault Systémes, Providence, RI, USA) was utilized in this study, both for modeling the nanoindentation experiment and for solving the finite element 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.

The contact between the indenter and the sample was assumed sliding contact with a friction coefficient of 0.25 and was defined as surface-to-surface contact. The indenter and the sample were assigned as the master and the slave surface, respectively. In ABAQUS surface-to-surface contact, master surface nodes can penetrate the slave surface (i.e. causing deformation to the slave surface), while the slave surface nodes cannot penetrate master surface. In indentation modeling using FE, it is generally assumed that the indenter is much stiffer than the sample surface; hence, deformation of the indenter by the sample surface is neglected. The element type was chosen from the eight-node brick element family (C3D8). The finite element problem consisted a total of 1323 elements and 1817 nodes. Figure 4.2 shows a schematic of the FE model generated in ABAQUS and Fig. 4.3 shows the typical ABAQUS simulation’s stress contour outputs at the end of a Berkovich indentation.

(a) Stress–strain relationship of bilinear plasticity model
(b) Schematic of ABAQUS finite element model for Berkovich indentation
Figure 4.2: ABAQUS finite element modeling details
(a) Von Mises equivalent stress
(b) X-axis stress component,
(c) Y-axis stress component,
(d) Z-axis stress component,
Figure 4.3: Stress contours ABAQUS output after unloading for elastic–plastic Berkovich indentation simulation

The elastic–plastic material behavior of the sample was incorporated in the FE software using a UMAT subroutine. The material behavior was chosen as isotropic elastic–plastic with linear hardening as shown in Fig. 4.2 (a). This model defines both the elastic and the plastic part of the stress–strain relationship as linear [147, 148]. Only four parameters are required to describe this particular material model, which are elastic modulus (E), Poisson’s ratio (), yield strength (), and hardening coefficient (h). For the numerical study, Poission’s ratio was kept constant at the known value of 0.345. According to prior FE-based studies, Poisson’s ratio does not affect the FE simulation of indentation experiment as much as the other model parameters [149] and hence, kept constant in most of the indentation modeling studies [150, 151]. Table 4.1 lists the range of values used in this study for the three parameters.

Model parameter Lowest level Highest level
Young’s modulus, E 60 75
Yield strength, 0.05 0.20
Hardening coefficient, h 0.4 0.7
Table 4.1: Range of values for the model parameters

4.1.3 Taguchi Design of Experiments and Sensitivity Analysis

When a system’s output is governed by two or more independent variables, information about each variable’s influence over the output may provide deeper insight into the optimization problem. In other words, it is important to know how the system output is affected by each given input parameter. By doing so the performance of the overall optimization routine could be greatly improved since this information could subsequently be utilized in reducing the computation expense of the meta-model development.

In this study, a Taguchi-based design of experiments methodology with ANOVA was adopted to quantify each input parameters contribution towards the overall output or the error function. Employing Taguchi orthogonal arrays instead of a full factorial experimental design help in reducing the number of finite element simulations required in assessing the sensitivity of model parameters.

The first step of performing a sensitivity analysis is to define an experimental design, which involved choosing an appropriate orthogonal array. This was achieved by first calculating the ‘degrees of freedom’ of the experiment, as


where, , is designated as the number of levels for the input parameter , and are from the interaction between model parameters. In this study four levels were considered for each of the three input parameters, as listed in Table 4.2. As a result, the degrees of freedom for each factor equaled 3. No interaction was considered among the model parameters. Hence, the total number of degrees of freedom for the experiment was found to be 9. The Taguchi orthogonal array which can successfully accommodate this number of degrees of freedom is modified . The experimental design for this study according to the modified orthogonal array is listed in Table 4.3.

Levels Elastic
E (GPa) Yield
(GPa) Hardening
Level 1 60 0.05 0.40
Level 2 65 0.10 0.50
Level 3 70 0.15 0.60
Level 4 75 0.20 0.70
Table 4.2: Levels of material model parameters
Experiment Elastic modulus Yield strength Hardening coefficient
1 1 1 1
2 1 2 2
3 1 3 3
4 1 4 4
5 2 1 2
6 2 2 1
7 2 3 4
8 2 4 3
9 3 1 3
10 3 2 4
11 3 3 1
12 3 4 2
13 4 1 4
14 4 2 3
15 4 3 2
16 4 4 1
Table 4.3: Experimental design based on the modified L orthogonal array

For each of these experiments, finite element simulation yielded results in terms of indenter displacement as a function of indentation load. The load increments for the simulation was chosen in such a manner so that it matched with the experimental loading data. The error function, δ for this study was defined by the following equation.


where, is the number of data points in the load–displacement plot. By following the Taguchi orthogonal array experimental design the relationship of three model parameters with the system output or the error function was formulated, which was then analyzed using analysis of variance (ANOVA).

4.1.4 POD–RBF Based Surrogate Model

The proper orthogonal decomposition (POD) theory, also known as principal component analysis (PCA), was developed to approximate a function over some domain of interest based on the known relationships between the input and the output [137, 136, 127]. This study followed the POD–RBF procedure outlined by Rogers et al. [152].

As per POD terminology, the relationship between the input and the output of a particular system is called a snapshot. In the context of this study, snapshots or training points were relation of material model parameters and the output tip displacement data. If number of simulations were carried out where in each of them at least one input variable was changed, then the snapshot matrix was formulated by combining number of displacement vectors. Moreover, if the output of the simulation (displacement vector) had data points, then snapshot matrix can be defined as,


Input material model parameters were collected in the input matrix, . The first step towards creating a reduced order model using POD was to generate snapshots of the system for the range of input parameters and subsequently combining all these appropriately in the and matrix. A brief outline of surrogate model training using POD–RBF technique is provided here without detailed mathematical derivations, which can be found in the literature [153, 154].

Step 1:

Develop the covariance matrix for the snapshot matrix , where = . .

Step 2:

Find the POD orthonormal basis vectors (for ) which would optimally represent . Here, POD basis matrix = . , and represents the eigenvectors of . can be found by solving the eigenvalue problem noted as .

Step 3:

Truncate the POD basis based on the energy of the POD modes and determine . The subsequent POD model would retain majority of the information about the system, while reducing the dimension of the problem considerably. The truncated POD basis, = . .

Step 4:

Once truncated POD basis matrix is known, the amplitude matrix can be computed as, . is defined as a nonlinear function of matrix. At the time is known, POD reduced order model of the system is ready, and data can be interpolated to find out the surrogate approximation for unknown input parameters.

Step 5:

Compute the coefficient matrix as, , where, is the matrix of interpolation functions or RBFs in the context of this study. is defined as–


where, and are input parameter vectors used to generate the i-th and j-th snapshots, respectively.

Step 6:

At an unknown point, in the parameter space, the system output can be computed using the relationship,