Exact Time Domain Solutions of 1-D Transient Dynamic Piezoelectric Problems with Nonlinear Damper Boundary Conditions

Naum M. Khutoryansky, Vladimir Genis

Department of Engineering Technology, Drexel University, Philadelphia, PA, USA.

DOI: 10.4236/jamp.2017.54077                

Abstract

Novel exact solutions of one-dimensional transient dynamic piezoelectric problems for thickness polarized layers and disks, or length polarized rods, are obtained. The solutions are derived using a time-domain Green’s function method that leads to an exact analytical recursive procedure which is applicable for a wide variety of boundary conditions including nonlinear cases. A nonlinear damper boundary condition is considered in more detail. The corresponding nonlinear relationship between stresses and velocities at a current time moment is used in the recursive procedure. In addition to the exact recursive procedure that is effective for calculations, some new practically important explicit exact solutions are presented. Several examples of the time behavior of the output electric potential difference are given to illustrate the effectiveness of the proposed exact approach.

Keywords

Piezoelectric Layer, Transient Dynamic Problems, Time Domain Solutions, Green’s Function Method, Nonlinear Boundary Conditions, Nonlinear Damper, Output Voltage

1. Introduction

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” [1] [2] [3] [4] .

Some newer relevant applications include (among others) the high voltage generation from transient dynamic impact processes in vehicles [5] .

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 [6] - [14] . 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 [15] . 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 $x_3$ -axis as the poling direction and the $x_1-x_2$ plane as the isotropic plane. Let this piezoelectric material occupy a disk (or a cylinder) $\Omega$ bounded in $x_3$ -direction by planes $x_3=0$ and $x_3=h$ where $h>0$ is the thickness of the disk (or the length of the cylinder). Consider a uniaxial strain state or a stress stress state in $x_3$ direction when there is only one non- zero component of strain ${\gamma }_{33}$ or stress ${\sigma }_{33}$ the others being zero. We assume that the non-zero stress and strain components, and also the displacement $u_3$ and electric displacement $D_3$ in the $x_3$ -direction, and the electric potential $\varphi$ , depend only on $x_3$ and $t$ which is usually the case for a longitudinal mode piezoelectric element [16] :

$\begin{array}{l}{\sigma }_{33}={\sigma }_{33}\left(x_3,t\right),D_3=D_3\left(x_3,t\right),\\ u_3=u_3\left(x_3,t\right),\varphi =\varphi \left(x_3,t\right)\text{in}\Omega \end{array}$ (1)

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):

${\sigma }_{33}=C{u}_{3,3}+e{\varphi }_{,3};D_3=e{u}_{3,3}-ϵ{\varphi }_{,3}$ (2)

where coefficients are different for the uniaxial strain and uniaxial stress cases.

Then the corresponding equations of motions can be written as

$\begin{array}{l}\rho {\ddot{u}}_3-C{u}_{3,33}-e{\varphi }_{,33}=b_3;\\ -D_{3,3}=ϵ{\varphi }_{,33}-e{u}_{3,33}=-q\text{in}\Omega \end{array}$ (3)

where $b_3=b_3\left(x_3,t\right)$ and $q=q\left(x_3,t\right)$ denote the body force in $x_3$ -direction and electric charge.

To simplify further notations we will denote $x_3$ and derivatives with respect to $x_3$ by $x$ 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

$\begin{array}{l}\rho \ddot{u}-C{u}^{″}-e{\varphi }^{″}=b;\\ -D^{\prime }=ϵ{\varphi }^{″}-e{u}^{″}=-q\text{in}\Omega \mathrm{.}\end{array}$ (4)

The Green’s functions for vector $\left\{u\mathrm{,}\varphi \right\}$ can be obtained using concentrated impulses instead of $b$ or $q$ in (4) when $\Omega$ is substituted by the infinite media.

Since ${\varphi }^{″}$ can be expressed through $u^{″}$ due to the second equation in (4), then the first equation in (4) can be presented as the one-dimensional wave equation for displacement $u$ :

$\rho \ddot{u}-C^D{u}^{″}=b-\frac{e}{ϵ}q$ (5)

where

$C^D=C+\frac{e^2}{ϵ}$

is the Young’s modulus measured at constant $D$ .

The wave speed corresponding to Equation (5) is denoted below by

$c=\sqrt{\frac{C^D}{\rho }}$

The Green’s function for $u$ corresponding to load $\left\{b=\delta \left(x\right)\delta \left(t\right),q=0\right\}$ is the well-known Green’s function for the one-dimensional wave Equation (5):

$U\left(x,t\right)=\frac{1}{2\rho c}H\left(t-\left|x\right|/c\right)$ (6)

