Piezoelectric Crystals

From grassman, 3 Years ago, written in Plain Text, viewed 73 times.

URL https://paste.bugabuse.net/view/49482205 Embed

Download Paste or View Raw

They can be made out of soda ash and cream of tartar and some water. Make thousands of kV just by hitting a crystal! 
 
More info on Piezoelectrics - 
Piezoelectricity
Jim Emery
4/3/97Contents
0.1  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3
0.2  The Linear Piezoelectric Equations  . . . . . . . . . . . . . . .  3
0.3 Piezoelectric Ceramics Constants and One Dimensional Equations  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  6
0.4  The Stress Form of The Piezoelectric Matrix  . . . . . . . . . .  9
0.5 Three Systems of Notation For The Piezoelectric Equations:
Physics (IRE), ABAQUS, and ANSYS  . . . . . . . . . . . . .  11
0.6 Deriving the Tensor Stress Form of the Piezoelectric Equations
From the Strain Form  . . . . . . . . . . . . . . . . . . . . . .  14
0.7  Orthotropic Material  . . . . . . . . . . . . . . . . . . . . . . .  16
0.8  Poling and Piezoelectric Ceramics  . . . . . . . . . . . . . . . .  18
0.9 The Equality of Direct and Converse Piezo Coeﬃcients . . . . 20
0.10  Typical Values of Piezoelectric Parameters  . . . . . . . . . . .  21
0.11 Piezoelectric Current and Impedance From A Finite Element
Calculation  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22
0.12  Piezoelectric Cube Example  . . . . . . . . . . . . . . . . . . .  24
0.13  Generating Ferroelectric Hysteresis Curves  . . . . . . . . . . .  26
0.14  The Strain Hysteresis Curve  . . . . . . . . . . . . . . . . . . .  26
0.15  Relaxors  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  26
0.16  Electrostriction  . . . . . . . . . . . . . . . . . . . . . . . . . .  27
0.17 The Conversion Program: pzansys.ftn . . . . . . . . . . . . .  27
0.18 Example: Parameter Conversion Using Compliance Components 28
0.19 Example: Parameter Conversion Using Four Independent Elasticity Components  . . . . . . . . . . . . . . . . . . . . . . . .  31
0.20  Example, Orthotropic Elastic Constants  . . . . . . . . . . . .  34
0.21 Program Listing of pzansys.ftn Program  . . . . . . . . . . .  37
0.22 Location of Piezoelectric Information In ANSYS Manuals . . . 58
i0.23 ANSYS Comparison For A Composite Piezoelectric Transducer: Example VM176. . . . . . . . . . . . . . . . . . . . . .  59
0.24  Piezoelectricity Bibliography . . . . . . . . . . . . . . . . . . .  61
ii0.1 Introduction
We shall show how to derive the stress form of the piezoelectric equations
from the strain form, and we shall show how to ﬁnd the clamped permittivity
matrix from the unclamped matrix. These conversions are required for ﬁnite
element programs because in such programs equations are solved for the
elastic displacements, and so the strains rather than the stresses are the
independent variables. We shall do the basic derivation three times in the
following sections. We shall ﬁrst treat the one dimensional equation. In
later sections we shall treat the matrix case and the full tensor case. These
derivations are conceptually the same.
0.2 The Linear Piezoelectric Equations
A piezoelectric ceramic is a ferroelectric material. It exhibits hysteresis and
has nonlinear behavior. When this material is poled (placed in a strong
electric ﬁeld), it attains a permanently polarized state. For a small variation
in the electric ﬁeld, it behaves approximately linearly near that state.
The linear piezoelectric equations are written in two diﬀerent forms. In
the more common form, which is called the strain form, the strain tensor εij
and the electric polarization vector Pi are each written as linear functions of
the electric ﬁeld vector Ei and the stress tensor σij
. These equations are
εij = Ekdkij + c
 
ijkl
σkl
,
and
Pi = 0χijEj + dijkσj k.
The dijk are the components of the piezoelectric tensor (strain coeﬃcients).
The c
 
ijkl
are the components of the inverse elastic tensor. The constant 0
is the permittivity of free space. The χij are the components of the electric
susceptibility tensor. In the direct piezoelectric eﬀect, when a piezoelectric
material is put under a stress, the material becomes electrically polarized
and surface charges appear. The direct piezoelectric eﬀect is
Pi = dijkσj k,
which is obtained from the second equation, when the external electric ﬁeld
is zero.
3In the converse piezoelectric eﬀect, when a piezoelectric material is put
in an external electric ﬁeld (voltages applied to electrodes), the material
experiences a strain. The converse piezoelectric eﬀect is
εij = Ekdkij
,
which is obtained from the ﬁrst equation, when there are no external forces
and the stress is zero.
The fact that the coeﬃcients for both the direct eﬀect and the converse
eﬀect are the same is a consequence of conservation of energy. This fact can
be established by a thermodynamic argument (See Nye). Because of various
symmetries in the tensor indices, the equations have matrix forms. Let us
deﬁne an index function f : (i, j)− > k that takes a pair of indices and
produces a single index. The function f takes values as follows: f(1, 1) =
1, f(2, 2) = 2, f(3, 3) = 3, f(2, 3) = 4, f(1, 3) = 5, f(1, 2) = 6 The function is
symmetric, f(i, j) = f(j, i), and so f(3, 2) = 4, f(3, 1) = 5, f(2, 1) = 6. This
is the Physics and The Institute of Radio Engineers assignment function. We
deﬁne this function with a table:
f 1 2 3
1 1 6 5
2 6 2 4
3 5 4 3
This index assignment is the one most often used in elasticity theory to
convert the second rank stress and strain tensors to six component vectors.
Although this is the usual convention, the following assignment deﬁned by
function g is used sometimes (ABAQUS).
g 1 2 3
1 1 4 5
2 4 2 6
3 5 6 3
A third index function is also used. The function, which we call h is used
by the ﬁnite element program ANSYS.
4h 1 2 3
1 1 4 6
2 4 2 5
3 6 5 3
Warning! These diﬀerences in index assignment will change the matrices
such as c, d, e.
The matrix components of the piezoelectric tensor are di,f(i,j) and are
either equal to the tensor component itself, when i = j, or to twice the
component. This will be explained below.
There is a second form of the piezoelectric equations, which is called the
stress form.
The stress tensor σij and the electric displacement vector Di
, are each
written as linear functions of the electric ﬁeld vector Ei and the strain tensor
εij
. These equations are
σij = cijklεkl − Ekekij
,
and
Di = eijkεj k + ijEj
.
The cijkl are the components of the elastic tensor. The eijk are the components of the piezoelectric tensor (stress coeﬃcients). The ij are the components of the permittivity tensor.
These equations are suitable for ﬁnite element programs because in these
programs the displacements and hence the strains are the independent variables, whereas the boundary conditions are given as force loadings, that is as
stresses. It will be shown below that the piezoelectric stress coeﬃcients and
the strain coeﬃcients are related as follows:
eijk = cj klmdilm.
Notation varies greatly in the ﬁeld of piezoelectricity. In a following
section we will summarize the notation used in physics, and in two ﬁnite
element programs ABAQUS and ANSYS.
50.3 Piezoelectric Ceramics Constants and One
Dimensional Equations
The principal references for this section are: Jaﬀe, Cook and Jaﬀe, Piezoelectric Ceramics, and the Morgan Matroc document Guide To Modern
Piezoelectric Ceramics.
The original research on piezoelectric materials was conducted with parallel plate capacitors, and the problem was treated as one dimensional. So
the strain form of the equations considered were
ε = Ed + c−1
σ
P = 0χE + dσ,
where the quantities are in the 3 direction. From the deﬁnition of the electric
displacement
D = P + 0E = 0(χ + 1)E + dσ = 
σ
E + dσ.
So we have an alternate form
D = 
σ
E + dσ.
Here 
σ
is the permittivity at constant stress (unclamped).
The second set of piezoelectric equations (stress form) may be obtained
from the ﬁrst set (strain form). Solving the ﬁrst equation for the stress, we
have
σ = cε − cdE = cε − eE.
Making this substitution for σ in the alternate form of the second equation
of the ﬁrst set, we get the second equation of the second set.
D = 
σ
E + d(cε − eE) = E + eε,
= (
σ
− de)E + dcε
= 
ε
E + eε,
where the piezoelectric coeﬃcient e = cd. Here 
ε
is the permittivity at
constant stress (unclamped).
6The coeﬃcients may be deﬁned as partial derivatives. For example,
d = (
∂ε
∂E
)σ,
where the subscript σ means that ε is a function of E and σ, where σ is held
constant in computing the partial derivative. Similarly we have
d = (
∂P
∂σ
)E
The constant g is deﬁned by
g = (−
∂E
∂σ
)ε = (
∂ε
∂D
)σ.
We have
E =
1
 
(D − dσ)
so
g = −(
∂E
∂σ
)ε =
d
 
.
The constant e is deﬁned by
e = −(
∂σ
∂E
)ε = (
∂D
∂ε
)E
The poled ceramic is an orthotropic material with a plane of symmetry
whose normal is in the poled direction. Further it is symmetric with respect
to any rotation about the poling direction. Performing a reﬂection through
the symmetry plane, some of the strain components and stress components
are maintained in sign and some others are changed in sign. Writing the
various components of stress as functions of strain, in the two coordinate
systems, we ﬁnd that the proper signs can be maintained only if some of the
elastic coeﬃcients are zero (see pp62-64 Sokolnikoﬀ Mathematical Theory
of Elasticity). Carrying out rotation transformations about the poling axis
(3 direction) we ﬁnd further conditions on the elasticity coeﬃcients. We have
only the following independent elastic coeﬃcients (Jaﬀe et al, p20)
c11, c33, c44, c12, c13.
7We have
c22 = c11, c55 = c44, c66 = 2(c11 − c12), c23 = c13.
Coeﬃcients which are not in the upper 3-dimensional submatrix or on the
main diagonal are zero.
The independent piezoelectric coeﬃcients are
d33, d31, d15.
We have
d32 = d31, d24 = d15.
All other coeﬃcients are zero.
The independent dielectric coeﬃcients are
K11, K33
We have
K22 = K11.
All other coeﬃcients are zero.
The permittivity is
ij = 0Kij
.
The electromechanical coupling constant is deﬁned in terms of the ratio
of stored electrical to stored mechanical energy. Let Qmech be the electrical
energy converted to mechanical energy and let Qinput be the electrical energy
input. Then the electromechanical coupling factor is
k
2
=
Qmech
Qinput
.
It is claimed (Jaﬀe et al p7) that the relation between the unclamped dielectric constant Kσ
and the clamped dielectric constant Kε
is
K
ε
= K
σ
(1 − k
2
).
Because 0 < k  < 1, the clamped dielectric constant is smaller than the
clamped constant
Kε
< Kσ
.
80.4 The Stress Form of The Piezoelectric Matrix
Given the strain form of the piezoelectric equations, we can convert to the
stress form and give expressions for the three independent e tensor parameters. In matrix form we have
ε = d
T
E + c−1
σ
P = 0χE + dσ
From the ﬁrst equation multiplying by the elasticity matrix c we have
cε = cd
T
E + σ,
so solving for the stress matrix σ we have
σ = cε − cd
T
E
= cε − eE,
where the matrix e is deﬁned as the product of c and the transpose of d
e = cd
T
.
We have
e
T
= (cd
T
)
T
= c
T
d = cd
The last equality follows because c is symmetric. From the second equation
D = P + 0E = 0χE + dσ + 0E
= 0(χ + I)E + dσ
= 
σ
E + dσ,
where 
σ
is the permittivity at constant stress (unclamped permittivity).
Substituting the expression for the stress σ from above
D = 
σ
E + d(cε − eE)
= (
σ
− de)E + dcε
9= 
ε
E + e
T
ε,
where
 
ε
= 
σ
− de
is the permittivity at constant strain (clamped permittivity). In the previous
section we saw that the elastic matrix c has 5 independent parameters and
d has three independent parameters. We have
c =








c11 c12 c13 0 0  0
c12 c11 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)








,
and
d =



0 0  0 0 d15 0
0 0  0 d15 0 0
d31 d31 d33 0 0 0


 .
Then the piezoelectric stress matrix is
e = cd
T
=






c11 c12 c13 0 0  0
c12 c11 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)












0 0 d31
0 0 d31
0 0 d33
0 d15 0
d15 0 0
0 0 0






=








0 0 c11d31 + c12d31 + c13d33
0 0 c12d31 + c11d31 + c13d33
0 0 c13d31 + c13d31 + c33d33
0 c44d15 0
c44d15 0 0
0 0  0








10=








0 0 e13
0 0 e13
0 0 e33
0 e42 0
e42 0 0
0 0 0








So the three independent piezo e parameters are e13, e33, e42. (Warning! The
e matrix is sometimes deﬁned as the transpose of this e matrix.)
Explicitly we have
e13 = (c11 + c12)d31 + c13d33
e33 = (c13 + c13)d31 + c33d33
and
e42 = c44d15.
These equations are easily solved for d33, d31 and d15 in terms of the three
independent e and ﬁve independent c coeﬃcients.
0.5 Three Systems of Notation For The Piezoelectric Equations: Physics (IRE), ABAQUS,
and ANSYS
The index function f(ij), which maps two indices to one, in the way that
is conventional in elasticity theory, was deﬁned in the previous section, as
was g(ij). Here is a table showing the notation most commonly used in
physics (References, Nye, Auld, Cady, IRE (The Institute of Radio Engineers Standard) ANSYS Procedures p8-15, ), and in the ﬁnite element systems ABAQUS (Reference: Theory manual, section 3.10.1), and ANSYS
(Reference: Theory manual, Piezoelectrics chapter).
11Vari abl e Physics ABAQUS ANSYS
Electric Field Vector Ei Ei Ei
Electric Polarization Vector Pi Pi Pi
Electric Displacement Vector Di
qi Di
Stress Tensor σij σij σij
Stress Vector σk = σf(ij)
- Tk
Strain Tensor εij εij
-
Strain Vector εk - Sk
Electric Permittivity Tensor ij D
φ
ij
ij
Electric Susceptibility Tensor χij
- -
Elasticity Tensor cijkl Dijkl
cijkl
Elasticity Matrix cmn = cf(ij)f(kl) − cmn
Inverse Elasticity Tensor c
 
ijkl
- c
 
ijkl
Piezoelectric Tensor (strain) di,jk d
φ
i,jk
di,jk
Piezoelectric Matrix (strain) di,m = di,f(j k) di,m = di,g(j k) di,m
Piezoelectric Tensor (stress) ei,jk e
φ
i,jk
ej k,i
Piezoelectric Matrix (stress) ei,m = ei,f(j k)
ei,m = ei,g(j k)
em,i
The engineering strain coeﬃcients are written as γij
. They diﬀer from
the εij
, in that the shear coeﬃcients are doubled in value.
Note: the physics notation for eijk is deﬁned in equation 8.40 of Auld. Nye
does not introduce the piezoelectric e coeﬃcients. The Cady book is very
old and one can obtain piezoelectric constant notation only by implication.
Note: The ANSYS deﬁnition of the e matrix and tensor is the transpose of
the physics and ABAQUS deﬁnition. Thus the ANSYS e matrix is 6 by 3,
whereas the physics e matrix is 3 by  6.
The two ABAQUS tensor equations are:
σij = Dijklεkl − e
φ
mijEm
(or alternately: σij = Dijklεkl − Dijkld
φ
mkl
Em)
qi = e
φ
ijk + D
φ
ijEj
.
The ANSYS matrix equations are:
T = cS − eE
D = e
T
S + E.
12The Physics matrix equations are:
σ = cε − e
T
E,
D = eε + E.
Again note that the ANSYS e deﬁnition and the physics e deﬁnition are
transposes of one another.
Here are the ANSYS matrices written out:
T = σ =





σ1
σ2
...
σ6





S = ε =





ε1
ε2
...
ε6





E =



E1
E2
E3



D =



D1
D2
D3



c =








c11 c12 c13 0 0 0
c21 c22 c23 0 0 0
c31 c32 c33 0 0 0
0 0  0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66








13e =






e11 e12 e13
e21 e22 e23
e31 e32 e33
e41 e42 e43
e51 e52 e53
e61 e62 e63






 =



11 0 0
0 22 0
0 0 22



We will show in the next section that (physics notation)
eijk = cj klmdilm.
0.6 Deriving the Tensor Stress Form of the
Piezoelectric Equations From the Strain
Form
We have already done this derivation in a previous section for the matrix
forms of the piezoelectric equations. Here we will repeat this for the tensor
forms of the equations. We start with the strain form of the piezoelectric
equations
εij = Ekdkij + c
 
ijkl
σkl
,
and
Pi = 0χijEj + dijkσj k.
Let cpqij be the elasticity tensor. Because c
 
is its inverse, we have
cpqij
c
 
ijkl = δ
p
k
δ
q
l
.
Then multiplying the ﬁrst equation by cpqij and summing we get
cpqijεij = cpqijEkdkij + cpqij
c
 
ijkl
σkl
,
14which becomes
cpqijεij = cpqijEkdkij + σpq.
We let
ekpq = cpqijdkij
.
Then
σpq = cpqrsεrs − empqEm.
The equation relating electric displacement D, electric polarization P, and
the macro electric ﬁeld E is
Di = Pi + 0Ei
.
Substituting for the polarization from the second piezoelectric equation above,
we have
Di = 0χijEj + dijkσj k + 0Ei
.
We write this as
Di = 0(χij + δij
)Ej + dijkσj k
= 
σ
ijEj + dipqσpq,
where 
σ
ij
is the zero stress (unclamped ) permittivity. Substituting our previously derived expression for the stress, we get
Di = 
σ
imEm + dipq(cpqrsεrs − empqEm)
= (
σ
im − dipqempq)Em + dipqcpqrsεrs
= eirsεrs + 
ε
imEm,
where
 
ε
im = 
σ
im − dipqempq
is the zero strain (clamped) permittivity. We have assumed that cpqrs = crspq.
This is the case when the internal energy of the material is a function of the
values σ, ε, D, and E, and not dependent on the path through which the
state was reached, that is,when the internal energy is an exact diﬀerential in
the thermodynamic sense. Therefore redeﬁning the dummy indices, we have
the stress version of the piezoelectric equations:
σij = cijklεkl − Ekekij
,
15and
Di = eijkεj k + 
ε
ij
Ej
.
The matrix forms of these equations are:
σ = cε − e
T
E,
D = eε + 
ε
E.
0.7 Orthotropic Material
For the basic equations of elasticity, see Jim Emery Elasticity, (See ﬁles:
elastic.tex, elastic.ps).
An orthotropic problem is one in which the material properties are symmetric with respect to 3 mutually orthogonal planes (See Sokolnikoﬀ p62). A
piezoelectric ceramic such as PZT is orthotropic, but has even more symmetry. It is invariant to any rotation about the poled axis. By applying various
symmetry transformations, and invariance properties, one ﬁnds that many
of the elastic constants must be zero. The elastic matrix becomes
c =








c11 c12 c13 0 0 0
c21 c22 c23 0 0 0
c31 c32 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66








.
The inverse matrix is
c
 
= c−1
=








1/E1 −ν21/E2 −ν31/E3 0 0  0
−ν12/E1 1/E2 −ν32/E3 0 0  0
−ν13/E1 −ν23/E2 1/E3 0 0  0
0 0 0 1/(2µ32) 0  0
0 0 0 0 1/(2µ31) 0
0 0 0 0 0 1/(2µ12)








.
The Poisson ratio, νij
, is the ratio of the negative of the strain in the
xj direction, to the strain in the xi direction, for a normal stress in the xi
direction. For example, if the only nonzero stress component is σ1, then
16





1/E1 −ν21/E2 −ν31/E3 0 0 0
−ν12/E1 1/E2 −ν32/E3 0 0 0
−ν13/E1 −ν23/E2 1/E3 0 0 0
0 0 0 1/(2µ32) 0  0
0 0 0 0 1/(2µ31) 0
0 0 0 0 0 1/(2µ12)












σ1
0
0
0
0
0






=








σ1/E1
(−σ1ν12)/E1
(−σ1ν13)/E1
0
0
0








=






ε1
ε2
ε3
0
0
0






.
So that
ν12 =
−ε2
ε1
,
and
ν13 =
−ε3
ε1
,
By the symmetry of c, we have
νij/Ei = νj i/Ej
.
So there are three independent values of Poisson’s ratio. We can take these
to be
ν12, ν13, ν23,
or
ν21, ν31, ν32.
17For the case of poled piezoelectric ceramics poled in the 3 direction, the best
choice is
ν21, ν31, ν32.
It is clear that in that case
ν31 = ν32,
so that there are only two independent poisson ratio’s, which we may take
as
ν21, ν31.
Also we see that ν31 is the poisson ratio stated for piezoceramics and is said
to be about
.31
If the engineering strain vector is used then the shear strains are twice the
physics shear strains, so that the 2 multiplying a shear modulus is suppressed.
The engineering shear strain is equal to the shear angle, the physics shear
strain is equal to one half the angle.
The notation for shear modulus µij depends upon which of the index
mappings f or g is used in assigning the stress-strain tensor indices to the
vector indices.
0.8 Poling and Piezoelectric Ceramics
See:
Newnham R E, Trolier-McKinstry S,Giniewicz J R
Piezoelectric, Pyroelectric and Ferroic Crystals
J. Material Education, 15, 189-223 (1993),
(reprint: Penn State MRL annual report to ONR 1994, Appendix 1)
For the nature of the strain vs electric ﬁeld curve, the nature of domains
in PZT, and the morphotropic phase boundary between tetragonal PbTiO3
and rhombohedral PbZrO3 ferroelectric phases, see:
Subbarao E C, Srikanth V,Cao W, Cross L E
Domain Switching And Microcracking During
Poling Of Lead Zirconate Titonate Ceramics
18Ferroelectrics 1993 V 145 pp. 271-281
(reprinted: Materials For
Adaptive Structural Acoustic Control, Office
Of Navel Research, Penn State Materials
Research Laboratory, L. Eric Cross ed, V1 1994, appendix 15).
The morphrotropic phase boundary is the boundary between two different structures of the piezoelectric material. We note that “morph” is a
root meaning form, and “trop” is a root meaning to turn or change. Hence
morphotropic means: a change in form.
Consider the equation
Pi = dijσj
.
Let the poling direction be the z axis. Then the domains will be randomly
directed in the xy plane. Hence the ceramic will be orthotropic relative to
the xy plane. Because of this symmetry, the piezoelectric matrix takes the
form (see Jaﬀe, Cook, Jaﬀe Piezoelectric Ceramics, p20)
d =



0 0  0 0 d15 0
0 0  0 d24 0 0
d31 d32 d33 0 0 0


 .
=



0 0  0 0 d15 0
0 0  0 d15 0 0
d31 d31 d33 0 0 0


 .
The elasticity tensor takes the form
c =








c11 c12 c13 0 0  0
c12 c22 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)








.
The dielectric tensor is 


K1 0 0
0 K1 0
0 0 K3


 .
19Suppose we have d31 = d32, d33, d15 = d24,and dij = 0,  otherwise.  Then  we
have
d311 = d31,
d322 = d32,
d333 = d33,
d113 = d15/2,
d131 = d15/2,
d223 = d24/2,
d232 = d24/2,
and dijk = 0  otherwise.
0.9 The Equality of Direct and Converse Piezo
Coeﬃcients
This equality may be shown by applying the ﬁrst and second laws of thermodynamics. Let the stress coeﬃcients σij
, the electric ﬁeld components Ei and
the temperature T be the state variables. Then write diﬀerential expressions
for the strain deij
, the electric displacement dDi
, and the entropy dS. Write
the internal energy change as
dU = σijdeij + EidDi + T dS.
Introduce the function
Φ = U − σij eij − EiDi − T S.
Compute the diﬀerential of Φ using the last two equations, and also write it
in the general partial derivative form. Then equate the partial derivatives to
coeﬃcients of the equations. For example
∂Φ
∂Ei
= −Di
.
Diﬀerentiate to equate derivatives of eij
, Di
, or S to second order partial
derivatives of Φ. This shows for example that the coeﬃcients in the direct
piezoelectric eﬀect are equal to the coeﬃcients of the converse piezoelectric
eﬀect. See Nye p179.
200.10 Typical Values of Piezoelectric Parameters
Dielectric constant
K = 200 =
 
