Approximate analytical solutions for piezoelectric rectangular beams by using Boley-Tolins method

The Airy stress function satisfies the bipotential equation which is a fourth order partial differential equation in two-dimensional elasticity for isotropic materials. In 1956 Bruno Boley published an iterative procedure for the solution of the Airy stress function. Nowadays this approximation method would be classified as homotopy analysis method (HAM). Boley presented analytical formulas for the deflections and the stresses of two-dimensional, isotropic, rectangular strips and beams. In this contribution Boley's method is extended by considering a piezoelectric material (PZT-5A). After defining the stress function and the electric potential and substituting the constitutive relations into the compatibility equation and Gauss’ law of electrostatics, one obtains two coupled partial differential equations. Adapting Boley's approximation method one finds an iteratively improving solution series where the higher iterations are considered as correction terms. In particular, it is interesting to note that the first iteration is identical to the one-dimensional piezoelectric beam theory within the framework of Bernoulli-Euler if the x-component of the electric displacement is neglected in the charge equation. Analytical formulas are presented for a piezoelectric rectangular beam (unimorph) and a piezoelectric bimorph in case of electrical actuation. Numerical results for the displacement and the electric potential are presented for a cantilever with various thickness-to-length ratios. The analytical outcome of this novel approach is compared to two-dimensional finite element results, showing a very good agreement even for very thick beams if the second iteration is considered.


