# Development of accurate piezoelectric beam models based on Boley's method

J. Schoeftner, Ayech Benjeddou

► To cite this version:

J. Schoeftner, Ayech Benjeddou. Development of accurate piezoelectric beam models based on Boley's method. Composite Structures, 2019, 223, pp.110970. 10.1016/j.compstruct.2019.110970. hal-03991023

HAL Id: hal-03991023

https://hal.science/hal-03991023v1

Submitted on 26 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Development of accurate piezoelectric beam models based on Boley's method

<figure><img src="image_2.png" /></figure>

J. Schoeftner<sup>a,*</sup>, A. Benjeddou<sup>b,c</sup>

## ARTICLE INFO

### Keywords:

Boley's iterative method

Piezoelectricity

Beam modeling

Thick beams

Two-dimensional plane stress analytic solutions

## ABSTRACT

The solution of the bi-potential equation gives analytical formulas for two-dimensional isotropic problems. In 1954, Boley published a fast-converging, iterative method for rectangular cross-section elastic beams, in order to calculate in-plane stresses and deflections. In this contribution, Boley's iterative method is extended to piezoelectric materials. After defining potential functions for the in-plane stresses and electric fields, two coupled partial differential equations for these two potential functions are obtained. Extending the Boley's iterative method, by including the piezoelectric effect, it is shown that the solution of the first iteration exactly yields the analytical results for thin piezoelectric beams that are already known in the open literature. The first iteration also shows that the distribution of the electric potential through the beam thickness is a second-order polynomial with respect to the thickness coordinate. The second and third iterations are found to be higher-order terms, which depend on the thickness-to-length ratio of the beam. They can be interpreted as correction terms to obtain more accurate results for moderately thick beams. In this work, analytical and numerical results for a hinged-hinged bimorph under a sinusoidal force (sensing problem) or under a sinusoidal voltage load (actuation problem) are compared to analytical 2D plane stress results. It is shown that for higher values of the thickness-to-length ratio (i.e. moderately thick or very thick beams), where the Bernoulli-Euler beam bending theory yields unsatisfactory results, the latter are modified by the higher order Boley's method iterations such that the results converge to analytical 2D results.

## 1. Introduction

Smart and intelligent structures are systems that are equipped with multi-functional materials. Piezoelectric materials are prominent members of such systems. Thanks to the piezoelectric effects, it is possible to influence and to manipulate the motion or the stress state of a flexible structure if the electric potential or charge is prescribed. This is denoted as the converse piezoelectric effect; it is used for actuation applications. The direct piezoelectric effect results when a mechanical deformation or a stress changes the polarization, which affects the electric field and the electric displacement; it is exploited in sensor applications. Classical review works on this subject are, among others, Crawley [1], Miu [2] and Tzou [3].

Many practical problems can be modelled as truss structures (e.g. aeroplane fuselage, vehicle body, frame bridge), which basic elements are beams or truss elements; this is an important group in the field of active and passive vibration control. In this contribution, the focus is

made on piezoelectric strips or beam-like structures. One crucial point of piezoelectric laminates is the choice of an *a priori* accurate polynomial approximation order of the mechanical displacement and electric potential as functions of the thickness coordinate. A power series approach in the thickness direction based on the equivalent single layer (ESL) theory was proposed by Reddy [4] for purely elastic materials. Reddy [5], Mitchell and Reddy [6] have developed ESL theories for piezoelectric laminated composite plates. In general, the easiest modelling assumption for the electric field in the thickness direction is to consider it equal to the ratio of the electric potential difference to the thickness, i.e. the influence of the direct piezoelectric effect is ignored. This is also denoted as *small piezoelectric coupling*, see Tiersten [7]. For very thin beams, this modelling assumption can be sufficient, see Krommer [8] for example, but better results are achieved if the concept of the *induced potential*, as introduced by Benjeddou et al. [9], is considered. Tiersten [10] and Yang [11] modelled the electric field as a power series up to the third-order in the thickness direction, whereas

<footnote><sup>a</sup> Johannes Kepler University Linz, Institute of Technical Mechanics, Altenberger Strasse 69, A-4040 Linz, Austria</footnote>

<footnote><sup>b</sup> Institut Supérieur de Mécanique de Paris, 3 rue Fernand Hainaut, 93407 Saint Ouen CEDEX, France</footnote>

<footnote><sup>c</sup> ROBERVAL FRE 2012 UTC/CNRS, Centre de Recherche de Royallieu, CS 60319, 60203 Compiègne CEDEX, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: juergen.schoeftner@jku.at (J. Schoeftner), ayech.benjeddou@supmeca.fr (A. Benjeddou).</footnote>
Batra and Vidoli [12] used Legendre polynomials as basic functions. For more advanced electromechanically coupled beam theories, including higher-order shear deformation theories or the influence of the axial component of the electric field/displacement, the reader is referred to Krommer and Irschik [13] and Benjeddou et al. [9,14]. These modelling approaches are not restricted to systems, where the electric voltage is applied on the electrodes; they are also valid for systems, where an electric circuit is attached on the electrodes, see also Schoeftner and Irschik [15] and Schoeftner and Buchberger [16]. Furthermore, the electrodes do not necessarily have to fulfil the equipotential area conditions; they can be also considered as finitely conductive, see Schoeftner et al. [17,18].

In this contribution, we restrict our attention to a piezoelectric bimorph. The motivation of this work is influenced by an elegant method developed by Boley and Tolins [19], where an iterative method is presented to solve the bi-harmonic equation. This is a fourth-order linear partial differential equation (PDE) for two dimensional (2D) isotropic problems. From the outcome of the first iteration, one may derive the classical Bernoulli-Euler (BE) thin beam bending solution. The higher iterations can be interpreted as correction terms of the BE beam theory. They become more dominant for thick beams and thus increase the order of accuracy for beams where the thin assumption is violated. A further outcome of this method is that one may find accurate solutions for the stress distributions; e.g., the first iteration yields the well-known linear stress distribution, while the second iteration yields a third-order polynomial that gives a better approximation if the thin beam assumption is not valid anymore. An extension of this method for thermo-elastic problems is given by Boley [20] where explicit formulas for the stresses and the displacement are presented. Krommer and Irschik [21] extended this method to the charge equation of electrostatics in order to find appropriate solutions for the electric potential when the mechanical displacement is given in advance.

Motivated by the latter contribution of Krommer and Irschik [21], here it is investigated how to find coupled solutions for the mechanical and the electrical equilibria. Hereby, no a priori assumptions are made on the displacement field (such as cross-sections of the beam remain plane) and the electric field (such as negligence of the axial component). Therefore, in the following section, the Boley's iterative method for isotropic materials is first recalled. Then, in Section 3, the extension to piezoelectric materials of the previous method is presented. After introducing the 2D piezoelectric constitutive relations in strain-charge form and the mechanical potential function for the stress components and the electrical potential for the electric field, the constitutive relations are inserted into the compatibility equations and into the charge equation of electrostatics. Two coupled PDE are found for the mechanical and the electrical potential functions. For both functions, the iterative procedure is applied and analytical solutions are presented. Next, and before the conclusion section, these iteratively obtained solutions are presented in non-dimensional form and compared to analytical ones for mechanically (sinusoidal unit-width vertical force) and electrically (sinusoidal voltage) loaded hinged-hinged (doubly simply supported-SS) piezoelectric bimorph in order to verify the obtained results.

## 2. Recapitulation of Boley's method

For didactical reasons, the basic idea of Boley and Tolins [19], who calculated the in-plane stresses and deflections for isotropic materials by approximately solving the Airy stress function, is hereafter shortly introduced.

Elementary courses in mechanics dealing with homogenous and elastic thin beams state that the curvature $1/\rho$, where $\rho$ denotes the radius of the curvature, is related to the bending moment $M(x)$ and bending stiffness $K_M$ by this relation

$$ \frac{1}{\rho} = -\frac{M(x)}{K_M} \qquad (1) $$

This exactly holds if the bending moment is constant. Karman [22] and Seewald [23] derived the following more general equation, by means of the integral transform method, for a transversely loaded rectangular cross-section beam

$$ \frac{1}{\rho} = -\frac{M(x)}{K_M} + \frac{M_{,xx}(x)}{K_M} \left( \frac{4c^2}{5} + \frac{\nu c^2}{2} \right) \qquad (2) $$

showing the influence of the Poisson's ratio $\nu$, thickness $c$ and second derivative of the bending moment, $M_{,xx}(x)$. It is clear that the last term of Eq. (2) is negligible for smooth load distributions (i.e. qualitatively speaking for slowly varying loads with respect to $x$) or thin beams; i.e. for $c \ll L$. Boley and Tolins [19] also obtained this result of Eq. (2), when the applied load is perpendicular to the beam axis. Furthermore, they managed to investigate other loading conditions and derived solutions for the in-plane stress distributions.

The x-z in-plane stress ($\sigma$)-strain ($\epsilon$) relations for an isotropic material, having Young's modulus $E$ and Poisson's ratio $\nu$, read

$$ \begin{aligned} \epsilon_{xx} &= \frac{1}{E}\sigma_{xx} - \frac{\nu}{E}\sigma_{zz} \\ \epsilon_{zz} &= -\frac{\nu}{E}\sigma_{xx} + \frac{1}{E}\sigma_{zz} \\ \gamma_{xz} &= 2\epsilon_{xz} = \frac{2(1+\nu)}{E}\sigma_{xz} \end{aligned} \qquad (3) $$

It is noted that the equations in (3) reflect the plane stress assumptions for two-dimensional problems, when the stress components $\sigma_{xy} = \sigma_{yz} = \sigma_{yy} = 0$ vanish, see also Fig. 1. Inserting Eq. (3) and the definition of the Airy stress function $\psi$ (which fulfills the 2D body force free x-z in-plane equilibrium equations, see Timoshenko and Goodier [24])

$$ \sigma_{xx} = \psi_{,zzz}, \quad \sigma_{zz} = \psi_{,xx}, \quad \sigma_{xz} = -\psi_{,xz} \qquad (4) $$

into the compatibility equation

$$ \epsilon_{xx,zz} + \epsilon_{zz,xx} - 2\epsilon_{xz,xz} = 0 \qquad (5) $$

one obtains the bi-harmonic equation, a PDE for the Airy stress function

$$ \psi_{,zzzz} + 2\psi_{,xxzz} + \psi_{,xxxx} = 0 \qquad (6) $$

Following Boley's iterative method, one seeks for a solution $\psi$ in the form

$$ \psi = \psi_1 + \psi_2 + \psi_3 + \psi_4 + \dots \qquad (7) $$

where the functions $\psi_i$ satisfy the equations

$$ \psi_{i,zzzz} = -2\psi_{i-1,xxzz} - \psi_{i-2,xxxx}, \quad \psi_i = 0 \text{ for } i < 1 \qquad (8) $$

This means that one has to solve in a first step

$$ \psi_{1,zzzz} = 0 \qquad (9) $$

In a second step, one takes advantage of the first solution $\psi_1$ and solves for $\psi_2$

$$ \psi_{2,zzzz} = -2\psi_{1,xxzz} \qquad (10) $$

Then, the higher order solutions can be calculated in a straight forward iterative manner, see Eq. (8).

The boundary conditions (BCs) are fulfilled by the first ($i = 1$)

<figure><img src="image_3.png" /><figcaption>Fig. 1. Vertical unit-width force $q(x)$ (in N/m) for a two-layer beam (see Eqs. (11) and (12)).</figcaption></figure>
<table><caption>**Table 1**

Solutions for the first and second iterations of the mechanically loaded elastic beam (see Fig. 1).</caption><thead><tr><th>Iteration i</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Airy stress function ψ</td><td><sup>M(x)</sup>&frasl;<sub>I</sub> (<sup>z<sup>3</sup></sup>&frasl;<sub>6</sub> &minus; <sup>z<sup>2</sup></sup>&frasl;<sub>2</sub> + <sup>c<sup>3</sup></sup>&frasl;<sub>3</sub>)</td><td><sup>M<sub>xx</sub>(x)</sup>&frasl;<sub>I</sub> (<sup>z<sup>5</sup></sup>&frasl;<sub>60</sub> &minus; <sup>c<sup>2</sup>z<sup>3</sup></sup>&frasl;<sub>30</sub> + <sup>z<sup>4</sup></sup>&frasl;<sub>60</sub>)</td></tr><tr><td>Longitudinal deflection u(x, 0)</td><td>0</td><td><sup>νM<sub>x</sub>(x)c<sup>3</sup></sup>&frasl;<sub>3EI</sub></td></tr><tr><td>Vertical deflection w(x, 0)</td><td>&minus;<sup>1</sup>&frasl;<sub>EI</sub> ∮ M(x)dxdx</td><td><sup>M(x)</sup>&frasl;<sub>EI</sub> (<sup>4c<sup>2</sup></sup>&frasl;<sub>5</sub> + <sup>νc<sup>2</sup></sup>&frasl;<sub>2</sub>)</td></tr></tbody></table>

