Proper Orthogonal Decomposition and Radial Basis Functions in material characterization based on instrumented indentation

Engineering Structures 33 (2011) 492–501

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Proper Orthogonal Decomposition and Radial Basis Functions in material characterization based on instrumented indentation Vladimir Buljak, Giulio Maier ∗ Department of Structural Engineering, Politecnico di Milano (Technical University), piazza Leondardo da Vinci 32, 20133, Italy

article

info

Article history: Received 20 November 2009 Received in revised form 26 October 2010 Accepted 2 November 2010 Available online 30 November 2010 Keywords: Instrumented indentation Inverse analysis Proper Orthogonal Decomposition Radial Basis Functions Material characterization Structural diagnosis

abstract For the mechanical characterization of structural materials non-destructive tests combined with computer simulations and inverse analyses are more and more frequently and advantageously employed in various engineering fields. The contribution to such development presented in this paper can be outlined as follows. With reference to isotropic elastic–plastic material models, indentation test simulations are done preliminarily, once-for-all, by a conventional finite element forward operator. Results of these simulations are employed in a procedure which is centered on Proper Orthogonal Decomposition and Radial Basis Functions approximation and is used for fast interpolations which replace further finite element analyses in the parameter identification process. Comparative computational exercises are presented in order to point out the consequent significant reduction of computing times in test simulations and, hence, in the minimization of the discrepancy function by the Trust Region Algorithm, namely by a traditional first-order mathematical programming method. Such a parameter identification procedure may be carried out routinely and economically on small computers for in situ structural diagnoses. Both the force–penetration relationship (provided by an instrumented indenter) and the average imprint profile (achievable by laser profilometer) are considered as sources of measurable response quantities or experimental data. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction The assessments of inelastic properties of materials are usually called ‘‘destructive’’ if laboratory tests are required on specimens extracted from structures. Indentation tests, which are becoming more and more frequent in engineering for structural diagnosis (and hardness tests which represent their historical origin) can be considered ‘‘non-destructive’’ (or at least ‘‘almostnondestructive’’). In fact, they usually imply neither service interruption nor a reduction of safety margins of the structure or plant involved. At micro- and nano-scales such an advantage is not exhibited by indentation but it is usually not practically important (e.g. for MEME and NEMS not considered here). Truly ‘‘nondestructive’’ experiments, like ultra-sonic tests, cannot provide data on material inelastic properties, a meaningful limitation from a structural mechanics standpoint. Inverse analyses based on experimental data collected from instrumented indentation tests and on their simulations have been widely used in recent years for material parameter identification in structural engineering. Several developments of such methodology

∗

Corresponding author. Tel.: +39 02 2399 4221; fax: +39 02 2399 4220. E-mail addresses: [email protected] (V. Buljak), [email protected], [email protected] (G. Maier). 0141-0296/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2010.11.006

oriented to different industrial applications can be found at present in the literature. For example, mechanical properties of functionally graded materials and of thin films have been assessed by such an approach in Refs. [1,2], respectively. In [3] the exploitation of two kinds of experimental data was proposed and shown to be useful, namely of both the load-versuspenetration digitalized curves provided by the instrumented indenter and the residual imprint geometry measured by a laser profilometer as an additional instrument. Possible perturbations on the imprint geometry due to the inhomogeneity of material microstructures are made irrelevant by averaging imprint profiles in various directions and by adopting a penetration much larger than the typical size of that microstructure. By means of such an amplification of available experimental data, nondestructive indentation tests, originally proposed to assess hardness (see e.g. [4,5]), at present can be successfully employed if associated to test simulation and inverse analysis, in order to assess also anisotropic material properties [6], tensorial residual stresses [7] and quasi-brittle fracture properties [8]. Clearly, the above mentioned properties cannot be assessed on the basis of conventional data provided by the instrumented indenter alone. However, even if directionindependent quantities are to be assessed in the indentation specimen (like parameters in isotropic elasto–plastic constitutive models to be considered herein), two, instead of one, sources of experimental data contribute to the ‘‘regularization’’ in the Tikhonov

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

