Hybrid asymptotic–direct approach to finite deformations of electromechanically coupled piezoelectric shells

We present a novel multistage hybrid asymptotic–direct approach to the modeling of the nonlinear behavior of thin shells with piezoelectric patches or layers, which is formulated in a holistic form for the first time in this paper. The key points of the approach are as follows: (1) the asymptotic reduction in the three-dimensional linear theory of piezoelasticity for a thin plate; (2) a direct approach to geometrically nonlinear piezoelectric shells as material surfaces, which is justified and completed by demanding the mathematical equivalence of its linearized form with the asymptotic solution for a plate; and (3) the numerical analysis by means of an FE scheme based on the developed model of the reduced electromechanically coupled continuum. Our approach is illustrated by examples of static and steady-state analysis and verified with three-dimensional solutions computed with the commercially available FE code ABAQUS as well as by comparison with results reported in the literature.


Introduction
Multifunctional materials with sensing and actuation capabilities and their integration into systems of structural mechanics are the basis for the development of so-called intelligent or smart structures, which are used in mechanical, aerospace or civil engineering; for reviews and applications, see, e.g., Crawley [1], Tani et al. [2] and Liu et al. [3].
The proper design of smart structures is strongly related to an accurate, but also numerically efficient modeling of the system. In many cases, these structures can be assumed as thin plates or shells with integrated multifunctional materials. In such problems, a three-dimensional modeling is often inefficient and mostly not necessary. Here, the theories of structural mechanics come into play, which must be extended with respect to the effect of the structurally integrated multifunctional materials. A prominent example for such materials with sensing and actuation capabilities is piezoelectric materials, for which the electromechanical coupling must be accounted for. Various approaches are reported in the literature, including equivalent single layer theories (Krommer [4], Batra and Vidoli [5]) or layer-wise and hybrid ones (Carrera and Boscolo [6], Moleiro et al. [7]). In these theories, a priori assumptions concerning the mechanical displacements and the electrical field in certain directions are made to reduce the three-dimensional continuum to a lower-order one. With respect to nonlinear formulations for piezoelectric plates and shells, we refer to the rich literature, e.g., Zheng et al. [8], Tan and Vu-Quoc [9], Klinkel and Wagner [10], Marinković et al. [11] and Lentzen et al. [12]. In this paper, we develop a novel hybrid approach to the modeling of the nonlinear behavior of piezoelectric shells, which is here presented in a holistic form for the first time; the approach is illustrated by two numerical examples.
The key point of the novel multistage hybrid asymptotic-direct approach for thin piezoelectric shells is as follows. The geometrically nonlinear shell theory, obtained with a direct approach to shells as material surfaces, is justified and completed by demanding the mathematical equivalence of its linearized version to the asymptotically reduced equations of the three-dimensional linear theory of piezoelectricity. This provides the constitutive relations of the dimensionally reduced continuum, allows for restoring the stressed state through the thickness with the solution of the reduced problem and, finally, leads to consistent solutions of challenging problems of the nonlinear behavior of thin shells with complex properties.
The first component of the hybrid approach is the procedure of asymptotic splitting in the form suggested by Eliseev [13], which we apply to the three-dimensional linear problem of deformation of piezoelectric plates. The approach has been systematically used for the dimensional reduction in the theories of curved and twisted thin rods with an inhomogeneous cross section [14], of thin-walled rods of open profile [15], of curved composite strips [16] and of thin piezoelectric plates with an inhomogeneous material structure through the thickness [17]. The main idea is that the leading-order terms of the series expansion of the solution are determined from the conditions of solvability for the minor terms.
The second component is the direct approach, based on the principles of Lagrangian mechanics extended to the electromechanically coupled continuum: The system of equations of a dimensionally reduced theory follows from an extended equation of virtual work. The procedure is compact and formal, as soon as the degrees of freedom of particles of the reduced continuum are determined and the constraint conditions are formulated.
The third component of the approach is the numerical analysis, which is based on the developed model of the reduced continuum. The kinematics of a material surface is approximated according to the theory, and the corresponding strain measures and the expression for the virtual work are used. We will compare the outcome of our numerical simulations against results computed with three-dimensional electromechanically coupled elements using ABAQUS software, 1 to semi-analytical solutions using the linearized version of the proposed shell theory and to results from the literature. The developed coupled shell theory is verified in the geometrically linear and the finite deformation regime by convergence studies.
The paper is structured as follows. After this introductory section, the governing equations of the threedimensional linear theory of piezoelectricity are shortly reviewed in Sect. 2. The application of the procedure of asymptotic splitting-the first component-in Sect. 3 results in a linear two-dimensional electromechanically coupled theory for thin piezoelectric plates. Section 4 presents the direct approach-the second componentfor geometrically nonlinear thin piezoelectric shells as material surfaces with mechanical and electrical degrees of freedom. The FE implementation-the third component-is described in Sect. 5. The paper is completed with two example problems in Sects. 6 and 7.

