Material model calibration through indentation test and stochastic inverse analysis

Indentation test is used with growing popularity for the characterization of various materials on different scales. Developed methods are combining the test with computer simulation and inverse analyses to assess material parameters entering into constitutive models. The outputs of such procedures are expressed as evaluation of sought parameters in deterministic sense, while for engineering practice it is desirable to assess also the uncertainty which affects the final estimates resulting from various sources of errors within the identification procedure. In this paper an experimental-numerical method is presented centered on inverse analysis build upon data collected from the indentation test in the form of force-penetration relationship (so-called indentation curve). Recursive simulations are made computationally economical by an a priori model reduction procedure. Resulting inverse problem is solved in a stochastic context using Monte Carlo simulations and non-sequential Extended Kalman filter. Obtained results are presented comparatively as for accuracy and computational efficiency.


Introduction
Indentation tests are employed since decades for mechanical characterization of materials at diverse scales and in various technological fields. The original form of the test was fairly simple, intended for the measurement of 'hardness' in metals. Later the assessment of mechanical parameters entering into constitutive models and necessary for overall inelastic structural analyses became the purpose of this test. The parameter identification at present is usually based on 'instrumented' indenters, namely those that are capable to provide a so-called 'indentation curve', a curve that correlates the indenter tip penetration versus force applied on it, see Figure 1.
Transition from indentation curves to the sought parameters in the early works was accomplished by semi-empirical formulae (see e.g. [1]). Recently the use of computational methods through 'inverse analysis' has been gaining popularity and is a subject of vast literature (see e.g. [2,3,4]). In this concept, the test is simulated, traditionally by the use of Finite Element Modeling (FEM), while the results of simulation are compared to those measured in the experiment in the form of a 'discrepancy function'. This function quantifies the difference between computed and measured data, and in a subsequent phase it is minimized with respect to sought parameters by a suitable mathematical programming technique.
Material parameter calibration based on inverse analysis offers improved accuracy and better flexibility with respect to semi-empirical formulae. This is particularly the case when uncertain or incomplete data, are used from the experiment.
In practical applications of material model calibration based on indentation test, this circumstance occurs quite frequently since there can be many factors affecting the measured indentation curve (e.g. rough surface of tested specimen, error in measured force and corresponding displacement, slight difference in local material properties etc, [5]).
In order to use measurable quantities endowed by a covariance matrix quantifying their uncertainties, it is required to design an inverse analysis procedure capable to estimate the propagation of this uncertainty up to the assessed parameters. To achieve such goal, a stochastic framework needs to be adopted for the resulting inverse problem.
In this paper Monte Carlo (MC) simulations and non-linear Kalman filters (KF) are employed to solve the parameter identification problem in a stochastic context, on the basis of indentation curves generated by an instrumented indenter. The indentation test referred to in this study considers conical indenter with 120 0 cone opening angle and a 200µm spherical radius on a tip, like in traditional Rockwell tests [6]. The material specimen response is assumed to be isotropic elastic-plastic with linear hardening, following Drucker-Prager yield criterion in a view of potential application of proposed approach to the characterization of metal-ceramic composites.
The developed procedure is verified within pseudo-experimental context (i.e. computer generated data are used instead of those collected from a real experiment) in order to comparatively check different strategies in terms of accuracy and computational efficiency in the present context. Computational burden connected to recursive simulations of an indentation test, particularly required by Monte Carlo simulations, is significantly reduced by the employment of a "model reduction" technique based on Proper Orthogonal Decomposition (POD) apt for the acceleration of non-linear simulations (see [7] for details).
The paper is organized as follows: Section 2 gives a brief outline of inverse analysis by mathematical programming and model reduction technique for fast and economical parameter identification. Section 3 is devoted to the stochastic approaches described in the form implemented in present context. A reference problem with several specific modeling features adopted in this study is given in Section 4, followed by comparative presentation of achieved numerical results. Conclusions and future research prospects considering discussed, and related items are outlined in Section 5.

