Nonlinear finite element modeling of vibration control of plane rod-type structural members with integrated piezoelectric patches

This paper addresses modeling and finite element analysis of the transient large-amplitude vibration response of thin rod-type structures (e.g., plane curved beams, arches, ring shells) and its control by integrated piezoelectric layers. A geometrically nonlinear finite beam element for the analysis of piezolaminated structures is developed that is based on the Bernoulli hypothesis and the assumptions of small strains and finite rotations of the normal. The finite element model can be applied to static, stability, and transient analysis of smart structures consisting of a master structure and integrated piezoelectric actuator layers or patches attached to the upper and lower surfaces. Two problems are studied extensively: (i) FE analyses of a clamped semicircular ring shell that has been used as a benchmark problem for linear vibration control in several recent papers are critically reviewed and extended to account for the effects of structural nonlinearity and (ii) a smart circular arch subjected to a hydrostatic pressure load is investigated statically and dynamically in order to study the shift of bifurcation and limit points, eigenfrequencies, and eigenvectors, as well as vibration control for loading conditions which may lead to dynamic loss of stability.


Introduction
In recent years, the subject area of theories and finite element tools for modeling and simulation of smart materials and structures have experienced tremendous growth in terms of research and development. Especially for thin-walled structures, integration of smart material layers or patches that can sense as well as be actuated is a promising possibility for shape and vibration control. In this context, modeling and numerical simulation of the static and dynamic response of adaptive beam, plate, and shell structures with integrated distributed control capabilities has attracted considerable research interest in recent years.
Geometrically nonlinear, large deflection problems of structures with integrated smart material layers have been treated considerably less in literature, although these structures are thin and flexible and the deflections often are much larger than the thickness. Therefore, consideration of the structural nonlinearity may be significant to account for phenomena like stress stiffening or softening, shift of eigenfrequencies, shift of bifurcation and limit points, etc., that may affect static as well as dynamic control problems of smart structures.
Concerning static analysis of smart structures, the von Kármán nonlinearity has been taken into account for interlaminar stress analysis of multilayered plates with induced-strain actuators by Icardi and Di Sciuva [29] in the framework of third-order shear deformation theory, and for the analysis of the sensed voltage of a PVDF bimorph cantilever beam by Mukherjee and Chaudhuri [30] adopting first-order shear deformation theory. Sensor output voltage analysis, deflection analysis, and shape control of smart isotropic or composite laminated beams, plates, and shells with integrated piezoelectric actuator and sensor layers have been considered in recent papers (Lentzen and Schmidt [31][32][33][34][35], Vu et al. [36], Nguyen et al. [37], Zhang et al. [38][39][40][41]) accounting for moderate rotations in the framework of first-and third-order shear deformation theories. Various types of nonlinearities were studied by Zhang et al. [42]. Note that the active buckling control and post-buckling analysis was considered by Krishna and Mei [43], Chandrashekhara and Bhatia [44], Wang and Varadan [45], Chróscielewski, Klosowski and Schmidt [46], and Lentzen and Schmidt [47]. Piezothermoelastic problems accounting for the von Kármán nonlinearity have been considered theoretically by Tzou et al. [48], Pai et al. [49], Reddy [50], and numerically by Oh et al. [51,52].
Concerning nonlinear dynamics and control of smart structures, the literature is rather scarce. Yi et al. [53] used solid elements for geometrically nonlinear analysis of the sensor voltage output of surface-bonded piezoelectric patches on beams, plates, and shells. Recently, Mukherjee and Chaudhuri [54] have analyzed the generation of sensed voltage due to nonlinear vibrations of a PVDF bimorph cantilever beam by adopting the von Kármán nonlinearity in the framework of first-order shear deformation theory. The sensor voltage output of surface-bonded piezoelectric patches on beams, plates, and shells has been analyzed in the framework of a nonlinear moderate rotation first-order shear deformation theory in recent papers (Lentzen and Schmidt [55][56][57][58], Rao et al. [60]).
Nonlinear vibration control problems based on Bernoulli-type beam theories have been treated by Shi and Atluri [61], Lee and Beale [62], and recently by Zhou and Wang [63] adopting the von Kármán nonlinearity, and by Chróscielewski, Klosowski and Schmidt [46,64,65] adopting a finite rotation theory. Tzou and Zhou [66,67], Zhou and Tzou [68] studied the active control of nonlinear piezoelectric circular plates and shells with von Kármán-type geometric nonlinearity. Recently, Mukherjee and Chaudhuri [54] adopted the von Kármán nonlinearity in the framework of first-order shear deformation theory for nonlinear vibration control of a PVDF bimorph cantilever beam. Lentzen and Schmidt [55][56][57][58] developed a finite element code for the simulation of nonlinear vibration control of smart isotropic or composite laminated beams, plates, and shells with integrated piezoelectric actuator and sensor layers based on a moderate rotation first-order shear deformation model. Large rotations were analyzed by Zhang and Schmidt [59]. Batra et al. [69,70] dealt with shape and vibration control of plates at finite deformations taking into account also nonlinear constitutive equations for piezoelectric patches. Nonlinear flutter suppression has been considered by Lai et al. [71], Zhou et al. [72,73] based on classical von Kármán large deflection plate theory, and by Shen and Sharpe [74] using discrete Kirchhoff theory finite elements. FE analysis for piezoelectric laminates under strong electric field was performed by Zhang et al. [75], Rao et al. [76]. The magneto-piezoelectric thin-walled structures were studied by Xin and Hu [77], Rao et al. [78]. Piezoelectric materials are commonly used either as sensors or as actuators, since the piezoelectric effect is reversible. However, due to this property, it is possible to make efficient use of these piezoelectric patches by using them simultaneously as sensors and actuators. In this way, it is possible to design controls with high performances and with fewer stability problems, due to the coincidence of the position of the two sensor and actuator functions (see for more detail [79][80][81][82][83]). Nonlinear control often is prone to many issues which require careful consideration. Among them, a very important point is to address loss of stability problem (see, e.g., [84,85]). For a careful modeling of the interfaces between the main mechanical structure and piezoelectric patches, we refer, for example, to [3,4,86,87].
In the present paper, a fully geometrically nonlinear theory of piezolaminated rod-type structures (e.g., plane curved beams, arches, ring shells) based on the Bernoulli hypothesis is developed (Chapter 2). In this framework, a nonlinear finite beam element for the static, stability, and transient analysis of piezolaminated structures is developed that is based on the assumptions of small strains and finite rotations of the normal (Chapter 3). The focus of this paper is on the FE simulation of the large-amplitude vibration response of thin rod-type structures and its control by integrated piezoelectric layers. Our ultimate goal is to expose the complexity of vibration suppression of nonlinear piezolaminated structures and of the related problem of associated finite element approximations. To this end, two problems are studied extensively. In Chapter 4, we consider a clamped semicircular ring shell with piezoelectric actuator and sensor patches bonded to its upper and lower surfaces that was originally treated by Tzou and Ye [23] and has been used as a benchmark problem for linear vibration control in several recent papers, for example, by Sze and Yao [88], Sze et al. [89], Balamurugan and Nayaranan [90]. The refereed results are critically reviewed and extended to account for the effects of structural nonlinearity. In Chapter 5, a circular arch consisting of a master structure with piezoelectric actuator patches bonded to its upper and lower surfaces is investigated statically and dynamically. Nonlinear static analysis of the response to a hydrostatic pressure load is performed in order to investigate the influence of piezoelectric control on the structural stability. The shift of both bifurcation and limit loads associated with unsymmetric and symmetric loss of stability, respectively, due to active stability control is investigated. The eigenvalue problem for small free vibrations about deformed configurations is solved in order to study the shift of eigenforms and eigenfrequencies due to the electric field applied to the actuator patches. Finally, numerical simulation of vibration control is performed, where particular interest is focused on transient loading conditions that may lead to symmetric or unsymmetric dynamic loss of stability. It is demonstrated that both bifurcation and snap-through detected by nonlinear static analysis yield corresponding phenomena under dynamic loading conditions. Finally, let us note that the results and techniques used can be of broader application, for example, one can use not only piezoelectric patches but also ionic polymer metal composites, see, for example, [91,92].