Piezoelastic continuum: electromechanically coupled theory
A coupled system of equations for mechanical and electrical field entities needs to be considered, when the material of the structure exhibits piezoelectric effects; see Nowacki [18]. Using the index '3' where appropriate to distinguish the three-dimensional fields from their counterparts in the two-dimensional plate or shell theory, we write the equations of balance in the volume of the body V and the corresponding condition at the boundary = ∂ V with the outer unit normal n: Here, τ 3 is the symmetric stress tensor, f 3 is the external body forces (which may include inertial terms in dynamics), p is surface tractions, and the differential (Hamilton) operator in three dimensions is denoted ∇ 3 .
In the linear theory, the strain tensor ε 3 is the symmetric part of the gradient of displacements u 3 : Strains, for which a suitable field of displacements exists, fulfill the condition of compatibility: Either displacements u 3 or tractions p are to be prescribed at as mechanical boundary conditions. The constitutive relations in Voigt's linearized theory of piezoelasticity involve also the electric displacement vector D and the electric field vector E: The elastic properties of the material are defined by the fourth-rank tensor 4 C, is the tensor of dielectric constants, and the third-rank tensor of piezoelectric constants 3 e is nonzero only in materials exhibiting the piezoelectric effect.
The electric field is defined as the negative gradient of the electric potential ϕ 3 , such that ∇ 3 × E = 0 holds, and the divergence of the electric displacements vanishes in the volume: Electric voltage, prescribed on a part of the surface, provides a boundary condition for ϕ 3 . The complementary boundary condition is written in terms of the free charge density σ per unit surface area: The constitutive relations (4) may be rewritten using the enthalpy function: Hamilton's principle for piezoelectric continua is formulated with the following counterpart to the total potential energy: This is a functional over the fields of displacements and electric potential: we seek the static equilibrium as a stationary point of (9): The variations δu 3 and δϕ 3 are independent, and the equations of balance (1), (6) as well as all kinds of boundary conditions follow from (11). In a purely mechanical problem, δ H Σ = 0 is the principle of virtual work-also denoted as Gibb's principle for a conservative problem; see, e.g., the textbook by Ziegler [19].

Asymptotic solution for a thin plate
Within this section, the first component of the hybrid approach-the procedure of asymptotic splitting-is discussed. The literature concerning the solution of equations for the three-dimensional model of a plate with the help of various asymptotic methods is exhaustive. The three-dimensional solution is typically sought as a series expansion with respect to a small parameter, which is related to the thickness of the plate; e.g., two-dimensional equations for the leading-order terms of the series expansion for a homogeneous plate were obtained in [20,21]. For an advanced study of piezoelectric plates with a periodic material non-homogeneity in all three dimensions, see [22]. Alternatively, the asymptotics can be based on variational principles as developed in [23]. For the present paper, the works [24,25] are relevant; the equations for the components of stresses and displacements are written for a specific material structure in a non-dimensional form, and the analysis proceeds by assigning particular orders of smallness to different components of the stress tensor a priori. The two-dimensional equations are obtained from the conditions of solvability for the minor terms of the series expansion of the solution.
In this study, we rely on the advantages of the procedure of asymptotic splitting in the form suggested by Eliseev [13,14]. All groups of equations of the theory of piezoelasticity are treated independently with the help of the conditions of compatibility, and no assumptions concerning the orders of smallness of different components of the stress tensor are involved. Leading terms in the solution follow from the conditions of solvability for the minor terms. The presented theory has been numerically validated by demonstrating the convergence to exact solutions for thin plates and curved strips in [16,17,26]. Below we summarize the notions and main steps of the asymptotic analysis of piezoelectric plates, first published in [17].

Small parameter and asymptotic expansions
We apply Eqs. (1)-(7) to a thin plate. Choosing the first two cartesian axes x, y in the plane of the plate and the third one z in the thickness direction, we write the position vector of a point of the structure: We have introduced unit basis vectors i, j , k in the directions of Cartesian axes, is the two-dimensional domain, and z ± are the coordinates of the upper and lower free surfaces of the plate. Its thinness is indicated by the formal small parameter λ; hence, the magnitudes of x and z are of the same order. We will analyze the asymptotic behavior of the solution of the problem as λ → 0. The advantage of the dimensionless formal small parameter is that it can be set equal to 1 as soon as the dominant terms in the solution are determined. Now, ∇ 3 r = I is the identity tensor, and the corresponding form of Hamilton's operator will be ∇ is the differential operator with respect to the in-plane position vector x. We seek solutions, which vary in the plane much slower than over the thickness (z is a "fast" variable), and therefore, the derivatives with respect to the in-plane coordinates x acquire a corresponding order of smallness in (13). Unknown fields are sought as asymptotic expansions in terms of the small parameter. Stresses grow when z + − z − → 0, which means that negative powers shall be included: We are interested in finding those terms in the solution, which dominate as the plate is getting thinner and λ → 0. In doing so, we generally rely on the existence of the solution of the original problem for any positive thickness of the plate. The leading power λ −2 in (14) is known from the results of the asymptotic analysis of a homogeneous elastic plate [27].
Owing to the structure of our problem, it is convenient to separate the in-plane and out-of-plane parts of vectors and tensors. The in-plane part will be denoted with a subscript ⊥: in which γ and τ are the out-of-plane shear strain and stress vectors.

