Existence, Uniqueness and Regularity of Piezoelectric Partial Differential Equations

TL;DR

This paper investigates the mathematical properties of the PDEs governing piezoelectric materials, establishing existence, uniqueness, regularity, and long-term behavior of solutions to improve understanding and simulation accuracy.

Contribution

It provides rigorous proofs of well-posedness and regularity for the coupled PDE system modeling piezoelectric behavior, including effects of Rayleigh damping.

Findings

Proved existence and uniqueness of solutions

Established regularity and long-term behavior of solutions

Applied Galerkin approximation and energy estimates

Abstract

Piezoelectric appliances have become hugely important in the past century and computer simulations play an essential part in the modern design process thereof. While much work has been invested into the practical simulation of piezoelectric ceramics there still remain open questions regarding the partial differential equations governing the piezoceramics. The piezoelectric behavior of many piezoceramics can be described by a second order coupled partial differential equation system. This consists of an equation of motion for the mechanical displacement in three dimensions and a coupled electrostatic equation for the electric potential. Furthermore, an additional Rayleigh damping approach makes sure that a more realistic model is considered. In this work we analyze existence, uniqueness and regularity of solutions to theses equations and give a result concerning the long-term behavior. The well-posedness of the initial boundary value problem in a bounded domain with sufficiently smooth boundary is proved by Galerkin approximation in the discretized weak version, followed by an energy estimation using Gronwall inequality and using the weak limit to show the results in the infinite dimensional space. Initial conditions are given for the mechanical displacement and the velocity.

Equations

$$H_{0,\Gamma }^1(\Omega )$$

$$H_{\mathcal{B}}^1(\Omega )$$

$$\mathcal{B}:=\left(\begin{array}{ccc}\frac{\partial }{\partial x}& 0& 0\\ 0& \frac{\partial }{\partial y}& 0\\ 0& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial z}& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}& 0\end{array}\right)$$

$$\mathcal{N}:=\left(\begin{array}{ccc}{n}_x& 0& 0\\ 0& n_y& 0\\ 0& 0& n_z\\ 0& n_z& n_y\\ n_z& 0& n_x\\ n_y& n_x& 0\end{array}\right).$$

$$\displaystyle c^{E}:=\left(\begin{array}[]{cccccc}c_{11}&c_{12}&c_{13}&0&0&0\\ c_{12}&c_{11}&c_{13}&0&0&0\\ c_{13}&c_{13}&c_{33}&0&0&0\\ 0&0&0&c_{44}&0&0\\ 0&0&0&0&c_{44}&0\\ 0&0&0&0&0&\frac{1}{2}(c_{11}-c_{12})\\ \end{array}\right)\in\mathbb{R}^{6\times 6}$$

$$\displaystyle\epsilon^{S}:=\left(\begin{array}[]{ccc}\epsilon_{11}&0&0\\ 0&\epsilon_{11}&0\\ 0&0&\epsilon_{33}\end{array}\right)\in\mathbb{R}^{3\times 3}$$

$$\displaystyle e:=\left(\begin{array}[]{cccccc}0&0&0&0&e_{15}&0\\ 0&0&0&e_{15}&0&0\\ e_{13}&e_{13}&e_{33}&0&0&0\\ \end{array}\right)\in\mathbb{R}^{6\times 3}.$$

$$\rho \ddot{\mathbf{u}}(t)+\alpha \rho \dot{\mathbf{u}}(t)-{\mathcal{B}}^T\left(c^E\mathcal{B}\mathbf{u}(t)+\beta c^E\mathcal{B}\dot{\mathbf{u}}(t)+e^T\nabla \varphi (t)\right)$$

$$-\nabla \cdot \left(e\mathcal{B}\mathbf{u}(t)-{ϵ}^S\nabla \varphi (t)\right)$$

$$\varphi (t)$$

$$\vec{n}\cdot \left(e\mathcal{B}\mathbf{u}(t)-{ϵ}^S\nabla \varphi (t)\right)$$

$${\mathcal{N}}^T\left(c^E\mathcal{B}\mathbf{u}(t)+\beta c^E\mathcal{B}\dot{\mathbf{u}}(t)+e^T\nabla \varphi (t)\right)$$

$$\mathbf{u}(0)$$

$$\dot{\mathbf{u}}(0)$$

$${\int }_{\Omega }{\left({\mathcal{B}}^T\sigma \right)}^{T}vd\Omega =-{\int }_{\Omega }{\sigma }^T\mathcal{B}vd\Omega +{\int }_{\partial \Omega }{\left({\mathcal{N}}^T\sigma \right)}^{T}vd\Omega .$$

$${\int }_{\Omega }\rho {\ddot{\mathbf{u}}}^{T}vd\Omega +\alpha {\int }_{\Omega }\rho {\dot{\mathbf{u}}}^{T}vd\Omega +{\int }_{\Omega }{\left(c^E\mathcal{B}\mathbf{u}\right)}^T\mathcal{B}vd\Omega +\beta {\int }_{\Omega }{\left(c^E\mathcal{B}\dot{\mathbf{u}}\right)}^T\mathcal{B}vd\Omega$$

$$+$$

$$=$$

$$ab\le \frac{a^p}{p}+\frac{b^q}{q}.$$

$${\int }_{\Omega }|uv|dx\le \Vert u{\Vert }_{L^p(\Omega )}\Vert v{\Vert }_{L^q(\Omega )}.$$

$$\eta (t)\le C_1{\int }_0^t\eta (s)ds+C_2$$

$$\eta (t)\le C_2\left(1+C_{1}t{e}^{C_{1}t}\right)$$

$$\eta (t)\le C_1{\int }_0^t\eta (s)ds$$

$$\eta (t)=0\text{ a.e. }$$