Introduction
Piezoelectric materials are prominent members of so-called smart or intelligent systems.By means of the piezoelectric effect the motion of structures can be manipulated, i.e. the deflection or the stress.Prescribing the electric charge or the electric potential causes a change of the mechanical state (e.g.strain, stress); this is called converse piezoelectric effect.Contrary to this, the direct effect is exploited in sensor applications when the mechanical deformation causes a change of the electrical properties (e.g. of the electric charge or the electric field).An extensive summary on piezoelectricity can be found in the books by Crawley [1], Miu [2], Tzou [3] and Mura [4].
Active and passive vibration control methods can be efficiently developed if an accurate structural model is present in analytical form.Here we focus on the modeling of beams, which are the core components of e.g.aeroplane fuselages, vehicle bodies and frame bridges.Most developed models follow the principle that one assumes kinematic restrictions on the mechanical displacement.In an analogous manner these restrictions can be imposed on the electric potential or the electric displacement.The proper choice of an accurate polynomial order is essential for the accuracy of the outcome.For laminates made of elastic materials Reddy [5] assumes a power series approach in the thickness direction for each layer, hence containing the classical Bernoulli-Euler (BE) and Kirchhoff hypothesis and shear deformation theories (Timoshenko and Mindlin assumptions).In later contributions, these theories are extended for piezoelectric materials, see Reddy [6], Mitchell and Reddy [7].
Focusing on piezoelectric models, the simplest assumption on the electric field is to approximate the electric field component in thickness direction (normally denoted as E z ) by the quotient of the electric potential difference and the distance of the electrodes (=the layer thickness).In such a case the direct piezoelectric effect is disregarded, because the electric field components do not have to be computed by means of the charge equation of electrostatics.Strictly speaking for beam structures, the replacement of the differential quotient by the difference quotient is only correct if bending is negligible and membrane forces are dominant.
If the direct piezoelectric effect is neglected, this also called small piezoelectric coupling, Tiersten [8].For very thin beams within the BE framework, this electric field assumption yields accurate results only for thin beams, but not for thick beams.Even for thin beams the concept of the induced potential yields a small additional quadratic potential distribution due to the bending strain which also increases the bending stiffness of the structure, see Krommer [9] and Benjeddou et al [10].Krommer [9] derives in a systematic way that the electric potential is a second order polynomial in the thickness coordinate if Bernoulli-Euler (BE) kinematics and a simplified version of the charge equation (i.e. the x-component of the electric displacement is neglected) holds.The outcome is a fully coupled piezoelectric beam theory, which is in very good agreement to twodimensional finite element (FE) results computed in ABAQUS (static and dynamic analysis) as long as thin beams are considered.A further research from Krommer includes the presence of higher-order shear deformation theories or the influence of the axial component of the electric field/displacement (this is the electrical analogon to the Timoshenko theory), see Krommer and Irschik [11].These modelling approaches can be also extended to so-called passive structures, when the electrodes are connected to an electric circuit which dissipates energy or another energy conversion takes place from mechanical to electrical domain.These analytical models allow for the derivation of model-based and intelligent vibration control algorithms, see also Schoeftner and Irschik [12,13] and Schoeftner and Buchberger [14].
For very thick or laminated beams with strongly varying elastic moduli (e.g.sandwich beams with attached piezoelectric transducers) the derivation of validated structural models for the development of position tracking (i.e.shape control) and position sensing methods are either rare or not available in the open literature.It is the goal of this contribution to close this gap by presenting an easy-to-understand and efficient mathematical technique.This contribution is motivated by the work of Boley and Tolins [15], who present an iterative solution method for the biharmonic equation in order to compute the stress and the deflection of rectangular beams made of isotropic materials.From a mathematical point of view, the Boley-Tolins approach is a special perturbation method which enables the derivation of analytical, but approximate solutions for partial differential equations, see Yildirim [16] and Liao [17] and [18].The outcome of Boley-Tolins work is that the solution of the first iteration yields the Bernoulli-Euler result (e.g.concerning the deflection curve and also automatically contains the assumption of the linear axial stress distribution over the cross-section), whereas the higher order iterations mainly play a role for thicker beams because these terms depend on the thickness-tolength ratio.Hence, this theory goes even beyond the Timoshenko beam theory because the normal stress in the thickness direction is also considered.In Irschikʼs contribution [19] this approach was extended up to the fourth iteration and the author introduced the notion of actuating moments and actuating shear forces for the higher iterations in order to compute also statically indeterminate beam structures.Krommer and Irschik [20] extended the approximation method to the charge equation of electrostatics and derived more accurate results for the electric potential for a prescribed mechanical displacement field.
This contribution is an extension of Schoeftner [21] in the sense that an analytical solution procedure for the deflection, the stress potential and the electric potential is derived for a finite polynomial voltage or load distribution, i.e. the iteration ends after a finite number of iterations.Non-dimensional parameters for the geometry and the material constants are introduced and analytical results for voltage-actuated piezoelectric unimorph and bimorph cantilevers are presented.Taking into account the first and the second iteration only, the numerical results for the benchmark problem perfectly match with finite element results (computed by a modified code in MATLAB based on the FE-code of Ferreira [22]).Results for the vertical deflection and the electric potential for a thin (thickness-to-length ratio λ = 1/20) and for a moderately thick beam (λ = 1/5) are compared.A convergence study shows the influence of the thickness ratio and demonstrates that for moderately thick beams the presented approach provides an accurate model for piezoelectric strips and beams if the second iteration is considered.

Constitutive relations
The constitutive relations relate the stress σ, the strain ε, the electric field E and the electric displacement D. According to the d-form (=the strain-charge form, see [23]) they read e s s In (1) the vectors for the strain, the stress, the electric displacement and the electric field read e e e g e = = , , 2 2 The compliance matrix at constant electric field is S E , the piezoelectric coupling matrix is d, and the permittivity matrix at constant stress is ò σ Within the linear domain, the three strain components (2) depend on the displacement in the x-and the zdirections u(x, z) and w(x, z) e e g ( ) A dielectric medium satisfies Gauss' law.This reads for piezoelectric material, where free body charges vanish The final form of the coupled partial differential equations is obtained by inserting the stress function ψ(x, z) from (10) and the electric potential j(x, z) from (11) into the compatibility equation (13) and Gauss' law (15) y y y j j y y j j With the exception of some special boundaries and geometries, analytical solutions for ψ(x, z) and j(x, z) can be hardly derived for the two coupled partial differential equations in (16).In the following an iterative solution method from Boley-Tolins is adapted for a piezoelectric material (PZT-5A).In the original work of Boley-Tolins [15] the compatibility equation is iteratively solved for an isotropic material (note: the bipotential equation for isotropic materials is obtained from the first equation in (16) if one sets ) and d ij = 0, see Timoshenko and Goodier [24]).