Asymptotic splitting in the equations of piezoelasticity
We use (13), (15) in (1) and rewrite the equations of equilibrium with the small parameter; see [17,26] for details. For the upper and lower surfaces free from tractions, which means we substitute the expansion (14) and begin the asymptotic procedure. At the first step, we balance the principal terms in the equations and boundary conditions, which have the order λ −2 . From the simple equalities, it follows that the shear and the transverse components vanish in the leading-order term of the stress tensor: The principal term of the in-plane stress tensor 0 τ ⊥ is yet unknown and shall be determined as a final result. At the second step, we proceed to the terms of the order λ −1 , taking (17) into account. We conclude that 1 σ z = 0, and The equation of balance of the in-plane force factor (stress resultant) τ(x) appears as a result of the condition of solvability for the shear stress 1 τ. A complete system of equations requires the third step of the asymptotic procedure, at which the external force factors come into play. Collecting the terms of the order λ 0 , substituting the results of the previous steps and writing the conditions of solvability for the minor terms 2 τ, we arrive at The equation of balance for the vector of transverse forces Q is the same as in the classical plate theory. Eliminating Q, we arrive at the second classical equation of equilibrium for the tensor of internal moments (stress couples) μ(x): Not only the known balance equations but rather the asymptotic structure of the components of the stress tensor is the outcome of the procedure. The boundary conditions at the side edges ∂ result from the asymptotics of the edge layer; see [17].
Being linearly related to stresses via (4), the remaining fields have the same asymptotic behavior with respect to the small parameter: the field of electric potential has the same asymptotic order as there is no λ in front of k∂ z in the expression of ∇ 3 (13). Now we use the asymptotic expansion of the strain in the condition of compatibility (3). Gathering the principal term of the in-plane component of the resulting equation, we conclude that the plane part of the principal term of the strain tensor is linearly distributed through the thickness: The plane tensors κ(x) and ε(x) are functions of the in-plane coordinates. An asymptotic analysis of displacements [17] shows that κ is the negative linearized curvature and ε is the in-plane strain of the plate: with u and w being the leading terms in u ⊥ and u z . This is the modern interpretation of Kirchhoff's kinematic hypotheses. Now, we consider the practically important case of a piezoelectric plate with two electrodes, one at the upper and one at the lower surface. The material parameters can be arbitrarily distributed through the plate thickness. The voltage between the electrodes is the potential difference: Using (6) and the differential operator (13), we find This corresponds to the assumption which is traditionally made in deriving the theory of piezoelectric plates with the method of hypotheses; see [4,28,29].
Writing the second constitutive relation in (4) for the principal terms, projecting on the thickness direction and using (5) and (13), we arrive at The principal terms in the solution shall be influenced by both the mechanical and the electrical loadings. The electric voltage λ −2 ϕ between the two electrodes and the total charge λ −2 Σ on the upper electrode need to be two orders of smallness larger than the volumetric mechanical force f z to have comparable effects on the behavior of the structure, and Here e is a two-dimensional domain, occupied by the electrodes. Now we proceed to the principal terms in the first of the constitutive relations (4). With (26), we write: As k · 0 τ = 0, we find 0 ε z and 0 γ as linear functions of 0 ε ⊥ = −κz + ε and D z0 from (28). We substitute them into the relation for the voltage, which follows by integration of (26) according to (27): This results in a constitutive relation for the distributed charge on the electrodes: Integrating further the in-plane part of (28) over the thickness, we finally obtain the internal plate force factors τ, μ as linear functions of κ, ε and ϕ: In these constitutive relations, the fourth-rank tensors 4 A, 4 B and 4 D determine the elastic properties of the plate at constant voltage. We also like to point out that this result accounts for electromechanical coupling, as one can see from the definition of the effective material parameters in (28) and (29). The coefficientsp,m are equivalent to p, m from (30) for orthotropic materials, a case to which the present paper is restricted to. A general proof of this fact for arbitrary material anisotropy, which is important for the use of these constitutive relations in the framework of a direct approach in Sect. 4 with a function of enthalpy per unit area of the plate, is yet to be found.
The above asymptotic procedure has resulted in the governing equations of a two-dimensional continuum, namely a piezoelastic plate. The equations are equivalent to the ones derived in [4] by making a priori assumptions concerning the distribution of the displacements and the electric field through the thickness. Hence, the asymptotic procedure justifies our previous results. We also note that the asymptotic procedure has been conducted for a plate with electrodes at both sides, for which, however, the material parameters are arbitrary functions of the thickness coordinate z. One can extend this procedure for plates with more than one piezoelectric layer and hence multiple sets of electrodes, which can be independently used for actuation and sensing. For the definition of material parameters in such a case, we refer to [4] as well. The asymptotic accuracy and simplicity in applications are the benefits of the presented approach in comparison with treating the fields in a laminate structure layer-wise [6,30].