0
0 = 8.85 × 10
−12
farad/meter.
KT
dielectric constant at constant stress. KS
dielectric constant at zero
strain (constrained piezoelectric coeﬃcient)
d33 = 200 × 10−12
coulomb/newton (or meter/volt).
d31 = −100 × 10−12
coulomb/newton (or meter/volt).
d15 = 400 × 10−12
coulomb/newton (or meter/volt).
gij voltage coeﬃcients.
kij = gij
Poisson’s ratio is approximately .31 for all ceramics. Curie temperature
Tc = 300.
fr = resonant frequency. Typical value is 70 kHz.
Panasonic motor.
Operating voltage 22 volts.
Rated output 1 watt.
Torque 250 gf.cm (gf gram force)
250 gf cm (pound/454 gf)(inch/2.54 cm)(16 ounce/pound) = 3.47 inch.ounce
Current .4 amp at 5 to 10 volts
Input power 2 to 4 watts
Voltage of control circuit 12 volts
250 g.cm (kg/1000 g)(9.8 m/s.s)(m/100 cm) = .0245 N.m
Angular velocity = 4002π/(min)(min/60sec) =  10.47(1/s)
Power = (.0245)(10.47) = .2565 watt
Speed 400 rpm, noload 600 rpm.
fa = antiresonant frequency (parallel resonance)
Electromagnetic coupling constant
k = (electrical energy stored)/(mechanical energy stored).
21kp = .2 to .8 = planar coupling factor of a thin disk.
Qm = quality factor.
Example: Thickness
l = .030inches = 7.63 × 10−4
m
v = 400volts.
E = v/l = 3.67 × 10
6
v/m.
e = dE = (200 × 10−12
)(3.67 × 10
6
) = 7.34 × 10−4
.
∆z = el = 22µinches
0.11 Piezoelectric Current and Impedance From
A Finite Element Calculation
Let D be the electric displacement and ρ the charge  density.  We can  use one
of Maxwell’s equations
∇ · D = ρ,
to show that the surface charge density on a conductor equals the normal
component Dn of D outside of the conductor. We enclose a conducting
surface with a thin pillbox of area A, volume V , and thickness ∆t. The
charge contained in the volume is Aτ where τ is the surface charge density.
Then letting ∆t go to zero, we have
Aτ =
 
V
ρdV =
 
V
∇ · DdV =
 
∂ V
D · nds → D · nA.
Then the charge density is
τ = D · n = Dn.
We compute the nodal strains εij
from the nodal displacements ui at a
face of an element.
From the ﬁnite element nodal potentials φn, we compute derivatives of
the shape functions, and then the gradient of the potential. Then we obtain
the electric ﬁeld Ei as the negative gradient of the potential,
E = −∇φ.
22We compute the electric displacement from the piezoelectric equation
Di = eijkεj k + KijEj
.
Then we compute the normal component Dn. This gives the surface charge
density. Multiplying by the area A, we get the charge q. We sum up the
charges on the boundaries of the elements to get the total electrode charge.
If we calculate the charge as a function of time, we can ﬁt it to a function
Q0 cos(ωt + ψ),
and thus to the real part of the exponential
Q0 exp(j(ωt + ψ)).
Diﬀerentiating with respect to time we get the current
I = I0 exp(j(ωt + ψ)).
Then if V is the applied voltage, the impedance is
Z = V /I .
We can also compute energy storage and dissipation in the ﬁnite element
model, and relate it to energy storage and dissipation in an equivalent circuit.
We can obtain or compute the energy density of the electric ﬁeld
D · E
2
,
the elastic strain energy, the kinetic energy, heat dissipation, and frictional
contact energy. The kinetic energy will correspond to inductor energy in
the equivalent circuit, the strain energy and the electric ﬁeld energy will
correspond to capacitive energy, and the dissipation energy will correspond
to a resistance loss. If the energies correspond, then the equivalent circuit is
valid.
230.12 Piezoelectric Cube Example
Consider a piezoelectric cube of side a. Let thin metal plates be located at
z = 0  and z = a. Let there be a voltage V on these plates. Let there be a
horizontal compressive forces F in the y direction acting on the y = 0,  and
the y = a faces of the block. The electric ﬁeld is
E =



E1
E2
E3


 =



0
0
V /a



The stress is
σ =






σ1
σ2
σ3
σ4
σ5
σ6






=






F /a
2
0
0
0
0
0






.
An average dielectric constant for PZT is
K = 200 =
 
o
= 1 + χ.
We have
0 = 8.85 × 10
10
(C
2
/N M2
).
Then the permittivity is
 = 1.77 × 10−9
.
If we take the material to be isotropic, the polarization equation becomes
P = 0χE +



0 0 0 0 d15 0
0 0 0 0 d25 0
d31 d32 d33 0 0  0











σ1
σ2
σ3
σ4
σ5
σ6








We have
d31 = d32 = −100 × 10−12
C/N or M/V
24d33 = 2100 × 10−12
d51 = d15 = 400 × 10−12
Because only σ1 is nonzero, we have
P = 0χE +



0
0
d31σ1


 .
Then
P3 = 0χE3 + d31σ1
D = 0E + P.
D3 = 0(1 + χ)E3 + d31σ1 = 0K E3 + d31σ1.
The surface charge on the z = a face is
Q = D3a
2
.
The strain is
e =








0 0 d31
0 0 d32
0 0 d33
0 d25 0
d15 0 0
0 0 0











E1
E2
E3


 +








c11 c12 c13 0 0 0
c21 c22 c33 0 0 0
c31 c32 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66
















σ1
σ2
σ3
σ4
σ5
σ6