$$(\mathbf{u},\varphi )\in L^{\infty }(0,T;H_{\mathcal{B}}^1(\Omega ))\times L^{\infty }(0,T;H_{0,\Gamma }^1(\Omega ))$$

$$\dot{\mathbf{u}}\in L^{\infty }(0,T;L^2(\Omega ))\text{ and }\ddot{\mathbf{u}}\in L^2(0,T;(H_{\mathcal{B}}^1(\Omega ){)}^{\prime })$$

$$\mathbf{u}(0)={\mathbf{u}}_0,\dot{\mathbf{u}}(0)={\mathbf{u}}_1\text{ on }\Omega$$

$$\Vert \mathbf{u}{\Vert }_{L^{\infty }(0,T;H_{\mathcal{B}}^1(\Omega ))}+\Vert \dot{\mathbf{u}}{\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}+\Vert \ddot{\mathbf{u}}{\Vert }_{L^2(0,T;(H_{\mathcal{B}}^1(\Omega ){)}^{\prime })}+\Vert \varphi {\Vert }_{L^{\infty }(0,T;H_{0,\Gamma }^1(\Omega ))}$$

$$\le$$

$$⟨\rho \ddot{\mathbf{u}}(t),v⟩+\alpha ⟨\rho \dot{\mathbf{u}}(t),v⟩+⟨c^E\mathcal{B}\mathbf{u}(t),\mathcal{B}v⟩+\beta ⟨c^E\mathcal{B}\dot{\mathbf{u}}(t),\mathcal{B}v⟩$$

$$+$$

$$=$$

$$⟨f(t),v⟩:=-{\varphi }^e(t){\int }_{\Omega }(e^T\nabla \chi {)}^T\mathcal{B}vd\Omega ,$$

$$⟨g(t),w⟩:={\varphi }^e(t){\int }_{\Omega }({ϵ}^S\nabla \chi {)}^T\nabla wd\Omega .$$

$$\Vert f{\Vert }_{H^1(0,T;(H_{\mathcal{B}}^1(\Omega ){)}^{\prime })}\le c_1(\Omega )\Vert {\varphi }^e{\Vert }_{H^1(0,T)},$$

$$\Vert g{\Vert }_{H^1(0,T;H^{-1}(\Omega ))}\le c_2(\Omega )\Vert {\varphi }^e{\Vert }_{H^1(0,T)}.$$

Full text

Existence, Uniqueness and Regularity of Piezoelectric Partial Differential Equations

Benjamin Jurgelucks1

Veronika Schulze1

Tom Lahmer2

(1Paderborn University, 2Bauhaus-University Weimar)

Abstract

Piezoelectric appliances have become hugely important in the past century and computer simulations play an essential part in the modern design process thereof. While much work has been invested into the practical simulation of piezoelectric ceramics there still remain open questions regarding the partial differential equations governing the piezoceramics.

The piezoelectric behavior of many piezoceramics can be described by a second order coupled partial differential equation system. This consists of an equation of motion for the mechanical displacement in three dimensions and a coupled electrostatic equation for the electric potential. Furthermore, an additional Rayleigh damping approach makes sure that a more realistic model is considered.

In this work we analyze existence, uniqueness and regularity of solutions to theses equations and give a result concerning the long-term behavior. The well-posedness of the initial boundary value problem in a bounded domain with sufficiently smooth boundary is proved by Galerkin approximation in the discretized weak version, followed by an energy estimation using Gronwall inequality and using the weak limit to show the results in the infinite dimensional space. Initial conditions are given for the mechanical displacement and the velocity.

1 Introduction

Piezoelectricity has become more and more important for technical purposes and innovations especially when high-frequency vibrations are to be measured or produced. Typical applications as actuators range from piezo-igniters over ultrasonic toothbrushes to diesel fuel injectors as well as many others, e.g., as part of intelligent sensory equipment. The piezoelectric effect describes the transaction between electrical and mechanical energy changes of a piezoelectric sample. The effect is caused by the structure of the material and its polarization. Therefor it is clear, that the effect and its usage is material based (cf. [6]).

There are two problems which can be solved regarding the piezoelectric equations, the forward and the inverse problem. For details regarding the inverse problem and optimization of sensitivities see e.g. [5], [7].

In order to design and analyse new piezoelectric devices, models are employed [14]. However for a reliable use existence, uniqueness and regularity of the solutions for these models need to be guaranteed.

The underlying application of the well-posedness result is a piezoelectric ceramic disc with top and bottom surface electrodes. The material parameters are extracted from real measurements for the forward simulation to compute the mechanical displacement and the electrical potential after electrical excitation.

The proof of the properties mentioned above assumes a bounded domain with sufficiently smooth boundary. Our piezoceramic and the electrodes on top and bottom fulfill these requirements.

The underlying model is linear, includes Rayleigh damping and neglects thermal effects. The behavior of the piezoelectric material can be described by a second order partial differential equation system, which defines the mechanical displacement and the electrical potential. By an appropriate choice of the Rayleigh damping parameters, the equation of motion of the mechanical displacement is a hyperbolic partial differential equation and the electrostatic equation of the electrical potential is an elliptic partial differential equation. The density, the elastic stiffness, the dielectric permittivity and the piezoelectric coupling matrices are the given material components in the standard Voigt notation. There are several existing works on the well-posedness of the piezoelectric initial boundary problem usually without any damping models. The proof structures used in this paper are similar. Parts of our work are based on the proof presented in [10]. Technical details are however elaborated in more detail and some derivations are developed in a more rigorous way. Proofs for the static and harmonic case can be found in [8] and [10].