Solution of the coupled structural problem: workflow
Two electrical unknowns exist for each pair of electrodes at the opposite sides of the plate: the total charge Σ or the voltage ϕ; see (27). Depending on the impedance of the electric circuit, which connects the electrodes, a relation between Σ and ϕ complements the structural equations. Two typical cases are as follows: -Actuation. The voltage ϕ is prescribed, and the corresponding terms in (31) may be interpreted as generalized forces, which result in the deformation of the structure. A proper distribution (design) of these forces can be used to assign a desired deformation to the plate (see [31,32] for applications to beam and frame structures). -Sensing. The measured voltage ϕ is interpreted by an observer in terms of structural entities: displacements, vibration amplitudes, etc. In the simulation, ϕ is an additional unknown, and this potential difference is the same for the whole pair of opposing electrodes. The system of equations is completed by a given value for the total charge Σ, which remains constant as the external electric circuit is open. A proper design of the sensor can be used to measure various kinematic entities of interest; see [32][33][34][35].
In the present paper, we solve a structural shell problem either according to the above linear equations, or using the nonlinear shell model from Sect. 4 and applying the finite element scheme from Sect. 5. Then, knowing the strains and the voltages in terms of the structural theory allows restoring the leading-order terms of the three-dimensional mechanical and electrical fields through the thickness of the shell.

Direct approach to finite deformations of a piezoelectric shell
We proceed with the second component of the hybrid approach-the direct approach to thin piezoelectric shells as material surfaces. As the curvature and geometrically nonlinear effects essentially couple the governing equations for the membrane and bending deformations of a shell, methods of Lagrangian mechanics are applied to enable the extension of the asymptotically validated theory of linear plates to the case of nonlinear shells with relative ease.

Finite deformations of piezoelastic shells as material surfaces
A classical Kirchhoff-Love shell is a material surface with five mechanical degrees of freedom of particles: three translations δr and two rotations δn; the variation of the unit normal to the surface n lies in the tangent plane. Each particle can be interpreted as a material unit normal (a needle), such that no virtual work is produced on its "drilling" rotations. This determines the equation of virtual work and the resulting theory in the purely elastic case; see [36,37]. The complete system of equations of the electromechanically coupled problem at hand is again provided by a variational principle after assigning an electrical variable ϕ, which has the meaning of the potential difference over the thickness, to each particle of the material surface. We extend the three-dimensional condition of stationarity of the total enthalpy (9) to a form, which is similar to the purely mechanical equation of virtual work: The term σ δϕ has "migrated" into the integral over the domain together with the free charge σ on the electrodes. The area change from the reference to the actual configuration appears from d = J d • , in which the small circle denotes variables in the undeformed reference state. The external forces inside the domain q and at the boundary P are counted per unit actual area and length of the deformed contour ∂ . The enthalpy per unit area in the reference configuration H will be determined in the course of the analysis, and the moment loadings are ignored for simplicity. At a position of static equilibrium, the functional H Σ has a stationary point.
Three parts of the domain need to be considered.
-At the part of the plate, which is free from piezoelectric patches/layers, H is the strain energy, and ϕ is not present. -If a patch/layer works as an actuator, then some voltage is prescribed on the electrodes. This shall be considered as a kinematic constraint for ϕ. -The electric circuit remains open at the electrodes of the patches/layers, which act as sensors, and ϕ is an additional unknown.
All electrodes are equipotential areas: ϕ = ϕ i , δϕ = δϕ i in i , and we rewrite the variational principle introducing total charges Σ i at the electrodes: i are those domains of the shell, at which sensors/actuators are present, and which are electrically independent.
In the following, we shall demonstrate that the equations of balance result as a consequence of (33). Again, τ and μ are the in-plane tensors of forces and moments in the shell, the in-plane vector Q determines the shear force, and a = ∇r = I − nn, b = −∇n (35) are the first and the second metric tensors of the surface. The differential operator ∇ can be determined by the expression of a differential of a (not necessarily scalar) field u, defined in the two-dimensional domain: du = dr · ∇u; dr = dr · ∇r = dr · a.
For a discussion of the basic notions of differential geometry of a surface, we refer to [38,39]. For a presentation with coordinates on the surface, see [13,36] as well as the finite element implementation in Sect. 5.2.
It is known (see "rigidity theorem for surfaces" in [38]) that a surface undergoes a rigid body motion if (and only if) the components of both metric tensors in a local materially fixed basis ("convected components") remain constant. In the invariant form these conditions read with the geometrically nonlinear strain tensors see (56) for a representation using components. The deformation gradient F relates an actual line element in the deformed configuration to its pre-image in the reference one: G plays the role of the "inverse" of F. In the absence of deformation, both E and K vanish owing to the metric tensors of the reference configuration • a and • b in (38). Now, for classical shells the unit normal n remains orthogonal to the surface in the course of deformation, which imposes the constraint δn + ∇δr · n = 0. (40) Indeed, any differential dr is orthogonal to n. Considering δ(dr · n) = 0 and using (36) for dδr, we arrive at (40).
Like in the purely mechanical case [26,36], we consider the variational formulation (33) for a shell, undergoing a rigid body motion with preserved electric potentials. The enthalpy shall remain constant under the following constraints: We introduce Lagrange multipliers for the above constraints, transformed with G, and rewrite (33): In addition to the stress resultants, which already appeared in (34), we have introduced a new Lagrange multiplier, namely the total charge at an electrodeΣ i (which certainly makes a difference only for open-circuit conditions with δϕ i = 0). Symmetry of the left-hand sides of the first and the second constraints in (41) implies that both τ and μ are symmetric.
After mathematical transformations and the use of the divergence theorem on the surface, the variational equation yields Here · · · 1 and · · · 2 are the left-hand sides of equilibrium Eq. (34), which are thus proved because the Lagrange multiplier Q allows treating the variations δr and δn as independent inside the domain. The contour integral leads to the consistent static and kinematic boundary conditions, and at the open-circuited piezoelectric patches/layers, the total electric charge needs to be equal to a prescribed value: Further analysis leads to the general form of the constitutive relations. Considering again (33) for arbitrary variations of unknown fields, we find The prescribed charges are equal to the actual ones (denoted previously with a tilde), and as δϕ = δϕ i at open-circuited electrodes and vanishes elsewhere. In a purely mechanical case, the equality of integrals means the equality of integrated expressions. This conclusion is commonly based on a consideration that (45) remains valid for an arbitrary sub-domain inside , which may be justified using the principles of Lagrangian mechanics [26]. In the present electromechanically coupled case, a non-controversial theory follows with a similar local variational relation, which is in correspondence with (45): the arbitrariness of variations at the right-hand side of a variational relation determines the list of arguments of a state function H . The constitutive relations follow from (47):