where $H\left(t\right)$ is the Heaviside step function (right-continuous), i.e. $H\left(t\right)=0$ for $t<0$ and $H\left(t\right)=1$ for $t\ge 0$ .

The corresponding Green’s function for $\varphi$ is calculated using the second equation in (4):

$U_{\varphi }\left(x,t\right)=\frac{e}{ϵ}U\left(x,t\right)$ (7)

Based on (6) and (7), the representation formula for the displacement vector in 3-D case described in [15] reduces to the following expression for the dis- placement component $u_3\left(x,t\right)=u\left(x,t\right)$ :

$\begin{array}{c}u\left(x,t\right)=\frac{1}{2}\left[u\left(0,t-x/c\right)+u\left(h,t-\left(h-x\right)/c\right)\right]\\ -\frac{1}{2\rho c}{\displaystyle {\int }_{-\infty }^{t-x/c}}\left[\sigma \left(0,\tau \right)+\frac{e}{ϵ}D\left(0,\tau \right)\right]\text{d}\tau \\ +\frac{1}{2\rho c}{\displaystyle {\int }_{-\infty }^{t-\left(h-x\right)/c}}\left[\sigma \left(h,\tau \right)+\frac{e}{ϵ}D\left(h,\tau \right)\right]\text{d}\tau \\ +\frac{1}{2\rho c}{\displaystyle {\int }_0^h}{\displaystyle {\int }_{-\infty }^{t-\left|\xi -x\right|/c}}\left[b\left(\xi ,\tau \right)-\frac{e}{ϵ}q\left(\xi ,\tau \right)\right]\text{d}\tau \text{d}\xi \text{in}\Omega \mathrm{,}\end{array}$ (8)

where ${\sigma }_{33}$ is denoted by $\sigma$ and it is taken into account that the outward normals to the lower and upper boundaries of the layer $0\le x\le h$ have opposite directions.

In many practical applications, the electric volume charges are absent. There- fore, we consider henceforth only the case when $q=0$ . Then the terms related to $D$ in the above expression can be simplified since, based on Equation (4) in this case, $D\left(x\mathrm{,}t\right)$ is spatially uniform:

$D\left(x,t\right)=D\left(t\right).$ (9)

Due to the property (9) the representation formula (8) can be rewritten as

$\begin{array}{c}u\left(x,t\right)=\frac{1}{2}\left[u\left(0,t-x/c\right)+u\left(h,t-\left(h-x\right)/c\right)\right]\\ -\frac{1}{2\rho c}{\displaystyle {\int }_{-\infty }^{t-x/c}}\left[\sigma \left(0,\tau \right)+\frac{e}{ϵ}D\left(\tau \right)\right]\text{d}\tau \\ +\frac{1}{2\rho c}{\displaystyle {\int }_{-\infty }^{t-\left(h-x\right)/c}}\left[\sigma \left(h,\tau \right)+\frac{e}{ϵ}D\left(\tau \right)\right]\text{d}\tau \\ +\frac{1}{2\rho c}{\displaystyle {\int }_0^h}{\displaystyle {\int }_{-\infty }^{t-\left|\xi -x\right|/c}}b\left(\xi ,\tau \right)\text{d}\tau \text{d}\xi \text{in}\Omega \mathrm{.}\end{array}$ (10)

To obtain a representation formula for $\varphi \left(x\mathrm{,}t\right)$ , let us consider an auxiliary function

$\psi \left(x,t\right)=\varphi \left(x,t\right)-\frac{e}{ϵ}u\left(x,t\right)$ (11)

that has the following connection to the electric displacement:

$D=-ϵ{\psi }^{\prime }\mathrm{.}$

According to (3), ${\psi }^{″}=-D^{\prime }/ϵ=0$ . Then, using the corresponding Green’s function $\left|x\right|/2$ and Equation (11), we get a representation formula for $\psi \left(x\mathrm{,}t\right)$ involving only boundary value of function $\psi \left(x\mathrm{,}t\right)$ and a spatially uniform electric displacement:

$\psi \left(x,t\right)=\frac{1}{2}\left[\psi \left(0,t\right)+\psi \left(h,t\right)\right]+\frac{h-2x}{2ϵ}D\left(t\right)\text{in}\Omega \mathrm{.}$ (12)

Formulas (12) and (11) lead to the following expression for $\varphi \left(x\mathrm{,}t\right)$ :

$\begin{array}{c}\varphi \left(x,t\right)=\frac{1}{2}\left[\varphi \left(0,t\right)+\varphi \left(h,t\right)\right]-\frac{e}{2ϵ}\left[u\left(0,t\right)+u\left(h,t\right)\right]\\ +\frac{e}{ϵ}u\left(x,t\right)+\frac{h-2x}{2ϵ}D\left(t\right)\text{in}\Omega \mathrm{.}\end{array}$ (13)

