- Email: [email protected]
piezoelectric constants of piezoelectric thin film by combining nanoindentation test with FEM
Computational Materials Science 63 (2012) 134–144
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Evaluation of engineering/piezoelectric constants of piezoelectric thin film by combining nanoindentation test with FEM S.T. Song a, X.J. Zheng a,b,⇑, H. Zheng a, W. Liu a a b
Faculty of Materials, Optoelectronics and Physics, Xiangtan University, Xiangtan, Hunan 411105, PR China School of Materials Science and Engineering, University of Shanghai for Science & Technology, Shanghai 200093, PR China
a r t i c l e
i n f o
Article history: Received 21 February 2012 Received in revised form 23 March 2012 Accepted 1 May 2012 Available online 30 June 2012 Keywords: Engineering constants Piezoelectric constants Finite element method Transversely isotropic Thin film
a b s t r a c t Based on the weight averaged method, the variety of the PZT thin film is determined by combining nanoindentation test with finite element method (FEM) simulation, and the piezoelectric constitution is used to establish the supplemental equation, in which the piezoelectric strain constants are related with the piezoelectric stress constants, so that all of the piezoelectric constants of piezoelectric thin film can be determined by combining nanoindentation test and FEM. In the forward analysis, the numerical loading curves at the broad spectrum of possible material combinations are simulated by using ABAQUS software to extract the numerical maximum indentation loads and loading curve exponents, and they are used to establish two dimensionless equations related with the engineering/piezoelectric constants of film/substrate system. In the reverse analysis, the loading part of experimental indentation curves performed on PZT thin film in nanoindentation test are fitted as the power function to obtain the maximum indentation loads and the loading curve exponents, and the engineering constants can be solved by using the simultaneity of dimensionless and supplemental equations. In order to verify the validity, the optimized engineering constants are inputted into ABAQUS software to obtain the numerical loading curves, and they are compared with the engineering constants of PZT thin film measured by the previous method to determine the variety of the PZT thin film. The piezoelectric constants of the certain variety are determined by combining forward analysis and inverse analysis. The engineering constants of the PZT thin film are determined as ET = 142 GPa, EL = 137 GPa, GL = 52 GPa, and the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined as e15 = 4.75 C/m2, e31 = 0.85 C/m2, e33 = 7.08 C/m2. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction The instrumented nanoindentation systems have been used to study the mechanical properties of thin film deposited on substrate effectively, and allow the application of a specified force or depth history [1–3], such that indentation load and the depth, are controlled and/or measured simultaneously and continuously over a complete loading cycle [4]. In general, the elastic modulus of isotropic thin film can be obtained from the indentation curve according to the method proposed by Oliver and Pharr [2]. The approaches cannot calibrate the whole set of the elastic parameters in anisotropic elasticity, even if indentations are performed in multiple directions and the shape of the contact area is exploited as supplementary information [5]. In nanoindentation test, the contact depth is generally less than one tenth of film thickness in order to obtain ⇑ Corresponding author at: Faculty of Materials, Optoelectronics and Physics, Xiangtan University, Xiangtan, Hunan 411105, PR China. Tel.: +86 731 58293648; fax: +86 731 58298119. E-mail address: [email protected] (X.J. Zheng). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.05.001
the intrinsic properties and avoid the so-called substrate effect [6,7]. However, the rule limiting the indentation depth may be only applicable to relatively thick films, and it is often not practical with the decrease of film thickness in many applications [8]. When the indentation depth is too shallow, measure error is induced by the small-scale effects, such as the adhesion effect between indenter tip and film, the strain gradient effect and the surface roughness effect [9]. Ma et al. [10] and Zhao et al. [9] considered the substrate effect during finite element method (FEM) simulation and nanoindentation test, and they proposed a method via combining nanoindentation with FEM to extract Young’s modulus and yield strength of isotropic thin films [11]. In our previous work, we have identified the elastic parameters of ZnO thin film when the substrate effect is considered for the transversely isotropic thin film [7]. Piezoelectric materials are a type of functional materials with mechanical and electrical coupling behaviors that are attracting increasingly wide engineering applications, such as actuators, sensors in smart structures and systems [12]. With the rapid development of nanoelectromechanicalsystem (NEMS) technologies, thin
135
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
films with the thickness of the film in the nanoscale have been extensively used as structural as well as electrical components of NEMS devices, especially ultrasensitive sensors for ultrafine resolution applications. Proper measurement of electromechanical coupling parameters of piezoelectric thin film accurately is essential to their applications [13]. It is also believed that indentation technique can be capable of probing the electrical properties of piezoelectric materials as the mechanical properties of thin film [14]. For the indentation responses of transversely isotropic piezoelectric thin film, there is a constant potential excited due to the piezoelectric effect, and the applied load will be affected not only by the substrate effect [9] but also by the excited electric potential during the indentation process [15]. Because of their inherent complexity, relatively few methods are available in literatures to evaluate the whole sets of the elastic parameters/piezoelectric constants. In this paper, the variety of the PZT thin film is determined by combining nanoindentation test with FEM simulation based on the weight averaged method, and the piezoelectric constitution is used to establish the supplemental equation, in which the piezoelectric strain constants are related with the piezoelectric stress constants, so that all of the piezoelectric constants of piezoelectric thin film can be determined by combining nanoindentation test with FEM. In the forward analysis, the numerical loading curves at the broad spectrum of possible material combinations are simulated by using ABAQUS software to extract the numerical maximum indentation loads and loading curve exponents, and they are used to establish two dimensionless equations related with the engineering/piezoelectric constants of film/substrate system. In the reverse analysis, the loading parts of experimental indentation curves performed on PZT thin film in nanoindentation tests are fitted as the power function to obtain the maximum indentation loads and the loading curve exponents, and the engineering constants can be solved by using the simultaneity of dimensionless and supplemental equations. In order to verify the validity, the optimized engineering constants are inputted into ABAQUS software to obtain the numerical loading curves and they are compared with the engineering constants of PZT thin film measured by the previous method to determine the variety of the PZT thin film. The piezoelectric constants of the certain variety of PZT thin film are determined by combining forward analysis and inverse analysis. The engineering constants of PZT thin film are determined and the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined and they are reasonable by comparing with the previous results. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. 2. Indentation model Generally, a piezoelectric thin film is perfectly bonded to a substrate and loaded by an axisymmetric rigid indenter as shown in Fig. 1, and the piezoelectric thin film is simplified as transversely isotropic [16]. The rectangular coordinate system is introduced, and the x1–x2 plane is assumed as the transversely isotropic plane
while the x3 axis is the rotational symmetry axis as shown in Fig. 1. Here, the subscripts 1 and 2 denote for T transverse directions while the subscript 3 denotes for L longitudinal direction. The film/substrate nanoindentation system is schematized in Fig. 1. The rigid axisymmetric indenter is a cone with the included semi-apex angle a and a perfect electric insulator with a known radial distribution of electric charge at its surface, and P, a and h are the indentation load, the radius of indentation contact region and the indentation depth. t is the thickness of the thin film. R and H are radius and thickness of substrate. The friction is usually weak during the indentation test therefore it is considered that there is no friction at the interface between the indenter and thin film [14]. 2.1. Governing equations In the absence of body and inertia forces, the equilibrium equations for the piezoelectric thin film in Fig. 1 can be expressed in the rectangular coordinate system as [17]
@ r11 @ s12 @ s31 þ þ ¼0 @x1 @x2 @x3 @ s12 @ r22 @ s23 þ þ ¼0 @x1 @x2 @x3 @ s31 @ s23 @ r33 þ þ ¼0 @x1 @x2 @x3
where rij (i, j = 1,2,3) is the normal stress, sij (i, j = 1, 2, 3) the shear stress. Correspondingly, the electrostatic equation without body charge is
@D1 @D2 @D3 þ þ ¼0 @x1 @x2 @x3
ð2Þ
where Di (i = 1, 2, 3) is the electric displacement component. The geometric conditions are related with the strains and displacements and they are expressed as follows
@u1 @x1 @u2 ¼ @x2 @u3 ¼ @x3
1 @u1 @u2 þ 2 @x2 @x1 1 @u1 @u3 ¼ þ 2 @x3 @x1 1 @u2 @u3 ¼ þ 2 @x3 @x2
e11 ¼
e12 ¼
e22
e13
e33
e23
ð3Þ
where eij (i, j = 1, 2, 3) is the strain, ui (i = 1, 2, 3) the displacement. The electric field can be written as
E1 ¼
@/ ; @x1
E2 ¼
@/ ; @x2
E3 ¼
@/ @x3
ð4Þ
where Ei (i = 1, 2, 3) is the electric field, / the electric potential. For transversely isotropic piezoelectric materials, the constitutive equations can be represented as [7]
ekl ¼ sklij rij þ dikl Ei ði; j ¼ 1; 2; 3; k; l ¼ 1; 2; 3Þ Dk ¼ dkij rij þ nki Ei
ð5Þ
where sklij is the elastic compliance constants, dkij the piezoelectric coefficients, nki the dielectric constants. The elastic constant matrix c is the inverse matrix of the full compliance matrix s [18]. For transverse isotropic materials, the elastic compliance matrix is [7,19]
2
Fig. 1. The schematic representation of film/substrate system in the selected coordinate system.
ð1Þ
1=ET
6 m =E T T 6 6 6 mTL =ET s¼6 6 0 6 6 4 0 0
3
mT =ET
mLT =EL
0
0
0
1=ET
mLT =EL
0
0
mTL =ET
1=EL
0
0
0
0
1=GT
0
0 0
0 0
0 0
1=GL 0
0 7 7 7 0 7 7 0 7 7 7 0 5 1=GL
ð6Þ
136
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
where ET is the transverse Young’s modulus, EL the longitudinal Young’s modulus, GT the transverse shear modulus, GL the longitudinal shear modulus, mT the Poisson’s ratio in the transversely isotropic (x1–x2) plane as shown in Fig. 1. mTL has the physical interpretation of the Poisson’s ratio that characterizes the strain in the x1–x2 plane resulting from stress normal to it, and mLT characterizes the longitudinal strain in the direction normal to the transversely isotropic plane resulting from stress in the x1–x2 plane. The Poisson’s ratios mTL and mLT are related as
mTL ET
¼
mLT
ð7Þ
EL
The transverse shear modulus GT can be expressed as
GT ¼
ET 2ð1 þ mT Þ
ð8Þ
The substrate is assumed as a linear elastic isotropic material characterized by Young’s modulus ES and the Poisson’s ratio mS [9].
Within the contact area of the conical indenter in Fig. 1, the downward displacement u2 is given as follow [14]
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 6 x21 þ x23 6 a ð9Þ
For the conical indenter, there is no friction between the indenter and the film, and it can be expressed as
s12 ¼ s23 ¼ 0; ðr P 0Þ
ð10Þ
There is no normal stress outside the contact region and it can be written as
r22 ¼ 0; ðr > aÞ
ð11Þ
The indenter is a perfect electrical insulator with a zero electric charge distribution. The electric charge distribution of the surface is zero and it can be given as
D2 ðx1 ; 0; x3 Þ ¼ 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 þ x23 > a
3. Evaluation of the engineering constants 3.1. Dimensionless analysis For an idealized transversely isotropic bulk material, the loading part of indentation curve from a conical indentation test obeys a quadratic relationship F = Ch2, where C is loading curvature, F indentation load, h indentation depth [20]. The quadratic relationship does not accurately describe the indentation response of thin film, so it is described by the power function F = Chx, in which the power exponent x is defined as the impact factor to reflect the substrate effect [7]. Then the relationship between indentation load F and indentation depth h can be rewritten as [7,21]
F ¼ F max
2.2. Boundary conditions
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 ðx1 ; 0; x3 Þ ¼ h x21 þ x23 cot a ;
and GL, and the simplified simulation results under PMI and PI modes can reflect the influence of the piezoelectric effect on characterization of engineering constants in the preliminary study.
ð12Þ
2.3. Simplification of the material model Generally, the mechanical and electrical properties of the transversely isotropic piezoelectric materials can be especially described by five independent engineering constants ET, EL, GL, mT and mTL, three piezoelectric coefficients e15, e31 and e33, and two dielectric constants n11 and n33[15]. As for the transversely isotropic materials without piezoelectric effect, the elastic properties can be only described by five independent engineering constants, and it is difficult to determine all the unknown engineering constants ET, EL, GL, mT and mTL [19]. A few of assumptions are conducted to simplify the material model, so the total number of independent mechanical variables should be reduced before determining the mechanical constants. When piezoelectric thin film is simplified as a kind of transversely isotropic materials, the Poisson’s ratio mT can be assumed as 0.3 [5]. The Poisson’s ratios mTL and mLT are related as mTL + mLT = 2mT [19]. When the piezoelectric and dielectric constants of the transversely isotropic thin film are taken as the known constants, the purely mechanical indentation (PMI) and piezoelectric indentation (PI) modes of ABAQUS software can be conveniently used to simulate the numerical loading curve in order to establish the dimensionless equations. Namely, the engineering constants of the transversely isotropic thin film can be characterized by the independent unknown engineering constants ET, EL
h hm
x ð13Þ
where Fmax and hm are the maximum indentation load and the corresponding maximum penetration. According to the previous results of transversely isotropic thin film indentation [22], the maximum indentation load Fmax and the loading curve exponent x should be the functions of the material constants, such as the engineering constants ET, EL and GL, the substrate elastic parameter ES, the geometric parameter of film thickness t, and the maximum indentation depth hm. Namely, they can be expressed as follows
F max ¼ /0 ðET ; EL ; GL ; ES ; t; hm Þ
ð14Þ
x ¼ u0 ðET ; EL ; GL ; ES ; t; hm Þ
ð15Þ
It is easy to know the parameters ES, t and hm by nanoindentation technique, therefore they are regarded as known constants in order to simulate the numerical loading curve, in which the numerical maximum indentation load and numerical loading curve exponent can be extracted. Applying Buckingham Pi Theorem to Eqs. (14) and (15), they can be written as dimensionless form
ET EL GL ¼ / ; ; 2 ES ES ES E S hm ET EL GL x¼u ; ; ES ES ES F max
ð16Þ ð17Þ
The elastic properties of transversely isotropic thin film can be described by the independent unknown engineering constants ET, EL and GL. In addition, the indentation modulus Em is the weighted average of the engineering constants of anisotropic materials along the different orientations and it is defined as [7,19]
Em ¼
EL þ ET 2
ð18Þ
Generally the unloading part of the indentation curve can reflect factual and comprehensive information of indented materials, such as plastic deformation, the expansion of surface and subsurface cracks during the indentation loading process, therefore Em is traditionally determined by the unloading part slope of indentation curve at the maximum indentation load according to the wellestablished Oliver and Pharr method. The elastic and plastic deformation processes are gained by the analysis of the indentation curve, and the widely studied parameters are the residual depth of penetration (hres), total work of contact deformation (Wtot), work from elastic recovery (Wel) and the work of inelastic deformation (Wpl). The numerical values of the quantities Wtot, Wel and Wpl can be directly obtained from the indentation curve, as the area under the loading curve, under the unloading curve and the area enclosed by the two curves, respectively. The residual depth ratio
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
hres/hm and plastic work ratio Wpl/Wtot characterize the extent of plastic deformation and thus are examined with respect to the deformation behavior of the materials. Although it is difficult to simulate the unloading part of the indentation curve by using commercial package, one can easily obtain the numerical loading curve under the simplified idealization through ABAQUS. As we all know the indentation curve can reflect ideal elasticity at low indentation load only, therefore we choose the simulation on loading curve in order to determine the elastic constants. However the surface and subsurface cracks may be expansive because of the maximum indentation loads at the beginning of the unloading process, and it indicates that it is very difficult to simulate the unloading part of the indentation. At the appropriate penetration depth the numerical maximum indentation loads and numerical loading curve exponents can be extracted from the numerical loading curves of the thin film to establish the dimensionless Eqs. (16) and (17), and the engineering constants EL, ET and GL can be determined by solving Eqs. (16)–(18). 3.2. Forward analysis 3.2.1. Finite element procedure Numerical indentation tests performed on the film/substrate nanoindentation system schematized as Fig. 1 are simulated by
137
using commercial numerical finite element package ABAQUS, and we adopt PMI and PI modes to obtain the numerical loading curves. For Berkovich indenter, the corresponding apex angle a of the equivalent cone was chosen as 70.3° [23]. The numerical indentation is modeled as a two-dimensional axisymmetric problem. A fine finite element mesh is used in the vicinity of the contact region to capture the localized deformation underneath the indenter and a gradually coarser mesh away from the contact region is used to cut down the computational work. As shown as Fig. 2, the substrate is meshed by using four-node reduced integrated linear axisymmetric solid element (CAX4R). The film is meshed by using CAX4R element for PMI mode and four-node reduced integrated linear axisymmetric piezoelectric element (CAX4E) for PI mode [16]. Eighty elements are used through the thin film thickness. Gradually, the influence of friction is usually weak on the nanoindentation test, and there is no friction at the interface between the indenter and the thin film [7,14]. The numerical indentation load is recorded as a function of the vertical displacement beneath the tip of the indenter. In order to consider the material constants accurately, a broad spectrum of possible material constants combinations 0.1 6 EL/ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 0.1 6 GL/ES 6 2.0 are chosen as input parmeters to simulate the indentation process, and then the numerical loading curves can be obtained. When the engineering constants of thin film ET, EL, GL and mT, the substrate
Fig. 2. Finite element mesh of an axisymmetric indentation.
138
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
elastic parameters ES and mS, and the maximum indentation depth hm are regarded as input parameters, the numerical loading curve in situ reflects the relationship between the indentation force on the rigid indenter and the indenter penetration into the thin film during FEM simulation. Generally, a well-known empirical rule is that the indentation depth should be smaller than one tenth of film thickness at micrometer to avoid substrate effect [21], however it is too shallow to extract the accurate materials properties of extreme thin film at nanometer. At first, it is very difficult to determine the contact area and depth accurately when the contact depth is of the same order as the surface roughness. Secondly, there is a little distinguishing with indentation force and adhesion between the indenter tip and specimen, in order that the latter may not be properly measured based on the present nanoindentation technique. Finally, even a brand new indenter tip is blunt at the nanoscale, with a radius of curvature typically tens of nanometers; the conventional sharp indentation analysis cannot be applied directly to the blunt tip at a small-scale. In a word, the indentation depth is generally larger than one tenth of the film thickness to obtain the exact properties [9,16]. On the other hand, the indentation depth must be limited for different materials to avoid excessive substrate deformation and cracking, and for simplicity it is chosen as one fifth of the film thickness as the previous report [7]. 3.2.2. Numerical loading curves According to the previous results on the materials constants, the broad spectrum of possible material combinations 0.1 6 EL/ ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 0.1 6 GL/ES 6 2.0 are chosen as input parameters to simulate the nanoindentation test, and von-Mises stress in contour and the numerical loading curves can be obtained. The former can indicate whether the indentation process is disturbed by the substrate, while the latter will give an insight into how the various engineering constants of thin film influence the shape of the numerical loading curves. For brevity, the different combinations of normalized engineering constants are fixed as EL/ ES = (0.6, 0.8, 1.0, 1.2), ET/ES = (0.5, 0.7, 0.8, 1.0) and GL/ES = (0.2, 0.3, 0.4, 0.5), and the indentation depth is chosen as one fifth of the film thickness as the previous report [7]. The typical von-Mises stress in contour at EL/ES = 1.0, ET/ES = 0.8 and GL/ES = 0.3 is given in Fig. 3.
Obviously, the stress decreases continuously away from the contact region in the thin film while it is discontinuousness at the interface of film/substrate, indicating that the indentation process is disturbed by the substrate. The typical numerical loading curves at ET/ES = 0.8 and GL/ ES = 0.3 can be obtained under PMI mode, and it is given in Fig. 4. 2 Both of the indentation load F=ES hm as a function of the indentation depth h/hm and the curvature of the numerical loading curve increase with the increase of EL/ES. At the fixed ratios of ET/ES and 2 GL/ES, the normalized maximum indentation load F max =ES hm increases with the ratio EL/ES. They give an insight into how the various engineering constants of thin film influence the shape of the numerical loading curves. According to the numerical loading curves in Fig. 4, the numer2 ical maximum indentation load F max =ES hm and the loading curve exponent x can be obtained after the numerical loading curves are fitted as the power function F = Fmax(h/hm)x with Origin soft2 ware. The F max =ES hm and x versus EL/ES and GL/ES are respectively plotted in Fig. 5a–f. With the increase of ET/ES, EL/ES and GL/ES, 2 F max =ES hm increases while x decreases. One can conclude that the 2 F max =ES hm and x regularly vary with the engineering constants, and the dimensionless functions given as Eqs. (16) and (17) will be established by fitting the engineering constants vs numerical
Fig. 4. The typical numerical loading curves obtained from FEM indentation test with ET/ES = 0.8 and GL/ES = 0.3 (a) for the PMI mode and (b) for the PI mode.
Fig. 3. The representative von–Mises stress in contour of numerical indentation.
139
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Fig. 5. Given the representative combinations of two normalized engineering constants, the curves of the other engineering constant related with the normalized maximum 2 load F max =ES hm (a)–(c) and the loading curve exponent x (d)–(f).
Table 1 Piezoelectric and dielectric constants of PZT-6B and PZT-4. Dielectric constants (1010 Fm1)
Materials
Piezoelectric constants (Cm2) e31
e33
e15
j11
j33
PZT-4 PZT-6B
13.44 4.6
6.98 0.9
13.84 7.1
60.0 36.0
54.7 34.0
2
maximum load F max =ES hm curves shown in Fig. 5(a)-(c) and the engineering constants vs numerical loading curve exponent x curves shown in Fig. 5d–e, which are used to solve the unknown engineering constants in the reverse analysis as the following.
3.2.3. Establish of dimensionless functions A plenty of combinations of the numerical maximum load 2 F max =ES hm and numerical loading curve exponent x extracted from the numerical loading curves will be used to determine the specific
form of dimensionless functions. According to Eqs. (16) and (17), letting g EL/ES, n ET/ES and c GL/ES, the dimensionless functions can be written as
/e ðn; g; cÞ ¼ p1 þ p2 n þ ðp3 þ p4 nÞc þ ðp5 þ p6 nÞg þ ðp7 þ p8 nÞcg þ ðp9 þ p10 nÞg2 þ ðp11 þ p12 nÞcg2
ð19Þ
ue ðn; g; cÞ ¼ q1 þ q2 n þ ðq3 þ q4 nÞg þ ðq5 þ q6 nÞc þ ðq7 þ q8 nÞgc þ ðq9 þ q10 nÞc2 þ ðq11 þ q12 nÞgc2
ð20Þ
140
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
The coefficients pi, qi (i = 1, 2, . . . , 12) fitted are listed in Table 2 for PMI mode and Table 3 for PI mode. 3.3. Reverse analysis In the reverse analysis, the nanoindentation tests were performed on PZT thin film laid on a silicon substrate. The experimental indentation curves, the maximum indentation load Pmax, the maximum indentation depth hm and the indentation modulus Em can be automatically obtained by the nanoindenter. According to the loading part of the indentation curve the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The reverse analysis algorithm is schematic as the flow chart in Fig. 6, and it includes the extraction of the experimental data, the setting on the combination range of unknown engineering constants, the calculation on the engineering constant errors, the selection combination of the smallest total error, and the output of the engineering constants. After determining the solutions of the engineering constants, the method validity is discussed. 3.3.1. Reverse analysis algorithm According to Fig. 6, at first ES is substrate elastic parameter, and Pmax, hm, xm and Em can be obtained from the nanoindentation test. When ES, Pmax, hm, xm and Em are known, the engineering constants ET, EL and GL can be determined by solving Eqs. (18)–(20). Secondly, the combination ranges of unknown engineering constants are assumed within the wide ranges, and the engineering constants EL/ES, ET/ES, GL/ES are setting as 0.1 6 EL/ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 2 0.1 6 GL/ES 6 2.0. Thirdly, the errors of Pmax =ES hm , xm and Em can be calculated, and the total error is calculated as the sum of absolute values. Finally, the combination of engineering constants is selected as the ultimate solution with the smallest total error, and the values of ET, EL, GL are outputted. Table 2 The fitting coefficients in Eqs. (19) and (20) for PMI mode. Coefficients in Eq. (19)
Coefficients in Eq. (20)
p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12
0.21389 0.41697 0.79622 0.58033 1.44269 0.59941 12.2136 6.05127 4.06122 5.47017 1.80912 10.4618
2.25659 0.0372 0.04018 0.06076 1.19663 0.26827 1.75892 1.11014 0.95761 0.54791 0.34978 1.98759
Table 3 The fitting coefficients in Eqs. (19) and (20) for PI mode. Coefficients in Eq. (19)
p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
Coefficients in Eq. (20)
PZT-4
PZT-6B
0.57831 0.14498 0.13407 0.15477 1.36684 0.08392 1.12744 0.47344 0.41528 0.07819 0.36699 0.13658
0.21943 0.13188 0.25928 0.37999 2.59463 0.95794 3.85829 0.99306 0.72688 0.25946 2.66842 1.09466
q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12
PZT-4
PZT-6B
2.28813 0.05121 0.09849 0.04084 0.04187 0.21079 0.01569 0.08768 0.58711 0.07475 0.06713 0.11082
2.01627 0.04739 0.04061 0.29019 0.37009 0.09921 0.42868 0.04421 0.56598 0.18855 0.03988 0.08139
Fig. 6. The flow chart of the reverse analysis for evaluating engineering constants.
3.3.2. Application to PZT thin film Using the nanoindenter (Hysitron Triboindenter, Hysitron Inc., USA), the nanoindentation test was performed on PZT thin film deposited on silicon substrate, and the maximum indentation depth was set as one fifth of film thickness, i.e. 70 nm. five indents were performed during the nanoindentation test, and the five experimental curves ‘‘1–5’’ and one averaged experimental curve ‘‘a’’ are given in Fig. 7a, and AFM image of nanoindentation impression is described as Fig. 7b. The experiment data Pmax and xm are determined as 1.3743 mN and 2.0336 from the averaged experimental curve in Fig. 7a, and the indentation modulus Em is obtained as 139.4 GPa by using the well-established Oliver and Pharr [2] method. According to the loading part of the experimental indentation curves, the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The Young’s modulus of silicon substrate is assumed to be known as ES = 130 GPa for simplicity according to the experimental report [24]. The related experimental data is listed in Table 4. 3.3.3. Determination of the engineering constants 3.3.3.1. Solutions for PMI mode. Substituting Pmax, hm, ES and xm into Eqs. (18)–(20), the multiple solutions will be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the engineering constants are selected as the solutions given in Table 5. In order to verify the validity of the solutions, the solutions in Table 5 are taken as inputted parameters of ABAQUS software, and the corresponding numerical loading curves S1, S2, and S3 for the PMI mode are obtained and shown in Fig. 8. One can see that the numerical loading curves are in good agreement with the averaged experimental loading curve, and the numerical loading curve S3 with the smallest total error is closest to the averaged experimental loading
141
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Fig. 7. (a) The five experimental load versus depth curves and the averaged load versus depth curve and (b) AFM image of indentation impression on PZT thin film.
Table 4 The related experimental data used in reverse analysis.
Table 6 Some representative engineering constants of PZT thin film [25].
t (nm)
hm (nm)
Fmax (mN)
exponent x
d33 (pC/ N)
Em (GPa)
ES (GPa)
350
70
1.3743
2.0336
60
139.4
130
Table 5 The multiple solutions for the PMI mode. Solutions
etotal (%)
ET (GPa)
EL (GPa)
GL (GPa)
S1 S2 S3
2.385 1.361 1.049
157 154 152
124 125 128
65 64 61
Materials
ET (GPa)
EL (GPa)
GL (GPa)
PZT-4 PZT-2 PZT-6B
81.3 86.2 137.0
64.1 67.5 131.6
56.8 33.4 54.1
of PZT-6B and PZT-4 are chosen as the input parameters listed in Table 1. The multiple solutions can be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the engineering constants are selected as the solutions given in Table 7. In order to verify the validity of the solutions, the multiple solutions are taken as inputted parameters of ABAQUS software to obtain the numerical loading curves S1, S2, and S3 shown in Fig. 9. Obviously, the total error of S3 is smallest, and S3 is taken as the ultimate engineering constants for the PZT thin film under the PI mode of PZT-6B. Comparing with the engineering constants of PZT-6B [25], the relative differences of ET, EL and GL are 3.65%, 3.94% and 4.04%, indicating that the PZT sample is identified as PZT-6B. 4. Evaluation of the piezoelectric constants The variety of PZT thin film has been determined as PZT-6B by evaluating the engineering constants, and in this section, the piezoelectric constants of PZT-6B will be evaluated by combining FEM simulation with the nanoindentation test. 4.1. Dimensionless equations
Fig. 8. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 5 for the PMI mode.
curve. The solution corresponding to the closest numerical loading curve S3 is taken as the ultimate engineering constants to compare with the previous results, such as PZT-6B, PZT-2 and PZT-4 listed in Table 6 [25]. The relative differences of ET, EL and GL are 43.3%, 47.3%, 45.2% for PZT-2, 59.7%, 58.4%, 65.4% for PZT-4 and 9.9%, 2.8%, 11.3% for PZT-6B. Obviously, the relative differences are least for PZT-6B therefore the PZT thin film is primarily determined as PZT-6B. 3.3.3.2. Solutions for PI mode. In order to identify the variety of PZT thin film further, the piezoelectric effect is considered and the engineering constants can be solve for PI mode. According to the results under PMI mode, the piezoelectric and dielectric constants
The engineering constants and the dielectric constants of the PZT-6B are known parameters. According to the previous results of transversely isotropic piezoelectric thin film indentation [7], the maximum indentation load Fmax and the loading curve exponent x should be the functions of the material parameters, such as the piezoelectric constants e15, e31, and e33, the substrate elastic parameter ES, the geometric parameter of film thickness t, and the maximum indentation depth hm. Namely, they can be expressed as follows
F max ¼ /p ðe15 ; e31 ; e33 ; ES ; t; hm Þ
ð21Þ
x ¼ up ðe15 ; e31 ; e33 ; ES ; t; hm Þ
ð22Þ
The moderate ratio of penetration depth to the thin film thickness is controlled as one fifth, and then applying the Buckingham PI Theorem Eqs. (21) and (22) can be written as
142
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Table 7 The multiple solutions for PI mode. Solutions
S1 S2 S3
PZT-6B
PZT-4
etotal (%)
ET (GPa)
EL (GPa)
GL (GPa)
etotal (%)
EL (GPa)
GL (GPa)
ET (GPa)
0.9 0.7 0.3
149 146 142
131 134 137
54 53 52
1.922 1.891 0.816
150 146 140
130 135 139
55 52 50
Fig. 9. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 7 for the PI mode.
e15 e31 ; 2 e33 e33 ES hm e15 e31 x ¼ up ; e33 e33 F max
¼ /p
ð23Þ ð24Þ
Based on the piezoelectric constitution, the longitudinal direction piezoelectric strain constant d33 can be written as [26]
d33 ¼ 2e31 s13 þ e33 s33
ð25Þ
where s13 and s33 are known according to Eq. (6). When d33 is the known parameter, the unknown piezoelectric constants e15, e31, and e33 will be determined by solving Eqs. (23)–(25). 4.2. Forward analysis 4.2.1. Numerical loading curves Numerical piezoelectric indentation responses are simulated using CAX4E elements. In order to cover a broad range of materials, the broad spectrum of possible material combinations 0 < e15/ e33 < 3.0 and 3.0 < e31/e33 < 0 are chosen as input parameters to obtain the numerical loading curves. For brevity, the different combinations of normalized piezoelectric constants are fixed as e15/ e33 = (0.1, 0.5, 1.0, 2.2) and e31/e33 = (0.01, 0.15, 1.4, 3.0), and they are used to obtain the numerical loading curves for evaluating piezoelectric constants. To evaluate the piezoelectric constants, the engineering constants and dielectric constants of PZT-6B thin film [25] are chosen as the input parameters, and the typical numerical loading curves are given in Fig. 10 under the combination e31/e33 = 1.4 and e15/e33 = (0.1, 0.5, 1.0, 2.2). In the figure, the indentation load increases with piezoelectric constants e15/e33 at the same indentation depth, which is attributed to the piezoelectric effect on the indentation responses. When the insulated indenter with zero distributed electric charge is pressed into the piezoelectric materials, the indenter creates a constant electric potential inside the contact area [15]. The electric field induced during indentation as a result of the electrical– mechanical coupling can resist or aid in the penetration of the indenter into the piezoelectric material depending on the electrical
Fig. 10. The typical numerical loading curves with different e15/e33 at e31/e33 = 1.4.
2
Fig. 11. The normalized maximum indentation load F max =ES hm as the function of the piezoelectric constant ratios e31/e33 at different ratios of e15/e33.
conductivity of the indenter and the surface boundary conditions of the indented material. After the numerical loading curves under the above four combinations are fitted as the power function F = Fmax(h/hm)x by Origin
143
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144 Table 8 The fitting coefficients in Eqs. (28) and (29). Coefficients in Eq. (28)
Coefficients in Eq. (29)
m1 m2 m3 m4 m5 m6 m7 m8 m9
n1 n2 n3 n4 n5 n6 n7 n8 n9
1.98558 0.12324 0.04294 0.02850 5.1170 0.00283 0.00178 0.00184 9.3870
2.06537 0.04794 8.8940 2.1750 0.01095 0.00397 6.8920 0.0020 7.7940
Fig. 12. The loading curve exponent x as the function of the piezoelectric constant ratios e31/e33 at different ratios of e15/e33.
software, the maximum indentation loads and the loading curve exponents can be obtained via extensive FEM simulations, the nor2 malized maximum indentation load F max =ES hm and loading curve exponent x as the function of the piezoelectric constant ratio e31/ e33 are given in Figs. 11 and 12 at different ratios of e15/e33. Obviously, with the increase of e31/e33, the normalized maximum 2 indentation load F max =ES hm decreases while the loading curve exponent x increases at the fixed e15/e33. The normalized maxi2 mum indentation load F max =ES hm and the loading curve exponent x continuously and regularly vary with e15/e33 and e31/e33, and the dimensionless functions given as Eqs. (23) and (24) will be established by fitting the piezoelectric constants vs numerical 2 maximum load F max =ES hm curve shown in Fig. 11 and the piezoelectric constants vs numerical loading curve exponent x curve shown in Fig. 12, which are used to solve the unknown piezoelectric constants in the reverse analysis as the following. 4.2.2. Dimensionless functions Similarly, according to Eqs. (23) and (24), letting a e15/e33 and b e31/e33 the dimensionless functions can be written as
F max 2
ES hm
¼ /p ða; bÞ
x ¼ up ða; bÞ
ð26Þ ð27Þ
The formulas of the dimensionless functions (26) and (27) are fitted by Origin software as follows
Fig. 13. The flow chart of the reverse analysis for evaluating the piezoelectric constants.
/p ða;bÞ ¼ m1 þm2 a þm3 a2 þðm4 þm5 a þm6 a2 Þbþðm7 þm8 a þm9 a2 Þb2 ð28Þ
up ða;bÞ ¼ n1 þn2 a þn3 a2 þðn4 þn5 a þn6 a2 Þbþðn7 þn8 a þn9 a2 Þb2
ð29Þ
The coefficients mi, ni (i = 1, 2, . . . , 9) in Eqs. (26) and (27) are listed in Table 8. 4.3. Reverse analysis In the reverse analysis, the nanoindentation tests were performed on PZT-6B thin film laid on a silicon substrate. The experimental indentation curves, the maximum indentation load Pmax and the maximum indentation depth hm can be automatically obtained by the nanoindenter. According to the loading part of the indentation curve the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The reverse analysis algorithm is schematic as the flow chart in Fig. 13, and it includes the extraction of the experimental data, the setting on the combination range of unknown piezoelectric constants, the calculation on the piezoelectric constants errors, the selection combination of the smallest
total error, and the output of the piezoelectric constants. After determining the solutions of the piezoelectric constants, the method validity is discussed. 4.3.1. Reverse analysis algorithm At first ES is substrate elastic parameter, and Pmax, hm and xm can be obtained from the nanoindentation test. The piezoelectric strain constant d33 can be determined by a PFM experiment. When Pmax, hm, xm and d33 are known, the piezoelectric constants e15, e31, and e33 can be determined by solving Eqs. (25)–(27). Secondly, the combination ranges of unknown engineering constants are assumed within the wide ranges, and the normalized piezoelectric constants e15/e33 and e31/e33 are setting as 0 < e15/e33 < 3.0 and 2 3.0 < e31/e33 < 0. Thirdly, the errors of P max =ES hm , xm and d33 can be calculated, and the total error is calculated as the sum of absolute values. Finally, the combination of engineering constants is selected as the ultimate solution with the smallest total error, and the values of e15, e31, e33 are outputted.
144
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
differences are 0.3–5.6% indicating that the piezoelectric constants are reasonable.
Table 9 Multiple solutions of piezoelectric constants. Solution
etotal (%)
e33 (C/m2)
e31 (C/m2)
e15 (C/m2)
S1 S2 S3
0.10 0.08 0.07
4.60 4.70 4.75
1.30 1.00 0.85
7.68 7.26 7.08
Fig. 14. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 9.
4.3.2. Application to PZT thin film Using the nanoindenter (Hysitron Triboindenter, Hysitron Inc., USA), the nanoindentation test was performed on PZT thin film deposited on silicon substrate, and the maximum indentation depth was set as one fifth of film thickness, i.e. 70 nm. five indents were performed during the nanoindentation test, and the five experimental curves ‘‘1–5’’ and one averaged experimental curve ‘‘a’’ are given in Fig. 7a, and AFM image of nanoindentation impression is described as Fig. 7b. The experiment data Pmax and xm are determined as 1.3743 mN and 2.0336 from the averaged experimental curve in Fig. 7a, and the piezoelectric strain constant d33 is obtained as d33 = 60 pC/N according to the experimental report [27]. The Young’s modulus of silicon substrate is assumed to be known as ES = 130 GPa for simplicity according to the experimental report [24].
4.3.3. Determination of the piezoelectric constants Substituting Pmax, hm, ES, xm and d33 into Eqs. (25)–(27), the multiple solutions will be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the piezoelectric constants are selected as the solutions given in Table 9. In order to verify the validity of the solutions, the solutions in Table 9 are taken as inputted parameters of ABAQUS software, and the corresponding numerical loading curves S1, S2, and S3 are obtained and shown in Fig. 14. One can see that the numerical loading curves are in good agreement with the averaged experimental loading curve, and the numerical loading curve S3 with the smallest total error is closest to the averaged experimental loading curve. The solution corresponding to the closest numerical loading curve S3 is taken as the ultimate piezoelectric constants of the PZT thin film. Similarly, comparing with the previous results e15 = 7.1, e31 = 0.9 and e33 = 4.6 [28], the relative
5. Conclusion Combining nanoindentation test with FEM simulation, a method is proposed to identify the variety of PZT thin film via evaluating engineering constants based on the weight average method, and then the piezoelectric constants e15, e31 and e33 of the PZT thin film are determined by combining forward and inverse analysis. The numerical loading curves are obtained from the extensive FEM simulations at the moderately deep indentation for the PMI and PI modes, and the nanoindentation tests were performed on PZT thin film to extract the maximum indentation load, the loading curve exponent and the maximum indentation depth. The engineering constants of PZT thin film are determined as ET = 142 GPa, EL = 137 GPa and GL = 52 GPa, and the variety of the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined as e15 = 4.75 C/m2, e31 = 0.85 C/m2 and e33 = 7.08 C/m2 based on piezoelectric constitution by combining forward analysis and inverse analysis. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. Acknowledgments This work was supported by PCSIRT (IRT1080), NNSF of China (10825209), Changjiang Scholar Incentive Program ([2009]17), Aid Program for Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan Province, and Shanghai Nano Special Foundation (11nm0502600). References [1] M.K. Khan, S.V. Hainsworth, M.E. Fitzpatrick, L. Edward, Comput. Mater. Sci. 49 (2010) 751–760. [2] W.C. Oliver, G.M. Pharr, J. Mater. Res. 7 (1992) 1564–1583. [3] S. Liu, Y.T. Gu, H. Huang, Mater. Sci. Eng. A 528 (2011) 7948–7951. [4] M.R. VanLandingham, J. Res. Natl. Inst. Stand Technol. 108 (2003) 249–265. [5] M. Bocciarelli, G. Bolzon, G. Maier, Mech. Mater. 37 (2005) 855–868. [6] R. Saha, W. Nix, Acta Mater. 50 (2002) 23–38. [7] J.S. Wang, X.J. Zheng, H. Zheng, et al., Comput. Mater. Sci. 49 (2010) 378–385. [8] F.Q. Yang, Mater. Sci. Eng. A 358 (2003) 226–232. [9] M.H. Zhao, X. Chen, Y. Xiang, J.J. Vlassk, et al., Acta Mater. 55 (2007) 6260– 6274. [10] D.J. Ma, K.W. Xu, J.W. He, et al., Surf. Coat. Technol. 116 (1999) 128–132. [11] Y.G. Liao, Y.C. Zhou, Y.L. Huang, et al., Mech. Mater. 41 (2009) 308–318. [12] K. Uchino, Piezoelectric Actuators and Ultrasonic Motors, Kluwer Academic Publishers, Boston, Dordrecht, London, 1997. [13] W.N. Sharpe, B. Yuan, R.L. Edwards, J. Microelectromech. Syst. 6 (1997) 193– 199. [14] J.H. Wang, C.Q. Chen, T.J. Lu, J. Mech. Phys. Solids 56 (2008) 3331–3351. [15] A.E. Giannakopoulos, S. Suresh, Acta Mater. 47 (1999) 2153–2164. [16] S.M. Han, Saha R, W.D. Nix, Acta Mater. 54 (2006) 1571–1581. [17] Z.K. Wang, B.L. Zheng, Int. J. Solids Struct. 32 (1995) 105–115. [18] S. Lee, B.C. Cho, H.C. Park, et al., Mater. Chem. Phys. 75 (2002) 174–177. [19] T. Nakamura, Y. Gu, Mech. Mater. 39 (2007) 340–356. [20] B. Poon, D. Rittel, G. Ravichandran, Int. J. Solids Struct. 45 (2008) 6399–6415. [21] H. Zheng, X.J. Zheng, J.S. Wang, et al., Mater. Des. 32 (2011) 1407–1413. [22] C. Gamonpilas, E.P. Busso, Mater. Sci. Eng. A 380 (2004) 52–61. [23] M. Dao, N. Chollacoop, K.J. Van Vlient, et al., Acta Mater. 49 (2001) 3899–3918. [24] J.J. Wortman, R.A. Evans, J. Appl. Phys. 36 (1965) 153–156. [25] M.H. Zhao, Y.P. Shen, Y.J. Liu, et al., Theor. Appl. Fract. Mech. 26 (1997) 141– 149. [26] D. Damjanovic, Rep. Prog. Phys. 61 (1998) 1267–1324. [27] S. Hiboux, P. Muralt, N. Setter, Mater. Res. Soc. 596 (2000) 499–504. [28] M.H. Zhao, N. Li, C.Y. Fan, Eng. Anal. Bound Elem. 33 (2008) 545–555.
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Evaluation of engineering/piezoelectric constants of piezoelectric thin film by combining nanoindentation test with FEM S.T. Song a, X.J. Zheng a,b,⇑, H. Zheng a, W. Liu a a b
Faculty of Materials, Optoelectronics and Physics, Xiangtan University, Xiangtan, Hunan 411105, PR China School of Materials Science and Engineering, University of Shanghai for Science & Technology, Shanghai 200093, PR China
a r t i c l e
i n f o
Article history: Received 21 February 2012 Received in revised form 23 March 2012 Accepted 1 May 2012 Available online 30 June 2012 Keywords: Engineering constants Piezoelectric constants Finite element method Transversely isotropic Thin film
a b s t r a c t Based on the weight averaged method, the variety of the PZT thin film is determined by combining nanoindentation test with finite element method (FEM) simulation, and the piezoelectric constitution is used to establish the supplemental equation, in which the piezoelectric strain constants are related with the piezoelectric stress constants, so that all of the piezoelectric constants of piezoelectric thin film can be determined by combining nanoindentation test and FEM. In the forward analysis, the numerical loading curves at the broad spectrum of possible material combinations are simulated by using ABAQUS software to extract the numerical maximum indentation loads and loading curve exponents, and they are used to establish two dimensionless equations related with the engineering/piezoelectric constants of film/substrate system. In the reverse analysis, the loading part of experimental indentation curves performed on PZT thin film in nanoindentation test are fitted as the power function to obtain the maximum indentation loads and the loading curve exponents, and the engineering constants can be solved by using the simultaneity of dimensionless and supplemental equations. In order to verify the validity, the optimized engineering constants are inputted into ABAQUS software to obtain the numerical loading curves, and they are compared with the engineering constants of PZT thin film measured by the previous method to determine the variety of the PZT thin film. The piezoelectric constants of the certain variety are determined by combining forward analysis and inverse analysis. The engineering constants of the PZT thin film are determined as ET = 142 GPa, EL = 137 GPa, GL = 52 GPa, and the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined as e15 = 4.75 C/m2, e31 = 0.85 C/m2, e33 = 7.08 C/m2. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction The instrumented nanoindentation systems have been used to study the mechanical properties of thin film deposited on substrate effectively, and allow the application of a specified force or depth history [1–3], such that indentation load and the depth, are controlled and/or measured simultaneously and continuously over a complete loading cycle [4]. In general, the elastic modulus of isotropic thin film can be obtained from the indentation curve according to the method proposed by Oliver and Pharr [2]. The approaches cannot calibrate the whole set of the elastic parameters in anisotropic elasticity, even if indentations are performed in multiple directions and the shape of the contact area is exploited as supplementary information [5]. In nanoindentation test, the contact depth is generally less than one tenth of film thickness in order to obtain ⇑ Corresponding author at: Faculty of Materials, Optoelectronics and Physics, Xiangtan University, Xiangtan, Hunan 411105, PR China. Tel.: +86 731 58293648; fax: +86 731 58298119. E-mail address: [email protected] (X.J. Zheng). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.05.001
the intrinsic properties and avoid the so-called substrate effect [6,7]. However, the rule limiting the indentation depth may be only applicable to relatively thick films, and it is often not practical with the decrease of film thickness in many applications [8]. When the indentation depth is too shallow, measure error is induced by the small-scale effects, such as the adhesion effect between indenter tip and film, the strain gradient effect and the surface roughness effect [9]. Ma et al. [10] and Zhao et al. [9] considered the substrate effect during finite element method (FEM) simulation and nanoindentation test, and they proposed a method via combining nanoindentation with FEM to extract Young’s modulus and yield strength of isotropic thin films [11]. In our previous work, we have identified the elastic parameters of ZnO thin film when the substrate effect is considered for the transversely isotropic thin film [7]. Piezoelectric materials are a type of functional materials with mechanical and electrical coupling behaviors that are attracting increasingly wide engineering applications, such as actuators, sensors in smart structures and systems [12]. With the rapid development of nanoelectromechanicalsystem (NEMS) technologies, thin
135
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
films with the thickness of the film in the nanoscale have been extensively used as structural as well as electrical components of NEMS devices, especially ultrasensitive sensors for ultrafine resolution applications. Proper measurement of electromechanical coupling parameters of piezoelectric thin film accurately is essential to their applications [13]. It is also believed that indentation technique can be capable of probing the electrical properties of piezoelectric materials as the mechanical properties of thin film [14]. For the indentation responses of transversely isotropic piezoelectric thin film, there is a constant potential excited due to the piezoelectric effect, and the applied load will be affected not only by the substrate effect [9] but also by the excited electric potential during the indentation process [15]. Because of their inherent complexity, relatively few methods are available in literatures to evaluate the whole sets of the elastic parameters/piezoelectric constants. In this paper, the variety of the PZT thin film is determined by combining nanoindentation test with FEM simulation based on the weight averaged method, and the piezoelectric constitution is used to establish the supplemental equation, in which the piezoelectric strain constants are related with the piezoelectric stress constants, so that all of the piezoelectric constants of piezoelectric thin film can be determined by combining nanoindentation test with FEM. In the forward analysis, the numerical loading curves at the broad spectrum of possible material combinations are simulated by using ABAQUS software to extract the numerical maximum indentation loads and loading curve exponents, and they are used to establish two dimensionless equations related with the engineering/piezoelectric constants of film/substrate system. In the reverse analysis, the loading parts of experimental indentation curves performed on PZT thin film in nanoindentation tests are fitted as the power function to obtain the maximum indentation loads and the loading curve exponents, and the engineering constants can be solved by using the simultaneity of dimensionless and supplemental equations. In order to verify the validity, the optimized engineering constants are inputted into ABAQUS software to obtain the numerical loading curves and they are compared with the engineering constants of PZT thin film measured by the previous method to determine the variety of the PZT thin film. The piezoelectric constants of the certain variety of PZT thin film are determined by combining forward analysis and inverse analysis. The engineering constants of PZT thin film are determined and the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined and they are reasonable by comparing with the previous results. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. 2. Indentation model Generally, a piezoelectric thin film is perfectly bonded to a substrate and loaded by an axisymmetric rigid indenter as shown in Fig. 1, and the piezoelectric thin film is simplified as transversely isotropic [16]. The rectangular coordinate system is introduced, and the x1–x2 plane is assumed as the transversely isotropic plane
while the x3 axis is the rotational symmetry axis as shown in Fig. 1. Here, the subscripts 1 and 2 denote for T transverse directions while the subscript 3 denotes for L longitudinal direction. The film/substrate nanoindentation system is schematized in Fig. 1. The rigid axisymmetric indenter is a cone with the included semi-apex angle a and a perfect electric insulator with a known radial distribution of electric charge at its surface, and P, a and h are the indentation load, the radius of indentation contact region and the indentation depth. t is the thickness of the thin film. R and H are radius and thickness of substrate. The friction is usually weak during the indentation test therefore it is considered that there is no friction at the interface between the indenter and thin film [14]. 2.1. Governing equations In the absence of body and inertia forces, the equilibrium equations for the piezoelectric thin film in Fig. 1 can be expressed in the rectangular coordinate system as [17]
@ r11 @ s12 @ s31 þ þ ¼0 @x1 @x2 @x3 @ s12 @ r22 @ s23 þ þ ¼0 @x1 @x2 @x3 @ s31 @ s23 @ r33 þ þ ¼0 @x1 @x2 @x3
where rij (i, j = 1,2,3) is the normal stress, sij (i, j = 1, 2, 3) the shear stress. Correspondingly, the electrostatic equation without body charge is
@D1 @D2 @D3 þ þ ¼0 @x1 @x2 @x3
ð2Þ
where Di (i = 1, 2, 3) is the electric displacement component. The geometric conditions are related with the strains and displacements and they are expressed as follows
@u1 @x1 @u2 ¼ @x2 @u3 ¼ @x3
1 @u1 @u2 þ 2 @x2 @x1 1 @u1 @u3 ¼ þ 2 @x3 @x1 1 @u2 @u3 ¼ þ 2 @x3 @x2
e11 ¼
e12 ¼
e22
e13
e33
e23
ð3Þ
where eij (i, j = 1, 2, 3) is the strain, ui (i = 1, 2, 3) the displacement. The electric field can be written as
E1 ¼
@/ ; @x1
E2 ¼
@/ ; @x2
E3 ¼
@/ @x3
ð4Þ
where Ei (i = 1, 2, 3) is the electric field, / the electric potential. For transversely isotropic piezoelectric materials, the constitutive equations can be represented as [7]
ekl ¼ sklij rij þ dikl Ei ði; j ¼ 1; 2; 3; k; l ¼ 1; 2; 3Þ Dk ¼ dkij rij þ nki Ei
ð5Þ
where sklij is the elastic compliance constants, dkij the piezoelectric coefficients, nki the dielectric constants. The elastic constant matrix c is the inverse matrix of the full compliance matrix s [18]. For transverse isotropic materials, the elastic compliance matrix is [7,19]
2
Fig. 1. The schematic representation of film/substrate system in the selected coordinate system.
ð1Þ
1=ET
6 m =E T T 6 6 6 mTL =ET s¼6 6 0 6 6 4 0 0
3
mT =ET
mLT =EL
0
0
0
1=ET
mLT =EL
0
0
mTL =ET
1=EL
0
0
0
0
1=GT
0
0 0
0 0
0 0
1=GL 0
0 7 7 7 0 7 7 0 7 7 7 0 5 1=GL
ð6Þ
136
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
where ET is the transverse Young’s modulus, EL the longitudinal Young’s modulus, GT the transverse shear modulus, GL the longitudinal shear modulus, mT the Poisson’s ratio in the transversely isotropic (x1–x2) plane as shown in Fig. 1. mTL has the physical interpretation of the Poisson’s ratio that characterizes the strain in the x1–x2 plane resulting from stress normal to it, and mLT characterizes the longitudinal strain in the direction normal to the transversely isotropic plane resulting from stress in the x1–x2 plane. The Poisson’s ratios mTL and mLT are related as
mTL ET
¼
mLT
ð7Þ
EL
The transverse shear modulus GT can be expressed as
GT ¼
ET 2ð1 þ mT Þ
ð8Þ
The substrate is assumed as a linear elastic isotropic material characterized by Young’s modulus ES and the Poisson’s ratio mS [9].
Within the contact area of the conical indenter in Fig. 1, the downward displacement u2 is given as follow [14]
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 6 x21 þ x23 6 a ð9Þ
For the conical indenter, there is no friction between the indenter and the film, and it can be expressed as
s12 ¼ s23 ¼ 0; ðr P 0Þ
ð10Þ
There is no normal stress outside the contact region and it can be written as
r22 ¼ 0; ðr > aÞ
ð11Þ
The indenter is a perfect electrical insulator with a zero electric charge distribution. The electric charge distribution of the surface is zero and it can be given as
D2 ðx1 ; 0; x3 Þ ¼ 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 þ x23 > a
3. Evaluation of the engineering constants 3.1. Dimensionless analysis For an idealized transversely isotropic bulk material, the loading part of indentation curve from a conical indentation test obeys a quadratic relationship F = Ch2, where C is loading curvature, F indentation load, h indentation depth [20]. The quadratic relationship does not accurately describe the indentation response of thin film, so it is described by the power function F = Chx, in which the power exponent x is defined as the impact factor to reflect the substrate effect [7]. Then the relationship between indentation load F and indentation depth h can be rewritten as [7,21]
F ¼ F max
2.2. Boundary conditions
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 ðx1 ; 0; x3 Þ ¼ h x21 þ x23 cot a ;
and GL, and the simplified simulation results under PMI and PI modes can reflect the influence of the piezoelectric effect on characterization of engineering constants in the preliminary study.
ð12Þ
2.3. Simplification of the material model Generally, the mechanical and electrical properties of the transversely isotropic piezoelectric materials can be especially described by five independent engineering constants ET, EL, GL, mT and mTL, three piezoelectric coefficients e15, e31 and e33, and two dielectric constants n11 and n33[15]. As for the transversely isotropic materials without piezoelectric effect, the elastic properties can be only described by five independent engineering constants, and it is difficult to determine all the unknown engineering constants ET, EL, GL, mT and mTL [19]. A few of assumptions are conducted to simplify the material model, so the total number of independent mechanical variables should be reduced before determining the mechanical constants. When piezoelectric thin film is simplified as a kind of transversely isotropic materials, the Poisson’s ratio mT can be assumed as 0.3 [5]. The Poisson’s ratios mTL and mLT are related as mTL + mLT = 2mT [19]. When the piezoelectric and dielectric constants of the transversely isotropic thin film are taken as the known constants, the purely mechanical indentation (PMI) and piezoelectric indentation (PI) modes of ABAQUS software can be conveniently used to simulate the numerical loading curve in order to establish the dimensionless equations. Namely, the engineering constants of the transversely isotropic thin film can be characterized by the independent unknown engineering constants ET, EL
h hm
x ð13Þ
where Fmax and hm are the maximum indentation load and the corresponding maximum penetration. According to the previous results of transversely isotropic thin film indentation [22], the maximum indentation load Fmax and the loading curve exponent x should be the functions of the material constants, such as the engineering constants ET, EL and GL, the substrate elastic parameter ES, the geometric parameter of film thickness t, and the maximum indentation depth hm. Namely, they can be expressed as follows
F max ¼ /0 ðET ; EL ; GL ; ES ; t; hm Þ
ð14Þ
x ¼ u0 ðET ; EL ; GL ; ES ; t; hm Þ
ð15Þ
It is easy to know the parameters ES, t and hm by nanoindentation technique, therefore they are regarded as known constants in order to simulate the numerical loading curve, in which the numerical maximum indentation load and numerical loading curve exponent can be extracted. Applying Buckingham Pi Theorem to Eqs. (14) and (15), they can be written as dimensionless form
ET EL GL ¼ / ; ; 2 ES ES ES E S hm ET EL GL x¼u ; ; ES ES ES F max
ð16Þ ð17Þ
The elastic properties of transversely isotropic thin film can be described by the independent unknown engineering constants ET, EL and GL. In addition, the indentation modulus Em is the weighted average of the engineering constants of anisotropic materials along the different orientations and it is defined as [7,19]
Em ¼
EL þ ET 2
ð18Þ
Generally the unloading part of the indentation curve can reflect factual and comprehensive information of indented materials, such as plastic deformation, the expansion of surface and subsurface cracks during the indentation loading process, therefore Em is traditionally determined by the unloading part slope of indentation curve at the maximum indentation load according to the wellestablished Oliver and Pharr method. The elastic and plastic deformation processes are gained by the analysis of the indentation curve, and the widely studied parameters are the residual depth of penetration (hres), total work of contact deformation (Wtot), work from elastic recovery (Wel) and the work of inelastic deformation (Wpl). The numerical values of the quantities Wtot, Wel and Wpl can be directly obtained from the indentation curve, as the area under the loading curve, under the unloading curve and the area enclosed by the two curves, respectively. The residual depth ratio
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
hres/hm and plastic work ratio Wpl/Wtot characterize the extent of plastic deformation and thus are examined with respect to the deformation behavior of the materials. Although it is difficult to simulate the unloading part of the indentation curve by using commercial package, one can easily obtain the numerical loading curve under the simplified idealization through ABAQUS. As we all know the indentation curve can reflect ideal elasticity at low indentation load only, therefore we choose the simulation on loading curve in order to determine the elastic constants. However the surface and subsurface cracks may be expansive because of the maximum indentation loads at the beginning of the unloading process, and it indicates that it is very difficult to simulate the unloading part of the indentation. At the appropriate penetration depth the numerical maximum indentation loads and numerical loading curve exponents can be extracted from the numerical loading curves of the thin film to establish the dimensionless Eqs. (16) and (17), and the engineering constants EL, ET and GL can be determined by solving Eqs. (16)–(18). 3.2. Forward analysis 3.2.1. Finite element procedure Numerical indentation tests performed on the film/substrate nanoindentation system schematized as Fig. 1 are simulated by
137
using commercial numerical finite element package ABAQUS, and we adopt PMI and PI modes to obtain the numerical loading curves. For Berkovich indenter, the corresponding apex angle a of the equivalent cone was chosen as 70.3° [23]. The numerical indentation is modeled as a two-dimensional axisymmetric problem. A fine finite element mesh is used in the vicinity of the contact region to capture the localized deformation underneath the indenter and a gradually coarser mesh away from the contact region is used to cut down the computational work. As shown as Fig. 2, the substrate is meshed by using four-node reduced integrated linear axisymmetric solid element (CAX4R). The film is meshed by using CAX4R element for PMI mode and four-node reduced integrated linear axisymmetric piezoelectric element (CAX4E) for PI mode [16]. Eighty elements are used through the thin film thickness. Gradually, the influence of friction is usually weak on the nanoindentation test, and there is no friction at the interface between the indenter and the thin film [7,14]. The numerical indentation load is recorded as a function of the vertical displacement beneath the tip of the indenter. In order to consider the material constants accurately, a broad spectrum of possible material constants combinations 0.1 6 EL/ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 0.1 6 GL/ES 6 2.0 are chosen as input parmeters to simulate the indentation process, and then the numerical loading curves can be obtained. When the engineering constants of thin film ET, EL, GL and mT, the substrate
Fig. 2. Finite element mesh of an axisymmetric indentation.
138
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
elastic parameters ES and mS, and the maximum indentation depth hm are regarded as input parameters, the numerical loading curve in situ reflects the relationship between the indentation force on the rigid indenter and the indenter penetration into the thin film during FEM simulation. Generally, a well-known empirical rule is that the indentation depth should be smaller than one tenth of film thickness at micrometer to avoid substrate effect [21], however it is too shallow to extract the accurate materials properties of extreme thin film at nanometer. At first, it is very difficult to determine the contact area and depth accurately when the contact depth is of the same order as the surface roughness. Secondly, there is a little distinguishing with indentation force and adhesion between the indenter tip and specimen, in order that the latter may not be properly measured based on the present nanoindentation technique. Finally, even a brand new indenter tip is blunt at the nanoscale, with a radius of curvature typically tens of nanometers; the conventional sharp indentation analysis cannot be applied directly to the blunt tip at a small-scale. In a word, the indentation depth is generally larger than one tenth of the film thickness to obtain the exact properties [9,16]. On the other hand, the indentation depth must be limited for different materials to avoid excessive substrate deformation and cracking, and for simplicity it is chosen as one fifth of the film thickness as the previous report [7]. 3.2.2. Numerical loading curves According to the previous results on the materials constants, the broad spectrum of possible material combinations 0.1 6 EL/ ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 0.1 6 GL/ES 6 2.0 are chosen as input parameters to simulate the nanoindentation test, and von-Mises stress in contour and the numerical loading curves can be obtained. The former can indicate whether the indentation process is disturbed by the substrate, while the latter will give an insight into how the various engineering constants of thin film influence the shape of the numerical loading curves. For brevity, the different combinations of normalized engineering constants are fixed as EL/ ES = (0.6, 0.8, 1.0, 1.2), ET/ES = (0.5, 0.7, 0.8, 1.0) and GL/ES = (0.2, 0.3, 0.4, 0.5), and the indentation depth is chosen as one fifth of the film thickness as the previous report [7]. The typical von-Mises stress in contour at EL/ES = 1.0, ET/ES = 0.8 and GL/ES = 0.3 is given in Fig. 3.
Obviously, the stress decreases continuously away from the contact region in the thin film while it is discontinuousness at the interface of film/substrate, indicating that the indentation process is disturbed by the substrate. The typical numerical loading curves at ET/ES = 0.8 and GL/ ES = 0.3 can be obtained under PMI mode, and it is given in Fig. 4. 2 Both of the indentation load F=ES hm as a function of the indentation depth h/hm and the curvature of the numerical loading curve increase with the increase of EL/ES. At the fixed ratios of ET/ES and 2 GL/ES, the normalized maximum indentation load F max =ES hm increases with the ratio EL/ES. They give an insight into how the various engineering constants of thin film influence the shape of the numerical loading curves. According to the numerical loading curves in Fig. 4, the numer2 ical maximum indentation load F max =ES hm and the loading curve exponent x can be obtained after the numerical loading curves are fitted as the power function F = Fmax(h/hm)x with Origin soft2 ware. The F max =ES hm and x versus EL/ES and GL/ES are respectively plotted in Fig. 5a–f. With the increase of ET/ES, EL/ES and GL/ES, 2 F max =ES hm increases while x decreases. One can conclude that the 2 F max =ES hm and x regularly vary with the engineering constants, and the dimensionless functions given as Eqs. (16) and (17) will be established by fitting the engineering constants vs numerical
Fig. 4. The typical numerical loading curves obtained from FEM indentation test with ET/ES = 0.8 and GL/ES = 0.3 (a) for the PMI mode and (b) for the PI mode.
Fig. 3. The representative von–Mises stress in contour of numerical indentation.
139
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Fig. 5. Given the representative combinations of two normalized engineering constants, the curves of the other engineering constant related with the normalized maximum 2 load F max =ES hm (a)–(c) and the loading curve exponent x (d)–(f).
Table 1 Piezoelectric and dielectric constants of PZT-6B and PZT-4. Dielectric constants (1010 Fm1)
Materials
Piezoelectric constants (Cm2) e31
e33
e15
j11
j33
PZT-4 PZT-6B
13.44 4.6
6.98 0.9
13.84 7.1
60.0 36.0
54.7 34.0
2
maximum load F max =ES hm curves shown in Fig. 5(a)-(c) and the engineering constants vs numerical loading curve exponent x curves shown in Fig. 5d–e, which are used to solve the unknown engineering constants in the reverse analysis as the following.
3.2.3. Establish of dimensionless functions A plenty of combinations of the numerical maximum load 2 F max =ES hm and numerical loading curve exponent x extracted from the numerical loading curves will be used to determine the specific
form of dimensionless functions. According to Eqs. (16) and (17), letting g EL/ES, n ET/ES and c GL/ES, the dimensionless functions can be written as
/e ðn; g; cÞ ¼ p1 þ p2 n þ ðp3 þ p4 nÞc þ ðp5 þ p6 nÞg þ ðp7 þ p8 nÞcg þ ðp9 þ p10 nÞg2 þ ðp11 þ p12 nÞcg2
ð19Þ
ue ðn; g; cÞ ¼ q1 þ q2 n þ ðq3 þ q4 nÞg þ ðq5 þ q6 nÞc þ ðq7 þ q8 nÞgc þ ðq9 þ q10 nÞc2 þ ðq11 þ q12 nÞgc2
ð20Þ
140
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
The coefficients pi, qi (i = 1, 2, . . . , 12) fitted are listed in Table 2 for PMI mode and Table 3 for PI mode. 3.3. Reverse analysis In the reverse analysis, the nanoindentation tests were performed on PZT thin film laid on a silicon substrate. The experimental indentation curves, the maximum indentation load Pmax, the maximum indentation depth hm and the indentation modulus Em can be automatically obtained by the nanoindenter. According to the loading part of the indentation curve the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The reverse analysis algorithm is schematic as the flow chart in Fig. 6, and it includes the extraction of the experimental data, the setting on the combination range of unknown engineering constants, the calculation on the engineering constant errors, the selection combination of the smallest total error, and the output of the engineering constants. After determining the solutions of the engineering constants, the method validity is discussed. 3.3.1. Reverse analysis algorithm According to Fig. 6, at first ES is substrate elastic parameter, and Pmax, hm, xm and Em can be obtained from the nanoindentation test. When ES, Pmax, hm, xm and Em are known, the engineering constants ET, EL and GL can be determined by solving Eqs. (18)–(20). Secondly, the combination ranges of unknown engineering constants are assumed within the wide ranges, and the engineering constants EL/ES, ET/ES, GL/ES are setting as 0.1 6 EL/ES 6 5.0, 0.1 6 ET/ES 6 5.0 and 2 0.1 6 GL/ES 6 2.0. Thirdly, the errors of Pmax =ES hm , xm and Em can be calculated, and the total error is calculated as the sum of absolute values. Finally, the combination of engineering constants is selected as the ultimate solution with the smallest total error, and the values of ET, EL, GL are outputted. Table 2 The fitting coefficients in Eqs. (19) and (20) for PMI mode. Coefficients in Eq. (19)
Coefficients in Eq. (20)
p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12
0.21389 0.41697 0.79622 0.58033 1.44269 0.59941 12.2136 6.05127 4.06122 5.47017 1.80912 10.4618
2.25659 0.0372 0.04018 0.06076 1.19663 0.26827 1.75892 1.11014 0.95761 0.54791 0.34978 1.98759
Table 3 The fitting coefficients in Eqs. (19) and (20) for PI mode. Coefficients in Eq. (19)
p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
Coefficients in Eq. (20)
PZT-4
PZT-6B
0.57831 0.14498 0.13407 0.15477 1.36684 0.08392 1.12744 0.47344 0.41528 0.07819 0.36699 0.13658
0.21943 0.13188 0.25928 0.37999 2.59463 0.95794 3.85829 0.99306 0.72688 0.25946 2.66842 1.09466
q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12
PZT-4
PZT-6B
2.28813 0.05121 0.09849 0.04084 0.04187 0.21079 0.01569 0.08768 0.58711 0.07475 0.06713 0.11082
2.01627 0.04739 0.04061 0.29019 0.37009 0.09921 0.42868 0.04421 0.56598 0.18855 0.03988 0.08139
Fig. 6. The flow chart of the reverse analysis for evaluating engineering constants.
3.3.2. Application to PZT thin film Using the nanoindenter (Hysitron Triboindenter, Hysitron Inc., USA), the nanoindentation test was performed on PZT thin film deposited on silicon substrate, and the maximum indentation depth was set as one fifth of film thickness, i.e. 70 nm. five indents were performed during the nanoindentation test, and the five experimental curves ‘‘1–5’’ and one averaged experimental curve ‘‘a’’ are given in Fig. 7a, and AFM image of nanoindentation impression is described as Fig. 7b. The experiment data Pmax and xm are determined as 1.3743 mN and 2.0336 from the averaged experimental curve in Fig. 7a, and the indentation modulus Em is obtained as 139.4 GPa by using the well-established Oliver and Pharr [2] method. According to the loading part of the experimental indentation curves, the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The Young’s modulus of silicon substrate is assumed to be known as ES = 130 GPa for simplicity according to the experimental report [24]. The related experimental data is listed in Table 4. 3.3.3. Determination of the engineering constants 3.3.3.1. Solutions for PMI mode. Substituting Pmax, hm, ES and xm into Eqs. (18)–(20), the multiple solutions will be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the engineering constants are selected as the solutions given in Table 5. In order to verify the validity of the solutions, the solutions in Table 5 are taken as inputted parameters of ABAQUS software, and the corresponding numerical loading curves S1, S2, and S3 for the PMI mode are obtained and shown in Fig. 8. One can see that the numerical loading curves are in good agreement with the averaged experimental loading curve, and the numerical loading curve S3 with the smallest total error is closest to the averaged experimental loading
141
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Fig. 7. (a) The five experimental load versus depth curves and the averaged load versus depth curve and (b) AFM image of indentation impression on PZT thin film.
Table 4 The related experimental data used in reverse analysis.
Table 6 Some representative engineering constants of PZT thin film [25].
t (nm)
hm (nm)
Fmax (mN)
exponent x
d33 (pC/ N)
Em (GPa)
ES (GPa)
350
70
1.3743
2.0336
60
139.4
130
Table 5 The multiple solutions for the PMI mode. Solutions
etotal (%)
ET (GPa)
EL (GPa)
GL (GPa)
S1 S2 S3
2.385 1.361 1.049
157 154 152
124 125 128
65 64 61
Materials
ET (GPa)
EL (GPa)
GL (GPa)
PZT-4 PZT-2 PZT-6B
81.3 86.2 137.0
64.1 67.5 131.6
56.8 33.4 54.1
of PZT-6B and PZT-4 are chosen as the input parameters listed in Table 1. The multiple solutions can be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the engineering constants are selected as the solutions given in Table 7. In order to verify the validity of the solutions, the multiple solutions are taken as inputted parameters of ABAQUS software to obtain the numerical loading curves S1, S2, and S3 shown in Fig. 9. Obviously, the total error of S3 is smallest, and S3 is taken as the ultimate engineering constants for the PZT thin film under the PI mode of PZT-6B. Comparing with the engineering constants of PZT-6B [25], the relative differences of ET, EL and GL are 3.65%, 3.94% and 4.04%, indicating that the PZT sample is identified as PZT-6B. 4. Evaluation of the piezoelectric constants The variety of PZT thin film has been determined as PZT-6B by evaluating the engineering constants, and in this section, the piezoelectric constants of PZT-6B will be evaluated by combining FEM simulation with the nanoindentation test. 4.1. Dimensionless equations
Fig. 8. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 5 for the PMI mode.
curve. The solution corresponding to the closest numerical loading curve S3 is taken as the ultimate engineering constants to compare with the previous results, such as PZT-6B, PZT-2 and PZT-4 listed in Table 6 [25]. The relative differences of ET, EL and GL are 43.3%, 47.3%, 45.2% for PZT-2, 59.7%, 58.4%, 65.4% for PZT-4 and 9.9%, 2.8%, 11.3% for PZT-6B. Obviously, the relative differences are least for PZT-6B therefore the PZT thin film is primarily determined as PZT-6B. 3.3.3.2. Solutions for PI mode. In order to identify the variety of PZT thin film further, the piezoelectric effect is considered and the engineering constants can be solve for PI mode. According to the results under PMI mode, the piezoelectric and dielectric constants
The engineering constants and the dielectric constants of the PZT-6B are known parameters. According to the previous results of transversely isotropic piezoelectric thin film indentation [7], the maximum indentation load Fmax and the loading curve exponent x should be the functions of the material parameters, such as the piezoelectric constants e15, e31, and e33, the substrate elastic parameter ES, the geometric parameter of film thickness t, and the maximum indentation depth hm. Namely, they can be expressed as follows
F max ¼ /p ðe15 ; e31 ; e33 ; ES ; t; hm Þ
ð21Þ
x ¼ up ðe15 ; e31 ; e33 ; ES ; t; hm Þ
ð22Þ
The moderate ratio of penetration depth to the thin film thickness is controlled as one fifth, and then applying the Buckingham PI Theorem Eqs. (21) and (22) can be written as
142
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
Table 7 The multiple solutions for PI mode. Solutions
S1 S2 S3
PZT-6B
PZT-4
etotal (%)
ET (GPa)
EL (GPa)
GL (GPa)
etotal (%)
EL (GPa)
GL (GPa)
ET (GPa)
0.9 0.7 0.3
149 146 142
131 134 137
54 53 52
1.922 1.891 0.816
150 146 140
130 135 139
55 52 50
Fig. 9. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 7 for the PI mode.
e15 e31 ; 2 e33 e33 ES hm e15 e31 x ¼ up ; e33 e33 F max
¼ /p
ð23Þ ð24Þ
Based on the piezoelectric constitution, the longitudinal direction piezoelectric strain constant d33 can be written as [26]
d33 ¼ 2e31 s13 þ e33 s33
ð25Þ
where s13 and s33 are known according to Eq. (6). When d33 is the known parameter, the unknown piezoelectric constants e15, e31, and e33 will be determined by solving Eqs. (23)–(25). 4.2. Forward analysis 4.2.1. Numerical loading curves Numerical piezoelectric indentation responses are simulated using CAX4E elements. In order to cover a broad range of materials, the broad spectrum of possible material combinations 0 < e15/ e33 < 3.0 and 3.0 < e31/e33 < 0 are chosen as input parameters to obtain the numerical loading curves. For brevity, the different combinations of normalized piezoelectric constants are fixed as e15/ e33 = (0.1, 0.5, 1.0, 2.2) and e31/e33 = (0.01, 0.15, 1.4, 3.0), and they are used to obtain the numerical loading curves for evaluating piezoelectric constants. To evaluate the piezoelectric constants, the engineering constants and dielectric constants of PZT-6B thin film [25] are chosen as the input parameters, and the typical numerical loading curves are given in Fig. 10 under the combination e31/e33 = 1.4 and e15/e33 = (0.1, 0.5, 1.0, 2.2). In the figure, the indentation load increases with piezoelectric constants e15/e33 at the same indentation depth, which is attributed to the piezoelectric effect on the indentation responses. When the insulated indenter with zero distributed electric charge is pressed into the piezoelectric materials, the indenter creates a constant electric potential inside the contact area [15]. The electric field induced during indentation as a result of the electrical– mechanical coupling can resist or aid in the penetration of the indenter into the piezoelectric material depending on the electrical
Fig. 10. The typical numerical loading curves with different e15/e33 at e31/e33 = 1.4.
2
Fig. 11. The normalized maximum indentation load F max =ES hm as the function of the piezoelectric constant ratios e31/e33 at different ratios of e15/e33.
conductivity of the indenter and the surface boundary conditions of the indented material. After the numerical loading curves under the above four combinations are fitted as the power function F = Fmax(h/hm)x by Origin
143
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144 Table 8 The fitting coefficients in Eqs. (28) and (29). Coefficients in Eq. (28)
Coefficients in Eq. (29)
m1 m2 m3 m4 m5 m6 m7 m8 m9
n1 n2 n3 n4 n5 n6 n7 n8 n9
1.98558 0.12324 0.04294 0.02850 5.1170 0.00283 0.00178 0.00184 9.3870
2.06537 0.04794 8.8940 2.1750 0.01095 0.00397 6.8920 0.0020 7.7940
Fig. 12. The loading curve exponent x as the function of the piezoelectric constant ratios e31/e33 at different ratios of e15/e33.
software, the maximum indentation loads and the loading curve exponents can be obtained via extensive FEM simulations, the nor2 malized maximum indentation load F max =ES hm and loading curve exponent x as the function of the piezoelectric constant ratio e31/ e33 are given in Figs. 11 and 12 at different ratios of e15/e33. Obviously, with the increase of e31/e33, the normalized maximum 2 indentation load F max =ES hm decreases while the loading curve exponent x increases at the fixed e15/e33. The normalized maxi2 mum indentation load F max =ES hm and the loading curve exponent x continuously and regularly vary with e15/e33 and e31/e33, and the dimensionless functions given as Eqs. (23) and (24) will be established by fitting the piezoelectric constants vs numerical 2 maximum load F max =ES hm curve shown in Fig. 11 and the piezoelectric constants vs numerical loading curve exponent x curve shown in Fig. 12, which are used to solve the unknown piezoelectric constants in the reverse analysis as the following. 4.2.2. Dimensionless functions Similarly, according to Eqs. (23) and (24), letting a e15/e33 and b e31/e33 the dimensionless functions can be written as
F max 2
ES hm
¼ /p ða; bÞ
x ¼ up ða; bÞ
ð26Þ ð27Þ
The formulas of the dimensionless functions (26) and (27) are fitted by Origin software as follows
Fig. 13. The flow chart of the reverse analysis for evaluating the piezoelectric constants.
/p ða;bÞ ¼ m1 þm2 a þm3 a2 þðm4 þm5 a þm6 a2 Þbþðm7 þm8 a þm9 a2 Þb2 ð28Þ
up ða;bÞ ¼ n1 þn2 a þn3 a2 þðn4 þn5 a þn6 a2 Þbþðn7 þn8 a þn9 a2 Þb2
ð29Þ
The coefficients mi, ni (i = 1, 2, . . . , 9) in Eqs. (26) and (27) are listed in Table 8. 4.3. Reverse analysis In the reverse analysis, the nanoindentation tests were performed on PZT-6B thin film laid on a silicon substrate. The experimental indentation curves, the maximum indentation load Pmax and the maximum indentation depth hm can be automatically obtained by the nanoindenter. According to the loading part of the indentation curve the modified power function described as Eq. (13) can be fitted by using Origin software to determine the experimental loading curve exponent xm. The reverse analysis algorithm is schematic as the flow chart in Fig. 13, and it includes the extraction of the experimental data, the setting on the combination range of unknown piezoelectric constants, the calculation on the piezoelectric constants errors, the selection combination of the smallest
total error, and the output of the piezoelectric constants. After determining the solutions of the piezoelectric constants, the method validity is discussed. 4.3.1. Reverse analysis algorithm At first ES is substrate elastic parameter, and Pmax, hm and xm can be obtained from the nanoindentation test. The piezoelectric strain constant d33 can be determined by a PFM experiment. When Pmax, hm, xm and d33 are known, the piezoelectric constants e15, e31, and e33 can be determined by solving Eqs. (25)–(27). Secondly, the combination ranges of unknown engineering constants are assumed within the wide ranges, and the normalized piezoelectric constants e15/e33 and e31/e33 are setting as 0 < e15/e33 < 3.0 and 2 3.0 < e31/e33 < 0. Thirdly, the errors of P max =ES hm , xm and d33 can be calculated, and the total error is calculated as the sum of absolute values. Finally, the combination of engineering constants is selected as the ultimate solution with the smallest total error, and the values of e15, e31, e33 are outputted.
144
S.T. Song et al. / Computational Materials Science 63 (2012) 134–144
differences are 0.3–5.6% indicating that the piezoelectric constants are reasonable.
Table 9 Multiple solutions of piezoelectric constants. Solution
etotal (%)
e33 (C/m2)
e31 (C/m2)
e15 (C/m2)
S1 S2 S3
0.10 0.08 0.07
4.60 4.70 4.75
1.30 1.00 0.85
7.68 7.26 7.08
Fig. 14. Comparison of the experimental loading curve with the numerical loading curves obtained via solutions in Table 9.
4.3.2. Application to PZT thin film Using the nanoindenter (Hysitron Triboindenter, Hysitron Inc., USA), the nanoindentation test was performed on PZT thin film deposited on silicon substrate, and the maximum indentation depth was set as one fifth of film thickness, i.e. 70 nm. five indents were performed during the nanoindentation test, and the five experimental curves ‘‘1–5’’ and one averaged experimental curve ‘‘a’’ are given in Fig. 7a, and AFM image of nanoindentation impression is described as Fig. 7b. The experiment data Pmax and xm are determined as 1.3743 mN and 2.0336 from the averaged experimental curve in Fig. 7a, and the piezoelectric strain constant d33 is obtained as d33 = 60 pC/N according to the experimental report [27]. The Young’s modulus of silicon substrate is assumed to be known as ES = 130 GPa for simplicity according to the experimental report [24].
4.3.3. Determination of the piezoelectric constants Substituting Pmax, hm, ES, xm and d33 into Eqs. (25)–(27), the multiple solutions will be obtained within a certain error range given in the flow chart. When the total error is least, three optimization combinations of the piezoelectric constants are selected as the solutions given in Table 9. In order to verify the validity of the solutions, the solutions in Table 9 are taken as inputted parameters of ABAQUS software, and the corresponding numerical loading curves S1, S2, and S3 are obtained and shown in Fig. 14. One can see that the numerical loading curves are in good agreement with the averaged experimental loading curve, and the numerical loading curve S3 with the smallest total error is closest to the averaged experimental loading curve. The solution corresponding to the closest numerical loading curve S3 is taken as the ultimate piezoelectric constants of the PZT thin film. Similarly, comparing with the previous results e15 = 7.1, e31 = 0.9 and e33 = 4.6 [28], the relative
5. Conclusion Combining nanoindentation test with FEM simulation, a method is proposed to identify the variety of PZT thin film via evaluating engineering constants based on the weight average method, and then the piezoelectric constants e15, e31 and e33 of the PZT thin film are determined by combining forward and inverse analysis. The numerical loading curves are obtained from the extensive FEM simulations at the moderately deep indentation for the PMI and PI modes, and the nanoindentation tests were performed on PZT thin film to extract the maximum indentation load, the loading curve exponent and the maximum indentation depth. The engineering constants of PZT thin film are determined as ET = 142 GPa, EL = 137 GPa and GL = 52 GPa, and the variety of the PZT thin film is determined as PZT-6B. The piezoelectric constants of the PZT-6B are determined as e15 = 4.75 C/m2, e31 = 0.85 C/m2 and e33 = 7.08 C/m2 based on piezoelectric constitution by combining forward analysis and inverse analysis. The proposed method is effective to determine the engineering/piezoelectric constants of piezoelectric thin film. Acknowledgments This work was supported by PCSIRT (IRT1080), NNSF of China (10825209), Changjiang Scholar Incentive Program ([2009]17), Aid Program for Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan Province, and Shanghai Nano Special Foundation (11nm0502600). References [1] M.K. Khan, S.V. Hainsworth, M.E. Fitzpatrick, L. Edward, Comput. Mater. Sci. 49 (2010) 751–760. [2] W.C. Oliver, G.M. Pharr, J. Mater. Res. 7 (1992) 1564–1583. [3] S. Liu, Y.T. Gu, H. Huang, Mater. Sci. Eng. A 528 (2011) 7948–7951. [4] M.R. VanLandingham, J. Res. Natl. Inst. Stand Technol. 108 (2003) 249–265. [5] M. Bocciarelli, G. Bolzon, G. Maier, Mech. Mater. 37 (2005) 855–868. [6] R. Saha, W. Nix, Acta Mater. 50 (2002) 23–38. [7] J.S. Wang, X.J. Zheng, H. Zheng, et al., Comput. Mater. Sci. 49 (2010) 378–385. [8] F.Q. Yang, Mater. Sci. Eng. A 358 (2003) 226–232. [9] M.H. Zhao, X. Chen, Y. Xiang, J.J. Vlassk, et al., Acta Mater. 55 (2007) 6260– 6274. [10] D.J. Ma, K.W. Xu, J.W. He, et al., Surf. Coat. Technol. 116 (1999) 128–132. [11] Y.G. Liao, Y.C. Zhou, Y.L. Huang, et al., Mech. Mater. 41 (2009) 308–318. [12] K. Uchino, Piezoelectric Actuators and Ultrasonic Motors, Kluwer Academic Publishers, Boston, Dordrecht, London, 1997. [13] W.N. Sharpe, B. Yuan, R.L. Edwards, J. Microelectromech. Syst. 6 (1997) 193– 199. [14] J.H. Wang, C.Q. Chen, T.J. Lu, J. Mech. Phys. Solids 56 (2008) 3331–3351. [15] A.E. Giannakopoulos, S. Suresh, Acta Mater. 47 (1999) 2153–2164. [16] S.M. Han, Saha R, W.D. Nix, Acta Mater. 54 (2006) 1571–1581. [17] Z.K. Wang, B.L. Zheng, Int. J. Solids Struct. 32 (1995) 105–115. [18] S. Lee, B.C. Cho, H.C. Park, et al., Mater. Chem. Phys. 75 (2002) 174–177. [19] T. Nakamura, Y. Gu, Mech. Mater. 39 (2007) 340–356. [20] B. Poon, D. Rittel, G. Ravichandran, Int. J. Solids Struct. 45 (2008) 6399–6415. [21] H. Zheng, X.J. Zheng, J.S. Wang, et al., Mater. Des. 32 (2011) 1407–1413. [22] C. Gamonpilas, E.P. Busso, Mater. Sci. Eng. A 380 (2004) 52–61. [23] M. Dao, N. Chollacoop, K.J. Van Vlient, et al., Acta Mater. 49 (2001) 3899–3918. [24] J.J. Wortman, R.A. Evans, J. Appl. Phys. 36 (1965) 153–156. [25] M.H. Zhao, Y.P. Shen, Y.J. Liu, et al., Theor. Appl. Fract. Mech. 26 (1997) 141– 149. [26] D. Damjanovic, Rep. Prog. Phys. 61 (1998) 1267–1324. [27] S. Hiboux, P. Muralt, N. Setter, Mater. Res. Soc. 596 (2000) 499–504. [28] M.H. Zhao, N. Li, C.Y. Fan, Eng. Anal. Bound Elem. 33 (2008) 545–555.