The proof is divided in four general steps. First, the system is transformed into the weak form and discretized, via Galerkin approximation. Then, via standard theory for ordinary differential equations there exist unique solutions. The finiteness of the finite dimensional solution is shown by the energy estimates via the Gronwall inequality. The weak limit of the discretized solutions provide the weak existence of a solution in infinite dimensional function spaces. The uniqueness of the solution is shown by applying the estimates to the homogeneous system and getting the trivial solution.

In the second part of the paper, Theorem 36 studies higher regularities for the solution of the system based on higher regularity requirements for the initial condition of the mechanical displacement, the velocity and the boundary value for the electrical potential. Finally, a remark about the long-term behaviour of an energy functional considered in the proof of Theorem 7 is stated.

2 Setting

Before we can begin to solve any partial differential equation we must first establish an exact setup - the geometry $\Omega$, the boundary $\partial\Omega$, the boundary conditions and initial values of the partial differential equations in question. We consider the case of a mechanically unclamped piezoceramic which is excited by prescribing a voltage on a part of the boundary. Let $\Omega\subseteq\mathbb{R}^{3}$ be an open domain describing the piezoelectric ceramic and let $\partial\Omega=:\Gamma$ be the nonempty boundary of $\Omega$. The boundary is divided into nonempty, disjunct, covering subsets of $\Gamma$ (see also Fig. 1) which are assumed to have a positive 2D measure. Let $\Gamma_{e}$ be the section of the boundary which is electrically excited, $\Gamma_{g}$ the section of the boundary which is grounded, $\Gamma_{r}=\Gamma\setminus\left(\Gamma_{e}\cup\Gamma_{g}\right)$ the remaining boundary section.

For the readers convenience the usual definitions of common function spaces which will be required later on are stated in the appendix A. Only the newly defined function spaces for the considered differential equation system are described now:

[TABLE]

where

[TABLE]

with $x,y,z$ refering to Cartesian coordinates. In this paper we denote derivatives with respect to time by the dot symbol e.g. $\dot{\sigma}$ and derivatives with respect to space by the nabla or $\mathcal{B}$ symbol, e.g. $\nabla\sigma$ or $\mathcal{B}\sigma$. Here $\mathcal{B}$ denotes the symmetric gradient in Voigt notation. It should be noted that the last three entries of the matrix vector product $\mathcal{B}u$ still contains the factor $2$, but for simplicity, no attention is paid here. The factor can be included in the definition of the linear strain vector $S$, where $S=\mathcal{B}u$.

All derivatives in the above are understood in the distributional sense. In addition, the dual space of a Hilbert space $X$ is denoted by $X^{\prime}$. In particular, $H^{-1}(\Omega)$ denotes the dual space of $H^{1}_{0}(\Omega)$. Note that in order to simplify the notation superscripts indicating the dimension of $\mathbf{u}$ or $\mathcal{B}\mathbf{u}$, which are 3 and 6 respectively, are omitted. This is reasonable as the vectorial scalar product inside $\int_{\Omega}\sigma^{T}\sigma\,d\Omega$ always returns a scalar no matter what dimensions $\sigma$ has.

Let $\vec{n}:=(n_{x},\,n_{y},\,n_{z})$ be the normal vector and

[TABLE]

Definition 1**.**

The material parameters $c^{E},\epsilon^{S}$ and $e$ $([c^{E}]=N\cdot m^{-2},\ [\epsilon^{S}]=F\cdot m^{-1},\ [e]=C\cdot m^{-2})$ are given by

[TABLE]

[TABLE]

[TABLE]

The material parameters are said to fulfill nonnegativity conditions if $c^{E}$ and $\epsilon^{S}$ are positive definite matrices.

The three dimensional transient linear piezoelectric equations with Rayleigh damping parameters $\alpha,\beta>0$ (chosen sufficiently large enough so that the system is parabolic) and density $\rho>0$ describing the mechanical displacement $\mathbf{u}\in\mathbb{R}^{3}$ and the electrical potential $\phi\in\mathbb{R}$ with given boundary conditions are stated as:

[TABLE]

The weak form of the equations above can easily be obtained [10] by testing with appropriate functions $v\in\mathbb{R}^{3}$ (for the first line) and $w\in\mathbb{R}$ (for the second line), integration by parts and using boundary conditions:

[TABLE]

First, we use a Dirchlet lift ansatz to homogenize the Dirichlet boundary condition for $\phi(t)$: Let $t\in[0,T]$ and let $\chi\in H^{1}(\Omega)$ where $\chi|_{\Gamma_{g}}=0$ and $\chi|_{\Gamma_{e}}=1$. Such a $\chi$ exists if we assume that $\Omega$ is at least a Lipschitz domain. Let $\phi(t)$ consist of two parts $\phi(t)=\phi_{0}(t)+\phi_{\phi^{e}}(t)$ where $\phi_{0}(t)\in H^{1}_{0}(\Omega)$ and $\phi_{\phi^{e}}(t)\in H^{1}(\Omega)$. We then rewrite $\phi_{\phi^{e}}(t)=\phi^{e}(t)\chi$. Therefor we set $\phi_{0}(t):=\phi(t)-\phi^{e}(t)\chi$.

As $\phi^{e}(t)$ is a given value $\phi^{e}(t)\chi$ can be taken out of the left hand side of the weak form and added to the right hand side. The weak form of the piezoelectric system for all $t\in[0,T]\text{ a.e.}$ and for all test functions $(v,w)\in H^{1}_{\mathcal{B}}(\Omega)\times H^{1}_{0}(\Omega)$ is given by

[TABLE]

Note that in light of [4, Thm. 2 in section 5.9.2] it makes sense to demand $\mathbf{u}(0)=\mathbf{u}_{0}$ and $\mathbf{\dot{u}}(0)=\mathbf{u}_{1}$. See also the only remark in [4, section 7.2.1].