After $u\left(x\mathrm{,}t\right)$ is calculated, $\varphi \left(x\mathrm{,}t\right)$ can be determined using this calculated value, $D\left(t\right)$ and boundary values of $\varphi \left(x\mathrm{,}t\right)$ .

The representation formula (10) allows us to get representation formulas for the velocity $v\left(x,t\right)=\stackrel{\dot{}}{u}\left(x,t\right)$ and stress $\sigma \left(x\mathrm{,}t\right)$ . Differentiating (10) with respect to time provides the following representation formula for the velocity:

$\begin{array}{c}v\left(x,t\right)=\frac{1}{2}\left[v\left(0,t-x/c\right)+v\left(h,t-\left(h-x\right)/c\right)\right]\\ -\frac{1}{2\rho c}\left[\sigma \left(0,t-x/c\right)+\frac{e}{ϵ}D\left(t-x/c\right)\right]\\ +\frac{1}{2\rho c}\left[\sigma \left(h,t-\left(h-x\right)/c\right)+\frac{e}{ϵ}D\left(t-\left(h-x\right)/c\right)\right]\\ +\frac{1}{2\rho c}{\displaystyle {\int }_0^h}b\left(\xi ,t-\left|\xi -x\right|/c\right)\text{d}\xi \text{in}\Omega .\end{array}$ (14)

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:

$\begin{array}{c}\sigma \left(x,t\right)=C{u}^{\prime }\left(x,t\right)+e\left[\frac{e}{ϵ}{u}^{\prime }\left(x,t\right)-\frac{1}{ϵ}D\left(t\right)\right]\\ =C^D{u}^{\prime }\left(x,t\right)-\frac{e}{ϵ}D\left(t\right).\end{array}$ (15)

After differentiating (10) with respect to $x$ and substituting the result into (15) we get

$\begin{array}{c}\sigma \left(x,t\right)=\frac{\rho c}{2}\left[v\left(0,t-x/c\right)+v\left(h,t-\left(h-x\right)/c\right)\right]\\ +\frac{1}{2}\left[\sigma \left(0,t-x/c\right)+\sigma \left(h,t-\left(h-x\right)/c\right)\right]\\ +\frac{e}{2ϵ}\left[D\left(t-x/c\right)+D\left(t-\left(h-x\right)/c\right)-2D\left(t\right)\right]\\ +\frac{1}{2}{\displaystyle {\int }_0^h}b\left(\xi ,t-\left|\xi -x\right|/c\right)\mathrm{sgn}\left(x-\xi \right)\text{d}\xi \text{in}\Omega .\end{array}$ (16)

A representation formula for $D\left(x\mathrm{,}t\right)$ is not needed under assumption that $q=0$ 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 $x$ tends to the upper and lower boundaries of the piezoelectric element, that is, when $x$ tends to $h$ or 0:

$\begin{array}{c}v\left(h,t\right)=v\left(0,t-\theta \right)+\frac{1}{\rho c}\left[\sigma \left(h,t\right)-\sigma \left(0,t-\theta \right)\right]\\ +\frac{e}{\rho cϵ}\left[D\left(t\right)-D\left(t-\theta \right)\right]+\frac{1}{\rho c}{\displaystyle {\int }_0^h}b\left(\xi ,t-\theta +\xi /c\right)\text{d}\xi \mathrm{,}\end{array}$ (17)

$\begin{array}{c}v\left(\mathrm{0,}t\right)=v\left(h\mathrm{,}t-\theta \right)+\frac{1}{\rho c}\left[\sigma \left(h\mathrm{,}t-\theta \right)-\sigma \left(\mathrm{0,}t\right)\right]\\ +\frac{e}{\rho cϵ}\left[D\left(t-\theta \right)-D\left(t\right)\right]+\frac{1}{\rho c}{\displaystyle {\int }_0^h}b\left(\xi \mathrm{,}t-\xi /c\right)\text{d}\xi \end{array}$ (18)

where $\theta$ denotes the time taken by the elastic wave to travel the thickness of the piezoelectric layer:

$\theta =\frac{h}{c}.$

Similarly, the stress representation formula (16) generates the following boun- dary equations:

$\begin{array}{c}\sigma \left(h\mathrm{,}t\right)=\sigma \left(\mathrm{0,}t-\theta \right)+\rho c\left[v\left(h\mathrm{,}t\right)-v\left(\mathrm{0,}t-\theta \right)\right]\\ +\frac{e}{ϵ}\left[D\left(t-\theta \right)-D\left(t\right)\right]-{\displaystyle {\int }_0^h}b\left(\xi \mathrm{,}t-\theta +\xi /c\right)\text{d}\xi \mathrm{,}\end{array}$ (19)