Enthalpy function for a shell
Solving particular problems, we need a function H (E, K, ϕ), which can be assumed as a quadratic form, because the local strains remain small even at large overall displacements and rotations. Moreover, the constitutive relations of the nonlinear shell theory (48) applied to the problem of small deformation of a plate shall be equivalent to the results of the asymptotic splitting of the three-dimensional electromechanically coupled problem in Sect. 3. Comparing with the constitutive relations of the linear plate theory (30), (31) for ε = E and κ = K, we arrive at the following expression of the enthalpy: The fourth-rank tensors 4 A, 4 B and 4 C comprise the elastic properties of the cross section when the voltage ϕ = const. The coupling between the voltage and the mechanical deformations is described by the second-rank tensors p, m and the capacity c. As discussed in Sect. 3.3, the solution of a structural problem allows restoring three-dimensional fields in the volume of the body of the shell because the local effect of the curvature is expected to be asymptotically negligible.
For a shell made of n rigidly/perfectly bonded transversally isotropic layers, the enthalpy can be written as Each piezoelectric layer is assumed to be electrically independent with the potential ϕ j ; non-piezoelectric layers are accounted for with p j = 0 and m j = 0. Explicit expressions of the coefficients and numerical values for such a layered shell are provided in Appendix A.

Finite element analysis of electromechanically coupled piezoelectric shells
Commercial finite element codes such as ABAQUS or ANSYS still do not include plate or shell finite elements that are capable of solving electromechanically coupled problems. Concerning the literature on piezoelectric thin-walled structures, the degenerated shell approach prevails; see, e.g., [40,41]. Here, we extend a numerical scheme, presented originally in [42], to the electromechanically coupled variational problem (33). This is the third and final component of the present hybrid approach.
The above modern version of the classical Kirchhoff-Love theory of shells requires in general C 1 continuity in the approximation of the deformed surface, which can be obtained in various ways [43][44][45]. We achieve this goal using a four-node finite element with the approximation scheme, suggested for linear plates in [46]: 16 The element has thus 48 mechanical degrees of freedom and belongs to the class of the so-called absolute nodal coordinate formulation (ANCF, [47,48]). A constant mass matrix and smooth approximation make it as well attractive for the transient analysis and non-material formulations for axially moving structures [49]. Despite inherent restrictions concerning the topology of the mesh and connections between shell segments, the present finite element has a relatively broad spectrum of potential applications with respect to both research and development.

Element kinematics and shape functions
The shell is modeled as a continuum of material normal vectors, and its configuration is defined by the field r(q 1 , q 2 ) ≡ r(q α ); −1 ≤ q α ≤ 1 are two material coordinates on a finite element. The vector of unit normal results from the natural basis: while the differential operator on the surface features the co-basis: Rewriting the expressions for the strain measures (38) with components (56), we see that the enthalpy H in (49), (50) will be integrable as long as the smoothness condition is fulfilled. We denote the four nodes of the finite element as i, j, k and l and write the position vector as The approximation features 12 nodal degrees of freedom: position vector r m , its derivatives with respect to the local coordinates on the element r m 1 and r m 2 , and the mixed second-order derivative r m 12 ; see Fig. 1. The conditions of smooth coupling with the approximation on the neighboring elements lead to a unique set of the 16 bi-cubic shape functions S m,n (q α ). The element itself is isoparametric: The reference geometry • r is also approximated by means of (53) and is thus C 1 continuous. The validity of the presented approximation requires that the coordinate lines q α = const are continuous across the element boundaries. This poses certain limitations on the topology of the mesh, which restricts the range of potential applications. Releasing this requirement is a non-trivial task, which will not be discussed here.