3 Existence, uniqueness and regularity of solutions

Before we attempt to show existence, uniqueness and regularity of solutions some additional tools are required:

Lemma 2**.**

*(Young inequality)

Let $1<p,q<\infty,\frac{1}{p}+\frac{1}{q}=1$. Then for $a,b>0$ the following inequality holds:*

[TABLE]

Proof.

See [4, Appendix B.2]. ∎

Lemma 3**.**

*(Hölder inequality)

Let $1\leq p,q\leq\infty,\frac{1}{p}+\frac{1}{q}=1$. Then for $u\in L^{p}(\Omega),v\in L^{q}(\Omega)$ the following inequality holds:*

[TABLE]

Proof.

See [4, Appendix B.2]. ∎

Remark 4**.**

The latter two inequalities are especially true for $p=q=2$. The latter inequality is then known as Cauchy–Schwarz (C.S.) inequality.

Lemma 5**.**

*(Gronwall inequality, integral form)

1. a)

Let $\eta:[0,T]\rightarrow\mathbb{R}^{\geq 0}$ be a nonnegative, summable function on $[0,T]$, which satisfies for almost every $t$ the differential inequality

[TABLE]

for constants $C_{1},C_{2}\geq 0$. Then

[TABLE]

for a.e. $0\leq t\leq T$. 2. b)

In particular, if

[TABLE]

for a.e. $t\in[0,T]$, then

[TABLE]

Proof.

See [4, Appendix B.2]. ∎

Remark 6**.**

(Sufficiently smooth boundary)

We say the boundary $\partial\Omega$ is sufficiently smooth if it permits application of the trace theorem (cf. [4]).

Thus, a $C^{1}-$boundary is sufficient. However, it is possible to utilize a variation of the trace theorem under less strict requirements (cf. [3]). We note that the boundary for our specific application (see Fig. 1) satisfies the special Lipschitz condition stated in Definition 5 of [3] and thus it appears that it can also be considered sufficiently smooth.

A proof for the following theorem was first given in [10]. The proof given there is also heavily oriented on work of [2] which itself is based on [13]. Here we present a proof with the similar essential steps as in the other given proofs, but with more technical details and necessary exact descriptions.

The proof follows the usual guideline as seen for many partial differential equations (e.g. [4, p. 353]): We get existence and uniqueness of a weak solution by the usual procedure:

Discretization via Galerkin approximation of infinite dimensional function spaces, 2. 2.

energy estimates via Gronwall inequality in discretized space which provide finiteness of the discretized solution, 3. 3.

weak limit of discretized solution provides weak existence of a solution in infinite dimensional function space, 4. 4.

uniqueness of the solution is shown by applying estimates to the difference $w:=w_{1}-w_{2}$ of two solutions $w_{1}$ and $w_{2}$. Thus, the only solution to the homogeneous case is the trivial solution.

Theorem 7**.**

*Let $\Omega\subseteq\mathbb{R}^{3}$ be a bounded domain with sufficiently smooth boundary as specified in Remark 6. Let the real valued material parameters $c^{E},e$ and $\epsilon^{S}$ be defined as in Def. 1 and let $c^{E}$ and $\epsilon^{S}$ be symmetric and positive definite. The Rayleigh coefficients $\alpha$ and $\beta$ are assumed to be nonnegative. Let $T>0$ and $\rho>0$.

Then there exists a $C>0$ such that for any $\mathbf{u}_{0}\in H^{1}_{\mathcal{B}}(\Omega),\mathbf{u}_{1}\in L^{2}(\Omega)$ and $\phi^{e}\in H^{1}(0,T;H^{1/2}(\Gamma_{e}))$ there exists a unique solution*

[TABLE]

with

[TABLE]

to Eq. (4) satisfying the initial conditions

[TABLE]

and the following estimate holds:

[TABLE]

Proof.

Note that many concepts of this proof are taken from [4, chapter 7] and information regarding involved spaces can be found in [1].

In the following constants denoted by the letter $C$ or $\tilde{C}$ are used. Unless explicitly specified otherwise we note that all these constants are positive $C_{i}>0\,,i\geq 1$.

Weak solutions are functions $\mathbf{u},\mathbf{\dot{u}},\mathbf{\ddot{u}}$ and $\phi_{0}$ as in Eq. (5) and Eq. (6) where $\phi_{0}=\phi+\phi_{\phi_{e}}$ such that for almost all $t\in[0,T]$ for all $(v,w)\in H^{1}_{\mathcal{B}}(\Omega)\times H^{1}_{0}(\Omega)$ the following equation holds:

[TABLE]

with

[TABLE]

and

[TABLE]

Note that by the Riesz representation theorem there exists a unique representation for the latter functionals as an inner product, i.e. $\left<f,\cdot\right>$ and $\left<g,\cdot\right>$. As is common in the field of partial differential equation for convenience we will also use the same symbols $f$ and $g$ to refer to the Riesz-representative as well as the functionals $\left<f,\cdot\right>$ and $\left<g,\cdot\right>$. Furthermore, we remember that $\chi\in H^{1}(\Omega)$ and that $\epsilon^{S}$, $e$ are constant. The integrals of the right hand side $\int_{\Omega}(e^{T}\nabla\chi)^{T}\mathcal{B}v\,d\Omega$, $\int_{\Omega}(\epsilon^{S}\nabla\chi)^{T}\nabla w\,d\Omega$ are finite, their values $c_{1}(\Omega),c_{2}(\Omega)<\infty$ depend, e.g.,  only on domain $\Omega$ but not on time $t$. Thus, by integrating this constant value over time we can estimate the Bochner-space norm of $f$ by

[TABLE]

and analogously we get

[TABLE]