$\begin{array}{c}\sigma \left(0,t\right)=\sigma \left(h,t-\theta \right)+\rho c\left[v\left(h,t-\theta \right)-v\left(0,t\right)\right]\\ +\frac{e}{ϵ}\left[D\left(t-\theta \right)-D\left(t\right)\right]+{\displaystyle {\int }_0^h}b\left(\xi \mathrm{,}t-\xi /c\right)\text{d}\xi \mathrm{.}\end{array}$ (20)

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 $v\left(x\mathrm{,}t\right)$ , and then calculate unknown boundary values of the displacement $u\left(x\mathrm{,}t\right)$ by integrating the boundary velocity over time (using also an initial condition for $u\left(x\mathrm{,}t\right)$ ).

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 $x$ tends to $h$ or to 0 are equivalent and, therefore, they are pre- sented below as one equation:

$\varphi \left(h\mathrm{,}t\right)-\varphi \left(\mathrm{0,}t\right)\mathrm{<mo>\; =\; </mo>}\frac{e}{ϵ}\left[u\left(h\mathrm{,}t\right)-u\left(\mathrm{0,}t\right)\right]-\frac{h}{ϵ}D\left(t\right)$ (21)

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 $F$ be a damping force acting on the lower end face which is defined by the following nonlinear relationship [17] :

$F=-k_{\alpha }{\left|v\left(0,t\right)\right|}^{\alpha }\mathrm{sgn}\left(v\left(0,t\right)\right)\mathrm{.}$ (22)

where $k_{\alpha }>0$ is the damping constant, $\alpha >0$ is the damping exponent, and $\mathrm{sgn}(.)$ is the signum function defined for all real numbers (including 0 where its value is also 0). If $v\left(\mathrm{0,}t\right)\ne 0$ , the direction of $F$ is opposite to $v\left(\mathrm{0,}t\right)$ . The exponent $\alpha$ has a value 1 for a linear damper, but may vary in practice in the interval $\left(\mathrm{0,2}\right]$ [17] creating a set of possible boundary conditions at $x=0$ . We assume that the force $F$ is uniformly distributed over the lower end face. Then (22) transforms into the following nonlinear (in general) boundary condition at the lower end face:

$\sigma \left(0,t\right)=\frac{k_{\alpha }}{A}{\left|v\left(0,t\right)\right|}^{\alpha }\mathrm{sgn}\left(v\left(0,t\right)\right)$ (23)

where $A$ 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 $u,\sigma ,b,\varphi ,D$ are defined for $-\infty <t<\infty$ . In addition, let us assume henceforth that

$u\left(x,t\right)=0,\varphi \left(x,t\right)=0\text{if}0<x<h,t<0$ (24)

which means, based on (2) and (3), that $\sigma \mathrm{,}D$ and $b$ are also zero inside the piezoelectric body at negative times. The next additional assumption is that

$b\left(x,t\right)=0$ (25)

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 $u\left(x\mathrm{,}t\right)$ are zero, as discussed in [15] . 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 $x=0$ and $x=h$ . 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 $\Delta \varphi =\varphi \left(0,t\right)-\varphi \left(h,t\right)$ , is of primary interest below.

The electric boundary condition at $x=0$ corresponds to the grounded electrode:

$\varphi \left(0,t\right)=0\text{if}t\ge 0$ (26)

At the upper end face, the following mechanical boundary condition is used:

$\sigma \left(h,t\right)=p\left(t\right)\mathrm{if}t\ge 0$ (27)

where $p\left(t\right)$ is an applied normal stress load which is assumed to be known and negative.

The electric boundary condition at the upper end face $x=h$ is formulated as follows:

$D\left(h,t\right)=0\text{if}t\ge 0.$ (28)

So, based on (9), $D\left(t\right)=0$ .

Using the above assumptions the representation formulas (14) and (16) for the velocity and stress take the following simplified forms:

$\begin{array}{c}v\left(x,t\right)=\frac{1}{2}\left[v\left(0,t-x/c\right)+v\left(h,t-\left(h-x\right)/c\right)\right]\\ +\frac{1}{2\rho c}\left[p\left(t-\left(h-x\right)/c\right)-\sigma \left(0,t-x/c\right)\right]\text{in}\Omega \mathrm{,}\end{array}$ (29)

