Piezoelectric materials and devices have been widely used in many technical applications. Nowadays, the coupling between electrical and mechanical beha- viors is used in different devices based both on the so-called “direct piezoelectric effect” or the “converse piezoelectric effect”     .
Some newer relevant applications include (among others) the high voltage generation from transient dynamic impact processes in vehicles  .
Analysis of operating electrical and mechanical parameters of such processes can be done by using various analytical and numerical methods. Although ana- lytical approaches are limited to rather simple geometries and other restrictions (homogeneous or piecewise-homogeneous bodies, linear governing equations, etc.), they often provide exact solutions.
Analytical methods have been successfully used for many transient dynamic one-dimensional piezoelectric problems  -  . Among analytical methods for transient dynamic piezoelectric problems, the Laplace transform methods play a very significant role. They solve boundary value problems in the frequency- domain, possibly for complex frequencies, using transformed boundary condi- tions for a piezoelectric body. After obtaining such solutions, the transformation back to the time-domain employs special methods for the inversion of Laplace transforms. However, using the Laplace transform methods is not instrumental even for one-dimensional problems if nonlinear boundary conditions are con- sidered. Time-domain numerical methods (e.g., finite element or finite diffe- rence methods) that can be used under such conditions usually lack precision associated with the use of analytical methods. Therefore, development of time- domain analytical or semi-analytical methods combining advantages of analyti- cal and numerical methods can be of interest for such problems.
In this paper, a time-domain Green’s function method is implemented for so- lution of one-dimensional transient dynamic piezoelectric problems for thick- ness polarized disks or length polarized rods. This method stems from a time- domain representation formulas approach for transient dynamic piezoe-lectric problems described in  . For one-dimensional problems with a variety of boun- dary conditions including nonlinear ones, this method produces exact solutions which are shown below. Such solutions can be used both for analyses of longitu- dinal mode, piezoelectric devices and as benchmark solutions for numerical me- thods of piezoelectricity.
2. Representation Formulas
Consider a transversely isotropic homogeneous piezoelectric material (piezoe- lectric element) with the -axis as the poling direction and the plane as the isotropic plane. Let this piezoelectric material occupy a disk (or a cylinder) bounded in -direction by planes and where is the thickness of the disk (or the length of the cylinder). Consider a uniaxial strain state or a stress stress state in direction when there is only one non- zero component of strain or stress the others being zero. We assume that the non-zero stress and strain components, and also the displacement and electric displacement in the -direction, and the electric potential , depend only on and which is usually the case for a longitudinal mode piezoelectric element  :
Under conditions (1), we can use the following one-dimensional constitutive equations (both for the uniaxial strain state and for the uniaxial stress state) that relate the mechanical and electrical fields in (1):
where coefficients are different for the uniaxial strain and uniaxial stress cases.
Then the corresponding equations of motions can be written as
where and denote the body force in -direction and electric charge.
To simplify further notations we will denote and derivatives with respect to by and the prime, respectively, and will skip subindex 3 for the elastic displacement, electric displacement and body force components presented in (3). Then system (3) becomes
The Green’s functions for vector can be obtained using concentrated impulses instead of or in (4) when is substituted by the infinite media.
Since can be expressed through due to the second equation in (4), then the first equation in (4) can be presented as the one-dimensional wave equation for displacement :
is the Young’s modulus measured at constant .
The wave speed corresponding to Equation (5) is denoted below by
The Green’s function for corresponding to load is the well-known Green’s function for the one-dimensional wave Equation (5):
where is the Heaviside step function (right-continuous), i.e. for and for .
The corresponding Green’s function for is calculated using the second equation in (4):
Based on (6) and (7), the representation formula for the displacement vector in 3-D case described in  reduces to the following expression for the dis- placement component :
where is denoted by and it is taken into account that the outward normals to the lower and upper boundaries of the layer have opposite directions.
In many practical applications, the electric volume charges are absent. There- fore, we consider henceforth only the case when . Then the terms related to in the above expression can be simplified since, based on Equation (4) in this case, is spatially uniform:
Due to the property (9) the representation formula (8) can be rewritten as
To obtain a representation formula for , let us consider an auxiliary function
that has the following connection to the electric displacement:
According to (3), . Then, using the corresponding Green’s function and Equation (11), we get a representation formula for involving only boundary value of function and a spatially uniform electric displacement:
Formulas (12) and (11) lead to the following expression for :
After is calculated, can be determined using this calculated value, and boundary values of .
The representation formula (10) allows us to get representation formulas for the velocity and stress . Differentiating (10) with respect to time provides the following representation formula for the velocity:
To get a representation formula for the stress we need to use the first contitu- tive equation from (2) (in the new notations introduced after equations (3)) and expression (13) which gives the following expression for the stress:
After differentiating (10) with respect to and substituting the result into (15) we get
A representation formula for is not needed under assumption that since the electric displacement is uniform in space in this case and deter- mined solely by the electric boundary conditions.
3. Boundary Equations
The velocity representation formula (14) generates two boundary equations when tends to the upper and lower boundaries of the piezoelectric element, that is, when tends to or 0:
where denotes the time taken by the elastic wave to travel the thickness of the piezoelectric layer:
Similarly, the stress representation formula (16) generates the following boun- dary equations:
It is easy to verify that Equations (17) and (19), though presented in different forms, are equivalent. The same is true for the pair of Equations (18) and (20). Therefore, we shall use the equations in these pairs interchangeably.
We will not work with boundary equations that can be obtained directly from the displacement representation formula (10), since it is computationally more effective to determine at first unknown boundary values of the velocity , and then calculate unknown boundary values of the displacement by integrating the boundary velocity over time (using also an initial condition for ).
We also need to consider boundary values of the expression (13) for the electric potential. It is important to emphasize that two equations obtained from (13) when tends to or to 0 are equivalent and, therefore, they are pre- sented below as one equation:
The boundary equations presented above will be used in the next section to create an exact time domain solution procedure in the case when nonlinear damper boundary conditions are sprecified.
4. Nonlinear Damper Boundary Conditions and Exact Solutions
Suppose that the lower end face of the piezoelectric element is fixed to a non- linear damper. Let be a damping force acting on the lower end face which is defined by the following nonlinear relationship  :
where is the damping constant, is the damping exponent, and is the signum function defined for all real numbers (including 0 where its value is also 0). If , the direction of is opposite to . The exponent has a value 1 for a linear damper, but may vary in practice in the interval  creating a set of possible boundary conditions at . We assume that the force is uniformly distributed over the lower end face. Then (22) transforms into the following nonlinear (in general) boundary condition at the lower end face:
where is the lower end face area.
Consider additional assumptions that will be used to get exact solutions for the damper boundary conditions based on the results of the previous section. We suppose that the values of are defined for . In addition, let us assume henceforth that
which means, based on (2) and (3), that and are also zero inside the piezoelectric body at negative times. The next additional assumption is that
inside the piezoelectric body at any time in the sense of generalized functions. This also includes the assumption that the initial conditions for the elastic displacement are zero, as discussed in  . These assumptions will simplify using boundary Equations (17)-(20) for particular problems considered below.
Regarding the design of the piezoelectric element, we assume that it is a cylinder (or a rod) with two coated electrodes at and . The elec- trodes are considered to be of negligible thickness (from the mechanical point of view) and their deformation is neglected. The output voltage, which is defined as the electric potential difference between the lower and upper electrodes , is of primary interest below.
The electric boundary condition at corresponds to the grounded electrode:
At the upper end face, the following mechanical boundary condition is used:
where is an applied normal stress load which is assumed to be known and negative.
The electric boundary condition at the upper end face is formulated as follows:
So, based on (9), .
Using the above assumptions the representation formulas (14) and (16) for the velocity and stress take the following simplified forms:
where all the time dependent functions are equal to zero for negative times.
In the representation formulas (29) and (30), there are three unknown boundary functions and first two of which are related by Equation (23). Two additional equations needed for determination of these three functions will be derived below based on (17) and (18).
After the the velocity is determined for any particular over time, the corresponding displacement can be obtained (due to the zero initial conditions) as
Boundary values of the displacement provide (according to (21) and (26)) the electric potential value at :
4.1. An Exact Recursive Procedure
The solution of the above problem will be obtained by using an exact recursive procedure based on the following equations obtained from (17) and (18) under the boundary conditions (23) (26) (27) (28):
There are two unknowns and at each time moment in these equations. The right-hand sides of the equations are known at each time point since they contain either or time-dalayed function values at that had to be determined at a previous step of the recursive process.
In order to simplify deriving next results, we need to introduce some addi- tional notations:
Then, Equation (34) reads as
Let the left-hand side of Equation (36) be denoted by . Since , is a continuous strictly monotonically increasing function on ranging from to . Therefore, for any real , there exists one and only one solution of Equation (36) in .
Denote by the operator that tranforms into this solution of equation (36). Thus, is the left inverse operator of the nonlinear operator acting on in the left-hand side of Equation (36). If , 1, or , the corres- ponding expressions of are very simple for computations:
The calculation of for other values of can effectively be imple- mented using a symbolic computation software like Maple  .
With help of the inverse operator Equation (34) can be rewritten in the following explicit form for calculating :
Equation (38) combined with (33) creates the recursive procedure that can be used directly for calculations or can lead to building explicit exact solution for vector step by step over consecutive time intervals . In doing so, it is helpful to substitute in (33) by its expression obtained from (38) which provides the following recursive equation for :
or, using the identity operator (that leaves unchanged the element on which it operates),
4.2. Explicit Exact Solutions
Now we derive some explicit exact solutions for and corres- ponding to three practically important ranges of the duration of the stress load at . Our goal is to present the boundary velocities directly through the transient stress load at that generates the dynamic process in the piezoelectric body.
4.2.1. Case 1:
So, if . Using the recursive Equation (39) under this condition for consecutive intervals , we finally obtain the following explicit expression for :
Substituting (40) into (38) we get the corresponding explicit expression for :
4.2.2. Case 2:
In this case, if . Acting similarly to section 4 we derive the following explicit exact solutions:
4.2.3. Case 3:
So, if , and the corresponding exact solutions have the following closed form:
Similar explicit formulas for are excessively cumbersome. In this case, it is easier to directly use the recursive procedure based on (33) and (38) which has the same simple form regardless of the transient load duration and also provides exact results.
5. Examples and Discussions
Consider some examples of using the results of the previous section for mathematical modeling of piezoelectric cylindrical devices installed in a car as proposed in  . These devices transform the mechanical energy of the moving pistons or crank-shafts into electrical energy, which will be stored in the capacitor or the battery charger. We consider the uniaxial stress state for a cylinder and assume that the material of the cylinder is PZT-5A  . In this case, parameters and in Equations (3) have the following values:
Then the elastic wave speed in the piezoelectric material is equal to . Next, we take into account that the total force instantaneously applied to the top of a piston in an internal combustion engine is around 6300 pounds, which corresponds to approximately 28,640 N  . Suppose that this force is applied downward to a piezoelectric cylinder with a length of and a diameter . So, the area of the upper end face of the cylinder . Assuming that is uniformly distributed over the upper end face, we get the amplitude of the pressure impulse acting on the top of the cylinder: . Let us assume that the applied normal stress load takes the form of the following rectangular compressive load (pressure) impulse:
where is the duration of the pressure impulse.
We assume first that and (see, e.g.,  ) in the damper boundary conditions (23). Consider the following three values of . Since the transit time of elastic waves between the upper and lower end faces , then these three durations correspond to the three cases of explicit exact solutions considered in 4. The operator is calculated according to (37). The calculated results for the output voltage (in kV) based on these exact solutions are presented for in Figures 1-3 below.
Comparing the graphs we can see that the maximum or peak value does not depend on the pressure impulse duration . However, the number of peaks in each figure depends on the . The time distance between two neighboring peaks is approximately equal to . After the pressure load is removed, there is an attenuation of the output voltage vibrations.
Now let us consider another set of the damper parameters: and . The parameters of the material and the impulse durations are the same as above. The operator is also calculated according to (37).
Figure 1. Output voltage for and .
Figure 2. Output voltage for and .
Figure 3. Output voltage for and .
The calculated results for the output voltage (in kV) are presented for in Figures 4-6.
Comparison of these graphs shows that the maximum value of the output voltage does not depend on the pressure impulse duration which similar to the case when . The difference is that now there is only one peak but its width depends on the . After the pressure load is removed, the attenuation of the output voltage vibrations is very pronounced: the amplitude of vibrations after the load removal is almost negligible.
Figure 4. Output voltage for and .
Figure 5. Output voltage for and .
One-dimensional transient dynamic piezoelectric problems for thickness polarized layers and disks, or length polarized rods, are considered here in the framework of a time-domain Green’s function method. As the result, a novel exact analytical recursive procedure is derived which is applicable for a wide variety of boundary conditions including the nonlinear damper case. Some new practically important explicit exact solutions are presented. The effectiveness of the proposed exact approach is demonstrated by examples of the time behavior of the output electric potential difference between two electrodes coated at the
Figure 6. Output voltage for and .
end faces of a piezoelectric cylinder fixed to a nonlinear damper at one end, and subjected to impulsive loading at the other.
The authors would like to thank the anonymous reviewers for their comments.
 Cook, E.G. (1956) Transient and Steady-State Response of Ultrasonic Piezoelectric Transducers. 1956 IRE National Convention, New York, 19-22 March 1956, 61-69.
 Zhang, H.L., Li, M.X. and Ying, C.F. (1983) Complete Solutions of the Transient Behavior of a Transmitting Thickness-Mode Piezoelectric Transducer and Their Physical Interpretations. The Journal of the Acoustical Society of America, 74, 1105.
 Ing, Y.S., Liao, H.F. and Huang, K.S. (2013) The Extended Durbin Method and Its Application for Piezoelectric Wave Propagation Problems. International Journal of Solids and Structures, 50, 4000-4009.
 Gazonas, G.A., Wildman, R.A., Hopkins, D.A. and Scheidler, M.J. (2016) Longitudinal Impact of Piezoelectric Media. Archive of Applied Mechanics, 86, 497-515.
 Ho, C., Lang, Z.Q., Sapiński, B. and Billings, S.A. (2013) Vibration Isolation Using Nonlinear Damping Implemented by a Feedback-Controlled MR Damper. Smart Materials and Structures, 22. https://doi.org/10.1088/0964-1726/22/10/105010