**Phase 1: Galerkin approximation

**The weak form is tested with test functions ${v}_{j}\in H^{1}_{\mathcal{B}}(\Omega)$ and ${w}_{j}\in H^{1}_{0}(\Omega),\quad j\in\mathbb{N}$, with

[TABLE]

and

[TABLE]

where ’$\approx$’ is to be understood in the sense of an orthogonal projection in the appropriate spaces. The finite dimensional spaces spanned by the test functions are defined as

[TABLE]

We can assume that the dimension of the test function spaces $\operatorname{dim}(V_{m})=\operatorname{dim}(W_{m})=m$ are the same, for $V_{m}$ in each vectorial component. So the test functions can be selected to be linearly independent. Furthermore, the functions can be chosen such that

[TABLE]

Then via standard theory for ordinary differential equations (see e.g. [4] or [12]) for all $m\in\mathbb{N}$ and for all $({v}_{m},w_{m})\in V_{m}\times W_{m}$ there exists a unique solution

[TABLE]

to the discretized version of Eq. (8) that fulfills the initial conditions $\mathbf{u}_{m}(0)=(\mathbf{u}_{0})_{m},$ $\mathbf{\dot{u}}_{m}(0)=(\mathbf{u}_{1})_{m}$. For more information on Sobolev spaces involving time and space see also [4, section 5.9.2].

**Phase 2: Energy estimates

**The aim of this phase is to use Gronwall inequality to show an energy estimate from which the finiteness of the finite dimensional solutions $(\mathbf{u}_{m}(t),\phi_{m}(t))$ in $L^{\infty}(0,T;H^{1}_{\mathcal{B}}(\Omega))\times L^{\infty}(0,T;H^{1}_{0}(\Omega))$, $\mathbf{\dot{u}}_{m}$ in $L^{\infty}(0,T;L^{2}(\Omega))$ and $\mathbf{\ddot{u}}_{m}$ in

$L^{2}(0,T;(H^{1}_{\mathcal{B}}(\Omega))^{\prime})$ can be deduced:

Let

[TABLE]

In order to use the Gronwall inequality we must show that there are constants $p,q\geq 0$ such that $\eta(t)\leq p\int_{0}^{t}\eta(s)\,ds+q$ holds. If this condition is true, then it can be shown that

[TABLE]

holds almost everywhere in $[0,T]$. Thus, this must also be true for the essential supremum over $0\leq t\leq T$ and we will get finiteness in the $L^{\infty}(0,T;X)$ norm for the appropriate sub-spaces $X$. In order to show the requirement we consider the following:

First, the discretized version of the weak form Eq. (8) is supposed to hold for all test functions $({v}_{m},w_{m})$. Thus, it should also hold for $(\mathbf{\dot{u}}_{m}(t),0)$:

[TABLE]

By transposing the inner product and direct computation it is easy to see that one can swap the placement of constant scalars or matrices such as $\rho,\epsilon^{S},c^{E}$ etc. (which are symmetric) in this bilinear form, e.g. the following holds:

[TABLE]

Thus, by bilinearity of the inner product

[TABLE]

Hence, the above equation simplifies to

[TABLE]

Now we differentiate the weak form Eq. (8) with respect to $t$ and test it with $(0,\phi_{m}(t))$, taking into account that the test functions ${v}_{m},w_{m}$ do not depend on time $t$, therefor the time derivatives $\dot{{v}}_{m},\dot{w}_{m}\equiv 0$:

[TABLE]

A subtraction of Eq. (10) and Eq. (11) gives

[TABLE]

The last equation Eq. (12) is integrated with respect to $t$.

[TABLE]

Hence, in short we can write

[TABLE]

Now the aim is to use this equation to show that the requirements for the Gronwall inequality are met.

We start by showing that the left-hand side $\mathcal{F}_{l}(t)$ of Eq. (13) has a lower bound. With $\lambda_{1,mech}$ the smallest eigenvalue of $c^{E}$ (which is strictly positive) one estimates

[TABLE]

With $\lambda_{1,elec}$ the smallest eigenvalue of $\epsilon^{S}$ (which is strictly positive) one estimates

[TABLE]

From the Poincaré inequality (see e.g. [15]), we obtain $c_{1},c_{2}\in\mathbb{R}$ such that

[TABLE]

By nonnegativity of $\rho,\alpha,\beta$ and the two inequalities Eq. (14) and Eq. (15) one can now estimate

[TABLE]

with a positive constant $C_{1}>0$. Furthermore, by the inequalities Eq. (14) and Eq. (15) and Cauchy–Schwarz and Young inequalities the right hand side $\mathcal{F}_{r}(t)$ can be bounded from above with $c,\tilde{c}>0$:

[TABLE]

Hence, we get

[TABLE]

with a positive constant $C_{2}>0$. As $\mathcal{F}_{l}(t)=\mathcal{F}_{r}(t)$ it is now clear that

[TABLE]

Utilizing the inequality $\|u(t)\|_{L^{2}(\Omega)}^{2}\leq 2\|u(0)\|_{L^{2}(\Omega)}^{2}+2T\int_{0}^{t}\|\dot{u}(s)\|^{2}_{L^{2}(\Omega)}\,ds$ (see [17] p. 425) for $T$ large enough, we can remove $c_{mech}\|\mathbf{u}_{m}(t)\|^{2}_{L^{2}(\Omega)}$ from the left hand side of the inequality to obtain:

[TABLE]

where $C_{4}>0$ now also depends on the fixed value $T$.

Let

[TABLE]

and let

[TABLE]

Then the above inequality simplifies to

[TABLE]

Hence, all requirements for Gronwall inequality have been shown to hold and it can now be safely applied and the result simplified to:

[TABLE]