Inverse analysis with mathematical programming and model reduction
Material model calibration based on inverse analysis relies on performing of an experiment (suitably chosen to provide measurable quantities sensitive to the sought parameters) and its subsequent simulation. Results gathered from the experiment are used together with their computed counterpart in order to form a 'discrepancy function' that quantifies the difference between the two. Resulting inverse analysis is centered on minimization of this function. Such problem can be briefly formulated as follows.
Let u e be the vector of experimentally measured data, and u(p) the computed counterpart resulting from a simulation, with sought parameters collected in vector p.
The discrepancy function in a purely deterministic formulation in its simplest form is given by: The solution to the formulated inverse problem is defined by a set of parameters collected in vector p SOL that is minimizing the above defined discrepancy function. This minimization problem is solved numerically as there is no analytical description of the objective function (1). To this purpose one can make a recourse to 'first order' algorithms which are relying on a quadratic approximation of the objective function.
Example of these are methods based on 'Trust Region Algorithm' (TRA). Description of various TRA types is widely available in the literature (see e.g. [8,9]) and a detailed presentation is omitted here for brevity. However, a brief outline of one of the most popular versions of TRA, used also in this study, is given in what follows for the sake of completeness.
In each iteration of TRA a quadratic approximation of the objective function is generated by an approximation of the Hessian matrix computed through Jacobian of the 'residual' vector, namely a vector that collects differences between measurable and computed quantities. Such quadratic approximation of the objective function represents a model function that is 'trusted' to be an adequate representation of the real function within defined region. Minimization of this model function is computationally inexpensive and its minimizer is searched for within a two-dimensional subspace Remarkable practical advantages can be achieved by invoking a model reduction technique based on Proper Orthogonal Decomposition in order to accelerate simulations. The POD procedure has been already adopted in researches and its description can be found in [10,11]. Here just a brief outline of stages involved in it will be given.
Within 'search domain' in the sought parameter space N points are selected that are uniformly covering this domain. Each one of this N 'nodes' is assumed as input to the simulation of test leading to computed outputs collected in vector u i . Exploiting The above model reduction procedure makes the parameter identification fast and economical, as the recursive simulations required by the optimization algorithm are performed using 'light' computational tool.
Using numerical optimization algorithms like TRA the output of resulting inverse problem is a vector of parameters with their estimated values. In those cases when experimental data are endowed by corresponding covariance matrix, denoted by C EXP , quantifying the experimental error of measured quantities, the discrepancy function can be written in weighted residual form such as: This modification is important in the situation where larger noise in measurements is present in certain zones (e.g. starting part of the indentation curve due to asperities of the unpolished surface). In this form of discrepancy function merely more importance is conferred to more accurate measurements, but the context remains deterministic, as the output of the minimization is still a single vector of assessed parameters without quantification of the uncertainty affecting them. In engineering applications however, it is highly desirable to estimate also the uncertainty of the identified parameters. To reach such result a stochastic procedure needs to be employed for solving the inverse problem. In the following section a brief outline of stochastic methods employed in this study is presented.

Stochastic inverse analysis for material parameter calibration
When uncertain measurable quantities, collected from the test, are used within inverse analysis procedure for material parameter calibration, it is important to quantify to which extent these uncertainties are affecting the outputs (namely sought parameters).
To this purpose a stochastic procedure needs to be employed in order to solve the inverse problem. In this study two particular methodologies are used in order to reach such result: Monte Carlo sampling and Kalman filters. In what follows an outline of the application of these methods in the present context of material characterization based on indentation test is given.

Monte Carlo sampling
Monte Carlo (MC) methods are quite popular stochastic methods and the simplest to implement (see e.g. [12]). Basic concept of MC methods is a generation of set of outputs corresponding to a given set of inputs. It is further possible to apply sampling methods to estimate statistical nature of the outputs for a given variety on inputs (that can also be statistically sampled). The main advantage of MC methods is that it is sufficient to solve the deterministic model. Major limitation concerns the convergence rate, as the number of pairs of inputs and outputs needs to be relatively high in order to represent a statistically meaningful sample. Generation of this large samples may be computationally expensive.
Applying MC methods to inverse analysis is computationally even more costly, It is worth noting that the computation burden connected to S inverse analysis (with S equal to 2000 or more) is significantly reduced by the use of reduced basis model specifically calibrated for the indentation test referred to in this study.

Extended Kalman filters
Kalman filtering is developed in control engineering as a tool that provides a recursive solution to the linear optimal filtering problem for dynamic systems [13]. The developed algorithm can be used to produce estimates of unknown state variables with higher accuracy relying on series of measurements observed over time, instead of using a single measurement.
By applying Kalman filters to solve inverse problems it is possible to generate the covariance matrix of estimated parameters, starting from measured quantities endowed by quantified uncertainty. The original form of Kalman filters is intended to be used for linear systems, and hence cannot be directly applied for the problem at hand, as the inverse analysis based on indentation test is highly non-linear problem. Thus, extended versions of Kalman filters will be employed, in its modified form to take into account non-sequential experimental data. The governing equations used for the updating procedure are derived by linearization about current estimate.
Theoretical fundamentals and detailed derivations considering Kalman filtering methodology can be found in [14,15]. In what follows a sequence of operations with proper updating equations adopted in this study is outlined.
Let d EXP be a vector of experimentally measured quantities, here representing R pairs of points (i.e. force and corresponding penetration depths) from the pseudoexperimental indentation curve, while C EXP is an R×R covariance matrix referring to the uncertainty in the measurements. Iterative procedure used for the estimates in the parameters and corresponding covariance matrix is following steps listed below.
(1) A starting value is attributed to all P sought parameters, collected in vector p GUESS , and a vector of residuals is computed using direct operator (i.e. numerical model for test simulation) according to