sense of the inverse problem (namely, to the convexity of the discrepancy function to be minimized), see e.g. [9]. However, even in the presence of well-posedness, the numerical solution of inverse problems, both by traditional and modern procedures, turns out to be computationally heavy. In fact, it generally implies a sequence of direct analyses, namely of computer simulations of the test with diverse inputs of sought parameters. Such circumstances frequently represent a severe burden and handicap in engineering, when parameter identification becomes the main tool for the diagnostic search of possible structural material damages and should be done repeatedly and routinely (hopefully, if possible, in situ rather than in a computing center). The above practical difficulty may be drastically mitigated or avoided by the procedure presented and numerically validated in this paper. A simple popular constitutive model is attributed to the (metallic, ductile) material considered for mechanical characterization (e.g. required by damage diagnosis on pipelines and power-plant components, which specifically have motivated the present study). The identification of possible deteriorated material parameters is pursued by a conventional deterministic approach, namely by minimization of a norm (called ‘‘discrepancy function’’) of the differences between measurable quantities and their counterparts computed as functions of the unknown parameters. Test simulations, namely ‘‘direct analyses’’, are performed by the commercial Finite Element (FE) code ABAQUS [10]. Both the loading–unloading force-versus-penetration relationship (denominated here ‘‘indentation curves’’) and the average profile of the residual imprint measured by a laser profilometer (‘‘imprint profile’’) are considered as sources of experimental data (or ‘‘pseudoexperimental’’, computer-generated data), as proposed in [3]. The constrained minimization of the discrepancy function is carried out by the traditional mathematical programming algorithm called ‘‘Trust Region’’ (TRA), see e.g. [11], which involves first derivatives only (computed by finite differences, of course) and is available in MATLAB [12]. In order to drastically reduce the computing effort for the sequence of simulations required by TRA, the following operative procedure is adopted and investigated herein (basically similar to the one recently developed and applied to thermal problems in [13,14]): (i) construction of a Proper Orthogonal Decomposition (POD) basis of system responses (‘‘snapshots’’) in terms of both indentation curves and imprint profiles; these ‘‘snapshots’’ (namely vectors of measurable quantities) are computed ‘‘a-priori’’, once-for-all, on the basis of ‘‘nodes’’ of a suitable grid selected in the space of the material parameters looked for; (ii) ‘‘truncation’’ as low-order approximation of the POD basis and relevant snapshot ‘‘amplitudes’’ (see e.g. [15]); (iii) replacement of FE direct analyses within each TRA iteration through computationally fast interpolation by means of Radial Basis Functions (RBF) (see e.g. [16,17]). Section 2 is devoted to a brief survey of POD and RBF in view of their specific employment as mathematical tools within the present structural diagnostic analysis. The various comparative computational exercises in Section 3 on direct problems, and Section 4 on inverse problems are intended to evidence potentialities and limitations of the proposed procedure in structural engineering. Section 5 is devoted to conclusions and future prospects. The usual notation of matrix algebra is adopted (e.g. bold-face letters for matrices and vectors). 2. On Proper Orthogonal Decomposition (POD) and Radial Basis Functions (RBF) in the present context

493

most real-life practical problems: (a) in the space of the parameters to identify, a domain where to confine the search can ‘‘a priori’’ be identified by physical considerations and/or by an ‘‘expert’’; (b) if M tests are simulated by assuming different parameter vectors pi (i = 1, . . . , M) included in the domain mentioned in (a), the N quantities which are measurable within responses to the tests can be gathered in a N-vector u (called a ‘‘snapshot’’ in the POD jargon); this vector u is ‘‘correlated’’ to all other snapshots included in the set of M vectors generated through the above simulations, since all the snapshots represent consequences due to the same external action in the same system, with only differences in material parameters. The correlation mentioned in (b) means mechanically that changes of parameters in the domain (a) give rise to moderate changes of snapshot and, hence, geometrically that the snapshots are ‘‘almost parallel’’ and equally oriented in their N-dimensional space. This geometrical interpretation naturally suggests to find a new reference system whose axes will be sorted in the descending order of norm of all M snapshot projections: namely, the first new axis will have in N-dimensional space the direction that maximizes a norm of projections of the snapshots, the second one will give the second largest norm etc. In this new basis a low-order approximation of high accuracy can be achieved by preserving first K < N prevailing components, i.e. by ‘‘truncating’’ the negligible ones. In the present context the above ‘‘modus operandi’’, which is central to POD conceptual kernel, implies the size reduction (or ‘‘compression’’) of the information contained in the snapshots and concerning the essential role played by the parameters in the material specimen response to the indentation test. 2.2. Proper Orthogonal Decomposition: an outline The mathematical procedure of POD adopted herein and briefly summarized below is consistent with the approach called Principal Component Analysis (PCA), described with details and proofs e.g. in [15]. Recurrent symbols in the sequel have the following meaning: the S × M matrix P gathers as columns M sets of S identifiable parameters (pi , i = 1, 2, . . . , M); the corresponding snapshots ui (i = 1, 2, . . . , M) are collected in the N × M matrix U; 8 = [ϕ1 , . . . , ϕN ] denotes the N × N (or N × M matrix in the case when M < N) basis matrix which defines in the Nspace of the measurable quantities the ‘‘optimal’’ reference system to be found according to the criterion suggested by the expected snapshot correlation as mentioned in Section 2.1; the N × M matrix A quantifies the ‘‘amplitudes’’ of the snapshots ui in the basis 8, namely the snapshot matrix U can be re-constructed when A is known by the relationship: U = 8 · A.

(1)

Clearly, if the basis 8 is known, in view of its orthonormality, the amplitude matrix A is computed easily, namely (I being the identity matrix):

8T · 8 = I,

A = 8T · U.

(2)

The central mathematical kernel, originally demonstrated in probability and information theories and dealt with in various publications (see e.g. [15,18]), consists of the sequential maximization of Euclidean norms of snapshot components, as mentioned in Section 2.1, whose result can be condensed into the formula which follows: −1/2

2.1. Preliminary remarks

ϕ¯ i = U · vi · λi

In the present indentation methodology for materials characterization the following two circumstances can be expected in

where λi is the ith positive eigenvalue and v represents the corresponding normalized eigenvector of the following matrix