holds almost everywhere in $[0,T]$.

We will return to this inequality shortly after considering the bilinear form

[TABLE]

and the continuous linear functional on $H^{1}_{0}(\Omega)$ for a fixed $\mathbf{u}_{m}(t)$

[TABLE]

which together form the weak form Eq. (8) tested by $(0,w)$. This bilinear form $A$ is coercive (inequality Eq. (15)) and continuous:

[TABLE]

Using the Lax–Milgram lemma and the Young inequality we get the estimate for $A(\phi_{m}(t),w)=b(w)$ $\forall w\in H^{1}_{0}(\Omega)$:

[TABLE]

Furthermore, for $t=0$ we get

[TABLE]

Hence, we obtain

[TABLE]

Finally, from the Gronwall inequality we can thus deduce

[TABLE]

Now knowing that all these values are finite we can deduce from Eq. (13) with $\beta>0$ that also

[TABLE]

It now remains to show that $\|\mathbf{\ddot{u}}_{m}\|_{L^{2}(0,T;(H^{1}_{\mathcal{B}}(\Omega))^{\prime})}$ is finite. We follow the general guideline given in e.g. [4, p. 355].

Fix any $\tilde{v}\in H^{1}_{\mathcal{B}}(\Omega)$ with $\|\tilde{v}\|_{H^{1}_{\mathcal{B}}(\Omega)}\leq 1$ and $\tilde{v}:=\tilde{v}^{1}+\tilde{v}^{2}$ with $\tilde{v}^{1}\in\operatorname{span}\{{v}_{i}\}_{i=1}^{m}$ and $\left<\tilde{v}^{2},{v}_{i}\right>=0$ for all $1\leq i\leq m$. Since $\{{v}_{i}\}_{i=0}^{m}$ can be assumed orthogonal in $H^{1}_{\mathcal{B}}(\Omega)$,

[TABLE]

Now with $\mathbf{u}_{m}=\sum_{i=0}^{m}u_{m}^{i}(t){v}_{i}$ the following holds almost everywhere in $[0,T]$:

[TABLE]

where the subscript $H^{1}_{\mathcal{B}}(\Omega)$ denotes the duality pairing between $\left(H^{1}_{\mathcal{B}}(\Omega)\right)^{\prime}$ and $H^{1}_{\mathcal{B}}(\Omega)$.

Using the Cauchy–Schwarz inequality we can deduce

[TABLE]

Using the Lax–Milgram lemma again on the form Eq. (20) for an arbitrary $t\in[0,T]$ we can further deduce with analogous arguments as in Eq. (21) that the following holds

[TABLE]

Thus by repetitive application of the Young inequality we get for the norm

[TABLE]

Now we have finiteness for all components, hence we can finally integrate inequality Eq. (25) over $[0,T]$.

We rearrange the terms and apply the estimates Eq. (23) and Eq. (24).

[TABLE]

Thus, it is now clear that

[TABLE]

**Phase 3: Weak limit

**Following e.g. [4, p. 384], [15, p. 239] from the energy estimates Eq. (23) and Eq. (26) we get the boundedness of the sequences

[TABLE]

Thus there exist subsequences

[TABLE]

with $(\mathbf{u},\phi)\in L^{\infty}(0,T;H^{1}_{\mathcal{B}}(\Omega))\times L^{\infty}(0,T;H^{1}_{0,\Gamma}(\Omega))$, $\mathbf{\dot{u}}\in L^{\infty}(0,T;L^{2}(\Omega))$ and $\mathbf{\ddot{u}}\in L^{2}(0,T;(H^{1}_{\mathcal{B}}(\Omega))^{\prime})$ such that

[TABLE]

We now proceed to show that the weak limit is a solution of the weak form. Following [4, p. 384] we fix a $N\in\mathbb{N}$ and choose functions ${v}\in C^{1}(0,T;H^{1}_{\mathcal{B}}(\Omega))$ and $w\in C^{1}(0,T;H^{1}_{0}(\Omega))$ having the form

[TABLE]

We choose $m\geq N$, multiply the discretized versions for each pair $({v}_{k},w_{k})$ of the weak form Eq. (8) with $(u^{k}_{m}(t),\phi^{k}_{m}(t))$, sum over $k=1,\dots,N$, integrate with respect to $t$. This yields

[TABLE]

Fixing $m=m_{l}$ and using Eq. (28) we obtain in the limit $m\rightarrow\infty$ along the subsequence $m_{l}$

[TABLE]

Noting that all functions of form Eq. (29) are dense in the according spaces this equality holds for all functions ${v}\in L^{2}(0,T;H^{1}_{\mathcal{B}}(\Omega)),\,w\in L^{2}(0,T;H^{1}_{0}(\Omega))$. In particular it follows that also

[TABLE]

almost everywhere $t\in[0,T]$ for all ${v}\in H^{1}_{\mathcal{B}}(\Omega)$ and $w\in H^{1}_{0}(\Omega)$.

Following [10] and [4, section 7.2.2 Thm. 3, p. 385] we confirm that the initial conditions are also met. Choose any function $({v},0)$ with ${v}\in C^{2}([0,T];H^{1}_{\mathcal{B}}(\Omega))$ and ${v}(T)=\dot{{v}}(T)=0$. By integrating by parts twice with respect to $t$ of Eq. (30) we get

[TABLE]

and analogously using Eq. (31) we get

[TABLE]

For Eq. (33) we set $m=m_{l}$ and recall Eq. (28) to deduce

[TABLE]

By equating coefficients of Eq. (34) and Eq. (35) (set either ${v}(0)$ or ${\dot{v}}(0)$ to zero) we conclude $\mathbf{u}(0)=\mathbf{u}_{0}$ and $\mathbf{\dot{u}}(0)=\mathbf{u}_{1}$.