Elastic forces and stiffness of the element
Solving the variational problem (33), we assemble the total enthalpy integrating over the reference configurations of all finite elements with H ext being the potential of external forces q and P. Unknowns are the nodal degrees of freedom and electric voltages ϕ k at open-circuit regions (sensors). This latter fact introduces additional coupling into the system for all elements, which belong to the same i . This coupling results in a non-local formulation over the domain i , as it has been discussed in detail in [4] for a linear plate. As in (50) c > 0 and the enthalpy H is not positive definite, we cannot speak about the minimality problem in the presence of electrical unknowns, but experience shows that small negative eigenvalues of the stiffness matrix do not influence the convergence of the quasi-Newton scheme for solving the equation in (54). Dividing the domain into finite elements, we keep in mind that the function H Σ is integrable due to the continuity of the approximation. Numerical results in Sect. 6 were obtained by using the Gaussian quadrature rule with 3 × 3 integration points for the integral over the element surface which also holds for K = K αβ The components of the metric tensors in the reference configuration • a αβ , • b αβ shall be precomputed in the integration points of the finite elements: They fully determine the undeformed geometry of the shell, and the elementary surface is We employ the Newton method for seeking the stationary points of H Σ , for which we need the derivatives to be computable for any configuration of the shell. There are N = 48 mechanical degrees of freedom e p for finite elements without piezoelectric properties or with prescribed electric voltage, and N = 49 for elements in open-circuit regions, in which the electrical unknown ϕ is shared among multiple elements. The global force vector F and stiffness matrix K result from assembling (summation) of F el and K el over the integration points, where the derivatives of the distributed enthalpy H (E, K, ϕ) are computed efficiently by a chain rule.

Boundary conditions
If an edge is free from kinematic constraints, then the external force factors acting on that edge need to be accounted for. In static problems, it is common to deal with conservative loads, which allows to speak about the potential energy of external force factors at the boundary. The most simple case of a conservative edge load is a force, which is distributed per unit length of the edge in the reference configuration. This means that the force vector P changes with the extension or contraction of the edge. The contribution of the potential energy of this force to H ext in (54) is easy to compute by integrating over the edges of the elements at the boundary. At a simply supported edge, the particles are fixed by appropriate penalty terms for the nodal positions r m and derivatives r m α (α corresponds to the direction along the edge). If the edge is clamped, then the direction of the normal vector n needs to be additionally constrained. For a straight edge, n = • n = const, and the constraint will be fulfilled exactly, if we demand • n · r m β = 0 and • n · r m 12 = 0, in which β corresponds to the direction pointing outwards of the domain. The same conditions can be applied at curved edges with satisfactory results.

Extension to dynamics
The equations of motion of a shell in the form of Lagrange equations of the second kind (in which the enthalpy replaces the potential energy) require the kinetic energy of the element T el . The effect of rotary inertia of a through-the-thickness element does not need to be accounted for classical shells, and we write here, (ρh) • is the material density per unit area in the reference configuration, and the dot means a time derivative.
Rewriting the approximation (53) as {r} = e T S, in which S is the 3 × 48 matrix of shape functions and e is a column matrix of degrees of freedom e p , we obtain The time integration is simple and straightforward due to the constant mass matrix M el , which is typical for the absolute nodal coordinate formulation [47,48].

Example 1: Cylindrical panel with piezoelectric patches
Some parts of the presented approach have been earlier demonstrated in applications to simple piezoelectric plates [17,26,28,50] and to shell problems [26,51]. In this section, we present a thorough study of the capabilities of the method on the practically relevant example of the static and dynamic behavior of a cylindrical panel; see Fig. 2.
The panel is clamped at one end and equipped with three piezoelectric patches, which are perfectly bonded at the outer surface of the panel (with larger radius). The structure of the panel is made of aluminum, and an orthotropic piezoelectric material PZT-5A is used for the patches.
We choose the middle of the substrate (structural layer) as a reference surface z = 0 and the thickness of both layers equals h = 10 −3 m, such that the substrate is at −h/2 ≤ z ≤ h/2 and the piezoelectric patch spans h/2 ≤ z ≤ 3h/2; see Fig. 3.
The structural properties of the through-the-thickness element z − ≤ z ≤ z + with z − = −h/2 and z + = 3h/2, which appear as coefficients in the quadratic form of the enthalpy (50), were computed according to the procedure, described in detail in Appendix A. The length of the panel is L = 0.4 m, and the radius of curvature of the middle surface of the structural layer is R = 0.08 m. The size of the piezoelectric patches in the length direction is 0.05 m, and the distance between them and the clamped edge is 0.075 m. The angular size of each patch in circumferential direction is π/6, such that the angular sizes of the two gaps between them are π/4. The panel is under the action of the field of gravity with the magnitude g = χg max , in which 0 ≤ χ ≤ 1 is a variable load factor and g max = 5000 ms −2 is the maximum loading. As it is indicated by an arrow in Fig. 2, the gravity force acts downwards along the line connecting two free corners of the undeformed panel.
In the following, we consider solutions, which were obtained with the above finite element scheme based on the hybrid approach to shells as material surfaces. The analysis is performed using the in-house software code ShellFE3, developed by the authors. The results are compared to the solution for an equivalent threedimensional finite element model, which was obtained using ABAQUS software. We used the 20-nodal brick elements with electromechanical coupling C3D20RE, and a regular mesh featuring 3 elements in the thickness direction for each layer, 33 elements per patch in the length direction and 28 elements per patch in the circumferential direction (136575 finite elements in total). This mesh provided solutions, which were converged to a sufficient degree for the comparisons below.