Modification of Boley-Tolins iterative solution method
Boley-Tolins assume that the solution of the stress function, or at least a very good approximation of it, is a series of iterative solutions where ψ i (x, z) is the solution of the i th iteration.In a similar manner, we assume that this holds for the electric potential It is noted that for elastic materials, it does not matter if isotropic, anisotropic or functionally-graded materials are considered, it is sufficient to stop after I = 2 or 3 iterations, see Gahleitner [25] and [26] where the accuracy of the iterative procedure is validated by finite element calculations.Taking advantage of the approximate series solutions proposed in (17) and ( 18) and inserting them into (16), it follows y y y j j It is assumed that j i = ψ i = 0 holds if the subscript i is smaller than one.Anyway, this kind of solution splitting requires a short remark: A necessary condition that this solution procedure converges at all is the existence of a small perturbation parameter.In this case, the introduction of non-dimensional coordinates yields that the thickness-to-length ratio λ can be considered as a perturbation parameter (see (A.1)-(A.4) in A).This technique is known under the notion homotopy perturbation method (HPM) in the scientific literature, if a small physical parameter exists, see Yildirim [16] for an approximate solution of Poissonʼs equation.If this is not the case, a more advanced mathematical approach has been developed recently, which is called homotopy analysis method (HAM) that is a promising mathematical procedure to compute approximate solutions for nonlinear differential equations, see Liao [17] and [18].Anyway, for the fully coupled problem in piezoelectricity (19) and ( 20), the HPM-approach is studied first.
and for the third iteration y y y j j y y j j It is noted that the boundary conditions have to be satisfied by the solutions j 1 and ψ 1 from the first iteration (21) only.Assuming a distributed load q(x) at the upper surface (at z = − c) the traction conditions read s s s s If the potential of the upper surface is prescribed and the lower surface is grounded, the electrical boundaries are see also the single rectangular strip presented in section 4 (figure 1a).Exemplarily, we assume both a mechanical and an electrical actuation here: on the one hand the load distribution over the upper surface is q(x), on the other hand the potential difference is V(x).The boundary conditions (24) and (25) are only fulfilled by the solutions ψ 1 and j 1 of the first iteration.For the solutions of the higher iteration ψ 2 , ψ 3 , K and j 2 , j 3 , K the boundary conditions must vanish.Considering (21), ( 24) and (25) one finds that the first iteration yields y y y y j j j j The left solutions hold if the beam is only electrically actuated (i.e.V(x) ≠ 0 and M(x) = 0), whereas the right solutions hold for a short-circuited beam subjected to the vertical load q(x), which is related to the bending moment distribution M(x) by the fundamental relation M ,xx (x) = − q(x).Here y z again.Six coefficients of the functions y z 2 ¯( ) and j z 2 ¯( ) have to be determined by the second iteration a 02 , a 12 , a 22 , a 32 , b 02 , b 12 , which are determined by the trivial boundary conditions at z = ± c.For the third iteration one proceeds in a similar manner: the solutions (26) and (28) are inserted into (23), and one finds The coefficients a 03 , a 13 , a 23 , a 33 , b 03 , b 13 are fixed by the trivial boundary conditions.It is noted that for the second and for the third iteration, the formal structure of the solution is identical in case of a loaded beam where the derivatives of the voltage actuation have to replaced by the load distribution (i.e.: V ,xx (x) → − q(x) and V ,xxxx (x) → − q ,xx (x)), only the coefficients a i2 , b i2 and a i3 , b i3 will be different.