**Phase 4: Uniqueness

**

Following e.g. [4, p. 385] it suffices to show that the only weak solution with

[TABLE]

is

[TABLE]

Notice that by property Eq. (24) $\|\mathcal{B}\mathbf{\dot{u}}_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}$ is finite. Hence, the remark in [4, remark below Thm 4, section 7.2.2 c), p. 385] does not apply to our case and we can continue in the fashion of [4, Thm. 4, section 7.1.2c), p. 358] instead. Passing to limits, we substitute ${v}=\mathbf{u}$ and $w=\phi_{0}$ in the original weak form. This is not prohibited as by property Eq. (24) all components exist also in the limit. Hence, we can deduce that the following non-discretized inequality holds

[TABLE]

In the case $t=0$ we get from Eq. (22) that $\|\phi(0)\|^{2}_{H^{1}_{0}(\Omega)}=0$. Hence, we now note that

[TABLE]

Finally, we can apply the second part of the Gronwall inequality to conclude that

[TABLE]

Thus, the only solution can be the trivial solution. ∎

Through the theorem we know what requirements we need to get existence of a solution of the weak form. Now prerequisites can be derived to achieve higher regularities of the solutions.

The following theorem is inspired by Thm. 5, chapter 7.2 in [4]. The proof uses ideas from [12] adapted for additional Rayleigh damping.

Theorem 8**.**

Let all requirements of Thm. 7 hold. If additionally $\mathbf{u}_{0}\in H^{2}(\Omega)$, $\mathbf{u}_{1}\in H^{1}(\Omega)$, $\beta\mathbf{u}_{1}\in H^{2}(\Omega)$, $\phi^{e}\in H^{2}(0,T;H^{1/2}(\Gamma_{e}))$, then

[TABLE]

Proof.

We differentiate the weak form Eq. (8) once with respect to time $t$ and test the result first with $(\mathbf{\ddot{u}}_{m}(t),0)$ to obtain

[TABLE]

Then we differentiate the weak form Eq. (8) twice with respect to time $t$ and test the result first with $(0,{\dot{\phi}}_{m}(t))$ to obtain

[TABLE]

Analogously to the proof of Thm. 7 we subtract these two results and integrate with respect to $t$ to obtain in analogy to Eq. (13)

[TABLE]

or, again, abbreviated as $\mathcal{F}_{l}=\mathcal{F}_{r}$. Analogously to inequality Eq. (18) we then can obtain

[TABLE]

for some $C_{1},C_{3},C_{4}>0$. Note that by deriving the weak form which we then test by $(0,\phi_{m}(t))$ we additionally obtain a bilinear form similar to Eq. (20) and can analogously deduce with the Lax–Milgram lemma that

[TABLE]

This is only possible because of the added requirement of increased regularity of $\mathbf{\dot{u}}(0)$ and $\dot{g}$. Furthermore, by the additional requirements on $\mathbf{u}_{m}(0)\in H^{2}(\Omega)$ we also obtain (estimating the $H^{2}$ norm by the norm of the Laplacian, see e.g. [16])

[TABLE]

In order to utilize the Gronwall lemma we are left to show finiteness of

$\|\mathbf{\ddot{u}}_{m}(0)\|^{2}_{L^{2}(\Omega)}$. Notice that by the increased regularity of $\mathbf{u}_{m}(0)\in H^{2}(\Omega)$ and $\phi_{m}(0)\in H^{2}(\Omega)$ the weak solution is also a strong solution, i.e., not quite a classical solution but solves the classical equations in $t=0$ almost everywhere, see e.g. [11, section 2.3 and 3.5 ]. Thus, by evaluating the strong system in $t=0$ and using the initial data and previously deduced inequalities we obtain

[TABLE]

Note that this $f$ is given by the Dirichlet ĺift ansatz for the strong system. Therefor we choose $\chi\in H^{2}(\Omega)$ where $\chi|_{\Gamma_{g}}=0$ and $\chi|_{\Gamma_{e}}=1$. With this requirement the right hand sight of the above inequality is bounded independently of $m$.

Since all components are finite, analogously to inequality Eq. (19) with

[TABLE]

we can apply the Gronwall lemma to obtain that

[TABLE]

holds almost everywhere in $[0,T]$.

Using results from Thm. 7 it is now clear that

[TABLE]

∎

Remark 9**.**

One may think that in order to achieve $\mathbf{\ddot{u}}\in$ $L^{\infty}(0,T;L^{2}(\Omega))$ it is only required that $\mathbf{\dot{u}}(0)=\mathbf{u}_{1}\in H^{1}(\Omega)$ instead of $\beta\mathbf{u}_{1}\in H^{2}(\Omega)$ (or more precisely $\|\mathcal{B}^{T}\beta c^{E}\mathcal{B}\mathbf{u}_{1}\|_{L^{2}}(\Omega))<\infty$). However, this is not the case.

The condition $\beta\mathbf{u}_{1}\in H^{2}(\Omega)$ is required to show that $\|\mathbf{\ddot{u}}_{m}(0)\|_{L^{2}(\Omega)}$ is finite, such that the Gronwall inequality can be applied.

Remark 10**.**

In e.g. [4, p. 390, Eq. (59)] a $H^{2}$ regularity for $\mathbf{u}$ is achieved by selecting the test functions for $\mathbf{u}$ to be the complete eigenfunction sequence of $-\Delta\mathbf{u}$ which, indirectly, allows an estimation of $\|\mathbf{u}\|_{H^{2}(\Omega)}$. A similar argument should also be possible for $\mathcal{B}^{T}\mathcal{B}$ (or more precisely the operator that works on the solution vector $(\mathbf{u},\phi)^{T}$ and contains $\mathcal{B}^{T}\mathcal{B}$). This would directly increase the regularity of $\phi$ so that not only $\mathbf{u}\in L^{\infty}(0,T;H^{2}(\Omega))$ but also $\phi\in L^{\infty}(0,T;H^{2}(\Omega)\cap H^{1}_{0,\Gamma}(\Omega))$.

