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 Coefficients . . . . 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 find the clamped permittivity matrix from the unclamped matrix. These conversions are required for finite 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 first 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 field), it attains a permanently polarized state. For a small variation in the electric field, it behaves approximately linearly near that state. The linear piezoelectric equations are written in two different 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 field 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 coefficients). 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 effect, when a piezoelectric material is put under a stress, the material becomes electrically polarized and surface charges appear. The direct piezoelectric effect is Pi = dijkσj k, which is obtained from the second equation, when the external electric field is zero. 3In the converse piezoelectric effect, when a piezoelectric material is put in an external electric field (voltages applied to electrodes), the material experiences a strain. The converse piezoelectric effect is εij = Ekdkij , which is obtained from the first equation, when there are no external forces and the stress is zero. The fact that the coefficients for both the direct effect and the converse effect 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 define 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 define 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 defined 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 finite element program ANSYS. 4h 1 2 3 1 1 4 6 2 4 2 5 3 6 5 3 Warning! These differences 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 field 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 coefficients). The ij are the components of the permittivity tensor. These equations are suitable for finite 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 coefficients and the strain coefficients are related as follows: eijk = cj klmdilm. Notation varies greatly in the field of piezoelectricity. In a following section we will summarize the notation used in physics, and in two finite element programs ABAQUS and ANSYS. 50.3 Piezoelectric Ceramics Constants and One Dimensional Equations The principal references for this section are: Jaffe, Cook and Jaffe, 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 definition 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 first set (strain form). Solving the first equation for the stress, we have σ = cε − cdE = cε − eE. Making this substitution for σ in the alternate form of the second equation of the first 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 coefficient e = cd. Here ε is the permittivity at constant stress (unclamped). 6The coefficients may be defined 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 defined by g = (− ∂E ∂σ )ε = ( ∂ε ∂D )σ. We have E = 1 (D − dσ) so g = −( ∂E ∂σ )ε = d . The constant e is defined 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 reflection 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 find that the proper signs can be maintained only if some of the elastic coefficients are zero (see pp62-64 Sokolnikoff Mathematical Theory of Elasticity). Carrying out rotation transformations about the poling axis (3 direction) we find further conditions on the elasticity coefficients. We have only the following independent elastic coefficients (Jaffe et al, p20) c11, c33, c44, c12, c13. 7We have c22 = c11, c55 = c44, c66 = 2(c11 − c12), c23 = c13. Coefficients which are not in the upper 3-dimensional submatrix or on the main diagonal are zero. The independent piezoelectric coefficients are d33, d31, d15. We have d32 = d31, d24 = d15. All other coefficients are zero. The independent dielectric coefficients are K11, K33 We have K22 = K11. All other coefficients are zero. The permittivity is ij = 0Kij . The electromechanical coupling constant is defined 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 (Jaffe 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 first 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 defined 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 defined 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 five independent c coefficients. 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 defined 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 finite 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 coefficients are written as γij . They differ from the εij , in that the shear coefficients are doubled in value. Note: the physics notation for eijk is defined in equation 8.40 of Auld. Nye does not introduce the piezoelectric e coefficients. The Cady book is very old and one can obtain piezoelectric constant notation only by implication. Note: The ANSYS definition of the e matrix and tensor is the transpose of the physics and ABAQUS definition. 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 definition and the physics e definition 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 first 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 field 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 differential in the thermodynamic sense. Therefore redefining 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 files: elastic.tex, elastic.ps). An orthotropic problem is one in which the material properties are symmetric with respect to 3 mutually orthogonal planes (See Sokolnikoff 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 finds 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 field 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 Jaffe, Cook, Jaffe 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 Coefficients This equality may be shown by applying the first and second laws of thermodynamics. Let the stress coefficients σij , the electric field components Ei and the temperature T be the state variables. Then write differential 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 differential of Φ using the last two equations, and also write it in the general partial derivative form. Then equate the partial derivatives to coefficients of the equations. For example ∂Φ ∂Ei = −Di . Differentiate to equate derivatives of eij , Di , or S to second order partial derivatives of Φ. This shows for example that the coefficients in the direct piezoelectric effect are equal to the coefficients of the converse piezoelectric effect. 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 coefficient) 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 coefficients. 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 finite element nodal potentials φn, we compute derivatives of the shape functions, and then the gradient of the potential. Then we obtain the electric field 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 fit it to a function Q0 cos(ωt + ψ), and thus to the real part of the exponential Q0 exp(j(ωt + ψ)). Differentiating 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 finite 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 field 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 field 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 field 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 field in the ferroelectric capacitor C1. The charge on the plates of C1 is proportional to the D field (actually the D field 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 field 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 field 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 field curve for a piezoelectric ceramic is a butterfly 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 field is removed, the material almost returns to an unpolarized state. It requires a bias field. 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 field 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 field 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 Jaffe: Piezoelectric Ceramics, Electrostriction is mentioned on pages 15 and 78. He relates electrostriction to domain reversal. Jaffe says that the effect is significant only for ferroelectrics above the Curie point. Cady: Piezoelectricity defines electrostriction to be that effect that is proportional to E 2 . See pages 4, 198, 199, 614. Cady states that electrostriction is only significant for fields 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 effect 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 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 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 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 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 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 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 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 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. • Defining 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 Stuff 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 Diffraction, Addison-Wesley, 1956. [19]Cummings James R, Predicting The Effect 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, Effects 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- 62ficial 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 Modified 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 Buffalo, 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]Jaffe Bernard, Cook William, Jaffe Hans, Piezoelectric Ceramics, Academic Press, 1971, Linda Hall, TK7871.15.C4 J3, Definitions, 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]Mansfield, E H, The Bending And Stretching Of Plates, (The Small Deflection Theory Of Plates Is Analogous To Beam Small, Deflection 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 Differential 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 Scientific 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 Influence 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 Coefficients 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]Sokolnikoff I S, Mathematical Theory Of Elasticity, 1956, (Biharmonic Equation p77, Two Dimensional Problems Ch. 5, Finite Differences, 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 Coeficients 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 Coefficient, 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 Differential 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