Calculation of the displacement
The mechanical state of stress and strain, the electric displacement and the electric field are known by (10) and (11).Integrating the first and the second equations in (9) with respect to x and z respectively, the deflection follows as ò ò where F(z) and G(x) are arbitrary functions.Assuming that the x-dependency of the load actuation V(x) is a polynomial, one expects the functions F(z) and G(x) also to be polynomials = + Inserting the expressions for the displacement (30) into the equation for the shear strain (=third equation in (9)) one finds ò ò ≔ ( ) Considering the compatibility equation (12) and the right-hand side of (32), one finds that h(x, z) is a univariate polynomial in x and z (but not a multivariate polynomial otherwise (12) would be violated) Substituting (31) and (33) into (32), one finds by comparison of coefficients From (34) all the coefficients of the function F(z) and G(x) in (31) are known, except for f 0 , g 0 and g 1 .These yet unknown three constants have to be specified by the kinematic boundary conditions, where three node restrictions for a statically determinate structure are required to prevent rigid body motions, see Timoshenko and Goodier [24].From (30) it follows for the axial deflection ò ò y j y j y Considering the piezoelectric unimorph in figure 1a, the remaining unknowns f 0 , g 0 , g 1 from F(z) and G(x) (35) and (36) are determined by the kinematic boundary conditions If the structure is not electrically loaded, but subjected to the lateral load q(x), V(x) needs to be simply replaced by the bending moment M(x) in (35) and (36).

Benchmark examples-Unimorph and bimorph piezoelectric beams
In this section the analytical results obtained by the Boley-Tolins method are verified.First a piezoelectric unimorph cantilever is investigated (Figure 1a, section 4.1), then a bimorph cantilever is investigated (figure 1b, section 4.2).
The material properties (PZT-5A) and the geometry for the examples are listed in table 1.In order to give analytical results in a reduced form it is useful to introduce non-dimensional parameters for the geometry (i.e. the thickness ratio λ = 2c/L) and for the material constants β, π s , π p , π d , π ò and ν 13 , whose meaning is given in table 2. For the computation of the three kinematic unknowns f 0 , g 0 and g 1 (for the functions F(z) and G(x) in (30)) and those for the stress function and the electric potential (a ji , b ji , see (26)-( 29)) the symbolic software MATHEMATICA is used for evaluating the analytical expressions.Then these results are compared to the results from our own finite element code in MATLAB.A Q8-element (eight-node quadrilateral element with quadratic ansatzfunctions for the displacement components and for the voltage) is used for the finite elements.The modified structure of this FE code is based on the book by Ferreira [22], who presents a huge variety of standard finite elements (beams, plates and two-dimensional elements) for elastic structures.In our calculation this code is adapted and piezoelectric properties are properly taken into account, the reader is referred to Piefort [27].

Example-unimorph piezoelectric cantilever
For the first numerical example the thickness-to-length ratio of the piezoelectric unimorph is λ = 1/5.At the upper surface the voltage actuation is j . Hence, according to (26)-( 29), the first, the second and the third iteration exist (V(x) ≠ 0, V ,xx (x) ≠ 0, V ,xxxx (x) ≠ 0), but the fourth iteration vanishes (∂ 6 V/∂x 6 (x) = 0).Taking advantage of the non-dimensional coordinates and the parameters from table 2 the results for the displacement field components u = u 1 + u 2 + K and w = w 1 + w 2 + K read (see ( 35) and (36)) The electrical potential reads