iteration only (other solutions fulfill trivial BCs). In case of a vertically loaded (at $z = -c$) beam (Fig. 1) with thickness $2c$, this means for iteration $i = 1$

$$ 
\begin{aligned}
\sigma_{zz}(x, -c) &= \psi_{1,xx}(x, -c) = -q(x)/b, & \sigma_{zz}(x, +c) &= \psi_{1,xx}(x, +c) = 0 \\
\sigma_{xz}(x, -c) &= -\psi_{1,xz}(x, -c) = 0, & \sigma_{xz}(x, +c) &= -\psi_{1,xz}(x, +c) = 0
\end{aligned}
\qquad (11) $$

While, for the higher iterations $i \ge 2$, these trivial BCs hold

$$ 
\begin{aligned}
\psi_{1,xx}(x, -c) &= 0, & \psi_{1,xx}(x, +c) &= 0 \\
\psi_{1,xz}(x, -c) &= 0, & \psi_{1,xz}(x, +c) &= 0 \quad i \ge 2
\end{aligned}
\qquad (12) $$

The results for the first and second iterations are given in Table 1, where the bending moment replaces the load distribution.

Note that the following relation holds between the second derivative of the moment and the unit-width vertical force (N/m) load (see Boley and Tolins [19])

$$ M_{xx}(x) = -q(x) \qquad (13) $$

It is noted that the signs in the solutions of Table 1 may differ from the results found by Boley and Tolins [19] since our coordinate system and the transverse load do not match with those chosen therein. One recognizes from Table 1 that the first iteration yields the well-known BE formula for the mid-plane deflection $w(x, 0)$ (second column-third row). The second iteration modifies this result (third column-third row), cf. Eq. (2). The second iteration also shows that the centerline (i.e. $z = 0$) is not the line of the neutral axis (third column-second row) unless $M_{xx}(x) = 0$ holds (e.g. pure bending) or in case of a vanishing Poisson's ratio.

Using the definition of the Airy stress function, see Eq. (4), one finds that only the fundamental solution, $i = 1$, yields a through-the-thickness linear distribution of the axial stress component but the second iteration yields a through-the-thickness third-order polynomial for $\sigma_{xx}$ as long as the second derivative of the bending moment is non-zero (third column-first row and consider $\sigma_{xx} = \psi_{zz}$).

### 3. Extension of Boley's method to piezoelectric materials

In this section, the Boley's procedure from Section 2 is extended by considering also the piezoelectric effects. In a first step, the constitutive equations are given and the stress and electric potential functions are defined. Then, two coupled PDE for these potentials are derived before presenting a similar iterative solution method to that given in Section 2.

#### 3.1. Electromechanical behavior and kinematics

The constitutive relations in strain-charge form relate the mechanical stresses $\sigma_{ij}$ and strains $\varepsilon_{ij}$ to the electric fields $E_i$ and displacements $D_i$, $i, j = 1, 2, 3$. In matrix notation, the corresponding constitutive equations (d-form) read

$$ 
\begin{aligned}
\varepsilon &= S^E \sigma + d^T E \\
D &= d \sigma + \varepsilon^\sigma E
\end{aligned}
\qquad (14) $$

here, the x-z in-plane vectors of the strain and the stress in Voigt's

notation and the electric displacement and electric field vectors read

$$ 
\begin{aligned}
\varepsilon &= \langle \varepsilon_{xx}, \varepsilon_{zz}, \gamma_{xz} \rangle^T & \sigma &= \langle \sigma_{xx}, \sigma_{zz}, \sigma_{xz} \rangle^T \\
D &= \langle D_x, D_z \rangle^T & E &= \langle E_x, E_z \rangle^T
\end{aligned}
\qquad (15) $$

For piezoceramics, the components of the short-circuit compliance matrix $S^E$ (at constant electric field), the strain piezoelectric coupling matrix $d$ and the electric permittivity matrix $\varepsilon^\sigma$ (at constant stress) are

$$ S^E = \begin{bmatrix} S_{11}^E & S_{13}^E & 0 \\ S_{13}^E & S_{33}^E & 0 \\ 0 & 0 & S_{55}^E \end{bmatrix}, \quad d = \begin{bmatrix} 0 & 0 & d_{15} \\ d_{31} & d_{33} & 0 \\ 0 & 0 & S_{55}^E \end{bmatrix}, \quad \varepsilon^\sigma = \begin{bmatrix} \varepsilon_{11}^\sigma & 0 \\ 0 & \varepsilon_{33}^\sigma \end{bmatrix} \qquad (16) $$

Inserting Eqs. (15) and (16) into Eq. (14), one finds the following six explicit relations

$$ 
\begin{aligned}
\varepsilon_{xx} &= S_{11}^E \sigma_{xx} + S_{13}^E \sigma_{zz} + d_{31} E_z \\
\varepsilon_{zz} &= S_{13}^E \sigma_{xx} + S_{33}^E \sigma_{zz} + d_{33} E_z \\
\gamma_{xz} &= S_{55}^E \sigma_{xz} + d_{15} E_x \\
D_x &= d_{15} \sigma_{xz} + \varepsilon_{11}^\sigma E_x \\
D_z &= d_{31} \sigma_{xx} + d_{33} \sigma_{zz} + \varepsilon_{33}^\sigma E_z
\end{aligned}
\qquad (17) $$

The strains in Eq. (17) are related to the displacement components $u(x, z)$ and $w(x, z)$ by

$$ \varepsilon_{xx} = u_{,x}, \quad \varepsilon_{zz} = w_{,z}, \quad \gamma_{xz} = u_{,z} + w_{,x} \qquad (18) $$

The gradient of the electric potential $\varphi$ represents the electric field vector. Hence, its components $E_x, E_z$ are

$$ E_x = -\varphi_{,x}, \quad E_z = -\varphi_{,z} \qquad (19) $$

#### 3.2. Compatibility equation and Gauss' law

Inserting the first three relations of Eq. (17) into the compatibility relation for the strains, Eq. (5), one finds

$$ (S_{11}^E \sigma_{xx,zz} + S_{13}^E \sigma_{zz,zz} + d_{31} E_{z,zz}) + (S_{13}^E \sigma_{xx,xx} + S_{33}^E \sigma_{zz,xx} + d_{33} E_{z,xx}) - (S_{55}^E \sigma_{xz,zz} + d_{15} E_{x,zz}) = 0 \qquad (20) $$

The displacement field in a dielectric medium satisfies

$$ D_{x,x} + D_{z,z} = \rho_f \qquad (21) $$

For piezoelectric materials, the free body charges do not exist; i.e., $\rho_f = 0$. Thus, substituting the last two relations of Eq. (17) into Eq. (21), one finds a second PDE relating the in-plane stresses to the electrical fields

$$ (d_{15} \sigma_{xz,x} + \varepsilon_{11}^\sigma E_{x,x}) + (d_{31} \sigma_{xx,z} + d_{33} \sigma_{zz,z} + \varepsilon_{33}^\sigma E_{z,z}) = 0 \qquad (22) $$

#### 3.3. Coupled PDE for mechanical and electric potential functions

Substituting the functions $\psi$ and $\varphi$ from Eqs. (4) and (19) into the coupled PDEs (20) and (22), one finds

$$ S_{11}^E \psi_{zzzz} + (2S_{13}^E + S_{55}^E) \psi_{xxxx} + S_{33}^E \psi_{xxxxx} - d_{31} \varphi_{zzz} - (d_{33} - d_{15}) \varphi_{zxx} = 0 \qquad (23) $$

$$ d_{31} \psi_{zzz} + (d_{33} - d_{15}) \psi_{xxz} - \varepsilon_{11}^\sigma \varphi_{xx} - \varepsilon_{33}^\sigma \varphi_{zz} = 0 \qquad (24) $$

These equations are the piezoelectric extensions of the bi-harmonic differential equation for the Airy stress function in the isotropic case, see Eq. (6). Seeking again for iterative solutions, the stress and electric potential functions, $\psi$ and $\varphi$, are written in this form

$$ 
\begin{aligned}
\psi &= \psi_1 + \psi_2 + \psi_3 + \psi_4 + \dots \\
\varphi &= \varphi_1 + \varphi_2 + \varphi_3 + \varphi_4 + \dots
\end{aligned}
\qquad (25) $$

Similar to Eq. (8), where the isotropic case is considered, iterative solutions have to satisfy the following relations if Eq. (25) is inserted into Eqs. (23) and (24)
$$ 
\begin{aligned}
& S_{11}^{E} \psi_{1,zzzz} + (2S_{13}^{E} + S_{55}^{E}) \psi_{1-1,xxzz} + S_{33}^{E} \psi_{1-2,xxxx} - d_{31} \varphi_{1,zzz} \\
& \quad - (d_{33} - d_{15}) \varphi_{1-1,xxz} = 0 \\
& \psi_{1} = 0 \text{ for } i < 1
\end{aligned} 
\quad (26) 
$$

$$ 
\begin{aligned}
& d_{31} \psi_{1,zzz} + (d_{33} - d_{15}) \psi_{1-1,xxz} - \epsilon_{11}^{\sigma} \varphi_{1-1,xx} - \epsilon_{33}^{\sigma} \varphi_{1,zz} = 0 \\
& \varphi_{1} = 0 \text{ for } i < 1
\end{aligned} 
\quad (27) 
$$

### 3.3.1. First iteration

The coupled differential equations (26) and (27) for the first iteration, $i = 1$, read

$$ S_{11}^{E} \psi_{1,zzzz} - d_{31} \varphi_{1,zzz} = 0 \quad (28) $$

$$ d_{31} \psi_{1,zzz} - \epsilon_{33}^{\sigma} \varphi_{1,zz} = 0 \quad (29) $$

These equations can be easily decoupled by inserting Eq. (29) into Eq. (28)

$$ S_{11}^{E} \left( 1 - \frac{d_{31}^{2}}{S_{11}^{E} \epsilon_{33}^{\sigma}} \right) \psi_{1,zzzz} = 0 \quad (30) $$

$$ \varphi_{1,zz} = \frac{d_{31}}{\epsilon_{33}^{\sigma}} \psi_{1,zzz} \quad (31) $$

It's worth mentioning that the modified compliance coefficient in Eq. (30) represents the compliance (at constant electric displacement) defined as

$$ S_{11}^{D} = S_{11}^{E} \left( 1 - \frac{d_{31}^{2}}{S_{11}^{E} \epsilon_{33}^{\sigma}} \right) = S_{11}^{E} (1 - k_{31}^{2}) \quad (32) $$

where, $k_{31}^2$ is the transverse-mode piezoelectric material electro-mechanical coupling factor (EMCF).

By introducing the product ansatz, one obtains that Eqs. (30) and (31) are fulfilled

$$ \psi_1(x, z) = f(x) \bar{\psi}_1(z), \quad \bar{\psi}_1(z) = a_{01} + a_{11}z + a_{21}z^2 + a_{31}z^3 \quad (33) $$

$$ \varphi_1(x, z) = f(x) \bar{\varphi}_1(z), \quad \bar{\varphi}_1(z) = g_{01} + g_{11}z + 3a_{31} \frac{d_{31} z^2}{\epsilon_{33}^{\sigma}} \quad (34) $$

where, the six coefficients $a_{01}$, $a_{11}$, $a_{21}$, $a_{31}$, $g_{01}$, $g_{11}$ and the function $f(x)$ depend on the BCs, see Section 4 and Eqs. (57)–(63).

### 3.3.2. Second iteration

Based on the outcome of the first iteration (see Eqs. (33) and (34)), the stress and electric potential functions $\psi_2$ and $\varphi_2$ for the second iteration, $i = 2$, can be calculated by solving these coupled PDEs

$$ S_{11}^{E} \psi_{2,zzzz} + (2S_{13}^{E} + S_{55}^{E}) \psi_{1,xxzz} - d_{31} \varphi_{2,zzz} - (d_{33} - d_{15}) \varphi_{1,xxz} = 0 \quad (35) $$

$$ d_{31} \psi_{2,zzz} + (d_{33} - d_{15}) \psi_{1,xxz} - \epsilon_{11}^{\sigma} \varphi_{1,xx} - \epsilon_{33}^{\sigma} \varphi_{2,zz} = 0 \quad (36) $$

Condensing Eq. (36) into Eq. (35) provides the following single PDE

$$ S_{11}^{E} \psi_{2,zzzz} + S_{s}^{E} \psi_{1,xxzz} - \bar{d}_{s} \varphi_{1,xxz} = 0 \quad (37) $$

with