Static analysis: sensing
Considering open electric circuits between the electrodes and zero total charges, we computed the static deformation of the structure and the electric voltages at the patches, which are thus acting as sensors, whose signals may be interpreted in terms of structural entities. The deformed model with 4 × 4 finite elements per patch for the final value of the load factor χ = 1 is presented in Fig. 4 from two points of view.
We notice a local bulge at the upper edge near the clamping, where the shell is softer than under the patch. The computed histories of the evolution of all three sensor signals ϕ k are shown in Fig. 5.  Even the coarse model with just 4 finite elements per patch produces accurate results, and refining the mesh we converge to a solution, which is close to the results of the three-dimensional analysis.
Owing to the symmetry of the problem, it is easy to conclude that the signal of the middle patch ϕ 2 (χ ) would vanish in a linear formulation, which is shown in Fig. 5 for small χ . Geometric nonlinearity results in ϕ 2 = 0 later on. One notices also a visible nonlinearity in the other two signals starting at χ ≈ 0.7, which corresponds to the beginning of the formation of the bulge. The local additional stiffness of the perfectly bonded patches slightly shifts the location of the bulge and increases the critical load, but the supercritical behavior remains nearly the same as in the structure with no patches. The location of the bulge also determines the fact that the nonlinearity is most prominent in the signal ϕ 1 , which can be used for drawing conclusions concerning the location of defects in the structure [34,35].

Evolution of eigenfrequencies
Having the stiffness matrix K and the mass matrix M of the finite element model at hand, we can solve a generalized eigenvalue problem and find eigenfrequencies ω k of lower vibration modes, which are roots of the equation For an undeformed structure at χ = 0 and with all three patches in open-circuit conditions, the results of computations with different meshes and in ABAQUS are compared in Table 1. The small difference between the converged eigenfrequencies in the shell and in the three-dimensional solutions is nevertheless higher than the effect of the electric circuit at the patches, which can be observed by comparison with the ABAQUS solution with prescribed voltages ϕ 1,2,3 = 0.
We also computed the variation of eigenfrequencies with respect to the above static solution. Increasing the load factor χ typically reduces the stiffness along with the eigenfrequencies due to the corresponding pre-stresses and pre-deformation, which can be observed in Fig. 6. Again, a qualitative change in the behavior is seen as χ > 0.7.

Steady-state vibrations: sensing and actuation
We performed a steady-state vibration analysis for the model, which was linearized in the vicinity of the undeformed state. Harmonic vibrations were actuated by the voltage at the first patch ϕ 1 =φ 1 sin ωt varying The computed variations of amplitudes of the measured signals are plotted against the excitation frequency ω in a logarithmic scale in Fig. 7. We again observe a good correspondence to the three-dimensional solution even at a coarse finite element mesh. The second resonance peak is missing in the first plot as the corresponding vibration mode is symmetric and not measured by the second sensor, placed in the middle.

Example 2: Static actuation of a hemispherical shell with an opening
In the second example, we study a hemispherical shell with a circular hole of 18 • , shown in the left plot of Fig. 8. The same problem has been studied by Klinkel et al. [52], using a solid shell finite element formulation. The sphere has an inner radius of 0.1 m and consists of a steel kernel sandwiched between two PZT-4 ceramic layers with the thickness values (0.25; 0.5; 0.25) × 10 −3 m; perfect bonding of the layers is assumed. The  upper edge of the shell (the boundary of the opening) is clamped, whereas the lower edge is free to deform under the action of an electrical voltage, applied to the outer and inner piezoelectric layer of the shell: The electrical voltage equals ϕ at the outer piezoelectric layer and −ϕ at the inner piezoelectric layer. This actuation contracts the outer layer and extends the inner one, which is comparable to a change in the curvature of the shell • b in the reference configuration; see (38) and (56). On the other hand, the reference configuration in-plane metric of the shell • a remains unchanged, such that the resulting undeformed state becomes incompatible, and the shell deforms as shown in the middle plot of Fig. 8.