,

(i = 1, . . . , M )

(3) i

494

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

(symmetric, positive semi-definite of order M): D = U · U. T

(4)

Orthonormality, Eq. (2), provides equality constraints in the sequential maximizations (with ϕi as variables), the solution of which by a Lagrangean approach leads to Eq. (3). The computational burden of POD turns out to be concentrated in solving the eigenvalue problem related to matrix D. Once eigenvalues and eigenvectors are computed, Eq. (3) is used to compute POD directions which are then sorted in the descending order of eigenvalues λi into matrix 8, that represents the POD basis for the snapshot matrix. Due to its optimality, the new coordinate system constructed as a POD basis makes possible a low-dimensional approximation of highest accuracy with respect to any other approximation of the same order. This low-dimensional approximation is computed ˆ , that is derived by in a ‘‘truncated’’ POD basis, denoted by 8 preserving the first K (K ≪ N) POD directions (i.e. columns of 8) that correspond to the largest eigenvalues. Consequently, a simple matrix multiplication straightforwardly generates, according to ˆ consistent with the new Eq. (2), the K × M matrix of amplitudes A truncated POD basis. Finally, a low-dimensional approximation of the snapshot matrix, or of a single snapshot, can be expressed as a linear transformation of amplitudes, namely:

ˆ · Aˆ , U≈8

ˆ · aˆ i , or ui ≈ 8

(i = 1, . . . , M ).

Eq. (5), of the snapshots due to grid nodes pi , by modifying the amplitude vector as a function of parameters, namely:

ˆ · aˆ (p). u(p) ≈ 8

aˆ k (p) =

2.3. Interpolation by Radial Basis Functions Let the considered system (specimen or structure) be characterized, as for material properties, by a vector p of constitutive parameters different from those selected as grid nodes pi , (i = 1, . . . , M), in the parameter domain, but still belonging to the same domain. The assessment of the snapshot u(p) of the system response, to the same considered external action, can be carried out by FEM once again (or by any method of computational mechanics apt to solve ‘‘direct’’ problems). An alternative approach, as proposed e.g. in [13], is centered on suitable interpolation of the available simulation results gathered in the snapshot matrix U and approximated by truncated POD according to the procedure outlined in Section 2.2, Eqs. (1)–(5). The approximation of the snapshot u corresponding to a new parameter vector p can be expressed like the approximation,

M −

bik · gi (p)

(7)

i =1

where bik are the coefficients of this combination to be evaluated. In the present context the functions gi (p) are chosen as Radial Basis Functions (RBF), i.e. they are functions of the ‘‘distance’’ (i.e. Euclidean norm of coordinate differences) between the argument p of the function and the ith ‘‘node’’ pi : gi (p) = gi (‖p − pi ‖).

(8)

In this case, the interpolation coefficients bik can be determined in a conventional way when RBF are used for data interpolation, namely by satisfying exactly Eq. (7) in all ‘‘nodes’’ [16]. Clearly, in each node pi (i = 1, . . . , M) the amplitudes a⌢k are known at this ⌢ stage and collected in the K × M matrix A . Let the values of chosen RBF in all nodes be gathered in the M × M matrix G, namely:

(5)

Recourse to a truncated POD basis instead of the original full one, naturally introduces some error. However, it turns out that it is quite easy to evaluate the amount of this error and make it negligible for practical purposes. In fact, an assessment of the error in the truncated formulation can be done by comparing the sum of the eigenvalues corresponding to the K kept POD ‘‘modes’’ (i.e. directions) with the summation of all N of them. In fluid mechanical problems where the snapshot matrix is related to velocity field, it can be shown [19] that the ith eigenvalue quantifies the kinetic energy related to the ith POD mode; then the sum of neglected eigenvalues gives information about the loss of energy introduced by the low-order approximation. In other physical contexts there is no direct analogy; however, still the ratio between the sum of the preserved eigenvalues and the sum of all of them gives global quantitative assessment about the accuracy of the approximation by the above specified ‘‘truncation’’. The approximation expressed by Eq. (5) concerns the nodes of the pre-selected grid in the space of parameters p. For the present purposes it is necessary to derive an approximation apt to quickly compute the measurable quantities in the structural response to the test for any arbitrary combination of parameters within the pre-selected domain. To achieve this goal, in what follows recourse is made to interpolation by Radial Basis Functions.

(6)

The new amplitude vector aˆ is a (generally non-linear) function of parameter vector p. Each component a⌢k (k = 1, . . . , K ) of aˆ can be expressed as a linear combination of M suitably chosen functions gi , namely:

 G=

g1 (|p1 − p1 |)

...

gM (|p1 − pM |)

... ... ...

g1 (|pM − p1 |)

...



gM (|pM − pM |)

.

(9)

Using the above matrices of data, in order to determine the coefficients bik gathered in the K × M matrix B, the linear equation to solve reads:

ˆ = B · G. A

(10)

After this computation, combining Eqs. (6) and (7) leads to the following equation which represents POD–RBF approximation for any arbitrary set p of parameters:

ˆ · B · g(p). u(p) ≈ 8

(11)