$$ 
\begin{aligned}
\bar{S}_{11}^{E} &= S_{11}^{D} = S_{11}^{E} - d_{31}^{2} / \epsilon_{33}^{\sigma}, \\
\bar{S}_{s}^{E} &= 2S_{13}^{E} + S_{55}^{E} - d_{31} (d_{33} - d_{15}) / \epsilon_{33}^{\sigma}, \\
\bar{d}_{s} &= (d_{33} - d_{15}) - d_{31} \epsilon_{11}^{\sigma} / \epsilon_{33}^{\sigma}
\end{aligned} 
\quad (38) 
$$

Assuming the solutions of the second iteration as

$$ \psi_2(x, z) = f_2(x) \bar{\psi}_2(z) \quad (39) $$

$$ \varphi_2(x, z) = f_2(x) \bar{\varphi}_2(z) \quad (40) $$

and inserting Eqs. (33), (34), (39) and (40) into Eq. (37), the latter turns into

$$ \bar{S}_{11}^{E} f_2 \bar{\psi}_{2,zzzz} + \bar{S}_{s}^{E} f_{xx} \bar{\psi}_{1,zz} - \bar{d}_{s} f_{xx} \bar{\varphi}_{1,z} = 0 \quad (41) $$

The functions $\bar{\psi}_{1,zz}(z)$ and $\bar{\varphi}_{1,z}(z)$ in Eq. (41) are linear functions in $z$, one finds that $\bar{\psi}_2(z)$ is a 5th order polynomial in $z$, depending on the first iteration, see Eqs. (33) and (34). Consequently, it follows from Eq. (36), that $\bar{\varphi}_2(z)$ is a 4th order polynomial in $z$. Additionally one concludes (either from Eq. (35), (36) or (41)) that

$$ f_2(x) = f_{xx}(x) \quad (42) $$

must hold. The second iteration solutions have then the following form

$$ 
\begin{aligned}
\bar{\psi}_2(z) &= (5^{\text{th}} \text{ order polynomial in } z) + a_{02} + a_{12}z + a_{22}z^2 + a_{32}z^3 \\
\bar{\varphi}_2(z) &= (4^{\text{th}} \text{ order polynomial in } z) + g_{02} + g_{12}z
\end{aligned} 
\quad (43) 
$$

In a similar manner as for the first iteration, the six integration constants $a_{02}$, $a_{12}$, $a_{22}$, $a_{32}$, $g_{02}$, $g_{12}$ have to be determined, so that the trivial BCs are fulfilled, see Section 4.

### 3.3.3. Third iteration

The potentials $\psi_3$ and $\varphi_3$ for the third iteration are calculated by finding the solution of the following two equations

$$ 
\begin{aligned}
S_{11}^{E} \psi_{3,zzzz} + (2S_{13}^{E} + S_{55}^{E}) \psi_{2,xxzz} + S_{33}^{E} \psi_{1,xxxx} - d_{31} \varphi_{3,zzz} - (d_{33} - d_{15}) \varphi_{2,xxz} \\
= 0
\end{aligned} 
\quad (44) 
$$

$$ d_{31} \psi_{3,zzz} + (d_{33} - d_{15}) \psi_{2,xxz} - \epsilon_{11}^{\sigma} \varphi_{2,xx} - \epsilon_{33}^{\sigma} \varphi_{3,zz} = 0 \quad (45) $$

The solutions $\psi_3$ and $\varphi_3$ depend on the outcome of the first and second iterations, $(\psi_1, \varphi_1)$ and $(\psi_2, \varphi_2)$, which are already known. Assuming a solution of the form

$$ \psi_3(x, z) = f_3(x) \bar{\psi}_3(z) \quad (46) $$

$$ \varphi_3(x, z) = f_3(x) \bar{\varphi}_3(z) \quad (47) $$

one finds a relation for the x-dependent function as

$$ f_3(x) = f_{2,xx}(x) = f_{xxxx}(x) \quad (48) $$

and for the z-dependent functions of the potentials as

$$ 
\begin{aligned}
\bar{\psi}_3(z) &= (7^{\text{th}} \text{ order polynomial in } z) + a_{03} + a_{13}z + a_{23}z^2 + a_{33}z^3 \\
\bar{\varphi}_3(z) &= (6^{\text{th}} \text{ order polynomial in } z) + g_{03} + g_{13}z
\end{aligned} 
\quad (49) 
$$

The six integration constants, $a_{03}$, $a_{13}$, $a_{23}$, $a_{33}$, $g_{03}$, $g_{13}$, can be calculated by the (trivial) BCs, see Section 4.

### 3.3.4. Iterations $i \ge 4$

The differential equations for the higher iterations are already stated in Eqs. (26) and (27). In a similar manner as for the second and third iterations, one finds the x- and z-dependent functions via this recurrence

$$ f_i(x) = f_{i-1,xx}(x) = \frac{\partial^{2(i-1)} f(x)}{\partial x^{2(i-1)}} \quad (50) $$

and

$$ 
\begin{aligned}
\bar{\psi}_i(z) &= ((2i + 1)^{\text{th}} \text{ order polynomial in } z) + a_{0i} + a_{1i}z + a_{2i}z^2 + a_{3i}z^3 \\
\bar{\varphi}_i(z) &= (2i^{\text{th}} \text{ order polynomial in } z) + g_{0i} + g_{1i}z
\end{aligned} 
\quad (51) 
$$

### 3.3.5. Mechanical displacements

For the development of a higher order beam theory, the displacement components are of special interest. Inserting the first strain-displacement relations (18) and the stress and electric potential functions, Eqs. (4) and (19), into the first relations of Eqs. (17) and (18), one finds for the axial displacement expression

$$ u(x, z) = \int (S_{11}^{E} \psi_{zz} - d_{31} \varphi_{z}) dx + S_{13}^{E} \psi_{x} + U_{0}(z) \quad (52) $$
Furthermore, a relation for the second derivative with respect to x of the vertical deflection is found from the third relations of Eqs. (17) and (18)

$$w_{,xx} = -u + S_{55}^{E} \sigma_{xz,x} + d_{15} E_{x,x} \quad (53)$$

Inserting the stress and electric potential functions Eqs. (4) and (19) and the axial deflection Eq. (52) into Eq. (53), one derives this expression for the vertical deflection

$$w(x, z) = - \iint (S_{11}^{E} \psi_{1,zz} - d_{31} \varphi_{zz}) dx \, dx - (S_{13}^{E} + S_{55}^{E}) \psi_{z} - d_{15} \varphi + W_{0}(z) + W_{1}(z)x \quad (54)$$

Inserting the results for the potentials obtained in each iteration step (see Eqs. (25) and (50)), one finds a series of solutions in these forms, $u(x, z) = u_1(x, z) + u_2(x, z) + \dots$ and $w(x, z) = w_1(x, z) + w_2(x, z) + \dots$, or explicitly as

$$\begin{align*}
u(x, z) &= (S_{11}^{E} \bar{\psi}_{1,zz} - d_{31} \bar{\varphi}_{1,z}) \iint f(x) dx + U_{01}(z) + \dots \\
&= u_1(x,z) \\
(S_{11}^{E} \bar{\psi}_{2,zz} - d_{31} \bar{\varphi}_{2,z}) &+ S_{13}^{E} \bar{\psi}_{1,zz} f_{,x}(x) + U_{02}(z) + \dots \\
&= u_2(x,z) \\
(S_{11}^{E} \bar{\psi}_{3,zz} - d_{31} \bar{\varphi}_{3,z}) &+ S_{13}^{E} \bar{\psi}_{2,zz} f_{,xxx}(x) + U_{03}(z) + \dots \\
&= u_3(x,z) \tag{55}
\end{align*}$$

and

$$\begin{align}
w(x, z) &= -(S_{11}^{E} \bar{\psi}_{1,zz} - d_{31} \bar{\varphi}_{1,z}) \iint \overline{f(x)} dx \, dx + W_{01}(z) + W_{11}(z)x + \dots \nonumber \\
&\quad - [(S_{11}^{E} \bar{\psi}_{2,zz} - d_{31} \bar{\varphi}_{2,z}) + (S_{13}^{E} + S_{55}^{E}) \bar{\psi}_{1,z} + d_{15} \bar{\varphi}_{1}] f(x) \nonumber \\
&\quad + W_{02}(z) + W_{12}(z)x \nonumber \\
&\quad + \dots \nonumber \\
&\quad - [(S_{11}^{E} \bar{\psi}_{3,zz} - d_{31} \bar{\varphi}_{3,z}) + (S_{13}^{E} + S_{55}^{E}) \bar{\psi}_{2,z} + d_{15} \bar{\varphi}_{2}] \overline{f_{,xx}(x)} \nonumber \\
&\quad + W_{03}(z) + W_{13}(z)x \nonumber \\
&\quad + \dots \nonumber \\
&= w_1(x,z) \nonumber \\
&= w_2(x,z) \nonumber \\
&= w_3(x,z) \tag{56}
\end{align}$$

The unknown functions $U_{0i}(z)$, $W_{0i}(z)$, $W_{1i}(z)$, which are constants if the kinematic BCs are located at $z = 0$, have to fulfill the BCs for every iteration $i$, see Section 4 for the application to a hinged-hinged piezoelectric bimorph beam.

# 4. Application to a hinged-hinged piezoelectric bimorph beam

## 4.1. Problem statement

For the verification of the previously presented theory, a two-layer piezoelectric (bimorph) beam, as shown in Fig. 2, is studied. The two identical piezoceramic (PZT-5A) layers are assumed perfectly glued together (no slip of their interfaces) at $z = 0$. Each of them has length, thickness and width of $L$, $c$ and $b$.

It has to be stressed here that the thickness-to-length ratio, $\lambda = 2c/L$ is an essential parameter in the present investigations. The material properties of this benchmark are detailed in Appendix A. The electric potential at $z = 0$ is assumed nil, $\varphi_0 = 0V$, while the $x$-dependent voltage distribution at $z = -c$ is prescribed. Two load cases are considered:

<figure><img src="image_4.png" /><figcaption>Fig. 2. Hinged-hinged piezoelectric bimorph benchmark.</figcaption></figure>

* Mechanical loading (sensing problem): the distributed load $q(x) = q_0 \sin(\pi x/L)$ acts on the upper layer surface (at $z = -c$). The electric potentials at the external surfaces are short-circuited; i.e., $\varphi_1 = \varphi_2 = 0$ (so, the potential drops are: $V_1 = V_2 = 0$).

* Electrical loading (actuation problem): no mechanical load is considered, $q(x) = 0$, but the electric potentials at the external surfaces are now: $\varphi_1 = \varphi_2 = V_0 \sin(\pi x/L)$ (so, the resulting potential drops are: $V_1 = V_0 \sin(\pi x/L)$ and $V_2 = -V_0 \sin(\pi x/L)$).

The stress and electric potential functions, which have been derived for every iteration, contain 6 unknowns that are: $a_{0i}$, $a_{1i}$, $a_{2i}$, $a_{3i}$, $g_{0i}$, $g_{1i}$. Furthermore, 3 unknown functions, which have to be evaluated at $z = 0$, have to be determined; they are: $U_{0i}(0)$, $W_{0i}(0)$, $W_{1i}(0)$, see Eqs. (55) and (56). This means that for the two-layer problem, depicted in Fig. 2, 18 unknowns have to be calculated for each iteration. In this benchmark, it will be seen that there are:

* 6 support conditions for the mechanical displacements (see Section 4.1.1),

* 4 mechanical continuity conditions at the bimorph interface, $z = 0$ (see Section 4.1.2),

* 4 mechanical BCs for the stress (at $z = \pm c$, see Section 4.1.3),

* 4 electrical BCs for the prescribed voltages at $z = \pm c$ and $z = 0$ (see Section 4.1.4 for the electrical loading).

It is noted that the shear stress, the axial stress and the electric field at the vertical faces $x = 0$ and $L$ are related to the distributed load $q(x)$ or to the applied voltage $V(x)$, hence the relations (107) and (110) have to be fulfilled. For arbitrary stress or electric field distributions at the vertical faces, the reader is referred to Ruan et al. [25], where the residual stress or the residual electric field from (107) and (110) has to be calculated: this is related to the Saint-Venant end effect problem for piezoelectric materials.

### 4.1.1. Mechanical displacements BCs

For a doubly SS bimorph, the following kinematic restrictions hold for each iteration step $i$, see Fig. 2

$$\begin{align}
u_i^{(1)}(0, 0) &= u_i^{(2)}(0, 0) = 0 \nonumber \\
w_i^{(1)}(0, 0) &= w_i^{(2)}(0, 0) = 0 \nonumber \\
w_i^{(1)}(L, 0) &= w_i^{(2)}(L, 0) = 0; \quad i = 1, 2, 3\dots \tag{57}
\end{align}$$

where, the superscript (1) denotes the lower layer variables and superscript (2) that of the upper layer ones.

### 4.1.2. Mechanical continuity conditions at the bimorph interface

At the bimorph interface, the normal and shear stresses have to be equal for the lower and upper layers