Rotationally symmetric deformation
We start with a semi-analytical solution of the rotationally symmetric problem of small deformation. The geometrically linear counterpart of the system of Eqs. (34), (44), (48), three kinematic boundary conditions at the upper edge (vanishing displacements and no rotation about the edge) and three static boundary conditions at the lower edge (no applied mechanical forces and no moment about the edge) provide a linear boundary value problem. The small thickness of the shell results in narrow boundary layers, in which the moment effects are concentrated, and the membrane stresses dominate in the interior part of the domain. This results in an ill-conditioned boundary value problem, such that particular care had to be taken in order to obtain numerically accurate results when solving it using the shooting method. For the sake of comparing the results to the literature, we used the material parameters from [52] when computing the structural properties of the through-the-thickness element of the shell according to the procedure presented in Appendix A.
A finite element simulation again featured the in-house finite element code ShellFE3. The entire hemisphere was modeled using 64 × 64 elements; we did not make use of the symmetry of the problem to preserve the generality for the subsequent buckling analysis. The applied electrical voltage ϕ = χϕ max was varied from zero to the maximum magnitude ϕ max = 1.4 × 10 5 V using a load factor χ .
At small voltages, the linear semi-analytical solution is nearly identical to the finite element one, as it is demonstrated by comparing the radial deflection of the lower edge for a particular load factor in Table 2.
The results are also sufficiently close to the value, obtained in [52] using a solid shell model. Increasing the load factor, we observed the geometrically nonlinear effects by computing the relative deflection difference in which the linear solution w lin is proportional to χ . The plot in Fig. 9 demonstrates the increasing role of the geometric nonlinearity, as the relative deflection difference increases steadily and reaches approx. 8% at a critical value of the load factor χ * , after which the solution changes qualitatively, the rotational symmetry is lost and the deflection w cannot be uniquely defined.  At χ = χ * ≈ 0.825 the equilibrium path eventually bifurcates into a buckling mode with waves in the circumferential direction. We studied this behavior with the help of the energy norm of the solution. In the finite element simulation, we integrated the mechanical part of the enthalpy H , i.e., the first three terms in the quadratic form (50) with the coefficients A 1,2 , B 1,2 , D 1,2 over the entire domain to compute the strain energy of the shell U FEM . Its variation with an increasing load factor is shown in Fig. 10 in comparison with the linear solution U linear , which is quadratic in χ .
The equilibrium path becomes slightly different from the linear solution at higher electrical voltages, but remains continuous until χ = χ * . Here, the stability of the rotationally symmetric solution is lost and the numerical solver jumps to a neighboring solution branch. We resolved the details of the behavior of the structure using very small load factor increments between sequential computations of the equilibrium states. Shortly after, another jump in the solution is observed (see the loading path in the zoomed area in Fig. 10), and the supercritical deformation mode with 8 waves in the circumferential direction is established. This branch of the solution remains stable as the loading increases further and the buckling mode is becoming more pronounced; the deformed shell for χ = 1 is shown in the right plot of Fig. 8. The unloading branch (dashed line in Fig. 10) is not identical with the loading one, which can be attributed to phenomena such as snap-through and snapback or snap-buckling; see Krommer and Irschik [53] and Krommer et al. [54] for the discussion of such phenomena for the simpler problem of a piezoelectric plate. In any case, we found multiple equilibria in the vicinity of χ * , and we observed a small hysteresis loop.

Conclusions
We have presented a novel multistage hybrid asymptotic-direct approach to the modeling of the nonlinear behavior of thin piezoelectric shells in a holistic form. All three components of the approach, namely asymptotic splitting, direct approach to geometrically nonlinear shells and numerical analysis with a FE scheme, result in a consistent method for the dimensional reduction and analysis of the coupled three-dimensional continuum. For two example problems, numerical results have been computed, verified with three-dimensional solutions obtained with ABAQUS and compared to results from the literature as well as semi-analytical solutions of a linearized version of the developed shell theory.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Material parameters of the shell model
For a material which belongs to the crystal class 2mm with the polarization in the 3 directions (such as PZT ), the linearized 3D constitutive relations can be written in matrix form as: here, the (1, 2) plane is the isotropic plane, and Q 66 = (Q 11 − Q 12 )/2 holds. All materials used in the two example problems in the paper belong to this class, and the thickness direction of the shell is the 3 directions.
In the following, we use effective material parameters defined as For isotropic materials, E and ν are actual material parameters and the piezoelectric parameter e vanishes. Expressing the electric displacement through the electric voltage according to (29), using the asymptotically justified plane stress condition and integrating the enthalpy for given ε, κ and D z0 over z − ≤ z ≤ z + , we arrive at the quadratic form (50). The coefficients A 1 , A 2 , B 1 , B 2 , D 1 , D 2 , p, m and c in the quadratic form are then integral properties of the through-the-thickness element, which are computed according to the procedure below.
We assume the shell to be made of layers i = 1, . . . , n, and the thickness coordinate z varies within each layer between z i−1 and z i ; we have z 0 = z − and z n = z + in terms of Sect. 3.1. The coefficients A 1 , A 2 , B 1 and B 2 are found as To compute D 1 and D 2 , we first introduce the parameters with the aid of which we have Finally, we compute p i , m i and c i , which must be done for each piezoelectric layer: In the examples, four materials were used. The elastic material parameters for aluminum and steel are For the dynamic simulations, we used the value of the mass density ρ = 2700 kg m −3 for aluminum and ρ = 7750 kg m −3 for PZT-5A.