I. INTRODUCTION
Properly designed multilayer protective thin films have been reported as having enhanced mechanical and tribological properties relative to the monolayers of equivalent thickness, with improvements attributed to a variety of mechanisms such as the superlattice effectReference Molina-Aldareguia, Lloyd, Oden, Joelsson, Hultman and Clegg1 and the modified Hall–Petch effect.Reference Kang, Kim, Kwon and Hong2 Considering materials deposited in multilayers may have considerably different microstructures compared with the same materials deposited as monolayers, characterizing the intrinsic mechanical properties of the individual constituent layers is of great importance to explore the dominant mechanism for the enhancement of the properties to design multilayer films with desired properties. As the most effective tool to investigate the mechanical properties of thin film structures, nanoindentation technique has gained significant advances in measuring the “film-only” properties of monolayer film structures either by experimentally suppressing the substrate effect (e.g., to limit the indentation depth to <10% of the film thickness) or theoretically modeling the substrate effect during interpreting the indentation data.Reference Yu, Sanday and Rath3–Reference Hay and Crawford11 However, due to the difficulty in separating the effects of the constituent layers, the indentation methods appropriate for monolayer films are generally not applicable to multilayered film structures. As a consequence, the film-only properties of the individual layers of the multilayered films are rarely characterized.
The indentation response of a multilayer thin film on a substrate is a function of the mechanical properties of the constituent layers and the indentation depth if the substrate and the multilayer configuration (i.e., the thickness and the laying order of the layers) are known.Reference Chima-Okereke, Bushby, Reece, Whatmore and Zhang12 In principle, the mechanical properties of the individual layer can be inferred from the indentation tests performed at different depths through an inverse analysis procedure. Unlike monolithic materials or monolayer films of which the inverse analysis procedures have been established and well validated, inverse analysis procedures appropriate for the stacking layers of the multilayered film structures are far less investigated due to the complicated interactions of the substrate and the multiple layers.
In fact, the inverse analysis of an indentation test can be viewed a fitting process during which the difference between the predicted indentation data and the measured indentation data is minimized. Attributed to the rapidly increasing computing power, more general inverse analysis methods combining finite element simulations and numerical optimizations have been developed in recent several years.Reference Peyrot, Bouchard, Ghisleni and Michler13–Reference Hamasaki, Shinbat and Yoshida15 Various optimization algorithms have been implemented with the deviation between the experimental indentation data and the corresponding simulation result as the objective functions. Peyrot et al.Reference Peyrot, Bouchard, Ghisleni and Michler13 proposed a coupling of semianalytical solution with an Evolution Strategy Algorithm to ensure both efficiency and global optimization while Rauchs and BardoReference Rauchs and Bardo14 adopted a gradient-based method, and the gradient of the objective function with respect to the material parameters is calculated using the direct differentiation method. To reduce computation time, an algorithm combing high- and low-fidelity models was developed to identify multiple parameters.Reference Hamasaki, Shinbat and Yoshida15 Considering various constitutive laws, loading schemes, and structural configurations can be conveniently defined in the finite element model, these methods are especially suitable with identifying complex materials and structures. Although multiple finite element calculations are needed to simulate the indentation test during the optimization process, the computation cost is manageable if an appropriate optimization algorithm is chosen and this is especially true when a parallel version of the optimization algorithm is implemented.
Very recently, Karimpour et al.Reference Karimpour, Balint, Rzepiejewska-Malyska, Szerling, Michler and Lin16 have proposed an inverse method to determine the elasto-plastic parameters of the constituent materials of a multilayer by least squares fitting of the finite element simulation indentation curves to the experimental data. The study demonstrated that the rule of mixtures, although a sensible starting point, is inappropriate for obtaining the mechanical properties of the bimaterial multilayers.Reference Karimpour, Balint, Rzepiejewska-Malyska, Szerling, Michler and Lin16 However, to extract the mechanical properties of the constituent materials of a multilayer film, the method requires multiple indentation curves of the possible bimaterial combinations, and for each material, an elastic perfectly plastic material model has to be assumed.Reference Karimpour, Balint, Rzepiejewska-Malyska, Szerling, Michler and Lin16 Considering this may not be well satisfied in practical applications, the present study proposes a simpler inverse method to determine the elastic modulus of the each layer of multilayered films on a substrate. The proposed method requires only multiple loading–unloading indentation cycles performed under different load levels. This inverse analysis procedure is first numerically verified against monolayer, bilayer, and trilayer films and then applied to identify the elastic modulus of monolayered tungsten film, aluminum film, and bilayered tungsten/aluminum from the available nanoindentation test. Uniqueness and depth independence of the extracted parameters suggest that this is a feasible approach to characterize the intrinsic elastic properties of the constituent layers of multilayered film structures by using nanoindentation techniques.
II. GENERIC INVERSE ANALYSIS PROCEDURE
As mentioned above, the proposed inverse analysis is treated as an optimization problem of which the objective function is the difference between the simulated indentation data and the measured indentation data. Considering the focus of the present study is to determine the elastic properties of the individual layers, the indentation contact stiffness (i.e., the initial unloading slope) or the nominal-reduced modulus collected at different indentation depths, are chosen as the quantities to be fitted since it is believed that the quantities are independent of the plastic properties of the constituent layers.Reference Karimpour, Balint, Rzepiejewska-Malyska, Szerling, Michler and Lin16 Considering the fitting nature of the present method, the quantities can be determined either by Oliver and Pharr's methodReference Oliver and Pharr17 or by Doerner and Nix's methodReference Doerner and Nix18 as long as a consistent method is used to analyze the experimental results and the corresponding numerical results.
Define the objective function to be minimized as
$$F\left( {{{\bf E}^k}} \right) = {{{{\mathop\sum\limits_{i = 1}^N {\left( {{{{f^{\exp }}\left( {{h_i}} \right) - {f^{{\rm{cal}}}}\left( {{{\bf E}^k},{h_i}} \right)} \over {{f^{\exp }}\left( {{h_i}} \right)}}} \right)} }^{\!2}}} \mathord{\left/ {\vphantom {{{{\sum\limits_{i = 1}^N {\left( {{{{f^{\exp }}\left( {{h_i}} \right) - {f^{{\rm{cal}}}}\left( {{E^k},{h_i}} \right)} \over {{f^{\exp }}\left( {{h_i}} \right)}}} \right)} }^2}} {F({E^0})}}} \right. \kern-\nulldelimiterspace} {F({{\bf E}^0})}}\quad ,$$where Ek is the list of the elastic moduli to be identified of the constituent layers at the kth iteration, f exp(h i) is the measured indentation response (the indentation contact stiffness or the effective-reduced modulus) at the maximum indentation depth of h i, and f cal(Ek, h i) is the calculated indentation response at the same indentation depth which is obtained from finite element simulations. In the present study, the Levenberg–Marquardt (L–M) algorithmReference Levenberg19,Reference Marquardt20 is chosen to minimize the objective function to mitigate the disadvantage of the Gauss–Newton method. In a typical L–M procedure, the parameters to be identified adjusted in an iterative way (Fig. 1):
$${{\bf E}^{k + {\bf{1}}}} = {{\bf E}^{k + {\bf{1}}}} + {{\bf g}^k}\quad ,$$where the correction term, i.e., gk, is determined by solving the equation:
$$\left( {{{\bf{A}}^T}{\bf{A}} + {\rm{\lambda }}{\bf{I}}} \right){{\bf{g}}^k} = - {{\bf{A}}^T}F\left( {{{\bf{E}}^k}} \right)\quad ,$$where 
${\bf{A}} \,\equiv\, {{{\partial {\bf{F}}} \mathord{\left/ {\vphantom {{\partial {\bf{F}}} {\partial {\bf{E}}}}} \right. \kern-\nulldelimiterspace} {\partial {\bf{E}}}}^k}$ is the Jacobian matrix which is evaluated by the finite difference method, I is the matrix identity and λ is a scalar. It is noticed that if λ = 0, the Gauss–Newton method is recovered and if λ → +∞ = 0, the steepest descent method is recovered. To gradually move from the steepest descent method to the Gauss–Newton method, λ can be evolved during the iteration. For example, λ is increased by a factor of 10 if the objective function decreases in an iteration step, and otherwise, λ is decreased by a factor of 10.
Flowchart of the inverse analysis procedure within the framework of the L–M method.
The axisymmetric finite element mesh used to calculate the simulated indentation response is shown in Fig. 2. The Berkovich indenter is replaced with an equivalent conical indenter with a half angle of 70.3°. Only three layers of films are included in the mesh for illustration and all the layers share the same thickness. Monolayer and bilayer films can also be modeled by the mesh by simply replacing the mechanical properties of the second and/or the third layer with those of the substrate.
The finite element mesh used to calculate the indentation response of a trilayer film system.
To minimize the effect of boundary conditions, the width and the height of the indented solid are both set to be 10 times of the potential contact radius. Presently, the contact between the indenter and the solid is assumed to be frictionless. 4960 quadrilateral axisymmetric elements are used to discretize the indented solid with elements refined in the potential contact region. Each layer contains six elements through its thickness. The meshing strategy has been tested to be a best compromise between computing time and accuracy of the numerical results. Convergence studies are conducted, and insensitivity of the results to the mesh is ensured.
The material size effect generally associated with the superlattice effect is not considered in the present study. The indented layers are assumed to obey the J 2 flow theory. A simple linear hardening is also incorporated for the constituent layers to check the uniqueness of the extracted moduli since it is found that inclusion of both yield stress and hardening could lead to nonuniqueness problem.Reference Luo and Lin21 Without loss of generality, the tangent modulus is set to be tenth of the elastic modulus throughout the study.
Both the layers and the substrate are considered to be initially stress-free, although residual stresses resulting from the deposition process may affect the measured indentation response to some degree. A displacement-control boundary condition with multiple loading and unloading segments is applied. The indentations penetrate the films layer by layer and it is suggested that for each layer, at least one indentation be performed at the depth deeper than 20% of the thickness of the layer. An open-source finite element solver, Code-Aster,22 is chosen to conduct all the finite element calculations.
III. VERIFICATION AND APPLICATION
A. Numerical verification
The proposed inverse analysis method was firstly verified by numerical indentation tests of several typical monolayer, bilayer, and trilayer films. The contact stiffness collected from the simulated indentation cycles was fed into the inverse analysis procedure as the indentation response to be fitted. Both the simulated indentation tests and the inverse analysis procedure share the same finite element mesh. Two, three, and four loading–unloading cycles were, respectively, chosen for the monolayer, bilayer, and trilayer films. The inverse analyses were repeated five times, and for each time, the initial values of the elastic moduli to be identified were chosen randomly. Both the mean value and the standard error of the identified parameter were reported.
The exact elastic moduli and the identified counterpart of typical monolayer, bilayer, and trilayer films are shown, respectively, in Tables I–III. 
$E_{\rm{f}}^i$, 
${\rm{\nu }}_{\rm{f}}^i$, 
${\rm{\sigma }}_{\rm{y}}^i$, and 
$H_{\rm{f}}^i$, respectively denote the Young's modulus, the Poisson's ratio, the yield strength, and the tangent modulus of the i-th layer, respectively. E s and νs denote the Young's modulus and the Poisson's ratio of the purely elastic substrate, respectively. Only for the purpose of verification, all the substrates were assumed to be purely elastic in the section. It is evident that all the elastic modulus of the films can be accurately and uniquely extracted from the numerical indentation measurements. This suggests that the proposed indentation scheme and the inverse analysis procedure are suitable for the multilayer films.
Comparison between the indentified elastic modulus and the exact elastic modulus of typical monolayer films.
Comparison between the indentified elastic moduli and the exact elastic moduli of typical bilayer films.
Comparison between the indentified elastic moduli and the exact elastic moduli of typical trilayer films.
In practical applications, the plastic properties of the constituent layers are usually not known a priori. To investigate the effect of the plastic properties on the extracted elastic parameters, typical yield strengths at the levels of 25%, 50%, 100%, 150%, and 200% of the “true” yield strengths were, respectively, set for the constituent layers to generate groups of different simulated indentation data which were then fed into the inverse analysis procedure. It is found that if the assumed yield strength is within the range of 50–150% of the true value, the deviation of the extracted elastic modulus from the known exact modulus is negligible. Considering the plastic properties of the bulk materials of most thin films usually fall into the scope, the unknown plastic properties of the constituent layers can be set to be those of their bulk counterparts and this would not affect the extraction of the elastic parameters very much.
B. Application
The inverse analysis method was applied to reanalyze the well-documented experimental data of several typical monolayer filmsReference Saha and Nix6 and a bilayer film.Reference Shang, Yang and Li23 Due to the measurement data of trilayer films or multilayer films is yet unavailable, experimental verification of the method against these even more complicated structures is left for our further work.
For the monolayer films, both soft films on hard substrates, i.e., 500-nm aluminum films on the sapphire substrate and on the silicon substrate, and hard films on soft substrates, i.e., 640-nm tungsten films on the silicon substrate and on aluminum substrates, are included. The elastic moduli of the substrates were determined by using Oliver and Pharr's method: for sapphire E s = 440 GPa, for silicon E s = 172 GPa,, and for aluminum E s = 75 GPa.Reference Saha and Nix6 The elastic moduli of the aluminum films and the tungsten films were identified at series of indentation depth no larger than the thickness of the films (Figs. 3 and 4). The plastic properties of the aluminum films and the tungsten films were set as those of the bulk counterparts. For the systems of hard films on soft substrates (i.e., the 640-nm tungsten films on the aluminum substrates), the substrates were assumed to be elasto-plastic to capture the possible plastic deformation of the substrates. Negligence of the plastic deformation of the substrates would lead to an underestimation of the elastic properties of the films.
The identified elastic moduli of the aluminum films on the sapphire substrate and on the silicon substrate.
The identified elastic moduli of the tungsten films on the aluminum substrate and on the silicon substrate.
For the tungsten films, the silicon and the aluminum substrates were assumed to be elasto-plastic, and the plastic properties were as those of the bulk counterparts. The nominal-reduced modulus determined by using Oliver and Pharr'sReference Oliver and Pharr17 method was chosen as the fitted quantity. Two nominal-reduced moduli collected at two neighboring depths were fed into the inverse analysis procedure, and identified elastic modulus is reported on the mean indentation depth.
It can be seen that the indentified moduli are insensitive to the substrates and match the expected film modulus reasonably well (Figs. 3 and 4). For the aluminum film, the elastic moduli extracted by the present method seem to be insensitive to the indentation depth. This is different from the results by Saha and NixReference Saha and Nix6 who applied the modified King's solution to identify the elastic moduli, whereas for the tungsten films, both our results and Saha and NixReference Saha and Nix6 results show a large variation with the indentation depth (Fig. 4). By randomly choosing different initial values for the identification procedure, it is confirmed that the identified modulus at each indentation depth is uniquely determined. The variation should not be attributed to the possible nonuniqueness problem and needs further investigation.
The tungsten–aluminum bilayer film was deposited on a glass substrate of which the elastic modulus is around 70.0 GPa and the thicknesses of the two layers are both around 400 nm.Reference Shang, Yang and Li23 Three nominal-reduced moduli, i.e., 120.5 GPa, 87.9 GPa, and 88.7 GPa, which were measured, respectively, at the indentation depths of 100 nm, 295 nm, and 467 nm, were chosen as the fitted quantities. The inverse calculations were repeated five times, and for each time, random initial values were chosen. Initial values-independent results were achieved, and the mean identified elastic moduli of the tungsten film and the aluminum film are 350.5 GPa and 80.1 GPa, respectively. The identified moduli of both the films match their typical values reasonably well.
IV. DISCUSSION
The solution of Yu et al.Reference Yu, Sanday and Rath3 to the axisymmetric contact problem for a layered elastic half-space has been applied to analyze elasto-plastic indentations subject to various restrictions. For an elastic monolayer film, the indentation contact stiffness can be written as a dimensionless function by analyzing Yu's solution:
$${{\rm{S}} \over {E_{\rm{s}}^{\rm{r}}}} = f\left( {{{E_{\rm{f}}^{\rm{r}}} \over {E_{\rm{s}}^{\rm{r}}}},{a \over H}} \right)\quad ,$$where 
$E_{\rm{s}}^{\rm{r}} \,\equiv\, {{{E_{\rm{s}}}} \mathord{\left/ {\vphantom {{{E_{\rm{s}}}} {1 - \nu _{\rm{s}}^{\rm{2}}}}} \right. \kern-\nulldelimiterspace} {1 - \nu _{\rm{s}}^{\rm{2}}}}$ and 
$E_{\rm{f}}^{\rm{r}} \,\equiv \,{{{E_{\rm{f}}}} \mathord{\left/ {\vphantom {{{E_{\rm{f}}}} {1 - \nu _{\rm{f}}^{\rm{2}}}}} \right. \kern-\nulldelimiterspace} {1 - \nu _{\rm{f}}^{\rm{2}}}}$ are the reduced modulus of the substrate and the film, respectively, a is the contact radius, and H is the film thickness. For elastic films, the film thickness is well defined and a single indentation contact stiffness is enough to solve 
$E_{\rm{f}}^{\rm{r}}$, whereas for elasto-plastic films, the film between the indenter and the substrate is thinned as a result of plastic flow and it is not immediately clear what thickness should be used. Li et al.Reference Li, Randall and Vlassak10 found that the local thinning effect is well captured if an effective film thickness is used
$${H_{{\rm{eff}}}} = H - {\rm{\eta }} \times h\quad ,$$where H eff is the effective film thickness, h is the indentation depth, and η is the dimensionless parameter quantifying the local thinning of the film caused by plastic flow in the film. It was shown that this parameter is nearly independent on the indentation depth.Reference Li, Randall and Vlassak10Equations (4) and (5) imply that at least two measured quantities should be fed into the inverse analysis procedure to uniquely obtain the modulus of a monolayer film. This is why two nominal-reduced moduli were used when identifying the moduli of the aluminum films and the tungsten films in the above Sec. III. B.
A similar rule may be applicable to multilayer films as well and therefore, at least three nominal-reduced moduli should be used to identify the two true elastic moduli of the tungsten–aluminum bilayer films and at least four nominal-reduced moduli for the trilayer films. The rule ensures the uniqueness of the extracted parameters as demonstrated by the monolayer films and the bilayer film.
V. SUMMARY
An inverse analysis method to identify the elastic moduli of the constituent layers of multilayer films is developed by fitting the finite element calculations with indentation measurements within the framework of the Levenberg–Marquardt algorithm. The method is verified against typical monolayer, bilayer, and trilayer film structures both numerically and experimentally. Substrate-independent elastic moduli of the constituent layers are obtained, and uniqueness of the identified parameters is ensured by conducting an indentation containing multiple loading–unloading cycles. The method provides a feasible way to characterize the intrinsic mechanical properties of the constituent layers of multilayered thin films and further to explore the dominant mechanism for the enhancement of their mechanical properties.
ACKNOWLEDGMENTS
Financial supports from the National Natural Science Foundation of China (Grant Nos. 11002057, 10902128, 10732100, 50802026, 10972239, and 11072271) are gratefully acknowledged.