$$\begin{align}
\sigma_{i,zz}^{(1)}(x, 0) &= \sigma_{i,zz}^{(2)}(x, 0) \nonumber \\
\sigma_{i,xz}^{(1)}(x, 0) &= \sigma_{i,xz}^{(2)}(x, 0); \quad i = 1, 2, 3\dots \tag{58}
\end{align}$$

Furthermore, perfect bonding is assumed between both layers (i.e. no relative motions)

$$\begin{align}
u_i^{(1)}(x, 0) &= u_i^{(2)}(x, 0) \nonumber \\
w_i^{(1)}(x, 0) &= w_i^{(2)}(x, 0); \quad i = 1, 2, 3\dots \tag{59}
\end{align}$$

### 4.1.3. Stress boundary conditions

At the external surface of the upper layer (at $z = -c$), the (negative) normal stress is $q(x)/b$, whereas the shear stress is zero (see Fig. 2)

$$\begin{align}
\sigma_{zz}^{(2)}(x, -c) &= -q(x)/b \rightarrow \psi_{1,xx}^{(2)}(x, -c) = -q(x)/b \nonumber \\
\psi_{1,xz}^{(2)}(x, -c) &= 0; \quad i \ge 2 \nonumber \\
\sigma_{xz}^{(2)}(x, -c) &= 0 \rightarrow -\psi_{i,xz}^{(2)}(x, -c) = 0; \quad i = 1, 2, 3\dots \tag{60}
\end{align}$$
The case of mechanical loading (this means Eq. (60) holds and instead of Eq. (62) the electrodes are short-circuited) requires further remarks: If the product ansatz Eq. (33) is inserted, the first equation in (60) yields $\bar{\psi}_1^{(2)}(-c) = 1/b$ and $f_{,xx}(x) = -q(x)$. Considering the discussion in Appendix C one can show that the unknown function $f(x)$ is equal to the bending moment distribution $M(x)$. For the second and third iterations Eq. (50) holds: $f_2(x) = f_{,xx}(x) = -q(x)$ and $f_3(x) = f_{,xx}(x) = -q_{,xx}(x)$.

It is noted that only the first iteration of the normal stress boundary condition in Eq. (60) fulfills the inhomogeneity $-q(x)/b$, while the higher iterations vanish. Hence, it is guaranteed that the sum of the iterative solutions $\psi_{,xx} = \psi_{1,xx} + \psi_{2,xx} + \dots$ fulfills this BC. This procedure seems plausible, because the first iteration mostly governs the solution: another kind of solution splitting or weighting may be advantageous for specific problems, but that is a matter of current research.

The stress components of the lower layer (at $z = c$) are traction-free

$$ 
\begin{aligned}
\sigma_{xx}^{(1)}(x, c) &= 0 \rightarrow -\psi_{1,xz}^{(1)}(x, c) = 0 \\
\sigma_{zz}^{(1)}(x, c) &= 0 \rightarrow \psi_{1,zz}^{(1)}(x, c) = 0; \quad i = 1, 2, 3\dots
\end{aligned}
\quad (61) 
$$

### 4.1.4. Electrical boundary conditions

At the external surface, non-trivial electrical BCs are allowed only for the first iteration, i.e.

$$ 
\begin{aligned}
\varphi_1^{(1)}(x, +c) &= V(x), \quad \varphi_1^{(1)}(x, +c) = 0; \quad i \ge 2 \\
\varphi_1^{(2)}(x, -c) &= V(x), \quad \varphi_1^{(2)}(x, -c) = 0; \quad i \ge 2
\end{aligned}
\quad (62) 
$$

Here, we also add a short remark: considering electrical loading only (this means all equations in (60) are zero) and inserting the product ansatz Eq. (33) into Eq. (62), one finds $\bar{\varphi}^{(1)}(c) = \bar{\varphi}^{(2)}(-c) = 1$ and $f(x) = V(x)$ for the first iteration. For the higher order iterations one obtains $f_3(x) = f_{,xxx}(x) = V_{,xxx}(x)$ and $f_3(x) = f_{,xxxx}(x) = V_{,xxxx}(x)$.

At the centerline, $z = 0$, the electric potential is kept at zero

$$ 
\begin{aligned}
\varphi_1^{(1)}(x, 0) &= 0 \\
\varphi_1^{(2)}(x, 0) &= 0; \quad i = 1, 2, 3\dots
\end{aligned}
\quad (63) 
$$

## 4.2. Analytical solutions

In order to give a condensed form of the solutions, the following non-dimensional parameters are introduced

$$ 
\begin{aligned}
k_{31}^2 &= d_{31}^2/(\epsilon_{33}^\sigma S_{11}^E), & \pi_S^E &= (2S_{13}^E + S_{35}^E)/S_{11}^E, & \nu_{13}^E &= -S_{13}^E/S_{11}^E \\
\pi_p &= (d_{33} - d_{15})/d_{31}, & \pi_d &= d_{15}/d_{31}, & \pi_\epsilon^\sigma &= \epsilon_{11}^\sigma/\epsilon_{33}^\sigma
\end{aligned}
\quad (64) 
$$

Additionally, the thickness influence is characterized by the ratio $\lambda = 2c/L$ and the non-dimensional coordinates are defined by $\bar{x} = x/L$ and $\bar{z} = z/c$. It will be seen that the electromechanical coupling factor $k_{31}^2$ and the thickness ratio $\lambda$ are the only non-dimensional parameters in the first iteration. The Poisson's ratio $\nu_{13}^E$, stiffness ratio $\pi_p^E$ (at constant electric field), piezoelectric ratios $(\pi_p, \pi_d)$ and permittivity ratio $\pi_\epsilon^\sigma$ (at constant stress) influence the solutions of the second iteration. Noticeably, more complex non-dimensional parameters can be defined for the third and higher iterations. Due to the lack of readability (even after introducing non-dimensional parameters), the results are not easy to interpret for higher iterations and, thus, the gain of scientific information is weak. Hence, hereafter, only the first and second iterations analytical results are presented. Nevertheless, numerical results (see Section 5) are shown up to order $i = 3$.

For the sake of completeness, the numerical values of the non-dimensional parameters are listed after inserting the PZT-5A material properties, see Table 3 (Appendix A), into Eq. (64)

$$ 
\begin{aligned}
k_{31}^2 &= 0.118, & \pi_S^E &= 2.016, & \nu_{13}^E &= 0.44 \\
\pi_p &= 1.228, & \pi_d &= -3.415, & \pi_\epsilon^\sigma &= 1.018
\end{aligned}
\quad (65) 
$$

$$ 
\begin{aligned}
\text{4.2.1. Mechanical loading: } q(\bar{x}) &= q_0 \sin(\pi \bar{x}), & \varphi_1(\bar{x}, \bar{z}) &= 1 \\
\varphi_2(\bar{x}, \bar{z}) &= -1
\end{aligned}
$$

The electric potential for the first iteration reads

$$ 
\varphi_1^{(1,2)}(\bar{x}, \bar{z}) = \frac{1}{\lambda b} \frac{L}{\pi^2} \frac{6 d_{31} (\pm 1 + k_{31}^2 \bar{z}) \bar{z}}{(4 - 3k_{31}^2) q_0 \sin(\pi \bar{x})} \quad (66) 
$$

The upper sign holds for the upper layer (i.e., superscript (2)) and the lower sign is for the lower layer (i.e., superscript (1)).

For the second iteration, the electric potential reads

$$ 
\varphi_2^{(1,2)}(\bar{x}, \bar{z}) = -\frac{\lambda L}{b} \frac{d_{31}}{\pi^3} \frac{(m_1 \bar{z} + m_2 \bar{z}^2 + m_3 \bar{z}^3 + m_4 \bar{z}^4)}{(4 - 3k_{31}^2)^2 (1 - k_{31}^2)} q_0 \sin(\pi \bar{x}) \quad (67) 
$$

with

$$ m_1 = \mp [(5k_{31}^2 - 4)\pi_s^E + (75k_{31}^4 - 177k_{31}^2 + 100)\pi_p + (21k_{31}^2 - 20)\pi_\epsilon^\sigma]/40 $$

$$ m_2 = -3[( - 5k_{31}^4 + 10k_{31}^2 - 4)\pi_s^E + (20k_{31}^4 - 42k_{31}^2 + 20)\pi_p + k_{31}^2\pi_\epsilon^\sigma]/20 $$

$$ m_3 = \mp (4 - 3k_{31}^2)[k_{31}^2(\pi_s^E - 2\pi_p) + \pi_\epsilon^\sigma]/4 $$

$$ m_4 = -(4 - 3k_{31}^2)[\pi_s^E - (k_{31}^2 + 1)\pi_p + \pi_\epsilon^\sigma]/8 $$

The vertical deflection yields for the first iteration

$$ w_1^{(1,2)}(\bar{x}, \bar{z}) = \frac{1}{\lambda^3 b} \frac{L}{\pi^4} S_{11}^E \frac{(1 - k_{31}^2)}{(4 - 3k_{31}^2)} q_0 \sin(\pi \bar{x}) \quad (68) $$

For the second iteration, one finds

$$ w_2^{(1,2)}(\bar{x}, \bar{z}) = \frac{1}{\lambda b} \frac{L}{\pi^2} S_{11}^E \frac{(n_0 + n_1 \bar{z} + n_2 \bar{z}^2)}{(4 - 3k_{31}^2)^2} q_0 \sin(\pi \bar{x}) \quad (69) $$

with

$$ n_0 = \frac{1}{5} [5\nu_{31}^E (3k_{31}^4 - 7k_{31}^2 + 4) + (10k_{31}^4 - 25k_{31}^2 + 16)\pi_s^E + (5k_{31}^4 - 7k_{31}^2)\pi_p + k_{31}^2 \pi_\epsilon^\sigma] $$

$$ n_1 = \mp k_{31}^2 (4 - 3k_{31}^2) (\nu_{31}^E + \pi_d + \pi_p) $$

$$ n_2 = -(4 - 3k_{31}^2)[\nu_{31}^E + k_{31}^2(\pi_d + \pi_p)] $$

The longitudinal deflection yields for the first iteration

$$ u_1^{(1,2)}(\bar{x}, \bar{z}) = -\frac{1}{\lambda^2 b} \frac{L}{\pi^3} S_{11}^E \frac{(1 - k_{31}^2) \bar{z}}{(4 - 3k_{31}^2)} q_0 \cos(\pi \bar{x}) \quad (70) $$

For the second iteration, one finds

$$ u_2^{(1,2)}(\bar{x}, \bar{z}) = \frac{L}{b} \frac{S_{11}^E}{\pi} q_0 (p_0 + p_1 \bar{z} + p_2 \bar{z}^2 + p_3 \bar{z}^3) \cos(\pi \bar{x}) \quad (71) $$

with

$$ p_0 = \frac{\nu_{31}^E (1 - \cos(\pi \bar{x}))}{2 \cos(\pi \bar{x})} $$

$$ p_1 = \frac{3}{5(4 - 3k_{31}^2)^2} [5\nu_{13}^E (3k_{31}^4 - 7k_{31}^2 + 4) + (5k_{31}^4 - 10k_{31}^2 + 4)\pi_s^E - (5k_{31}^4 - 7k_{31}^2)\pi_p - k_{31}^2\pi_\epsilon^\sigma] $$

$$ p_2 = \mp \frac{3}{2(4 - 3k_{31}^2)} k_{31}^2 (\nu_{13}^E + \pi_s^E - \pi_p) $$

$$ p_3 = -\frac{1}{(4 - 3k_{31}^2)} (\nu_{13}^E + \pi_s^E - k_{31}^2 \pi_p) $$

It's worth noticing that the electromechanical solutions under mechanical loading (sensing problem), Eqs. (66)–(71), depend on the length-to-thickness ($1/\lambda$) and the length-to-width ($L/b$) ratios. The latter ratio is a direct consequence of the applied load $q(\bar{x})/b$ in the stress BC, see Eq. (60), where $q(\bar{x})$ has the dimension N/m. We used this load since the first iteration results can be directly compared to beam theory results, where $q(\bar{x})$ represents the distributed beam load (in N/m). If we used $\sigma_{zz}(x, -c) = q^*(\bar{x})$ (where $q^*(\bar{x})$ has the dimension N/m$^2$), Eqs. (66)–(71) would not depend on the length-to-width ratio $(L/b)$, but on length $L$ only.
### 4.2.2. Electrical loading: $q(\bar{x}) = 0$, $\varphi_1(\bar{x}, \bar{z} = 1) = -\varphi_2(\bar{x}, \bar{z} = -1) = V_0 \sin(\pi\bar{x})$

The electric potential for the first iteration reads

$$ \varphi_1^{(1,2)}(\bar{x}, \bar{z}) = - \frac{(\pm 4 + 3k_{31}^2 \bar{z}) \bar{z}}{4 - 3k_{31}^2} V_0 \sin(\pi \bar{x}) \quad (72) $$

