In linear piezoelectricity the equations of elasticity are coupled to the charge equation of electrostatics by means of piezoelectric constants (IEEE Standard on Piezoelectricity([90])):

(10–51)

(10–52)

or equivalently

(10–53)

where:

Equation 10–51 and Equation 10–52 are the usual constitutive equations for structural and electrical fields, respectively, except for the coupling terms involving the piezoelectric matrix [e].

The elasticity matrix [c] is the usual [D] matrix described in Structural Fundamentals (input using the MP commands). It can also be input directly in uninverted form [c] or in inverted form [c]^-1 as a general anisotropic symmetric matrix (input using TB,ANEL):

(10–54)

The piezoelectric stress matrix [e] (input using TB,PIEZ with TBOPT = 0) relates the electric field vector {E} in the order X, Y, Z to the stress vector {T} in the order X, Y, Z, XY, YZ, XZ and is of the form:

(10–55)

The piezoelectric matrix can also be input as a piezoelectric strain matrix [d] (input using TB,PIEZ with TBOPT = 1). Mechanical APDL will automatically convert the piezoelectric strain matrix [d] to a piezoelectric stress matrix [e] using the elasticity matrix [c] at the first defined temperature:

(10–56)

The orthotropic dielectric matrix [ε^S] uses the electrical permittivities (input as PERX, PERY and PERZ on the MP commands) and is of the form:

(10–57)

The anisotropic dielectric matrix at constant strain [ε^S] (input using TB,DPER,,,,0 command) is used by PLANE222, PLANE223, SOLID225, SOLID226, and SOLID227 and is of the form:

(10–58)

The dielectric matrix can also be input as a dielectric permittivity matrix at constant stress [ε^T] (input using TB,DPER,,,,1). The program will automatically convert the dielectric matrix at constant stress to a dielectric matrix at constant strain:

(10–59)

where:

For the LINK228 element, the anisotropic matrices defined by the TB command (TB,ANEL; TB,DPER) are excluded. Only the axial component of the piezoelectric matrix [e₁₁] (input using TB,PIEZ with TBOPT = 0) is considered. Only the axial component of the electrical permittivity (input as MP,PERX) is used to form the dielectric matrix.

Structural and electrical losses can be introduced in dynamic piezoelectric analyses using various material properties. This section describes losses that are included in the piezoelectric constitutive equations of elements PLANE222, PLANE223, SOLID225, SOLID226, and SOLID227 and result in heat generation.

Viscous structural damping can be included in a transient or harmonic piezoelectric analysis using a stiffness matrix multiplier or in the form of a viscosity matrix:

(10–60)

where:

βm = stiffness matrix multiplier (input as BETD on MP or SDAMP on TB)

[η] = anisotropic viscosity matrix (input using TB,AVIS with TBOPT = 0)

= strain rate vector

The stiffness matrix multiplier β^m is described in Damping Matrices. The symmetric viscosity matrix [η] is of the form:

(10–61)

Anisotropic damping can also be in the form of fluency matrix [ζ] (input using TB,AVIS with TBOPT = 1):

(10–62)

where:

The program will convert the fluency matrix [ζ] to viscosity matrix [η] using the elastic stiffness matrix [c^E]:

(10–63)

Hysteretic structural damping can be included in harmonic and damped modal piezoelectric analyses using several structural damping coefficients or in the form of a loss tangent matrix:

(10–64)

where:

The damping coefficients β^ξ, m, and g are described in Damping Matrices. The symmetric elastic loss tangent matrix [ψ] (input using TB,ELST) is of the form:

(10–65)

In a harmonic and modal piezoelectric analysis, the contributions from all types of viscous and hysteretic damping can be summarized as follows:

(10–66)

Note that constant structural damping coefficients input using BETAD and DMPSTR commands are not included in the above piezoelectric constitutive relations since they are not directly associated with any loss in the material.

Electrical losses can be included in a harmonic piezoelectric analysis by specifying electrical resistivity, electrical loss tangent coefficient, or a dielectric loss tangent matrix:

(10–67)

where:

j = imaginary unit

electrical conductivity matrix

ρxx, ρyy, ρzz = coefficients of electrical resistivity (input as RSVX (Y, Z) on MP command)