permittivity ratio
It is noted that the results from the third iteration do not vanish (i.e.u 3 ≠ 0, w 3 ≠ 0 and j 3 ≠ 0 hold), but these are lengthy expressions that are not explicitly written down due to reasons of readability.From a practical point of view the results will show that the influence is small (see figure 2 and 3).The lateral deflection of the unimorph is depicted in figure 2 (cross section A-A at x = L/2 and cross section B-B at x = L), the two-dimensional deformation plot is shown in figure 3.
The first iteration yields (blue) w 1 (ξ, η) = 0, see (39) (figure 2(a), (b)), i.e. bending deflections can not be realized with a piezoelectric unimorph.The second iteration w 2 causes a low order deflection and is nearly a linear function in z or η, see (39).One recognizes from figure 2(a) and b that there is no significant deviation between the FE result (black) and the analytical result if the second iteration (red) is included.Considering also the third iteration (cyan), one finds an even better agreement.At ξ = 0.5 (cross-section A-A, see figure 2(a)) one observes an almost linear distribution in the thickness direction: − 1.625 × 10 −10 m (at the upper surface η = − 1) and 1.881 × 10 −10 m (at the lower surface η = 1).This is very close to the FE outcome: − 1.624 × 10 −10 m (at η = − 1) and 1.886 × 10 −10 m (at η = 1).The third iteration still improves the results from the extended Boley-Tolins procedure: − 1.623 × 10 −10 m (at η = − 1) and 1.883 × 10 −10 m (at η = 1).At ξ = 1 (cross-section B-B, see figure 2b) one recognizes that the lateral deflection does not remain constant over the thickness (note that w 2 does not depend on η if ξ = 1, see (39)): there is an additional higher order distribution from the third iteration w 3 (ξ, η) ≠ 0, but with a low amplitude.Figure 3 shows four twodimensional plots: the FE-result (a) and the analytical results up to the first, the second and the third Boley iteration (b, c, d).One observes that the first iteration only reflects the axial deformation, which is uniformly distributed across the cross-section u 1 (ξ = 1, η) ≈ − 6.84 • 10 −10 m (figure 3(b), one observes that the horizontal deflections are the same for locations 3, 4 and 5).This result also follows if one integrates the simplified relation for the axial strain u ,x ≈ d 31 E ,z by neglecting the axial and the compressive stresses (i.e.σ xx ≈ σ zz ≈ 0).Hence the axial deflection at the free end follows as . It is evident from the deformation plot that the first iteration result does not correspond to the FE result (compare figure 3(a) and (b)), hence the higher order iterations have to be taken into account (see ( 38) and (39)).Considering also the second iteration one finds a very good correlation to the finite element result.One observes that the vertical deflection remains constant at the free end 2.05 × 10 −10 m (cf the deflections for the locations 3, 4 and 5 and the result B-B in figure 2(b) (red) do not depend onη).As already mentioned, the third iteration yields a small nonlinear distribution, so that the analytical solution is very close to the FE result (see also figure 2(b)).
Summarizing the results of this section in a few words one may conclude that the third iteration is practically negligible although a very thick beam (λ = 1/5) is considered.Considering the first and the second Boley-Tolins iteration, one may qualitatively find the same outcome as that from the numerical simulation (FE).

Example-bimorph piezoelectric cantilever
In this section analytical results are presented for a more practical example: a piezoelectric bimorph which is actuated by the voltage load 1b).According to section 2.3, six unknowns have to be determined for each layer and each iteration: six boundary conditions are given by ( 24) and (25), the remaining six equations are due to the continuity relations for the displacement, for the stress and for the electric voltage.The continuity relations require the stress continuity and the no-slip condition, if perfect bonding is realized: i.e. on the one hand, the shear and the transversal stress components of the upper and lower layer must match.On the other hand, for the axial and the lateral deflection at z = 0, interlaminar slip does not occur (i.e.Δu = u (low) (ξ, 0) − u (upp) (ξ, 0) = 0 and Δw = w (low) (ξ, 0) − w (upp) (ξ, 0) = 0).Furthermore, the voltage actuation over the upper and over the lower surface is j , but vanishes at the interface j(x, z = 0) = 0.
In particular, the influence of the thickness ratio is investigated, numerical results are presented for λ = 1/20 (thin beam, figure 4) and λ = 1/5 (thick beam, figure 5).Finally a convergence study compares the various calculations (FE results, Boley 1st iteration and Boley 1st and 2nd iteration, in particular for the tip deflection and the electric potential) as a function of the thickness-to-length ratio λ = {1/2,K,1/40}, ranging from very thick to very thin beams (figure 6).