where, here also, the upper sign holds for the upper layer (i.e., superscript (2)) and the lower sign for the lower layer (i.e., superscript(1)).

For the second iteration one finds

$$ \varphi_2^{(1,2)}(\bar{x}, \bar{z}) = \lambda^2 \pi^2 \frac{(r_1 \bar{z} + r_2 \bar{z}^2 + r_3 \bar{z}^3 + r_4 \bar{z}^4)}{(4 - 3k_{31}^2)(1 - k_{31}^2)} V_0 \sin(\pi \bar{x}) \quad (73) $$

with

$$ r_1 = \pm \frac{1}{30(4 - 3k_{31}^2)} \left\{ -k_{31}^2 [\pi_s^E + (3k_{31}^2 - 5)\pi_p] + \pi_\sigma^E (9k_{31}^4 - 30k_{31}^2 + 20) \right\} $$

$$ r_2 = - \frac{k_{31}^2}{80(4 - 3k_{31}^2)} \left[ (36 - 25k_{31}^2) \pi_s^E + (-15k_{31}^4 + 93k_{31}^2 - 100) \pi_p + (-9k_{31}^2 + 20) \pi_\sigma^E \right] $$

$$ r_3 = \mp \frac{1}{6} [k_{31}^2 (\pi_s^E - 2\pi_p) + \pi_\sigma^E] $$

$$ r_4 = - \frac{k_{31}^2}{16} [\pi_s^E - (1 + k_{31}^2)\pi_p + \pi_\sigma^E] $$

The vertical deflection yields for the first iteration

$$ w_1^{(1,2)}(\bar{x}, \bar{z}) = - \frac{1}{\lambda^2} \frac{24}{\pi^2} \frac{(1 - k_{31}^2)}{(4 - 3k_{31}^2)} d_{31} V_0 \sin(\pi \bar{x}) \quad (74) $$

For the second iteration, one finds

$$ w_2^{(1,2)}(\bar{x}, \bar{z}) = \frac{(s_0 + s_1 \bar{z} + s_2 \bar{z}^2)}{(4 - 3k_{31}^2)} d_{31} V_0 \sin(\pi \bar{x}) \quad (75) $$

with

$$ s_0 = \frac{1}{10(4 - 3k_{31}^2)} [10\nu_{13}^E(4 - 3k_{31}^2) + (-5k_{31}^2 + 4)\pi_s^E + (15k_{31}^4 - 63k_{31}^2 + 60)\pi_p + (9k_{31}^2 - 20)\pi_\sigma^E] $$

$$ s_1 = \pm 4(\nu_{13}^E + \pi_d + \pi_p) $$

$$ s_2 = 3[\nu_{13}^E + k_{31}^2(\pi_d + \pi_p)] $$

The longitudinal deflection yields for the first iteration

$$ u_1^{(1,2)}(\bar{x}, \bar{z}) = \frac{1}{\lambda} \frac{12}{\pi} \frac{(1 - k_{31}^2) \bar{z}}{(4 - 3k_{31}^2)} d_{31} V_0 \cos(\pi \bar{x}) \quad (76) $$

For the second iteration, one finds

$$ u_2^{(1,2)}(\bar{x}, \bar{z}) = -\lambda \pi \frac{(t_1 \bar{z} + t_2 \bar{z}^2 + t_3 \bar{z}^3)}{(4 - 3k_{31}^2)^2} d_{31} V_0 \cos(\pi \bar{x}) \quad (77) $$

with

$$ t_1 = -\frac{1}{20} [10\nu_{13}^E (3k_{31}^2 - 4) + (25k_{31}^2 - 36)\pi_s^E + (15k_{31}^4 - 63k_{31}^2 + 60)\pi_p + (9k_{31}^2 - 20)\pi_\sigma^E] $$

$$ t_2 = \mp (3\pi_p - 4)(\nu_{13}^E + \pi_s^E - \pi_p) $$

$$ t_3 = -(3\pi_p - 4)(\nu_{13}^E + \pi_s^E - k_{31}^2 \pi_p)/2 $$

It's worth noticing here that, in contrary to the sensing problem, the electromechanical solutions of the present actuation problem, Eqs. (71)–(76), are independent of the aspect ratio (L/b).

# 4.3. Solutions discussion & comparison with existing theories

## 4.3.1. Vertical deflection

Within the framework of BE assumptions, one finds the following differential equation of the vertical deflection if the beam is mechanically loaded

$$ w_{\text{beam}}(\bar{x}) = \frac{L^4}{K_M} \iiint q(\bar{x}) d\bar{x} d\bar{x} d\bar{x} \quad (78) $$

while, if the beam is electrically loaded (voltage actuation), see Eq. (4) in Schoeftner et al [17]

$$ w_{\text{beam}}(\bar{x}) = \frac{e_{31} \text{ beam} c b L^2}{K_M} \iint V(\bar{x}) d\bar{x} d\bar{x} \quad (79) $$

where, $e_{31\text{beam}}$ is the transverse mode response stress piezoelectric coupling coefficient defined in Appendix A

Considering $q(\bar{x}) = q_0 \sin(\pi\bar{x})$ and $V(\bar{x}) = V_0 \sin(\pi\bar{x})$, one finds that Eqs. (68), (78) and Eqs. (74), (79) are identical for the hinged-hinged piezoelectric bimorph. Note that the bending stiffness of a piezoelectric bimorph is (dependent on thickness, $1/\lambda$, and aspect, L/b, ratios, as well as on axial compliance and EMCF)

$$ K_M(N. m^2) = c^3 b \left( \frac{2C_{11}^E}{3} + \frac{e_{31}^2}{6 \epsilon_{33}^E} \right) = \lambda^3 \left( \frac{b}{L} \right) \left( \frac{L^4}{S_{11}^E} \right) \frac{(4 - 3k_{31}^2)}{48(1 - k_{31}^2)} \quad (80) $$

where, $C_{11}^E$ and $\epsilon_{33}^E$ are the axial stiffness (at constant electric field) and the transverse dielectric (at constant strain) coefficients, see Eq. (93) in Appendix A.

## 4.3.2. Electric potential

Within the framework of BE assumptions, one finds the following relation between the potential difference of the piezoelectric layers and their vertical deformation, see Schoeftner et al [17,18] or Krommer [8]

$$ \varphi_{\text{beam}}(\bar{x}, \bar{z}) = \mp V(\bar{x})\bar{z} + \frac{\lambda^2 e_{31} \text{ beam}}{8 \epsilon_{33}^E} (\mp 1 - \bar{z})\bar{z} w_{\text{beam},\bar{x}\bar{z}}(\bar{x}) \quad (81) $$

where, the first term on the right-hand side denotes the through-the-thickness linear potential distribution, while the second through-the-thickness quadratic contribution represents the so-called induced potential [9]; it reflects the influence of the bimorph bending deformation (Section 4.3.1).

Inserting the beam deflections Eq. (68) (mechanical loading) or Eq. (74) (electrical loading) into Eq. (81), one finds that the outcome is equal to Eq. (66) (mechanical excitation) or Eq. (72) (electrical actuation). Hence the BE solutions for the electric potential as well as for the vertical deflection are identical to the Boley's ones, when only the first iteration is considered. In the following, they are used as synonyms.

## 4.3.3. Neutral fiber

The mechanical loading solution of the axial deflection, Eq. (71), shows that the neutral fiber is not at the geometric centroid (at $\bar{z} = 0$)

$$ u_2^{(1,2)}(\bar{x}, 0) = \frac{L}{b} \nu_{13}^E S_{11}^E q_0 \frac{1 - \cos(\pi \bar{x})}{2\pi} \quad (82) $$

After dividing Eq. (82) by the vertical displacement, Eq. (68), one finds that the order of magnitude is related to the cubic power of the thickness-to-length ratio $\lambda$

$$ O[u_2^{(1,2)}(\bar{x}, 0)]/O[w_1^{(1,2)}(\bar{x}, 0)] = \lambda^3 \nu_{13}^E \frac{\pi^3 (4 - 3k_{31}^2)}{96 (1 - k_{31}^2)} \tan(\pi\bar{x}/2) \quad (83) $$

Thus, the relative location of the neutral fiber is very close to the beam axis unless for very thick beams, i.e. $\lambda \gg 1/5$ (if the location $\bar{x} = 1$ is not accounted for).

For the electric loading, the location of the neutral fiber is unchanged, see Eq. (77), where $t_0 = 0$ holds.

## 4.3.4. Thickness deformation

The relative thickness deformation is $\Delta w_2^{(1)}(\bar{x}) = w_2^{(1)}(\bar{x}, 1) - w_2^{(1)}(\bar{x}, 0)$ for the lower layer and $\Delta w_2^{(2)}(\bar{x}) = w_2^{(2)}(\bar{x}, 0) - w_2^{(2)}(\bar{x}, -1)$ for the upper layer. It follows from Eq. (69) for the mechanical loading

$$ \Delta w_2^{(1,2)}(\bar{x}) = \pm \frac{1}{\lambda} \frac{L}{b} \frac{6}{\pi^2} \nu_{13}^E S_{11}^E \frac{(1 - k_{31}^2)}{(4 - 3k_{31}^2)} q_0 \sin(\pi \bar{x}) \quad (84) $$
One observes that this influence is of lower order than the first iteration of the vertical deformation (see Eq. (68)), since the ratio yields

$$O(\Delta w_2^{(1,2)})/O(w_1^{(1,2)}) = \pm \frac{\lambda^2 \nu_{13}^E \pi^2}{8} \qquad (85)$$

For the electrical loading, the relative thickness deformation is

$$\Delta w_2^{(1,2)}(\bar{x}) = \pm \frac{\nu_{13}^E + (4 - 3k_{31}^2)(\pi_d + \pi_p)}{(4 - 3k_{31}^2)} d_{31} V_0 \sin(\pi \bar{x}) \qquad (86)$$

Hence, it is also a lower order influence in comparison to the first iteration (see Eq. (74)), since the ratio yields

$$O(\Delta w_2^{(1,2)})/O(w_1^{(1,2)}) = \pm \lambda^2 \pi^2 \frac{\nu_{13}^E + (4 - 3k_{31}^2)(\pi_d + \pi_p)}{24(1 - k_{31}^2)} \qquad (87)$$

### 4.4. Numerical benchmarking results

In this section, numerical benchmarking results are presented for a SS piezoelectric bimorph. First (Section 4.4.1) considers a sinusoidal distributed load at the upper piezoelectric layer (see Fig. 2), $q(\bar{x}) = q_0 \sin(\pi \bar{x})$, where $q_0 = 1N/m$ holds. This is denoted as mechanical loading. Then (Section 4.4.2), a sinusoidal voltage distribution is applied to the outer electrodes, $V_1(\bar{x}) = -V_2(\bar{x}) = V_0 \sin(\pi \bar{x})$, where $V_0 = 1V$ holds (while the inner electrodes are kept at zero potential, see Fig. 2); this is the electric loading case. The outcome of the extended Boley's method (up to third iteration) is compared to analytical solutions (2D plane stress, see Appendix B). For the convergence study (Sections 4.4.1.3 and 4.4.2.3), when the relative error of the vertical deflection and the potential is computed, the thickness-to-length ratio $\lambda = 2c/L$ is varied from $\lambda = 1/2$ (very thick beam) to $\lambda = 1/80$ (very thin beam). The material properties of the benchmark are listed in Tables 3 and 4 of Appendix A. The length is $L = 1m$ and the layer thickness depends on the $\lambda$-ratio (note: $c = L\lambda/2$). It is noted that these geometric parameters cancel out for the depicted non-dimensional results, they are only relevant for the electric field in Figs. 6 and 7.

#### 4.4.1. Mechanical loading

For the mechanical loading, the load on the upper layer is assumed to be $q(\bar{x}) = q_0 \sin(\pi \bar{x})$, while the electrodes of the piezoelectric layer are short-circuited. Results are shown for the thickness-to-length ratio $\lambda = 1/20$ (thin beam) and $1/5$ (very thick beam). We show the non-dimensional displacement $\bar{w} = [(\lambda^3 b)/(q_0 S_{11}^E L)]w$ along the non-dimensional $\bar{x}$-axis (Figs. 3a, 4a), the non-dimensional electrical potential $\bar{\phi} = -[(\lambda \epsilon_{33}^g b)/(d_{31} L q_0)]\phi$ (see Figs. 3b, 4b) and the non-dimensional deflection of the cross section $\bar{u} = [(\lambda^2 b)/(q_0 S_{11}^E L)]u$ (see Figs. 3c, d, 4c and 4d) in the non-dimensional thickness direction $\bar{z}$. The proper choice of the non-dimensional coefficients for the electric potential and the deflections are extracted from the analytical results for the mechanical loading, see Eqs. (66), (68) and (70). Finally, a convergence study is performed for a wide range of thickness-to-length ratios $\lambda = [1/2 - 1/80]$, where the relative errors of the lateral centerline deflection and of the electric potential are plotted against the inverse of the thickness-to-length ratio (Fig. 5). For this study, the plane stress closed-form solution (CFS) is used as the reference solution, see Eq. (88).

