The powerful extended Kantorovich method (EKM) originally proposed by Kerr in 1968 is generalized to obtain a three-dimensional coupled piezoelasticity solution of smart piezoelectric laminated plates in cylindrical bending. Such solutions are needed to accurately predict the edge effects in these laminates under electromechanical loading. The Reissner-type mixed variational principle extended to piezoelasticity is used to develop the governing equations in terms of displacements, electric potential as well as stresses and electric displacements. It allows for exact satisfaction of the boundary conditions, including the non-homogeneous ones at all points. An n-term solution generates a set of 11n algebraic ordinary differential equations in the inplane direction and a similar set in the thickness direction for each lamina, which are solved in closed form. The multi-term EKM is shown to predict the coupled electromechanical response, including the edge effects, of single-layer piezoelectric sensors as well as hybrid laminated panels accurately, for both pressure and electric potential loadings. This work will facilitate development of accurate semi-analytical solutions of many other unresolved problems in three-dimensional piezoelasticity, such as the free-edge stresses in hybrid laminates under bending, tension and twisting.
The presence of widely different mechanical, electric and thermal material properties in adjacent layers in piezoelectric laminated structures and the geometric discontinuity at the boundaries result in a stress field near the edges, which is truly three dimensional in nature and decays rapidly away from the edges. Because of piezoelectric coupling effects, the edge effect is a more complex phenomenon in the hybrid laminates than their elastic counterparts. These localized stresses are of major concern because they may result in premature interlaminar failure, and reduction or complete loss of actuation/sensing authority of the piezoelectric layers. The need for accurate analytical/semianalytical solutions based on three-dimensional piezoelasticity/elasticity for such stress fields in laminated structures has been felt for a long time, but has remained largely unsolved so far , understandably so because of the inherent complexities involved in the problem. Here, three-dimensional piezoelasticity solutions refer to those wherein the governing equations of three-dimensional piezoelasticity are used to solve for the distributions of the field variables along the thickness direction, in contrast to the two-dimensional theories, which make a priori assumptions on the through-thickness distributions of the field variables, and the governing equations are obtained in terms of the inplane directions only.
Analytical three-dimensional solutions for laminates with arbitrary edge conditions and electric boundary conditions can provide valuable insight into the edge effects and the influence of electromechanical coupling on them. The limitations of the conventional three-dimensional finite-element (FE) method in accurately predicting the stress field near some edges are well known in the literature [1,2]. A three-dimensional semi-analytical FE solution based on an FE discretization in inplane coordinates has been presented by Qing et al.  for static and free vibration analysis of hybrid piezoelectric plates. The method, however, inherits the limitations of the FE method in predicting the stress field with sharp gradients, leading to possible numerical instabilities . Vel & Batra used the Eshelby–Stroh formalism to obtain a series solution for the cylindrical bending  and general bending  of hybrid plates with arbitrary edge conditions. However, the interface continuity and the boundary conditions were satisfied in the sense of Fourier series, resulting in an infinite system of equations in an infinite number of unknowns. The truncation in the series led to unwanted oscillations in the variations of transverse stresses at the boundaries. The last two remain the only three-dimensional approximate analytical solutions for bending of hybrid piezoelectric plates presented to date.
In the year 1968, Kerr  proposed a very elegant method, called the extended Kantorovich method (EKM), to obtain semi-analytical solutions of two-dimensional elasticity problems involving bivariate partial differential equations. This method assumes the solution in the form of a sum of products of separable functions in the two directions and alternatively considers the functions in one direction as known (from an initial guess or previous iteration) and reduces the problem to a set of ordinary differential equations (ODEs) for the functions in the other directions. The iterations over the two directions are repeated until the convergence is achieved, which occurs very rapidly. Unlike other approximate methods such as Ritz’s and Galerkin’s methods, the accuracy of the final solution is independent of the initial guess, whether or not the guess satisfies the boundary conditions. The method has been successfully used to analyse several problems of single-layer [8–11] and multi-layer [12–14] elastic plates and shells. Lucy & Abramovich  have employed the method to obtain the bending response of rectangular symmetrically laminated piezoelectric plates using the classical laminate theory, but without considering two-way electromechanical coupling.
All the above problems are, however, based on two-dimensional theories and involve only homogeneous boundary conditions. Even after 43 years of its inception, the EKM had not been applied to three-dimensional elasticity problems until recently when the authors extended it to the three-dimensional elasticity solution for cylindrical bending of laminated plates [16,17]. The solution involved satisfaction of non-homogeneous boundary conditions due to applied traction at the top and bottom surfaces. It was shown that a one-term solution, which is known to yield reasonably accurate results for two-dimensional plate/shell problems, fails to accurately predict the stress field near the boundaries in the three-dimensional elasticity solution. However, just one or two additional terms led to a very accurate prediction of all entities, reconfirming the efficacy of the method. Tahani & Andakhshideh  have subsequently presented a three-dimensional elasticity solution for general bending of rectangular elastic laminated plates. But, since a displacement formulation is adopted, all boundary and interface conditions involving stresses, including those of applied pressure at the top surface, are satisfied only in an average sense over the surface and not at all points. Clearly, the stress field near the edges will be highly sensitive to such approximations and hence will not be reliable. Kim et al.  presented a two-dimensional analysis of free-edge effects in a composite laminate using the EKM, which also involved only homogeneous boundary conditions.
To the best of the authors’ knowledge, no EKM solution has been presented so far for the coupled electromechanical problem of piezoelectric single- or multi-layered plates, even using two-dimensional plate theories. In this work, the multi-term EKM for the three-dimensional elasticity solution  is extended to the three-dimensional piezoelasticity case of cylindrical bending, considering full electromechanical coupling. Using the method, a semi-analytical solution is obtained for the angle-ply hybrid panels, integrated with piezoelectric fibre reinforced composite (PFRC) actuators and sensors. An extended Reissner-type variational principle for the piezoelastic medium is used to develop a mixed formulation with displacements, stresses, electric potential and electric displacements as independent variables. This allows for the exact satisfaction of the boundary and interface conditions at all points, and ensures the same degree of accuracy for all response entities. Unlike in the elastic case wherein the same form of constitutive equations was used in both inplane and thickness directions, the present formulation requires two different forms of the constitutive equations involving the electric variables. Unlike the pressure load, which is applied only at the top and bottom surfaces, the electric potential is prescribed also at the interfaces, calling for a layerwise solution for the part satisfying the non-homogeneous boundary conditions. It also leads to discontinuity in the transverse electric displacement at the layer interfaces, which is accounted for. The formulation allows for both open and closed circuit electric boundary conditions at the outer surfaces of the piezoelectric layers so that they can be modelled as sensors as well as actuators.
2. Extended Reissner-type variational principle for piezoelasticity
Consider an infinitely long, hybrid angle-ply composite panel (figure 1) with span length a along direction x and total thickness h. The hybrid laminate is made of L perfectly bonded generally orthotropic laminas, some of which can be of monolithic piezoelectric or PFRC materials with poling along the thickness direction z, which can be used as distributed sensors and actuators. The piezoelectric materials are considered to be of class mm2 symmetry with respect to the principal material axes xi (i=1,2,3), x3 being oriented along the z-direction. The mid-plane of the plate is chosen as the reference xy-plane. The principal material axis x1 of the kth layer from the bottom is at an angle θk to the x-axis. The thickness of the kth layer is t(k), and the z-coordinate of its bottom surface is denoted as zk−1. The interface between the kth and the (k+1)th layer is denoted as the kth interface. The layer superscript is omitted unless required for clarity.
The three-dimensional piezoelasticity constitutive equations relating stress components (σi and τij), electric displacement components (Di) to the corresponding strain components (εi and γij) and electric field components (Ei) with respect to the structural coordinate system (x,y,z) are given by 2.1where and denote the transformed elastic compliances, piezoelectric strain constants and dielectric permittivities at constant stress field, respectively. They can be expressed in terms of the engineering material constants, in the material coordinate system, namely, Young’s moduli Y i, shear moduli Gij, major Poisson ratios νij and constant strain dielectric permittivities ηij. All entities are non-dimensionalized as follows such that the dimensionless forms of the governing equation are exactly identical to the original equations: 2.2where Y 0 and d0 are, respectively, the values of Young’s modulus and the piezoelectric coefficient, used for non-dimensionalization. This non-dimensionalization is performed to ensure numerical stability in the computations. In the following, the dimensionless form of the equations is used. However, since the resulting equations are identical to their original counterparts, the superscript asterisks (*) are omitted from the dimensionless entities for convenience.
For the generalized plane strain case of cylindrical bending, considering that all variables are independent of y, the strain–displacement and electric field–potential relations reduce to 2.3where a subscript comma denotes differentiation. Using equations (2.3)2, and (2.1)2 for ϵy in equations (2.1)1, (2.1)3, (2.1)6 and (2.1)9, εx,εz,γxy and Ez can be expressed in terms of σx,τxy,σz and Dz as 2.4where 2.5Here, for example, equation (2.3)2 denotes the second equation in the set of equations given by equation (2.3). The strains εx,εz and γxy can alternatively be expressed in terms of σx,τxy,σz and Ez as 2.6where 2.7While equation (2.4) is used for obtaining equations in the z-direction, equation (2.6) is used for the x-direction. Using equations (2.1)4, (2.1)5 and (2.1)7, γyz, γzx and Ex can be expressed in terms of τyz, τzx and Dx as 2.8where 2.9
The Reissner-type mixed variational principle  can be extended for a piezoelectric medium without body force and internal charge source for the cylindrical bending case as 2.10∀δui, δσi, δτij, δϕ and δDi, where V denotes the volume of the panel per unit width in the y-direction. Equation (2.10) implies that the associated variationally consistent boundary conditions, as stated below, are to be satisfied exactly, 2.11where are the components of the prescribed surface traction , and are the components of the prescribed electric displacement on a surface with outward normal , with (i=1,2,3) being the unit vectors along the x,y,z-directions. AT and Au are the surface boundaries where surface traction and displacements are prescribed, respectively. Similarly, AD and Aϕ are the surface boundaries where electric displacement and potential are prescribed, respectively.
3. Boundary and interface conditions
A dimensionless inplane coordinate ξ is defined as ξ=x/a, which takes values 0, 1 at x=0,a; and a local thickness coordinate ζ(k) for the kth layer is defined as ζ(k)=(z−zk−1)/t(k), taking values 0, 1 at z=zk−1,zk. The plate is subjected to uniform pressures p1 and p2 at the bottom and top surfaces, respectively. These surfaces may be under either prescribed electric potential ϕα or prescribed transverse electric displacement (charge density) Dα, which are uniformly distributed over x. Thus, the boundary conditions at z=∓h/2 are 3.1The interfaces between the piezoelectric layers and the elastic layers are taken as grounded for effective actuation/sensing. For such interfaces nq (q=1,…,La), the electric potential is prescribed as 3.2The equilibrium and the compatibility conditions at the kth interface between adjacent layers are 3.3for k=1,…,L−1, except for Dz for the interfaces k=nq, q=1,…,La. The Dz is discontinuous for the interfaces nq. For such surfaces, the continuity condition for Dz is to be replaced by equation (3.2). The boundary conditions at the edges ξ=0 and 1 can be arbitrary, for example, 3.4Electrically, the ends may be under the open circuit (OC) condition with Dx=0 or under the closed circuit (CC) condition with prescribed potential.
4. The generalized extended Kantorovich method
The solution of the field variables X=[u v w σx σz τxy τyz τzx ϕ Dx Dz]T is assumed in terms of an n-term series of the products of separable functions of the two coordinates ξ and ζ. To take care of the non-homogeneous boundary conditions for σz and ϕ as given in equation (3.1), a solution satisfying the non-homogeneous part is superimposed to the above solutions for σz and ϕ, so that the latter are required to satisfy homogeneous boundary conditions only. Thus, the general solution of the lth variable Xl of X for the kth layer takes the following form: 4.1where and are unknown univariate functions of ξ and ζ to be determined iteratively, satisfying homogeneous boundary conditions. The repeated index l does not mean summation here. The second term in the general solution is added to satisfy the non-homogeneous boundary conditions on σz, where δlm is Kronecker’s delta, pa=−(p1+p2)/2 and pd=−(p2−p1)/h. The solution for the non-homogeneous boundary condition of ϕ is given by 4.2While functions are defined for the kth layer, ’s are valid for all layers.
(a) First iteration step
In this step, functions are assumed, for which the variation δXi is given by 4.3Functions are partitioned into a column vector of eight independent variables that appear in the boundary and interface conditions (3.1) and (3.3), and a column vector consisting of the remaining three dependent variables, 4.4Now, εi, γij and Ei in the variational equation (2.10) are substituted with their expressions from equations (2.4) and (2.8). Using the solution given by equations (4.1) and (4.3) in the resulting equation, integrating over the ξ-direction and considering that the variations are arbitrary, the coefficients of are set to zero. This results in the following set of 11n differential–algebraic equations for for the kth layer: 4.5and 4.6where and are 8n×8n, 8n×8n, 8n×3n, 3n×3n and 3n×8n matrices, respectively. Their non-zero terms are defined using the notation for integration over the span length a as 4.7where ip=n(p−1)+i and jq=n(q−1)+j for p,q=1,2,…,11. , and are load vectors of size 8n, 3n and 8n, respectively, whose non-zero terms are 4.8where . Since the functions are known analytical functions, elements of the matrices defined in equations (4.7) and (4.8) can be evaluated in closed form. Solving the system of algebraic equations (4.6) for and substituting back the solution into equation (4.5) yields 4.9with and . Equation (4.9) represents a set of 8n non-homogeneous first-order ODEs with constant coefficients. Following the procedure given by Kapuria & Kumari , its general solution can be expressed in terms of 8n real constants as 4.10where Fi(ζ) are column vectors of functions expressed in closed form using exponential and trigonometric functions in terms of the ith eigenvalue and eigenvector of matrix A. U0 and U1 correspond to the particular solution.
The boundary and interface conditions in equations (3.1)–(3.3) can be written in terms of functions as 4.11and 4.12where k=1,2,…,L−1, except for at the interfaces k=nq (q=1,…La) where ϕ is prescribed. At these interfaces, the continuity condition on is replaced by 4.13The 8n×L constants ’s for L layers are obtained from the 8n boundary conditions and 8n×(L−1) interface continuity conditions given by equations (4.11)–(4.13). The initial trial functions for for the first iteration are assumed as 4.14which correspond to simply supported boundary conditions at both ends.
(b) Second iteration step
In this step, the solution of the previous step is taken as the known approximate solution for , while functions are considered unknown. The variation δX for this case is given by 4.15The functions are partitioned into vector of eight independent variables that appear in the boundary conditions (3.4) and a vector of the remaining variables that are dependent, 4.16In this case, εi, γij and Ei in the variational equation (2.10) are substituted with their expressions from equations (2.6) and (2.8). Substituting equations (4.1) and (4.15) in the resulting equation, performing integration over the ζ-direction on the known functions of ζ, applying integration by parts wherever necessary, and equating the coefficient of to zero individually, yields the following system of differential–algebraic equations for : 4.17and 4.18where and are 8n×8n, 8n×8n, 8n×3n, 3n×3n and 3n×8n matrices, respectively, and and are column vectors of size 8n and 3n, respectively, comprising the loading terms. Using the notation for integration across the thickness, the non-zero elements of the above-mentioned matrices are given by 4.19and 4.20Since the functions are known in closed form from equation (4.10), the integrations for the above terms are also evaluated in closed form. Substituting from equation (4.18) into equation (4.17) yields the following set of first-order ODEs for : 4.21with and . Equation (4.21) is of the same nature as equation (4.9) and is solved similarly satisfying the prescribed boundary conditions (3.4), expressed in terms of . This completes the second iterative step.
These two iteration steps (§4a,b) of alternating directions for approximate and analytical solutions are repeated until a prescribed level of convergence is achieved. For higher-order terms in the trial function, the matrices A and B in equations (4.7) and (4.19) often have large positive eigenvalues, whose exponentials, being very large, may cause numerical instability to the solution. To prevent this, for positive eigenvalues, the exponential term eαζ is replaced with eα(ζ−1), which limits the maximum value of the exponential to one at ζ=1 and decays towards ζ=0. Obviously, this does not alter the final solution, since the constants associated with these eigenvectors get modified accordingly. This method has been used by many researchers for improved numerical stability in such solutions .
5. Numerical results
Numerical results are presented for a single-layer piezoelectric panel (a) of PZT-5A, which can be used in sensory applications, and a hybrid soft-core sandwich panel (b) integrated with PFRC sensor and actuator layers at the top and bottom surfaces, as shown figure 2. The properties of the materials are listed in table 1.
The response of the PZT-5A panel (a) is obtained for a uniform pressure loading applied at the top surface. For the hybrid panel (b), the response is obtained for the following pressure and electric potential load cases:
— uniform pressure applied on the top surface, p2=p0 and
— uniform actuation potentials applied to top and bottom, ϕ1=ϕ0, ϕ2=ϕ0.
The results are non-dimensionalized with Y 0=10.3 GPa, and d0=374.0×10−12 CN−1 for panel (a) and d0=100.0×10−12 CN−1 for panel (b), as
pressure load case:
electric potential load case:
(a) Single-layer piezoelectric panel
Figure 3 presents the response of a thick PZT-5A panel (a) (S=5) with simply supported ends that are electrically grounded (ϕ=0) (i.e. under CC condition). The longitudinal variations of deflection , and stresses and , at z-locations where they are large, are plotted for both single-term (n=1) and multi-term (n>1) trial functions. For this boundary condition, the exact analytical solution of the three-dimensional piezoelasticity equations was presented by Kapuria et al.  using a Fourier series expansion in the span direction. The present results are compared with this exact solution in figure 3. Since the results for uniform loading were not presented by Kapuria et al. , these results have been obtained here using the computer program developed by them for their numerical study. Converged results are obtained using 49 non-zero terms in the Fourier series solution, whereas the EKM solution converges with only one term in the expansion in just two iterations. The EKM solution is thus much more computationally efficient than the Fourier series solution, apart from the fact that the latter is valid only for simply supported panels with CC edge conditions, whereas the present EKM is applicable to arbitrary mechanical as well as electrical boundary conditions. The figure indicates that the EKM solution is in excellent agreement with three-dimensional exact solution for the displacements and stresses.
Figure 4 shows the through-thickness distributions of the electric potential and stresses and at various locations along the span of a clamped simply supported (CS) PZT-5A panel. Converged results for this case are obtained with n=3 and two iterations. These results are compared with an FE solution obtained using ABAQUS software. The eight-node quadratic serendipity plane strain piezoelectric (CPE8RE) element with reduced integration is used for the FE model. Converged results are achieved with a mesh size of 100 (length) ×40 (thickness). An excellent match is observed between the two results at all ξ-locations, for the stresses as well as the induced electric potential.
The effect of electric boundary conditions is studied in figure 5, by plotting results for both CC and OC conditions at the top for the piezoelectric panel with CS boundary conditions. It is revealed that for the single-layer piezoelectric plate, the electric boundary condition at the top (or bottom) surface does not have any significant effect on the displacements and stresses. In figure 5, the results obtained without considering the electromechanical constants (dij=0) are also plotted. It is observed that the electromechanical coupling has a significant effect on the displacements over the entire span, but the shear stress gets affected only near the clamped edge, where has a large value (in the coupled case).
(b) Hybrid sandwich panel under pressure loading
The results for a thick (S=5) simply supported hybrid sandwich panel under pressure loading with CC conditions at all boundaries are compared in figure 6 with the three-dimensional exact piezoelectric solution . The results are obtained with up to three terms in the trial functions (n=1,2 and 3). In each case, convergence was achieved in two iterations. It is seen that a two-term solution with two iterations is sufficient to accurately predict the response of the simply supported hybrid panel under the pressure load case, including for the electric variables. The longitudinal variations of the sensory potential induced at the top surface under OC conditions and the electric displacement are depicted in figure 7. Once again, the two-term solution is an excellent match with the three-dimensional exact solution, even as the single-term solution is unable to predict it accurately.
The longitudinal variations of displacements and , normal stress in the piezoelectric layer and in the elastic substrate , shear stress and induced electric potential of the panel under CS boundary conditions are shown in figure 8. The results are presented for n=1,2 and 3 with two iterations in each case, which yielded converged results, even though the initial trial functions did not satisfy the given boundary conditions. For all response entities, the EKM results converge with two terms (n=2). For comparison, a detailed FE analysis is performed in ABAQUS, using the eight-node plane strain piezoelectric (CPE8RE) and elastic (CPE8R) elements for the piezoelectric and elastic layers with a mesh size of 100 (length) ×36 (thickness). It is observed that the single-term solution yields very inaccurate results for the displacement and induced sensory potential , and also for the stresses , and near the clamped support, for the thick hybrid laminated plate. The error in is larger than that in . Overall, the error in the single-term solution for the hybrid plate is higher than that observed by Kapuria & Kumari  for elastic sandwich plates. The two-term solution results are in excellent agreement with the FE solution for all variables, even in the vicinity of the clamped support. There is, however, a mismatch with the FE solution for stresses and at the clamped support. But, this mismatch is because the FE solution does not satisfy the conditions of applied normal and (zero) shear tractions at the top and bottom surfaces and the conditions of continuity of transverse stresses at the layer interfaces at the clamped edge. This is seen from the through-thickness distributions presented in figure 9. However, the FE predictions for the transverse stresses satisfy the boundary and interface conditions, and match well with the EKM solution at some distance away from the clamped support.
(c) Hybrid sandwich panel under potential loading
Results for a simply supported hybrid sandwich panel under electric potential loading are presented in figure 10. The edges of all the panels studied for the potential loading case are considered to be under CC conditions. In this case, there is a sharp increase of from zero at the support to near maximum within a distance of a/25 from the support. While a two-term solution is generally good for other variables, more terms (n=4) are needed to predict the sharp variation in near the edge, for this load case. Also, contrary to the pressure load case for which a single-term solution gives fairly accurate results for the simply supported panels, its results are way out for the potential load case, and hence are not presented here for comparison. It is also found that the solution for any value of n converges in the first iteration itself, for the potential load case.
The longitudinal variations of displacements, stresses and electric displacement for the clamped-free (CF) boundary condition are presented in figure 11. Similar to the simply supported case, the EKM solution converges with n=4 and matches well with the FE solution, except for the stresses and electric displacement at the clamped and free edges. This is contrary to the pressure loading case for which the two results differ only at the clamped edge. It is clear from the figure that the FE solution does not satisfy the conditions of σx=0 and τzx=0 at the free edge in the potential load case, and hence the mismatch.
Results for the CS boundary condition are presented in figure 12 for the thick (S=5) panel. In this figure, stresses are presented at different locations across the thickness in the piezoelectric and elastic layers. It is observed that the nature of variations of and the boundary-layer effect significantly depend on the location across the thickness. The through-thickness distributions of and for the panel are plotted in figure 13. For this case too, the FE fails to satisfy the interfacial continuity condition of at the support (x=0) and hence cannot predict even accurately at this location. The two distributions, however, match away from the support.
An accurate three-dimensional piezoelasticity solution is presented for the coupled electromechanical response of piezoelectric laminated plates in cylindrical bending, using the powerful multi-term EKM. It represents the first application of the EKM to a coupled multi-field problem. The convergence and accuracy of the solution has been established by comparison with available exact three-dimensional piezoelasticity solutions and detailed FE analyses. The numerical study reveals the following.
— The method yields accurate response of piezoelectric sensors as well as piezolaminated plates for both pressure and electric potential loadings.
— For pressure loading, a two-term solution is good enough to yield accurate results for all variables, including the sharp variations of the induced electric potential ϕ and electric displacement Dz observed near a clamped support.
— Unlike the pressure loading for which the boundary layer is observed only at the clamped support, for the electric potential loading, very sharp variations in stresses σx and τzx occur at the simply supported and free edges.
— A higher number of terms (up to n=4) is usually required to predict these sharp variations of stresses and electric variables in the potential load case.
The present study will facilitate the development of accurate semi-analytical solutions of many other unresolved problems in three-dimensional piezoelasticity such as the free-edge stresses in hybrid laminates under bending, tension and twisting.
- Received September 26, 2012.
- Accepted December 20, 2012.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.