As for structural diagnoses by indentation, once a series of M simulations are performed and the resulting data are used in the above presented way, Eq. (11) provides approximate measurable quantities in the response to the indentation test, without the need to perform any further FE simulation. It is reasonable to expect that the above interpolation by linear matrix computations is by far more computationally economical than non-linear FEM as a tool for test simulations. Approximation defined by Eq. (11) involves a set of FE simulations a priori ˆ and B. performed in order to determine the entries in matrices 8 This process called ‘‘training’’ of the procedure is done once-for-all. The operations that it involves are summarized by the flow-chart in Fig. 1 and the main symbols of recurrent vectors and matrices are gathered and clarified in Table 1. The above procedure of condensing a relatively large number of numerical simulations within a preliminary phase of training is somehow similar to procedures centered on Artificial Neural Networks (ANN). Some advantages of the present method with respect to ANN consist of easier ways to control computational errors, and a better capability of dealing with noisy experimental data. This is consequence of regularization provided by the use of POD in the first phase of the procedure. Such advantages also occur with respect to the Lagrange interpolation, proposed by Aoki et al. [20]. It can be shown that the effectiveness of the latter technique drops with input of noisy data. A further advantage of the POD–RBF procedure presented here arises from the fact that it doesn’t have to use a regular (i.e. equally spaced) grid of points in parameter space for the training phase;

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

495

Table 1 Meaning of the main symbols in this paper. Vectors computed for any set of parameters (vector p)

u g P

Matrices computed once-for-all in the phase of ‘‘training’’

U

ˆ A ˆ 8 G B

N-vector of computed quantities representing the response of a system to an indentation test (e.g. penetration depths for given load levels) M-vector of the values acquired by the RBF as functions of the parameter vector p S × M matrix of parameters; each column collects one parameter combination p corresponding to one ‘‘node’’ of the grid preselected over the chosen domain in the parameter space N × M matrix called the ‘‘snapshot matrix’’; each column represents one response u to the test with the parameters placed in the corresponding column in matrix P K × M matrix of ‘‘amplitudes’’ of ‘‘snapshots’’ in the ‘‘truncated’’ basis N × K matrix of ‘‘truncated’’ POD basis M × M matrix collecting values of the adopted RBF in all ‘‘nodes’’ (i.e. columns of matrix P) K × M matrix collecting the coefficients for RBF interpolation apt to provide amplitude vector aˆ .

Fig. 2. Finite element model used herein for indentation test simulations.

Unlike classical, local interpolation techniques, the present procedure approximates the response of the system by means of just one function (Eq. (11)) over the whole domain of interest. This approach is therefore effective for the repetitive test simulations required by iterative optimization algorithms, as will be shown by numerical exercises in the next sections. 3. Indentation testing and its fast simulation 3.1. Model formulation

Fig. 1. Flow chart of the operations to be performed in a POD–RBF procedure (Snumber of parameters to identify; N-quantities to be measured in the test; Msnapshots, namely test simulations based on pre-selected parameters).

an important feature when the number of parameters to identify is large.

In what follows the proposed procedure for material characterization will be illustrated and evaluated by means of computational exercises with reference to a practical situation specified here below for its pertinent features. A Rockwell conical indenter is considered with: 120° opening angle, 200 µm tip radius, Young modulus 1170 GPa, Poisson ratio 0.07. The maximum load is equal to 150 N. The isotropy of the material to be tested suggests axisymmetrical modeling. The FE model adopted for the test simulation by the commercial code ABAQUS [10] is visualized in Fig. 2. It exhibits the following features: rigid horizontal alignment of the top boundary of the modeled part of the indenter subjected to the imposed load; zero displacements on the boundaries common with the surrounding practically undisturbed zone (adopted here for simplicity instead of a more accurate alternative such as statical condensation of the surrounding solids); 1715 four-nodes quadrilateral elements (hence 3556 degrees of freedom); large elastic–plastic strains and Coulomb friction unilateral contacts (with

496

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

Fig. 3. Quantities entering into the snapshots: (a) from indentation curves; (b) from residual imprint.

coefficient 0.15) according to their implementation in ABAQUS [10] and not discussed here. To the material of the laboratory specimen (or of the structure to be examined in situ) remarkable ductility is attributed as in many structural metals and, hence, the classical Huber–Mises elastic–plastic model is adopted. Whereas the Poisson ratio is a priori assumed, ν = 0.3, the Young modulus E and the yield stress σY are the parameters to identify within the following reasonable lower and upper bounds: 140 GPa ≤ E ≤ 220 GPa;

300 MPa ≤ σY ≤ 500 MPa.

(12)