tanδ = electric loss tangent coefficient (input as LSST on MP command)

[φ] = dielectric loss tangent matrix (input using TB,DLST with TBOPT = 0)

Electrical resistivity and loss tangent are described in Quasistatic Electric Analysis. The symmetric dielectric loss tangent matrix [φ] (input using TB,DLST) is of the form:

(10–68)

The electric loss tangent coefficient (MP,LSST) and the dielectric loss tangent matrix (TB,DLST) can also be included in a piezoelectric damped modal analysis (Equation 10–72).

Electrical resistivity (MP,RSVX/Y/Z) is the only type of electric loss that can be included in a transient piezoelectric analysis. To account for this loss, the transient piezoelectric analysis must be current-based (Equation 10–73).

For the LINK228 element, these anisotropic material properties are not included in the structural damping:

For the same element, the electrical losses include the axial component of the electrical resistivity coefficients (input as MP,RSVX) and the electric loss tangent coefficient (input as MP,LSST).

Different types of losses and their availability in dynamic piezoelectric analyses are summarized in the table below.

After the application of the variational principle and finite element discretization (Allik([82])), the coupled finite element matrix equation derived for a one element model is:

(10–69)

where:

[K] = element stiffness matrix (defined by Equation 2–58)

[M] = element mass matrix (defined by Equation 2–58)

[C] = element structural damping matrix (discussed in Damping Matrices)

[Cm] = element anisotropic structural damping matrix

{F} = vector of nodal and surface forces (defined by Equation 2–56 and Equation 2–58)

[Kd] = element dielectric permittivity coefficient matrix ([Kvs] in Equation 5–117 or [Kvh] in Equation 5–116)

{L} = vector of nodal, surface, and body charges (defined by Equation 5–117)

[B] = strain-displacement matrix (see Equation 2–44)

[Ce] = element electric damping matrix

{εth} = thermal strain vector (as defined by equation Equation 2–3)

{N} = element shape functions

---

Note:  In a strongly coupled thermo-piezoelectric analysis (see Equation 10–24), the electric potential and temperature degrees of freedom are coupled by:

where:

{α} = vector of coefficient of thermal expansion.

---

In a harmonic analysis, the structural [C^m] and electric [C^e] damping matrices are formed as follows:

(10–70)

(10–71)

where:

= element viscous structural damping matrix

= element hysteretic structural damping

= element resistive damping matrix

= element hysteretic electric damping

In a modal analysis, when the frequency ω is not known initially, the hysteretic damping terms are included in the finite element matrix equation as imaginary stiffness:

(10–72)

A current-based piezoelectric analysis is obtained by differentiating the electric charge equation of the piezoelectric coupled system (Equation 10–69) with respect to time:

(10–73)

where:

The current-based piezoelectric analysis can be harmonic or transient. Resistive losses in a piezoelectric material can be simulated using a transient current-based analysis.

This section describes the results calculation in a piezoelectric analysis using elements PLANE222, PLANE223, SOLID225, SOLID226, and SOLID227.

The total stress {T} (output as S) and total electric flux density {D} (output as D) are computed from elastic strain {S} and electric field {E} as shown in Equation 10–51 and Equation 10–52. In dynamic piezoelectric analyses, structural and electrical damping is included in the values of stress and electric flux according to Equation 10–60, Equation 10–64, and Equation 10–67. Both the stress {T} and electric flux density {D} are evaluated at the integration point locations and then moved or extrapolated to the corner nodes for output.

In addition to stress {T} and electric flux {D}, piezoelectric results include the total electric current density {J} (output as JS) calculated as a sum of the total (conduction and displacement) electric current {J^e} (see Quasistatic Electric Analysis) and the piezoelectric displacement current {J^p}:

(10–74)

The total electric current density {J} is calculated at the element centroidal location and output in the global Cartesian coordinate system. It can be used as a source current density in a subsequent magnetic analysis (LDREAD,JS).

In a current-based piezoelectric analysis, the conduction current density {J^c} = [σ]{E} is output (as JC) instead of the electric flux {D}.

In static and transient piezoelectric analyses, the element instantaneous energies are calculated as:

(10–75)

(10–76)

(10–77)