Geometrically nonlinear theory of the 3-D elastic piezoelectric body
In this section, we consider the geometrically nonlinear boundary value problem of the classical elastic piezoelectric body. In Lagrangian description, the basic field equations read: Mechanical equilibrium equations: Div(FS) + f 0 = ρ 0ü ; (2.1) Electrical equilibrium equations: Coupled mechanical constitutive equations: Coupled electrical constitutive equations: Here, S denotes the second Piola-Kirchhoff stress tensor, ε is the Green-Lagrange strain tensor, E and D stand for the electric field and the electric displacement vector, respectively, ρ 0 and f 0 are the mass and body force vectors per unit undeformed volume, respectively, u = x−X is the displacement vector with x and X denoting the position vector in the deformed and undeformed configuration, and C, e, ξ are the elastic moduli, dielectric permittivity, and dielectric tensors, respectively. Here, we also assumed that free charges are vanished. Furthermore, in Eq. (2.1), F is the deformation gradient tensor F = ∇x = ∇(X + u) = 1 + ∇u, (2.5) where the gradient operator is defined as ∇(·) = ∂(·)/∂X. These equations are the starting point for the reduction to a 1D theory of plane curved rods with integrated piezoelectric layers.

Geometrically nonlinear theory of plane piezolaminated rods
In this section, we shall consider the planar deformation of a naturally curved plane rod with integrated piezoelectric layers. We take the (x 1 , x 2 )-plane spanned by the orthonormal basis {e 1 , e 2 } to be the plane of the rod deformation (see Fig. 1).
The subsequent analysis is based on the Bernoulli theory of plane deformations of rods, which is rich enough to accommodate the longitudinal extension and flexure neglecting the transverse shear strains. In the framework of the present theory, the deformation of the rod is completely determined by two components of the displacement vector at the rod axis. The rotation of the normal is not independent variable, but is expressed in terms of the displacement components at the rod axis. In Bernoulli-type rod theory, it is convenient to resolve the displacement vector at the rod axis into the tangential and normal components u t and u n , respectively, relative to the undeformed axis. Based on the Bernoulli kinematical hypothesis, the displacement vector in the 3D rod space at a distance ζ from the axis, restricted to a planar deformation, takes the form where u(s) = u t (s) t(s) + u n (s) n(s) (2.8) is the displacement vector at the rod axis in this 1D problem, with s denoting the arc length coordinate along the undeformed rod axis. Here, we introduce two linearized parameters ϕ and υ, defined by where K is the curvature of the rod, given by K = R −1 ≡ b 11 in terms of the local radius of curvature R and the curvature tensor component b 11 , respectively, with the index 1 denoting the arc length coordinate s. Furthermore, K = , where (s) is the angle between the x 1 -axis and the tangent to the undeformed rod axis. A prime, i.e., (.) , stands for the derivative with respect to the arc length parameter s along the undeformed rod axis, so that u n and u t are dimensionless quantities.
The strain state associated with the above kinematical hypothesis consists of extensional and bending strains. Starting from (2.6), applying the hypothesis (2.7), and using the standard transformations for curved rods, the strain-displacement relations and the stretch of the axis take the form Equations (2.10) and (2.11) are valid for the geometrically nonlinear theory of rods accounting for small strains but large displacements and rotations. The work-conjugate stresses are obtained according to (2.3), consequently yielding the normal stress resultant N and the bending couple M. The reduction of the 3D problem (2.1)-(2.6) to the theory of curved plane rods by invoking the Bernoulli hypothesis leads to a 1D formulation where only the stresses S 11 = 0 and the only independent strains are ε 11 .
We assume that the axis of polarization of the piezoelectric layers is perpendicular to the axis s of the rod and that the electric field is applied to the piezoelectric layers only in normal direction ζ (denoted by the index 3). Under this assumption, Eqs. (2.3) and (2.4) reduce locally to the form Integration of Eq. (2.12) through the cross-sectional area by using (2.10) 1 leads to the stress resultants N and M incorporating a mechanical (m) and an electrical (e) part according to with (l) as layer number. Since a uniform voltage is applied to the piezoelectric layers and these layers are thin, so that the strains can be assumed constant through the layer thickness, Eq. (2.2) is satisfied identically. Finally, taking into account (2.10) and (2.13), the stress-strain relations for a thin linearly elastic piezolaminated rod take the form where In the present study, only the actuation and control problems are investigated. So the leading terms are the mechanical ones as by Zhang et al. [42]. For general formulations of the principle of virtual work for piezoelectric materials, we refer to classic works by Tiersten [93], Allik and Hughes [94], and McMeeking et al. [95], see also Chapter 3 in the recent book by Abali [96]. In the case of shells, the complete form of the principle of virtual work was used, for example, by Lammering [20], Zhang and Schmidt [59], Lammering and Mesecke-Rischmann [97]. Using (2.12) and (2.13), the weak formulation of Eq. (2.1) is obtained as is the virtual work if the inertia forces, where stands for the virtual work of the external loads of surface forces t 0 and body forces f 0 , which must be specified for each type of loading. For the terms G e (2.19), the transition from the 3D to 1D formulation is performed in the standard way, based on assumption (2.7) like in (2.17). The effect of viscous damping can be introduced into the constitutive equations (2.12), or can be directly included in (2.16) by taking into account additional damping forces in the external loads, i.e., by assuming surface forces t 0 as (t 0 − c t 0u ) and/or body forces f 0 as 0 and c f 0 are the respective viscous damping ratios. The principle of virtual displacements (2.16) serves as the basis for a finite element formulation.