$\begin{array}{c}\sigma \left(x,t\right)=\frac{\rho c}{2}\left[v\left(0,t-x/c\right)+v\left(h,t-\left(h-x\right)/c\right)\right]\\ +\frac{1}{2}\left[\sigma \left(0,t-x/c\right)+p\left(t-\left(h-x\right)/c\right)\right]\text{in}\Omega \end{array}$ (30)

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 $v\left(\mathrm{0,}t\right)\mathrm{,}\sigma \left(\mathrm{0,}t\right)$ and $v\left(h\mathrm{,}t\right)$ 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 $v\left(x\mathrm{,}t\right)$ is determined for any particular $x$ over time, the corresponding displacement $u\left(x\mathrm{,}t\right)$ can be obtained (due to the zero initial conditions) as

$u\left(x\mathrm{,}t\right)={\displaystyle {\int }_0^t}v\left(x\mathrm{,}\tau \right)\text{d}\tau \mathrm{.}$ (31)

Boundary values of the displacement provide (according to (21) and (26)) the electric potential value at $x=h$ :

$\varphi \left(h\mathrm{,}t\right)=\frac{e}{ϵ}\left[u\left(h\mathrm{,}t\right)-u\left(\mathrm{0,}t\right)\right]$ (32)

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):

$v\left(h,t\right)=2v\left(0,t-\theta \right)-v\left(h,t-2\theta \right)+\frac{1}{\rho c}\left[p\left(t\right)-p\left(t-2\theta \right)\right],$ (33)

$v\left(0,t\right)+\frac{k_{\alpha }}{\rho cA}{\left|v\left(0,t\right)\right|}^{\alpha }\mathrm{sgn}\left(v\left(0,t\right)\right)=v\left(h,t-\theta \right)+\frac{1}{\rho c}p\left(t-\theta \right).$ (34)

There are two unknowns $v\left(h\mathrm{,}t\right)$ and $v\left(\mathrm{0,}t\right)$ at each time moment $t$ in these equations. The right-hand sides of the equations are known at each time point since they contain either $p\left(t\right)$ or time-dalayed function values at $t-\theta$ 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:

$\gamma =\frac{k_{\alpha }}{\rho cA},\xi =v\left(0,t\right),r=v\left(h,t-\theta \right)+\frac{1}{\rho c}p\left(t-\theta \right).$ (35)

Then, Equation (34) reads as

$\xi +\gamma {\left|\xi \right|}^{\alpha }\mathrm{sgn}\left(\xi \right)=r.$ (36)

Let the left-hand side of Equation (36) be denoted by $f\left(\xi \right)$ . Since $\alpha >0$ , $f\left(\xi \right)$ is a continuous strictly monotonically increasing function on $\left(-\infty \mathrm{,}\infty \right)$ ranging from $-\infty$ to $\infty$ . Therefore, for any real $r$ , there exists one and only one solution of Equation (36) in $\left(-\infty \mathrm{,}\infty \right)$ .

Denote by $Q_{\alpha }$ the operator that tranforms $r$ into this solution of equation (36). Thus, $Q_{\alpha }$ is the left inverse operator of the nonlinear operator acting on $\xi$ in the left-hand side of Equation (36). If $\alpha =2$ , 1, $\text{1}/\text{2}$ or $\text{1}/\text{3}$ , the corres- ponding expressions of $Q_{\alpha }r$ are very simple for computations:

$\{\begin{array}{l}{Q}_{2}r=\frac{1}{2\gamma }\left(-1+\sqrt{1+4\left|r\right|\gamma }\right)\mathrm{sgn}\left(r\right),\\ Q_{1}r=\frac{r}{1+\gamma },Q_{1/2}r=\frac{1}{4}{\left(-\gamma +\sqrt{{\gamma }^2+4\left|r\right|}\right)}^2\mathrm{sgn}\left(r\right),\\ Q_{1/3}r=-\gamma \left[\frac{1}{6}{\left(108r+12\sqrt{12{\gamma }^3+81r^2}\right)}^{1/3}-\frac{2\gamma }{{\left(108r+12\sqrt{12{\gamma }^3+81r^2}\right)}^{1/3}}\right]+r.\end{array}$ (37)

The calculation of $Q_{\alpha }r$ for other values of $\alpha$ can effectively be imple- mented using a symbolic computation software like Maple [18] .

With help of the inverse operator $Q_{\alpha }$ Equation (34) can be rewritten in the following explicit form for calculating $v\left(\mathrm{0,}t\right)$ :

$v\left(0,t\right)=Q_{\alpha }\left[v\left(h,t-\theta \right)+\frac{1}{\rho c}p\left(t-\theta \right)\right].$ (38)

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 $\left\{v\left(\mathrm{0,}t\right)\mathrm{,}v\left(h\mathrm{,}t\right)\right\}$ step by step over consecutive time intervals $j\theta \le t<\left(j+1\right)\theta \left(j=0,1,2,\cdots \right)$ . In doing so, it is helpful to substitute $v\left(\mathrm{0,}t-\theta \right)$ in (33) by its expression obtained from (38) which provides the following recursive equation for $v\left(h\mathrm{,}t\right)$ :

$\begin{array}{c}v\left(h,t\right)=2Q_{\alpha }\left[v\left(h,t-2\theta \right)+\frac{1}{\rho c}p\left(t-2\theta \right)\right]\\ -\left[v\left(h,t-2\theta \right)+\frac{1}{\rho c}p\left(t-2\theta \right)\right]+\frac{1}{\rho c}p\left(t\right),\end{array}$

or, using the identity operator $I$ (that leaves unchanged the element on which it operates),

$v\left(h,t\right)=\left(2Q_{\alpha }-I\right)\left[v\left(h,t-2\theta \right)+\frac{1}{\rho c}p\left(t-2\theta \right)\right]+\frac{1}{\rho c}p\left(t\right).$ (39)

4.2. Explicit Exact Solutions

Now we derive some explicit exact solutions for $v\left(h\mathrm{,}t\right)$ and $v\left(\mathrm{0,}t\right)$ corres- ponding to three practically important ranges of the duration $t_1$ of the stress load at $x=h$ . Our goal is to present the boundary velocities directly through the transient stress load at $x=h$ that generates the dynamic process in the piezoelectric body.

4.2.1. Case 1: $t_1<2\theta$

So, $p\left(t\right)=0$ if $t\notin \left[0,2\theta \right)$ . Using the recursive Equation (39) under this condition for consecutive intervals $\left[2k\theta ,2\left(k+1\right)\theta \right),k=0,1,2,\cdots$ , we finally obtain the following explicit expression for $v\left(h\mathrm{,}t\right)$ :

$v\left(h,t\right)=\{\begin{array}{l}\frac{1}{\rho c}p\left(t\right)\text{if}0\le t<2\theta ,\\ {\left(2Q_{\alpha }-I\right)}^k\left[\frac{2}{\rho c}p\left(t-2k\theta \right)\right]\\ \text{if}2k\theta \le t<2\left(k+1\right)\theta ,k=1,2,\cdots .\end{array}$ (40)

Substituting (40) into (38) we get the corresponding explicit expression for $v\left(\mathrm{0,}t\right)$ :

$v\left(0,t\right)=\{\begin{array}{l}0\text{if}0\le t<\theta ,\\ Q_{\alpha }{\left(2Q_{\alpha }-I\right)}^{k-1}\left[\frac{2}{\rho c}p\left(t-2k\theta \right)\right]\\ \text{if}\left(2k-1\right)\theta \le t<\left(2k+1\right)\theta ,k=1,2,\cdots .\end{array}$ (41)

4.2.2. Case 2: $t_1<4\theta$

In this case, $p\left(t\right)=0$ if $t\notin \left[0,4\theta \right)$ . Acting similarly to section 4 we derive the following explicit exact solutions:

$v\left(h,t\right)=\{\begin{array}{l}\frac{1}{\rho c}p\left(t\right)\text{if}0\le t<2\theta ,\\ {\left(2Q_{\alpha }-I\right)}^{k-1}[\frac{1}{\rho c}p\left(t-2\left(k-1\right)\theta \right)\\ +\left(2Q_{\alpha }-I\right)\left(\frac{2}{\rho c}p\left(t-2k\theta \right)\right)]\text{if}2k\theta \le t<2\left(k+1\right)\theta ,k=1,2,\cdots ;\end{array}$ (42)

$v\left(0,t\right)=\{\begin{array}{l}0\text{if}0\le t<\theta ,\\ Q_{\alpha }\left[\frac{2}{\rho c}p\left(t-\theta \right)\right]\text{if}\theta \le t<3\theta ,\\ Q_{\alpha }{\left(2Q_{\alpha }-I\right)}^{k-1}[\frac{1}{\rho c}p\left(t-2k\theta +\theta \right)\\ +\left(2Q_{\alpha }-I\right)\left(\frac{2}{\rho c}p\left(t-2k\theta -\theta \right)\right)]\\ \text{if}\left(2k-1\right)\theta \le t<\left(2k+1\right)\theta ,k=2,3,\cdots .\end{array}$ (43)

4.2.3. Case 3: $t_1<6\theta$

So, $p\left(t\right)=0$ if $t\notin \left[0,6\theta \right)$ , and the corresponding exact solutions have the following closed form:

$v\left(h,t\right)=\{\begin{array}{l}\frac{1}{\rho c}p\left(t\right)\text{if}0\le t<2\theta ,\\ \frac{1}{\rho c}p\left(t\right)+\left(2Q_{\alpha }-I\right)\left[\frac{2}{\rho c}p\left(t-2\theta \right)\right]\text{if}2\theta \le t<4\theta ,\\ {\left(2Q_{\alpha }-I\right)}^{k-1}[\frac{1}{\rho c}p\left(t-2\left(k-1\right)\theta \right)+\left(2Q_{\alpha }-I\right)[\frac{2}{\rho c}p\left(t-2k\theta \right)\\ +\left(2Q_{\alpha }-I\right)\left(\frac{2}{\rho c}p\left(t-2\left(k+1\right)\theta \right)\right)]]\\ \text{if}2\left(k+1\right)\theta \le t<2\left(k+2\right)\theta ,k=1,2,\cdots ;\end{array}$ (44)

$v\left(0,t\right)=\{\begin{array}{l}0\text{if}0\le t<\theta ,\\ Q_{\alpha }\left[\frac{2}{\rho c}p\left(t-\theta \right)\right]\text{if}\theta \le t<3\theta ,\\ Q_{\alpha }\left[\frac{1}{\rho c}p\left(t-\theta \right)+\left(2Q_{\alpha }-I\right)\left(\frac{2}{\rho c}p\left(t-3\theta \right)\right)\right]\text{if}3\theta \le t<5\theta ,\\ Q_{\alpha }{\left(2Q_{\alpha }-I\right)}^{k-1}[\frac{1}{\rho c}p\left(t-2k\theta +\theta \right)\\ +\left(2Q_{\alpha }-I\right)[\frac{2}{\rho c}p\left(t-2k\theta -\theta \right)\\ +\left(2Q_{\alpha }-I\right)\left(\frac{2}{\rho c}p\left(t-2k\theta -3\theta \right)\right)]]\\ \text{if}\left(2k+3\right)\theta \le t<\left(2k+5\right)\theta ,k=1,2,\cdots .\end{array}$ (45)

Similar explicit formulas for $t_1\ge 6\theta$ 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 [5] . 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 [4] . In this case, parameters $C\mathrm{,}e\mathrm{,}ϵ$ and $\rho$ in Equations (3) have the following values:

$\begin{array}{l}C={\mathrm{5.32\; <mo>\; *\; </mo>10}}^{10}\text{N}/{\text{m}}^{\text{2}}\mathrm{<mo>\; ;\; </mo>}e=19.89\text{N}/\left(\text{V}\cdot \text{m}\right)\mathrm{<mo>\; ;\; </mo>}\\ ϵ={\mathrm{76.12\; <mo>\; *\; </mo>10}}^{-10}\text{farad}/\text{m}\mathrm{<mo>\; ;\; </mo>}\rho =7750.0\text{kg}/{\text{m}}^{\text{3}}\end{array}$

Then the elastic wave speed $c$ in the piezoelectric material is equal to $3684.06\text{m}/\text{s}$ . 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 [19] . Suppose that this force $F_a$ is applied downward to a piezoelectric cylinder with a length of $h=1\text{cm}$ and a diameter $d=1\text{cm}$ . So, the area of the upper end face of the cylinder $A=\text{π}{d}^2/4=0.785{\text{cm}}^{\text{2}}$ . Assuming that $F_a$ is uniformly distributed over the upper end face, we get the amplitude of the pressure impulse acting on the top of the cylinder: $p_a=364.66\text{MPA}$ . Let us assume that the applied normal stress load takes the form of the following rectangular compressive load (pressure) impulse:

$p\left(t\right)=-p_a\left[H\left(t\right)-H\left(t-t_1\right)\right]$ (46)

where $t_1$ is the duration of the pressure impulse.

We assume first that $\alpha =0.5$ and $k_{\alpha }=1000\text{N}\cdot {\text{s}}^{\text{1}/\text{2}}\cdot {\text{m}}^{-\text{1}/\text{2}}$ (see, e.g., [20] ) in the damper boundary conditions (23). Consider the following three values of $t_1:5,10,15\text{μs}$ . Since the transit time of elastic waves between the upper and lower end faces $\theta =2.71\text{μs}$ , then these three durations correspond to the three cases of explicit exact solutions considered in 4. The operator $Q_{\alpha }$ is calculated according to (37). The calculated results for the output voltage $V=\Delta \varphi$ (in kV) based on these exact solutions are presented for $t\le 50\text{μs}$ 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 $t_1$ . However, the number of peaks in each figure depends on the $t_1$ . The time distance between two neighboring peaks is approximately equal to $2\theta$ . 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: $\alpha =2.0$ and $k_{\alpha }=250\text{N}\cdot {\text{s}}^{\text{2}}\cdot {\text{m}}^{-\text{2}}$ . The parameters of the material and the impulse durations are the same as above. The operator $Q_{\alpha }$ is also calculated according to (37).

Figure 1. Output voltage for $\alpha =0.5,k_{\alpha }=1000\text{N}\cdot {\text{s}}^{\text{1}/\text{2}}\cdot {\text{m}}^{-\text{1}/\text{2}}$ and $t_1=5\text{μs}$ .

Figure 2. Output voltage for $\alpha =0.5,k_{\alpha }=1000\text{N}\cdot {\text{s}}^{\text{1}/\text{2}}\cdot {\text{m}}^{-\text{1}/\text{2}}$ and $t_1=10\text{μs}$ .

Figure 3. Output voltage for $\alpha =0.5,k_{\alpha }=1000\text{N}\cdot {\text{s}}^{\text{1}/\text{2}}\cdot {\text{m}}^{-\text{1}/\text{2}}$ and $t_1=15\text{μs}$ .

The calculated results for the output voltage $V=\Delta \varphi$ (in kV) are presented for $t\le 50\text{μs}$ 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 $t_1$ which similar to the case when $\alpha =0.5$ . The difference is that now there is only one peak but its width depends on the $t_1$ . 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 $\alpha =2,k_{\alpha }=250\text{N}\cdot {\text{s}}^{\text{2}}\cdot {\text{m}}^{-\text{2}}$ and $t_1=5\text{μs}$ .

Figure 5. Output voltage for $\alpha =2,k_{\alpha }=250{\text{Ns}}^{\text{2}}{\text{m}}^{-\text{2}}$ and $t_1=10\text{μs}$ .

6. Conclusion

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 $\alpha =2,k_{\alpha }=250\text{N}\cdot {\text{s}}^{\text{2}}\cdot {\text{m}}^{-\text{2}}$ and $t_1=15\text{μs}$ .

end faces of a piezoelectric cylinder fixed to a nonlinear damper at one end, and subjected to impulsive loading at the other.

Acknowledgements

The authors would like to thank the anonymous reviewers for their comments.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Mason, W.P. (1981) Piezoelectricity, Its History and Applications. The Journal of the Acoustical Society of America, 70, 1561. https://doi.org/10.1121/1.387221
[2]	Yang, J. (2006) Analysis of Piezoelectric Devices. World Scientific, Singapore. https://doi.org/10.1142/6156
[3]	Fang, D., Wang, J. and Chen, W., Eds. (2013) Analysis of Piezoelectric Structures and Devices. De Gruyter, Berlin, Boston.
[4]	Vijaya, M.S. (2013) Piezoelectric Materials and Devices: Applications in Engineering and Medical Sciences. CRC Press, Boca Raton.
[5]	Genis, V. and Soukhomlinoff, A. (2010) United States Patent: 7679271—Piezoelectric Powered Vehicles and Motors, March 2010.
[6]	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. https://doi.org/10.1109/IRENC.1956.199155
[7]	Filipczynski, L. (1960) Transients and the Equivalent Electrical Circuit of the Piezoelectric Transducer. Acta Acustica united with Acustica, 10, 149-154.
[8]	Redwood, M. (1961) Transient Performance of a Piezoelectric Transducer. The Journal of the Acoustical Society of America, 33, 527-536. https://doi.org/10.1121/1.1908709
[9]	Stuetzer, O.M. (1967) Multiple Reflections in a Free Piezoelectric Plate. The Journal of the Acoustical Society of America, 42, 502. https://doi.org/10.1121/1.1910607
[10]	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. https://doi.org/10.1121/1.390034
[11]	Le, K.C. (1999) Vibrations of Shells and Rods. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-59911-8
[12]	Ma, C.C., Chen, X.H. and Ing, Y.S. (2007) Theoretical Transient Analysis and Wave Propagation of Piezoelectric Bi-Materials. International Journal of Solids and Structures, 44, 7110-7142.
[13]	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.
[14]	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. https://doi.org/10.1007/s00419-015-1042-3
[15]	Khutoryansky, N.M. and Sosa, H. (1995) Dynamic Representation Formulas and Fundamental Solutions for Piezoelectricity. International Journal of Solids and Structures, 32, 3307-3325.
[16]	Yuan, F.G. (2016) Structural Health Monitoring (SHM) in Aerospace Structures. Elsevier Science, Amsterdam.
[17]	Dixon, J. (2008) The Shock Absorber Handbook. John Wiley & Sons, Hoboken.
[18]	Monagan, M.B., Geddes, K.O., Heal, K.M., Labahn, G., Vorkoetter, S.M., McCarron, J. and DeMarco, P. (2005) Maple 10 Programming Guide. Maplesoft, Waterloo, ON, Canada.
[19]	Pulkrabek, W.W. (1997) Engineering Fundamentals of the Internal Combustion Engine. Prentice Hall, Upper Saddle River.
[20]	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