Analytical solutions for the deflection curve and the electric potential
The deflection curve (=lateral deflection at centerline z = 0) reads for the first and for the second iteration Unlike as for the piezoelectric unimorph, the first iteration yields a non-zero deflection (cf the unimorph solution (39) and the bimorph solution (41)).Furthermore, it is noted that the result from the first iteration does not depend on the thickness-coordinate η, but the second iteration of the adapted Boley-Tolins method also includes the thickness deformation ).The thickness deformation is small, but for multi-layered beams like sandwich-structures with a low core stiffness it might become significant.As already stated the first iteration yields the result of a piezoelectric beam based on Bernoulli-Euler assumptions.From literature, the deflection curve for a voltage-loaded piezoelectric bimorph reads (see [9] and [29]) where A 0 and A 1 depend on the kinematic boundary conditions.In Appendix B it is shown that the results from the first Boley iteration (41) and that obtained by considering BE-assumptions (43) really match.) and (d) electric potential j(ξ = 0.75, in the thickness direction. The results for the electrical potentials from the first and the second iteration are given by (44) and (45)  ) and (d) electric potential j(ξ = 0.75, η) in the thickness direction.As expected Boleyʼs first iteration (blue) hardly deviates from the other solutions, the deviation from w 1 2 (red), when the second iteration is considered, is only 0.5%.The FE result (black) coincides with w 1 2 .Figure 4(c) shows the absolute deflection over the cross section at half of the beam length ξ = 0.5.One recognizes that the through-thickness variation of the deflection is almost identical for Boleyʼs 2nd iteration (red) and the FE result (black), whereas the first Boley iteration (blue) does not include the thickness deformation at all.The thickness variation is a lower order effect, which is practically insignificant for thin beams.At ξ * = 0.5 the difference Δw(ξ * , η = 1) = w(ξ * , η = 1) − w(ξ * , η = 0) ≈ 3.1 × 10 −10 m is only 0.6% of the midline deflection.Alternatively, a rough approximation of the thickness deformation is found by neglecting the normal stresses in the axial and the thickness direction .The analytical results from the extended Boley-Tolins method and the FE-results (black) are almost indistinguishable at ξ = 0.75 and η = 0.5: 0.303 V computed by FE and second order Boley and 0.3032 V by the first order Boley (=BE result).

4.2.3.
Comparison of the Boley-Tolins method to finite element results for a thick beam (λ = 1/5) For the thick beam (λ = 1/5) the results are shown in figure 5.One observes the good agreement between the FE results (black) and the outcome of the second order Boley results (red): the tip deflection difference is less than 0.1% (figure 5(b)).The variation over the thickness for upper and lower layer is almost linear again which is the outcome from the second iteration.For comparison, the difference between first order Boley and FE results is 0.8% for the tip deflection.In the middle of the layer at η = ± 0.5, the result for the electrical potential (FE and second order Boley) are almost equal j 12 (0.75, 0.5) ≈ j FE (0.75, 0.5) ≈ 0.2999 V, whereas the value from the first Boley iteration yields j 1 (0.75, 0.5) = 0.3032 V.
4.2.4.Convergence study-vertical deflection and electric potential as a function of the thickness ratio Finally, a convergence study is presented in figure 6 showing the relative deviation of the first and the second order Boley-results as a function of the inverse of the thickness-to-length ratio 1/λ.The relative deviation for the tip deflection reads

=
´- and for the electric potential ´- Here a thin beam is classified if 1/λ > 30, a moderate thick beam if 10 < 1/λ 30 and a thick beam if 1/λ 10 holds.One observes that for 1/λ = 10, the relative deviation is 0.2 % for the first iteration.Even for a very thick beam (1/λ = 3) the error is only − 0.34 % if the second Boley iteration is considered.For thinner beams the results from the second Boley-Tolins iteration converge to the basic solution.It is noted that in our calculation the analytical results slightly deviate from the FE-result, but it is only 0.02 % for 1/λ = 40.This is expected because it is well-known that FE-results underestimate the analytical results (depending on the number of finite elements, ansatzfunctions, aspect ratio, etc used in the calculation).It is noted that the accuracy of twodimensional FE results (note that a standard Q8-element is used which is also available in commercial software packages like ANSYS and ABAQUS) strongly depends on the number of elements in thickness direction.In particular, this becomes a problem for very thin beams if FE results are used for verification: even for a thin beam (i.e.λ = 1/40) 120 elements in thickness direction are used in order to achieve the desired accuracy.A lower number of elements means an underestimation of the deflection.The authors are aware of this discrepancy between numerical and analytical results which is present for very thin beams, but in this study we are mainly interested in verifying the Boley-Tolins procedure for thick structures and not in the development of highly convergent finite elements.For the electric potential the relative difference has the same order of magnitude 0.27 % for a moderately thick beam (i.e.1/λ = 10).The FE calculation, whose main purpose is to verify the presented theory, indeed verifies the results.Even for a very thick beam (1/λ = 3) the error is only 0.02 % for the electric potential.