However, the authors did not follow that argumentation. Note that the here occurring differential operators are slightly different from the Laplacian. Hence, this leads to rather unpleasant changes due to the now very technical arguments and spaces. In that case, it would be possible to reduce the regularity requirements, however this would also change the resulting spaces and increase the cost of technical proof steps.

With the estimations and equations in Theorem 7 and the corresponding proof, a long-time behavior of the energy function $\eta$ and in particular of each component can be derived.

Corollary 11**.**

Let all requirements of Thm. 7 hold and let $\alpha,\beta>0$ strictly. If additionally there exists a $t_{0}\in\mathbb{R},\,t_{0}\geq 0$ such that $\phi^{e}(t)=0$ for $t\geq t_{0}$, then

[TABLE]

*for $t\rightarrow\infty$.

Furthermore the energy of the system*

[TABLE]

converges to a constant $\eta(t)\rightarrow c\in\mathbb{R}^{+}$ for $t\rightarrow\infty$.

Proof.

The right hand side $\mathcal{F}_{r}(t)$ of the energy balance Eq. (13) is constant for $t\geq t_{0}$ as no new energy is given into the system starting from time $t_{0}$, i.e. $\mathcal{F}_{r}(t)=c_{1}\in\mathbb{R}^{\geq 0}$ for $t\geq t_{0}$. Let

[TABLE]

and let

[TABLE]

Then Eq. (13) implies that $\tilde{\eta}(t)+\gamma(t)=c_{1}$ for $t\geq t_{0}$. As $\gamma(t)$ is monotonically increasing, it follows that $\tilde{\eta}(t)$ is monotonically decreasing. Both $\tilde{\eta}(t),\gamma(t)$ are bounded below and above by zero and $c_{1}$, respectively. Hence, $\gamma(t)$ and $\tilde{\eta}(t)$ must converge. Based on these results we get $0\leq\tilde{\eta}(t)\rightarrow c_{2}\leq c_{1}<\infty$ for $t\rightarrow\infty$. We know that $\gamma(t)$ converges, thus the occurring integrands must converge towards 0, i.e.,

[TABLE]

Through these results and by utilizing positive definiteness of $c^{E}$ and estimations similar to Eq. (14) we also get the following convergence result

[TABLE]

For the next steps it is already known that

[TABLE]

We have to show that one of the other two summands converges and determine the limit values.

We get the convergence of $\left<c^{E}\mathcal{B}\mathbf{u}_{m}(t),\mathcal{B}\mathbf{u}_{m}(t)\right>$ by taking advantage of Eq. (41) and using the characteristic of the time derivative of $\mathcal{B}\mathbf{u}_{m}(t)$. As all other summands converge this implies that also $\left<\epsilon^{S}\nabla\phi_{m}(t),\nabla\phi_{m}(t)\right>$ must converge. In order to specify the limit values we test the weak form Eq. (8) first with $(\mathbf{u}_{m}(t),0)$ and get

[TABLE]

to obtain

[TABLE]

and then a second time with $(0,\phi_{m}(t))$

[TABLE]

yielding

[TABLE]

From Eq. (43) and Eq. (44) we get

[TABLE]

Note that from Eq. (21) we are aware that $\|e\mathcal{B}\mathbf{u}_{m}(t)\|_{L^{2}(\Omega)}\rightarrow 0$ and using the requirement for $\phi^{e}(t)$ implies that also $\|\phi_{m}(t)\|_{H^{1}_{0}(\Omega)}\rightarrow 0.$

It is clear that $\tilde{\eta}(t)\rightarrow 0$ and with the characteristics of the material parameters

[TABLE]

we conclude

[TABLE]

for $t\rightarrow\infty$. Finally, we know that the derivatives in time and space of $\mathbf{u}_{m}(t)$ converge to zero, so $\|\mathbf{u}_{m}(t)\|_{L^{2}(\Omega)}\rightarrow\tilde{c}\in\mathbb{R}$.

Then we can conclude that $\eta(t)\rightarrow c\in\mathbb{R}^{+}$ for $t\rightarrow\infty$. ∎

As expected, we also find this theorized behavior in our numerical simulation results, see also Fig. 3. There, the monotonically decreasing energy term $\tilde{\eta}$ is shown. The electrode on the top of the piezoceramic disk is excited by the potential pulse as shown in Fig. 2. The time integration is given by a HHT-method, which is commonly used for piezoceramics (see [9]). These results were obtained by applying our simulation tool which will be focused on in upcoming publications. Note that small inaccuracies can occur due to numerical reasons.

Remark 12**.**

By using similar techniques as in the second part of the proof of Thm. 7, the last Corollary can be extended to non-discretized solutions of the partial differential equations.

4 Conclusion

Piezoelectric materials are widely diversified in their applications. Since measurements on real specimens are very expensive, computer simulations are used instead.

However, in order to confidently use these computer simulations, the underlying damped partial differential equation must be analyzed. In this paper, we prove existence, uniqueness and regularity of weak solutions of the governing partial differential equations and show some results on the long term behavior of solutions.

The obtained theoretical results are consistent with numerical results gained from a computational simulation of the model. With this, the basis is formed for ongoing design optimization of piezoelectric transducers.

Appendix A Definitions

Let $d,k\in\mathbb{N}$ be integers and let $\alpha$ be a multi-index. Then we define the functional spaces

[TABLE]

[TABLE]