2251
A mathematical model in three-dimensional
piezoelectric continuum to predict non-linear
responses of piezoceramic materials
M K Samal1∗, P Seshu2, U von Wagner4, P Hagedorn3, B K Dutta1,and H S Kushwaha1
1Reactor Safety Division, Bhabha Atomic Research Centre, Trombay, Mumbai, India
2Department of Mechanical Engineering, Indian Institute of Technology Bombay, Mumbai, India
3Department of Applied Mechanics, Technical University of Darmstadt, Darmstadt, Germany
4Department of Applied Mechanics, Technical University of Berlin, Berlin, Germany
The manuscript was received on 10 January 2008 and was accepted after revision for publication on 28 April 2008.
DOI: 10.1243/09544062JMES1002
Abstract: It has been experimentally observed that the piezoceramic materials exhibit different
types of non-linearities under different combinations of electrical and mechanical fields. When
excited near resonance in the presence of weak electric fields, they exhibit typical non-linearities
similar to a Duffing oscillator such as jump phenomena and the presence of superharmonics in
the response spectra. In this work, these non-linearities have been modelled for a generalized
three-dimensional piezoelectric continuum using higher-order quadratic and cubic terms in the
electric enthalpy density function and the virtual work. The identification of the parameters
of the model requires a closed form solution for non-linear response of a simplified geometry.
A simple proportional damping formulation has been used in the model. Experiments have
been conducted on rectangular and cylindrical geometries of piezoceramic PIC 181 at different
magnitudes of applied electric fields and results have been compared with those of simulation.
Keywords: analytical solution, piezoceramic, non-linear modelling, perturbation method
1 INTRODUCTION
Day by day, the usage of piezoceramic materials is
increasing and these materials are widely applied
in many different areas of modern technology, for
instance, aerospace structures [1]. Prediction of per-
formance of the structures with piezoelectric sensors
and actuators (for various control applications such as
active shape, vibration, noise control, etc.) requires a
thorough understanding of the piezoelectric material
behaviour at different applied electrical and mechan-
ical fields. Though the non-linearities under strong
electric fields such as dielectric and butterfly hystere-
sis have drawn the attention of many authors recently
[2–7], the modelling of non-linear behaviour of piezo-
ceramics under weak electric fields is limited in liter-
ature. With a number of piezoelectric finite-element
∗Corresponding author: Reactor Safety Division, Bhabha Atomic
Research Centre, Hall 7, Trombay, Mumbai 400085, Maharashtra,
(FE) models existing in the literature to cope with
complicated boundary problems of smart structures,
it is very important to develop a simple procedure to
obtain analytical solutions for a few simple geome-
tries for verification of finite-element formulations
and also to gain confidence in the procedure. More-
over, it is also necessary to develop analytical solutions
to determine various material parameters including
those of the non-linear model to be used in the finite-
element analysis. The coupling of mechanical and
electrical properties of piezoelectric materials com-
plicates the analytical solution. Piezoelectric plates
comprising orthotropic or anisotropic layers are usu-
ally analysed using laminated composite plate theory.
For the analytical solution of these plates, Pagano
[8] first presented the exact solution of a laminated
plate under cylindrical bending. Ray et al. [9–11]
and Heyliger and Brooks [12] extended this method-
ology to develop exact solutions for the piezoelec-
tric laminated plates using Eshelby–Stroh formalism.
Later, considering thermal effect, Xu et al. [13]pre-
sented a three-dimensional analysis for the coupled
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2252 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
thermo-electro-elastic response of multi-layered
hybrid composite plates with four simply supported
edges by using the state space method. Dube et al. [14]
developed a series solution method for a single-layer
thermo-electro-elastic plate subjected to cylindrical
bending using the approach of Pagano [8].
To solve the dynamic problem of the piezoelec-
tric array elements, several simplified solutions were
tried for modelling the linear piezoelectric behaviour
in one-, two-, and three-dimensional domains.
Ray et al. [15] presented an exact solution for analysing
the dynamics of composite plates with piezoelec-
tric layers bonded to the top and bottom surfaces.
Later, Ray [16] also presented the closed form solu-
tion for the optimal control of flexural vibration of
a simply supported symmetric laminated plate. Wolf
and Gottlieb [17] derived the non-linear equations of
motion for a cantilever beam covered by piezoelec-
tric ceramic material (PZT) layers in symmetric and
asymmetric configurations. By defining a non-linear
enthalpy function with non-linear strain and applying
the Hamilton’s principle, they derived the equations
of motion and solved by perturbation analysis and
multiple scales method.
Beige and Schmidt [18] modelled the non-
linearities under weak electric fields using higher-
order quadratic and cubic elastic terms in the electric
enthalpy function. von Wagner and Hagedorn [19]
studied the non-linearity in a piezobeam system.
Later von Wagner [20] and Neumann [21] predicted
the amplitude of displacement and current responses
of piezorods of different diameters using a one-
dimensional non-linear formulation. Neumann [21]
used the non-linear damping formulation to derive
the virtual work done by the dissipative damping
forces. Samal et al. [22] presented a generalized three-
dimensional formulation for modelling the non-linear
behaviour of piezoelectric materials under weak elec-
tric fieldsandit was used toderive athree-dimensional
finite-element model. The application of the formu-
lation was demonstrated in Samal et al. [23] using
different types of examples. The parameter identifica-
tion procedure used by vonWagner [20] and Neumann
[21] is restricted to those of one-dimensional contin-
uum. In order to predict the non-linear response of
actual three-dimensional piezoelectric systems, the
material properties of the three-dimensional piezo-
electric continuum are required.This requires a closed
form solution for a simple geometry. In this paper,
such closed form solutions for predicting the displace-
ment and current responses of simple rectangular
and cylindrical geometries have been developed using
a simple proportional damping formulation for the
material damping.
In an earlier work, Samal et al. [24] have derived a
closed form solution for a rectangular piezoceramic
plate of a soft piezoceramic PIC 255. However, a
non-linear damping formulation was used in the
above model that requires identification of a large
number of non-linear model parameters, which is very
time-consuming, complicated, and cumbersome. On
the other hand, Usher and Sim [25], von Wagner [20],
Wang et al. [26], etc. demonstrated that the simple pro-
portional damping theory satisfactorily predicts the
response of wide varieties of piezoceramics. Hence,
the objective of this work is to develop an efficient
non-linear model using a simple damping theory, so
that the electric field dependent damping coefficient
can be evaluated easily from experimental measure-
ments and at the same time both the mechanical
and electrical responses of the piezoceramic are sat-
isfactorily predicted. Another important aspect of this
model is the use of an approach where the numerical
shape functions (obtained from a prior FE analysis)
can be used in place of analytical shape functions
in the Rayleigh–Ritz discretization. This approach is
especially helpful for geometries where the evaluation
of analytical shape functions is difficult. Experiments
have been conducted on plate and cylindrical geome-
tries of a PIC 181 material (which is a hard piezoce-
ramic) and the predicted results have been compared
with those of the experiments.
2 GENERALIZED NON-LINEAR CONSTITUTIVE
EQUATION OF THE PIEZOELECTRIC
CONTINUUM
The electric enthalpy density function His generally
used to derive the governing equations of the coupled
piezoelectric continuum and for the linear piezoelec-
tric behaviour, it is given as (IEEE standard [27])
H=1
2CE
ijklSijSkl −ekijEkSij −1
2εS
ijEiEj(1)
where CE
ijkl is the fourth-order elastic tensor under con-
stant electric field, ekij is the third-order piezoelectric
tensor, εS
ij is the second-order dielectric tensor, Sij is
the strain tensor, and Eiis the electric field vector.
The stress tensor Tij and the electric displacement vec-
tor Dican be derived from Husing the expressions
Tij =∂H/∂Sij and Di=−∂H/∂Ei, respectively. It can
be easily seen that the expressions for Tij and Diare
coupled with the secondary field variables Sij (strain
tensor) and Ei(electric field vector), respectively. The
strain tensor, Sij is expressed in terms of mechanical
displacement vector uias Sij =(ui,j+uj,i)/2 and Eiis
expressed in terms of φas Ei=−φ,i. Keeping in mind
the fact that stress and strain tensor are symmetric, H
can be simplified and written in terms of second order
tensors as [28]
H=1
2CE
ij SiSj−eijEiSj−1
2εS
ijEiEj(2)
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2253
In order to model the non-linearities in a coupled
piezoelectric medium, the linear electric enthalpy
density function is modified in the following way. It
was experimentally observed [19,20] that second- and
third-order harmonics were present in the response
spectra of PZT piezoceramic plates and rod geome-
tries. Hence, the generalized expression for non-linear
electric enthalpy density function can be written
as [22]
Hnon-lin =1
2CE
ij SiSj−emiEmSi−1
2εS
mnEmEn
+1
3CE
ijkSiSjSk+1
4CE
ijklSiSjSkSl
−1
3εS
mnoEmEnEo−1
4εS
mnopEmEnEoEp
−1
2emijEmSiSj−1
2e∗
mniEmEnSi
−1
3e∗∗
mijkEmSiSjSk−1
4e∗∗∗
mnijEmEnSiSj
−1
3e∗∗∗
mnoiEmEnEoSi(3)
where the coefficients with 3 or 4 indices are higher-
order cubic and quadratic coefficients, respectively. A
matrix expression for Hnon-lin can be obtained from
equation (3) as follows
Hnon-lin =1
2{S}T[C]{S}−{E}T[d][C]{S}
−1
2{E}T[εT]−[d][C][d]T{E}
+1
3{S}T[C1]{S2}+1
4{S2}T[C21]{S2}
+1
4{S}T[C22]{S3}−1
2{E}T[γ11]{S2}
−1
2{E2}T[γ12]{S}−1
3{E}T[υ1]{E2}
−1
3{E}T[γ21]{S3}−1
2{E2}T[γ22]{S2}
−1
3{E3}T[γ23]{S}−1
4{E}T[υ21]{E3}
−1
4{E2}T[υ22]{E2}(4)
where {S}is the strain vector, [C]is the linear elas-
ticity matrix, [d]is the linear piezoelectric coefficient
matrix, {E}is the electric field vector, [εT]is the dielec-
tric coefficient matrix, [C1]is the quadratic elasticity
matrix, and [C21]and [C22]are the cubic elasticity
matrices at constant electric field. The superscript
E and S (for the matrices at constant electric field
and constant strain) are omitted for clarity. The terms
[γ11]and [γ12]are quadratic piezoelectric matrices,
[γ21],[γ22], and [γ23]are cubic piezoelectric matrices,
[υ0]=[[εT]−[d][C][d]T],[υ1],[υ21], and[υ22]are linear,
quadratic, and cubic dielectric matrices at constant
strain field, S, respectively. The vectors with super-
scripts 2 and 3 are defined as the vectors with square
and cube of individual terms of the vectors. The first
three terms of Hnon-lin correspond to the linear elec-
tric enthalpy density function. Using this expression
of Hnon-lin and the virtual work done by damping
forces (described later in this paper), the closed form
solutions for simple rectangular plate and cylindrical
geometries have been derived in the next section using
the Hamilton’s principle.
3 HAMILTON’S PRINCIPLE
The dynamic equations of a piezoelectric contin-
uum can be derived from the Hamilton’s principle, in
which the Lagrangian and the virtual work is properly
adapted to include the electrical, mechanical as well
as the coupled electro-mechanical terms. The poten-
tial energy density of a piezoelectric material includes
contributions from the strain energy and from the
electrostatic energy. The Hamilton’s principle can be
written as
δt1
t0V
LdVdt+t1
t0
δWdt=0 (5)
where t0and t1define the time interval, Lis the
Lagrangian and defined in terms of the kinetic energy
density Tke, electrical enthalpy density Has L=
(Tke −H), and δWis the virtual work done by external
mechanical and electrical forces. The kinetic energy is
given as
Tke =1
2ρ{˙
u}T{˙
u}(6)
where ρis the mass density, {u}is the general-
ized displacement field and {˙
u}is the generalized
velocity field.
The non-linear electric enthalpy density Hnon-lin is
given in equation (4). The virtual work done by the
external mechanical forces Fand the applied electric
charges Qfor an arbitrary variation of the displace-
ment field {δu}and electrical potential {δφ}both satis-
fying the essential boundary conditions (i.e. {δu}=0
on surface As3and {δφ}=0 on surface As4) can be
written as
δW=V{δu}T{FV}dV+As1 {δu}T{FAs}dAs
+{δu}T{FP}−As2
δφς dAs−δφQ+δWD(7)
where {FV}is the applied body force vector, {FAs}is
the applied surface force vector (defined on the sur-
face As1),{FP}is the applied point load vector, φis
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2254 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
the electric potential, ςis the surface charge brought
on surface As2,Qis the applied concentrated elec-
tric charge, and δWDis the work done by damping
forces. Figure 1 shows the various types of mechanical
forces and electrical charges operating in the domain
V. Substituting the non-linear electric enthalpy den-
sity function Hnon-lin and kinetic energy density Tke
from equations (4) and (6) in equation (5), yields the
following expression for the Hamilton’s principle
−t1
t0Vρ{δu}T{¨
u}−t1
t0V{δS}T[C]{S}
−[C]T[d]T{E}−[diag({S})][γ11]T{E}−1
2[γ12]T{E2}
+diag({S})[C21]diag({S}){S}+1
4[C22]{S3}
+3
4diag({S2})[C22]T{S}+1
3[C1]{S2}
+2
3[diag({S})][C1]T{S}−[diag({S2})][γ21]T{E}
−[diag({S})][γ22]T{E2}−1
3[γ23]T{E}3dVdt
−t1
t0V{δE}Td[C]{S}+[υ0]{E}+1
2[γ11]{S2}
+[diag({E})][γ12]{S}+1
3[γ13]{S2}
+2
3[diag({E})][γ13]T{E}+1
3[γ21]{S3}
+[diag({E})][γ22]{S2}+[diag({E2})][γ23]{S}
+1
4[υ1]{E3}+3
4[diag({E2})][υ1]T{E}
+[diag({E})][υ2][diag({E})]{E}dVdt
+t1
t0V{δu}T{Fv}dVdt+t1
t0As1{δu}T{FAs}dAsdt
+t1
t0{δu}T{FP}dt−t1
t0As2
δφζ dAsdt
−t1
t0
δφQdt−t1
t0[δWD]dt=0 (8)
Fig. 1 Generalized mechanical and electrical forces of
various types acting in the domain V
where the matrix with ‘diag’ refers to a diagonal matrix
with the main diagonal terms being the terms of its
vector. For deriving the expression for the virtual work
done by damping forces, the proportional damping
or Rayleigh damping model is considered and the
formulation is given in next section.
4 PROPORTIONAL DAMPING FORMULATION
In this model, modal proportional damping is con-
sidered for the piezoceramic material and hence
the virtual work done by viscous damping forces is
given as (using the principle of Rayleigh proportional
damping)
δWD=V{δu}T[Cdamp]{˙
u}dV(9)
where [Cdamp]is the proportional damping coefficient
matrix and is expressed in terms of mass [M]and stiff-
ness [K]matrices with the help of constants αand βas
[Cdamp]=α[M]+β[K](10)
5 THREE-DIMENSIONAL ANALYSIS OF A
RECTANGULAR PIEZOCERAMIC PLATE
The rectangular plate analysed in this work is shown
in Fig. 2. One needs to perform three-dimensional
analysis in order to determine all the elements of the
non-linear material property matrices so that it can
be used in the finite-element analysis to analyse other
complex geometries made up of the same material.
The form of solution for the non-linear displacement
amplitude using Rayleigh–Ritz method and perturba-
tion technique is described in the next section where
three-dimensional matrices and vectors are used in
the electric enthalpy function, kinetic energy, and
virtual work.
Fig. 2 Rectangular piezoceramic plate excited with
electric field applied across its thickness
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2255
Using the Hamilton’s principle, one gets three partial
differential equations for motions in three co-ordinate
directions (i.e. x,y, and z). The displacement degrees
of freedom used in the analysis are u,v, and walong
x,y, and zdirections. The governing equations are
solved along with the boundary conditions to get
the shape functions of the problem so that it can be
used in Rayleigh–Ritz method to discretize the equa-
tions of motion and to solve by perturbation analysis.
In this section, the solution of the problem will be
discussed for deriving the shape functions for the
three-dimensional problem corresponding to the first
in-plane mode of free-free vibration of the rectangu-
lar plate. Again for simplification, it will be assumed
that the electric field is constant in z-direction and the
expression is given as follows
Ez=−U0
hcos(t)(11)
In this three-dimensional analysis, all the stress and
strain components exist in the formulation. The
kinematic relations between strain and the displace-
ment for the three-dimensional continuum are as
follows
Sxx =S1=∂u
∂x,Syy =S2=∂v
∂y,Szz =S3=∂w
∂z
Syz =S4=∂v
∂z+∂w
∂y,Sxz =S5=∂u
∂z+∂w
∂x
Sxy =S6=∂u
∂y+∂v
∂x
(12)
The three-dimensional linear material property matrix
[c]relating strains to stresses for three-dimensional
continuum are given as
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
T1
T2
T3
T4
T5
T6
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
c11 c12 c13 000
c21 c22 c23 000
c31 c32 c33 000
000c44 00
0000c55 0
00000c66
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
S1
S2
S3
S4
S5
S6
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
(13)
5.1 Derivation of the shape function of the
three-dimensional linear eigenvalue problem
The linear electric enthalpy function Hlinear corre-
sponding to equation (2) has been used to derive
the shape functions. After substituting the electric
enthalpy density function in the Hamilton’s principle,
i.e. equation (5), and performing variation with respect
to u,v, and w, the three equations of motion yielded,
respectively, are
ρ¨
u=c11
∂2u
∂x2+c66
∂2u
∂y2+c55
∂2u
∂z2+(c12 +c66)∂2v
∂x∂y
+(c13 +c55)∂2w
∂x∂z
ρ¨
v=(c21 +c66)∂2u
∂x∂y+c66
∂2v
∂x2+c22
∂2v
∂y2+c44
∂2v
∂z2
+(c23 +c44)∂2w
∂y∂z
ρ¨
w=(c31 +c55)∂2u
∂x∂z+(c32 +c44)∂2v
∂y∂z+c55
∂2w
∂x2
+c44
∂2w
∂y2+c33
∂2w
∂z2
(14)
and the associated boundary conditions are given
in Appendix 1. Using variable separable method, the
solutions for u,v, and ware as follows
u(x,y,z,t)=U(x,y,z)p(t);v(x,y,z,t)
=V(x,y,z)p(t);w(x,y,z,t)
=W(x,y,z)p(t)(15)
Substituting u(x,y,z,t),v(x,y,z,t), and w(x,y,z,t)in
equation(14)and separatingthe termsspace andtime,
yields
¨
p(t)
p(t)=constant =−ω2
c11
∂2U(x,y,z)
∂x2+c66
∂2U(x,y,z)
∂y2+c55
∂2U(x,y,z)
∂z2
+(c12 +c66)∂2V(x,y,z)
∂x∂y+(c13 +c55)∂2W(x,y,z)
∂x∂z
=−ρω2U(x,y,z)
(c21 +c66)∂2U(x,y,z)
∂x∂y+c66
∂2V(x,y,z)
∂x2
+c22
∂2V(x,y,z)
∂y2+c44
∂2V(x,y,z)
∂z2
+(c23 +c44)∂2W(x,y,z)
∂y∂z=−ρω2V(x,y,z)
(c31 +c55)∂2U(x,y,z)
∂x∂z+(c32 +c44)∂2V(x,y,z)
∂y∂z
+c55
∂2W(x,y,z)
∂x2+c44
∂2W(x,y,z)
∂y2
+c33
∂2W(x,y,z)
∂z2=−ρω2W(x,y,z)
(16)
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2256 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
The solutions can be assumed as
U(x,y,z)=u0sin(λ1x)cos(λ2y)cos(λ3y)
V(x,y,z)=v0cos(λ1x)sin(λ2y)cos(λ3y)
W(x,y,z)=w0cos(λ1x)cos(λ2y)sin(λ3y)
p(t)=sin(ωt)or cos(ωt)
(depending upon the initial condition)
(17)
Substituting the assumed solutions into the equation
(14), yields the eigenvalue equation as
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
c11λ2
1+c66λ2
2(c12 +c66)λ1λ2(c13 +c55)λ1λ3
+c55λ2
3
(c21 +c66)λ1λ2c66λ2
1+c22λ2
2(c23 +c44)λ2λ3
+c44λ2
3
(c31 +c55)λ1λ3(c32 +c44)λ2λ3c55λ2
1+c44λ2
2
+c33λ2
3
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
×⎧
⎨
⎩
u0
v0
w0
⎫
⎬
⎭=ρω2⎧
⎨
⎩
u0
v0
w0
⎫
⎬
⎭
(18)
Solving the eigenvalue equations (i.e. putting the
determinant of the matrix to be zero), yields
the relationship between the wave numbers (i.e.
λ1,λ2, and λ3)and the eigenfrequencies ω. The
assumed solutions must satisfy all the boundary con-
ditions given in Appendix 1 (for short-circuited elec-
trodes, i.e. Ez=0). Substituting the expressions for u,
w, and vin the boundary conditions yields various
combinations of the five unknowns (i.e. v0/u0,w0/u0,
λ1,λ2, and λ3). Hence, finally u,v, and wbecome
summation of various combinations of the terms con-
taining (λ1)i,(λ2)j, and (λ3)k, where i,j, and krefer to
each possible solution. Hence u(x,y,z,t),v(x,y,z,t),
and w(x,y,z,t)can be written as
u(x,y,z,t)=
i=1,n
j=1,n
k=1,n
u0ijk sin(λ1)i
×xcos(λ2)jycos(λ3)kzpijk(t)
v(x,y,z,t)=
i=1,n
j=1,n
k=1,n
v0ijk cos(λ1)i
×xsin(λ2)jycos(λ3)kzpijk (t)
w(x,y,z,t)=
i=1,n
j=1,n
k=1,n
w0ijk cos(λ1)i
×xcos(λ2)jysin(λ3)kzpijk (t)
(19)
For the first longitudinal mode of the plate (when
length (l)>width (b)>height (h)), the solutions for
Fig. 3 Mode-shape of the plate in the first free-free
in-plane longitudinal vibration mode
U(x,y,z),V(x,y,z), and W(x,y,z)are given as
λ1=λ2=λ3=π
l;v0
u0=−c21c33 −c31c23
c22c33 −c32c23
w0
u0=−c21c32 −c31c22
c22c32 −c22c33
U(x,y,z)=u0sinπx
lcosπy
lcosπz
l
V(x,y)=v0cosπx
lsinπy
lcosπz
l
W(x,y)=w0cosπx
lcosπy
lsinπz
l(20)
Using these shape functions, the non-linear equations
ofmotion for the maximumamplitude Ausing thepro-
cedure can be solved as described in the next section.
The mode shape of the plate for the first free-free in-
plane mode is shown in Fig. 3 and the three shape
functions U(x,y,z),V(x,y,z), and W(x,y,z)for the
same mode are shown in Figs 4(a) to (c), respectively.
5.2 Solution for the non-linear analysis with
proportional damping formulation
The Rayleigh–Ritz method has been used to solve the
non-linear equations represented by the variational
equation (8) and the displacement functions used
in the analysis represented by equation (15). These
functions are substituted into the variational equation
form of the Hamilton’s principle containing non-linear
enthalpy function Hnon-lin fromequation(4).Theshape
functions corresponding to the mode knearest to the
excitation frequency alone are used. For the first lon-
gitudinal vibration mode of the plate, the mode shape
functions as given by equation (20) are used in this
analysis. After performing the variation of the expres-
sion in the Hamilton’s principle with respect to pk(t),
yields the discretized non-linear equation of motion as
mk¨
pk(t)+ckpk(t)+Bkp2
k(t)+εkp3
k(t)
=f(1)
kcos t+f(2)
kpk(t)cos t+f(3)
kp2
k(t)cos t
+f(4)
kcos2t+f(5)
kpk(t)cos2t+f(6)
kcos3t
(21)
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2257
Fig. 4 Contour plot of three-dimensional shape func-
tions for the first in-plane free-free longitudi-
nal vibration mode of the plate on the origi-
nal plate geometry: (a) U(x,y,z); (b) V(x,y,z);
(c) W(x,y,z)
where k=1, 2, 3, ... and the coefficients are given in
Appendix 2.
5.3 Solution by perturbation analysis for the
formulation with proportional damping
An approximate solution for the non-linear equation,
i.e. equation (21), can be found using perturbation
technique [29]. Linear modal proportional damp-
ing with the non-dimensional damping coefficient
ξhas been used in this analysis. One can non-
dimensionalize time as
τ=ω0t(22)
where ω0=√ck/mkis the natural frequency of the
kth mode of the system. Equation (21) can now
be written as
p +εD∗p+p+εχp2+εα∗p3
=ε(h1cos ητ +h2sin ητ)(1+ψ0p+β∗p2)
+ε2γ(h1cos ητ +h2sin ητ)2(1+ψ1p)
+ε3δ∗(h1cos ητ +h2sin ητ)3(23)
where
()
=∂
∂τ ,η=
ω0
,D∗=2ξ
ε,χ=Bk
mkω2
0ε
h2
1+h2
2=h2
0,h0=f(1)
k
mkω2
0ε,ψ0=f(2)
k
f(1)
k
β∗=f(3)
k
f(1)
k
,γ=f(4)
k
mkω2
0ε2h2
0
,ψ1=f(5)
k
f(4)
k
δ∗=f(6)
k
mkω2
0ε3h3
0
(24)
The coefficients h1and h2are introduced in order to
allow for a phase shift between the excitation and
the response of the system due to the presence of
damping. The parameter εis called the perturbation
parameter and can be chosen arbitrarily. But, it should
be very small, e.g. 0.001. Now, the Lindstedt–Poincaré
method is used with the following expansion for p
and η.
p=p0+εp1+···,η=η0+εη1+··· (25)
Substituting equation (25) in equation (23) and col-
lecting the coefficients for ε0and ε1, yields
p
0+η2p0=0
p
1+η2p1=−D∗p
0+2η1p0−χp2
0−α∗p3
0
+(h1cos ητ +h2sin ητ)
×(1+ψ0p0+β∗p2
0)
(26)
The above expressions were obtained by using the
following expansion
1=η−εη1(taking only the first two terms
of the expansion of η)
⇒12=(η −εη1)2(27)
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2258 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
The solution of first part of equation (26) can be
written as
p0=Acos ητ (28)
where the amplitude Ais unknown. By substituting
equation (28) for p0in second part of equation (26),
one gets
p
1+η2p1=ξ(ηAsin ητ) +2η1Acos ητ
−χA2cos2ητ −α∗A3cos3ητ
+h1cos ητ +h1ψ0Acos2ητ
+h1β∗A2cos3ητ +h2sin ητ
+h2ψ0Asin ητ cos ητ
+h2β∗A2sin ητ cos2ητ (29)
Rewriting equation (29) in terms of sin ητ, cos ητ,
sin 2ητ, cos 2ητ, sin 3ητ, and cos 3ητ, one gets
p
1+η2p1
=ξηA+h2+1
4h2β∗A2sin ητ
+2η1A−3
4α∗A3+h1+3
4h1β∗A2cos ητ
+−1
2ξA2+1
2h1ψ0A+1
2h2ψ0Asin 2ητ
+1
2(h1ψ0A−ξA2)cos 2ητ +1
4h2β∗A2sin 3ητ
+1
4(h1β∗A2−α∗A3)cos 3ητ (30)
For the solution of equation (30) to remain finite, the
secular terms (terms with sin ητ and cos ητ) must be
zero [29]. Hence, one gets
h1=−2η1A+(3/4)α∗A3
1+(3/4)β∗A2,h2=−D∗Aη
1+(1/4)β∗A2
(31)
Combining h1and h2, one gets h0which is known and
this relation reduces to a polynomial equation in the
unknown stationary amplitude Awhich is expressed as
a5A10 +a4A8+a3A6+a2A4+a1A2+a0=0 (32)
where the coefficients are written as
a5=9
256α∗2β∗2
a4=9
32α∗2β∗−9
256β∗4h2
0−3
16β∗2η1α∗
a3=9
16α∗2−3
8β∗3h2
0
+9
16β∗2D∗2η2−3
2α∗β∗η1+1
4β∗2η2
1
a2=−11
8β∗2h2
0+3
2β∗D∗2η2−3α∗η1+2β∗η2
1
a1=−2β∗h2
0+D∗2η2+4η2
1,a0=−h2
(33)
The roots of the polynomial represent the solutions
for the stationary amplitude A. In the region of the
response curve where jump phenomenon is observed,
it will give three real roots and two imaginary roots.
The three real roots correspond to sweep-up frequency
response, sweep-down frequency response, and the
unstable responses, respectively. At other regions, the
polynomial will have only one real root and four imag-
inary roots and hence there is a unique response for
both sweep-up and sweep-down analyses. The total
response corresponding to the series expansion of p
is given by (neglecting higher-order terms like sin 2ητ,
cos 2ητ, sin 3ητ, and cos 3ητ)
p=p0+εp1=Acos ητ +ε
2η2(h1ψ0A−χA2)(34)
After the amplitude of displacement is calculated,
current flowing through the piezoceramic can be cal-
culated following the method for current calculation
as described in next section.
5.4 Calculation of current flowing through the
piezoceramic plate
The current flowing through the piezoelectric plate
can be calculated using the expression for the electric
displacement as follows
I(t)=As˙
Dz(t)dAs(35)
where Asis the area through which current is flowing in
the piezoceramic plate and zis the direction perpen-
dicular to the area Asand ˙
Dzis the time derivative of
the z-component of the electric displacement vector
{D}.The electric displacement vector {D}can be calcu-
lated from the generalized non-linear electric enthalpy
density function as
{D}=−∂Hnon-lin
∂{E}(36)
Using equation (4) for Hnon-lin in the equation (36) and
taking time derivative, one can write {˙
D}as
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2259
{˙
D}=[d][c]{˙
S}+[υ0]{˙
E}+[γ11]{S˙
S}
+[diag({˙
E})][γ12]{S}+[diag({E})][γ12]{˙
S}
+2
3[γ13]{E˙
E}+2
3[diag({˙
E})][γ13]T{E}
+2
3[diag({E})][γ13]T{˙
E}+[γ21]{S2˙
S}
+[diag({˙
E})][γ22]{S2}+2[diag({E})][γ22]
×{S˙
S}+2[diag({E˙
E})][γ23]{S}
+[diag({E2})][γ23]{˙
S}+3
4[υ1]{E2˙
E}
+3
2[diag({E˙
E})][υ1]{E}+3
4[diag({E2})]
×[υ1]T{˙
E}+1
2[diag({˙
E})][υ2][diag({E})]
×{E}+[diag({E})][υ2][diag({E})]{˙
E}
+1
2[diag({˙
E})][υ2]T[diag({E})]{E}
+[diag({E})][υ2]T[diag({E})]{˙
E}(37)
6 RESULTS AND DISCUSSION
6.1 Prediction of non-linear response of the PIC
181 rectangular plate in the first in-plane mode
Experiments have been conducted on a rectangu-
lar piezoelectric plate made of PIC 181 piezoceramic
manufactured by PI Ceramic, Lederhose, Germany.
Fig. 5 Displacement response of the PIC 181 plate for
the applied sinusoidal electric field of 7.5V/mm
peak amplitude (first in-plane mode)
Figure 2 shows the plate with the wires attached on
two faces of the plate coated with electrode mate-
rial. The electric field is applied along the thickness,
and the vibration amplitude of the plate is measured
at the end of its length in the in-plane free-free mode
by using a laser vibrometer (velocity amplitude from
the laser vibrometer is converted to the displacement
amplitude by dividing it with the circular frequency ω
of the corresponding mode of excitation of the plate).
The current flowing through the piezoelectric plate is
measured by connecting a resistance in series with
the plate and measuring the potential drop across
the resistance (the potential difference is converted
to the current by dividing it with the known value
Fig. 6 Current response of the PIC 181 plate for the
applied sinusoidal electric field of 7.5V/mm peak
amplitude (first in-plane mode)
Fig. 7 Displacement response of the PIC 181 plate for
the applied sinusoidal electric field of 0.75V/mm
peak amplitude (first in-plane mode)
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2260 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
of resistance of the resistor). The material properties
used in the analysis are given in Appendix 3 (Fig. 25
for the damping coefficient). These material proper-
ties have been obtained by optimization of the above
closed form solution with experimental results and
the details are given elsewhere [30,31]). The length,
width, and thickness of the plate are 70, 20, and 4 mm,
respectively.
In order to study the non-linear response of the
above-said piezoceramic plate, it was excited with
applied sinusoidal electric fields of different ampli-
tudes. One of the non-linear phenomena observed is
the jump behaviour, which appears when the plate
is excited with different frequencies near the reso-
nance region (e.g. the resonance of the linear sys-
tem is at 23.65 kHz for the first in-plane longitudinal
mode). When the frequency of excitation is gradu-
ally increased from 21 to 25 kHz, the amplitude of
displacement increases gradually (Fig. 5) and then
the amplitude gets increased suddenly (jump for the
Fig. 8 Current response of the PIC 181 plate for the
applied sinusoidal electric field of 0.75V/mm
peak amplitude (first in-plane mode)
Fig. 9 Mode-shape of the piezoceramic plate for the
second free-free in-plane longitudinal vibration
mode
sweep-up mode) at a frequency that depends upon the
excitation field (e.g. jump is at 23.51 kHz for electric
field amplitude of 7.5V/mm). Similarly, when the fre-
quency of excitation is decreased gradually from 25 to
21 kHz, the amplitude of displacement increases grad-
ually (Fig. 5) and then the amplitude gets decreased
suddenly (jump for the sweep-down mode) at a fre-
quency that depends again upon the excitation field
(e.g. jump is at 23.23 kHz for electric field amplitude
of 7.5V/mm). The above jump phenomenon reflects
Fig. 10 Contour plot of three-dimensional shape func-
tions for the second in-plane free-free longitu-
dinal vibration mode of the plate on the original
plate geometry: (a) U(x,y,z); (b) V(x,y,z); (c)
W(x,y,z)
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2261
that the piezoelectric plate behaves like a softening
spring (with cubic non-linear behaviour [29] where
the spring stiffness changes with applied frequency
of excitation), where the response depends upon the
path (i.e. increasing or decreasing frequency ampli-
tude) of excitation. In order to capture this effect, one
must include cubic terms in the material constitutive
equation for the electrical, mechanical, and coupled
piezoelectric domains and the hence the same are
included in the generalized three-dimensional electric
enthalpy function of equation (4).
The results of analysis with a high electric field (viz.
7.5V/mm) andwithalow electricfield (viz.0.75V/mm)
Fig. 11 Displacement response of the PIC 181 plate for
the applied sinusoidal electric field of 7.5V/mm
peak amplitude (second in-plane mode)
Fig. 12 Current response of the PIC 181 plate for the
applied sinusoidal electric field of 7.5V/mm
peak amplitude (second in-plane mode)
have been presented here and compared with those of
experiment. Figures 5 and 6 show the variation of dis-
placement and current amplitudes with frequency at
the electric field amplitude of 7.5V/mm, and Figs 7
and 8 show the same at the electric field amplitude of
0.75V/mm. Results have been shown for both sweep-
up and sweep-down modes of frequency response. As
this piezoceramic material has a low damping prop-
erty, the jump phenomenon can be clearly observed
in Figs 5 and 6 respectively for both displacement
and current responses of the plate. Again, the non-
linear behaviour (i.e. jump phenomenon) is dominant
at higher electric field amplitudes, whereas at the very
Fig. 13 Displacement response of the PIC 181 plate
for the applied sinusoidal electric field of
0.75V/mm peak amplitude (second in-plane
mode)
Fig. 14 Current response of the PIC 181 plate for the
applied sinusoidal electric field of 0.75V/mm
peak amplitude (second in-plane mode)
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2262 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
Fig. 15 (a) Piezoceramic cylinder excited in the longitu-
dinal mode with electric field applied across the
length and the vibration amplitude measured
using laser vibrometer; (b) boundary conditions
used for the axisymmetric modes of the cylinder
small values of applied electric field, the responses of
both sweep-up and sweep-down modes are same and
hence the behaviour is almost linear.
The other non-linearity observed in the frequency
response is the dependence of the resonance fre-
quency on applied electric field amplitude. The linear
resonance frequency for the first in-plane mode of the
plate is 23.65 kHz as can be seen from frequencies cor-
responding to maximum amplitudes of Figs 7 and 8
where it has linear response. However, the frequency
corresponding to maximum amplitude is different for
sweep-up and sweep-down modes (Figs 5 and 6).
This frequency is also less compared with the lin-
ear natural frequency, which corresponds to softening
phenomenon in the piezoceramic and as the applied
electric field increases, this frequency decreases, and
hence for this type of non-linearity, the resonance fre-
quencyisdependent upon theamplitudeof excitation.
It can also be observed that the closed form solution
has been able to represent the different experimental
responses quite satisfactorily.
6.2 Prediction of non-linear response of the PIC
181 rectangular plate in the second in-plane
mode
Experiments have also been carried on the same PIC
181 rectangular plate with electric field excitation near
its second in-plane longitudinal vibration mode (the
deformed mode-shape is shown in Fig. 9). The lin-
ear natural frequency corresponds to 66.2 kHz. For
the analytical solution near the second in-plane mode,
the wave numbers λ1,λ2, and λ3are taken as 3π/l(the
corresponding three-dimensional shape functions, i.e.
U(x,y,z),V(x,y,z), and W(x,y,z), are plotted in
Figs 10(a) to (c), respectively, as contour plots on the
original geometry) and hence the non-linear displace-
ment and current responses are obtained using the
same solution strategy as discussed for the first in-
plane longitudinal vibration mode. Results for applied
electric fields of 7.5 and 0.75V/mm are discussed
here. Figures 11 and 12 show the comparison of
theoretical and experimental displacement and cur-
rent responses, respectively, for the PIC 181 plate
when excited near the second in-plane longitudinal
vibration mode with an electric field of 7.5V/mm.
Experiment and analyses were carried out for both
sweep-up and sweep-down frequency modes of exci-
tation and at this applied electric field, the non-
linearities such as jump phenomena and dependence
of resonance frequency on applied electric field are
clearly observed. The responses at an electric field
of 0.75V/mm is almost linear for both displacement
Fig. 16 Plot of (a) radial and (b) longitudinal displacements corresponding to the eigenvector of
the first axisymmetric longitudinal mode shape of the piezocylinder
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2263
and current responses (Figs 13 and 14). The compar-
ison between experimental and theoretical results is
satisfactory.
6.3 Prediction of non-linear response of the PIC
181 cylinders (10 and 25 mm diameter) in the
first longitudinal mode
After demonstrating the prediction capability of the
solution technique proposed in this work for the non-
linear response of piezoceramic materials with the
help of rectangular piezoplates, the authors intend to
demonstrate the same for another specimen geome-
try (i.e. a solid piezocylinder). They have carried out
Fig. 17 Displacement response of PIC 181 cylinder
(10 mm diameter) for the applied sinusoidal
electric potential of 30V peak amplitude
Fig. 18 Current response of PIC 181 cylinder (10 mm
diameter) for the applied sinusoidal electric
potential of 30V peak amplitude
experiments on piezocylinders of PIC 181 (of 20 mm
length and diameters ranging from 6 to 25 mm). The
electric potential is applied across the length of the
piezocylinder as shown in Fig. 15(a). For the analy-
sis, the authors consider the axisymmetric boundary
conditions (Fig. 15(b)) as the modes under consider-
ation are axisymmetric. The radial and longitudinal
co-ordinates are denoted as rand zand the corre-
sponding displacement are U(r,z)and V(r,z), respec-
tively. In order to derive the mode shape function
corresponding to the first longitudinal vibration
mode of the cylinder numerical simulation (i.e.
finite-element method) is used. The finite-element
linear shape functions corresponding to this mode
are plotted in Fig. 16. Once the shape functions are
Fig. 19 Displacement response of PIC 181 cylinder
(10 mm diameter) for the applied sinusoidal
electric potential of 5V peak amplitude
Fig. 20 Current response of PIC 181 cylinder (10 mm
diameter) for the applied sinusoidal electric
potential of 5V peak amplitude
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2264 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
obtained, the procedure to obtain the displacement
and current response is same as discussed for the plate
geometry. The piezocylinders have been excited with
sinusoidal electric fields of different amplitudes (in
both sweep-up and sweep-down frequency modes).
The amplitude of sinusoidal electric potential across
the length ranges from 5 to 30V.
The results for displacement amplitude (at the end
of the cylinder which is measured with laser vibrom-
eter as shown in Fig. 15(a)) and current responses at
5 and 30V of applied electric potential are presented
andcomparedwith those of theexperiment. Figures 17
and 18 show the non-linear displacement and current
responses for the sweep-up and sweep-down modes
Fig. 21 Displacement response of PIC 181 cylinder
(25 mm diameter) for the applied sinusoidal
electric potential of 30V peak amplitude
Fig. 22 Current response of PIC 181 cylinder (25 mm
diameter) for the applied sinusoidal electric
potential of 30V peak amplitude
of frequency response of the cylinder of 10 mm diam-
eter. It can be observed that non-linearities similar to
those of the plate geometry are also observed here
and predicted results compare very well with those
of experiment. For 5V electric potential amplitude,
the response of the piezocylinder is linear as can be
seen from Figs 19 and 20. The results of non-linear
responses of the piezocylinder of 25 mm diameter are
shown in Figs 21 and 22 corresponding to 30V of
applied electric potential, whereas at 5V of applied
electric potential, the responses are linear (Figs 23
and 24). Again, the predictions compare very well
with those of experiment for the piezocylinder with
different diameters.
Fig. 23 Displacement response of PIC 181 cylinder
(25 mm diameter) for the applied sinusoidal
electric potential of 5V peak amplitude
Fig. 24 Current response of PIC 181 cylinder (25 mm
diameter) for the applied sinusoidal electric
potential of 5V peak amplitude
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2265
The linear natural frequency of the piezocylinders
depends upon their diameter (for the same length
as stiffness and inertia properties change) and hence
the linear natural frequency of the piezocylinder with
10 mm diameter is 66.1 kHz whereas it is 76.91 kHz
for the piezocylinder with 25 mm diameter. When
they are excited with an electric potential of 30V,
it can be observed (Figs 17 and 21) that the fre-
quencies corresponding to resonance peaks in both
sweep-up and sweep-down modes decrease and the
extent of decrease depends upon the cylinder diame-
ter. For example, the sweep-down resonance peak for
the cylinder with 10 mm diameter (Fig. 17) decreased
to 65.4 kHz whereas the same peak for the cylinder
with 25 mm diameter (Fig. 21) decreased to 65.08 kHz.
These results demonstrate the fact that the non-linear
response depends upon the geometrical dimensions
of piezomaterial and the model described here is
able to predict both the displacement and current
responses satisfactorily for the piezomaterial with
different geometries (i.e. plate and cylinders with
different diameters).
7 CONCLUSIONS
Piezoelectric continuum exhibits different types of
non-linearities under weak electric fields (when the
system is operating near resonance frequency) such
as jump phenomena, dependence of resonance fre-
quency on vibration amplitude, etc. In this work,
a generalized non-linear electric enthalpy density
function incorporating higher-order non-linear terms
(quadratic and cubic) in the energy expression of the
coupled piezoelectric medium has been formulated.
Damping has been formulated using proportional
damping method. This non-linear electric enthalpy
density function, as well as the virtual work due to
proportional damping, has been used in the Hamilton
principle to get the variational form of the equa-
tions of motion. This form has been discretized using
Rayleigh–Ritz method.
The shape function of the linear eigenvalue problem
corresponding to the desired mode of vibration of the
rectangular piezoelectric plate has been obtained by
solving the governing equations of motion along with
the boundary conditions for the three-dimensional
domain. These shape functions are used in the
Rayleigh–Ritz method to get the non-linear equation.
The non-linear equation is expanded using the per-
turbation technique of Lindstedt and Poincaré. The
solutionfor the vibrationamplitude has been obtained
finally as the solution of a polynomial equation
where the coefficients are functions of the material
properties. Experiments have been conducted on a
rectangular plate geometry of PIC 181 piezoceramic in
both first and second in-plane modes.The closed form
solutions obtained from three-dimensional analysis
have been compared with the experimental results. It
is observed that the closed form solutions were able
to reproduce the experimental non-linear responses
quite satisfactorily.
Experiments have also been conducted on piezo-
cylinders of different diameters of PIC 181 with electric
field excitations along the length of the cylinder corre-
sponding to the axisymmetric modes. Mode-shapes of
the cylinder for the Rayleigh–Ritz procedure have been
obtained numerically using finite-element method.
Following the same procedure as the piezoplate, the
non-linear displacement and current responses of
the piezocylinder have been predicted for different
amplitudes of electric field excitations. The predicted
results have been compared with experiments and
it was observed that the method could predict sat-
isfactory responses for different types of geometries.
These closed form solutions should form a basis to
validate the finite-element formulations for this non-
linear problem as well as determination of material
parameters using optimization procedure.
REFERENCES
1 Crawley, E. F. Intelligent structures for aerospace: a
technology overview and assessment. AIAA J., 1994,
32(8), 1689–1699.
2 Viswamurthy, S. R., Rao, A. K., and Ganguli, R. Dynamic
hysteresis of piezoceramic stack actuators used in heli-
copter vibration control: experiments and simulations.
Smart Mater. Struct., 2007, 16, 1109–1119.
3 Viswamurthy, S. R. and Ganguli, R. Modeling and anal-
ysis of piezoceramic actuator hysteresis for helicopter
vibration control. Sens. Actuators A, Phys., 2007, 35(2),
801–810.
4 Hwang, S. C. and McMeeking, R. M. A finite element
model of ferroelectric polycrystals. Ferroelectrics, 1998,
211, 177–194.
5 Thakkar, D. and Ganguli, R. Induced shear actuation
of helicopter rotor blade for active twist control. Thin-
Walled Struct., 2007, 45(1), 111–121.
6 Thakkar,D. and Ganguli,R. Use of single crystal and soft
piezoceramics for alleviation of flow separation induced
vibration in smart helicopter rotor. Smart Mater. Struct.,
2006, 15, 331–341.
7 Thakkar, D. and Ganguli, R. Helicopter vibration reduc-
tion in forward flight with induced shear based piezo-
ceramic actuation. Smart Mater. Struct., 2004, 13(3),
599–608.
8 Pagano,N.J. Exact solutionsfor rectangularbidirectional
composites and sandwich plates. J. Compos. Mater., 1970,
4(1), 20–34.
9 Ray, M. C., Rao, K. M., and Samanta, B. Exact analysis of
coupled electro-elastic behaviour of a piezoelectric plate
under cylindrical bending. Comput. Struct., 1992, 45(4),
667–677.
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2266 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
10 Ray, M. C., Rao, K. M., and Samanta, B. Exact solution
for static analysis of an intelligent structure under cylin-
drical bending. Comput. Struct., 1993, 47(6), 1031–1042.
11 Ray,M. C.,Bhattacharya,R.,and Samanta,B. Exact solu-
tions for static analysis of intelligent structures. AIAA J.,
1993, 31(1), 1684–1691.
12 Heyliger, P. and Brooks, S. Exact solutions for laminated
piezoelectric plates in cylindrical bending. Trans. ASME,
J. Appl. Mech., 1996, 63(4), 903–910.
13 Xu, K. M., Noor, A. K., and Tang, Y. Three-dimensional
solutions for coupled thermoelectroelastic response of
multilayered plates. Comput. Methods Appl. Mech. Eng.,
1995, 126(3–4), 355–371.
14 Dube,G. P.,Kapuria,S.,and Dumir,P. C. Exact piezother-
moelastic solution of simply supported orthotropic flat
panel in cylindrical bending. Int. J. Mech. Sci., 1996, 38(3),
1161–1177.
15 Ray, M. C., Bhattacharya, R., and Samanta, B. Exact
solutions for dynamic analysis of composite plates with
distributed piezoelectric layers. Comput. Struct., 1998,
66(6), 737–743.
16 Ray, M. C. Closed-form solution for optimal control of
laminated plate. Comput. Struct., 1998, 69(2), 283–290.
17 Wolf, K. and Gottlieb, O. Nonlinear dynamics of a can-
tilever beam actuated by piezoelectric layers in sym-
metric and asymmetric configuration. Technical report
ETR-2001-02, Faculty of Mechanical Engineering, Israel
Institute of Technology, Israel, 2001.
18 Beige, H. and Schmidt, G. Electromechanical reso-
nances for investigating linear and nonlinear properties
of dielectrics. Ferroelectrics, 1982, 41(1), 39–49.
19 von Wagner, U. and Hagedorn, P. Piezo-beam-systems
subjected to weak electric field: experiment and mod-
elling of nonlinearities. J. Sound Vibr., 2002, 256(5),
861–872.
20 von Wagner, U. Non-linear longitudinal vibrations of
piezoceramics excited by weak electric fields. Int. J.
Non-Linear Mech., 2003, 38(4), 565–574.
21 Neumann, N. Nichtlineare Effekte bei Längsschwingun-
gen axial polarisierter piezokeramischer Staabe: Experi-
mentelle Untersuchungen und Parameteridentifikation.
Diplomarbeit, Institut für Mechanik, Technische Univer-
sität Darmstadt, 2002.
22 Samal, M. K., Seshu, P., Parashar, S. K., von Wagner, U.,
Hagedorn, P., Dutta, B. K., and Kushwaha, H. S. Non-
linear behaviour of piezoceramics under weak electric
fields, part 1: 3-D finite element formulation. Int. J. Solids
Struct., 2006, 43(6), 1422–1436.
23 Samal, M. K., Seshu, P., Parashar, S. K., von Wagner,
U., Hagedorn, P., Dutta, B. K., and Kushwaha, H. S.
Nonlinear behaviour of piezoceramics under weak elec-
tric fields, part 2: numerical results and validation with
experiments. Int. J. Solids Struct., 2006, 43(6), 1437–1458.
24 Samal, M. K., Seshu, P., and Dutta, B.K. An analytical
formulation in 3D domain for the nonlinear response of
piezoelectric slabs under weak electric fields. Int. J. Solids
Struct., 2007, 44, 4656–4672.
25 Usher,T.and Sim,A.Nonlineardynamics of piezoelectric
high displacement actuators in cantilever mode. J. Appl.
Phys., 2005, 98(6), 064102, pp. 1–7.
26 Wang, Q.-M., Zhang, Q., Xu, B., Liu, R., and Cross, L. E.
Nonlinear piezoelectric behavior of ceramic bending
mode under strong electric fields. J. Appl. Phys., 1999,
86(6), 3352–3360.
27 IEEE Standards on Piezoelectricity. ANSI/IEEE Stan-
dard,The Institute of Electrical and Electronic Engineers,
New York, 1988.
28 Maugin, G. A. Nonlinear electromechanical effects and
applications, 1985 (World Scientific Publishing Com-
pany, Singapore).
29 Nayfeh, A. H. and Mook, D. Nonlinear oscillations, 1979
(Willey publications, New York).
30 Parashar, S. K. and von Wagner, U. Nonlinear longitudi-
nal vibrations of transversally polarized piezoceramics:
experiments and modeling. Nonlinear Dyn., 2004, 37(1),
51–73.
31 Samal, M. K. Modelling of nonlinear behaviour of piezo-
electric materials under weak electric fields. MSc Thesis,
Department of Mechanical Engineering, Indian Institute
of Technology, Bombay, 2003.
APPENDIX 1
Boundary conditions for equation (14)
The boundary conditions are derived as
+h/2
−h/2 +b/2
−b/2 c11
∂u
∂x+c12
∂v
∂y+c13
∂w
∂z
+Ez(d31c11 +d32c12 ++d33c13)x=±l/2
dydz=0
+h/2
−h/2 +l/2
−l/2 c21
∂u
∂x+c22
∂v
∂y++c23
∂w
∂z
+Ez(d31c21 +d32c22 ++d33c23)y=±b/2
dxdz=0
+b/2
−b/2 +l/2
−l/2 c31
∂u
∂x+c32
∂v
∂y++c33
∂w
∂z
+Ez(d31c31 +d32c32 ++d33c33)z=±h/2
dxdy=0
+h/2
−h/2 +b/2
−b/2 c66 ∂u
∂y+∂v
∂xx=±l/2
dydz=0
+h/2
−h/2 +l/2
−l/2 c66 ∂u
∂y+∂v
∂xy=±b/2
dxdz=0
+h/2
−h/2 +l/2
−l/2 c44 ∂v
∂z+∂w
∂yy=±b/2
dxdz=0
+h/2
−h/2 +l/2
−l/2 c44 ∂v
∂z+∂w
∂yz=±h/2
dxdy=0
+h/2
−h/2 +b/2
−b/2 c55 ∂u
∂z+∂w
∂xx=±l/2
dydz=0
+b/2
−b/2 +l/2
−l/2 c55 ∂u
∂z+∂w
∂xz=±h/2
dxdy=0
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
Mathematical model to predict non-linear responses of piezoceramic materials 2267
APPENDIX 2
Coefficients of equation (21)
mk=ρ+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [U2(x,y,z)+V2(x,y,z)
+W2(x,y,z)]dxdydz
ck=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{S∗}T[c]{S∗}]dxdydz
Bk=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{S∗}T[c1]{S∗2}]dxdydz
εk=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{S∗2}T[c21]{S∗2}
+{S∗}T[c22]{S∗3}]dxdydz
f(1)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{E∗}T[d][c]{S∗}]dxdydz
f(2)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{E∗}T[γ11]{S∗2}]dxdydz
f(3)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{E∗}T[γ21]{S∗3}]dxdydz
f(4)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 1
2{E∗2}T[γ12]{S∗}dxdydz
f(5)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 [{E∗2}T[γ22]{S∗2}]dxdydz
f(6)
k=+h/2
−h/2 +b/2
−b/2 +l/2
−l/2 1
3{E∗3}T[γ23]{S∗}dxdydz
where
{S∗}=
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
∂U(x,y,z)
∂x
∂V(x,y,z)
∂y
∂W(x,y,z)
∂z
∂W(x,y,z)
∂y+∂V(x,y,z)
∂z
∂W(x,y,z)
∂x+∂U(x,y,z)
∂z
∂V(x,y,z)
∂x+∂U(x,y,z)
∂y
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
;{S∗2}=
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
∂U(x,y,z)
∂x2
∂V(x,y,z)
∂y2
∂W(x,y,z)
∂z2
∂W(x,y,z)
∂y+∂V(x,y,z)
∂z2
∂W(x,y,z)
∂x+∂U(x,y,z)
∂z2
∂V(x,y,z)
∂x+∂U(x,y,z)
∂y2
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
{S∗3}=
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
∂U(x,y,z)
∂x3
∂V(x,y,z)
∂y3
∂W(x,y,z)
∂z3
∂W(x,y,z)
∂y+∂V(x,y,z)
∂z3
∂W(x,y,z)
∂x+∂U(x,y,z)
∂z3
∂V(x,y,z)
∂x+∂U(x,y,z)
∂y3
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
;{E∗}=⎧
⎪
⎨
⎪
⎩
0
0
U0
h
⎫
⎪
⎬
⎪
⎭
;{E∗2}=⎧
⎪
⎪
⎨
⎪
⎪
⎩
0
0
U0
h2⎫
⎪
⎪
⎬
⎪
⎪
⎭
;{E∗3}=⎧
⎪
⎪
⎨
⎪
⎪
⎩
0
0
U0
h3⎫
⎪
⎪
⎬
⎪
⎪
⎭
JMES1002 ©IMechE 2008 Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science
2268 M K Samal, P Seshu, U von Wagner, P Hagedorn, B K Dutta, andHSKushwaha
APPENDIX 3
Fig. 25 Variation of stiffness proportional damping con-
stant with the applied electric field for the
piezoceramic PIC 181
The material properties of PIC 181 are as follows
(Fig. 25)
[SE]=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
1.1600 −0.4070 −0.4996 0 0 0
−0.4070 1.1750 −0.4996 0 0 0
−0.4996 −0.4996 1.4110 0 0 0
0 0 0 3.5330 0 0
0 0 0 0 3.5330 0
0 0 0 0 0 3.1640
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
×10−11 m2/N
[d]=⎡
⎢
⎣
0 0 0 0 3.89 0
0 0 03.8900
−1.148 −1.148 2.53 0 0 0
⎤
⎥
⎦10−10 m/V
[C21]=−
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
4.7000 2.7914 2.6525 0 0 0
2.7914 4.7000 2.6525 0 0 0
2.6525 2.6525 3.9978 0 0 0
0 0 0 0.8465 0 0
00000.8465 0
000000.9452
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
×1016 N/m2
Proc. IMechE Vol. 222 Part C: J. Mechanical Engineering Science JMES1002 ©IMechE 2008