GUESS EXP COM GUESS
(2) Sensitivity matrix denoted by L is computed in the following form: resulting matrix K is of the order P×R.
(4) Updating equation is used for the evaluation of 'new' parameter values using experimentally measured data, and computed gain matrix by the following expression: Expression in the square brackets represents already computed vector of residuals, so no further computational 'cost' is required for the calculation of (6) (5) A R×R covariance matrix of estimated parameters is evaluated using the equation given by: ( ) Presented updating procedure results from a full derivation by attributing to the initial covariance matrix of the parameters a large value, causing its inverse to vanish.
The presented scheme therefore can be used in the absence of the a priori knowledge of the uncertainty affecting initial 'guess' value of the parameters. In case when also the covariance matrix of initial parameters is specified at the beginning of identification, equation (7) is modified so that also this matrix is updated throughout the procedure (see [16] for details). In this study the most general procedure requiring least amount of data as inputs is tested.
Iterative procedure given by equations (3) to (7) is terminated by reaching convergence criteria, previously specified on minimal changes in parameters between two consecutive iterations. Result of this calculation is vector of resulting parameter p RES endowed by appropriate covariance matrix C P , with elements on the main diagonal being a square roots of standard deviation of assessed parameters (i.e. The isotropy of the material to be tested suggests axially symmetrical model. Finite Element Model (FEM) adopted for the test simulation by the commercial code ABAQUS [17] is visualized in Figure 2. It exhibits the following features: unilateral contact between surface of a rigid indenter and modeled specimen surface with Coulomb friction coefficient equal to 0.05 (considered as known parameter here, alternatively can be identified together with other sought parameters, see e.g. [18]); large strain regime according to their implementation in ABAQUS [17], not discussed here; zero displacements on the boundaries common with the surrounding practically undisturbed zone (adopted here for simplicity instead of a more accurate alternative as statical condensation of the surrounding solid); 4365 four-node quadrilateral elements (hence 8878 degrees of freedom).
Material response of tested specimen is assumed to be elastic with plasticity following Drucker-Prager (DP) yield criterion, in a view of potential applications of tested procedures to metal-ceramic composite materials. The use of more general DP formulation instead of Hencky-Huber-Mises (HHM), which is included in former as a special case of it, is suited better for describing mechanical behavior of ceramic materials and composites with higher content of it (see [19], [20]). DP yield criterion therefore can be seen as a modification of HHM criterion in order to take into account pressure sensitivity. Yield surface according to DP criterion with linear hardening has the following form: where s is deviatoric stress tensor and J 2 its second invariant, I 1 is the first invariant of stress tensor σ σ σ σ, α is the internal friction angle, k is the cohesion, while the last term represent contribution of linear strain hardening, with h being a material parameter and λ plastic strain-rate multiplier. By assuming associative plasticity (and therefore adopting that dilatancy angle is equal to internal friction angle) the number of constitutive parameters that govern plastic response following DP criterion is equal to 3.
The internal friction angle α and cohesion k depend on initial tensile and compressive yield limits, 0 t σ and 0 c σ as follows: For the material under investigation in this study constitutive parameters to assess by designed inverse analysis procedure, are Young's modulus (E), and additional three plastic parameters: α, 0 c σ and h, while remaining parameters can be calculated using equation (9).
In order to computationally verify outlined procedures, a pseudo-experimental indentation curve is generated using designed FEM model attributing to the material parameters the following reference values: E=80GPa, , α=10 0 and h=1000MPa. These values can be regarded as typical and representative for ceramicmetal composites ( [19], [21]), and will be treated as "targets" for the further identification procedure. Maximal adopted indentation force is taken to be equal to 500N, which is a value large enough for this type of materials, to ensure global response of the tested specimen.

Stochastic inverse analysis by Monte Carlo sampling
A reference indentation curve is generated through test simulation attributing to parameters above listed values. In order to produce additional curves starting from a reference one, a procedure outlined in section 3.1 is used. Here, in particular, the following level of noise was used: nf=±2%, nd=±5%, with M=15 different random noise extractions resulting in 15 different curves. These curves are further sampled in order to obtain a mean curve and a corresponding standard deviation for each of 50 levels of force along loading part and 50 levels of force along unloading part.
In a subsequent step, S=2500 different indentation curves are generated using Gauss distribution based on previously calculated statistical values. For each of these curves an inverse analysis is solved cantered on minimization of a discrepancy function defined by function (1). The solution of this minimization problem is here achieved by means of traditional mathematical programming algorithm, an iterative "first-order" Trust Region Algorithm, in its version available in MATALB [22], and briefly outlined in section 2. Recursive test simulations are performed by a fast computational tool based on POD and RBF interpolation, summarized in section 2, and described in details in [7], [11]. Designed Reduced Basis (RB) model is "trained" once-for-all by a set of 2520 simulations confined by the following upper and lower bounds of sought parameters, defining a reasonable range for the type of materials referred to in this study, namely: The use of RB model significantly reduced computational burden involved in solving 2500 inverse problems, as it provides results to the test simulation in a realtime, with controllable error that is kept on the level as of FEM calculation. Resulting parameters from these inverse analyses are sampled in order to obtain mean value and corresponding standard deviation and quantitative values of these are given in Table 1.
By observing results in the table in may be noticed that the mean value of all the assessed parameters is rather close to the target one with the largest error less than 1%, regardless, rather significant "noise" level used as perturbation of the reference curve.
This corroborates the fact that inverse problem is well-posed, and the stochastic solution of it provides results with limited influence of level of perturbation on the mean values of estimates. On the other hand standard deviation of all estimated parameters are affected to a larger extent by the perturbed inputs. Still however, they may be regarded as reasonably small (i.e. all comparable to considered level of noise which was ±5%) proving stability of the proposed stochastic solution.

Extended Kalman filters for stochastic solution of inverse problem
Present identification problem is further solved by an updating procedure of extended Kalman filters, outlined in sub-section 3.2. Reference curve, used as input, is accompanied also by a covariance matrix C EXP constructed using the same noise level as in the previous case. The matrix is calculated by attributing to all diagonal members values that correspond to standard deviation of particular level of force. Indentation curve here is divided into 50 points along loading part and 50 points along unloading part, therefore the resulting matrix C EXP is 100×100 with diagonal elements equal to: where S Fi is previously calculated standard deviation of i th level of force from the indentation curve. Stochastic solution of the inverse problem by KF leads to results that are summarized in Table 2.
Adopted updating procedure turned out to be quite efficient leading to mean values that are rather close to "target" one, with the largest error smaller than 1%, like in MC case. As pointed out in sub-section 3. Values of standard deviation for all parameters are higher than those obtained by Monte Carlo sampling. This is due to the selected updating procedure in which initial covariance matrix for parameters is assumed as unknown, attributing to the 0 P C large value, causing its inverse to vanish. Adopting this scheme inevitably led to larger values of assessed parameter's standard deviation. An improvement of this outcome can be reached by using different Kalman filter updating scheme which takes into account initial value of parameter covariance matrix.
It is worth commenting that somewhat larger values of standard deviation obtained for friction angle and hardening coefficient are due to the limited influence of these parameters to the test outcome (i.e. indentation curve). Sensitivity analysis conducted in this study showed that by perturbing separately each of four parameters with same relative amount with respect to their reference values, more pronounced effect on the indentation curve is provoked with perturbation of Young's modulus and yield limit than with remaining two parameters. Such circumstance led to somewhat larger values of standard deviation for these parameters, with mean values still accurately identified corroborating the robustness of the identification procedure.

Conclusions and closing remarks
Indentation instruments are presently equipped with tools capable to provide the relationship between applied force and the consequent penetration depth. The output of such test is a so-called 'indentation curve' which can be used for characterization of diverse materials, see e.g. [1,2,3]. These characterization procedures provide constitutive parameter estimates only. In engineering practice it is advantageous to assess also the uncertainty affecting the final estimates as a consequence of various sources of errors embedded within identification procedure. Updating procedure based upon the assumption of the absence of reliable quantification of 'guess' parameter uncertainty, described by a large initial covariance, obviously leads to somewhat larger values of standard deviation of estimated parameters, when compared to values obtained by MC, here taken as reference.
Future promising developments in progress concern the use of conceptually similar updating procedure that takes into account initial parameter guess covariance matrix. Further on an alternative, so called 'unscented' Kalman filter [23], [24], that doesn't require Taylor's expansion to approximate the non-linear forward operator, and its potential application in the present context is also subjection of current research.