Conclusion
In this contribution an approximation method is presented to solve the compatibility equation and the charge equation for two-dimensional piezoelectric structures.This method was originally presented by Bruno Boley in 1956 to find an iteratively converging solution series of the bipotential equation for an isotropic material.Here the method is adapted if piezoelectric structures are considered in order to find a series solution for both the Airy stress function and the electric potential.Inserting the constitutive relations into the compatibility equation and the charge equation of electrostatics, one finds two coupled partial differential equations.Solving these differential equations according to the adapted Boley-Tolins method simultaneously, one may find iteratively converging series of solutions.Analytical solutions are presented for a piezoelectric cantilever (unimorph and bimorph) and compared to the solutions obtained by two-dimensional finite element analysis.It is found that the basic solution (from the first iteration) corresponds to analytical solutions already known from the literature, if Bernoulli-Euler assumptions hold for a piezoelectric bimorph.In case of thick beams, the second iteration is a correction term which can be interpreted as a low order correction for piezoelectric materials.By varying the thickness of the beam a convergence study shows that the presented analytical procedure provides accurate, higher order corrections, which may be significant for thick structures in position control applications and health monitoring.From (A.4) one observes that the small-perturbation parameter is the thickness-to-length ratio λ.For a rectangular strip the relation λ < 1 holds, which is a necessary, but not a sufficient condition that the presented perturbation method in section 2.3 might work.In case of very thin structures λ → 0, one immediately finds the first iteration from (A.4), cf.

1 ¯( ) and j z 1 ¯equation ( 27 )
( ) are polynomials in z.Inserting(26) into(21) it follows for the polynomials y z 1 ¯( ) and j z 1 holds for a mechanical or for an electrical load.The second subscript in a i1 or b i1 denotes the first iteration.For the loaded beam with q(x) ≠ 0, the voltage may be simply replaced by the bending moment: V (x) → M(x).The unknown six coefficients a 01 , a 11 , a 21 , a 31 , b 01 , b 11 are computed by the six boundary conditions(24) and(25).For the second iteration the solutions of the first iteration (26) are substituted into(22) and it follows

Figure 3 .
Figure 3. Undeformed (black mesh grid) and deformed unimorph.The axial and the vertical deflection are given for the location 1-5.((a): FE result, (b) only 1st Boley-Tolins iteration, (c) 1st and 2nd iterations (d) 1st, 2nd and 3rd iterations are considered, note: the scale factor is chosen such that the maximal absolute deflection is 0.05).

Figure 6 .
Figure 6.Convergence study-relative deviations of the tip deflection and the electric potential to FE-results (see (50) and (51)).
differential equations(16), which directly follow from the compatibility and the charge equation Analytical results from literature for a piezoelectric bimorph Appendix B.1.Deflection curve w(x) In Schoeftner and Benjeddou[21] the deflection curve for a voltage-loaded piezoelectric bimorph can be evaluated by integrating e 31 ˜and K M are the piezoelectric coefficient on beam level (in stress-charge form) and the bending stiffness.In Krommer[9] the same equation is written in an alternative form is the piezoelectric moment and D ¯is related to the bending stiffness.The bending stiffness is Between the material parameters in strain-charge form and those parameters on beam level (i.e.denoted by the tilde sign), the following relations hold (B.1), one obtains the same analytical outcome for the lateral deflection as the result obtained by the first Boley iteration (41).
the stress function ψ(x, z), whose second partial derivatives are related to the stress components, and the electric potential j(x, z), whose partial derivatives are related to the electric field, are introduced

Table 2 .
Non-dimensional parameters for the numerical examples.
Evaluating the last part of (49) at ξ = 0.5 yields Δw approx ≈ 3.27 × 10 −10 m, which agrees quite well with the aforementioned result.The electric potential over the cross-section is shown in figure 4d.Additionally, the linear voltage distribution for the electric potential is drawn (black dash-dot).The linear distribution is given in (46),