In view of the snapshot generation and interpolation envisaged in Sections 1 and 2, within the above domain in the (here twodimensional) parameter space, let a grid of M = 189 ‘‘nodes’’ pi be defined by subdividing the intervals (12) into steps of 10 GPa for E and of 10 MPa for σY . The coordinates of all the above nodes are gathered in the parameter matrix P ≡ [p1 , p2 , . . . , pM ] of S × M size, with S = 2. The structural response (‘‘snapshot’’) to the above specified indentation test here varies only with variation of the two material parameters. It consists of many quantities, which, to the present purpose, are classified and collected into the following submatrices of the snapshot matrix U with M columns: (i) NL × M submatrix UL for measurable quantities from the loading–unloading curves, namely NL = 2 × 50 penetrations under given loads at equal intervals (see Fig. 3(a)) and NI × M submatrix UI for measurable quantities from the imprint profile, namely NI = 82 heights from the horizontal specimen surface (Fig. 3(b)); (ii) computable (but not all measurable) response quantities such as: nodal displacements in the (3556 × 189) submatrix UD ; von Mises equivalent stresses σM and principal stress components σ11 and σ22 in all Gauss points, gathered in submatrices UM , U11 and U22 , respectively, each one with dimensionality 6860 × 189 according to the FE mesh shown in Fig. 2. 3.2. Comparative computational exercises of direct analysis The response quantities contained in the snapshot submatrices UD , UM , U11 and U22 , above defined at (ii), result herein from 189 traditional FE simulations of the experiment described in the preceding subsection, each simulation being based on a parameter pair pi = {Ei , σYi }T , i = 1, . . . , 189. Such a routine computational effort, required about 8 h of computing time (an average of 2.5 min for a single simulation) on a Pentium Dual Core computer, 2.5 GHz with 2 GB of RAM, by using the commercial code ABAQUS [10]. The POD procedure outlined in Section 2.2 has been applied to the snapshot sets generated as mentioned above. Using the generated snapshot matrices, UD , UM , U11 and U22 , the eigenvalues and eigenvectors of positive semi-definite matrix D, Eq. (1), are computed and the orthonormal basis vectors are determined, Eq. (3). For the truncation the criterion was adopted of reducing to

less than 10−6 the ratio between the summation of neglected eigenvalues and the summation of all of them. This criterion led to a different number of ‘‘modes’’ preserved for each one of the four above specified snapshot matrices. The most correlated data turn out to be those collecting nodal displacements, since for them only the first four new reference axes were sufficient to satisfy the above specified accuracy criterion. The basis for the snapshot matrix of the indentation curve was truncated after the seventh new reference direction, while those for von Mises stress σM and for principal stresses σ11 and σ22 after the 14th, 10th and 12th new basis directions, respectively. It is worth noting that the above truncation implies a large reduction of the dimensionality (e.g. the stress snapshot size was originally N = 6860). Now a new parameter vector p is considered, specifically E = 195 GPa and σY = 405 MPa, diverse from the 189 vectors employed as grid nodes for the training phase. Once again a traditional FE simulation of the test provides the snapshot vector u, subdivided in subvectors uL , uD , uM , u11 and u22 according to the previously adopted subdivision of the response quantities to be considered for evaluative comparisons. Alternatively, the same quantities have been computed by interpolation exploiting the truncated POD results through RBF using Eq. (6). In this study, ‘‘inverse multiquadric’’ RBF are employed, expressed by the following formula: gi (p) = 

1

‖p − p i ‖2 + r 2

(13)

where r is a ‘‘smoothing coefficient’’, which has been assumed herein equal to 0.65, while parameters are normalized between 0 and 1 over the adopted intervals for training. When the determination of interpolation coefficients involves a large number of linear equations to be solved, here Eq. (10), high values (anyway within the 0–1 range) of smoothing coefficients are computationally advantageous but can cause numerical instabilities with matrix G being almost singular, as shown e.g. in [17]. Since values of r closer to 1 generally lead to better results, the usual practice is to adopt an r value that is as large as possible [17]. In the present case, Fig. 4 visualizes results for the indentation curve obtained for r = 0.1 and for r = 0.65 and compares them to their counterpart achieved by conventional FE analysis rather than by POD–RBF procedures. For practical applications of parameter identification by discrepancy function minimization, crucial are comparisons between the following alternative approaches to the test simulation sequence (i.e. generations of measurable quantities) within the minimization procedure: (a) conventional FE analyses; (b) the above described POD–RBF procedure performed after conventional FE analyses carried out ‘‘a priori’’ with parameters suitably selected over a domain in their space. Such comparisons (which, clearly, require same FE model in both approaches) are meaningful in terms of: (i) additional errors implied by POD truncation and RBF interpolation with respect

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

497

of an indentation test. The inverse problem is solved using TRA implemented in MATLAB and the repeated test simulations have been performed by the POD–RBF procedure described in the preceding section. In this procedure the preliminary training was carried out by direct analyses using commercial code ABAQUS and the provisions implemented in it for large strains [10]. All features of the FE model are here the same as in Section 3. The material model adopted here is isotropic elasto–plastic with exponential hardening. Its classical tensorial description (omitted here for brevity, see e.g. [10]) is based on the traditional Ramberg–Osgood model, which in its original holonomic uniaxial formulation reads:

Fig. 4. Effects of different smoothing coefficient r assumed in the multiquadric Radial Basis Functions (RBF), Eq. (13).