So the nonzero elements of the strain are
e1 = d31E3 + σ1c11 = d31E3 +
σ1
ξ
,
e2 = d32E3 + σ1c21 = d32E3 −
σ1ν
ξ
,
e3 = d33E3 + σ1c31 = d33E3 −
σ1ν
ξ
,
where ξ is Young’s modulus, and ν is Poisson’s ratio.
250.13 Generating Ferroelectric Hysteresis Curves
Suppose we arrange two capacitors in series, one C1 containing a ferroelectric
dielectric, and the second C2 a nonferroelectric dielectric. Let us connect in
series an alternating high voltage source. Let V1 be the voltage across capacitor C1. This is a measure of the E ﬁeld in the ferroelectric capacitor C1. The
charge on the plates of C1 is proportional to the D ﬁeld (actually the D ﬁeld
is equal to the surface charge density). But the charge on C1 is equal to the
charge on C2 and moves from C1 to C2 and back again during a cycle of the
applied alternating voltage. But the second capacitor is nonferroelectric and
its dielectric is linear, so that the voltage V2 across the second capacitor is
proportional to charge and hence to the D ﬁeld of the ferroelectric capacitor.
Therefore, if we connect V1 across the horizontal plates of an oscilloscope,
and V2 across the vertical plates, we will see the hysteresis curve. The polarization ﬁeld will be approximately equal to D, since the permittivity  for
the ferroelectric is much larger than the permittivity of free space 0.
(I think this is close to the method given by Sawyer and Tower, see
bibliography.)
0.14 The Strain Hysteresis Curve
The strain vs electric ﬁeld curve for a piezoelectric ceramic is a butterﬂy
curve with the strain always of one sign, and the curve lying in the positive
upper half plane.
0.15 Relaxors
A relaxor ferroelectric has some of the properties of a normal ferroelectric,
but is not so permanently polarized and has a very narrow hysteresis loop.
When the ﬁeld is removed, the material almost returns to an unpolarized
state. It requires a bias ﬁeld.
See:
Cross L Eric
Relaxor Ferroelectrics: An Overview
Ferroelectrics, 1994 v151 pp. 305-320
26(reprinted: annual report 1994 from MRL at Penn State to ONR v1 p305)
0.16 Electrostriction
Electrostriction is an internal stress caused by the force of the electric ﬁeld on
charges. It occurs in all materials, not just piezoelectric ones. A computation
shows it to be related to the gradient of the dielectric constant, and in practice
the stress is proportional to the square of the electric ﬁeld intensity. It is very
hard to measure, since the force can be masked by the force due to the free
charges on the capacitor plates.
Notes and references to electrostriction: Wert, Thomson Physics of
Solids, p. 321. Julius Adams Stratton, Electromagnetic Theory, pp149-
151. Panofsky and Philips Classical Electricity and Magnetism p.95,
As a term of an equation. Maxwell stress tensor. Classical derivation. The
treatment in Panofsky and Phillips is based on the treatment in Becker:
(formerly called Becker and Abraham) Electromagnetic Fields and Interactions p. 134. In Jaﬀe: Piezoelectric Ceramics, Electrostriction is
mentioned on pages 15 and 78. He relates electrostriction to domain reversal. Jaﬀe says that the eﬀect is signiﬁcant only for ferroelectrics above the
Curie point. Cady: Piezoelectricity deﬁnes electrostriction to be that effect that is proportional to E
2
. See pages 4, 198, 199, 614. Cady states
that electrostriction is only signiﬁcant for ﬁelds above 20,000 volts per cm.
See also Auld, Acoustic Fields and Waves in Solids, volumes I and II,
who presents a nice one dimensional model of the piezoelectric eﬀect on p265.
0.17 The Conversion Program: pzansys.ftn
This program does the conversion calculations derived in the previous sections. It requires as input the d parameters, and either the elasticity matrix,
or its inverse. It computes the e parameters, and computes the clamped
permittivity (strain constant) from the unclamped permittivity. It generates
material property cards for the ANSYS program. We list the program in a
later section. The next sections illustrate the operation of the program.
270.18 Example: Parameter Conversion Using
Compliance Components
yonada:/users/u51195 $ pzansys
Reference: Jim Emery "Piezoelectricity",
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
Enter the orthotropic elastic constants
for a coordinate system in the principle directions
Enter Young’s Modulus, in x [ 0.12800E+12]
28Enter Young’s Modulus, in y [ 0.12800E+12]
Enter Young’s Modulus, in z [ 0.11000E+12]
Enter ANSYS Poisson ratios (Theory sec. 2.1)
Enter Major Poisson ratio prxy[ 0.30000 ]
Enter Major Poisson ratio prxz[ 0.41976 ]
Enter Major Poisson ratio pryz[ 0.41976 ]
Enter shear modulus g23 [ 0.52000E+11]
Enter shear modulus g13 [ 0.52000E+11]
Enter shear modulus g12 [ 0.12200E+12]
Compliance matrix (elasticity inverse) cinv=
0.7813E-11 -0.2344E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.2344E-11 0.7813E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3279E-11 -0.3279E-11 0.9091E-11 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Note: In ANSYS components 4,5,6 are permuted
(A different mapping from tensor to vector is used):
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
The five independent compliance coeficients are:
s11= 0.78125000E-11
s33= 0.90909091E-11
s44= 0.19230769E-10
s12=-0.23437500E-11
s13=-0.32793388E-11
Computing inverse matrix
Elasticity matrix =
0.2104E+12 0.1119E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1119E+12 0.2104E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1163E+12 0.1163E+12 0.1939E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.12800000E+12
Young’s Modulus 2= 0.12800000E+12
Young’s Modulus 3= 0.11000000E+12
Major Poisson ratio prxy= 0.30000000
Major Poisson ratio prxz= 0.41975537
Major Poisson ratio pryz= 0.41975537
Shear Modulus g23= 0.52000000E+11
Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
29Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -8.977066554171117
e33= 15.52345455897376
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.4900E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permitivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permitivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.3950E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.12800000E+12
EY, 2, 0.12800000E+12
EZ, 2, 0.11000000E+12
30/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.30000000
PRXZ, 2, 0.41975537
PRYZ, 2, 0.41975537
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.39498958E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -8.9770666
TBDATA, 6, -8.9770666
TBDATA, 9, 15.523455
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.19 Example: Parameter Conversion Using
Four Independent Elasticity Components
yonada:/users/u51195 $ pzansys
Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals, and so on)
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
31c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
2
Enter c11 [ 0.13200E+12]
Enter c33 [ 0.11500E+12]
Enter c44 (ANSYS c66) [ 0.52000E+11]
Enter c12 [ 0.71000E+11]
Enter c13 [ 0.73000E+11]
Elasticity matrix c=
0.1320E+12 0.7100E+11 0.7300E+11 0.0000E+00 0.0000E+00 0.0000E+00
0.7100E+11 0.1320E+12 0.7300E+11 0.0000E+00 0.0000E+00 0.0000E+00
0.7300E+11 0.7300E+11 0.1150E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Note: In ANSYS components 4,5,6 are permuted
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
Compliance matrix=
0.1273E-10 -0.3665E-11 -0.5754E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3665E-11 0.1273E-10 -0.5754E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.5754E-11 -0.5754E-11 0.1600E-10 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.78561263E+11
Young’s Modulus 2= 0.78561263E+11
Young’s Modulus 3= 0.62497537E+11
Major Poisson ratio prxy= 0.28788955
Major Poisson ratio prxz= 0.45203533
Major Poisson ratio pryz= 0.45203533
Shear Modulus g23= 0.52000000E+11
Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
32Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 8.400
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -5.700000000000001
e33= 8.400000000000000
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 8.400
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.2820E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permitivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permitivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.6030E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.78561263E+11
EY, 2, 0.78561263E+11
EZ, 2, 0.62497537E+11
33/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.28788955
PRXZ, 2, 0.45203533
PRYZ, 2, 0.45203533
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.60300000E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -5.7000000
TBDATA, 6, -5.7000000
TBDATA, 9, 8.4000000
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.20 Example, Orthotropic Elastic Constants
Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals, and so on)
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
34(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
Enter the orthotropic elastic constants
for a coordinate system in the principle directions
Enter Young’s Modulus, in x [ 0.12800E+12]
Enter Young’s Modulus, in y [ 0.12800E+12]
Enter Young’s Modulus, in z [ 0.11000E+12]
Enter ANSYS Poisson ratios (Theory sec. 2.1)
Enter Major Poisson ratio prxy[ 0.30000 ]
Enter Major Poisson ratio prxz[ 0.41976 ]
Enter Major Poisson ratio pryz[ 0.41976 ]
Enter shear modulus g23 [ 0.52000E+11]
Enter shear modulus g13 [ 0.52000E+11]
Enter shear modulus g12 [ 0.12200E+12]
Compliance matrix (elasticity inverse) cinv=
0.7813E-11 -0.2344E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.2344E-11 0.7813E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3279E-11 -0.3279E-11 0.9091E-11 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Note: In ANSYS components 4,5,6 are permuted
(A different mapping from tensor to vector is used):
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
The five independent compliance coeficients are:
s11= 0.78125000E-11
s33= 0.90909091E-11
s44= 0.19230769E-10
s12=-0.23437500E-11
s13=-0.32793388E-11
Computing inverse matrix
Elasticity matrix =
0.2104E+12 0.1119E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1119E+12 0.2104E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1163E+12 0.1163E+12 0.1939E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.12800000E+12
Young’s Modulus 2= 0.12800000E+12
Young’s Modulus 3= 0.11000000E+12
Major Poisson ratio prxy= 0.30000000
Major Poisson ratio prxz= 0.41975537
Major Poisson ratio pryz= 0.41975537
Shear Modulus g23= 0.52000000E+11
35Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -8.977066554171117
e33= 15.52345455897376
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.4900E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permittivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permittivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.3950E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.12800000E+12
EY, 2, 0.12800000E+12
EZ, 2, 0.11000000E+12
/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.30000000
36PRXZ, 2, 0.41975537
PRYZ, 2, 0.41975537
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.39498958E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -8.9770666
TBDATA, 6, -8.9770666
TBDATA, 9, 15.523455
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.21 Program Listing of pzansys.ftn Program
Here is a listing of the program. The program works through tensor summations, matrix calculations, and matrix inversions. Each subroutine has a one
line description and a description of the parameters.
% cat pzansys.ftn
c pzansys.ftn stress and strain versions of the piezoelectric tensor
c Version 2/6/97
c Author: Jim Emery (AlliedSignal Kansas City)
c for the ansys finite element program.
column 11111111112222222222333333333344444444445555555555666666666677777777778
c2345678901234567890123456789012345678901234567890123456789012345678901234567890
c Reference: Jim Emery "Piezoelectricity" (File: piezoelc.tex)
implicit real*8 (a-h,o-z)
dimension d(3,3,3)
dimension e(3,3,3)
dimension c(3,3,3,3)
dimension dc(3,3)
dimension cmat(6,6)
dimension cmatinv(6,6)
dimension a(6,6),b(6,6),cc(6,6)
dimension dmat(3,6)
dimension emat(6,3)
dimension epsilon(3,3)
dimension clamped(3,3)
dimension ain(10)
dimension pr(3)
dimension ym(3)
real*8 nu12,nu23,nu13,nu21,nu32,nu31
character*1 cm
character*1 reply
character*53 cs1
zero=0.
c piezoelectric constants (coulomb/newton)
37do i=1,3
do j=1,6
dmat(i,j)=0.
enddo
enddo
dmat(3,3)=200e-12
dmat(3,1)=-100e-12
dmat(1,5)=400e-12
write(*,*)’Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)’
write(*,*)’Document location:/u51195/ps/piezoelc.ps,piezoelc.tex ’
write(*,*)’Anonymous FTP: ftp.os.kcp.com /pub/emery/ ’
write(*,*)’Units are metric: Newton, Meter, Kilogram’
write(*,*)’stress: Newton/M^2 (Pascals)’
write(*,*)’Young’’s Modulus: Pascals’
write(*,*)’Elasticity coefficients: Pascals, and so on)’
write(*,*)’Enter piezoelectric d constants (coulomb/newton)’
write(*,*)’The IRE convention for piezoelectric labeling is used’
write(*,*)’Type <return> for default value’
write(*,*)’Enter d33 [’,dmat(3,3),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(3,3)=ain(1)
endif
write(*,*)’Enter d31 [’,dmat(3,1),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(3,1)=ain(1)
endif
dmat(3,2)=dmat(3,1)
write(*,*)’Enter d15 [’,dmat(1,5),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(1,5)=ain(1)
endif
dmat(2,4)=dmat(1,5)
write(*,*)’ d matrix= ’
do i=1,3
write(*,’(6(g11.4,1x))’)(dmat(i,j),j=1,6)
enddo
c construct the components of the d tensor
do i=1,3
do j=1,3
do k=1,3
call t2vec(j,k,jk)
if(j .eq. k)then
d(i,j,k)=dmat(i,jk)
else
d(i,j,k)=dmat(i,jk)/2.
endif
enddo
enddo
enddo
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 27 d tensor components? [n]’
read(*,’(a)’)reply
38if(lenstr(reply) .eq. 0)then
reply=’n’
endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
if(d(i,j,k) .ne. zero)then
write(*,’(a,i1,i1,i1,a,g15.8)’)’d(’,i,j,k,’)=’,d(i,j,k)
endif
enddo
enddo
enddo
endif
c orthotropic elastic constants
write(*,*)’Define the elasticity matrix or compliance matrix:’
write(*,*)
write(*,*)’(1)Use elastic constants to define compliance’
write(*,*)’for orthotropic material, using Young’’s moduli,’
write(*,*)’Poissons ratios, and Shear Moduli.’
write(*,*)
write(*,*)’(2)Use the five components of the elastic matrix’
write(*,*)’c11,c33,c44 (c44 is ANSYS c66),c12,c13, ’
write(*,*)’for a z-poled ceramic.’
write(*,*)’(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,’
write(*,*)’and c_{23}=c_{13}. Coefficients which are’
write(*,*)’not in the upper 3-dimensional submatrix,’
write(*,*)’not on the main diagonal, are zero.’
write(*,*)
write(*,*)’(3)Use the five components of the compliance matrix’
write(*,*)’s11,s33,s44 (s44 is ANSYS s66),s12,s13,’
write(*,*)’for a z-poled ceramic.’
write(*,*)’(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,’
write(*,*)’and s_{23}=s_{13}. Coefficients which are’
write(*,*)’not in the upper 3-dimensional submatrix,’
write(*,*)’not on the main diagonal, are zero.’
write(*,*)’type <return> for default selection: [1]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’1’
endif
c2345678901234567890123456789012345678901234567890123456789012345678901234567890
c begin case1
if(reply .eq. ’1’)then
write(*,*)’Enter the orthotropic elastic constants’
write(*,*)’for a coordinate system in the principle directions’
e1=12.8e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in x [’,e1,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
39e1=ain(1)
endif
e2=12.8e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in y [’,e2,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
e2=ain(1)
endif
e3=11.0e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in z [’,e3,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
e3=ain(1)
endif
c write(*,’(a)’)’NOTE. The poisson ratio nuij is the ’
c write(*,’(a)’)’the strain ratio -epsilon_j/epsilon_i ’
c write(*,’(a)’)’under a normal stress sigma_i’
write(*,*)
c nu31=.31
write(*,*)’Enter ANSYS Poisson ratios (Theory sec. 2.1)’
prxy=.3
prxz=.41975537
pryz=.41975537
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio prxy[’,prxy,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxy=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio prxz[’,prxz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxz=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio pryz[’,pryz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxz=ain(1)
endif
c g23=(e1*e2)/(e1+e2+2.*nu12*e1)
g23=5.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g23 [’,g23,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g23=ain(1)
endif
c g13=(e1*e2)/(e1+e2+2.*nu12*e1)
40g13=5.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g13 [’,g13,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g13=ain(1)
endif
c g12=(e1*e2)/(e1+e2+2.*nu12*e1)
g12=12.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g12 [’,g12,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g12=ain(1)
endif
c build the compliance matrix
do i=1,6
do j=1,6
cmatinv(i,j)=0.
enddo
enddo
cmatinv(1,1)=1./e1
cmatinv(2,2)=1./e2
cmatinv(3,3)=1./e3
cmatinv(1,2)=-prxy/e1
cmatinv(2,1)=cmatinv(1,2)
cmatinv(1,3)=-prxz/e1
cmatinv(3,1)=cmatinv(1,3)
cmatinv(2,3)=-pryz/e2
cmatinv(3,2)=cmatinv(2,3)
cmatinv(4,4)=1/g23
cmatinv(5,5)=1/g13
cmatinv(6,6)=1/g12
write(*,*)’Compliance matrix (elasticity inverse) cinv= ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’(A different mapping from tensor to vector is used):’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
write(*,*)’The five independent compliance coeficients are:’
write(*,’(a,g15.8)’)’s11=’,cmatinv(1,1)
write(*,’(a,g15.8)’)’s33=’,cmatinv(3,3)
write(*,’(a,g15.8)’)’s44=’,cmatinv(4,4)
write(*,’(a,g15.8)’)’s12=’,cmatinv(1,2)
write(*,’(a,g15.8)’)’s13=’,cmatinv(1,3)
do i=1,6
do j=1,6
a(i,j)=cmatinv(i,j)
enddo
41enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
write(*,*)’Computing inverse matrix’
call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
endif
do i=1,6
do j=1,6
cmat(i,j)=b(i,j)
enddo
enddo
write(*,*)’ Elasticity matrix = ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmatinv,ia,m,n,b,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end of case1
c begin case2
if(reply .eq. ’2’)then
do i=1,6
do j=1,6
cmat(i,j)=0.
enddo
enddo
cmat(1,1)=13.2e10
write(*,’(a,g12.5,a)’)’Enter c11 [’,cmat(1,1),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,1)=ain(1)
endif
cmat(2,2)=cmat(1,1)
cmat(3,3)=11.5e10
write(*,’(a,g12.5,a)’)’Enter c33 [’,cmat(3,3),’]’
42call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(3,3)=ain(1)
endif
cmat(4,4)=5.2e10
write(*,’(a,g12.5,a)’)’Enter c44 (ANSYS c66) [’,cmat(4,4),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(4,4)=ain(1)
endif
cmat(5,5)=cmat(4,4)
cmat(1,2)=7.1e10
write(*,’(a,g12.5,a)’)’Enter c12 [’,cmat(1,2),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,2)=ain(1)
endif
cmat(2,1)=cmat(1,2)
cmat(1,3)=7.3e10
write(*,’(a,g12.5,a)’)’Enter c13 [’,cmat(1,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,3)=ain(1)
endif
cmat(3,1)=cmat(1,3)
cmat(2,3)=cmat(1,3)
cmat(3,2)=cmat(2,3)
cmat(6,6)=2.*(cmat(1,1)-cmat(1,2))
write(*,*)’Elasticity matrix c=’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
do i=1,6
do j=1,6
a(i,j)=cmat(i,j)
enddo
enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
43call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
endif
do i=1,6
do j=1,6
cmatinv(i,j)=b(i,j)
enddo
enddo
write(*,*)’Compliance matrix= ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmat,ia,m,n,cmatinv,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end case2
c begin case3
if(reply .eq. ’3’)then
do i=1,6
do j=1,6
cmatinv(i,j)=0.
enddo
enddo
cmatinv(1,1) = 1.65e-11
cmatinv(3,3) = 2.07e-11
cmatinv(4,4) = 4.35e-11
cmatinv(1,2) = -4.78e-12
cmatinv(1,3) = -8.45e-12
write(*,’(a,g12.5,a)’)’Enter s11 [’,cmatinv(1,1),’]’
write(*,*)’Reference for default values:’
write(*,*)’Auld, Appendix A.4, PZT-5H’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,1)=ain(1)
endif
cmatinv(2,2)=cmatinv(1,1)
44write(*,’(a,g12.5,a)’)’Enter s33 [’,cmatinv(3,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(3,3)=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter s44 (ANSYS c66) [’,cmatinv(4,4),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(4,4)=ain(1)
endif
cmatinv(5,5)=cmatinv(4,4)
write(*,’(a,g12.5,a)’)’Enter s12 [’,cmatinv(1,2),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,2)=ain(1)
endif
cmatinv(2,1)=cmatinv(1,2)
write(*,’(a,g12.5,a)’)’Enter s13 [’,cmatinv(1,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,3)=ain(1)
endif
cmatinv(3,1)=cmatinv(1,3)
cmatinv(2,3)=cmatinv(1,3)
cmatinv(3,2)=cmatinv(2,3)
cmatinv(6,6)=2.*(cmatinv(1,1)-cmatinv(1,2))
write(*,*)’Compliance matrix s= cinv =’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
do i=1,6
do j=1,6
a(i,j)=cmatinv(i,j)
enddo
enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
write(*,*)’Computing inverse matrix’
call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
45endif
do i=1,6
do j=1,6
cmat(i,j)=b(i,j)
enddo
enddo
write(*,*)’ Elasticity matrix = ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmatinv,ia,m,n,b,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end case3
write(*,*)’Orthotropic constants from compliance matrix:’
e1=1./cmatinv(1,1)
e2=1./cmatinv(2,2)
e3=1./cmatinv(3,3)
write(*,’(a,g15.8)’)’Young’’s Modulus 1= ’,e1
write(*,’(a,g15.8)’)’Young’’s Modulus 2= ’,e2
write(*,’(a,g15.8)’)’Young’’s Modulus 3= ’,e3
prxy=-cmatinv(1,2)*e1
write(*,’(a,g15.8)’)’Major Poisson ratio prxy=’,prxy
prxz=-cmatinv(1,3)*e1
write(*,’(a,g15.8)’)’Major Poisson ratio prxz=’,prxz
pryz=-cmatinv(2,3)*e2
write(*,’(a,g15.8)’)’Major Poisson ratio pryz=’,pryz
g23= 1./cmatinv(4,4)
g13= 1./cmatinv(5,5)
g12= 1./cmatinv(6,6)
write(*,’(a,g15.8)’)’Shear Modulus g23= ’,g23
write(*,’(a,g15.8)’)’Shear Modulus g13= ’,g13
write(*,’(a,g15.8)’)’Shear Modulus g12= ’,g12
write(*,*)
c orthotropic elasticity tensor
do i=1,3
do j=1,3
do k=1,3
e(i,j,k)=0.
do l=1,3
c(i,j,k,l)=0.
enddo
enddo
enddo
46enddo
do i=1,3
do j=1,3
call t2vec(i,j,ij)
do k=1,3
do l=1,3
call t2vec(k,l,kl)
if(k .eq. l)then
c(i,j,k,l)=cmat(ij,kl)
else
c(i,j,k,l)=cmat(ij,kl)/2.
endif
enddo
enddo
enddo
enddo
c
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 81 elasticity tensor components of c? [n]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’n’
endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
do l=1,3
if(c(i,j,k,l) .ne. zero)then
write(*,’(a,i1,i1,i1,i1,a,g15.8)’)’c(’,i,j,k,l,’)=’,c(i,j,k,l)
endif
enddo
enddo
enddo
enddo
endif
write(*,*)’ piezoelectric (stress) tensor coefficients’
do i=1,3
do j=1,3
do k=1,3
e(i,j,k)=0.
do l=1,3
do m=1,3
e(i,j,k) = e(i,j,k) + d(i,l,m)*c(l,m,j,k)
enddo
enddo
c write(*,’(a,i1,i1,i1,a,g15.8)’)’e(’,i,j,k,’)=’,e(i,j,k)
enddo
enddo
enddo
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 27 piezoelectric tensor components of e? [n]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’n’
47endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
if(e(i,j,k) .ne. zero)then
write(*,’(a,i1,i1,i1,a,g15.8)’)’e(’,i,j,k,’)=’,e(i,j,k)
endif
enddo
enddo
enddo
endif
c compute the e piezoelectric matrix
do i=1,3
do jk=1,6
call vec2t(jk,j,k)
if(j .eq. k)then
emat(jk,i) = e(i,j,k)
else
emat(jk,i) = 2.*e(i,j,k)
endif
enddo
enddo
c
write(*,*)’Piezoelectric e matrix= ’
do i=1,6
write(*,’(3(g11.4,1x))’)(emat(i,j),j=1,3)
enddo
write(*,*)’Second calculation:’
e13=cmat(1,1)*dmat(3,1)+cmat(1,2)*dmat(3,1)+cmat(1,3)*dmat(3,3)
write(*,*)’e13=’,e13
e33=cmat(1,3)*dmat(3,1)+cmat(1,3)*dmat(3,1)+cmat(3,3)*dmat(3,3)
write(*,*)’e33=’,e33
e42=cmat(4,4)*dmat(1,5)
write(*,*)’e42=’,e42
write(*,*)’Note: ANSYS e rows are permuted’
write(*,*)’ANSYS row 4 is row 6 ’
write(*,*)’ANSYS row 5 is row 4 ’
write(*,*)’ANSYS row 6 is row 5 ’
write(*,*)’ANSYS Piezoelectric e matrix= ’
do i=1,3
write(*,’(3(g11.4,1x))’)(emat(i,j),j=1,3)
enddo
write(*,’(3(g11.4,1x))’)(emat(6,j),j=1,3)
write(*,’(3(g11.4,1x))’)(emat(4,j),j=1,3)
write(*,’(3(g11.4,1x))’)(emat(5,j),j=1,3)
c dielectric constants (farads/meter)
48c dc(1,1)=7.12e-9
c dc(2,2)=7.12e-9
c dc(3,3)=5.84e-9
c density (kilograms/meter^3)
c density=7.5e3
cm=’,’
idmat=3
m=3
n=6
iemat=6
l=3
icc=6
call matm(dmat,idmat,m,n,emat,iemat,l,cc,icc)
write(*,*)’d times e = ’
do i=1,3
write(*,’(3(g11.4,1x))’)(cc(i,j),j=1,3)
enddo
c compute permittivity tensor
eps0=8.85e-12
do i=1,3
do j=1,3
epsilon(i,j)=0.
enddo
enddo
dielx=1300.
write(*,’(a)’)’Enter the unclamped dielectric constant K11’
write(*,’(a,g15.8,a)’)’[’,dielx,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
dielx=ain(1)
endif
epsilon(1,1)=dielx*eps0
diely=dielx
c cs1=’Enter the dielectric constant K22’
c write(*,’(a,a,g15.8,a)’)cs1,’ [’,diely,’]’
c call readr(0,ain,nr)
c if(nr .eq. 1)then
c diely=ain(1)
c endif
epsilon(2,2)=diely*eps0
dielz=1000.
write(*,’(a)’)’Enter the unclamped dielectric constant K33’
write(*,’(a,g15.8,a)’)’[’,dielz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
dielz=ain(1)
endif
epsilon(3,3)=dielz*eps0
write(*,*)’Unclamped electric permittivity tensor=’
do i=1,3
49write(*,’(3(g11.4,1x))’)(epsilon(i,j),j=1,3)
enddo
do i=1,3
do j=1,3
clamped(i,j)= epsilon(i,j) - cc(i,j)
enddo
enddo
write(*,*)’Clamped electric permittivity tensor=’
do i=1,3
write(*,’(3(g11.4,1x))’)(clamped(i,j),j=1,3)
enddo
do i=1,3
do j=1,3
sum=0.
do k=1,3
do l=1,3
sum=sum+d(i,k,l)*e(j,k,l)
c if((i .eq. 1) .and. (j .eq. 1))then
c write(*,*)k,l,’ d= ’,d(i,k,l),’ e= ’,e(j,k,l)
c endif
enddo
enddo
epsilon(i,j)=epsilon(i,j)-sum
enddo
enddo
c write(*,*)’Clamped electric permittivity tensor=’
c do i=1,3
c write(*,’(3(g11.4,1x))’)(epsilon(i,j),j=1,3)
c enddo
rho= 7.500000E-03
write(*,*)’Enter density Kg/m^3 [7500]’
call readr(0,ain,nr)
if(nr .ge. 1)then
rho=ain(1)
endif
write(*,’(a)’)’/COM ’
write(*,’(a)’)’/COM Material properties’
c material number
mn=2
write(*,’(a)’)’/COM Young’’s Moduli’
write(*,’(a,i3,a,g15.8)’)’EX,’,mn,cm,e1
write(*,’(a,i3,a,g15.8)’)’EY,’,mn,cm,e2
write(*,’(a,i3,a,g15.8)’)’EZ,’,mn,cm,e3
write(*,’(a)’)’/COM ANSYS "Major" Poisson ratios’
c reference ANSYS Theory 2.1 structural fundamentals
prxy=-cmatinv(1,2)*e1
write(*,’(a,i3,a,g15.8)’)’PRXY,’,mn,cm,prxy
50prxz=-cmatinv(1,3)*e1
write(*,’(a,i3,a,g15.8)’)’PRXZ,’,mn,cm,prxz
pryz=-cmatinv(2,3)*e2
write(*,’(a,i3,a,g15.8)’)’PRYZ,’,mn,cm,pryz
write(*,’(a)’)’/COM Density Kg/m^3’
write(*,’(a,i3,a,g15.8)’)’/COM DENS,’,mn,cm,rho
write(*,’(a)’)’/COM clamped Permitivities’
c epsilon(1,1)=8.7969e-12
c epsilon(2,2)=8.7969e-12
c epsilon(3,3)=8.7969e-12
write(*,’(a,i3,a,g15.8)’)’MP,PERX,’,mn,cm,clamped(1,1)
write(*,’(a,i3,a,g15.8)’)’MP,PERY,’,mn,cm,clamped(2,2)
write(*,’(a,i3,a,g15.8)’)’MP,PERZ,’,mn,cm,clamped(3,3)
write(*,’(a)’)’/COM ANSYS Piezoelectric "e" matrix’
write(*,’(a)’)’TB,PIEZ,3’
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,3,cm,emat(1,3)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,6,cm,emat(2,3)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,9,cm,emat(3,3)
c ansysemat(5,2)=emat(4,2)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,14,cm,emat(4,2)
c ansysemat(6,1)=emat(5,1)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,16,cm,emat(5,1)
end
c+ readr read a row of floating point numbers
subroutine readr(nf, a, nr)
implicit real*8(a-h,o-z)
c numbers are separated by spaces
c examples of valid numbers are:
c 12.13 34 45e4 4.78e-6 4e2,5.6D-23,10000.d015
c nf=file number, 0 for standard input file
c a=array of returned numbers
c nr=number of values in returned array,
c or 0 for empty or blank line,
c or -1 for end of file on unit nf.
c requires functions val and length
dimension a(*)
character*200 b
character*200 c
character*1 d
c=’ ’
if(nf.eq.0)then
read(*,’(a)’,end=99)b
else
read(nf,’(a)’,end=99)b
endif
nr=0
l=lenstr(b)
if(l.ge.200)then
write(*,*)’ error in readr subroutine ’
write(*,*)’ record is too long ’
endif
do 1 i=1,l
d=b(i:i)
if (d.ne.’ ’) then
51k=lenstr(c)
if (k.gt.0)then
c=c(1:k)//d
else
c=d
endif
endif
if( (d.eq.’ ’).or.(i.eq.l)) then
if (c.ne.’ ’) then
nr=nr+1
call valsub(c,a(nr),ier)
c=’ ’
endif
endif
1 continue
return
99 nr=-1
return
end
c
c+ lenstr nonblank length of string
function lenstr(s)
c length of the substring of s obtained by deleting all
c trailing blanks from s. thus the length of a string
c containing only blanks will be 0.
character s*(*)
lenstr=0
n=len(s)
do 10 i=n,1,-1
if(s(i:i) .ne. ’ ’)then
lenstr=i
return
endif
10 continue
return
end
c+ valsub converts string to floating point number (r*8)
subroutine valsub(s,v,ier)
implicit real*8(a-h,o-z)
c examples of valid strings are: 12.13 34 45e4 4.78e-6 4E2
c the string is checked for valid characters,
c but the string can still be invalid.
c s-string
c v-returned value
c ier- 0 normal
c 1 if invalid character found, v returned 0
c
logical p
character s*(*),c*50,t*50,ch*15
character z*1
data ch/’1234567890+-.eE’/
v=0.
ier=1
l=lenstr(s)
if(l.eq.0)return
52p=.true.
do 10 i=1,l
z=s(i:i)
if((z.eq.’D’).or.(z.eq.’d’))then
s(i:i)=’e’
endif
p=p.and.(index(ch,s(i:i)).ne.0)
10 continue
if(.not.p)return
n=index(s,’.’)
if(n.eq.0)then
n=index(s,’e’)
if(n.eq.0)n=index(s,’E’)
if(n.eq.0)n=index(s,’d’)
if(n.eq.0)n=index(s,’D’)
if(n.eq.0)then
s=s(1:l)//’.’
else
t=s(n:l)
s=s(1:(n-1))//’.’//t
endif
l=l+1
endif
write(c,’(a30)’)s(1:l)
read(c,’(g30.23)’)v
ier=0
return
end
c+ str floating point number to string
subroutine str(x,s)
implicit real*8(a-h,o-z)
character s*25,c*25,b*25,e*25
zero=0.
if(x.eq.zero)then
s=’0’
return
endif
write(c,’(g11.4)’)x
read(c,’(a25)’)b
l=lenstr(b)
do 10 i=1,l
n1=i
if(b(i:i).ne.’ ’)go to 20
10 continue
20 continue
if(b(n1:n1).eq.’0’)n1=n1+1
b=b(n1:l)
l=l+1-n1
k=index(b,’E’)
if(k.gt.0)e=b(k:l)
if(k.gt.0)then
s=b(1:(k-1))
k1=index(b,’E+0’)
if(k1.gt.0)then
e=’E’//b((k1+3):l)
53else
k1=index(b,’E+’)
if(k1.gt.0)e=’E’//b((k1+2):l)
endif
k1=index(b,’E-0’)
if(k1.gt.0)e=’E-’//b((k1+3):l)
l=k-1
else
s=b
endif
j=index(s,’.’)
n2=l
if(j.ne.0)then
do 30 i=1,l
n2=l+1-i
if(s(n2:n2).ne.’0’)go to 40
30 continue
endif
40 continue
s=s(1:n2)
if(s(n2:n2).eq.’.’)then
s=s(1:(n2-1))
n2=n2-1
endif
if(k.gt.0)s=s(1:n2)//e
return
end
c+ t2vec stress tensor index to vector index
subroutine t2vec(i,j,k)
c input:
c i,j
c output:
c k
c so that sigma(k)=sigma(i,j)
integer m(3,3)
data m/1,6,5,6,2,4,5,4,3/
k=m(i,j)
return
end
c+ vec2t stress vector index to tensor index
subroutine vec2t(k,i,j)
c input:
c k
c output:
c i,j
c so that sigma(i,j)=sigma(k)
integer mi(6),mj(6)
data mi/1,2,3,2,1,1/
data mj/1,2,3,3,3,2/
i=mi(k)
j=mj(k)
return
end
c+ gaussr solution of linear equations, inverse, determinant (real) new version
c gaussr2.ftn 4/17/96 modernization of gaussr under construction
54subroutine gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
c solves the equation a*c=b for c, where a is an n by n matrix
c c and b are n row by m column matrices. c is returned as b.
c algorithm -gaussian elimination with partial pivoting.
c parameters a-n by n matrix containing the coefficients of
c the linear system.
c ia-row dimension of a in defining routine
c e.g. in the routine where a is defined,
c a might be dimensioned as:
c dimension a(nr,nc)
c then ia must be set to nr. we may have n < ia, but
c must not have n > ia, or n*n > nr*nc.
c ia is needed for proper addressing of matrix a.
c fortran stores by column first: a(i,j)=a(i+(j-1)*ia))
c b-n by m matrix containing the m right sides
c of the equations, on entering. on returning, b
c contains the solutions. the inverse of a is
c returned in b when inv=1
c ib-row dimension of b in defining routine
c n-row and column dimension of a.
c m-column dimension of b (usually 1)
c the program changes m to n when inv=1
c inv-the inverse of a is calculated
c and returned in b when inv=1
c b must be large enough to hold the inverse
c eps-each equation is normalized so that the
c coefficients are <= 1 in magnitude.
c when a pivot is less than eps the matrix is
c considered nearly singular, and ier is set to 1
c if a pivot is zero the matrix is singular, and
c ier is set to 2. one may set eps=1.e-5 for
c single precision, and 1.e-12 for double.
c eps does not effect any calculation.
c normalization may also prevent exponent overflow.
c idet-compute determinant only if idet = 1
c determinants are products of n numbers.
c overflow can occur if the elements of the
c matrix have large exponents.
c set idet=0 if the determinant is not needed.
c det-determinant of a.
c ier-return parameter,
c ier=0 normal return
c ier=1 matrix is nearly singular
c ier=2 matrix is singular
c
c warning!! the subroutine changes a and b. if they need to be
c saved, copies must be made before calling the subroutine.
c the subroutine can be converted to different number type
c by uncommenting the appropriate implicit statement.
implicit real*8(a-h,o-z)
c implicit complex(a-h,o-z)
c implicit complex*16(a-h,o-z)
logical cdet
dimension a(ia,*),b(ib,*)
cdet = idet .eq. 1
55zero=0.
ier=0
det=1.
if(m.le.0)m=1
if(inv.eq.1)then
c set b equal to the identity
do i=1,n
do j=1,n
b(i,j)=0.
if(i.eq.j)b(i,j)=1.
enddo
enddo
m=n
endif
c normalize rows
do 20 i=1,n
bigest=a(i,1)
do 16 j=2,n
ab=a(i,j)
if(abs(ab).gt.abs(bigest))bigest=ab
16 continue
if(bigest.eq.zero)then
ier=2
det=0.
return
endif
if(cdet)det=det*bigest
do 18 j=1,n
a(i,j)=a(i,j)/bigest
18 continue
do 19 j=1,m
b(i,j)=b(i,j)/bigest
19 continue
20 continue
j=1
do while( j .lt. n)
kk=j+1
l=j
c find row l with largest pivot
do 32 i=kk,n
if(abs(a(i,j)).gt.abs(a(l,j)))l=i
32 continue
if(abs(a(l,j)).eq.zero)then
ier=2
det=0.
return
endif
if(abs(a(l,j)).le.abs(eps))ier=1
if(l.ne.j)then
c interchange rows l and j
do 37 k=1,n
c=a(l,k)
a(l,k)=a(j,k)
a(j,k)=c
37 continue
56do 39 k=1,m
c=b(l,k)
b(l,k)=b(j,k)
b(j,k)=c
39 continue
if(cdet)det=det*(-1.)
endif
if(cdet)det=det*a(j,j)
c divide row by pivot
c=a(j,j)
do 50 k=j,n
a(j,k)=a(j,k)/c
50 continue
do 55 k=1,m
b(j,k)=b(j,k)/c
55 continue
c add multiple of row j to lower rows
c to eliminate jth coefficients
jj=j+1
do 80 i=jj,n
am=a(i,j)
do 60 k=1,n
a(i,k)=a(i,k)-am*a(j,k)
60 continue
do 70 k=1,m
b(i,k)=b(i,k)-am*b(j,k)
70 continue
80 continue
j=j+1
enddo
am=a(n,n)
if(abs(am).eq.zero)then
ier=2
det=0.
return
endif
if(abs(am).le.abs(eps))ier=1
if(cdet)det=det*am
c a is now in triangular form
c compute nth component of solution
do 90 k=1,m
b(n,k)=b(n,k)/am
90 continue
c back substitute to compute n-i component
c i=1,2,3,...
nn=n-1
do 120 i=1,nn
ni=n-i
do 110 j=1,m
nj=ni+1
do 100 ki=nj,n
b(ni,j)=b(ni,j)-a(ni,ki)*b(ki,j)
100 continue
110 continue
120 continue
57return
end
subroutine matm(a,ia,m,n,b,ib,l,c,ic)
implicit real*8(a-h,o-z)
c arguments
c a-matrix
c ia-row dimension of a in calling program
c m-number of rows of a
c n-number of columns of a
c b-matrix
c ib-row dimension of b in calling program
c l-number of columns of b
c c-product matrix: c=a*b
c ic-row dimension of c in calling program
c
dimension a(ia,*),b(ib,*),c(ic,*)
c c=a*b
do 10 i=1,m
do 10 j=1,l
c(i,j)=0.
do 5 k=1,n
5 c(i,j)=c(i,j)+a(i,k)*b(k,j)
10 continue
return
end
0.22 Location of Piezoelectric Information In
ANSYS Manuals
• Index assignment for stress-strain tensor to six-vector: Theory 2.1.
• Major and Minor Poisson ratios, orthotropic compliance matrices: Theory 2.1.
• Piezoelectric Solid Element SOLID5: Elements manual 4-39.
• Piezoelectric Plane Element PLANE13: Elements manual 4-85.
• Piezoelectric Tetrahedral Coupled Solid Element SOLID98: Elements
manual 4-653.
• Piezoelectric Analysis: Procedures Coupled-Field Analysis 8.1, piezoelectric analysis p8-13 to p8-15.
• Deﬁning matrices: Commands TB and TBDATA, Commands manual.
58• Piezoelectrics Theory: Theory Manual 11.1 to page 11-17.
• Example VM175: Natural Frequency of a Piezoelectric Transducer
...(ANSYS Stuﬀ Notebook). Reference: Boucher D, Lagier M, Maerfeld C, IEEE Transactions on Sonics and Ultrasonics, Volume SU-28,
No. 5 September 1981 pp318-330.
• Example VM176: Frequency Response of Electrical Admittance for a
Piezoelectric Transducer, Reference: Kagawa and Yamabuchi, IEEE
Trans. Sonics and Ultrasonics, V. SU-26, No. 2, March 1979.
• Example PM13, ANSYS Examples Supplement, Piezoelectric Beam
Resonator
0.23 ANSYS Comparison For A Composite
Piezoelectric Transducer: Example VM176.
Reference for the problem: Kagawa and Yamabuchi Finite Element Simulation of a Composite Piezoelectric Transducer IEEE Transactions
on Sonics and Ultrasonics, Vol. SU-26, No 2 march 1979.
The transducer consists of cylindrical disks of aluminum, adhesive, PZT,
adhesive, and aluminum.
PZT material: NEPEC 6
Elasticity matrix
c = 10
10








12.8 6.8 6.6 0  0  0
6.8 12.8 6.6 0  0  0
6.6 6.6 11  0  0  0
0 0 0 2.1 0  0
0 0 0 0 2.1 0
0  0  0 0 0 2.1








59Piezoelectric matrix (C oulomb/M  eter
2
)
e =






0 0 −6.1
0 0 −6.1
0 0  15.7
0 0  0
0 0  0
0 0  0






Dielectric Matrix (F arad/M eter)
 = 10−9



8.7969 0 0
0 8.7969 0
0 0 8.7969



See ANSYS notebook for more details.
600.24 Piezoelectricity Bibliography
[1]Agarwal Bhagwan D, Broutman Lawrence, Analysis And Performance
Of Fiber Composites.
[2]Ahmed F R, Crystallographic Computing, Proceedings International
Summer School, Linda Hall Library.
[3]Allik H, Hughes J R, Finite Element Method For Piezoelectric Vibration, International Journal of Numerical Methods For Engineering, V.
2, pp 151-157, 1970, (This is the formulation of the problem that is implemented in ANSYS.).
[4]Ashton J E, Whitney J M, Theory Of Laminated Plates, Linda Hall
Library, 1979.
[5]Ashton T E, Halpin T C, Petit P H, Primer On Composite Materials
Analysis, Technomic, Stamford Conn., 1969.
[6]Auld B A, Acoustic Fields And Waves In Solids, Volumes I and II.
[7]Billington E W, Introduction To The Mechanics And Physics Of
Solids, AlliedSignal Library, QA808 .2, 1986, A Mathematical Book. Von
Mises Yield Criterion. Von Mises Stress.
[8]Bottom Virgil E, The Theory And Design Of Quartz Crystal Units,
2nd Ed., 1974, (Self Published), Linda Hall, (Piezoelectric Transducers.).
[9]Boucher D, Lagier M, Maerfeld, Computation Of The Vibrational
Modes For Piezoelectric Array Transducers Using A Mixed Finite
Element-Perturbation Method, IEEE Transactions on Sonics And Ultrasonics, V. Su28, No. 8, Sept 81.
[10]Brebbia, Dewilde,Blain, Computer Aided Design In Composite Material Technology, TA 418.9 .c6 I585, 1988, Linda Hall.
[11]Buchanon Relva C, Editor, Ceramic Materials For Electronics, 2nd
Ed., 1991, AlliedSignal Library, (Piezoelectricity.).
[12]Buerger, Martin, Introduction To Crystal Geometry.
[13]Cady Walter Guyton, Piezoelectricity, Volumes I and II, Dover, 1964,
61KCMO Library.
[14]Cady Walter Guyton, The Piezoelectric Resonator, In: Physical
Acoustics, R. Bruce Lindsey Editor, (Benchmark Papers In Acoustics).
[15]Clough R W, Penzien J, Dynamics Of Structures, McGraw-Hill 1972,
AlliedSignal Library, (Rayleigh Damping, Piezoelectricity.).
[16]Cracknell A P, Ultrasonics, Wykeham Publications, 1980 AlliedSignal
Library.
[17]Crandall B C, Lewis James, Editors, Nanotechnology, Article: John
Foster, p17, Scanning Tunneling Microscope (STM), (Uses Piezoelectric Control,), 1992, MIT Press.
[18]Cullity, Elements Of X-Ray Diﬀraction, Addison-Wesley, 1956.
[19]Cummings James R, Predicting The Eﬀect Of Thickness Of Piezoelectric Plates On The Performance Of The Ultrasonic Traveling
Wave Motor, U. Mo. (Rolla,), May 1993.
[20]Cummings James R, Stutts Daniel S, Dharani L R, Eﬀects Of Material And Geometric Characteristics On The Performance Of A
Traveling-wave Motor, U. MO. Rolla Report, 11-9-93.
[21]Cummings Robert James, Analysis Of The Performance And Contact Behavior Of The Piezoelectric Traveling-wave Motor, Masters
Theses, University Of Missouri-Rolla MRTC Manufacturing Research And
Training Center, 1994.
[22]Emery James D, Elasticity, Unpublished Document.
[23]Emery James D, Electomagnetic Theory, Unpublished Document.
[24]Emery James D, Transducer Equivalent Circuits, Unpublished Document.
[25]Emery James D, Piezo Notes I, Unpublished Document.
[26]Emery James D, Piezo Notes II, Unpublished Document.
[27]Flynn Anita, Torque Production In Ultrasonic Motors, MIT Arti-
62ﬁcial Intelligence Laboratory, January 9, 1994, (Piezoelectric Motor.).
[28]Flynn Anita M, Piezoelectric Ultrasonic Micromotors, PhD Thesis,
MIT, June 1995.
[29]Fulton Langdon H, Techniques And Applications Of The Gas Squeeze
Film Bearing, Masters Thesis, Moore School, U. Pennsylvania, Dec 1967,
(Piezoelectric Application. Now At Martin-Marietta Communications Systems, Camden NJ).
[30]Gooberman G L, Ultrasonics, 1968, Acoustic Impedance, Transmission,
Lines, Transducers, Equivilent Circuit As Transmission Line, Q Of Transducer.
[31]Goree, Pagano-Goree Free Edge Stress, (Discussion, Goree: Singular
Elements In Abaqus Are Needed, Disagreement With Whitney.).
[32]Hagedorn P, Wallashek J, Traveling Wave Ultrasonic Motors, Part
I: Working Principle And The Mathematical Modeling Of The Stator, J. Of Sound And Vibration, V. 155, No. 1, pp 31-46, 1992.
[33]Hagood N, Chung W,Von Flotow A, Modeling Of Piezoelectric Actuator Dynamics For Active Structural Control, J. Intell.  Mat.  Syst.
And Structures, V. 1, pp327-354, July 1990, (Hamilton’s Principal Modiﬁed
For General Electromechanical Systems.).
[34]Hagood Nesbit, Mcfarland Andrew, Modeling Of A Piezoelectric Rotary Ultrasonic Motor, IEEE Trans. Ultrasonics Ferroelectrics And Frequency Control, V. 42, No. 2, March 1995, Nesbit was Anita Flynn Collegue
at MIT. Rayleigh-Ritz Method, which includes, Electrical Energy.
[35]Hauptman, Herbert, Crystal Structure Determination, Linda Hall,
QD911.h33, 1972, Mathematical Biophysics Department, Medical Foundation Of Buﬀalo, Cosine Semiinvariants.
[36]Heising Raymond A, Quartz Crystals For Electrical Circuits: Their
Design And Manufacture, Van Nostrand, 1946, Reprinted by Electronic
Industries Association, 2001 Eye St, NW, Washington DC 20006, 1978, 1982,
(Piezoelectric Equivalent Circuits.).
[37]Hirose, Takahashi, Uchino, Aoyagi, Tomakawa, Electromechanical Prop-
63erties Of PbZro3-pbzrtio3-pb(mn(1/3)sb(2/3)o3 Ceramics Under
Vibration Level Change, Mat. Res. Soc. Symp. Proc., V. 360, 1995,
Materials Research Society.
[38]Hirose, Takahashi, Uchino, Aoyagi, Tomakawa, Measuring Methods
For High-Power Characteristics Of Piezoelectric Materials, Mat.
Res. Soc. Symp. Proc., V. 360, 1995, Materials Research Society.
[39]Hossack John A, Hayward Gordon, Finite-Element Analysis of 1-3
Composite Transducers, IEEE Transactions on Ultrasonics, Ferroelectrics,
and Frequency Control, V. 38, N0. 6, November 1991, (Uses ANSYS.).
[40]Inman Daniel, Engineering Vibration, MatLab Disk, 1994, Chapter:
Experimental Modal Analysis.
[41]Institute of Radio Engineers, IRE Standards on Piezoelectric Crystals: Determination of the Elastic, Piezoelectric, and Dielectric
Constants, The Electromechanical Coupling Factor, 1958, Proc. IRE,
V. 46, pp764-78.
[42]Jaﬀe Bernard, Cook William, Jaﬀe Hans, Piezoelectric Ceramics, Academic Press, 1971, Linda Hall, TK7871.15.C4 J3, Deﬁnitions, Symmetry,
Constants, Measurements, Hysteresis.
[43]Johnson K L, Contact Mechanics, Linda Hall Library.
[44]Jones, Acoustic And Electromagnetic Waves, Oxford.
[45]Kagawa, Yamabuchi, Finite Element Simulation of a Composite
Piezoelectric Transducer, IEEE Transactions on Sonics and, Ultrasonics,
V. SU-26, No. 2 march 1979, (ANSYS example VM176.).
[46]Kikuchi And Oden, Contact Problems In Elasticity: A Study Of
Variational Inequalities And Finite Element Methods, Linda Hall,
QA931.k47 1988, (SIAM), Society For Industrial And Applied Mathematics.
[47]Kittle, Solid State Physics, (Piezoelectric Equations.).
[48]Koehler D R, Radial Vibration Modes In Thin Plates Of BaT iO3
And LiNBO3, Memo, Sandia National Laboratory, October 16, 1989.
64[49]Koiter W T , Miknailor G K, Theory Of Shells.
[50]Kolsky, Stress Waves In Solids, Dover.
[51]Kompaneyts A S, Theoretical Physics, Dover, 1961.
[52]Lass Victor, Vector And Tensor Analysis, (Relation Between Mechanical, Electric, And Magnetic Energy.).
[53]Lazan B J, Damping Of Materials And Members In Structural
Mechanics, Oxford 1968.
[54]Lerch, R, Simulation Of Piezoelectric Devices By Two And Three
Dimensional Finite Elements, IEEE Trans. On Ultrasonics Ferroelectrics
And Freq. Control, V. 37, No. 2, May 1990.
[55]Levich B G, Theoretical Physics, Volumes I, Ii, Iii, Wiley-Interscience
1971.
[56]Lur’e A L, Three-dimensional Problems Of Elasticity, (Basis Of
James Cummings Rolla Analysis Of Tooth Contact In Piezo Motor).
[57]Mansﬁeld, E H, The Bending And Stretching Of Plates, (The Small
Deﬂection Theory Of Plates Is Analogous To Beam Small, Deﬂection Theory. Gaussian And Mean Curvature, Where Curvature Is Taken To, Be
Approximated By The 2nd Order Derivatives. Twist Is The Cross Derivative, Transformation Matrix For Curvatures Handled As In JRC (above).
[58]Marsden Jerrold, The Mathematical Theory Of Elasticity, (Contains Much Diﬀerential Geometry).
[59]Martin, Transformation Geometry, (Symmetry Groups, Transformations).
[60]Mason Warren P, Piezoelectric Crystals and Their Application to
Ultrasonics, Van Nostrand, 1950.
[61]Mason Warren P, Physical Acoustics, 1964 Academic, AlliedSignal Library, QC225.m3, (Piezoelectric Actuators.).
[62]Mcfarland, Smith, Bernhart, Analysis Of Plates, Linda Hall Library.
65[63]McLachlan N W, Theory Of Vibrations, Circular Membrane, Circular
Plate, Bessel Functions p 129-144.
[64]Meirovitch Leonard, Elements Of Vibration Anlysis, Mcgraw-Hill,
1975.
[65]Meirovitch, Leonard, Analytical Methods Of Vibration, Vibration
Of Plates p 179, Eigenvalues, Rayleigh Quotient, The Biharmonic Equation.
[66]Mentesana C P, Optical Stronglink Development And Piezoelectric Motor Evaluation, AlliedSignal KCD, Kcp-613-4573, May 1992.
[67]Mercer C D, Reddy B D, Eve R A, Finite Element Method For Piezoelectric Media, UCT/CSIR Applied Mechanics Research Unit, Technical
Report No. 92, April, 1987.
[68]Merhaut Josef, Theory Of Electroacoustics, 1976, Czech Book, Johnson County Library, (Piezoelectric Transducers.).
[69]Mikhlin, S.g, Variational Methods In Mathematical Physics, Linda
Hall Library, 1964.
[70]Mizushima Masataka, Theoretical Physics, John Wiley, 1972, AlliedSignal, (From Classical Mechanics To General Relativity, Relativistic Quantum Mechanics And Group Theory.).
[71]Moore And Christopher Davis And Coplan, Building Scientiﬁc Apparatus, (Experimental Techniques.).
[72]Morgan Matroc Incorporated, Guide To Modern Piezoelectric Ceramics, Catalog, Vernitron Division.
[73]Morse, Philip M, Vibration And Sound.
[74]Musikant Solomon, What Every Engineer Should Know About Ceramics, 1991, Marcel Dekker, AlliedSignal Book.
[75]Newnham R E, Trolier-McKinstry S, Giniewicz J R, Piezoelectric, Pyroelectric And Ferroic Crystals, J. Material Education, 15, pp 189-223
(1993), (Reprint: Penn State MRL Annual Report To ONR 1994, Appendix
661).
[76]Nye J F, Physical Properties Of Crystals, QD931 .N8, 1957, 1987.
[77]Ostergard D, Pawlak T, Three-Dimensional Finite Elements For
Analyzing Piezoelectric Structures, Proceedings IEEE Ultrasonics Symp.,
Williamsburg, Va, 1986, pp 639-642.
[78]Paganno N J, Pipes R B, Interlaminar Stresses In Composite Laminates Under Uniform Axial Extension, Journal Of Composite Materials, V. 4, 1970, pp 538-548.
[79]Paganno N J, Pipes R B, The Inﬂuence Of Stacking Sequence On
Laminate Strength, Journal Of Composite Materials, V. 5, 1971, pp 50-57.
[80]Panasonic, Ultrasonic Motor, Panasonic Electric Motor Division, Matsushita Electric Industrial Co. Ltd., (Equivalent Circuit.).
[81]Parker Sybyl, Solid State Physics Sourcebook.
[82]Press, Flannery, Teukolsky, Numerical Recipes in Fortran, (Subroutines For Computing Bessel Functions.).
[83]Pressly Robert B, Mentesana Charles P, Piezoelectric Motor Development At AlliedSignal Kansas City Division, KCP-613-5515, November, 1994.
[84]Puppo A. H. Evenson H A, Journal Of Composite Materials, V. 4,
1970, pp 201-214.
[85]Raju P. Narayana, Vibration Of Annular Plates, Journal Of The
Aeronautical Society Of India, V. 14, No. 2, May, 1962, (Boundary Conditions For Maple Motor Code).
[86]Ramachandran G N , Srinivasan R, Fourier Methods In Crystallography, Linda Hall Library, 1970.
[87]Rayleigh Lord, The Theory Of Sound, Volumes I and II, Dover, (Vibrations Of Plates: Chapter 10, p 352, Origin of Rayleigh Damping Concept).
67[88]Roiguiskis, Bansevicius, Barauskas, Kulviets, Vibromotors For Precision Microrobots, AlliedSignal Library, (Many Kinds Of Piezoelectric
Motor Designs.).
[89]Saigoh Hiroaki, Kawasaki Mayumi, Maruko Nobuhiro, Kanayama Kouichi,
Multilayer Piezoelectric Motor Using The First Longitudinal And
Second Bending Vibrations, Japanese J Applied Physics, V. 34, 1995,
pp2760-2764.
[90]Sands Donald E, Introduction To Crystallography, 1975, Dover.
[91]Sashida Toshiiku, Kenjo Takashi, An Introduction To Ultrasonic
Motors, Oxford U. Press 1993, Piezoelectricity.
[92]Sawyer C B, Tower C H, Phys. Rev. 35, (1930) pp269-73, (Technique For Displaying Ferroelectric-Piezoelectric Hysteresis Curves.).
[93]Sherrit S, Gauther N, Wiederick H D, Mukherjee B K, Accurate Evaluation Of The Real And Imaginary Material Constants For A Piezoelectric Resonator In The Radial Mode, Ferroelectrics, 1991, V. 119,
pp17-32.
[94]Shields John Potter, Basic Piezoelectricity, 1966, Sams, QC595 .S4,
Linda Hall Library.
[95]Smits Jan G, Iterative Method For The Accurate Determination
Of The Real And Imaginary Parts Of The Materials Coeﬃcients
Of Piezoelectric Ceramics, Trans. Sonics Ultraacoustics. V. Su-23, No.
6, Nov. 1976, pp393-402.
[96]Soedel Werner, Vibrations Of Shells And Plates, 1981. Second Edition 1993, University Illinois At Chicago.
[97]Sokolnikoﬀ I S, Mathematical Theory Of Elasticity, 1956, (Biharmonic Equation p77, Two Dimensional Problems Ch. 5, Finite Diﬀerences,
Relaxation Methods, pp442-454.).
[98]Southwell R V, Theory Of Elasticity.
[99]Stoker, Nonlinear Vibrations.
68[100]Stutts Daniel, Loss Factors, FAX: Feb 23, 1996, Various Materials
Loss Coeﬁcients vs Young Modulus, Comparison Of Damping for a Variety
of Metals, Magnetoeleastic Damping In Ferromagnetic Metals, Internal Friction Vs Strain Amplitude, Experimental Setup With Pucot Drive Crystal
And Specimin.).
[101]Stutts Daniel, Modal Analysis Of The Free-Free Rod With Damping, FAX: April, 9 1996, (See Piezo96 Notebook.).
[102]Stutts Daniel, A Simple Example Of The Relationship Between
Electromechanical Coupling And Measured Impedance, Unpublished
Document, University of Missouri at Rolla, Mechanical Engineering Department, March 1, 1995.
[103]Stutts Daniel, The Free-free Annular Plate, Maple Program And
Document, U. Mo. Rolla Mech. Eng. Dept, January,1994.
[104]Stutts Daniel S, Friend James R (Formerly Cummings), Quarterly Report: Piezoelectric Motor Project, Oct. 2 1995 (tooth Motion, Asymmetry, Maple Conference, Contact Mech, Web Page).
[105]Takahashi Sadayuki, Herose Seiji, Uchino Kenji, Oh Ki-Young, Electromechanical Characteristics Of Lead-Zirconate-Titanate Ceramics
Under Vibration Level Change, 1995 IEEE.
[106]Takuro Ikeda, Fundamentals Of Piezoelectricity,, (Electromechanical Coupling Coeﬃcient, X-cut, Y-cut Plate, Equivalent Circuit. p84.).
[107]Termin, Radio Engineering, (Equivalent Circuit Of Piezo Oscillator).
[108]Tiersten H F, Linear Piezoelectric Plate Vibration, Plenum Press,
1969, pp33-99, (Fundamental Equations And Equivalent Circuit.).
[109]Timoshenko And Goodier, Theory Of Elasticity, Third Edition.
[110]Timoshenko, And Woinowsky-Krieger, Theory Of Plates And Shells,
2nd Ed.
[111]Tomikawa Y, Ueha S, Ultrasonic Motors: Theory And Application, Oxford 1993.
69[112]Tuinenga Paul, Spice, A Guide To Circuit Simulation And Analysis Using Pspice.
[113]Tzou H S, Piezoelectric Shells, Kluwer Academic Publishers, 1993.
[114]Urick Robert J, Principles Of Underwater Sound, 2nd Ed., McGrawHill, 1975, (Piezoelectric Transducers.).
[115]Van Randeraat J, Setterington R E, Piezoelectric Ceramics, London,
Mullard, Linda Hall Library, 1974.
[116]Voyiadjis, Karamanlidis, Editors, Advances In The Theory Of Plates
And Shells, Linda Hall Library.
[117]Wallashek Jorg, Piezoelectric Ultrasonic Motors, Conference On
Mechatronics And Robotics, Duisburg/Moers, Germany, Sept 27-29, 1993
pp107-125.
[118]Webb S, The Physics Of Medical Imaging, 1988, Ultrasound, NMR,
CT, Tomography.
[119]Weinberger, Hans, Partial Diﬀerential Equations, (Sturm-Liouville
Problems, pp182-187, Circular Membrane, Bessel Functions.).
[120]Weinberger, H F, Variational Methods In Boundary Value Problems, 1961.
[121]Wert, Thomson, Physics Of Solids.
[122]Whitney J M, Browning C E, Journal Of Composite Materials, V.
6, April 1972, pp 50-57.
[123]Whitney J M, Stansbarger D L, Howell H B, Journal Of Composite
Materials, V. 5, April 1971, pp24-34.
[124]Whitney James M, Structural Analysis Of Laminated Isotropic
Plates, 1987, Technomic, (Includes Lampcal Software.).
[125]Whitney, Pipes, Experimental Mechanics Of Fiber Reinforced
Composite Materials, 1982, T-C Publishers.
[126]Willmore, T J, Tensors As Tensor Products, Linda Hall Library,
701959.
[127]Woolett Ralph S, Theory Of Piezoelectric Flexural Disk Transducers With Applications To Underwater Sound, USL Report No.
490, U. S. Navy Underwater Sound Laboratory, Fort Trumbull, New London
Conn, 1960, S-F001 03 04-1.
[128]Wooster, W A, Breton A, Experimental Crystal Physics, KCMO
Library.
[129]Yang, Variational Formulation Of Piezoelectric Elements.
[130]Ziman J M, Principles Of The Theory Of Solids, Cambridge, 1964,
(Second Edition 1979).
71

Reply to "Piezoelectric Crystals"

Here you can reply to the paste above

Author What's your name?

Title Give your paste a title.

Language What language is your paste written in?

Your paste Paste your paste here

They can be made out of soda ash and cream of tartar and some water. Make thousands of kV just by hitting a crystal!

More info on Piezoelectrics - 
Piezoelectricity
Jim Emery
4/3/97Contents
0.1  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3
0.2  The Linear Piezoelectric Equations  . . . . . . . . . . . . . . .  3
0.3 Piezoelectric Ceramics Constants and One Dimensional Equations  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  6
0.4  The Stress Form of The Piezoelectric Matrix  . . . . . . . . . .  9
0.5 Three Systems of Notation For The Piezoelectric Equations:
Physics (IRE), ABAQUS, and ANSYS  . . . . . . . . . . . . .  11
0.6 Deriving the Tensor Stress Form of the Piezoelectric Equations
From the Strain Form  . . . . . . . . . . . . . . . . . . . . . .  14
0.7  Orthotropic Material  . . . . . . . . . . . . . . . . . . . . . . .  16
0.8  Poling and Piezoelectric Ceramics  . . . . . . . . . . . . . . . .  18
0.9 The Equality of Direct and Converse Piezo Coeﬃcients . . . . 20
0.10  Typical Values of Piezoelectric Parameters  . . . . . . . . . . .  21
0.11 Piezoelectric Current and Impedance From A Finite Element
Calculation  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22
0.12  Piezoelectric Cube Example  . . . . . . . . . . . . . . . . . . .  24
0.13  Generating Ferroelectric Hysteresis Curves  . . . . . . . . . . .  26
0.14  The Strain Hysteresis Curve  . . . . . . . . . . . . . . . . . . .  26
0.15  Relaxors  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  26
0.16  Electrostriction  . . . . . . . . . . . . . . . . . . . . . . . . . .  27
0.17 The Conversion Program: pzansys.ftn . . . . . . . . . . . . .  27
0.18 Example: Parameter Conversion Using Compliance Components 28
0.19 Example: Parameter Conversion Using Four Independent Elasticity Components  . . . . . . . . . . . . . . . . . . . . . . . .  31
0.20  Example, Orthotropic Elastic Constants  . . . . . . . . . . . .  34
0.21 Program Listing of pzansys.ftn Program  . . . . . . . . . . .  37
0.22 Location of Piezoelectric Information In ANSYS Manuals . . . 58
i0.23 ANSYS Comparison For A Composite Piezoelectric Transducer: Example VM176. . . . . . . . . . . . . . . . . . . . . .  59
0.24  Piezoelectricity Bibliography . . . . . . . . . . . . . . . . . . .  61
ii0.1 Introduction
We shall show how to derive the stress form of the piezoelectric equations
from the strain form, and we shall show how to ﬁnd the clamped permittivity
matrix from the unclamped matrix. These conversions are required for ﬁnite
element programs because in such programs equations are solved for the
elastic displacements, and so the strains rather than the stresses are the
independent variables. We shall do the basic derivation three times in the
following sections. We shall ﬁrst treat the one dimensional equation. In
later sections we shall treat the matrix case and the full tensor case. These
derivations are conceptually the same.
0.2 The Linear Piezoelectric Equations
A piezoelectric ceramic is a ferroelectric material. It exhibits hysteresis and
has nonlinear behavior. When this material is poled (placed in a strong
electric ﬁeld), it attains a permanently polarized state. For a small variation
in the electric ﬁeld, it behaves approximately linearly near that state.
The linear piezoelectric equations are written in two diﬀerent forms. In
the more common form, which is called the strain form, the strain tensor εij
and the electric polarization vector Pi are each written as linear functions of
the electric ﬁeld vector Ei and the stress tensor σij
. These equations are
εij = Ekdkij + c

ijkl
σkl
,
and
Pi = 0χijEj + dijkσj k.
The dijk are the components of the piezoelectric tensor (strain coeﬃcients).
The c

ijkl
are the components of the inverse elastic tensor. The constant 0
is the permittivity of free space. The χij are the components of the electric
susceptibility tensor. In the direct piezoelectric eﬀect, when a piezoelectric
material is put under a stress, the material becomes electrically polarized
and surface charges appear. The direct piezoelectric eﬀect is
Pi = dijkσj k,
which is obtained from the second equation, when the external electric ﬁeld
is zero.
3In the converse piezoelectric eﬀect, when a piezoelectric material is put
in an external electric ﬁeld (voltages applied to electrodes), the material
experiences a strain. The converse piezoelectric eﬀect is
εij = Ekdkij
,
which is obtained from the ﬁrst equation, when there are no external forces
and the stress is zero.
The fact that the coeﬃcients for both the direct eﬀect and the converse
eﬀect are the same is a consequence of conservation of energy. This fact can
be established by a thermodynamic argument (See Nye). Because of various
symmetries in the tensor indices, the equations have matrix forms. Let us
deﬁne an index function f : (i, j)− > k that takes a pair of indices and
produces a single index. The function f takes values as follows: f(1, 1) =
1, f(2, 2) = 2, f(3, 3) = 3, f(2, 3) = 4, f(1, 3) = 5, f(1, 2) = 6 The function is
symmetric, f(i, j) = f(j, i), and so f(3, 2) = 4, f(3, 1) = 5, f(2, 1) = 6. This
is the Physics and The Institute of Radio Engineers assignment function. We
deﬁne this function with a table:
f 1 2 3
1 1 6 5
2 6 2 4
3 5 4 3
This index assignment is the one most often used in elasticity theory to
convert the second rank stress and strain tensors to six component vectors.
Although this is the usual convention, the following assignment deﬁned by
function g is used sometimes (ABAQUS).
g 1 2 3
1 1 4 5
2 4 2 6
3 5 6 3
A third index function is also used. The function, which we call h is used
by the ﬁnite element program ANSYS.
4h 1 2 3
1 1 4 6
2 4 2 5
3 6 5 3
Warning! These diﬀerences in index assignment will change the matrices
such as c, d, e.
The matrix components of the piezoelectric tensor are di,f(i,j) and are
either equal to the tensor component itself, when i = j, or to twice the
component. This will be explained below.
There is a second form of the piezoelectric equations, which is called the
stress form.
The stress tensor σij and the electric displacement vector Di
, are each
written as linear functions of the electric ﬁeld vector Ei and the strain tensor
εij
. These equations are
σij = cijklεkl − Ekekij
,
and
Di = eijkεj k + ijEj
.
The cijkl are the components of the elastic tensor. The eijk are the components of the piezoelectric tensor (stress coeﬃcients). The ij are the components of the permittivity tensor.
These equations are suitable for ﬁnite element programs because in these
programs the displacements and hence the strains are the independent variables, whereas the boundary conditions are given as force loadings, that is as
stresses. It will be shown below that the piezoelectric stress coeﬃcients and
the strain coeﬃcients are related as follows:
eijk = cj klmdilm.
Notation varies greatly in the ﬁeld of piezoelectricity. In a following
section we will summarize the notation used in physics, and in two ﬁnite
element programs ABAQUS and ANSYS.
50.3 Piezoelectric Ceramics Constants and One
Dimensional Equations
The principal references for this section are: Jaﬀe, Cook and Jaﬀe, Piezoelectric Ceramics, and the Morgan Matroc document Guide To Modern
Piezoelectric Ceramics.
The original research on piezoelectric materials was conducted with parallel plate capacitors, and the problem was treated as one dimensional. So
the strain form of the equations considered were
ε = Ed + c−1
σ
P = 0χE + dσ,
where the quantities are in the 3 direction. From the deﬁnition of the electric
displacement
D = P + 0E = 0(χ + 1)E + dσ = 
σ
E + dσ.
So we have an alternate form
D = 
σ
E + dσ.
Here 
σ
is the permittivity at constant stress (unclamped).
The second set of piezoelectric equations (stress form) may be obtained
from the ﬁrst set (strain form). Solving the ﬁrst equation for the stress, we
have
σ = cε − cdE = cε − eE.
Making this substitution for σ in the alternate form of the second equation
of the ﬁrst set, we get the second equation of the second set.
D = 
σ
E + d(cε − eE) = E + eε,
= (
σ
− de)E + dcε
= 
ε
E + eε,
where the piezoelectric coeﬃcient e = cd. Here 
ε
is the permittivity at
constant stress (unclamped).
6The coeﬃcients may be deﬁned as partial derivatives. For example,
d = (
∂ε
∂E
)σ,
where the subscript σ means that ε is a function of E and σ, where σ is held
constant in computing the partial derivative. Similarly we have
d = (
∂P
∂σ
)E
The constant g is deﬁned by
g = (−
∂E
∂σ
)ε = (
∂ε
∂D
)σ.
We have
E =
1

(D − dσ)
so
g = −(
∂E
∂σ
)ε =
d

.
The constant e is deﬁned by
e = −(
∂σ
∂E
)ε = (
∂D
∂ε
)E
The poled ceramic is an orthotropic material with a plane of symmetry
whose normal is in the poled direction. Further it is symmetric with respect
to any rotation about the poling direction. Performing a reﬂection through
the symmetry plane, some of the strain components and stress components
are maintained in sign and some others are changed in sign. Writing the
various components of stress as functions of strain, in the two coordinate
systems, we ﬁnd that the proper signs can be maintained only if some of the
elastic coeﬃcients are zero (see pp62-64 Sokolnikoﬀ Mathematical Theory
of Elasticity). Carrying out rotation transformations about the poling axis
(3 direction) we ﬁnd further conditions on the elasticity coeﬃcients. We have
only the following independent elastic coeﬃcients (Jaﬀe et al, p20)
c11, c33, c44, c12, c13.
7We have
c22 = c11, c55 = c44, c66 = 2(c11 − c12), c23 = c13.
Coeﬃcients which are not in the upper 3-dimensional submatrix or on the
main diagonal are zero.
The independent piezoelectric coeﬃcients are
d33, d31, d15.
We have
d32 = d31, d24 = d15.
All other coeﬃcients are zero.
The independent dielectric coeﬃcients are
K11, K33
We have
K22 = K11.
All other coeﬃcients are zero.
The permittivity is
ij = 0Kij
.
The electromechanical coupling constant is deﬁned in terms of the ratio
of stored electrical to stored mechanical energy. Let Qmech be the electrical
energy converted to mechanical energy and let Qinput be the electrical energy
input. Then the electromechanical coupling factor is
k
2
=
Qmech
Qinput
.
It is claimed (Jaﬀe et al p7) that the relation between the unclamped dielectric constant Kσ
and the clamped dielectric constant Kε
is
K
ε
= K
σ
(1 − k
2
).
Because 0 < k  < 1, the clamped dielectric constant is smaller than the
clamped constant
Kε
< Kσ
.
80.4 The Stress Form of The Piezoelectric Matrix
Given the strain form of the piezoelectric equations, we can convert to the
stress form and give expressions for the three independent e tensor parameters. In matrix form we have
ε = d
T
E + c−1
σ
P = 0χE + dσ
From the ﬁrst equation multiplying by the elasticity matrix c we have
cε = cd
T
E + σ,
so solving for the stress matrix σ we have
σ = cε − cd
T
E
= cε − eE,
where the matrix e is deﬁned as the product of c and the transpose of d
e = cd
T
.
We have
e
T
= (cd
T
)
T
= c
T
d = cd
The last equality follows because c is symmetric. From the second equation
D = P + 0E = 0χE + dσ + 0E
= 0(χ + I)E + dσ
= 
σ
E + dσ,
where 
σ
is the permittivity at constant stress (unclamped permittivity).
Substituting the expression for the stress σ from above
D = 
σ
E + d(cε − eE)
= (
σ
− de)E + dcε
9= 
ε
E + e
T
ε,
where

ε
= 
σ
− de
is the permittivity at constant strain (clamped permittivity). In the previous
section we saw that the elastic matrix c has 5 independent parameters and
d has three independent parameters. We have
c =








c11 c12 c13 0 0  0
c12 c11 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)








,
and
d =



0 0  0 0 d15 0
0 0  0 d15 0 0
d31 d31 d33 0 0 0


 .
Then the piezoelectric stress matrix is
e = cd
T
=






c11 c12 c13 0 0  0
c12 c11 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)












0 0 d31
0 0 d31
0 0 d33
0 d15 0
d15 0 0
0 0 0






=








0 0 c11d31 + c12d31 + c13d33
0 0 c12d31 + c11d31 + c13d33
0 0 c13d31 + c13d31 + c33d33
0 c44d15 0
c44d15 0 0
0 0  0








10=








0 0 e13
0 0 e13
0 0 e33
0 e42 0
e42 0 0
0 0 0








So the three independent piezo e parameters are e13, e33, e42. (Warning! The
e matrix is sometimes deﬁned as the transpose of this e matrix.)
Explicitly we have
e13 = (c11 + c12)d31 + c13d33
e33 = (c13 + c13)d31 + c33d33
and
e42 = c44d15.
These equations are easily solved for d33, d31 and d15 in terms of the three
independent e and ﬁve independent c coeﬃcients.
0.5 Three Systems of Notation For The Piezoelectric Equations: Physics (IRE), ABAQUS,
and ANSYS
The index function f(ij), which maps two indices to one, in the way that
is conventional in elasticity theory, was deﬁned in the previous section, as
was g(ij). Here is a table showing the notation most commonly used in
physics (References, Nye, Auld, Cady, IRE (The Institute of Radio Engineers Standard) ANSYS Procedures p8-15, ), and in the ﬁnite element systems ABAQUS (Reference: Theory manual, section 3.10.1), and ANSYS
(Reference: Theory manual, Piezoelectrics chapter).
11Vari abl e Physics ABAQUS ANSYS
Electric Field Vector Ei Ei Ei
Electric Polarization Vector Pi Pi Pi
Electric Displacement Vector Di
qi Di
Stress Tensor σij σij σij
Stress Vector σk = σf(ij)
- Tk
Strain Tensor εij εij
-
Strain Vector εk - Sk
Electric Permittivity Tensor ij D
φ
ij
ij
Electric Susceptibility Tensor χij
- -
Elasticity Tensor cijkl Dijkl
cijkl
Elasticity Matrix cmn = cf(ij)f(kl) − cmn
Inverse Elasticity Tensor c

ijkl
- c

ijkl
Piezoelectric Tensor (strain) di,jk d
φ
i,jk
di,jk
Piezoelectric Matrix (strain) di,m = di,f(j k) di,m = di,g(j k) di,m
Piezoelectric Tensor (stress) ei,jk e
φ
i,jk
ej k,i
Piezoelectric Matrix (stress) ei,m = ei,f(j k)
ei,m = ei,g(j k)
em,i
The engineering strain coeﬃcients are written as γij
. They diﬀer from
the εij
, in that the shear coeﬃcients are doubled in value.
Note: the physics notation for eijk is deﬁned in equation 8.40 of Auld. Nye
does not introduce the piezoelectric e coeﬃcients. The Cady book is very
old and one can obtain piezoelectric constant notation only by implication.
Note: The ANSYS deﬁnition of the e matrix and tensor is the transpose of
the physics and ABAQUS deﬁnition. Thus the ANSYS e matrix is 6 by 3,
whereas the physics e matrix is 3 by  6.
The two ABAQUS tensor equations are:
σij = Dijklεkl − e
φ
mijEm
(or alternately: σij = Dijklεkl − Dijkld
φ
mkl
Em)
qi = e
φ
ijk + D
φ
ijEj
.
The ANSYS matrix equations are:
T = cS − eE
D = e
T
S + E.
12The Physics matrix equations are:
σ = cε − e
T
E,
D = eε + E.
Again note that the ANSYS e deﬁnition and the physics e deﬁnition are
transposes of one another.
Here are the ANSYS matrices written out:
T = σ =





σ1
σ2
...
σ6





S = ε =





ε1
ε2
...
ε6





E =



E1
E2
E3



D =



D1
D2
D3



c =








c11 c12 c13 0 0 0
c21 c22 c23 0 0 0
c31 c32 c33 0 0 0
0 0  0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66








13e =






e11 e12 e13
e21 e22 e23
e31 e32 e33
e41 e42 e43
e51 e52 e53
e61 e62 e63






 =



11 0 0
0 22 0
0 0 22



We will show in the next section that (physics notation)
eijk = cj klmdilm.
0.6 Deriving the Tensor Stress Form of the
Piezoelectric Equations From the Strain
Form
We have already done this derivation in a previous section for the matrix
forms of the piezoelectric equations. Here we will repeat this for the tensor
forms of the equations. We start with the strain form of the piezoelectric
equations
εij = Ekdkij + c

ijkl
σkl
,
and
Pi = 0χijEj + dijkσj k.
Let cpqij be the elasticity tensor. Because c

is its inverse, we have
cpqij
c

ijkl = δ
p
k
δ
q
l
.
Then multiplying the ﬁrst equation by cpqij and summing we get
cpqijεij = cpqijEkdkij + cpqij
c

ijkl
σkl
,
14which becomes
cpqijεij = cpqijEkdkij + σpq.
We let
ekpq = cpqijdkij
.
Then
σpq = cpqrsεrs − empqEm.
The equation relating electric displacement D, electric polarization P, and
the macro electric ﬁeld E is
Di = Pi + 0Ei
.
Substituting for the polarization from the second piezoelectric equation above,
we have
Di = 0χijEj + dijkσj k + 0Ei
.
We write this as
Di = 0(χij + δij
)Ej + dijkσj k
= 
σ
ijEj + dipqσpq,
where 
σ
ij
is the zero stress (unclamped ) permittivity. Substituting our previously derived expression for the stress, we get
Di = 
σ
imEm + dipq(cpqrsεrs − empqEm)
= (
σ
im − dipqempq)Em + dipqcpqrsεrs
= eirsεrs + 
ε
imEm,
where

ε
im = 
σ
im − dipqempq
is the zero strain (clamped) permittivity. We have assumed that cpqrs = crspq.
This is the case when the internal energy of the material is a function of the
values σ, ε, D, and E, and not dependent on the path through which the
state was reached, that is,when the internal energy is an exact diﬀerential in
the thermodynamic sense. Therefore redeﬁning the dummy indices, we have
the stress version of the piezoelectric equations:
σij = cijklεkl − Ekekij
,
15and
Di = eijkεj k + 
ε
ij
Ej
.
The matrix forms of these equations are:
σ = cε − e
T
E,
D = eε + 
ε
E.
0.7 Orthotropic Material
For the basic equations of elasticity, see Jim Emery Elasticity, (See ﬁles:
elastic.tex, elastic.ps).
An orthotropic problem is one in which the material properties are symmetric with respect to 3 mutually orthogonal planes (See Sokolnikoﬀ p62). A
piezoelectric ceramic such as PZT is orthotropic, but has even more symmetry. It is invariant to any rotation about the poled axis. By applying various
symmetry transformations, and invariance properties, one ﬁnds that many
of the elastic constants must be zero. The elastic matrix becomes
c =








c11 c12 c13 0 0 0
c21 c22 c23 0 0 0
c31 c32 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66








.
The inverse matrix is
c

= c−1
=








1/E1 −ν21/E2 −ν31/E3 0 0  0
−ν12/E1 1/E2 −ν32/E3 0 0  0
−ν13/E1 −ν23/E2 1/E3 0 0  0
0 0 0 1/(2µ32) 0  0
0 0 0 0 1/(2µ31) 0
0 0 0 0 0 1/(2µ12)








.
The Poisson ratio, νij
, is the ratio of the negative of the strain in the
xj direction, to the strain in the xi direction, for a normal stress in the xi
direction. For example, if the only nonzero stress component is σ1, then
16





1/E1 −ν21/E2 −ν31/E3 0 0 0
−ν12/E1 1/E2 −ν32/E3 0 0 0
−ν13/E1 −ν23/E2 1/E3 0 0 0
0 0 0 1/(2µ32) 0  0
0 0 0 0 1/(2µ31) 0
0 0 0 0 0 1/(2µ12)












σ1
0
0
0
0
0






=








σ1/E1
(−σ1ν12)/E1
(−σ1ν13)/E1
0
0
0








=






ε1
ε2
ε3
0
0
0






.
So that
ν12 =
−ε2
ε1
,
and
ν13 =
−ε3
ε1
,
By the symmetry of c, we have
νij/Ei = νj i/Ej
.
So there are three independent values of Poisson’s ratio. We can take these
to be
ν12, ν13, ν23,
or
ν21, ν31, ν32.
17For the case of poled piezoelectric ceramics poled in the 3 direction, the best
choice is
ν21, ν31, ν32.
It is clear that in that case
ν31 = ν32,
so that there are only two independent poisson ratio’s, which we may take
as
ν21, ν31.
Also we see that ν31 is the poisson ratio stated for piezoceramics and is said
to be about
.31
If the engineering strain vector is used then the shear strains are twice the
physics shear strains, so that the 2 multiplying a shear modulus is suppressed.
The engineering shear strain is equal to the shear angle, the physics shear
strain is equal to one half the angle.
The notation for shear modulus µij depends upon which of the index
mappings f or g is used in assigning the stress-strain tensor indices to the
vector indices.
0.8 Poling and Piezoelectric Ceramics
See:
Newnham R E, Trolier-McKinstry S,Giniewicz J R
Piezoelectric, Pyroelectric and Ferroic Crystals
J. Material Education, 15, 189-223 (1993),
(reprint: Penn State MRL annual report to ONR 1994, Appendix 1)
For the nature of the strain vs electric ﬁeld curve, the nature of domains
in PZT, and the morphotropic phase boundary between tetragonal PbTiO3
and rhombohedral PbZrO3 ferroelectric phases, see:
Subbarao E C, Srikanth V,Cao W, Cross L E
Domain Switching And Microcracking During
Poling Of Lead Zirconate Titonate Ceramics
18Ferroelectrics 1993 V 145 pp. 271-281
(reprinted: Materials For
Adaptive Structural Acoustic Control, Office
Of Navel Research, Penn State Materials
Research Laboratory, L. Eric Cross ed, V1 1994, appendix 15).
The morphrotropic phase boundary is the boundary between two different structures of the piezoelectric material. We note that “morph” is a
root meaning form, and “trop” is a root meaning to turn or change. Hence
morphotropic means: a change in form.
Consider the equation
Pi = dijσj
.
Let the poling direction be the z axis. Then the domains will be randomly
directed in the xy plane. Hence the ceramic will be orthotropic relative to
the xy plane. Because of this symmetry, the piezoelectric matrix takes the
form (see Jaﬀe, Cook, Jaﬀe Piezoelectric Ceramics, p20)
d =



0 0  0 0 d15 0
0 0  0 d24 0 0
d31 d32 d33 0 0 0


 .
=



0 0  0 0 d15 0
0 0  0 d15 0 0
d31 d31 d33 0 0 0


 .
The elasticity tensor takes the form
c =








c11 c12 c13 0 0  0
c12 c22 c13 0 0  0
c13 c13 c33 0 0  0
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 2(c11 − c12)








.
The dielectric tensor is 


K1 0 0
0 K1 0
0 0 K3


 .
19Suppose we have d31 = d32, d33, d15 = d24,and dij = 0,  otherwise.  Then  we
have
d311 = d31,
d322 = d32,
d333 = d33,
d113 = d15/2,
d131 = d15/2,
d223 = d24/2,
d232 = d24/2,
and dijk = 0  otherwise.
0.9 The Equality of Direct and Converse Piezo
Coeﬃcients
This equality may be shown by applying the ﬁrst and second laws of thermodynamics. Let the stress coeﬃcients σij
, the electric ﬁeld components Ei and
the temperature T be the state variables. Then write diﬀerential expressions
for the strain deij
, the electric displacement dDi
, and the entropy dS. Write
the internal energy change as
dU = σijdeij + EidDi + T dS.
Introduce the function
Φ = U − σij eij − EiDi − T S.
Compute the diﬀerential of Φ using the last two equations, and also write it
in the general partial derivative form. Then equate the partial derivatives to
coeﬃcients of the equations. For example
∂Φ
∂Ei
= −Di
.
Diﬀerentiate to equate derivatives of eij
, Di
, or S to second order partial
derivatives of Φ. This shows for example that the coeﬃcients in the direct
piezoelectric eﬀect are equal to the coeﬃcients of the converse piezoelectric
eﬀect. See Nye p179.
200.10 Typical Values of Piezoelectric Parameters
Dielectric constant
K = 200 =

0
0 = 8.85 × 10
−12
farad/meter.
KT
dielectric constant at constant stress. KS
dielectric constant at zero
strain (constrained piezoelectric coeﬃcient)
d33 = 200 × 10−12
coulomb/newton (or meter/volt).
d31 = −100 × 10−12
coulomb/newton (or meter/volt).
d15 = 400 × 10−12
coulomb/newton (or meter/volt).
gij voltage coeﬃcients.
kij = gij
Poisson’s ratio is approximately .31 for all ceramics. Curie temperature
Tc = 300.
fr = resonant frequency. Typical value is 70 kHz.
Panasonic motor.
Operating voltage 22 volts.
Rated output 1 watt.
Torque 250 gf.cm (gf gram force)
250 gf cm (pound/454 gf)(inch/2.54 cm)(16 ounce/pound) = 3.47 inch.ounce
Current .4 amp at 5 to 10 volts
Input power 2 to 4 watts
Voltage of control circuit 12 volts
250 g.cm (kg/1000 g)(9.8 m/s.s)(m/100 cm) = .0245 N.m
Angular velocity = 4002π/(min)(min/60sec) =  10.47(1/s)
Power = (.0245)(10.47) = .2565 watt
Speed 400 rpm, noload 600 rpm.
fa = antiresonant frequency (parallel resonance)
Electromagnetic coupling constant
k = (electrical energy stored)/(mechanical energy stored).
21kp = .2 to .8 = planar coupling factor of a thin disk.
Qm = quality factor.
Example: Thickness
l = .030inches = 7.63 × 10−4
m
v = 400volts.
E = v/l = 3.67 × 10
6
v/m.
e = dE = (200 × 10−12
)(3.67 × 10
6
) = 7.34 × 10−4
.
∆z = el = 22µinches
0.11 Piezoelectric Current and Impedance From
A Finite Element Calculation
Let D be the electric displacement and ρ the charge  density.  We can  use one
of Maxwell’s equations
∇ · D = ρ,
to show that the surface charge density on a conductor equals the normal
component Dn of D outside of the conductor. We enclose a conducting
surface with a thin pillbox of area A, volume V , and thickness ∆t. The
charge contained in the volume is Aτ where τ is the surface charge density.
Then letting ∆t go to zero, we have
Aτ =

V
ρdV =

V
∇ · DdV =

∂ V
D · nds → D · nA.
Then the charge density is
τ = D · n = Dn.
We compute the nodal strains εij
from the nodal displacements ui at a
face of an element.
From the ﬁnite element nodal potentials φn, we compute derivatives of
the shape functions, and then the gradient of the potential. Then we obtain
the electric ﬁeld Ei as the negative gradient of the potential,
E = −∇φ.
22We compute the electric displacement from the piezoelectric equation
Di = eijkεj k + KijEj
.
Then we compute the normal component Dn. This gives the surface charge
density. Multiplying by the area A, we get the charge q. We sum up the
charges on the boundaries of the elements to get the total electrode charge.
If we calculate the charge as a function of time, we can ﬁt it to a function
Q0 cos(ωt + ψ),
and thus to the real part of the exponential
Q0 exp(j(ωt + ψ)).
Diﬀerentiating with respect to time we get the current
I = I0 exp(j(ωt + ψ)).
Then if V is the applied voltage, the impedance is
Z = V /I .
We can also compute energy storage and dissipation in the ﬁnite element
model, and relate it to energy storage and dissipation in an equivalent circuit.
We can obtain or compute the energy density of the electric ﬁeld
D · E
2
,
the elastic strain energy, the kinetic energy, heat dissipation, and frictional
contact energy. The kinetic energy will correspond to inductor energy in
the equivalent circuit, the strain energy and the electric ﬁeld energy will
correspond to capacitive energy, and the dissipation energy will correspond
to a resistance loss. If the energies correspond, then the equivalent circuit is
valid.
230.12 Piezoelectric Cube Example
Consider a piezoelectric cube of side a. Let thin metal plates be located at
z = 0  and z = a. Let there be a voltage V on these plates. Let there be a
horizontal compressive forces F in the y direction acting on the y = 0,  and
the y = a faces of the block. The electric ﬁeld is
E =



E1
E2
E3


 =



0
0
V /a



The stress is
σ =






σ1
σ2
σ3
σ4
σ5
σ6






=






F /a
2
0
0
0
0
0






.
An average dielectric constant for PZT is
K = 200 =

o
= 1 + χ.
We have
0 = 8.85 × 10
10
(C
2
/N M2
).
Then the permittivity is
 = 1.77 × 10−9
.
If we take the material to be isotropic, the polarization equation becomes
P = 0χE +



0 0 0 0 d15 0
0 0 0 0 d25 0
d31 d32 d33 0 0  0











σ1
σ2
σ3
σ4
σ5
σ6








We have
d31 = d32 = −100 × 10−12
C/N or M/V
24d33 = 2100 × 10−12
d51 = d15 = 400 × 10−12
Because only σ1 is nonzero, we have
P = 0χE +



0
0
d31σ1


 .
Then
P3 = 0χE3 + d31σ1
D = 0E + P.
D3 = 0(1 + χ)E3 + d31σ1 = 0K E3 + d31σ1.
The surface charge on the z = a face is
Q = D3a
2
.
The strain is
e =








0 0 d31
0 0 d32
0 0 d33
0 d25 0
d15 0 0
0 0 0











E1
E2
E3


 +








c11 c12 c13 0 0 0
c21 c22 c33 0 0 0
c31 c32 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66
















σ1
σ2
σ3
σ4
σ5
σ6








So the nonzero elements of the strain are
e1 = d31E3 + σ1c11 = d31E3 +
σ1
ξ
,
e2 = d32E3 + σ1c21 = d32E3 −
σ1ν
ξ
,
e3 = d33E3 + σ1c31 = d33E3 −
σ1ν
ξ
,
where ξ is Young’s modulus, and ν is Poisson’s ratio.
250.13 Generating Ferroelectric Hysteresis Curves
Suppose we arrange two capacitors in series, one C1 containing a ferroelectric
dielectric, and the second C2 a nonferroelectric dielectric. Let us connect in
series an alternating high voltage source. Let V1 be the voltage across capacitor C1. This is a measure of the E ﬁeld in the ferroelectric capacitor C1. The
charge on the plates of C1 is proportional to the D ﬁeld (actually the D ﬁeld
is equal to the surface charge density). But the charge on C1 is equal to the
charge on C2 and moves from C1 to C2 and back again during a cycle of the
applied alternating voltage. But the second capacitor is nonferroelectric and
its dielectric is linear, so that the voltage V2 across the second capacitor is
proportional to charge and hence to the D ﬁeld of the ferroelectric capacitor.
Therefore, if we connect V1 across the horizontal plates of an oscilloscope,
and V2 across the vertical plates, we will see the hysteresis curve. The polarization ﬁeld will be approximately equal to D, since the permittivity  for
the ferroelectric is much larger than the permittivity of free space 0.
(I think this is close to the method given by Sawyer and Tower, see
bibliography.)
0.14 The Strain Hysteresis Curve
The strain vs electric ﬁeld curve for a piezoelectric ceramic is a butterﬂy
curve with the strain always of one sign, and the curve lying in the positive
upper half plane.
0.15 Relaxors
A relaxor ferroelectric has some of the properties of a normal ferroelectric,
but is not so permanently polarized and has a very narrow hysteresis loop.
When the ﬁeld is removed, the material almost returns to an unpolarized
state. It requires a bias ﬁeld.
See:
Cross L Eric
Relaxor Ferroelectrics: An Overview
Ferroelectrics, 1994 v151 pp. 305-320
26(reprinted: annual report 1994 from MRL at Penn State to ONR v1 p305)
0.16 Electrostriction
Electrostriction is an internal stress caused by the force of the electric ﬁeld on
charges. It occurs in all materials, not just piezoelectric ones. A computation
shows it to be related to the gradient of the dielectric constant, and in practice
the stress is proportional to the square of the electric ﬁeld intensity. It is very
hard to measure, since the force can be masked by the force due to the free
charges on the capacitor plates.
Notes and references to electrostriction: Wert, Thomson Physics of
Solids, p. 321. Julius Adams Stratton, Electromagnetic Theory, pp149-
151. Panofsky and Philips Classical Electricity and Magnetism p.95,
As a term of an equation. Maxwell stress tensor. Classical derivation. The
treatment in Panofsky and Phillips is based on the treatment in Becker:
(formerly called Becker and Abraham) Electromagnetic Fields and Interactions p. 134. In Jaﬀe: Piezoelectric Ceramics, Electrostriction is
mentioned on pages 15 and 78. He relates electrostriction to domain reversal. Jaﬀe says that the eﬀect is signiﬁcant only for ferroelectrics above the
Curie point. Cady: Piezoelectricity deﬁnes electrostriction to be that effect that is proportional to E
2
. See pages 4, 198, 199, 614. Cady states
that electrostriction is only signiﬁcant for ﬁelds above 20,000 volts per cm.
See also Auld, Acoustic Fields and Waves in Solids, volumes I and II,
who presents a nice one dimensional model of the piezoelectric eﬀect on p265.
0.17 The Conversion Program: pzansys.ftn
This program does the conversion calculations derived in the previous sections. It requires as input the d parameters, and either the elasticity matrix,
or its inverse. It computes the e parameters, and computes the clamped
permittivity (strain constant) from the unclamped permittivity. It generates
material property cards for the ANSYS program. We list the program in a
later section. The next sections illustrate the operation of the program.
270.18 Example: Parameter Conversion Using
Compliance Components
yonada:/users/u51195 $ pzansys
Reference: Jim Emery "Piezoelectricity",
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
Enter the orthotropic elastic constants
for a coordinate system in the principle directions
Enter Young’s Modulus, in x [ 0.12800E+12]
28Enter Young’s Modulus, in y [ 0.12800E+12]
Enter Young’s Modulus, in z [ 0.11000E+12]
Enter ANSYS Poisson ratios (Theory sec. 2.1)
Enter Major Poisson ratio prxy[ 0.30000 ]
Enter Major Poisson ratio prxz[ 0.41976 ]
Enter Major Poisson ratio pryz[ 0.41976 ]
Enter shear modulus g23 [ 0.52000E+11]
Enter shear modulus g13 [ 0.52000E+11]
Enter shear modulus g12 [ 0.12200E+12]
Compliance matrix (elasticity inverse) cinv=
0.7813E-11 -0.2344E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.2344E-11 0.7813E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3279E-11 -0.3279E-11 0.9091E-11 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Note: In ANSYS components 4,5,6 are permuted
(A different mapping from tensor to vector is used):
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
The five independent compliance coeficients are:
s11= 0.78125000E-11
s33= 0.90909091E-11
s44= 0.19230769E-10
s12=-0.23437500E-11
s13=-0.32793388E-11
Computing inverse matrix
Elasticity matrix =
0.2104E+12 0.1119E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1119E+12 0.2104E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1163E+12 0.1163E+12 0.1939E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.12800000E+12
Young’s Modulus 2= 0.12800000E+12
Young’s Modulus 3= 0.11000000E+12
Major Poisson ratio prxy= 0.30000000
Major Poisson ratio prxz= 0.41975537
Major Poisson ratio pryz= 0.41975537
Shear Modulus g23= 0.52000000E+11
Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
29Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -8.977066554171117
e33= 15.52345455897376
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.4900E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permitivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permitivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.3950E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.12800000E+12
EY, 2, 0.12800000E+12
EZ, 2, 0.11000000E+12
30/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.30000000
PRXZ, 2, 0.41975537
PRYZ, 2, 0.41975537
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.39498958E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -8.9770666
TBDATA, 6, -8.9770666
TBDATA, 9, 15.523455
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.19 Example: Parameter Conversion Using
Four Independent Elasticity Components
yonada:/users/u51195 $ pzansys
Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals, and so on)
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
31c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
2
Enter c11 [ 0.13200E+12]
Enter c33 [ 0.11500E+12]
Enter c44 (ANSYS c66) [ 0.52000E+11]
Enter c12 [ 0.71000E+11]
Enter c13 [ 0.73000E+11]
Elasticity matrix c=
0.1320E+12 0.7100E+11 0.7300E+11 0.0000E+00 0.0000E+00 0.0000E+00
0.7100E+11 0.1320E+12 0.7300E+11 0.0000E+00 0.0000E+00 0.0000E+00
0.7300E+11 0.7300E+11 0.1150E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Note: In ANSYS components 4,5,6 are permuted
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
Compliance matrix=
0.1273E-10 -0.3665E-11 -0.5754E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3665E-11 0.1273E-10 -0.5754E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.5754E-11 -0.5754E-11 0.1600E-10 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.78561263E+11
Young’s Modulus 2= 0.78561263E+11
Young’s Modulus 3= 0.62497537E+11
Major Poisson ratio prxy= 0.28788955
Major Poisson ratio prxz= 0.45203533
Major Poisson ratio pryz= 0.45203533
Shear Modulus g23= 0.52000000E+11
Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
32Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 8.400
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -5.700000000000001
e33= 8.400000000000000
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 -5.700
0.0000E+00 0.0000E+00 8.400
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.2820E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permitivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permitivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.6030E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.78561263E+11
EY, 2, 0.78561263E+11
EZ, 2, 0.62497537E+11
33/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.28788955
PRXZ, 2, 0.45203533
PRYZ, 2, 0.45203533
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.60300000E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -5.7000000
TBDATA, 6, -5.7000000
TBDATA, 9, 8.4000000
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.20 Example, Orthotropic Elastic Constants
Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)
Document location:/u51195/ps/piezoelc.ps,piezoelc.tex
Anonymous FTP: ftp.os.kcp.com /pub/emery/
Units are metric: Newton, Meter, Kilogram
stress: Newton/M^2 (Pascals)
Young’s Modulus: Pascals
Elasticity coefficients: Pascals, and so on)
Enter piezoelectric d constants (coulomb/newton)
The IRE convention for piezoelectric labeling is used
Type <return> for default value
Enter d33 [ 2.0000000000000001E-10]
Enter d31 [ -1.0000000000000000E-10]
Enter d15 [ 4.0000000000000001E-10]
d matrix=
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.4000E-09 0.0000E+00 0.0000E+00
-0.1000E-09 -0.1000E-09 0.2000E-09 0.0000E+00 0.0000E+00 0.0000E+00
Do you want to see the nonzero components
of the 27 d tensor components? [n]
Define the elasticity matrix or compliance matrix:
(1)Use elastic constants to define compliance
for orthotropic material, using Young’s moduli,
Poissons ratios, and Shear Moduli.
(2)Use the five components of the elastic matrix
c11,c33,c44 (c44 is ANSYS c66),c12,c13,
for a z-poled ceramic.
(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,
and c_{23}=c_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
34(3)Use the five components of the compliance matrix
s11,s33,s44 (s44 is ANSYS s66),s12,s13,
for a z-poled ceramic.
(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,
and s_{23}=s_{13}. Coefficients which are
not in the upper 3-dimensional submatrix,
not on the main diagonal, are zero.
type <return> for default selection: [1]
Enter the orthotropic elastic constants
for a coordinate system in the principle directions
Enter Young’s Modulus, in x [ 0.12800E+12]
Enter Young’s Modulus, in y [ 0.12800E+12]
Enter Young’s Modulus, in z [ 0.11000E+12]
Enter ANSYS Poisson ratios (Theory sec. 2.1)
Enter Major Poisson ratio prxy[ 0.30000 ]
Enter Major Poisson ratio prxz[ 0.41976 ]
Enter Major Poisson ratio pryz[ 0.41976 ]
Enter shear modulus g23 [ 0.52000E+11]
Enter shear modulus g13 [ 0.52000E+11]
Enter shear modulus g12 [ 0.12200E+12]
Compliance matrix (elasticity inverse) cinv=
0.7813E-11 -0.2344E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.2344E-11 0.7813E-11 -0.3279E-11 0.0000E+00 0.0000E+00 0.0000E+00
-0.3279E-11 -0.3279E-11 0.9091E-11 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1923E-10 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8197E-11
Note: In ANSYS components 4,5,6 are permuted
(A different mapping from tensor to vector is used):
ANSYS 44 is 66
ANSYS 55 is 44
ANSYS 66 is 55
The five independent compliance coeficients are:
s11= 0.78125000E-11
s33= 0.90909091E-11
s44= 0.19230769E-10
s12=-0.23437500E-11
s13=-0.32793388E-11
Computing inverse matrix
Elasticity matrix =
0.2104E+12 0.1119E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1119E+12 0.2104E+12 0.1163E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.1163E+12 0.1163E+12 0.1939E+12 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.5200E+11 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1220E+12
Orthotropic constants from compliance matrix:
Young’s Modulus 1= 0.12800000E+12
Young’s Modulus 2= 0.12800000E+12
Young’s Modulus 3= 0.11000000E+12
Major Poisson ratio prxy= 0.30000000
Major Poisson ratio prxz= 0.41975537
Major Poisson ratio pryz= 0.41975537
Shear Modulus g23= 0.52000000E+11
35Shear Modulus g13= 0.52000000E+11
Shear Modulus g12= 0.12200000E+12
Do you want to see the nonzero components
of the 81 elasticity tensor components of c? [n]
piezoelectric (stress) tensor coefficients
Do you want to see the nonzero components
of the 27 piezoelectric tensor components of e? [n]
Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
Second calculation:
e13= -8.977066554171117
e33= 15.52345455897376
e42= 20.80000000000000
Note: ANSYS e rows are permuted
ANSYS row 4 is row 6
ANSYS row 5 is row 4
ANSYS row 6 is row 5
ANSYS Piezoelectric e matrix=
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 -8.977
0.0000E+00 0.0000E+00 15.52
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 20.80 0.0000E+00
20.80 0.0000E+00 0.0000E+00
d times e =
0.8320E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.8320E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.4900E-08
Enter the unclamped dielectric constant K11
[ 1300.0000 ]
Enter the unclamped dielectric constant K33
[ 1000.0000 ]
Unclamped electric permittivity tensor=
0.1151E-07 0.0000E+00 0.0000E+00
0.0000E+00 0.1151E-07 0.0000E+00
0.0000E+00 0.0000E+00 0.8850E-08
Clamped electric permittivity tensor=
0.3185E-08 0.0000E+00 0.0000E+00
0.0000E+00 0.3185E-08 0.0000E+00
0.0000E+00 0.0000E+00 0.3950E-08
Enter density Kg/m^3 [7500]
/COM
/COM Material properties
/COM Young’s Moduli
EX, 2, 0.12800000E+12
EY, 2, 0.12800000E+12
EZ, 2, 0.11000000E+12
/COM ANSYS "Major" Poisson ratios
PRXY, 2, 0.30000000
36PRXZ, 2, 0.41975537
PRYZ, 2, 0.41975537
/COM Density Kg/m^3
/COM DENS, 2, 0.75000000E-02
/COM clamped Permitivities
MP,PERX, 2, 0.31850000E-08
MP,PERY, 2, 0.31850000E-08
MP,PERZ, 2, 0.39498958E-08
/COM ANSYS Piezoelectric "e" matrix
TB,PIEZ,3
TBDATA, 3, -8.9770666
TBDATA, 6, -8.9770666
TBDATA, 9, 15.523455
TBDATA, 14, 20.800000
TBDATA, 16, 20.800000
0.21 Program Listing of pzansys.ftn Program
Here is a listing of the program. The program works through tensor summations, matrix calculations, and matrix inversions. Each subroutine has a one
line description and a description of the parameters.
% cat pzansys.ftn
c pzansys.ftn stress and strain versions of the piezoelectric tensor
c Version 2/6/97
c Author: Jim Emery (AlliedSignal Kansas City)
c for the ansys finite element program.
column 11111111112222222222333333333344444444445555555555666666666677777777778
c2345678901234567890123456789012345678901234567890123456789012345678901234567890
c Reference: Jim Emery "Piezoelectricity" (File: piezoelc.tex)
implicit real*8 (a-h,o-z)
dimension d(3,3,3)
dimension e(3,3,3)
dimension c(3,3,3,3)
dimension dc(3,3)
dimension cmat(6,6)
dimension cmatinv(6,6)
dimension a(6,6),b(6,6),cc(6,6)
dimension dmat(3,6)
dimension emat(6,3)
dimension epsilon(3,3)
dimension clamped(3,3)
dimension ain(10)
dimension pr(3)
dimension ym(3)
real*8 nu12,nu23,nu13,nu21,nu32,nu31
character*1 cm
character*1 reply
character*53 cs1
zero=0.
c piezoelectric constants (coulomb/newton)
37do i=1,3
do j=1,6
dmat(i,j)=0.
enddo
enddo
dmat(3,3)=200e-12
dmat(3,1)=-100e-12
dmat(1,5)=400e-12
write(*,*)’Reference: Jim Emery "Piezoelectricity", (piezoelc.ps)’
write(*,*)’Document location:/u51195/ps/piezoelc.ps,piezoelc.tex ’
write(*,*)’Anonymous FTP: ftp.os.kcp.com /pub/emery/ ’
write(*,*)’Units are metric: Newton, Meter, Kilogram’
write(*,*)’stress: Newton/M^2 (Pascals)’
write(*,*)’Young’’s Modulus: Pascals’
write(*,*)’Elasticity coefficients: Pascals, and so on)’
write(*,*)’Enter piezoelectric d constants (coulomb/newton)’
write(*,*)’The IRE convention for piezoelectric labeling is used’
write(*,*)’Type <return> for default value’
write(*,*)’Enter d33 [’,dmat(3,3),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(3,3)=ain(1)
endif
write(*,*)’Enter d31 [’,dmat(3,1),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(3,1)=ain(1)
endif
dmat(3,2)=dmat(3,1)
write(*,*)’Enter d15 [’,dmat(1,5),’]’
call readr(0,ain,nr)
if(nr .ge. 1)then
dmat(1,5)=ain(1)
endif
dmat(2,4)=dmat(1,5)
write(*,*)’ d matrix= ’
do i=1,3
write(*,’(6(g11.4,1x))’)(dmat(i,j),j=1,6)
enddo
c construct the components of the d tensor
do i=1,3
do j=1,3
do k=1,3
call t2vec(j,k,jk)
if(j .eq. k)then
d(i,j,k)=dmat(i,jk)
else
d(i,j,k)=dmat(i,jk)/2.
endif
enddo
enddo
enddo
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 27 d tensor components? [n]’
read(*,’(a)’)reply
38if(lenstr(reply) .eq. 0)then
reply=’n’
endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
if(d(i,j,k) .ne. zero)then
write(*,’(a,i1,i1,i1,a,g15.8)’)’d(’,i,j,k,’)=’,d(i,j,k)
endif
enddo
enddo
enddo
endif
c orthotropic elastic constants
write(*,*)’Define the elasticity matrix or compliance matrix:’
write(*,*)
write(*,*)’(1)Use elastic constants to define compliance’
write(*,*)’for orthotropic material, using Young’’s moduli,’
write(*,*)’Poissons ratios, and Shear Moduli.’
write(*,*)
write(*,*)’(2)Use the five components of the elastic matrix’
write(*,*)’c11,c33,c44 (c44 is ANSYS c66),c12,c13, ’
write(*,*)’for a z-poled ceramic.’
write(*,*)’(Note: c22=c11, c55=c44, c66=2(c11-c12),c23=c13,’
write(*,*)’and c_{23}=c_{13}. Coefficients which are’
write(*,*)’not in the upper 3-dimensional submatrix,’
write(*,*)’not on the main diagonal, are zero.’
write(*,*)
write(*,*)’(3)Use the five components of the compliance matrix’
write(*,*)’s11,s33,s44 (s44 is ANSYS s66),s12,s13,’
write(*,*)’for a z-poled ceramic.’
write(*,*)’(Note: s22=s11, s55=s44, s66=2(s11-s12),s23=s13,’
write(*,*)’and s_{23}=s_{13}. Coefficients which are’
write(*,*)’not in the upper 3-dimensional submatrix,’
write(*,*)’not on the main diagonal, are zero.’
write(*,*)’type <return> for default selection: [1]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’1’
endif
c2345678901234567890123456789012345678901234567890123456789012345678901234567890
c begin case1
if(reply .eq. ’1’)then
write(*,*)’Enter the orthotropic elastic constants’
write(*,*)’for a coordinate system in the principle directions’
e1=12.8e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in x [’,e1,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
39e1=ain(1)
endif
e2=12.8e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in y [’,e2,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
e2=ain(1)
endif
e3=11.0e10
write(*,’(a,g12.5,a)’)’Enter Young’’s Modulus, in z [’,e3,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
e3=ain(1)
endif
c write(*,’(a)’)’NOTE. The poisson ratio nuij is the ’
c write(*,’(a)’)’the strain ratio -epsilon_j/epsilon_i ’
c write(*,’(a)’)’under a normal stress sigma_i’
write(*,*)
c nu31=.31
write(*,*)’Enter ANSYS Poisson ratios (Theory sec. 2.1)’
prxy=.3
prxz=.41975537
pryz=.41975537
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio prxy[’,prxy,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxy=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio prxz[’,prxz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxz=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter Major Poisson ratio pryz[’,pryz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
prxz=ain(1)
endif
c g23=(e1*e2)/(e1+e2+2.*nu12*e1)
g23=5.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g23 [’,g23,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g23=ain(1)
endif
c g13=(e1*e2)/(e1+e2+2.*nu12*e1)
40g13=5.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g13 [’,g13,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g13=ain(1)
endif
c g12=(e1*e2)/(e1+e2+2.*nu12*e1)
g12=12.2e10
write(*,’(a,g12.5,a)’)’Enter shear modulus g12 [’,g12,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
g12=ain(1)
endif
c build the compliance matrix
do i=1,6
do j=1,6
cmatinv(i,j)=0.
enddo
enddo
cmatinv(1,1)=1./e1
cmatinv(2,2)=1./e2
cmatinv(3,3)=1./e3
cmatinv(1,2)=-prxy/e1
cmatinv(2,1)=cmatinv(1,2)
cmatinv(1,3)=-prxz/e1
cmatinv(3,1)=cmatinv(1,3)
cmatinv(2,3)=-pryz/e2
cmatinv(3,2)=cmatinv(2,3)
cmatinv(4,4)=1/g23
cmatinv(5,5)=1/g13
cmatinv(6,6)=1/g12
write(*,*)’Compliance matrix (elasticity inverse) cinv= ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’(A different mapping from tensor to vector is used):’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
write(*,*)’The five independent compliance coeficients are:’
write(*,’(a,g15.8)’)’s11=’,cmatinv(1,1)
write(*,’(a,g15.8)’)’s33=’,cmatinv(3,3)
write(*,’(a,g15.8)’)’s44=’,cmatinv(4,4)
write(*,’(a,g15.8)’)’s12=’,cmatinv(1,2)
write(*,’(a,g15.8)’)’s13=’,cmatinv(1,3)
do i=1,6
do j=1,6
a(i,j)=cmatinv(i,j)
enddo
41enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
write(*,*)’Computing inverse matrix’
call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
endif
do i=1,6
do j=1,6
cmat(i,j)=b(i,j)
enddo
enddo
write(*,*)’ Elasticity matrix = ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmatinv,ia,m,n,b,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end of case1
c begin case2
if(reply .eq. ’2’)then
do i=1,6
do j=1,6
cmat(i,j)=0.
enddo
enddo
cmat(1,1)=13.2e10
write(*,’(a,g12.5,a)’)’Enter c11 [’,cmat(1,1),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,1)=ain(1)
endif
cmat(2,2)=cmat(1,1)
cmat(3,3)=11.5e10
write(*,’(a,g12.5,a)’)’Enter c33 [’,cmat(3,3),’]’
42call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(3,3)=ain(1)
endif
cmat(4,4)=5.2e10
write(*,’(a,g12.5,a)’)’Enter c44 (ANSYS c66) [’,cmat(4,4),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(4,4)=ain(1)
endif
cmat(5,5)=cmat(4,4)
cmat(1,2)=7.1e10
write(*,’(a,g12.5,a)’)’Enter c12 [’,cmat(1,2),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,2)=ain(1)
endif
cmat(2,1)=cmat(1,2)
cmat(1,3)=7.3e10
write(*,’(a,g12.5,a)’)’Enter c13 [’,cmat(1,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmat(1,3)=ain(1)
endif
cmat(3,1)=cmat(1,3)
cmat(2,3)=cmat(1,3)
cmat(3,2)=cmat(2,3)
cmat(6,6)=2.*(cmat(1,1)-cmat(1,2))
write(*,*)’Elasticity matrix c=’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
do i=1,6
do j=1,6
a(i,j)=cmat(i,j)
enddo
enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
43call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
endif
do i=1,6
do j=1,6
cmatinv(i,j)=b(i,j)
enddo
enddo
write(*,*)’Compliance matrix= ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmat,ia,m,n,cmatinv,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end case2
c begin case3
if(reply .eq. ’3’)then
do i=1,6
do j=1,6
cmatinv(i,j)=0.
enddo
enddo
cmatinv(1,1) = 1.65e-11
cmatinv(3,3) = 2.07e-11
cmatinv(4,4) = 4.35e-11
cmatinv(1,2) = -4.78e-12
cmatinv(1,3) = -8.45e-12
write(*,’(a,g12.5,a)’)’Enter s11 [’,cmatinv(1,1),’]’
write(*,*)’Reference for default values:’
write(*,*)’Auld, Appendix A.4, PZT-5H’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,1)=ain(1)
endif
cmatinv(2,2)=cmatinv(1,1)
44write(*,’(a,g12.5,a)’)’Enter s33 [’,cmatinv(3,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(3,3)=ain(1)
endif
write(*,’(a,g12.5,a)’)’Enter s44 (ANSYS c66) [’,cmatinv(4,4),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(4,4)=ain(1)
endif
cmatinv(5,5)=cmatinv(4,4)
write(*,’(a,g12.5,a)’)’Enter s12 [’,cmatinv(1,2),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,2)=ain(1)
endif
cmatinv(2,1)=cmatinv(1,2)
write(*,’(a,g12.5,a)’)’Enter s13 [’,cmatinv(1,3),’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
cmatinv(1,3)=ain(1)
endif
cmatinv(3,1)=cmatinv(1,3)
cmatinv(2,3)=cmatinv(1,3)
cmatinv(3,2)=cmatinv(2,3)
cmatinv(6,6)=2.*(cmatinv(1,1)-cmatinv(1,2))
write(*,*)’Compliance matrix s= cinv =’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmatinv(i,j),j=1,6)
enddo
write(*,*)’Note: In ANSYS components 4,5,6 are permuted’
write(*,*)’ANSYS 44 is 66’
write(*,*)’ANSYS 55 is 44’
write(*,*)’ANSYS 66 is 55’
do i=1,6
do j=1,6
a(i,j)=cmatinv(i,j)
enddo
enddo
ia=6
ib=6
n=6
m=6
inv=1
eps=1.e-10
idet=0
write(*,*)’Computing inverse matrix’
call gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
if(ier .ne. 0)then
write(*,*)’ matrix is singular ’
45endif
do i=1,6
do j=1,6
cmat(i,j)=b(i,j)
enddo
enddo
write(*,*)’ Elasticity matrix = ’
do i=1,6
write(*,’(6(g11.4,1x))’)(cmat(i,j),j=1,6)
enddo
m=6
n=6
l=6
ic=6
call matm(cmatinv,ia,m,n,b,ib,l,cc,ic)
c write(*,*)’Compliance matrix times elasticity matrix= ’
c do i=1,6
c write(*,’(6(g11.4,1x))’)(cc(i,j),j=1,6)
c enddo
endif
c end case3
write(*,*)’Orthotropic constants from compliance matrix:’
e1=1./cmatinv(1,1)
e2=1./cmatinv(2,2)
e3=1./cmatinv(3,3)
write(*,’(a,g15.8)’)’Young’’s Modulus 1= ’,e1
write(*,’(a,g15.8)’)’Young’’s Modulus 2= ’,e2
write(*,’(a,g15.8)’)’Young’’s Modulus 3= ’,e3
prxy=-cmatinv(1,2)*e1
write(*,’(a,g15.8)’)’Major Poisson ratio prxy=’,prxy
prxz=-cmatinv(1,3)*e1
write(*,’(a,g15.8)’)’Major Poisson ratio prxz=’,prxz
pryz=-cmatinv(2,3)*e2
write(*,’(a,g15.8)’)’Major Poisson ratio pryz=’,pryz
g23= 1./cmatinv(4,4)
g13= 1./cmatinv(5,5)
g12= 1./cmatinv(6,6)
write(*,’(a,g15.8)’)’Shear Modulus g23= ’,g23
write(*,’(a,g15.8)’)’Shear Modulus g13= ’,g13
write(*,’(a,g15.8)’)’Shear Modulus g12= ’,g12
write(*,*)
c orthotropic elasticity tensor
do i=1,3
do j=1,3
do k=1,3
e(i,j,k)=0.
do l=1,3
c(i,j,k,l)=0.
enddo
enddo
enddo
46enddo
do i=1,3
do j=1,3
call t2vec(i,j,ij)
do k=1,3
do l=1,3
call t2vec(k,l,kl)
if(k .eq. l)then
c(i,j,k,l)=cmat(ij,kl)
else
c(i,j,k,l)=cmat(ij,kl)/2.
endif
enddo
enddo
enddo
enddo
c
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 81 elasticity tensor components of c? [n]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’n’
endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
do l=1,3
if(c(i,j,k,l) .ne. zero)then
write(*,’(a,i1,i1,i1,i1,a,g15.8)’)’c(’,i,j,k,l,’)=’,c(i,j,k,l)
endif
enddo
enddo
enddo
enddo
endif
write(*,*)’ piezoelectric (stress) tensor coefficients’
do i=1,3
do j=1,3
do k=1,3
e(i,j,k)=0.
do l=1,3
do m=1,3
e(i,j,k) = e(i,j,k) + d(i,l,m)*c(l,m,j,k)
enddo
enddo
c write(*,’(a,i1,i1,i1,a,g15.8)’)’e(’,i,j,k,’)=’,e(i,j,k)
enddo
enddo
enddo
write(*,*)’Do you want to see the nonzero components’
write(*,*)’of the 27 piezoelectric tensor components of e? [n]’
read(*,’(a)’)reply
if(lenstr(reply) .eq. 0)then
reply=’n’
47endif
if(reply .eq. ’y’)then
do i=1,3
do j=1,3
do k=1,3
if(e(i,j,k) .ne. zero)then
write(*,’(a,i1,i1,i1,a,g15.8)’)’e(’,i,j,k,’)=’,e(i,j,k)
endif
enddo
enddo
enddo
endif
c compute the e piezoelectric matrix
do i=1,3
do jk=1,6
call vec2t(jk,j,k)
if(j .eq. k)then
emat(jk,i) = e(i,j,k)
else
emat(jk,i) = 2.*e(i,j,k)
endif
enddo
enddo
c
write(*,*)’Piezoelectric e matrix= ’
do i=1,6
write(*,’(3(g11.4,1x))’)(emat(i,j),j=1,3)
enddo
write(*,*)’Second calculation:’
e13=cmat(1,1)*dmat(3,1)+cmat(1,2)*dmat(3,1)+cmat(1,3)*dmat(3,3)
write(*,*)’e13=’,e13
e33=cmat(1,3)*dmat(3,1)+cmat(1,3)*dmat(3,1)+cmat(3,3)*dmat(3,3)
write(*,*)’e33=’,e33
e42=cmat(4,4)*dmat(1,5)
write(*,*)’e42=’,e42
write(*,*)’Note: ANSYS e rows are permuted’
write(*,*)’ANSYS row 4 is row 6 ’
write(*,*)’ANSYS row 5 is row 4 ’
write(*,*)’ANSYS row 6 is row 5 ’
write(*,*)’ANSYS Piezoelectric e matrix= ’
do i=1,3
write(*,’(3(g11.4,1x))’)(emat(i,j),j=1,3)
enddo
write(*,’(3(g11.4,1x))’)(emat(6,j),j=1,3)
write(*,’(3(g11.4,1x))’)(emat(4,j),j=1,3)
write(*,’(3(g11.4,1x))’)(emat(5,j),j=1,3)
c dielectric constants (farads/meter)
48c dc(1,1)=7.12e-9
c dc(2,2)=7.12e-9
c dc(3,3)=5.84e-9
c density (kilograms/meter^3)
c density=7.5e3
cm=’,’
idmat=3
m=3
n=6
iemat=6
l=3
icc=6
call matm(dmat,idmat,m,n,emat,iemat,l,cc,icc)
write(*,*)’d times e = ’
do i=1,3
write(*,’(3(g11.4,1x))’)(cc(i,j),j=1,3)
enddo
c compute permittivity tensor
eps0=8.85e-12
do i=1,3
do j=1,3
epsilon(i,j)=0.
enddo
enddo
dielx=1300.
write(*,’(a)’)’Enter the unclamped dielectric constant K11’
write(*,’(a,g15.8,a)’)’[’,dielx,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
dielx=ain(1)
endif
epsilon(1,1)=dielx*eps0
diely=dielx
c cs1=’Enter the dielectric constant K22’
c write(*,’(a,a,g15.8,a)’)cs1,’ [’,diely,’]’
c call readr(0,ain,nr)
c if(nr .eq. 1)then
c diely=ain(1)
c endif
epsilon(2,2)=diely*eps0
dielz=1000.
write(*,’(a)’)’Enter the unclamped dielectric constant K33’
write(*,’(a,g15.8,a)’)’[’,dielz,’]’
call readr(0,ain,nr)
if(nr .eq. 1)then
dielz=ain(1)
endif
epsilon(3,3)=dielz*eps0
write(*,*)’Unclamped electric permittivity tensor=’
do i=1,3
49write(*,’(3(g11.4,1x))’)(epsilon(i,j),j=1,3)
enddo
do i=1,3
do j=1,3
clamped(i,j)= epsilon(i,j) - cc(i,j)
enddo
enddo
write(*,*)’Clamped electric permittivity tensor=’
do i=1,3
write(*,’(3(g11.4,1x))’)(clamped(i,j),j=1,3)
enddo
do i=1,3
do j=1,3
sum=0.
do k=1,3
do l=1,3
sum=sum+d(i,k,l)*e(j,k,l)
c if((i .eq. 1) .and. (j .eq. 1))then
c write(*,*)k,l,’ d= ’,d(i,k,l),’ e= ’,e(j,k,l)
c endif
enddo
enddo
epsilon(i,j)=epsilon(i,j)-sum
enddo
enddo
c write(*,*)’Clamped electric permittivity tensor=’
c do i=1,3
c write(*,’(3(g11.4,1x))’)(epsilon(i,j),j=1,3)
c enddo
rho= 7.500000E-03
write(*,*)’Enter density Kg/m^3 [7500]’
call readr(0,ain,nr)
if(nr .ge. 1)then
rho=ain(1)
endif
write(*,’(a)’)’/COM ’
write(*,’(a)’)’/COM Material properties’
c material number
mn=2
write(*,’(a)’)’/COM Young’’s Moduli’
write(*,’(a,i3,a,g15.8)’)’EX,’,mn,cm,e1
write(*,’(a,i3,a,g15.8)’)’EY,’,mn,cm,e2
write(*,’(a,i3,a,g15.8)’)’EZ,’,mn,cm,e3
write(*,’(a)’)’/COM ANSYS "Major" Poisson ratios’
c reference ANSYS Theory 2.1 structural fundamentals
prxy=-cmatinv(1,2)*e1
write(*,’(a,i3,a,g15.8)’)’PRXY,’,mn,cm,prxy
50prxz=-cmatinv(1,3)*e1
write(*,’(a,i3,a,g15.8)’)’PRXZ,’,mn,cm,prxz
pryz=-cmatinv(2,3)*e2
write(*,’(a,i3,a,g15.8)’)’PRYZ,’,mn,cm,pryz
write(*,’(a)’)’/COM Density Kg/m^3’
write(*,’(a,i3,a,g15.8)’)’/COM DENS,’,mn,cm,rho
write(*,’(a)’)’/COM clamped Permitivities’
c epsilon(1,1)=8.7969e-12
c epsilon(2,2)=8.7969e-12
c epsilon(3,3)=8.7969e-12
write(*,’(a,i3,a,g15.8)’)’MP,PERX,’,mn,cm,clamped(1,1)
write(*,’(a,i3,a,g15.8)’)’MP,PERY,’,mn,cm,clamped(2,2)
write(*,’(a,i3,a,g15.8)’)’MP,PERZ,’,mn,cm,clamped(3,3)
write(*,’(a)’)’/COM ANSYS Piezoelectric "e" matrix’
write(*,’(a)’)’TB,PIEZ,3’
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,3,cm,emat(1,3)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,6,cm,emat(2,3)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,9,cm,emat(3,3)
c ansysemat(5,2)=emat(4,2)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,14,cm,emat(4,2)
c ansysemat(6,1)=emat(5,1)
write(*,’(a,i3,a,g15.8)’)’TBDATA,’,16,cm,emat(5,1)
end
c+ readr read a row of floating point numbers
subroutine readr(nf, a, nr)
implicit real*8(a-h,o-z)
c numbers are separated by spaces
c examples of valid numbers are:
c 12.13 34 45e4 4.78e-6 4e2,5.6D-23,10000.d015
c nf=file number, 0 for standard input file
c a=array of returned numbers
c nr=number of values in returned array,
c or 0 for empty or blank line,
c or -1 for end of file on unit nf.
c requires functions val and length
dimension a(*)
character*200 b
character*200 c
character*1 d
c=’ ’
if(nf.eq.0)then
read(*,’(a)’,end=99)b
else
read(nf,’(a)’,end=99)b
endif
nr=0
l=lenstr(b)
if(l.ge.200)then
write(*,*)’ error in readr subroutine ’
write(*,*)’ record is too long ’
endif
do 1 i=1,l
d=b(i:i)
if (d.ne.’ ’) then
51k=lenstr(c)
if (k.gt.0)then
c=c(1:k)//d
else
c=d
endif
endif
if( (d.eq.’ ’).or.(i.eq.l)) then
if (c.ne.’ ’) then
nr=nr+1
call valsub(c,a(nr),ier)
c=’ ’
endif
endif
1 continue
return
99 nr=-1
return
end
c
c+ lenstr nonblank length of string
function lenstr(s)
c length of the substring of s obtained by deleting all
c trailing blanks from s. thus the length of a string
c containing only blanks will be 0.
character s*(*)
lenstr=0
n=len(s)
do 10 i=n,1,-1
if(s(i:i) .ne. ’ ’)then
lenstr=i
return
endif
10 continue
return
end
c+ valsub converts string to floating point number (r*8)
subroutine valsub(s,v,ier)
implicit real*8(a-h,o-z)
c examples of valid strings are: 12.13 34 45e4 4.78e-6 4E2
c the string is checked for valid characters,
c but the string can still be invalid.
c s-string
c v-returned value
c ier- 0 normal
c 1 if invalid character found, v returned 0
c
logical p
character s*(*),c*50,t*50,ch*15
character z*1
data ch/’1234567890+-.eE’/
v=0.
ier=1
l=lenstr(s)
if(l.eq.0)return
52p=.true.
do 10 i=1,l
z=s(i:i)
if((z.eq.’D’).or.(z.eq.’d’))then
s(i:i)=’e’
endif
p=p.and.(index(ch,s(i:i)).ne.0)
10 continue
if(.not.p)return
n=index(s,’.’)
if(n.eq.0)then
n=index(s,’e’)
if(n.eq.0)n=index(s,’E’)
if(n.eq.0)n=index(s,’d’)
if(n.eq.0)n=index(s,’D’)
if(n.eq.0)then
s=s(1:l)//’.’
else
t=s(n:l)
s=s(1:(n-1))//’.’//t
endif
l=l+1
endif
write(c,’(a30)’)s(1:l)
read(c,’(g30.23)’)v
ier=0
return
end
c+ str floating point number to string
subroutine str(x,s)
implicit real*8(a-h,o-z)
character s*25,c*25,b*25,e*25
zero=0.
if(x.eq.zero)then
s=’0’
return
endif
write(c,’(g11.4)’)x
read(c,’(a25)’)b
l=lenstr(b)
do 10 i=1,l
n1=i
if(b(i:i).ne.’ ’)go to 20
10 continue
20 continue
if(b(n1:n1).eq.’0’)n1=n1+1
b=b(n1:l)
l=l+1-n1
k=index(b,’E’)
if(k.gt.0)e=b(k:l)
if(k.gt.0)then
s=b(1:(k-1))
k1=index(b,’E+0’)
if(k1.gt.0)then
e=’E’//b((k1+3):l)
53else
k1=index(b,’E+’)
if(k1.gt.0)e=’E’//b((k1+2):l)
endif
k1=index(b,’E-0’)
if(k1.gt.0)e=’E-’//b((k1+3):l)
l=k-1
else
s=b
endif
j=index(s,’.’)
n2=l
if(j.ne.0)then
do 30 i=1,l
n2=l+1-i
if(s(n2:n2).ne.’0’)go to 40
30 continue
endif
40 continue
s=s(1:n2)
if(s(n2:n2).eq.’.’)then
s=s(1:(n2-1))
n2=n2-1
endif
if(k.gt.0)s=s(1:n2)//e
return
end
c+ t2vec stress tensor index to vector index
subroutine t2vec(i,j,k)
c input:
c i,j
c output:
c k
c so that sigma(k)=sigma(i,j)
integer m(3,3)
data m/1,6,5,6,2,4,5,4,3/
k=m(i,j)
return
end
c+ vec2t stress vector index to tensor index
subroutine vec2t(k,i,j)
c input:
c k
c output:
c i,j
c so that sigma(i,j)=sigma(k)
integer mi(6),mj(6)
data mi/1,2,3,2,1,1/
data mj/1,2,3,3,3,2/
i=mi(k)
j=mj(k)
return
end
c+ gaussr solution of linear equations, inverse, determinant (real) new version
c gaussr2.ftn 4/17/96 modernization of gaussr under construction
54subroutine gaussr(a,ia,b,ib,n,m,inv,eps,idet,det,ier)
c solves the equation a*c=b for c, where a is an n by n matrix
c c and b are n row by m column matrices. c is returned as b.
c algorithm -gaussian elimination with partial pivoting.
c parameters a-n by n matrix containing the coefficients of
c the linear system.
c ia-row dimension of a in defining routine
c e.g. in the routine where a is defined,
c a might be dimensioned as:
c dimension a(nr,nc)
c then ia must be set to nr. we may have n < ia, but
c must not have n > ia, or n*n > nr*nc.
c ia is needed for proper addressing of matrix a.
c fortran stores by column first: a(i,j)=a(i+(j-1)*ia))
c b-n by m matrix containing the m right sides
c of the equations, on entering. on returning, b
c contains the solutions. the inverse of a is
c returned in b when inv=1
c ib-row dimension of b in defining routine
c n-row and column dimension of a.
c m-column dimension of b (usually 1)
c the program changes m to n when inv=1
c inv-the inverse of a is calculated
c and returned in b when inv=1
c b must be large enough to hold the inverse
c eps-each equation is normalized so that the
c coefficients are <= 1 in magnitude.
c when a pivot is less than eps the matrix is
c considered nearly singular, and ier is set to 1
c if a pivot is zero the matrix is singular, and
c ier is set to 2. one may set eps=1.e-5 for
c single precision, and 1.e-12 for double.
c eps does not effect any calculation.
c normalization may also prevent exponent overflow.
c idet-compute determinant only if idet = 1
c determinants are products of n numbers.
c overflow can occur if the elements of the
c matrix have large exponents.
c set idet=0 if the determinant is not needed.
c det-determinant of a.
c ier-return parameter,
c ier=0 normal return
c ier=1 matrix is nearly singular
c ier=2 matrix is singular
c
c warning!! the subroutine changes a and b. if they need to be
c saved, copies must be made before calling the subroutine.
c the subroutine can be converted to different number type
c by uncommenting the appropriate implicit statement.
implicit real*8(a-h,o-z)
c implicit complex(a-h,o-z)
c implicit complex*16(a-h,o-z)
logical cdet
dimension a(ia,*),b(ib,*)
cdet = idet .eq. 1
55zero=0.
ier=0
det=1.
if(m.le.0)m=1
if(inv.eq.1)then
c set b equal to the identity
do i=1,n
do j=1,n
b(i,j)=0.
if(i.eq.j)b(i,j)=1.
enddo
enddo
m=n
endif
c normalize rows
do 20 i=1,n
bigest=a(i,1)
do 16 j=2,n
ab=a(i,j)
if(abs(ab).gt.abs(bigest))bigest=ab
16 continue
if(bigest.eq.zero)then
ier=2
det=0.
return
endif
if(cdet)det=det*bigest
do 18 j=1,n
a(i,j)=a(i,j)/bigest
18 continue
do 19 j=1,m
b(i,j)=b(i,j)/bigest
19 continue
20 continue
j=1
do while( j .lt. n)
kk=j+1
l=j
c find row l with largest pivot
do 32 i=kk,n
if(abs(a(i,j)).gt.abs(a(l,j)))l=i
32 continue
if(abs(a(l,j)).eq.zero)then
ier=2
det=0.
return
endif
if(abs(a(l,j)).le.abs(eps))ier=1
if(l.ne.j)then
c interchange rows l and j
do 37 k=1,n
c=a(l,k)
a(l,k)=a(j,k)
a(j,k)=c
37 continue
56do 39 k=1,m
c=b(l,k)
b(l,k)=b(j,k)
b(j,k)=c
39 continue
if(cdet)det=det*(-1.)
endif
if(cdet)det=det*a(j,j)
c divide row by pivot
c=a(j,j)
do 50 k=j,n
a(j,k)=a(j,k)/c
50 continue
do 55 k=1,m
b(j,k)=b(j,k)/c
55 continue
c add multiple of row j to lower rows
c to eliminate jth coefficients
jj=j+1
do 80 i=jj,n
am=a(i,j)
do 60 k=1,n
a(i,k)=a(i,k)-am*a(j,k)
60 continue
do 70 k=1,m
b(i,k)=b(i,k)-am*b(j,k)
70 continue
80 continue
j=j+1
enddo
am=a(n,n)
if(abs(am).eq.zero)then
ier=2
det=0.
return
endif
if(abs(am).le.abs(eps))ier=1
if(cdet)det=det*am
c a is now in triangular form
c compute nth component of solution
do 90 k=1,m
b(n,k)=b(n,k)/am
90 continue
c back substitute to compute n-i component
c i=1,2,3,...
nn=n-1
do 120 i=1,nn
ni=n-i
do 110 j=1,m
nj=ni+1
do 100 ki=nj,n
b(ni,j)=b(ni,j)-a(ni,ki)*b(ki,j)
100 continue
110 continue
120 continue
57return
end
subroutine matm(a,ia,m,n,b,ib,l,c,ic)
implicit real*8(a-h,o-z)
c arguments
c a-matrix
c ia-row dimension of a in calling program
c m-number of rows of a
c n-number of columns of a
c b-matrix
c ib-row dimension of b in calling program
c l-number of columns of b
c c-product matrix: c=a*b
c ic-row dimension of c in calling program
c
dimension a(ia,*),b(ib,*),c(ic,*)
c c=a*b
do 10 i=1,m
do 10 j=1,l
c(i,j)=0.
do 5 k=1,n
5 c(i,j)=c(i,j)+a(i,k)*b(k,j)
10 continue
return
end
0.22 Location of Piezoelectric Information In
ANSYS Manuals
• Index assignment for stress-strain tensor to six-vector: Theory 2.1.
• Major and Minor Poisson ratios, orthotropic compliance matrices: Theory 2.1.
• Piezoelectric Solid Element SOLID5: Elements manual 4-39.
• Piezoelectric Plane Element PLANE13: Elements manual 4-85.
• Piezoelectric Tetrahedral Coupled Solid Element SOLID98: Elements
manual 4-653.
• Piezoelectric Analysis: Procedures Coupled-Field Analysis 8.1, piezoelectric analysis p8-13 to p8-15.
• Deﬁning matrices: Commands TB and TBDATA, Commands manual.
58• Piezoelectrics Theory: Theory Manual 11.1 to page 11-17.
• Example VM175: Natural Frequency of a Piezoelectric Transducer
...(ANSYS Stuﬀ Notebook). Reference: Boucher D, Lagier M, Maerfeld C, IEEE Transactions on Sonics and Ultrasonics, Volume SU-28,
No. 5 September 1981 pp318-330.
• Example VM176: Frequency Response of Electrical Admittance for a
Piezoelectric Transducer, Reference: Kagawa and Yamabuchi, IEEE
Trans. Sonics and Ultrasonics, V. SU-26, No. 2, March 1979.
• Example PM13, ANSYS Examples Supplement, Piezoelectric Beam
Resonator
0.23 ANSYS Comparison For A Composite
Piezoelectric Transducer: Example VM176.
Reference for the problem: Kagawa and Yamabuchi Finite Element Simulation of a Composite Piezoelectric Transducer IEEE Transactions
on Sonics and Ultrasonics, Vol. SU-26, No 2 march 1979.
The transducer consists of cylindrical disks of aluminum, adhesive, PZT,
adhesive, and aluminum.
PZT material: NEPEC 6
Elasticity matrix
c = 10
10








12.8 6.8 6.6 0  0  0
6.8 12.8 6.6 0  0  0
6.6 6.6 11  0  0  0
0 0 0 2.1 0  0
0 0 0 0 2.1 0
0  0  0 0 0 2.1








59Piezoelectric matrix (C oulomb/M  eter
2
)
e =






0 0 −6.1
0 0 −6.1
0 0  15.7
0 0  0
0 0  0
0 0  0






Dielectric Matrix (F arad/M eter)
 = 10−9



8.7969 0 0
0 8.7969 0
0 0 8.7969



See ANSYS notebook for more details.
600.24 Piezoelectricity Bibliography
[1]Agarwal Bhagwan D, Broutman Lawrence, Analysis And Performance
Of Fiber Composites.
[2]Ahmed F R, Crystallographic Computing, Proceedings International
Summer School, Linda Hall Library.
[3]Allik H, Hughes J R, Finite Element Method For Piezoelectric Vibration, International Journal of Numerical Methods For Engineering, V.
2, pp 151-157, 1970, (This is the formulation of the problem that is implemented in ANSYS.).
[4]Ashton J E, Whitney J M, Theory Of Laminated Plates, Linda Hall
Library, 1979.
[5]Ashton T E, Halpin T C, Petit P H, Primer On Composite Materials
Analysis, Technomic, Stamford Conn., 1969.
[6]Auld B A, Acoustic Fields And Waves In Solids, Volumes I and II.
[7]Billington E W, Introduction To The Mechanics And Physics Of
Solids, AlliedSignal Library, QA808 .2, 1986, A Mathematical Book. Von
Mises Yield Criterion. Von Mises Stress.
[8]Bottom Virgil E, The Theory And Design Of Quartz Crystal Units,
2nd Ed., 1974, (Self Published), Linda Hall, (Piezoelectric Transducers.).
[9]Boucher D, Lagier M, Maerfeld, Computation Of The Vibrational
Modes For Piezoelectric Array Transducers Using A Mixed Finite
Element-Perturbation Method, IEEE Transactions on Sonics And Ultrasonics, V. Su28, No. 8, Sept 81.
[10]Brebbia, Dewilde,Blain, Computer Aided Design In Composite Material Technology, TA 418.9 .c6 I585, 1988, Linda Hall.
[11]Buchanon Relva C, Editor, Ceramic Materials For Electronics, 2nd
Ed., 1991, AlliedSignal Library, (Piezoelectricity.).
[12]Buerger, Martin, Introduction To Crystal Geometry.
[13]Cady Walter Guyton, Piezoelectricity, Volumes I and II, Dover, 1964,
61KCMO Library.
[14]Cady Walter Guyton, The Piezoelectric Resonator, In: Physical
Acoustics, R. Bruce Lindsey Editor, (Benchmark Papers In Acoustics).
[15]Clough R W, Penzien J, Dynamics Of Structures, McGraw-Hill 1972,
AlliedSignal Library, (Rayleigh Damping, Piezoelectricity.).
[16]Cracknell A P, Ultrasonics, Wykeham Publications, 1980 AlliedSignal
Library.
[17]Crandall B C, Lewis James, Editors, Nanotechnology, Article: John
Foster, p17, Scanning Tunneling Microscope (STM), (Uses Piezoelectric Control,), 1992, MIT Press.
[18]Cullity, Elements Of X-Ray Diﬀraction, Addison-Wesley, 1956.
[19]Cummings James R, Predicting The Eﬀect Of Thickness Of Piezoelectric Plates On The Performance Of The Ultrasonic Traveling
Wave Motor, U. Mo. (Rolla,), May 1993.
[20]Cummings James R, Stutts Daniel S, Dharani L R, Eﬀects Of Material And Geometric Characteristics On The Performance Of A
Traveling-wave Motor, U. MO. Rolla Report, 11-9-93.
[21]Cummings Robert James, Analysis Of The Performance And Contact Behavior Of The Piezoelectric Traveling-wave Motor, Masters
Theses, University Of Missouri-Rolla MRTC Manufacturing Research And
Training Center, 1994.
[22]Emery James D, Elasticity, Unpublished Document.
[23]Emery James D, Electomagnetic Theory, Unpublished Document.
[24]Emery James D, Transducer Equivalent Circuits, Unpublished Document.
[25]Emery James D, Piezo Notes I, Unpublished Document.
[26]Emery James D, Piezo Notes II, Unpublished Document.
[27]Flynn Anita, Torque Production In Ultrasonic Motors, MIT Arti-
62ﬁcial Intelligence Laboratory, January 9, 1994, (Piezoelectric Motor.).
[28]Flynn Anita M, Piezoelectric Ultrasonic Micromotors, PhD Thesis,
MIT, June 1995.
[29]Fulton Langdon H, Techniques And Applications Of The Gas Squeeze
Film Bearing, Masters Thesis, Moore School, U. Pennsylvania, Dec 1967,
(Piezoelectric Application. Now At Martin-Marietta Communications Systems, Camden NJ).
[30]Gooberman G L, Ultrasonics, 1968, Acoustic Impedance, Transmission,
Lines, Transducers, Equivilent Circuit As Transmission Line, Q Of Transducer.
[31]Goree, Pagano-Goree Free Edge Stress, (Discussion, Goree: Singular
Elements In Abaqus Are Needed, Disagreement With Whitney.).
[32]Hagedorn P, Wallashek J, Traveling Wave Ultrasonic Motors, Part
I: Working Principle And The Mathematical Modeling Of The Stator, J. Of Sound And Vibration, V. 155, No. 1, pp 31-46, 1992.
[33]Hagood N, Chung W,Von Flotow A, Modeling Of Piezoelectric Actuator Dynamics For Active Structural Control, J. Intell.  Mat.  Syst.
And Structures, V. 1, pp327-354, July 1990, (Hamilton’s Principal Modiﬁed
For General Electromechanical Systems.).
[34]Hagood Nesbit, Mcfarland Andrew, Modeling Of A Piezoelectric Rotary Ultrasonic Motor, IEEE Trans. Ultrasonics Ferroelectrics And Frequency Control, V. 42, No. 2, March 1995, Nesbit was Anita Flynn Collegue
at MIT. Rayleigh-Ritz Method, which includes, Electrical Energy.
[35]Hauptman, Herbert, Crystal Structure Determination, Linda Hall,
QD911.h33, 1972, Mathematical Biophysics Department, Medical Foundation Of Buﬀalo, Cosine Semiinvariants.
[36]Heising Raymond A, Quartz Crystals For Electrical Circuits: Their
Design And Manufacture, Van Nostrand, 1946, Reprinted by Electronic
Industries Association, 2001 Eye St, NW, Washington DC 20006, 1978, 1982,
(Piezoelectric Equivalent Circuits.).
[37]Hirose, Takahashi, Uchino, Aoyagi, Tomakawa, Electromechanical Prop-
63erties Of PbZro3-pbzrtio3-pb(mn(1/3)sb(2/3)o3 Ceramics Under
Vibration Level Change, Mat. Res. Soc. Symp. Proc., V. 360, 1995,
Materials Research Society.
[38]Hirose, Takahashi, Uchino, Aoyagi, Tomakawa, Measuring Methods
For High-Power Characteristics Of Piezoelectric Materials, Mat.
Res. Soc. Symp. Proc., V. 360, 1995, Materials Research Society.
[39]Hossack John A, Hayward Gordon, Finite-Element Analysis of 1-3
Composite Transducers, IEEE Transactions on Ultrasonics, Ferroelectrics,
and Frequency Control, V. 38, N0. 6, November 1991, (Uses ANSYS.).
[40]Inman Daniel, Engineering Vibration, MatLab Disk, 1994, Chapter:
Experimental Modal Analysis.
[41]Institute of Radio Engineers, IRE Standards on Piezoelectric Crystals: Determination of the Elastic, Piezoelectric, and Dielectric
Constants, The Electromechanical Coupling Factor, 1958, Proc. IRE,
V. 46, pp764-78.
[42]Jaﬀe Bernard, Cook William, Jaﬀe Hans, Piezoelectric Ceramics, Academic Press, 1971, Linda Hall, TK7871.15.C4 J3, Deﬁnitions, Symmetry,
Constants, Measurements, Hysteresis.
[43]Johnson K L, Contact Mechanics, Linda Hall Library.
[44]Jones, Acoustic And Electromagnetic Waves, Oxford.
[45]Kagawa, Yamabuchi, Finite Element Simulation of a Composite
Piezoelectric Transducer, IEEE Transactions on Sonics and, Ultrasonics,
V. SU-26, No. 2 march 1979, (ANSYS example VM176.).
[46]Kikuchi And Oden, Contact Problems In Elasticity: A Study Of
Variational Inequalities And Finite Element Methods, Linda Hall,
QA931.k47 1988, (SIAM), Society For Industrial And Applied Mathematics.
[47]Kittle, Solid State Physics, (Piezoelectric Equations.).
[48]Koehler D R, Radial Vibration Modes In Thin Plates Of BaT iO3
And LiNBO3, Memo, Sandia National Laboratory, October 16, 1989.
64[49]Koiter W T , Miknailor G K, Theory Of Shells.
[50]Kolsky, Stress Waves In Solids, Dover.
[51]Kompaneyts A S, Theoretical Physics, Dover, 1961.
[52]Lass Victor, Vector And Tensor Analysis, (Relation Between Mechanical, Electric, And Magnetic Energy.).
[53]Lazan B J, Damping Of Materials And Members In Structural
Mechanics, Oxford 1968.
[54]Lerch, R, Simulation Of Piezoelectric Devices By Two And Three
Dimensional Finite Elements, IEEE Trans. On Ultrasonics Ferroelectrics
And Freq. Control, V. 37, No. 2, May 1990.
[55]Levich B G, Theoretical Physics, Volumes I, Ii, Iii, Wiley-Interscience
1971.
[56]Lur’e A L, Three-dimensional Problems Of Elasticity, (Basis Of
James Cummings Rolla Analysis Of Tooth Contact In Piezo Motor).
[57]Mansﬁeld, E H, The Bending And Stretching Of Plates, (The Small
Deﬂection Theory Of Plates Is Analogous To Beam Small, Deﬂection Theory. Gaussian And Mean Curvature, Where Curvature Is Taken To, Be
Approximated By The 2nd Order Derivatives. Twist Is The Cross Derivative, Transformation Matrix For Curvatures Handled As In JRC (above).
[58]Marsden Jerrold, The Mathematical Theory Of Elasticity, (Contains Much Diﬀerential Geometry).
[59]Martin, Transformation Geometry, (Symmetry Groups, Transformations).
[60]Mason Warren P, Piezoelectric Crystals and Their Application to
Ultrasonics, Van Nostrand, 1950.
[61]Mason Warren P, Physical Acoustics, 1964 Academic, AlliedSignal Library, QC225.m3, (Piezoelectric Actuators.).
[62]Mcfarland, Smith, Bernhart, Analysis Of Plates, Linda Hall Library.
65[63]McLachlan N W, Theory Of Vibrations, Circular Membrane, Circular
Plate, Bessel Functions p 129-144.
[64]Meirovitch Leonard, Elements Of Vibration Anlysis, Mcgraw-Hill,
1975.
[65]Meirovitch, Leonard, Analytical Methods Of Vibration, Vibration
Of Plates p 179, Eigenvalues, Rayleigh Quotient, The Biharmonic Equation.
[66]Mentesana C P, Optical Stronglink Development And Piezoelectric Motor Evaluation, AlliedSignal KCD, Kcp-613-4573, May 1992.
[67]Mercer C D, Reddy B D, Eve R A, Finite Element Method For Piezoelectric Media, UCT/CSIR Applied Mechanics Research Unit, Technical
Report No. 92, April, 1987.
[68]Merhaut Josef, Theory Of Electroacoustics, 1976, Czech Book, Johnson County Library, (Piezoelectric Transducers.).
[69]Mikhlin, S.g, Variational Methods In Mathematical Physics, Linda
Hall Library, 1964.
[70]Mizushima Masataka, Theoretical Physics, John Wiley, 1972, AlliedSignal, (From Classical Mechanics To General Relativity, Relativistic Quantum Mechanics And Group Theory.).
[71]Moore And Christopher Davis And Coplan, Building Scientiﬁc Apparatus, (Experimental Techniques.).
[72]Morgan Matroc Incorporated, Guide To Modern Piezoelectric Ceramics, Catalog, Vernitron Division.
[73]Morse, Philip M, Vibration And Sound.
[74]Musikant Solomon, What Every Engineer Should Know About Ceramics, 1991, Marcel Dekker, AlliedSignal Book.
[75]Newnham R E, Trolier-McKinstry S, Giniewicz J R, Piezoelectric, Pyroelectric And Ferroic Crystals, J. Material Education, 15, pp 189-223
(1993), (Reprint: Penn State MRL Annual Report To ONR 1994, Appendix
661).
[76]Nye J F, Physical Properties Of Crystals, QD931 .N8, 1957, 1987.
[77]Ostergard D, Pawlak T, Three-Dimensional Finite Elements For
Analyzing Piezoelectric Structures, Proceedings IEEE Ultrasonics Symp.,
Williamsburg, Va, 1986, pp 639-642.
[78]Paganno N J, Pipes R B, Interlaminar Stresses In Composite Laminates Under Uniform Axial Extension, Journal Of Composite Materials, V. 4, 1970, pp 538-548.
[79]Paganno N J, Pipes R B, The Inﬂuence Of Stacking Sequence On
Laminate Strength, Journal Of Composite Materials, V. 5, 1971, pp 50-57.
[80]Panasonic, Ultrasonic Motor, Panasonic Electric Motor Division, Matsushita Electric Industrial Co. Ltd., (Equivalent Circuit.).
[81]Parker Sybyl, Solid State Physics Sourcebook.
[82]Press, Flannery, Teukolsky, Numerical Recipes in Fortran, (Subroutines For Computing Bessel Functions.).
[83]Pressly Robert B, Mentesana Charles P, Piezoelectric Motor Development At AlliedSignal Kansas City Division, KCP-613-5515, November, 1994.
[84]Puppo A. H. Evenson H A, Journal Of Composite Materials, V. 4,
1970, pp 201-214.
[85]Raju P. Narayana, Vibration Of Annular Plates, Journal Of The
Aeronautical Society Of India, V. 14, No. 2, May, 1962, (Boundary Conditions For Maple Motor Code).
[86]Ramachandran G N , Srinivasan R, Fourier Methods In Crystallography, Linda Hall Library, 1970.
[87]Rayleigh Lord, The Theory Of Sound, Volumes I and II, Dover, (Vibrations Of Plates: Chapter 10, p 352, Origin of Rayleigh Damping Concept).
67[88]Roiguiskis, Bansevicius, Barauskas, Kulviets, Vibromotors For Precision Microrobots, AlliedSignal Library, (Many Kinds Of Piezoelectric
Motor Designs.).
[89]Saigoh Hiroaki, Kawasaki Mayumi, Maruko Nobuhiro, Kanayama Kouichi,
Multilayer Piezoelectric Motor Using The First Longitudinal And
Second Bending Vibrations, Japanese J Applied Physics, V. 34, 1995,
pp2760-2764.
[90]Sands Donald E, Introduction To Crystallography, 1975, Dover.
[91]Sashida Toshiiku, Kenjo Takashi, An Introduction To Ultrasonic
Motors, Oxford U. Press 1993, Piezoelectricity.
[92]Sawyer C B, Tower C H, Phys. Rev. 35, (1930) pp269-73, (Technique For Displaying Ferroelectric-Piezoelectric Hysteresis Curves.).
[93]Sherrit S, Gauther N, Wiederick H D, Mukherjee B K, Accurate Evaluation Of The Real And Imaginary Material Constants For A Piezoelectric Resonator In The Radial Mode, Ferroelectrics, 1991, V. 119,
pp17-32.
[94]Shields John Potter, Basic Piezoelectricity, 1966, Sams, QC595 .S4,
Linda Hall Library.
[95]Smits Jan G, Iterative Method For The Accurate Determination
Of The Real And Imaginary Parts Of The Materials Coeﬃcients
Of Piezoelectric Ceramics, Trans. Sonics Ultraacoustics. V. Su-23, No.
6, Nov. 1976, pp393-402.
[96]Soedel Werner, Vibrations Of Shells And Plates, 1981. Second Edition 1993, University Illinois At Chicago.
[97]Sokolnikoﬀ I S, Mathematical Theory Of Elasticity, 1956, (Biharmonic Equation p77, Two Dimensional Problems Ch. 5, Finite Diﬀerences,
Relaxation Methods, pp442-454.).
[98]Southwell R V, Theory Of Elasticity.
[99]Stoker, Nonlinear Vibrations.
68[100]Stutts Daniel, Loss Factors, FAX: Feb 23, 1996, Various Materials
Loss Coeﬁcients vs Young Modulus, Comparison Of Damping for a Variety
of Metals, Magnetoeleastic Damping In Ferromagnetic Metals, Internal Friction Vs Strain Amplitude, Experimental Setup With Pucot Drive Crystal
And Specimin.).
[101]Stutts Daniel, Modal Analysis Of The Free-Free Rod With Damping, FAX: April, 9 1996, (See Piezo96 Notebook.).
[102]Stutts Daniel, A Simple Example Of The Relationship Between
Electromechanical Coupling And Measured Impedance, Unpublished
Document, University of Missouri at Rolla, Mechanical Engineering Department, March 1, 1995.
[103]Stutts Daniel, The Free-free Annular Plate, Maple Program And
Document, U. Mo. Rolla Mech. Eng. Dept, January,1994.
[104]Stutts Daniel S, Friend James R (Formerly Cummings), Quarterly Report: Piezoelectric Motor Project, Oct. 2 1995 (tooth Motion, Asymmetry, Maple Conference, Contact Mech, Web Page).
[105]Takahashi Sadayuki, Herose Seiji, Uchino Kenji, Oh Ki-Young, Electromechanical Characteristics Of Lead-Zirconate-Titanate Ceramics
Under Vibration Level Change, 1995 IEEE.
[106]Takuro Ikeda, Fundamentals Of Piezoelectricity,, (Electromechanical Coupling Coeﬃcient, X-cut, Y-cut Plate, Equivalent Circuit. p84.).
[107]Termin, Radio Engineering, (Equivalent Circuit Of Piezo Oscillator).
[108]Tiersten H F, Linear Piezoelectric Plate Vibration, Plenum Press,
1969, pp33-99, (Fundamental Equations And Equivalent Circuit.).
[109]Timoshenko And Goodier, Theory Of Elasticity, Third Edition.
[110]Timoshenko, And Woinowsky-Krieger, Theory Of Plates And Shells,
2nd Ed.
[111]Tomikawa Y, Ueha S, Ultrasonic Motors: Theory And Application, Oxford 1993.
69[112]Tuinenga Paul, Spice, A Guide To Circuit Simulation And Analysis Using Pspice.
[113]Tzou H S, Piezoelectric Shells, Kluwer Academic Publishers, 1993.
[114]Urick Robert J, Principles Of Underwater Sound, 2nd Ed., McGrawHill, 1975, (Piezoelectric Transducers.).
[115]Van Randeraat J, Setterington R E, Piezoelectric Ceramics, London,
Mullard, Linda Hall Library, 1974.
[116]Voyiadjis, Karamanlidis, Editors, Advances In The Theory Of Plates
And Shells, Linda Hall Library.
[117]Wallashek Jorg, Piezoelectric Ultrasonic Motors, Conference On
Mechatronics And Robotics, Duisburg/Moers, Germany, Sept 27-29, 1993
pp107-125.
[118]Webb S, The Physics Of Medical Imaging, 1988, Ultrasound, NMR,
CT, Tomography.
[119]Weinberger, Hans, Partial Diﬀerential Equations, (Sturm-Liouville
Problems, pp182-187, Circular Membrane, Bessel Functions.).
[120]Weinberger, H F, Variational Methods In Boundary Value Problems, 1961.
[121]Wert, Thomson, Physics Of Solids.
[122]Whitney J M, Browning C E, Journal Of Composite Materials, V.
6, April 1972, pp 50-57.
[123]Whitney J M, Stansbarger D L, Howell H B, Journal Of Composite
Materials, V. 5, April 1971, pp24-34.
[124]Whitney James M, Structural Analysis Of Laminated Isotropic
Plates, 1987, Technomic, (Includes Lampcal Software.).
[125]Whitney, Pipes, Experimental Mechanics Of Fiber Reinforced
Composite Materials, 1982, T-C Publishers.
[126]Willmore, T J, Tensors As Tensor Products, Linda Hall Library,
701959.
[127]Woolett Ralph S, Theory Of Piezoelectric Flexural Disk Transducers With Applications To Underwater Sound, USL Report No.
490, U. S. Navy Underwater Sound Laboratory, Fort Trumbull, New London
Conn, 1960, S-F001 03 04-1.
[128]Wooster, W A, Breton A, Experimental Crystal Physics, KCMO
Library.
[129]Yang, Variational Formulation Of Piezoelectric Elements.
[130]Ziman J M, Principles Of The Theory Of Solids, Cambridge, 1964,
(Second Edition 1979).
71

Create Shorturl Create a shorter url that redirects to your paste?

Private Private paste aren't shown in recent listings.

Delete After When should we delete your paste?

Bugabuse.Net Pastebin

Piezoelectric Crystals

Reply to "Piezoelectric Crystals"