where:

= stored elastic strain energy (output as an NMISC element item UE).

= dielectric energy (output as an NMISC element item UD)

= electromechanical or mutual energy (output as an NMISC element item UM)

In a transient analysis, the kinetic energy (output as KENE) and the damping energy (output as DENE) are calculated as shown in Equation 14–342 and Equation 14–354, respectively.

In harmonic and modal piezoelectric analyses, the time-averaged element energies are calculated as:

(10–78)

(10–79)

(10–80)

where:

= complex conjugate of the elastic strain

= complex conjugate of the electric flux density

The time-averaged kinetic and damping energies are calculated as shown in Equation 14–344 and Equation 14–355, respectively.

The real parts of Equation 10–78 and Equation 10–79 represent the average stored elastic and dielectric energies, respectively. The imaginary parts represent the average elastic and electric losses. Therefore, the quality factor can be calculated from the total stored energy as:

(10–81)

where:

= number of piezoelectric elements

The total stored energy + is output as SENE. The factor can therefore be derived from the real and imaginary records of SENE summed over the piezoelectric elements.

In a damped modal analysis, the quality factor can also be calculated using the real and imaginary part of the complex eigenvalue:

(10–82)

The mutual energy can be used to calculate the electromechanical coupling coefficient as:

(10–83)

The piezoelectric Poynting vector is calculated as a sum of mechanical and electrical (quasi-static approximation) Poynting vectors:

(10–84)

where:

= stress matrix

{v} = velocity

φ = electric potential

{J} = total electric current density

In a harmonic analysis, the time-averaged Poynting vector is:

(10–85)

where:

The calculated Poynting vector is output as P. In a harmonic analysis, the real part of the Poynting vector represents the average power flow.

In a transient piezoelectric analysis, the heat generation rate produced by viscous losses is calculated as follows:

(10–86)

In a harmonic piezoelectric analysis, the generated heat is due to structural and electrical losses. The structural part of the time-averaged heat rate per unit volume is calculated from the complex values of mechanical stress {T_m} and strain {S} as:

(10–87)

The time-averaged electrical heat, or Joule heat, rate is calculated as:

(10–88)

The total heat rate {Q} = {Q^m} + {Q^e} is stored in both the real and imaginary sets of the JHEAT output item and can be used as a heat generation load in a subsequent thermal analysis (LDREAD,HGEN).

Perfectly matched layers are artificial anisotropic materials that absorb all incoming waves without any reflections, except for the graze wave that parallels the PML interface in the propagation direction. PMLs are currently constructed for the propagation of elastic waves in the elastic media in harmonic response analysis.

Considering two Cartesian coordinate systems, (1) a global system with respect to an orthogonal basis and (2) a local system with respect to another orthogonal basis , the displacements and voltage of the piezoelectric PML medium in the local coordinate system are governed by the following equations:

(10–89)

(10–90)

(10–91)

(10–92)

(10–93)

(10–94)

where:

= components of the infinitesimal strain tensor in the local Cartesian coordinate system

= components of the material stiffness tensor in the local Cartesian coordinate system

= components of the piezoelectric stress matrix in the local Cartesian coordinate system

= components of the dielectric matrix in the local Cartesian coordinate system

= components of the electric flux in the local Cartesian coordinate system

= components of the electric field in the local Cartesian coordinate system

= non-zero complex-valued coordinate stretching functions in the direction of the local Cartesian coordinate system

= density of the elastic medium

The tensors and vectors in the local coordinate system are transformed into the global coordinate system using the rotation-of-basis transformation matrix with the component , such as vector and matrix . On multiplying Equation 10–89 and Equation 10–90 with s = s₁s₂s₃, and taking the fact that s_j is a function of only, the governing equations of the displacements and voltage in the global coordinate system are written as:

(10–95)

(10–96)

(10–97)

(10–98)

(10–99)

(10–100)

The PML material is defined using coupled-field elements (PLANE222, PLANE223, SOLID225, SOLID226, or SOLID227) with KEYOPT(15) = 1. The PML element coordinate system (PSYS command) uniquely identifies each PML region. For more information about using PMLs, refer to Perfectly Matched Layers (PML) in the Acoustic Analysis Guide.