to the results provided by another conventional FE analysis; (ii) computing times in (a) compared to those in (b), but disregarding the computational efforts required by preliminary ‘‘a priori’’ training. The following circumstances are worth noting: the comparative criterion (i) is autonomous from computational comparisons among minimization algorithms (here first-order TRA, but zeroorder procedures like Nelder–Mead or genetic algorithms might be employed elsewhere); comparative criterion (ii) becomes crucial when repeated practical diagnostic assessments are needed routinely far from computing centers. The traditional and popular TRA employed herein is explained with details and implemented in MATLAB (see e.g. [11,12]); therefore only the following main features of it are reminded here: (i) at the beginning of each step the objective function ω is replaced by its quadratic approximation; (ii) its Hessian is approximated by means of its gradient and of the Jacobian of the measurable quantities, so that first derivatives only are computed (by forward finite differences); (iii) each iteration is reduced to quadratic programming in two variables along gradient and Newton vector, under modifiable ‘‘trust region’’ constraints. The above considerations motivate the present computational exercises, some results of which are visualized in Fig. 5. From such exercises it turns out that the results obtained by lowdimensional approximation of the indentation test and those results achieved by FE conventional analysis are practically the same. In fact, the differences in resulting nodal displacements mapped in Fig. 5(a) can be compared to, say, the maximum pile-up on the deformed surface; for the differences in resulting stresses visualized in Fig. 5(b)–(d) the yield stress σY = 405 MPa provides a reference value. In both kinds of response to an indentation test the computed differences between POD–RBF results and FEM turn out to be negligible in practical purposes. In terms of computing times, after ‘‘training’’, the POD–RBF generation of the response to an indentation test turns out to be about 20000 times faster than FE simulation. 4. Identification of elastic–plastic material parameters by TRA combined with POD and RBF 4.1. Validation by pseudo-experimental data In order to validate computationally the practical usefulness of the POD–RBF procedure for diagnostic structural analysis based on indentation, a numerical exercise is discussed here below on the identification of constitutive parameters for metallic materials starting from pseudo-experimental data provided by a simulation

σ (ε) =

  E · ε,

if ε ≤

 σ01−n · E n · εn ,

if ε >

σ0 E

σ0 E

.

(14)

In this case, there were three parameters to identify, namely: Young modulus E, yield stress σY and hardening exponent n. For the inverse analysis both indentation curves and residual imprint profiles are used and, therefore, two different snapshot matrices are constructed. As for the former snapshots concerning the indentation curves, the same procedure as in Section 3 leads to the size N = 100. The latter snapshot set consists of vectors collecting profile heights (with respect to the specimen surface) consequent of nodal displacements of the material points on the surface of the specimen. The connections between imprint profile heights and material particle displacements in a large-strain context are implemented in the commercial code ABAQUS [10] and are not discussed here. Since there is a total of 82 nodes in this region and the model is two-dimensional, in view of axial symmetry, the size of each snapshot is N = 82. In the first step, the snapshot matrices are constructed by a set of FE direct analyses, varying σ0 from 75 to 210 MPa with the step of 15 MPa, E from 50 to 120 GPa with the step of 10 GPa and the hardening exponent n from 0.2 to 0.5 with the step of 0.01. The total number of nodes in the parameter space turns out to be M = 2480. As a result of these simulations, the following snapshot matrices are constructed: U1 [82 × 2480] with surface displacements, and U2 [200 × 2480] with data concerning the indentation curve. Both POD bases are truncated after the sixth direction (namely after the sixth eigenvalue in order of magnitude). For the interpolation the same type of RBF as in the previous section, (Eq. (13)) is employed. A typical ‘‘pseudo experimental’’ exercise for the validation of the parameter identification method proposed herein is centered on the operative stages which follow. A parameter set pN , say E = 95 GPa, σ0 = 140 MPa and n = 0.295, is assumed different from those above considered as nodes of the grid for the generation of the snapshot matrices U1 and U2 . Once again, an indentation test is simulated by FEM for the selected parameters pN . The measurable quantities thus computed are gathered in snapshot vectors u1 and u2 (concerning imprint profile and indentation curve, respectively) and employed as a pseudo-experimental input in the inverse analysis procedure (POD–RBF–TRA) described in Section 3. The resulting estimates are compared to those in vector p earlier assumed and now acting as ‘‘targets’’: clearly, the differences quantify the accuracy of the estimation. The ‘‘discrepancy function’’ ω to be minimized is formulated here as the simplest one in inverse problems, i.e. the Euclidean norm of the differences between (here pseudo-) experimental data and their counterparts computed on the basis of the parameters as current variables, namely:

ω(p) = R(p)T R(p),

R(z) = ue − u(p).

(15)

The present numerical exercise led to the following results. The convergence of the TRA iterative procedure on the expected values

498

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

a

b

c

d

e Fig. 5. Differences between results computed by FE simulation and by POD–RBF approximation: (a) nodal displacements; (b) von Mises equivalent stress; (c) normal stress σ11 in axis direction 1 (d) normal stress σ22 in axis direction 2; (e) indentation curves.

has been achieved starting from diverse initializations: Figs. 6 and 7 visualize two numerical tests, with normalizations of the variables parameters by means of the relevant targets. Fig. 8 provides a typical visualization example of the discrepancy function. One of the sought parameters, namely, the exponent of hardening, is kept fixed, specifically n = 0, while the other two change over the domain of interest. The ‘‘regularity’’ of discrepancy function, i.e. its convexity (and smoothness), indicates