#### 4.4.1.1. Thin piezoelectric bimorph: $\lambda = 1/20$. The results for the thin beam ($\lambda = 1/20$) are shown in Fig. 3. The centerline deflection in the middle of the beam (Fig. 3a) is $\bar{w}(0.5, 0) = 0.1192$ for the BE-beam theory (blue, iteration $i = 1$). The second iteration (red) improves the results and matches quite well with the 2D plane stress CFS (black) $\bar{w}(0.5, 0) = 0.120$. No deviations can be recognized from the second- and third-order (grey) Boley's solutions. The voltage distribution in the thickness direction (Fig. 3b) shows the typical parabolic distribution in case of short-circuited electrodes $\bar{\phi}(0.5, -0.5) = 0.0417$, see Eq. (66). Additionally, the cross-section deformation in the $\bar{x}$-direction $\bar{u}(0.4, \bar{z})$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Distributions of non-dimensional vertical deflection $\bar{w}(\bar{x}, 0)$ (a), voltage $\bar{\phi}(0.5, \bar{z})$ (b), axial deflection $\bar{u}(0.4, \bar{z})$ (c) and deviation of the axial deflection from the linear distribution $\Delta\bar{u}(0.4, \bar{z})$ (d) for a thin beam $\lambda = 1/20$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Distributions of non-dimensional vertical deflection $\bar{w}(\bar{x}, 0)$ (a), voltage $\bar{\phi}(0.5, \bar{z})$ (b), axial deflection $\bar{u}(0.4, \bar{z})$ (c) and deviation of the axial deflection from the linear distribution $\Delta\bar{u}(0.4, \bar{z})$ (d) for a very thick beam $\lambda = 1/5$.</figcaption></figure>

is shown in Fig. 3c, showing the predominantly linear distribution for thin beams, i.e. $u(x, z) \approx -zw_x(x, z = 0)$ or $\bar{u}(\bar{x}, \bar{z}) \approx -\bar{z}\bar{w}_x(\bar{x}, \bar{z} = 0)\lambda/2$ in non-dimensional form: the “cross-section remains plane” BE assumption outperforms the higher order effects. In order to show the latter, one can subtract the linear distribution, see Fig. 3d: one recognizes also that the neutral fiber does not cross the beam axis: the axial deformation is $\bar{u}(\bar{x}, \bar{z} = 0) \approx 1.2 \times 10^{-4}$.

#### 4.4.1.2. Very thick piezoelectric bimorph: $\lambda = 1/5$. The very thick beam
<figure><img src="image_7.png" /><figcaption>Fig. 5. Convergence study – relative errors of non-dimensional vertical deflection $\bar{w}(0.5, 0)$ (a, b) and voltage $\bar{\phi}(0.5, -0.5)$ (c, d) as a function of the biomorph inverse of the thickness-to-length ratio (mechanical loading).</figcaption></figure>