Numerical method
In what follows, we restrict our considerations to plane rods with constant curvature. From the numerical point of view, the Bernoulli model (2.16) requires C 1 interelement continuity (Hermitian interpolation) to obtain a conforming finite element method (Fig. 2). To ensure this continuity requirement, we consider the 1D 2-node element with four degrees of freedom at each node resulting in the interpolation scheme for the components of the displacement vector according to

1) with the nodal element parameters
in each element e and node (a), respectively.
Here, α ξ (ξ ) denotes the Lame coefficient, ξ ∈ [0, +1] is the natural coordinate, while H (a) (ξ ) and H (a) (ξ ), a = 1, 2, are the cubic Hermitian shape functions, respectively, given for the 2-node element by In contrast to the standard finite element approximation, in this paper, the initial geometry of the arch is represented exactly within the element. Substituting the interpolation scheme (3.1) into the principle of virtual work (2.16) and following again the standard finite element procedure, one gets where (.) t ≡ (.)(t) denotes the value of the respective function at time t. Here, q T = {u (1) u (2) . . . u (a) . . .}, q ≡ ∂ t q,q ≡ ∂ tq denote the global vectors of nodal displacements, velocities, and accelerations, respectively, and is the discrete form of equilibrium equations (2.16). In (3.5), M = const and C t are the mass and damping matrices, is the reference load and λ(t) the load function. Let us note that in the piezoelectricity these matrices may have special properties which can be used for numerics, see, e.g., [96,[98][99][100][101]. The element matrices in (3.5) and those resulting from the consistent linearization of (2.16) can next be calculated in the standard way. In contrast to the finite element class C 0 , this element does not exhibit locking effects and therefore no special techniques to avoid them, for example reduced integration, are needed. For the evaluation of element matrices, we apply the Gauss-Legendre numerical integration with 3 sampling points. The linearization process and the use of the standard FEM approximation procedure (3.1) yield the classical incremental form of the equations of motion where β and γ are the parameters of the Newmark method (β = 1 4 , γ = 1 2 ). Sufficiently small time steps t have to be employed to assure convergence of the iteration process to be performed at each time step. However, this does not always guarantee stability of the numerical procedure in case of nonlinear problems. This fact must be taken into account while investigating the techniques of vibration control.
The nonlinear static analysis and the analysis of small free vibrations around the current, geometrically nonlinear configuration are based on the equations where ω i and a i are the i-th eigenfrequency and eigenvector, respectively.
To trace the static equilibrium paths, we consider the extended space R N +1 of size (N + 1) with elements being the increments q T = { q, λ} ∈ R N +1 . Then, the computation of nonlinear equilibrium paths may be described by the compact set of equations (see, e.g., Waszczyszyn [102]): Here, equation (3.9) 1 is supplemented with a constraint equation that is imposed on the vector q. By specifying the vectort, one obtains various control techniques: [103,104], Wempner [105]. (3.10) Here, λ * , q * , and s * are given values of the control parameters and where s i are weighting coefficients that allow choosing certain components of the generalized displacement increment vector q∈R N (and also j). We introduce the matrix S, since different units and orders of magnitudes of components of q (and also j) may cause problems with a reasonable definition of the arc length. To overcome this problem, Chróścielewski and Nolte [106] proposed the technique of selective elimination or scaling.
Here, for convergence control, the criterion of selective displacements and selective unbalanced forces was where e q and e j are, respectively, the assumed error range for displacements and unbalanced forces, ||Sδ( q)|| is the norm of the iterative change δ( q) of the increment q of the selected components of the displacement vector, ||S q|| is the norm of the total increment of the selected components of the displacement vector, ||Sj|| is the norm of the selected unbalanced force components, and p ref is the given reference value chosen appropriately for type of load under consideration. The secondary path (see Fig. 3) is obtained using load-displacement control along with the control window technique developed by Chróscielewski and Schmidt [107].
Other possible ways to treat the problem under study, particularly suited for nonlinear curved beams, are, for example, isogeometric analysis or Hencky-type lumped discretization (see, e.g., [108][109][110][111]). For nonlinear dynamic problems and its FEM analysis, we also refer to [112,113].   Figure 4 shows a laminated semicircular ring shell, which was considered in Tzou and Ye [23] in the context of geometrically linear finite element simulation of small-amplitude vibration control. The structure consists of a steel master structure with segmented piezoelectric layers bonded to the top and bottom surfaces.
The outer piezoelectric layer serves as a distributed actuator and the inner one as a distributed sensor, both being made of lead zirconate titanate piezoceramics (PZT). The steel ring is b = 5.08 cm wide, h = 0.635 cm thick, and L = 100 cm long, i.e., the curvature radius of the shell midsurface is R = 31.831 cm. The thickness of the PZT layers is h PZT = 254 µm, and that of the steel layer is h steel = 0.5842 cm. The material parameters are given in Table 1.
Our numerical analysis concerns three different problems: 1. Convergence study, linear eigenvalue problem for small free vibrations and static analysis, 2. Control of quasi-free small-amplitude vibrations, 3. Control of nonlinear dynamic problems. Figure 5 shows a convergence study for the first eigenfrequency of the ring shell, where the structure was discretized using CAM shell elements of C 0 -class (see Chróścielewski, Makowski, and Stumpf [112][113][114][115]), and C 1 -class beam finite elements (LUK), respectively, described in Chapter 3. In Table 2, results for the first eight eigenfrequencies are given using a 4 × 50 mesh of 16-node C 0 -class CAM shell elements with full integration (FI, 4 × 4 integration points), and 100 2-node C 1 -class beam finite elements. In the former analysis also, out-of-plane vibration modes are detected, while the results for the in-plane modes are in excellent agreement with those of the beam finite element analysis. In Fig. 5, N denotes the number of nodal points along the ring length in the respective mesh and in the symbol describing the CAM element type e4, e8, e9, e16 denotes 4-, 8-, 9-, and 16-node elements, while FI, URI, and SRI stand for full, uniform reduced, and selective reduced integration, respectively. The first eigenfrequency given in Tzou and Ye [23] differs from our converged result by approximately 25% and so do the higher eigenfrequencies. It can be seen from Fig. 5 that the results of Tzou and Ye [8] are probably obtained with an element which exhibits some locking effect and with a mesh which is too coarse.  Table 2 Eigenfrequencies Additionally, we have included in Table 2 results which can be found in Sze and Yao [88] and Sze et al. [89], obtained by a hybrid-stress-assumed natural strain 8-node solid-shell element and by the 4-node generalpurpose (thin and thick) curved shell element of ABAQUS S4. Based on these results, these two papers criticize the element of Tzou and Ye [23] as too stiff. Unfortunately, following an error in Table 2 of Tzou and Ye [23], Sze and Yao [74] and Sze et al. [89] have used a wrong value of Young's modulus for steel, E = 68.95 GPa instead of E = 210 GPa. Additionally, they do not realize that in Tzou and Ye [23], only in-plane modes are considered, i.e., they compare, for example, their second eigenfrequency (for an out-of-plane mode) with the second in-plane eigenfrequency of Tzou and Ye [23], and so forth. Therefore, we have recalculated their results for E = 210 GPa by multiplying their numbers by the factor √ 210/68.95 and put them in the right order (see Table 2). After these modifications, the results of Sze and Yao [88] and Sze et al. [89] are in reasonable agreement with our results. It can be also observed (Fig. 5) that the conforming elements (CAM (FI), LUK) converge monotonically in the typical way for eigenfrequencies.

Convergence study
Because of the observed differences described above, now we perform an additional convergence study, in a local sense, for the internal stress resultants. We consider a static loading of the shell by a tip normal force P n(a) = 1 N at point (a). Since the deformations are small and the system is statically determined, the displacements as well as the distribution of the normal force and bending moments, etc., can be computed easily. In Fig. 6, the results of the convergence study for the internal stress resultants are displayed. These results refer to the integration points of the finite elements.
It can be observed that a good approximation for the bending moment M(s) can be obtained already with a mesh of ten elements, which yields nearly the same results as a mesh of 100 elements. For the normal force N (s), however, the discretization by only ten elements is insufficient. Even if only the results for the element midpoints are considered, the error for the normal force is not acceptable. Additionally, in Fig. 6, the results for the discretization by 20 elements, corresponding to the mesh used in Tzou and Ye [23], and for 100 elements are given. Since the latter results are satisfactory, in the following a mesh of at least 100 elements will be used.  It should be noted that the present convergence study is of h-type, being performed using 2-node elements of the class C 1 . Elements of the class C 0 or with a different number of nodes may show different convergence behavior. Figures 7 and 8 show the first five in-plane vibration modes of the structure, where u n and u t denote the normal and tangential components of the displacement vector u, respectively, and ||u t (s)|| = u t (s)/u max , ||u n (s)|| = u n (s)/u max , with u max = max s∈L (u n (s), u t (s)). It can be observed that the first vibration mode is dominated by the tangential displacement components u t , while the other modes shown in these graphs are dominated by the normal displacement components u n . It can be seen that the distribution of the normal force N (s) is more complicated than that of the bending moment M(s)t, which is typically for curved beams and shells. These forms can provide some useful information about suitable positions and lengths of actuator and sensor patches depending on the modes to be suppressed.

Control of quasi-free small-amplitude vibrations
The following investigations concern the uncontrolled and controlled transient response of the ring shell to impulsively applied loads. A negative velocity proportional feedback algorithm is used in the analysis. The geometry of the ring and used notations are shown in Fig. 4. This study differs from the one reported by Tzou and Ye [23], where an initial velocity is imparted to the free end of the structure (u t (a) = −1ms −1 ). However, in the framework of structural mechanics, rather an initial velocity field instead of a velocity in a point should be prescribed, because otherwise in finite element analysis, the structural response is mesh dependent. This is even more pronounced when the consistent mass matrix approach is used, which is demonstrated in Fig. 11.
In the following finite element analyses, the ring shell is discretized by 100 C 1 -class 2-node beam elements (LUK) described in Chapter 3, i.e., ten elements per segment, which is sufficient to assure FE convergence. A viscous damping of 10N s/m 4 is assumed in analogy to the 0.2% initial damping used by Tzou and Ye [23]. For the time integration of the equations of motions, Newmark method (3.7) with the time step t = 2 ×10 −6 s and convergence criterion (3.12) for the iteration process with e q = 0.0001 and e j = 0.0001 were used.
First, we study the transient uncontrolled and controlled response of the ring shell for small-amplitude vibrations due to a force impulsively applied at the tip point (a). The time dependency is defined by the load function λ(t), which is linearly increased in a short time interval of 0.001s to the value 1, then kept constant, and decreased linearly in the time interval between 0.019 and 0.02 s (see Fig. 12). Figure 13 shows the uncontrolled tip point response to a tangential force P t (a) (t) = P t,re f (a) λ(t), P t,re f (a) = −1N impulsively applied at the free end of the structure. It can be seen from Fig. 13 that this load history leads to small-amplitude vibrations.
Next, distributed vibration control of the ring shell subjected to the same load history is studied with different lengths of actuators and sensors. The number of piezoelectric actuator and sensor patches varies from 1 to 10 (see Fig. 4) starting with one patch located near the clamped end (i.e., 10% of the upper and lower surface is covered) to ten patches (i.e., fully covered surfaces). All actuator patches are charged with the same voltage, but such that opposite action of the respective actuators at the upper and lower surface is obtained. A negative velocity proportional feedback algorithm is applied, in which the actuator voltage V is taken proportional to the current average normal velocity v n,average ≡u n,average of all sensor segments, i.e., V = pu n,average . Here, p is a proportionality coefficient chosen in the example under consideration as p = 10000V s/m, andu n,average is determined as the average value of the normal velocities in the center points of all finite elements located in the area covered by sensor patches. Figures 14, 15, and 16 show the controlled transient response of the ring shell with a 2-patch, 5-patch, and 10-patch actuator pair, respectively. In all cases, a rapid decay of the vibration amplitude is observed in the graphs of both the tangential and normal tip point displacements. It can be seen that the control effect increases strongly from the 2-patch to the 5-patch actuator, but the difference between the damping actions of the 5-and  A different problem is encountered when instead of a tangential force P t (a) (t), a normal force P n(a) (t) is applied in an impulsive manner at the free end (a) of the ring shell. In what follows, the time history for the force P n(a) (t) is chosen in the same way as described above for P t (a) (t) (see Fig. 12). Figure 17 shows the uncontrolled tip point response to a normal force P n,re f (a) = 1 N impulsively applied at the free end (a). It can be observed from Fig. 17 that this impulsive normal force leads to high-frequency small-amplitude vibrations. Figure 18 shows the result for the numerical simulation of the controlled tip point response of the ring shell with a 10-patch actuator pair, i.e., when the surface is 100% covered. One can observe that in this case, a damping of both the normal and tangential tip point displacements is achieved, which is good, but not as rapid as for the tangential impulsive force P t (a) (t) (see Fig. 16). Furthermore, for a smaller number of patches, the result is even worse: For the 5-patch actuator pair, the damping action is not evident (Fig. 19) and for the 2-patch actuator pair instead of damping, an excitation of the vibration was observed (Fig. 20).
In order to investigate the observed differences of the control effectiveness in more detail, we analyze the evolution of the uncontrolled deformations during time. Figures 21, 22, 23, and 24 show the results for the numerical simulation of the uncontrolled response for both impulsive forces P t (a) (t) and P n(a) (t) (see Fig. 12).  It can be observed that the deformations and velocity distributions caused by the impulse P n(a) (t) are more complicated than for P n(a) (t). This corresponds to the observations already made in Figs. 13 and 17. Both results illustrate that for the impulse P n(a) (t), the participation of higher modes in the motion is strong. This can be the reason for the low effectiveness (Fig. 18) or even failure (Figs. 19 and 20) of the applied control method.

Control of nonlinear dynamic problems
Finally, the transient uncontrolled and controlled response of the structure is studied in the geometrically nonlinear range of deformation. To this end, a tangential force P t (a) (t) = P t,re f (a) λ(t), P t,re f (a) = −300N with higher load intensity is applied at the free end (a) of the structure. Additionally, the pulse duration is increased to 0.2s (see Fig. 25).
Like in Sect. 4.1.2, the ring shell is discretized by 100 C 1 -class 2-node beam elements described, i.e., ten elements per segment. A viscous damping of 10 Ns/m 4 is assumed in analogy to the 0.2% initial damping used in Tzou and Ye [23]. For the time integration of the equations of motions, the Newmark method (3.7) with the time step t = 2 × 10 −6 s and the convergence criterion (3.12) for the iteration process with e q = 0.0001 and e j = 0.0001 were used.  For the investigation of the control effect in the geometrically nonlinear range of deformation with different numbers of piezoelectric actuator patches, we follow the procedure described above: A negative velocity proportional feedback algorithm is applied, in which the actuator voltage V is taken proportional to the current average normal velocity v n,average ≡u n,average of all sensor segments, i.e., V = pu n,average . Here, p is a proportionality coefficient chosen in the example under consideration as p = 20,000V s/m, andu n,average is determined as the average value of the normal velocities in the center points of all finite elements located in the area covered by sensor patches. Figures 28, 29, and 30 show the results for the numerical simulation of the controlled transient largeamplitude response of the ring shell with a 2-patch, 5-patch, and 10-patch actuator pair, respectively. It can be observed that only for the 10-patch actuator pair, the effectiveness of the control technique is acceptable. For this case, the corresponding sequence of the deformed configurations shows rapid suppression of the vibration within ∼ 2s (see Fig. 31).
Additionally, in Figs. 32, 33, and 34, the control effect is illustrated by the sequences of the decaying normal velocities, bending moments, and normal forces.       Next, we check the control technique for the case when the impulsively applied force acts in normal direction, i.e., P n,re f (a) = 300 N with the same load function (see Fig. 25) is applied at the free end (a) of the ring shell. As already shown for the small vibration analysis, in this case, the structural response is related to higher eigenmodes. Figure 35 shows the resulting uncontrolled large-amplitude high-frequency tip point vibrations versus time. Figure 36 shows the corresponding sequence of the deformed configurations every 0.5 s. Figures 37, 38, and 39 show the controlled large-amplitude response of the ring shell with a 2-, 5-, and 10patch actuator pair, respectively. Like in the small-amplitude vibration control only with the 10-patch actuator pair, a damping effect was achieved, while here for the 2-patch actuator pair, an excitation of the vibration is observed.
The result for the 5-patch actuator pair shows that in the nonlinear range, the additional problem of numerical stability of the time integration algorithm requires special attention. This is due to the fact that even algorithms that are unconditionally stable for linear dynamics lose this property in the nonlinear case (see, e.g., Kuhl and Crisfield [116], Chróścielewski et al. [117]. In this context, reduction of the time step can be helpful, but does not solve this problem. For the case of the 10-patch actuator pair, the corresponding sequence of the deformed configurations shows suppression of the vibration (see Fig. 40). Additionally, for the case of the 10-patch actuator, in Figs. 41, 42, and 43, the control effect is illustrated by the sequences of the decaying normal velocities, bending moments, and normal forces in time. It can be observed that the control effect for the impulsive normal force (P n,re f (a) = 300 N) is much worse than for the tangential impulse (P t,re f (a) = −300 N). This is especially visible in the graph for the internal normal force, which remains in the order of magnitude of 100 N during the entire time range.
Next, we repeat the investigation of the control effect as above except the fact that the voltage applied to the actuator patches in the negative velocity proportional feedback algorithm is limited to a maximum value of V max = 250 V. First, for the tangential impulse (P t,re f (a) = −300 N), Figs. 44, 45, and 46 show the results for the controlled transient large-amplitude response of the ring shell with a 2-patch, 5-patch, and 10-patch actuator pair, respectively. It can be observed that in all cases, the control effect is slower, but in contrast to the case of unlimited actuator voltage, here for the 2-patch and 5-patch actuator pairs, the control method leads to a suppression of the vibrations. Now, for the normal impulse (P n,re f (a) = 300 N), the results for the controlled transient large-amplitude response with limited actuator voltage V max = 250 V are presented in Figs. 47, 48, and 49. As in the case of  The results show that the relatively simple negative velocity proportional feedback control law can be successfully applied when the vibration of the structure is primarily dominated by the first mode. When the response is dominated by higher modes, more sophisticated control techniques are required (see, e.g., Mukherjee and Joshi [118,119] for adaptive shape design of continuously distributed actuator and sensor profiles or Lentzen and Schmidt [120] for the application of modal sensor arrays for the active damping of higher modes of vibrations around a geometrically nonlinear deformed equilibrium state). Figure 50 shows a laminated circular arch clamped at both ends and subjected to a uniformly distributed hydrostatic pressure loading.   The arch consists of a steel master structure with two layers of piezoelectric PVDF polymer patches bonded to its upper and lower surfaces, respectively. Each piezoelectric layer consists of ten equally segmented distributed actuators. The geometry of the axis of the arch is determined by the distance B = 120 mm and the angle α = 30 • . Accordingly, the radius of the arch is R = 1 2 B sin 1 2 α = 231.182 mm, and its height is H = 1 2 B(sin 1 2 α− tan 1 2 α) = 7.899 mm. The width of the arch is b = 10 mm, the thickness of the steel layer is h s = 0.2 mm, and the thickness of the piezoelectric patches is h PVDF = 0.028 mm each. The material parameters are given in Table 3.

Circular arch with piezoelectric patches
We assume that each piezoelectric patch can be charged independently. The voltages applied to the upper and lower piezoelectric patch, respectively, of an individual segment are chosen equal in magnitude but with different sign.
The arch is subjected to a hydrostatic pressure, t 0 = −q 0n (u), where q 0 = qdS/d S is the value of the pressure measured per unit area of the undeformed configuration, q(λ) = q n,ref λ where λ is the load parameter, andn(u(λ)) is the unit normal vector of the deformed axis.
The arch was discretized by 100 2-node C 1 -class finite elements, i.e., ten elements per segment, which is sufficient to assure FE convergence. For the time integration of the equations of motions, Newmark method (3.7) with the time step t = 2 × 10 −6 s and convergence criterion (3.12) for the iteration process with e q = 0.0001 and e j = 0.0001 were used.  The static as well as the nonlinear dynamic analyses concern both the symmetric (snap-through) and asymmetric (bifurcation) loss of stability of the arch.

Nonlinear static analysis
In the framework of the geometrically nonlinear static analysis, the load-deflection behavior of the arch was investigated first without any piezoelectric effects. For this case, Figs. 51 and 52 show a load-deflection curve with a limit point of the primary path at q lim = 5300.7N/m 2 , u n,lim(c) = −0.00012m (associated with symmetric snap-through) and bifurcation points of a secondary path (associated with asymmetric buckling forms) at q cr1 = 3338.5 N/m 2 , u n,cr1(c) = − 0.00000998 m and q cr2 = 684.3 N/m 2 , w n,cr2(c) = − 0.009826 m, respectively. For the primary path, the Riks-Wempner constant arc length method was applied. The secondary path (see Figs. 52 and 53) was obtained using load-displacement control along with the control window technique developed by Chróścielewski and Schmidt [107]. In order to determine the associated bifurcation points, a perturbation force of the intensity P t,impf(c) = 1 N was applied in tangential direction at the middle point (c) of the arch. This perturbation force was subsequently removed when the tangential displacement |u t (c) | in the middle point of the arch became greater than u t,rem(c) = 5 × 10 −5 m. The critical value of the hydrostatic pressure is determined for zero tangential displacement at the middle point (c) of the arch (see Fig. 53).
Next, the influence of the piezoelectric effect on the load-deflection response of the arch was investigated. To this end, the piezoelectric patches of the six sections located symmetrically in the middle of the arch were charged with a negative voltage for the top layers and a positive voltage for the bottom layers. The two sections located at the left and right clamped edge, respectively, were charged in the opposite way, i.e., with a positive voltage for the top layers and a negative voltage for the bottom layers. This voltage distribution was chosen in order to cause a state of prestress in the arch that acts against the applied loading in the initial (prebuckling) range of deformation.
The value of the voltage was V = 1000 V for each patch. Figures 51 and 52 show that the primary path changes qualitatively: Additional loops with respective limit points appear. For the secondary path, only quantitative changes are observed: The values of the bifurcation points increase with the applied voltage (see     actuator voltage decreases strongly. For dynamics, one has to take into account also the greater mass density of the PZT ceramic: ρ PZT ρ PVDF = 7600 kg/m 3 2840 kg/m 3 ≈ 2.676.

Eigenvalue analysis
Next, an analysis of the eigenvalue problem for small free vibrations about geometrically nonlinear deformed configurations was carried out. Figure 54 shows the change of the first three eigenfrequencies of the arch as function of the load parameter q for different values of the applied voltage V . The distribution of the voltage applied to the piezoelectric patches was the same as described in Sect. 4.1. In Table 4, results for the first five eigenfrequencies for q = 0 and V = 0 are given using 100 2-node C 1 -class beam finite elements. They are used as reference values in Fig. 54.

Control of nonlinear dynamic problems
For vibration control, the voltage V applied to the piezoelectric patches in each individual segment was chosen equal in magnitude and with different sign for the top and bottom patches, but now proportional to the current average normal velocityu n,average of the respective segment, i.e., V = pu n,average . Here, p is a proportionality coefficient chosen as p = 10,000 Vs/m. The maximum voltage was restricted to the breakdown voltage assumed to be V max = 1000 V. Velocity feedback has been used, for example, for dynamic instability control of piezolaminated columns by Mukherjee and Saha Chaudhuri [121][122][123].
The aim of the present analysis is the simulation of vibration control for the strongly nonlinear arch problem involving bifurcation and snap-through phenomena, which should play an important role also in the dynamic response of the arch.
To this end, three qualitatively different ranges of the hydrostatic pressure are investigated: (a) less than the bifurcation load: q(t) < q cr1 = 3338.5N /m 2 , (b) between the bifurcation load and the limit load: q cr1 = 3338.5N /m 2 < q(t) < q lim = 5300.7N /m 2 , (c) greater than the limit load: q(t) > q lim = 5300.7N /m 2 .
In the first load case, the hydrostatic pressure is linearly increased in a short time interval of 0.0001s to the value q and then kept constant as shown by the step load function in Figure 55 where λ(t) is the load function and q(t) = q n,ref λ(t).
Concerning the case (a), a hydrostatic pressure q ≡ 3300N /m 2 < q cr 1 = 3338.5N /m 2 was chosen. In order to check whether a bifurcation phenomenon occurs, a perturbation load like in the static analysis was applied (P t,impf(c) = 1N with u t,rem(c) = 5 × 10 −5 m). It can be seen in Figs. 56, 57, and 58 that under this loading condition, no snap-through phenomenon occurs. If the piezoelectric control is applied, the biggest components of the nonlinear vertical and horizontal vibrations are damped rapidly. However, since in this case, the horizontal displacement |u t (c) | remains smaller than the chosen value of u t,rem(c) = 5 × 10 −5 m, the perturbation load is not removed and the arch reaches a state of small asymmetric vibrations. For completeness, the same problem was investigated with a smaller value of u t,rem(c) = 5 × 10 −6 m. In this case, the perturbation load was removed quickly. In both cases, the velocities are much less than 0.1m/s and the applied voltage does not reach V = 1000V (see Figs. 56 and 57 for the case u t,rem(c) = 5 × 10 −6 m).
Concerning the case (b), a hydrostatic pressure q cr1 = 3338.5N/m 2 < q ≡ 3400N/m 2 < q lim = 5300.7N/m 2 was chosen, i.e., a value only slightly above the bifurcation load. Bearing in mind the results of the static analysis, a bifurcation phenomenon should be expected in this case. Two different perturbation loads When piezoelectric vibration control is applied, a rapid damping is observed before the loss of stability occurs and the time until loss of stability is substantially prolonged. In this range, the velocities are relatively small (less than 0.1 m/s) and the applied voltage does not reach V =1000 V. However, when the loss of stability occurs, the displacements, velocities, and accelerations increase very rapidly and the applied voltage would be always much greater than V =1000 V. As a consequence, the results of the control method (|V | ≤ 1000 V) applied here are not sufficient. However, the bifurcation phenomenon cannot be avoided. It should be remarked that also if no perturbation load is applied and an asymmetric movement is admitted, after a certain period of time, an asymmetric loss of stability is predicted. The reason for this is the accumulation of numerical errors, which play the role of a perturbation. This can be avoided if symmetry conditions are imposed which prohibit any asymmetric movement (see Figs. 61, 62 and 63 for q = 5000N/m 2 . Concerning the case (c), a hydrostatic pressure q ≡ 5310N/m 2 > q lim = 5300.7 N/m 2 was chosen. In order to avoid the possibility of asymmetric loss of stability, symmetry conditions were imposed. In Figs. 64, 65, and 66, the symmetric snap-through phenomenon can be observed. The remarks concerning the efficiency of the control method in the case (b) are also valid here. Next, an impulsive hydrostatic pressure q(t) = q n,ref λ(t) was used: The load was applied as mentioned before, but it was removed linearly in the interval between 0.0099 and 0.01 s (see Fig. 55). Two different values of loading were applied, q = 5500 N/m 2 and q = 6000 N/m 2 , both above the static limit load.
In the first case (see Fig. 67), the maximum amplitude of vibrations in the middle point of the arch is less than 1mm and no dynamic snap-through is reached. The piezoelectric control effectively damps the vibrations. In the second case (see Fig. 68), although the loading is not much bigger than in the first one, dynamic snapthrough occurs and the maximum amplitudes are more than ten times bigger than in Fig. 67. Even after the loading is removed (q = 0), large-amplitude vibrations can be observed. Also for the case of piezoelectrically controlled vibrations, dynamic snap-through is reached. However, the control method rapidly decreases the amplitude of the vibrations. For this type of impulsive loading, the proposed control method turned out to be especially efficient for the damping of both vibrations without and with dynamic loss of stability, because in the latter case, the amplitudes of vibrations return into the range of the prebuckling state after the loading is removed.

Conclusions
In the present paper, a fully geometrically nonlinear theory of piezolaminated rod-type structures based on the Bernoulli hypotheses and a corresponding finite element method for the simulation of static, stability, and vibration control have been presented. The numerical algorithm has been applied to two problems: (i) vibration control of a clamped piezolaminated semicircular ring shell subjected to impulsive tangential or radial tip loads, (ii) static and dynamic analysis of a piezolaminated circular arch under hydrostatic pressure loading conditions which may lead to static or dynamic loss of stability.
The conclusions drawn from the finite element analysis of the piezolaminated semicircular ring shell may be summarized as follows: • It was shown that lack of proper convergence studies for the eigenfrequencies is the source of disagreements between results published in the recent literature. It was also shown that additional convergence studies, in a local sense, for the internal stress resultants are of outmost importance in order to avoid significant errors. • For the case of small-amplitude vibrations of the ring shell due to a tangential force impulsively applied at the tip, the present study of the control effect performed with different lengths of actuators and sensors confirms the results of Tzou and Ye [23] concerning the damping ratio versus number of actuator patches. • The results presented for vibration control in case of an impulsively applied radial tip force, leading to high-frequency small-amplitude vibrations, are unprecedented in the literature. In this case, a lower control effectiveness is observed and explained by a detailed investigation of the evolution of the uncontrolled deformations during time. • For both loading cases (i.e., tangential and radial impulsive tip load, respectively), vibration suppression is studied also in the geometrically nonlinear range of deformation and the control effect is illustrated by the sequences of the decaying normal velocities, bending moments, and normal forces.
The example of a piezolaminated circular arch under hydrostatic pressure loading is used to study (i) the nonlinear static response taking into account the effect of the voltage applied to the piezoelectric patches, (ii) the eigenvalue problem for small free vibrations about geometrically nonlinear deformed configurations with different voltages applied to the piezoelectric patches, and (iii) the vibration control of nonlinear transient problems. The static as well as the dynamic analyses concern both the symmetric (snap-through) and asymmetric (bifurcation) loss of stability of the arch.
• It was shown that phenomena like bifurcation and snap-through occurring in nonlinear static analysis have corresponding counterparts also in dynamic analysis. Therefore, from the computational point of view, the static analysis of stability may be an important tool for the prognosis of the dynamic behavior of structures. • Concerning the nonlinear static analysis of structures, it was shown that the position of both bifurcation and limit points can be shifted using appropriate active piezoelectric stability control. • Likewise, by appropriate piezoelectric control, eigenfrequencies and eigenvectors of the structure can be changed. • For structures under transient loading, the considered vibration control method can suppress the vibrations, or delay the loss of dynamic stability, which can be important from the technical point of view. • For step loading conditions, the proposed control method is effective for the damping of vibrations in the prebuckling range. For dynamic impulsive loading conditions, it was also efficient for free vibrations after a dynamic loss of stability.
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.