well-posedness of the inverse problems and, at least qualitatively, sensitivity of the set of experimental data with respect to the remaining two searched parameters, E and σ0 . Finally, the computing times on the same computer as specified in Section 3, in the above typical validation exercise turned out to amount as follows (with equal initialization and equal stopping criteria): ∼3 s by means of the above POD–RBF–TRA; ∼4 h by means of the TRA implying FE simulations at each iteration step.

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

499

Fig. 6. Normalized values of parameters along the TRA iteration sequence with the first initialization.

Fig. 9. Measured residual imprint from the hardness test with maximum force of 1500 N; (a) 3D surface of the imprint; (b) average profile together with standard deviation. Table 2 Assessed parameters for copper subjected to indentation test. Fig. 7. Normalized values of parameters along the TRA iteration sequence with the second initialization. Initialization 1 Initialization 2 Initialization 3 Initialization 4

Fig. 8. Discrepancy function ω over the parameter domain with hardening set to zero (n = 0).

4.2. Validation with true experimental data In this subsection, the proposed POD–RBF procedure is applied on real experimental data collected from indentation tests on a copper specimen. As a first test, instrumented indentation was performed with maximum force equal to 200 N, and the indentation curve was used directly as input to the inverse analyses. A second test carried out was a classical hardness test (the indenter used was not instrumented) with applied force reaching 1500 N. The residual imprint after the removal of the indenter was measured by a contact profilometer with results visualized in Fig. 9. In view of the material isotropy and consequent expected axial symmetry of the imprint, a two-dimensional profile was computed as an average of 8 profiles (rotated by 45° from each other), to be used as inputs for the inverse analyses (Fig. 9(b)).

E (MPa)

σY (MPa)

n

ω

100,634 103,115 105,101 112,947

229 226 215 212

0.046 0.041 0.055 0.059

1.06 0.794 0.793 0.869

The system responses to both tests were computed by the proposed POD–RBF procedure, based on 832 analyses. These analyses were previously performed varying the three sought parameters in the following ranges (determined by ‘‘expert’s conjectures’’ on properties of the copper considered): yield limit between 150 and 255 MPa, step 15 MPa; Young modulus between 90 and 160 GPa, step 10 GPa; exponent of hardening between 0 and 0.12, step 0.01. The same RBF type as in the previous exercise was used, namely the one given by Eq. (13), with smoothing coefficient r = 0.5. The inverse problem was solved by means of TRA implemented in MATLAB [12] four times starting from randomly selected initialization vectors, in order to check whether the resulting estimates correspond to the global minimum of the discrepancy function. The assessed parameters are given in Table 2; they have been employed to construct the uniaxial stress–strain relationships visualized in Fig. 10. Such plots turn out to be satisfactorily close to each other, a circumstance which corroborates the validity of the parameter identification. Fig. 11 shows that computed system responses (namely indentation curve and residual imprint) are very close to their experimentally measured counterparts. These comparisons further check the accuracy of the inverse problem solution, but they also show that the trained POD–RBF procedure was capable of capturing the system responses resulting from the indentation tests. Finally it should be emphasized that using the described procedure based on POD–RBF–TRA, the identification process takes about 3 s. This means that, the procedure can be routinely used on a small computer, since, once it is trained, it does not require any further ‘‘heavy’’ and time consuming computations.

500

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501

Fig. 10. Stress–strain curves computed for the assessed parameters.

Fig. 11. Experimental results together with those computed by the POD–RBF–TRA procedure: (a) indentation curves; (b) residual imprints.

5. Closing remarks The structural diagnosis method proposed and preliminarily investigated in this paper for the mechanical characterization of ductile structural materials exploits two circumstances which usually occur in real-life diagnostic problems: search for constitutive parameters confined to a reasonably predictable domain in their space; correlation of the measurable responses to tests with diverse parameters within that domain. The non-destructive experiment considered herein is indentation (rooted in classical hardness test) combined with finite element simulation and inverse analysis based not only on force–penetration curves digitalized by an instrumented indenter, but also on residual imprint geometry (usually average profile) measured by a laser or contact instrument. Such an increase of experimental data was shown to be practically advantageous in preceding research [3,6,7]. The procedure proposed herein for deterministic practical diagnosis of structural material deteriorations implies the following sequence of operative stages:

(a) As a preliminary phase, over the above mentioned domain numerous parameter vectors (or ‘‘nodes’’ in a suitably selected grid) are used as inputs into test simulations by some conventional finite element code and the resulting vectors of measurable response quantities (‘‘snapshots’’) are computed once-for-all and stored. (b) The accumulated large set of snapshots is approximated by a ‘‘truncated’’ Proper Orthogonal Decomposition (POD) basis, namely by a classical tool of numerical mathematics apt to exploit correlation in a large set of vectors, in order to drastically simplify (i.e. reduce dimensionality), with optimal approximation, the informative contents of that set. The main computational effort consists of the solution (again once for all) of a single problem concerning eigenvalue analysis of a large, but symmetric and positive semi-definite, matrix. (c) For any new vector of sought parameters in the earlier specified domain, the relevant response snapshots of measurable quantities are computed (rather than by further finite element simulation of the test) by interpolation of the results achieved in the two preceding phases and stored. Such interpolation is done by means of suitably chosen Radial Basis Functions (RBF). Identification of interpolation coefficients requires the solution of a linear algebraic problem, after which the system responses are obtained straightforward by matrix multiplication. (d) The customary minimization, with respect to the sought parameters, of a suitable norm of the differences between experimental data and their computed counterparts has been performed here by a recent version of the first-order Trust Region Algorithm (TRA). The important provision (c) adopted for the many direct analyses required by TRA at each iteration (consisting of two-variable quadratic programming) leads to a drastic reduction of computing time for each parameter. Therefore the investigated procedure, after the above specified training in a computing center, can be carried out routinely, economically, quickly and even ‘‘in situ’’ on a portable computer. Clearly, the above conclusions, reached in this paper with reference to inverse analyses by TRA, hold for other iterative firstorder optimization algorithms and for the zero-order techniques like genetic algorithms. These soft-computing techniques are frequently suitable for the present estimation problems in view of a possible lack of convexity and local minima exhibited by the discrepancy function to be minimized. Research in progress concerns extensions of the present results to other diagnostic techniques at present frequently employed in engineering (such as flat jack tests in civil engineering) and to stochastic approaches (Monte Carlo, Bayesian techniques, Kalman filter [20,21]) to parameter identifications where the assessment of estimate uncertainties (i.e. their covariance matrix) due to experimental ‘‘noises’’ generally requires a significant increase of the number of correlated simulations. Acknowledgements The authors acknowledge with thanks the support from ENI and VETEC to a research project related to the results here summarized and the inspiration for POD–RBF adoption from a colleague, Ostrowski, during his visit. References [1] Nakamura T, Wang T, Sampath S. Determination of properties of graded materials by inverse analysis and instrumented indentation. Acta Mater 2000; 48:4293–306. [2] Van Vliet KJ, Gouldstone A. Mechanical properties of thin films quantified via instrumented indentation. Surf Eng 2001;17(2):140–5. [3] Bolzon G, Maier G, Panico M. Material model calibration by indentation, imprint mapping and inverse analysis. Internat J Solids Structures 2004;41: 2957–75.

V. Buljak, G. Maier / Engineering Structures 33 (2011) 492–501 [4] Oliver WC, Pharr GM. An improved techniques for determining hardness elastic modulus using load and displacement sensing indentation experiments. J Mater Res 1992;7:176–81. [5] Van Landingham MR. Review of instrumented indentaion. J Res Natl Inst Stand Technol 2003;108:249–65. [6] Bocciarelli M, Bolzon G, Maier G. Mechanical characterization of anisotropic materials by indentation test, imprint mapping and inverse analysis. Mech Mater 2005;37:855–68. [7] Bocciarelli M, Maier G. Indentation and imprint mapping method for identification of residual stresses. Comput Mater Sci 2007;39:381–92. [8] Maier G, Bocciarelli M, Bolzon G, Fedele R. Inverse analyses in fracture mechanics. Int J Fract 2006;138:47–73. [9] Tikhonov AN, Arsenin VR. Solution of ill-posed problems. Washington: Winston and Sons; 1977. [10] ABAQUS/standard, theory and user’s manuals, release 6.2-1. Pawtucket (RI, USA): HKS Inc; 1998. [11] Nocedal J, Wright SJ. Numerical optimization. New York: Springer-Verlag; 2006. [12] Matlab. User’s guide and optimization toolbox, release 6.13. USA: The Math Works Inc; 2002.

501

[13] Ostrowski Z, Bialecki RA, Kassab AJ. Solving inverse heat conduction problems using trained POD–RBF network. Inverse Probl Sci Eng 2008;16(1): 705–14. [14] Ostrowski Z, Bialecki RA, Kassab AJ. Estimation of constant thermal conductivity by use of Proper Orthogonal Decomposition. Comput Mech 2005; 37:52–9. [15] Liang YC, Lee HP, Lim SP, Lee KH, Lin WZ. Proper Orthogonal Decomposition and its applications — part i: theory. J Sound Vibration 2002;252(3): 527–44. [16] Kansa EJ. Motivations for using radial basis functions to solve PDEs; 2001. p. 1–8. http://rbf-pde.uah.edu/kansaweb.pdf. [17] Buhmann DM. Radial basis functions. Cambridge University Press; 2004. [18] Jolliffe IT. Principal component analysis. New York: Springer-Verlag; 1986. [19] Holmes P, Lumley JL, Berkoz G. The proper orthogonal decomposition in the analysis of turbulent flows. Annu Rev Fluid Mech 1993;25:539–75. [20] Aoki S, Amaya K, Sahashi M, Nakamura T. Identification of Gurson’s material constants by using Kalman filter. Comput Mech 1997;19:501–6. [21] Bittanti S, Maier G, Nappi A. Inverse problems in structural elasoplasticity: a Kalman filter approach. In: Sawaczuk A, Bianchi G, editors. Plasticity Today. Applied Science Publications; 1984. p. 311–29.