results ($\lambda = 1/5$) are shown in Fig. 4. The vertical deflection (Fig. 4a) is $\bar{w}(0.5, 0) = 0.1192$ for the BE-beam (blue, iteration $i = 1$), which underestimates the exact result 0.1311 (black). The second Boley iteration yields 0.1312 (red), whereas the third iteration yields 0.1311(grey) as the CFS result. The voltage distributions (Fig. 4b) are symmetrical with respect to the $\bar{z}$-axis, i.e. they are even functions, see Eqs. (66) and (67), for the first and second Boley iterations (e.g. BE-beam (blue) with $\bar{\phi}(0.5, \pm 0.5) = 0.0417$ or $\bar{\phi}(0.5, \pm 0.5) = 0.04338$ (red, second Boley iteration). The latter result is in good agreement with the analytical solution. Only the third Boley iteration (grey) shows that the potential distributions for upper and lower layers are not perfectly symmetric, and it matches best with the analytical result of course (e.g. $\bar{\phi}(0.5, -0.5) = 0.04348$ and $\bar{\phi}(0.5, 0.5) = 0.04332$ (grey, third iteration) and $\bar{\phi}(0.5, -0.5) = 0.04345$ and $\bar{\phi}(0.5, 0.5) = 0.04329$ (black, analytical). Fig. 4c and Fig. 4d indicate that the cross-section does not remain plane anymore and the neutral fiber does not cross the centerline. Indeed, the second Boley term yields $\bar{u}(\bar{x}, \bar{z} = 0) = 0.00194$, which yields the absolute midline displacement $u(x, z = 0) = 7.94 \times 10^{-13}m$. The latter dimensional value does not depend on the thickness ratio $\lambda$, see Eq. (71), where the thickness-to-length ratio cancels out. The Boley iteration up to order 3 yields $\bar{u}(\bar{x}, \bar{z} = 0) = 0.00202$, which is very close to the analytical result $\bar{u}(\bar{x}, \bar{z} = 0) = 0.002017$.

It is noted that only thanks to the proper choices of the non-dimensional coefficients (see Section 4.4.1) that BE solutions are independent of the thickness ratio $\lambda$; so, the second order corrections are clearly visible, cf. Figs. 3 and 4.

### 4.4.1.3. Convergence study.
A convergence study is performed to measure the convergence behavior of the presently extended Boley's method (EBM). The percent (%) relative (Rel.) error of an arbitrary indicator function $e$, obtained by the EBM, is calculated as

$$ \text{Rel. error (\%)} = 100 \frac{e_{\text{EBM}} - e_{\text{ref.}}}{e_{\text{ref.}}} \qquad (88) $$

The function $e$ represents here either the deflection in the middle $\bar{w}(0.5, 0)$ (Fig. 5a, and the zoom Fig. 5b) or the electric potential in the middle of the lower layer $\bar{\phi}(0.5, -0.5)$ (Fig. 5c, and the zoom Fig. 5d). The reference (ref.) is the analytical 2D plane stress CFS presented in Appendix B. The relative errors for the Boley iterations up to order three are plotted as functions of the inverse of the thickness ratio $\lambda$.

As expected, one observes that for very thin beams, the error converges to zero. For both the deflection and the electric potential, the second Boley iteration improves the results: e.g., for a moderately thin beam $\lambda = 1/10$, the BE result yields an error of about 2.5%, two Boley iterations cause an error of 0.0047% and three Boley iterations cause an

error of only 0.0001%. From a practical point of view, thickness ratios $\lambda > 1/10$ are not considered in general as beam structures, and it is well-known that BE formulas (i.e. first Boley iteration) can only be used as a very rough approximation. For an extremely thick piezoelectric two-layer beam, $\lambda = 1/2$, the errors are 37.5%, (blue, BE), 1.9% (red, Boley up to order $i = 2$) and 0.02% (grey, Boley up to order $i = 3$). For $\lambda = 1/5$, the error is only 0.07% and 0.0006% for the second and third order iterations.

For the electrical potential, the second iteration also yields a huge improvement of the result: e.g., for $\lambda = 1/8$ (moderately thick beam), the error is 0.02% instead of 1.6% (BE). The third iteration slightly yields better results as the error is only 0.01%.

In general, one can conclude from Fig. 5 that the presently extended Boley's method converges for thickness ratios $\lambda > 1$; this follows from the fact that the order of magnitude increases by a factor $\lambda^2$ for every higher iteration, see e.g. Eqs. (85) and (87).

### 4.4.2. Electrical loading

For the electrical loading, the external mechanical load vanishes, $q(\bar{x}) = 0$, while the electric voltage is $V_1(\bar{x}) = -V_2(\bar{x}) = \sin(\pi\bar{x})$. Results are plotted for the thickness-to-length ratio $\lambda = 1/20$ (thin beam) and $\lambda = 1/5$ (very thick beam). We show the non-dimensional lateral deflection $\bar{w} = -[(\lambda^2/(d_{31}V_0))w$ along the $\bar{x}$-axis (Figs. 6a, 7a), the non-dimensional potential distribution $\bar{\phi} = \phi/V_0$ in the non-dimensional thickness direction (see Figs. 6b, 7b) and the electric field in the thickness (Figs. 6c, 7c) and axial (Figs. 6d, 7d) directions. The amplitude of the sinusoidal voltage distribution is denoted by $V_0$, see Section 4.1. Notice that the electric field is intentionally not made non-dimensional as the other results. This choice is made in order to assess the validity of the common piezoelectric beam modeling assumption of the negligence of $E_x$ in comparison to $E_z$. Besides, it's worth mentioning that the proper choice of the non-dimensional coefficients for the electric potential and deflections are extracted from the analytical results for the electrical loading, see Eqs. (72), (74) and (76). Finally, a convergence study is performed for various thickness-to-length ratios, $\lambda = [1/2 - 1/80]$. For this purpose, the relative errors, with regards to analytical 2D plane stress CFS (see Appendix B), of the lateral centerline deflection and of the electric potential are plotted against the inverse of

<figure><img src="image_8.png" /><figcaption>Fig. 6. Distributions of non-dimensional vertical deflection $\bar{w}(\bar{x}, 0)$ (a), voltage $\bar{\phi}(0.5, \bar{z})$ (b) and real thickness and axial electric fields $E_z(L/4, z)$ (c) and $E_x(L/4, z)$ (d) for a thin beam $\lambda = 1/20$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Distributions of non-dimensional vertical deflection $\bar{w}(\bar{x}, 0)$ (a), voltage $\bar{\phi}(0.5, \bar{z})$ (b) and real thickness and axial electric fields $E_z(L/4, z)$ (c) and $E_x(L/4, z)$ (d) for a very thick beam $\lambda = 1/5$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Convergence study – relative errors of the non-dimensional vertical deflection $\bar{w}(0.5, 0)$ (a,b) and voltage $\bar{\phi}(0.5, -0.5)$ (c,d) as a function of the bimorph inverse of the thickness-to-length ratio (electrical loading).</figcaption></figure>

the thickness-to-length ratio, $1/\lambda$ (Fig. 8).

#### 4.4.2.1. Thin piezoelectric bimorph: $\lambda = 1/20$.
The results for the thin beam ($\lambda = 1/20$) are shown in Fig. 6. The vertical bimorph deflection (Fig. 6a) is $\bar{w}(0.5, 0) = 0.5882$ for the BE-beam theory (blue, iteration $i = 1$). The Boley iterations up to order three yield $\bar{w}(0.5, 0) = 0.5869$ (red, grey), which is also the result found by the analytical 2D plane stress CFS calculations (black). The potential distribution in the non-dimensional thickness direction (Fig. 6b) is shown to be dominated by a linear part, see Eq. (72) and the discussion in Section 4.3.2 and Eq. (81). The higher order influence, which is dominated by the quadratic part, is hardly recognizable and stems mainly due to the thickness-dependent deformation field, in particular due to $\epsilon_{xx} \neq 0$ and $w_{,xx} \neq 0$, see Eq. (81). The electric potential in the middle of the layer is $\bar{\phi}(0.5, -0.5) = 0.5244$ for the BE-theory (blue, Boley iteration $i = 1$) and 0.5239 for the three other cases under consideration (red, grey, black). Comparing the electric field components (Fig. 6c and d), one observes that the

absolute values of $E_z$ are in general larger than those of $E_x$. The electric field $E_z$ varies almost linearly in the thickness direction and depends on the piezoelectric EMCF, $k_{31}^2$, see Eq. (72).

#### 4.4.2.2. Very thick piezoelectric bimorph: $\lambda = 1/5$.
The vertical deflection (Fig. 7a) is $\bar{w}(0.5, 0) = 0.5882$ for the BE-beam (blue, iteration $i = 1$) again (cf. Fig. 6a, blue). It overestimates the analytical 2D plane stress CFS result $0.5679$ (note that for the mechanically loaded case, the analytical 2D plane stress CFS deflection was underestimated by BE, see Fig. 4a). The higher order Boley iterations are in very good agreement with the analytical 2D plane stress CFS results. At $\bar{x} = 0.5$ and $\bar{z} = -0.5$, the electric potential is $0.5244$ (blue), the same as for the thin beam (see Fig. 6b in Section 4.4.2.1). It can be shown that the first iteration does not consider the thickness ratio, see Eq. (72). But, the second iteration (see Eq. (73)) depends on $\lambda^2$, lowering the electric voltage to $0.5173$: this is very close to Boley's third iteration result and the analytical 2D plane stress CFS value, which are both equal to $0.5174$. Fig. 7c shows that the electric field in thickness direction increases linearly from $6.38V/m$ to $7.76V/m$ for the BE result (blue), see Eq. (72) by considering $E_z = -\phi_z$: it is a better result than $E_z \approx -7.071V/m = -V_1/c = -\sin(0.25\pi)/c$ (lower layer), when the electric field is simply approximated by the difference quotient; the latter is an often seen approximation in the open literature when modeling piezoelectric beams: it is noted that this simple approach yields reliable results only for very thin beams, see also Eq. (81), when this assumption is justified. Then both the $x$- and also the $z$-derivatives of the electric field can be neglected in Gauss' law of electrostatics Eq. (22). The second- and third-order Boley iterations match perfectly with the analytical results: both approximations show the higher order dependency on the thickness direction, which is only relevant for the thicker beam, see Eq. (73), but not for the thin beam (compare Figs. 7c to 6c).

#### 4.4.2.3. Convergence study.
Finally, a convergence study is also presented for the non-dimensional electrical loading (for the vertical deflection $\bar{w}(0.5, 0)$ (Fig. 8a and b) and the electric potential $\bar{\phi}(0.5, -0.5)$ in the upper layer (Fig. 8c and d) again. The relative error, with regards to the analytical plane stress CFS as in Appendix B, is calculated via Eq. (88).

For thin beams (e.g. $\lambda < 1/10$), the BE approach (Boley first iteration) yields sufficient results, where the relative error is lower than 1% for both the deflection as well as the electric potential. For a very thick beam ($\lambda = 1/5$, see also Fig. 7), it is found that the relative error of the displacement is reduced from 3.6% (blue, BE) to 0.07% (red, Boley second iteration) and 0.007% (grey, Boley third iteration), see Fig. 8a and b. For the electric potential, the error remains practically zero (rel. error <0.1%) for inverse thickness-to-length ratio higher than 5 (red, Boley second iteration) or higher than 3 (grey, Boley third iteration), see Fig. 8c and d.

## 5. Conclusions

In the present work, an iterative solution procedure has been presented for the explicit analytical and numerical evaluation of electro-mechanical states of piezoelectric beams under static mechanical and electric loadings. For this purpose, the Boley's method, presented initially for deriving analytically moderately thick and thick beam higher order solutions, was here extended by taking into account the piezoelectric inverse and direct effects without any *a priori* (axiomatic) assumptions. Two coupled PDE were derived, from the compatibility condition and Gauss' law of electrostatics, and solved iteratively. Analytical formulas and numerical results were given for a hinged-hinged (SS) piezoelectric bimorph, which is loaded either mechanically (sensing problem) or electrically (actuation problem). It was found that the first iterative solution perfectly represents the piezoelectric BE beam solutions, which have already been derived in the literature; as
expected, the higher order solutions were found to depend on the bi-morph thickness-to-length ratio. Results obtained from the first, second and third Boley iterations were compared to those from analytical plane stress CFS given in Appendix B. It was found that the second iteration corrects the BE beam solution (i.e. the outcome of the first iteration) in such a manner that it considers the higher order thickness variations of the deformation, the stress fields and the electric potentials.

In summary, the comparison with analytical results showed that the first iteration corresponds to the BE beam solution for piezoelectric material (BE<sub>piezo</sub>) where the induced potential concept is considered. For thin beams, the higher order solutions are negligible. But for beams with a thickness-to-length ratio λ larger than 1/20, the second iteration significantly increases the accuracy. Convergence studies for very thick beams (λ > 1/3) show that the iterations up to the third order result into relative errors lower than 0.5% for the deflection. The main ad-

vantage of the presented analytical expressions is that beam theories for piezoelectric material can be systematically derived for various kinds of loading and boundary conditions without assuming restrictions for the deflections, the electric fields or the electric displacements.

## Acknowledgements

This contribution is an extension of the work presented at ICCS19 (19th International Conference on Composite Structures) in Porto on 8th September 2016. The first author acknowledges support from the Austrian Science Fund FWF via the project P 26762-N30 and from the Johannes Kepler University Linz (Linz Institute of Technology–LIT-Project LIT-2017-4-YOU-004 RBMBM) and the second author acknowledges support from Linz Center of Mechatronics via COMET K2 Symbiotic Mechatronics project n° 9204849.

# Appendix A

## PZT-5A material data

In matrix notations, the stress-charge form of the 3D piezoelectric constitutive equations read

$$ 
\begin{aligned}
\sigma &= C^E \varepsilon - e^T E \\
D &= e \varepsilon + \varepsilon^T E
\end{aligned}
 $$

(89)

In strain-charge form, the former read

$$ 
\begin{aligned}
\varepsilon &= S^E \sigma + d^T E \\
D &= d \sigma + \varepsilon^T E
\end{aligned}
 $$

(90)

The meanings of the variables used in Eqs. (89) and (90) are given in Table 2.

<table><caption>**Table 2**
Nomenclature.</caption><thead><tr><th>Quantity</th><th>Meaning</th></tr></thead><tbody><tr><td>σ</td><td>Vector of stress components</td></tr><tr><td>ε</td><td>Vector of strain components</td></tr><tr><td>D</td><td>Electric displacement vector</td></tr><tr><td>E</td><td>Electric field vector</td></tr><tr><td>C<sup>E</sup></td><td>Stiffness matrix at constant electric field (short-circuit)</td></tr><tr><td>S<sup>E</sup></td><td>Compliance matrix at constant electric field (short-circuit)</td></tr><tr><td>e</td><td>Stress piezoelectric coupling matrix</td></tr><tr><td>d</td><td>Strain piezoelectric coupling matrix</td></tr><tr><td>ε<sup>ε</sup></td><td>Dielectric (permittivity) matrix at constant strain</td></tr><tr><td>ε<sup>σ</sup></td><td>Dielectric (permittivity) matrix at constant stress</td></tr></tbody></table>

One can switch between these two forms by applying the following transformations rules

$$ S^E = (C^E)^{-1}, \quad d = e(C^E)^{-1}, \quad \varepsilon = \varepsilon^\varepsilon + e(C^E)^{-1}e^T $$

(91)

or

$$ C^E = (S^E)^{-1}, \quad e = d(S^E)^{-1}, \quad \varepsilon = \varepsilon^\sigma - d(S^E)^{-1}d^T $$

(92)

Inserting the values from Table 3, showing PZT-5A data in strain-charge form, into Eq. (92), one obtains the corresponding stress-charge values listed in Table 3.

In order to find the relevant parameters on beam-level (when the normal and shear stress components and the electric field in x-direction are neglected, see Schoeftner and Irschik [15]), which are used for comparing the first Boley iteration results to the BE results from the literature (see the discussion in Sections 4.3.1 and 4.3.2 for the vertical deflection and for the potential distribution), one finds the following 1D transverse mode response relations and values

$$ C_{11}^E = 1/S_{11}^E = 61.0 \cdot 10^9 \text{Nm}^{-2} $$

$$ e_{31}^E = d_{31}/S_{11}^E = 10.43 \text{Asm}^{-2} (\text{Cm}^{-2}) $$

$$ \varepsilon_{33}^E = \varepsilon_{33}^\sigma - d_{31}^2/S_{11}^E = 13.31 \cdot 10^{-9} \text{AsV}^{-1}\text{m}^{-1} (\text{Fm}^{-1}) $$

(93)
<table><caption>**Table 3**

Strain-charge form of PZT-5A 3D parameters.</caption><thead><tr><th>Parameter</th><th>value</th><th>unit</th></tr></thead><tbody><tr><td>S<sup>E</sup><sub>11</sub></td><td>16.4 × 10<sup>-12</sup></td><td>N<sup>-1</sup>m<sup>2</sup></td></tr><tr><td>S<sup>E</sup><sub>12</sub></td><td>-5.74 × 10<sup>-12</sup></td><td>N<sup>-1</sup>m<sup>2</sup></td></tr><tr><td>S<sup>E</sup><sub>13</sub></td><td>-7.22 × 10<sup>-12</sup></td><td>N<sup>-1</sup>m<sup>2</sup></td></tr><tr><td>S<sup>E</sup><sub>33</sub></td><td>18.8 × 10<sup>-12</sup></td><td>N<sup>-1</sup>m<sup>2</sup></td></tr><tr><td>S<sup>E</sup><sub>55</sub></td><td>47.5 × 10<sup>-12</sup></td><td>N<sup>-1</sup>m<sup>2</sup></td></tr><tr><td>d<sub>31</sub></td><td>-171 × 10<sup>-12</sup></td><td>mV<sup>-1</sup>(CN<sup>-1</sup>)</td></tr><tr><td>d<sub>33</sub></td><td>374 × 10<sup>-12</sup></td><td>mV<sup>-1</sup>(CN<sup>-1</sup>)</td></tr><tr><td>d<sub>15</sub></td><td>584 × 10<sup>-12</sup></td><td>mV<sup>-1</sup>(CN<sup>-1</sup>)</td></tr><tr><td>ε<sup>σ</sup><sub>11</sub></td><td>15.32 × 10<sup>-9</sup></td><td>AsV<sup>-1</sup>m<sup>-1</sup>(Fm<sup>-1</sup>)</td></tr><tr><td>ε<sup>σ</sup><sub>33</sub></td><td>15.05 × 10<sup>-9</sup></td><td>AsV<sup>-1</sup>m<sup>-1</sup>(Fm<sup>-1</sup>)</td></tr></tbody></table>

<table><caption>**Table 4**

Stress-charge form of PZT-5A 3D parameters.</caption><thead><tr><th>Parameter</th><th>value</th><th>unit</th></tr></thead><tbody><tr><td>C<sup>E</sup><sub>11</sub></td><td>120.35 × 10<sup>9</sup></td><td>Nm<sup>-2</sup></td></tr><tr><td>C<sup>E</sup><sub>12</sub></td><td>75.18 × 10<sup>9</sup></td><td>Nm<sup>-2</sup></td></tr><tr><td>C<sup>E</sup><sub>13</sub></td><td>75.09 × 10<sup>9</sup></td><td>Nm<sup>-2</sup></td></tr><tr><td>C<sup>E</sup><sub>33</sub></td><td>110.87 × 10<sup>9</sup></td><td>Nm<sup>-2</sup></td></tr><tr><td>C<sup>E</sup><sub>55</sub></td><td>21.05 × 10<sup>9</sup></td><td>Nm<sup>-2</sup></td></tr><tr><td>e<sub>31</sub></td><td>-5.35</td><td>Asm<sup>-2</sup>(Cm<sup>-2</sup>)</td></tr><tr><td>e<sub>33</sub></td><td>15.78</td><td>Asm<sup>-2</sup>(Cm<sup>-2</sup>)</td></tr><tr><td>e<sub>15</sub></td><td>12.29</td><td>Asm<sup>-2</sup>(Cm<sup>-2</sup>)</td></tr><tr><td>ε<sup>σ</sup><sub>11</sub></td><td>8.14 × 10<sup>-9</sup></td><td>AsV<sup>-1</sup>m<sup>-1</sup>(Fm<sup>-1</sup>)</td></tr><tr><td>ε<sup>σ</sup><sub>33</sub></td><td>7.32 × 10<sup>-9</sup></td><td>AsV<sup>-1</sup>m<sup>-1</sup>(Fm<sup>-1</sup>)</td></tr></tbody></table>

# Appendix B

## Analytical 2D plane CFS for SS piezoelectric bimorph

The two piezoelectric coupled PDE were given in Eqs. (23) and (24); they read

$$S_{11}^{E} \psi_{zzzz} + (2S_{13}^{E} + S_{55}^{E}) \psi_{xxzz} + S_{33}^{E} \psi_{xxxx} - d_{31} \varphi_{zzz} - (d_{33} - d_{15}) \varphi_{xxz} = 0 \quad (94)$$

and

$$d_{31} \psi_{zzz} + (d_{33} - d_{15}) \psi_{xxz} - \epsilon_{11}^{\sigma} \varphi_{xx} - \epsilon_{33}^{\sigma} \varphi_{zz} = 0 \quad (95)$$

For sinusoidal loads, one can extend the method described by Timoshenko and Goodier [24] for isotropic material (see p.53ff., where a 2D analytical plane stress CFS in the form of Fourier series is presented) by considering the piezoelectric effects, see Eqs. (94) and (95). Hence, assuming a solution for the stress and electric potential functions as

$$\begin{aligned} \psi(x, z) &= h(z)\sin \alpha x \\ \varphi(x, z) &= g(z)\sin \alpha x \end{aligned} \quad (96)$$

One finds two coupled ordinary differential equations (ODE)

$$\begin{aligned} S_{11}^{E} h_{zzzz} - \alpha^2 (2S_{13}^{E} + S_{55}^{E}) h_{,zz} + \alpha^4 S_{33}^{E} h - d^{31} g_{zzz} + \alpha^2 (d_{33} - d_{15}) g_{,z} = 0 \\ d_{31} h_{zzz} - \alpha^2 (d_{33} - d_{15}) h_{,z} + \alpha^2 \epsilon_{11}^{\sigma} g - \epsilon_{33}^{\sigma} g_{,zz} = 0 \end{aligned} \quad (97)$$

For isotropic elastic materials, it is shown in Timoshenko and Goodier [24] that one has to treat a bi-quadratic ODE. Combining the relations in Eq. (97), one may find the following first-order matrix differential equation

$$\underline{\underline{y}}' = \underline{\underline{A}} \underline{\underline{y}}, \quad \underline{\underline{y}}^T = [h_{,zzz}, h_{,zz}, h_{,z}, h, g_{,z}, g]^T \quad (98)$$

A general solution of Eq. (98) is

$$y(z) = \sum_{j=1}^{6} C_j y_j(z) \quad (99)$$

where, $y_j$ are suitable eigenfunctions. For each layer, the six constants $C_j$ must fulfill the BCs, see Eqs. (57)–(63). This means that, for a bimorph, one has to calculate twelve constants from Eq. (99), and additional six constants for the displacement field, see Eqs. (55) and (56).
Appendix C

The stress and the electric field components in case of mechanical and electrical loading

C.1 Stress and the electric field as a function of the bending moment M(x) (mechanical loading due to q(x))

The elementary relations between the bending moment M(x), shear force Q(x) and the distributed load q(x) are given by

$$
\begin{align}
Q_x(x) &= -q(x) &\Leftrightarrow Q(x) &= -\int q(x) dx + A_1 \nonumber \\
M_x(x) &= Q(x) &\Leftrightarrow M(x) &= \int Q(x) dx + A_0 = -\int \int q(x) dx dx + A_1 x + A_0 \tag{100}
\end{align}
$$

Here $A_0$, $A_1$ are constants that are known from the boundary conditions at $x = 0$ and $L$ for a statically determinate beam. For statically in-determinate beams, the reader is referred to Irschik [26]. According to Eq. (25) the stress function $\psi$ is

$$
\begin{align}
\psi(x, z) &= \psi_1(x, z) + \psi_2(x, z) + \psi_3(x, z) + \dots \nonumber \\
&= \bar{\psi}_1(z)f(x) + \bar{\psi}_2(z)f_{,xx}(x) + \bar{\psi}_3(z)f_{,xxx}(x) + \dots \tag{101}
\end{align}
$$

The bending moment M(x), which is the sum of the bending moment for the upper M⁽²⁾(x) and lower layer M⁽¹⁾(x), is related to the axial stress by

$$
M(x) = M^{(2)}(x) + M^{(1)}(x) = b \int_{-c}^{0} \sigma_{xx}^{(2)} z \, dz + b \int_{0}^{c} \sigma_{xx}^{(1)} z \, dz \quad (102)
$$

In the case of mechanical loading (see Eq. (60)), one finds from the traction condition for the first iteration $\psi_1^{(1,2)}(z)$ (see also Eq. (33) when the unknowns $a_{01}$, $a_{11}$, $a_{21}$, $a_{31}$ are calculated)

$$
\bar{\psi}_1^{(2)}(-c) = 1/b, \quad \bar{\psi}_1^{(1)}(c) = \bar{\psi}_{1,z}^{(1)}(c) = \psi_{1,z}^{(2)}(-c) = 0 \quad (103)
$$

For the higher iterations, trivial traction conditions have to be satisfied for $\psi_2^{(1,2)}(z)$

$$
\bar{\psi}_1^{(1)}(c) = \bar{\psi}_1^{(2)}(-c) = \bar{\psi}_{1,z}^{(1)}(c) = \psi_{1,z}^{(2)}(-c) = 0 \quad i \ge 2 \quad (104)
$$

Considering Eqs. (58), (60) and (103), one finds that the function $f(x)$ equals the bending moment

$$
\begin{align}
M(x) &= f(x)b \left[ \int_{-c}^{0} \bar{\psi}_{1,z}^{(2)} z \, dz + \int_{0}^{c} \bar{\psi}_{1,z}^{(1)} z \, dz \right] = -f(x)b \left[ \int_{-c}^{0} \bar{\psi}_{1,z}^{(2)} dz - \bar{\psi}_{1,z}^{(2)} z |_{-c}^{0} + \int_{0}^{c} \bar{\psi}_{1,z}^{(1)} dz - \bar{\psi}_{1,z}^{(1)} z |_{0}^{c} \right] \nonumber \\
&= -f(x)b \left[ \begin{aligned}[t]
                 & [\bar{\psi}_1^{(2)}(0) - \bar{\psi}_1^{(1)}(0) - \bar{\psi}_1^{(2)}(-c) + \bar{\psi}_1^{(1)}(c) - \bar{\psi}_{1,z}^{(2)}(-c)c - \bar{\psi}_{1,z}^{(1)}(c)c] = f(x) \\
                 &= 0, \text{ see Eq. (58)} \end{aligned} \right. \nonumber \\
&= \frac{1}{b}, \text{ see Eq. (103)} \tag{105}
\end{align}
$$

For the higher iterations, it follows from Eqs. (42), (48), (50) and (100) that

$$
\begin{align}
f_2(x) &= f_{,xx}(x) = -q(x) \nonumber \\
f_3(x) &= f_{2,xx}(x) = f_{,xxxx}(x) = -q_{,xx}(x) \nonumber \\
\vdots \tag{106}
\end{align}
$$

It is noted that this relation holds only for mechanical loads (i.e. when the electrodes are short-circuited). For the axial stress, the shear stress and the electric field one finds by considering the product ansatz and substituting Eq. (105) into Eqs. (4) and (19)

$$
\left.
\begin{aligned}
\sigma_{xx}^{(1,2)} &= \psi_{,zz}^{(1,2)} = M(x)\bar{\psi}_{1,zz}^{(1,2)} + M_{,xx}(x)\bar{\psi}_{2,zz}^{(1,2)} + \dots \\
\sigma_{xz}^{(1,2)} &= -\psi_{,xz}^{(1,2)} = -M_{,x}(x)\bar{\psi}_{1,z}^{(1,2)} - M_{,xxx}(x)\bar{\psi}_{3,z}^{(1,2)} - \dots \\
E_x^{(1,2)} &= -\varphi_{,x}^{(1,2)} = -M_{,x}(x)\bar{\varphi}_1^{(1,2)} - M_{,xxx}(x)\bar{\varphi}_1^{(1,2)} - \dots \\
E_z^{(1,2)} &= -\varphi_{,z}^{(1,2)} = -M(x)\bar{\varphi}_{1,z}^{(1,2)} - M_{,xx}(x)\bar{\varphi}_{1,z}^{(1,2)} - \dots
\end{aligned}
\right\}
\quad \text{mechanical loading}
\tag{107}
$$

Eq. (107) shows the relation between the derived variables from both potentials and the distributed load $q(x)$, which is linked to the bending moment $M(x)$ or the shear force $Q(x)$ by Eq. (100).

C.2 Stress and the electric field as a function of the voltage distribution V (x) (electrical loading)

If the voltage is prescribed, then it follows from the product ansatz Eq. (34) and Eq. (62)

$$
f(x) = V(x) \tag{108}
$$

Inserting Eq. (108) into Eqs. (42), (48), (50), it follows for the first three iterations

$$
\begin{align}
f_1(x) &= V(x) \nonumber \\
f_2(x) &= f_{,xx}(x) = V_{,xx}(x) \nonumber \\
f_3(x) &= f_{2,xx}(x) = f_{,xxxx}(x) = V_{,xxxx}(x) \tag{109}
\end{align}
$$

and hence for the stress and the electric field (in a similar manner as for Eq. (107))
$$
\left.
\begin{aligned}
\sigma_{xx}^{(1,2)} &= \psi_{,zz}^{(1,2)} = V(x)\bar{\psi}_{1,zz}^{(1,2)} + V_{,xx}(x)\bar{\psi}_{2,zz}^{(1,2)} + \dots \\
\sigma_{xz}^{(1,2)} &= -\psi_{,xz}^{(1,2)} = -V_{,x}(x)\bar{\psi}_{1,z}^{(1,2)} - V_{,xxx}(x)\bar{\psi}_{3,z}^{(1,2)} - \dots \\
E_x^{(1,2)} &= -\varphi_{,x}^{(1,2)} = -V_{,x}(x)\bar{\varphi}_1^{(1,2)} - V_{,xxx}(x)\bar{\varphi}_1^{(1,2)} - \dots \\
E_z^{(1,2)} &= -\varphi_{,z}^{(1,2)} = -V(x)\bar{\varphi}_{1,z}^{(1,2)} - V_{,xx}(x)\bar{\varphi}_{1,z}^{(1,2)} - \dots
\end{aligned}
\right\}
\quad \text{electrical loading}
\tag{110}
$$

Eq. (110) shows the relation of the stress and the electric field as a function of the prescribed voltage V(x).

References

[1] Crawley EF. Intelligent structures for aerospace: a technology overview and assessment. AIAA J 1994;32(8):1689-99.

[2] Miu DK. Mechatronics: electro-mechanics and control mechanics. New York: Springer; 1993.

[3] Tzou HS. Multi-field transducers, devices, mechatronic systems and structronic systems with smart materials. Shock Vib Dig 1998;30:282-94.

[4] Reddy JN. On the generalization of displacement-based laminate theories. Appl Mech Rev 1989;42(11-Part 2):213-22.

[5] Reddy JN. On laminated composite plates with integrated sensors and actuators. Eng Struct 1999;21:568-93.

[6] Mitchell JA, Reddy JN. A refined hybrid plate theory for composite laminates with piezoelectric laminae. Int J Solids Struct 1995;32(16):2345-67.

[7] Tiersten HF. Linear piezoelectric plate vibrations. New York: Plenum Press; 1969.

[8] Krommer M. On the correction of the Bernoulli-Euler beam theory for smart piezoelectric beams. Smart Mater Struct 2001;10:668-80.

[9] Benjeddou A, Trindade MA, Ohayon R. A unified beam finite element model for extension and shear piezoelectric actuation mechanisms. J Intell Mater Syst Struct 1997;8:1012-25.

[10] Tiersten HF. Equations for the extension and flexure of relatively thin electroelastic plates undergoing large electric fields AMD 161, MD 42 In: Lee JS, Maugin GA, Shindo Y, editors. Mechanics of electromagnetic materials and structuresNew York: ASME; 1993. p. 21-34.

[11] Yang JS. Equations for the extension and flexure of electroelastic plates under strong electric fields. Int J Solids Struct 1999;36:3171-92.

[12] Batra RC, Vidoli S. Higher-order piezoelectric plate theory derived from a three-dimensional variational principle. AIAA J 2002;40(1):91-104.

[13] Krommer M, Irschik H. An electromechanically coupled theory for piezoelectric beams taking into account the charge equation of electrostatics. Acta Mech

2002;154:141-58.

[14] Benjeddou A, Trindade MA, Ohayon R. New shear actuated smart structure beam finite element. AIAA J 1999;37:378-83.

[15] Schoeftner J, Irschik H. A comparative study of smart passive piezoelectric structures interacting with electric networks: timoshenko beam theory versus finite element plane stress. Smart Mater Struct 2011;20:13.

[16] Schoeftner J, Buchberger G. Active shape control of a cantilever by resistively interconnected piezoelectric patches. Smart Struct Syst 2013;12:501-21.

[17] Schoeftner J, Buchberger G, Benjeddou A. Slender piezoelectric beams with resistive-inductive electrodes - modeling and axial wave propagation. Smart Struct Syst 2016;18:335-54.

[18] Schoeftner J, Buchberger G, Benjeddou A. Transverse dynamics of slender piezoelectric bimorphs with resistive-inductive electrodes. Smart Struct Syst 2016;18:355-74.

[19] Boley BA, Tolins IS. On the stresses and deflections of rectangular beams. ASME J Appl Mech 1956;23:339-42.

[20] Boley BA. The determination of temperature, stresses, and deflections in two-dimensional thermoelastic problems. New York Columbia University, Institute of Air Flight Structures; 1954. AF WADC Technical report No. 54-424.

[21] Krommer M, Irschik H. Boley's method for two-dimensional thermoelastic problems applied to piezoelastic structures. Int J Solids Struct 2004;41:2121-31.

[22] von Karman T. Über die Grundlagen der Balkentheorie 1927; vol. 7:3-10.

[23] Seewald F. Die Spannungen und Formänderungen von Balken mit rechteckigem Querschnitt 1927; vol. 7:11-33.

[24] Timoshenko SP, Goodier JN. Theory of elasticity. Singapore: McGraw-Hill; 1970.

[25] Ruan X, Danforth SC, Safari A, Chou TW. Saint-Venant end effect in piezoceramic materials. Int J Solids Struct 2000;37:2625-37.

[26] Irschik H. Enhancement of elementary beam theories in order to obtain exact solutions for elastic rectangular beams. Mech Res Commun 2015;68:46-51.