Geometrically nonlinear finite element simulation of smart laminated shells using a modified first-order shear deformation theory

Abstract

This article investigates geometrically nonlinear and linear analysis of multilayered shells with integrated piezoelectric materials. An efficient nonlinear shell element is developed to solve piezoelastic response of laminated structure with embedded piezoelectric actuators and sensors. A modified first-order shear deformation theory is introduced in the present method to remove the shear correction factor and improve the accuracy of transverse shear stresses. The electric potential is assumed to be a linear function through the thickness of each active sub-layer. Several numerical tests for different piezolaminated geometries are conducted to highlight the reliability and efficiency of the proposed implementation in linear and geometrically nonlinear finite element analysis.

Keywords

Modified first-order shear deformation theory geometrically nonlinear analysis sensors/actuators piezoelectric material four-node shell

Introduction

Smart structures have attracted intensive research interests because of their potential benefits in a wide range of applications, such as shape control, vibration suppression, noise attenuation, and damage detection (Gabbert et al., 2017; Gandhi and Thompson, 1992). The coupled mechanical and electrical properties of piezoelectric materials make them well suited for use as sensors and actuators in smart structures, named also intelligent structures. These structures offer great potential for use in advanced aerospace as well as hydrospace, nuclear and automotive structural applications due to their excellent electromechanical properties, easy fabrication, design flexibility, and efficiency to convert electrical energy into mechanical energy (Fereidoon et al., 2016; Foda et al., 2010; Ren and Azar, 2006; Zhang et al., 2007, 2011).

Many researches deal with analytical (Benjeddou et al., 2002; Bohlooly and Mirzavand, 2015; Kulikov and Plotnikova, 2013; Kulikov et al., 2015; Lam and Ng, 1999; Torres and Paulo de Tarso, 2010; Xiang and Shi, 2008) and semi-analytical (Daneshmehr and Shakeri, 2007) techniques to predict piezoelastic response of laminated structure with integrated piezoelectric layers. Nevertheless, analytical approach does not allow analysis of complicated structures, it can only be used for simple structures. To analyze complicated structures with integrated piezoactive layers, finite element method (FEM) is required. Three different approaches can be used in order to investigate the kinematic of shell structures, which are classical laminated theory (CLT), first-order shear deformation theory (FSDT), and high-order shear deformation theory (HSDT).

A large number of developed elements adopted the CLT that implements Kirchhoff–Love kinematics. Moita et al. (2002) studied the geometrically nonlinear analysis of laminated structures with embedded piezoelectric layers, using the classical Kirchhoff–Love hypothesis. The obtained results exhibit that the transverse displacements are significantly influenced by the geometrically nonlinear analysis. He et al. (2001) established a generic finite element code to model a curved functionally graded material (FGM) shell with piezoelectric sensors and actuators for the static and dynamic control of shells. However, the formulation based on the CLT theory does not take into account the transverse shear strains (Jrad et al., 2018a). Gabbert et al. (2002) selected the SemiLoof element, originally proposed by Irons (1976), for an extension to a piezoelectric shell element to predict the active behavior of a smart composite structure.

Although CLT theories have proved their efficiency to describe the behavior of thin shells by their simple implementation in the most finite element codes and low computational cost, they neglect the effect of transverse shear deformations. To overcome the limitation of CLT theory, Mindlin (1951) proposed the FSDT theory considering the effects of a constant transverse shear deformation. The formulation of FSDT theory is the simplest and has been employed by many researchers for developing finite element for the analysis of smart laminated shells. Recently, a three-node shell element was developed by Marinković and Rama (2017) and Rama (2017) in order to predict the static and dynamic piezoelastic response of piezoelectric-laminated composite shells. In this work, geometrically effects may easily become important in the behavior of shell structures which is related to kinematic boundary conditions and load level. The three-node finite shell element was also proposed by Neto et al. (2012) in order to predict piezoelastic static and dynamic response of smart laminated structures. Using a four-node degenerated shell element based on the coupled FSDT, an analysis of smart laminated piezoelectric plate and shell is conducted in the work of Kogl and Bucalem (2005). Yet, a nine-node piezoelectric shell element was established by Marinković and Nestorović (Marinković et al., 2006, 2009, 2012; Marinkovic and Zehn, 2015; Nestorović et al., 2012, 2014) and Balamurugan and Narayanan (2008) since the degenerated eight-node shell element is prone to notable locking effects, and so, the addition of the ninth node offers better behavior simulation.

Recently, some research works have also been done, using FSDT theory, to analyze the behavior of FGM shells with integrated piezoactive layers. In fact, FGMs have gained intensive consideration as special composites whose composition changes continuously through the thickness of structure, due to their high performance, ensuring smooth transition of stress distributions and novel thermo-mechanical properties (Kidane and Shukla, 2008; Müller et al., 2003; Pompea et al., 2003). The static and dynamic piezo-thermo-elastic analyses of FGM plates with integrated piezoelectric sensor and actuator layers were investigated by Liew et al. (2001). In this study, torsional and bending vibrational controls were carried out on the FGM plates under thermal loading, using self-monitoring sensors and actuators. Moreover, Ng et al. (2002) formulated a flat-shell element to predict the frequency response using a constant gain displacement and velocity feedback control algorithm. In another work, Liew et al. (2004) presented a generic finite element formulation for analyzing vibration control of a shell laminate using self-monitoring sensor and self-controlling actuator layers. The obtained results demonstrate that mode shapes and resonance frequencies could be controlled appropriately by adjusting the displacement control gain. Another type of the FGM called functionally graded piezoelectric materials (FGPM) obtained even more attention in recent years. These types of materials can be integrated with the host structure to avoid deterioration in interlayers bonding strength at high temperature. Studies dealing with static, dynamic, and free vibration responses of functionally graded piezoelectric shell structures can be found in the literature (Bodaghi et al., 2012, 2014; Joshi et al., 2003; Lee, 2005).

The FSDT theory seems to be insufficient to examine thick shell structures and it is necessary to introduce the shear correction coefficients which can be restrictive for such applications. Indeed, shear correction coefficients can be simply acquired for linear isotropic material (5/6), but their determination is more complex for composites, especially laminated structures (Hajlaoui et al., 2015). Hence, HSDT theories were proposed to overcome the shear coefficient limitations (Wali et al., 2015; Zghal et al., 2018, Mallek et al., 2018). This class of kinematic theories has been recently equipped with piezoelectric capability. For example, Correia et al. (2002) developed a finite shell element based on the high-order displacement field to study active control of axisymmetric shells with piezoelectric layers. The electromechanical behavior of piezoelectric generic shells with graded material properties in the thickness direction was examined by Wu et al. (2002). Furthermore, parametric studies are presented to exhibit the effects of graded material properties on the piezoelectrically induced displacements, stresses, electric potential, and electric displacements distributions. Sudhakar and Kamal (2003) also used HSDT theory to examine the active vibration control performance of the curved beam with integrated sensors and actuators. Recently, a unified formulation (developed by Cinefra and Carrera, 2013; Cinefra and Valvano, 2016; Cinefra et al., 2010, 2012), which can generate any refined theory, is developed in static and free vibration for laminate composites and FGM shells.

Despite HSDT theories provide a refined approximation of the displacements and deformations of the structure, the number of used kinematic variables is mainly high which leads to a high computational time effort. Therefore, attention has turned to modified FSDT shell type finite elements, which can provide a satisfying accuracy with acceptable numerical effort (Mellouli et al., 2018). Various studies have been made to correct the shear strain distribution. In addition, Tanov and Tabiei (2000) implemented a shear function in the FSDT formulation leads to a parabolic distribution of transverse shear stresses across the thickness. Inspired from the investigation of Shi (2007), the present formulation presents an improved FSDT theory by imposing a parabolic shear strain distribution across the shell thickness and a zero shear stress on the top and bottom faces of the shell.

Geometrically, nonlinear problems of structures with integrated piezolayers have been treated considerably less in the literature, although these structures undergo moderate deformation because of their flexible nature, and hence, geometrical nonlinearity has to be accounted for in the analysis (Jrad et al., 2018b). It has been shown in the literature (Kulkarni and Bajoria, 2007; Lee, 2005; Lentzen et al., 2007; Schulz and Klinkel, 2008) that in many cases, the geometrical nonlinearities cannot be neglected.

This article investigates the static piezoelastic response of laminated structure with integrated piezoelectric layers under mechanical and electrical loads. A very important aspect of the article is the extension of formulation to geometrically nonlinear analysis based on the improved FSDT. This new theory has been used by imposing a parabolic through thickness distribution for the transverse shear strains and zero transverse shear stresses requirements at the shell surfaces. The analysis is based on the FEM by using four-node shell element and the electric potential is assumed to be linear through the thickness of the piezoelectric layer. In addition, the performance of the element in modeling the actuator performance and the sensor response for different piezolaminated geometries are investigated. A set of numerical analyses is performed in order to highlight the applicability and effectiveness of the present finite element model notably for smart structures.

Geometrically nonlinear piezoelectric shell formulation

The objective of this work is to examine the linear and geometrically nonlinear behavior of piezolaminated shells under mechanical and electrical loads. The present theory is based on the modified FSDT. The reference surface of the shell is assumed to be smooth, continuous, and differentiable. The initial and the deformed configurations are denoted by $C_{0}$ and $C_{t}$ , respectively. Variables associated to initial configuration $C_{0}$ (respectively, $C_{t}$ ) are denoted by upper-case letters (respectively, lower-case letters). Vectors and tensors are expressed using bold letters.

Piezoelectric modified FSDT shell kinematic assumptions

Curvilinear coordinates $ξ = (ξ^{1}, ξ^{2}, ξ^{3} = z)$ are used in this work to define the curved surface of the shell. As shown in Figure 1, $X_{q}$ and $x_{q}$ are the position vectors of point (q) in the reference and current configurations, respectively. $X_{p}$ and $x_{p}$ are the position vectors of point (p) in the mid-surface, which is the reference surface of the shell. The displacement field in both initial and the deformed configurations are defined as follows

{\begin{matrix} X_{q} (ξ^{1}, ξ^{2}, z) = X_{p} (ξ^{1}, ξ^{2}) + z D (ξ^{1}, ξ^{2}) \\ x_{q} (ξ^{1}, ξ^{2}, z) = x_{p} (ξ^{1}, ξ^{2}) + z d (ξ^{1}, ξ^{2}) \end{matrix}, z \in [- \frac{h}{2}, \frac{h}{2}]

(1)

where h is the thickness, $D$ and $d$ are the director vectors of the reference and current configurations, respectively, as illustrated in Figure 1. The vector $D$ is perpendicular to the undeformed mid-surface.

The covariant vectors in the initial state $C_{0}$ are given by

G_{α} = A_{α} + z D_{, α}; G_{3} = D α = 1, 2

(2)

The surface element dA in the initial state is given by

dA = \sqrt{A} d A_{ξ}, \sqrt{A} = ‖ A_{1} \land A_{2} ‖, d A_{ξ} = d ξ^{1} d ξ^{2}

(3)

The covariant reference metric tensor $G$ at the material point is defined by

G = [G_{i} \cdot G_{j}], i, j = 1, 2, 3

(4)

The volume element $dV$ in the initial state is given by

dV = \sqrt{G} d ξ^{1} ξ^{2} dz \sqrt{G} = [G_{1} G_{2} G_{3}] = \sqrt{| G_{ij} |}

(5)

The metric tensor in the deformed configuration $C_{t}$ can be divided into in plane and out of plane parts. Hence, it is expressed as

g_{i j} = g_{i} . g_{j}, {\begin{matrix} g_{α β} = a_{α β} + z b_{α β} + o (z^{2}) \\ g_{α 3} = c_{α} \end{matrix} α, β = 1, 2

(6)

where $a_{α β}$ , $b_{α β}$ , and $c_{α}$ represent kinematic variables, which are expressed as

{\begin{matrix} \begin{matrix} a_{α β} = a_{α} . a_{β} \\ c_{α} = a_{α} . d \end{matrix} \\ {b_{α β}}_{} = a_{α} . d_{, β} + a_{β} . d_{, α} \end{matrix}

(7)

The Green–Lagrangian strain tensor $ε$ can be described by means of membrane, bending, and transverse shear strains using the following expression (Frikha et al., 2018b; Mars et al., 2017)

ε = \frac{1}{2} (g - G) ε_{ij} = \frac{1}{2} (g_{ij} - G_{ij}), {\begin{matrix} ε_{α β} = e_{α β} + z χ_{α β} \\ 2 ε_{α 3} = γ_{α} \end{matrix}

(8)

$e_{α β}$ , $χ_{α β}$ , and $χ_{α β}$ denote, respectively, membrane, bending, and transverse shear strains. The vectors of membrane, bending, and shear strains are summarized in matrix notation in the following way (Frikha et al., 2016)

\begin{array}{l} δ e = {\begin{matrix} δ e_{11} \\ δ e_{22} \\ 2 δ e_{12} \end{matrix}} = B_{m} . δ x, \\ δ χ = {\begin{matrix} δ χ_{11} \\ δ χ_{22} \\ 2 δ χ_{12} \end{matrix}} = B_{b m} δ x + B_{b b} δ d, \\ δ γ = {\begin{matrix} δ γ_{1} \\ δ γ_{2} \end{matrix}} = B_{s m} δ x + B_{s b} δ d \end{array}

(9)

where matrix differential operators $C_{0}$ are detailed in the following subsections.

Membrane strains

The membrane strains can be computed as

e_{α β} = \frac{1}{2} (a_{α β} - A_{α β}), α, β = 1, 2

(10)

The virtual membrane strains in $C_{t}$ can be written in the following way

δ e_{α β} = \frac{1}{2} (a_{α} . δ x_{, β} + a_{β} . δ x_{, α}), α, β = 1, 2

(11)

In matrix form, the membrane strain vector is given by

δ e = B_{m} . δ x, B_{m} = [\begin{matrix} a_{1}^{T} \frac{\partial}{\partial ξ^{1}} \\ a_{2}^{T} \frac{\partial}{\partial ξ^{2}} \\ a_{1}^{T} \frac{\partial}{\partial ξ^{2}} + a_{2}^{T} \frac{\partial}{\partial ξ^{1}} \end{matrix}]

(12)

Bending strains

The bending strains are given by

χ_{α β} = \frac{(b_{α β} - B_{α β})}{2}, α, β = 1, 2

(13)

The variation of the bending strains can be written in $C_{t}$ as

\begin{array}{l} δ χ_{α β} = \frac{1}{2} (a_{α} . δ d_{, β} + a_{β} . δ d_{, α} + δ x_{, α} . d_{, β} + δ x_{, β} . d_{, α}), \\ α, β = 1, 2 \end{array}

(14)

In matrix notation, the bending strain vector is given by

δ χ = B_{bm} δ x + B_{bb} δ d

(15)

where

\begin{array}{l} B_{b b} = [\begin{matrix} a_{1}^{T} \frac{\partial}{\partial ξ^{1}} \\ a_{2}^{T} \frac{\partial}{\partial ξ^{2}} \\ a_{1}^{T} \frac{\partial}{\partial ξ^{2}} + a_{2}^{T} \frac{\partial}{\partial ξ^{1}} \end{matrix}] a n d \\ B_{b m} = [\begin{matrix} d_{, 1}^{T} \frac{\partial}{\partial ξ^{1}} \\ d_{, 2}^{T} \frac{\partial}{\partial ξ^{2}} \\ d_{, 1}^{T} \frac{\partial}{\partial ξ^{2}} + d_{, 2}^{T} \frac{\partial}{\partial ξ^{1}} \end{matrix}], k = 1, 2 \end{array}

(16)

Shear strains

The shear strains can be expressed as

γ_{α} = c_{α} - C_{α}

(17)

The variation of the shear strains is defined as follows

δ γ_{α} = a_{α} . δ d + δ x_{, α} . d

(18)

In matrix notation, the shear strain vector is given by

δ γ = B_{sm} δ x + B_{sb} δ d

(19)

where

B_{sm} = [\begin{matrix} d^{T} \frac{\partial}{\partial ξ^{1}} \\ d^{T} \frac{\partial}{\partial ξ^{2}} \end{matrix}] and B_{sb} = [\begin{matrix} a_{1}^{T} \\ a_{2}^{T} \end{matrix}]

(20)

A modified FSDT shell model

It should be mentioned that based on FSDT theory, the transverse shear strains is assumed to be linear across the thickness (Mindlin, 1951). Nevertheless, it is already well known that the shear (and hence strain) stress distribution is parabolic across the thickness vanishing at a point ( $ξ^{1}$ , $ξ^{2}$ , ±h/2) on the top and bottom surfaces of the shell (Levinson, 1980; Reddy, 1984; Tanov and Tabiei, 2000). Inspired by the investigation of Shi (2007) and Murthy (1981), the present formulation constitutes an improved FSDT theory by assuming a parabolic shear strain through the thickness using a function (see Table 1)

f (z) = \frac{5}{4} (1 - \frac{4 z^{2}}{h^{2}})

(21)

Table 1.

The shear theories functions.

Theories	Shear function $f (z)$
FSDT: Mindlin (1951)	$1$
HSDT: Levinson (1980), Reddy (1984), Frikha et al. (2018a), Zghal et al. (2018a)	$1 - \frac{4 z^{2}}{h^{2}}$
FSDT: Tanov and Tabiei (2000)	$1 - \frac{4 z^{2}}{h^{2}}$
HSDT: Shi (2007), Murthy (1981)	$\frac{5}{4} (1 - \frac{4 z^{2}}{h^{2}})$
Modified FSDT: Present theory	$\frac{5}{4} (1 - \frac{4 z^{2}}{h^{2}})$

FSDT: first-order shear deformation theory; HSDT: high-order shear deformation theory.

The shear strains vectors becomes

δ γ^{Z} = f (z) δ γ

(22)

Electric field

The electrical field $E$ is evaluated based on the gradient of the electric potential $φ$ . Its expression is given by

E = - φ_{, α}, α = 1 \dots 3

(23)

In this work, since the piezoelectric layer is considered as thin shell structure with polarization in the thickness direction, the in-plane electric field $E_{1}$ and $E_{2}$ are assumed to be zero and only transverse electric field $E_{3}$ is considered. The electric field of an element is expressed in terms of the electrical potential as

E = - B_{e} . φ B_{e} = [\begin{matrix} 0 \\ 0 \\ \frac{\partial}{\partial z} \end{matrix}]

(24)

Denoting $E_{m}$ as the opposite of the electrical field $E$ , the relationship between $E_{m}$ and $φ$ is defined as

E_{m} = - E = B_{e} . φ

(25)

Weak form

The weak form of equilibrium equations can take the following form

\begin{matrix} G (Φ, δ Φ) = \int_{V} (S_{α β} δ ε_{α β} + S_{3 α} δ γ_{α}^{Z} + q_{3} δ E_{3}) dV - G_{ext} (Φ, δ Φ) = 0; \\ α, β = 1, 2 \end{matrix}

(26)

where $δ Φ = (δ x, δ d, δ φ)$ is an arbitrary variation, $δ ε_{α β}$ are the covariant components of the virtual Green–Lagrange strain tensor, $S_{α β}$ are the contravariant components of the second Piola–Kirchhoff stress tensor, $q_{α}$ is the electric displacement, $δ E_{α}$ is the virtual electric vector, and $G_{ext}$ is the external virtual work.

Using equation (9), the weak form becomes

G = \int_{A} (N . δ e + M . δ χ + T . δ γ + \bar{q} . δ E) dA - G_{ext} = 0

(27)

where $δ e$ , $δ χ$ , and $δ γ$ are the variations of the shell strains; $N$ , $M$ , and $T$ are the membrane, bending, and shear stress resultants; and $\bar{q}$ is the dielectric displacement which can be summarized in the following vectors as

N = {\begin{matrix} N^{11} \\ N^{22} \\ N^{12} \end{matrix}}, M = {\begin{matrix} M^{11} \\ M^{22} \\ M^{12} \end{matrix}}, T = {\begin{matrix} T^{1} \\ T^{2} \end{matrix}}, \bar{q} = {\begin{matrix} q^{1} \\ q^{2} \\ q^{3} \end{matrix}}

(28)

Their components are defined as follows

\begin{matrix} N^{α β} = \int_{- h / 2}^{h / 2} S_{α β} \sqrt{\frac{G}{A}} dz M^{α β} = \int_{- h / 2}^{h / 2} z S_{α β} \sqrt{\frac{G}{A}} dz \\ T^{α} = \int_{- h / 2}^{h / 2} f (z) S_{α 3} \sqrt{\frac{G}{A}} dz \bar{q} = \int_{- h / 2}^{h / 2} q \sqrt{\frac{G}{A}} dz \end{matrix}

(29)

The expressions of generalized resultant of stress $R$ and strain $Σ$ vectors are given by

R = {\begin{matrix} N \\ M \\ T \\ \bar{q} \end{matrix}}_{11 \times 1}, Σ = {\begin{matrix} e \\ χ \\ γ \\ E_{m} \end{matrix}}_{11 \times 1}

(30)

where $E_{m}$ represents the opposite of the electric field $E$ .

Hence, the weak form of the equilibrium equation can be rewritten as

G (Φ, δ Φ) = \int_{A} δ Σ^{T} . R dA - G_{ext} (Φ, δ Φ) = 0

(31)

The Newton–Raphson iterative method is used for the solution of nonlinear equation. Indeed, the consistent tangent operator is built up by providing the directional derivative of the weak form in the direction of the increment $Δ Φ = (Δ u, Δ d, φ)$ . To analyze large displacement, the tangent operator is decomposed into material and geometric parts, denoted by $D_{G} G . Δ Φ$ and $D_{M} G . Δ Φ$ , respectively

DG . Δ Φ = D_{G} G . Δ Φ + D_{M} G . Δ Φ

(32)

Material part

The material part of the tangent operator, which results from the variation in the stress resultants, can be written as follows

D_{M} G . Δ Φ = \int_{A} δ Σ^{T} Δ R dA

(33)

The coupling between the mechanical and electrical material behavior can be described through the constitutive relations, taking into account the direct and converse piezoelectric effect

{\begin{matrix} S = C ε - p^{T} E \\ q = p ε + kE \end{matrix}

(34)

where $C$ , $p$ , and $q$ represent, respectively, the elastic matrix for a constant electric field, the piezoelectric coupling matrix, and the dielectric permittivity matrix. In the three-dimensional case, the elastic matrix is given by

C = [\begin{matrix} c_{11} & c_{11} & c_{12} & 0 & 0 & 0 \\ c_{22} & c_{23} & 0 & 0 & 0 \\ c_{33} & 0 & 0 & 0 \\ c_{44} & 0 & 0 \\ Sym & c_{55} & 0 \\ c_{66} \end{matrix}]

(35)

The components of the third-order piezoelectric constant tensor $p$ and the second-order dielectric constant tensor $k$ in equation (34) can be arranged, respectively, in matrix form as

\begin{array}{l} p = [\begin{matrix} 0 & 0 & 0 & 0 & e_{14} & 0 \\ 0 & 0 & 0 & 0 & 0 & e_{25} \\ e_{31} & e_{32} & e_{33} & 0 & 0 & 0 \end{matrix}] and \\ k = [\begin{matrix} k_{11} & 0 & 0 \\ 0 & k_{22} & 0 \\ 0 & 0 & k_{33} \end{matrix}] \end{array}

(36)

For thin shell structures, a vanishing stress component $S_{33}$ is considered and the components in the x- and y-directions for the electric field can be neglected. Elimination of the strain component $ε_{33}$ finally yields for the in-plane stresses and the transverse dielectric displacement

[\begin{matrix} S_{11} \\ S_{22} \\ S_{12} \\ q_{3} \end{matrix}] = [\begin{matrix} {c^{*}}_{11} & {c^{*}}_{12} & 0 & {e^{*}}_{31} \\ {c^{*}}_{22} & 0 & {e^{*}}_{32} \\ {c^{*}}_{44} & 0 \\ Sym & {k^{*}}_{33} \end{matrix}] [\begin{matrix} ε_{11} \\ ε_{22} \\ ε_{12} \\ E_{m 3} \end{matrix}]

(37)

and for the transversal shear stresses and in-plane dielectric displacements

[\begin{matrix} S_{13} \\ S_{23} \\ q_{1} \\ q_{2} \end{matrix}] = [\begin{matrix} c_{55} & 0 \\ 0 & c_{66} \\ e_{14} & 0 \\ 0 & e_{25} \end{matrix}] [\begin{matrix} ε_{13} \\ ε_{23} \end{matrix}]

(38)

The plane stress components of the material constants are computed as follows

\begin{matrix} {c^{*}}_{11} = c_{11} - \frac{{(c_{13})}^{2}}{c_{33}} = \frac{Y (z)}{1 - ν^{2} (z)} \\ {c^{*}}_{22} = c_{22} - \frac{{(c_{23})}^{2}}{c_{33}} = \frac{Y (z)}{1 - ν^{2} (z)} \\ {c^{*}}_{12} = c_{12} - \frac{c_{13} c_{23}}{c_{33}} = \frac{Y (z) ν (z)}{1 - ν^{2} (z)} \\ {c^{*}}_{44} = c_{44} = \frac{Y (z)}{2 (1 + ν (z))} = c_{55} = c_{66} \\ {e^{*}}_{31} = e_{31} - \frac{c_{13}}{c_{33}} e_{33} = e_{31} - \frac{ν (z)}{1 - ν^{2} (z)} e_{33} \\ {e^{*}}_{32} = e_{32} - \frac{c_{23}}{c_{33}} e_{33} = e_{32} - \frac{ν (z)}{1 - ν^{2} (z)} e_{33} \\ {k^{*}}_{33} = k_{33} - \frac{{(e_{33})}^{2}}{c_{33}} = k_{33} - \frac{(1 + ν (z)) (1 - 2 ν (z))}{Y (z) (1 - ν (z))} {(e_{33})}^{2} \end{matrix}

(39)

with $Y (z)$ is the Young modulus and $ν (z)$ is Poisson’s ratio.

In the case of an elastic isotropic constitutive model, the variation in stress resultants $Δ R$ is related to the variation in strain field.

\begin{matrix} Δ R = H_{T} Δ Σ \end{matrix}

(40)

with $H_{T}$ is the linear coupling elastic and electric matrix expressed as

\begin{matrix} H_{T} & = [\begin{matrix} H_{11} & H_{12} & 0 & H_{14} \\ H_{22} & 0 & H_{24} \\ H_{33} & H_{34} \\ Sym & H_{44} \end{matrix}]; \\ {\begin{matrix} (H_{11}, H_{12}, H_{22}) = \int_{- h / 2}^{h / 2} (1, z, z^{2}) C dz \\ H_{33} = \int_{- h / 2}^{h / 2} {(f (z))}^{2} C_{τ} dz \\ (H_{14}, H_{24}) = \int_{- h / 2}^{h / 2} (1, z) p_{1}^{T} dz \\ H_{34} = \int_{- h / 2}^{h / 2} f (z) p_{2}^{T} dz \\ H_{44} = \int_{- h / 2}^{h / 2} k dz \end{matrix} \end{matrix}

(41)

where $C$ and $C_{τ}$ are in-plane and out-of-plane linear elastic sub-matrices, which can be written as

C = [\begin{matrix} {c^{*}}_{11} & {c^{*}}_{12} & 0 \\ {c^{*}}_{22} & 0 \\ Sym & {c^{*}}_{44} \end{matrix}], C_{τ} = [\begin{matrix} c_{44} & 0 \\ 0 & c_{55} \end{matrix}]

(42)

where $p_{1}^{T}$ , $p_{2}^{T}$ , and $k$ represent the in-plane and out-of-plane piezoelectric coupling sub-matrices and dielectric permittivity matrix, respectively, that can be expressed as

\begin{matrix} p_{1}^{T} = [\begin{matrix} 0 & 0 & {e^{*}}_{31} \\ 0 & 0 & {e^{*}}_{32} \\ 0 & 0 & 0 \end{matrix}], p_{2}^{T} = [\begin{matrix} e_{14} & 0 & 0 \\ 0 & e_{25} & 0 \end{matrix}], \\ k = [\begin{matrix} k_{11} & 0 & 0 \\ 0 & k_{22} & 0 \\ 0 & 0 & {k^{*}}_{33} \end{matrix}] \end{matrix}

(43)

The material part of the tangent operator becomes

D_{M} G . Δ Φ = \int_{A} δ Σ^{T} Δ R dA = \int_{A} δ Σ^{T} H_{T} Δ Σ dA

(44)

Using equations (12), (15), (19), and (24), the internal virtual work becomes

G_{int} = δ {Φ_{n}}^{T} \int_{A} B^{T} R dA; B = [\begin{matrix} B_{m} \\ B_{b} \\ B_{s} \\ B_{e} \end{matrix}]

(45)

where $B$ represents the stain–displacement matrix.

The material tangent operator is deduced from equation (45)

K_{m} = [\begin{matrix} K_{m}^{uu} & K_{m}^{u φ} \\ K_{m}^{φ u} & K_{m}^{φ φ} \end{matrix}] = \int_{A} B^{T} H_{T} B dA

(46)

where

\begin{array}{l} K_{m}^{u u} = \int_{A} {[\begin{matrix} B_{m} \\ B_{b} \\ B_{s} \end{matrix}]}^{T} [\begin{matrix} C & 0 \\ 0 & C_{τ} \end{matrix}] [\begin{matrix} B_{m} & B_{b} & B_{s} \end{matrix}] d A; \\ K_{m}^{u φ} = \int_{A} {[\begin{matrix} B_{m} \\ B_{b} \\ B_{s} \end{matrix}]}^{T} {[\begin{matrix} p_{1} & 0 \\ 0 & p_{2} \end{matrix}]}^{T} [B_{e}] d A \\ K_{m}^{φ u} = \int_{A} {[B_{e}]}^{T} [\begin{matrix} p_{1} & 0 \\ 0 & p_{2} \end{matrix}] [\begin{matrix} B_{m} & B_{b} & B_{s} \end{matrix}] d A; \\ K_{m}^{φ φ} = \int_{A} {[B_{e}]}^{T} k [B_{e}] d A \end{array}

(47)

Geometrical part

The geometric part results from the variation of the virtual strains while holding stress resultants constant

D_{G} G . Δ Φ = \int_{A} (Δ δ Σ^{T} . R) dA

(48)

This expression can be divided into membrane, bending, and shear terms as follows

D_{G} G . Δ Φ = D_{G} G_{m} . Δ Φ + D_{G} G_{b} . Δ Φ + D_{G} G_{s} . Δ Φ

(49)

The weak form equation (31), the material part equation (44), and the geometrical part equation (49) will be presented in the following section.

Finite element formulation

In this section, the finite element implementation of the proposed nonlinear modified FSDT shell to predict the electromechanical responses of the structure is detailed. A four-node element with an isoparametric shape function $N^{I}$ is considered (Figure 2)

N^{I} (ξ^{1}, ξ^{2}) = \frac{1}{4} (1 - ξ^{1^{I}} ξ^{1}) (1 - ξ^{2^{I}} ξ^{2}) I = 1 \dots 4

(50)

Each node possesses five mechanical and one electrical degrees of freedom (DOF): $u, v, w, θ_{x}, θ_{y}, φ$ .

The interpolation of the displacement vector $(U = x - X)$ can be written as

U = \sum_{I = 1}^{4} N^{I} U_{I}, Δ U = \sum_{I = 1}^{4} N^{I} Δ U_{I}

(51)

where $U_{I} \in R^{3}$ is the nodal displacement.

The variation and increment of the director vector $d$ are approximated with the same functions:

δ d = \sum_{I = 1}^{4} N^{I} δ d_{1 I}, Δ d = \sum_{I = 1}^{4} N^{I} Δ d_{1 I}

(52)

The electric potential $φ$ is approximated with the same shape functions as shown

φ = \sum_{I = 1}^{4} N^{I} φ_{I}

(53)

Local Cartesian coordinate

In order to move from curvilinear system to the Cartesian one, the relation between the derivations of $N^{I}$ in the local Cartesian and the local elementary systems can be written as

[\begin{matrix} {\bar{N}}_{, 1}^{I} \\ {\bar{N}}_{, 2}^{I} \end{matrix}] = {[J]}^{- 1} [\begin{matrix} N_{, 1}^{I} \\ N_{, 2}^{I} \end{matrix}]

(54)

where $J$ is the Jacobian transformation from the local Cartesian base ${n_{1}^{0}, n_{2}^{0}, n^{0}}$ to the covariant base ${A_{1}, A_{2}, A_{3}}$ as shown in Figure 3. It can be defined by

J = [\begin{matrix} n_{1}^{0} . A_{1} & n_{2}^{0} . A_{1} \\ n_{1}^{0} . A_{2} & n_{2}^{0} . A_{2} \end{matrix}]

(55)

The normal field $n^{0}$ to the mid-surface, in the initial state $C_{0}$ , can be evaluated as

n^{0} = \frac{A_{1} \land A_{2}}{‖ A_{1} \land A_{2} ‖}

(56)

Approximation of membrane and bending strain field

The discretization of the membrane and bending parts of the strain field are given by

δ e = B_{m} . δ Φ_{n}, δ χ = B_{b} . δ Φ_{n}

(57)

where $B_{m}$ and $B_{b}$ are the membrane and the bending strain–displacement matrices, respectively. Their expressions are given at nodal level I in Table 2.

Table 2.

The discrete membrane and bending strain–displacement operators.

B_{m}^{I} = [\begin{matrix} B_{mm}^{I} & 0 & 0 \end{matrix}]

B_{mm}^{I} = [\begin{matrix} n_{1}^{T} {\bar{N}}_{, 1}^{I} \\ n_{2}^{T} {\bar{N}}_{, 2}^{I} \\ n_{1}^{T} {\bar{N}}_{, 2}^{I} + n_{2}^{T} {\bar{N}}_{, 1}^{I} \end{matrix}]

I = 1, . . ., 4

B_{b}^{I} = [\begin{matrix} B_{bm}^{I} & B_{bb}^{I} & 0 \end{matrix}]

B_{bm}^{I} = [\begin{matrix} d_{, 1}^{T} {\bar{N}}_{, 1}^{I} \\ d_{, 2}^{T} {\bar{N}}_{, 2}^{I} \\ d_{, 1}^{T} {\bar{N}}_{, 2}^{I} + d_{, 2}^{T} {\bar{N}}_{, 1}^{I} \end{matrix}]

B_{bb}^{I} = [\begin{matrix} n_{1}^{T} {\bar{N}}_{, 1}^{I} \\ n_{2}^{T} {\bar{N}}_{, 2}^{I} \\ n_{1}^{T} {\bar{N}}_{, 2}^{I} + n_{2}^{T} {\bar{N}}_{, 1}^{I} \end{matrix}]

d_{, a} = \sum_{I = 1}^{4} {\bar{N}}_{, a}^{I} d_{I}, a = 1, 2, I = 1, . . ., 4

The unit vectors of the actual basis $(n_{1}, n_{2})$ in Table 2 are defined as

\begin{matrix} n_{, 1} = \sum_{I = 1}^{4} {\bar{N}}_{, 1}^{I} {(X + U)}_{I} \\ n_{, 2} = \sum_{I = 1}^{4} {\bar{N}}_{, 2}^{I} {(X + U)}_{I} \end{matrix}

(58)

Approximation of shear strain field

It should be noted that the interpolation of the shear strain is obtained by using Assumed Natural transverse Shear strain method (ANS Method). This advanced technique is introduced in order to avoid shear locking of the developed element. This has been presented for a purely mechanical element formulation in work of Bathe and Dvorkin (1985). Accordingly, the shear strain is expressed in the middle of each edge of the element by

δ γ_{ξ} = [\begin{matrix} δ γ_{1} \\ δ γ_{2} \end{matrix}] = [\begin{matrix} \frac{1}{2} (1 - ξ^{2}) δ γ_{1} (B) + \frac{1}{2} (1 + ξ^{2}) δ γ_{1} (D) \\ \frac{1}{2} (1 - ξ^{1}) δ γ_{2} (A) + \frac{1}{2} (1 + ξ^{1}) δ γ_{2} (C) \end{matrix}]

(59)

in which $γ_{2} (A)$ , $γ_{1} (B)$ , $γ_{2} (C)$ , and $γ_{1} (D)$ represent the strains at points A, B, C, and D, respectively, as defined in Figure 2.

The variation of shear strain can be rewritten as

δ γ_{ξ} = B_{s ξ} . δ Φ_{n}

(60)

where $B_{s ξ}$ is the discrete shear strain–displacement operator

\begin{matrix} B_{s ξ} = [\begin{matrix} N_{, 1}^{1} d_{B}^{T} & N_{, 1}^{2} A_{1 B}^{T} & 0 & N_{, 1}^{2} d_{B}^{T} & N_{, 1}^{2} A_{1 B}^{T} & 0 \\ N_{, 2}^{1} d_{A}^{T} & N_{, 2}^{4} A_{2 A}^{T} & 0 & N_{, 2}^{2} d_{C}^{T} & N_{, 2}^{3} A_{2 C}^{T} & 0 \end{matrix} \\ \begin{matrix} N_{, 1}^{3} d_{D}^{T} & N_{, 1}^{3} A_{1 D}^{T} & 0 & N_{, 1}^{4} d_{D}^{T} & N_{, 1}^{3} A_{1 D}^{T} & 0 \\ N_{, 2}^{3} d_{C}^{T} & N_{, 2}^{3} A_{2 C}^{T} & 0 & N_{, 2}^{4} d_{A}^{T} & N_{, 2}^{4} A_{2 A}^{T} & 0 \end{matrix}] \end{matrix}

(61)

The transverse shear strain can, hence, be expressed in the local Cartesian system as

δ γ = B_{s} . δ Φ_{n}

(62)

where $B_{s}$ denotes the discrete shear strain–displacement matrix in the local Cartesian system. Its expression is given as follows

B_{s} = J^{- 1} B_{s ξ}

(63)

Approximation of electric field

The electric field is interpolated as follows

δ E_{m} = B_{e} . δ Φ_{n}

(64)

where $B_{e}$ represents the discrete electric field–displacement vector

{B_{e}}^{I} = [\begin{matrix} 0 & 0 & B_{ee}^{I} \end{matrix}]; {B_{ee}}^{I} = {\begin{matrix} 0 \\ 0 \\ \frac{N^{I}}{t} \end{matrix}}

(65)

where t represents the thickness of the piezoelectric layer.

Linearization of weak form

The virtual and incremental generalized strain can be represented in approximate form as follows

Δ Σ = B . Δ Φ, δ Σ = B . δ Φ

(66)

where $B = [\begin{matrix} B_{m} \\ B_{b} \\ B_{s} \\ B_{e} \end{matrix}]$

Using equation (45), the internal virtual work becomes

G_{int} = δ {Φ_{n}}^{T} \int_{A} B^{T} R dA

(67)

The discrete material tangent operator is deduced from equation (44)

K_{m} = \int_{A} B^{T} H_{T} B dA

(68)

In matrix form, the geometric tangent operator, based on equation (61), becomes

D_{G} G . Δ Φ = δ {Φ_{n}}^{T} . (K_{Gm} + K_{Gf} + K_{Gc}) . Δ Φ_{n}

(69)

Membrane, bending, and shear contributions can be grouped to form the geometric tangent operator

D_{G} G . Δ Φ = δ {Φ_{n}}^{T} . K_{G} . Δ Φ_{n}

(70)

The global geometric tangent operator is detailed in the Appendix. Note that $K_{G}$ is symmetric.

The global stiffness matrix is expressed as

K = K_{m} + K_{G}

(71)

Nodal transformation

The variation of the directors $δ d$ and their derivatives $δ d_{, α}$ can be expressed either in spatial description, in all equations, as

δ d = δ θ \land d = Λ δ θ, \bar{Λ} = - \tilde{d}

(72)

where $\tilde{d}$ is the skew-symmetric tensor such that $\tilde{d} . d = 0$ , or in material description

δ d = Q δ \bar{Θ} E_{3} = \bar{Λ} δ Θ, \bar{Λ} = Q {\tilde{E}}_{3}

(73)

where $d = Q E_{3}$ , $E_{3} = [{\begin{matrix} 0 & 0 & 1] \end{matrix}}^{t}$ and $Q = [\begin{matrix} t_{1} & t_{2} & t_{3} \end{matrix}]$ .

A spatial description leads to a shell problem with seven DOF/node and the material description leads to a shell problem with six DOF/node. The transformation $\bar{Λ}$ takes the following form

\bar{Λ}_{n} = {[\begin{matrix} - t_{2 n} & t_{1 n} \end{matrix}]}_{3 \times 2}; n = 1, \dots, 4

(74)

The relation between the generalized displacement vector $δ Φ_{n} = (δ x, δ d, δ φ)_{n}$ and the nodal displacement $δ Γ_{n} = (δ x, δ Θ, δ φ)_{n}$ vector is defined as

Φ_{n} = Π_{n} δ Γ_{n}; n = 1, \dots, 4

(75)

where $Π_{n} = [\begin{matrix} I & 0 & 0 \\ 0 & \bar{Λ_{n}} & 0 \\ 0 & 0 & 1 \end{matrix}]$

The global stiffness matrix becomes

K_{T} = Π^{T} K Π

(76)

where $Π = diag (Π_{1}, Π_{2}, Π_{3}, Π_{4})$ .

Numerical examples and discussion

In order to validate the numerical performance of the formulation, the results of linear and geometrically nonlinear static problems are presented. Comparisons are made between the present obtained results and existing solutions available in literature for different piezolaminated geometries. The standard quadrilateral four-node shell element with three translational, two rotational, and one electrical DOF per node is used to model all proposed geometries.

Table 3 lists the electromechanical material properties of all materials used in this section according to Lammering and Yang (2009), Tzou and Ye (1996), Nestorović et al. (2012), Marinković et al. (2006); Marinkovic and Rama, 2017; Marinkovic and Zehn, 2015) and Moita et al. (2002).

Table 3.

Electromechanical material properties.

	T300/976	S-glass/Epoxy	PXE-52	PVDF	PIC 151	PZT-G1195	PZT
Elastic properties
Y11 (GPa)	150.0	55.0	62.5	2.0	61.3	63.0	63.0
Y22 (GPa)	9.0	16.0	62.5	2.0	61.3	63.0	63.0
$ν_{12}$	0.30	0.28	0.30	0.29	0.38	0.30	0.30
Piezoelectric coupling coefficients $(C m^{- 2})$
$e_{31}^{} = e_{32}^{}$	–	–	25	0.046	9.6	22.86	18.46
Dielectric permittivity coefficient $(\times 1 0^{- 8})$ $(F m^{- 1})$
$k_{33}^{*}$	–	–	1.71	0.01062	1.710	0.0254	1.53

PVDF: polyvinylidene fluoride; PZT: Lead Zirconate Titanate.

Pinched hemisphere shell

In order to prove the applicability of the developed nonlinear finite element formulation for curved shell structures, a pinched hemisphere, with an 18° hole at its pole, is tested. The pinched hemisphere shell is subjected to purely mechanical loads. Alternating radial forces are applied at points A and B at 90° interval (two inward and two outward forces). The spherical geometry is defined by the radius R = 10 and the thickness of t = 0.04, the material properties are the Young’s modulus Y = 6.82510 × 10⁷, and the Poisson’s ratio ν = 0.3. The symmetry of the model allows the consideration of only one quarter of the shell for the finite element simulation using $16 \times 16$ elements, as depicted in Figure 4.

The analysis considers large displacements and rotations. The loads are increased by a factor of 400 to compare with results by Kim et al. (2008), Simo et al. (1990), and Buechter and Ramm (1992). Figure 5 depicts load-deflection curves at points A and B of the pinched hemispherical shell and Figure 6 represents the initial and deformed configurations for the hemispherical shells under the maximum load. As may be seen in Figure 5, results obtained by the current approach are in good agreement with those acquired by Kim et al. (2008) and Simo et al. (1990). The small difference between the present and the reference results (Buechter and Ramm, 1992) can be explained by the constant shear correction factor used in the reference model (Buechter and Ramm, 1992).

Bimorph beam

The standard test problem of cantilever beam, proposed by Tzou and Ye (1996), is invoked. The structure consists of two layers with opposite polarities, the length $L = 100 mm$ , the width $b = 5 mm$ , and the thickness $h = 1 mm$ as shown in Figure 7. Two different cases of actuation are considered here: The beam is considered once to be made of polyvinylidene fluoride (PVDF) and then again to be made of piezoelectric ceramic PIC151 in order to demonstrate the influence of the piezoelectric material factor. The material properties are described in Table 3.

The cantilevered beam is modeled using four shell elements of equal length and two discrete layers. The structure will act as an actuator with an applied voltage of 1 V across the thickness direction. When an unidirectional electric field is applied, one piezoelectric layer expands while the other contracts. This creates a bending moment which causes deflection of the beam. Comparisons of deflection predicted using the present element, finite element presented by Nestorović et al. (2012), and exact solution obtained using Bernoulli approach (Tzou, 1993) are depicted in Figure 7. The analytical solution of deflection is a quadratic function with respect to the coordinate along the length of the beam (x axis), as defined in work of Tzou (1993)

W (x) = \frac{3}{2} \frac{e_{31}^{*} Δ φ}{Y h^{2}} x^{2}

(77)

As may be seen, the present solutions are in good agreement with analytical and numerical ones for both materials (PVDF and PIC 151 materials). The PIC151 beam is deflected more than the PVDF beam as expected because PIC151 has a higher elastic material property and piezoelectric coupling coefficients compared to PVDF and its bending stiffness is, therefore, greater. This implies that piezoelectric material coefficients $e_{31}^{*}$ , $e_{32}^{*}$ , and $k_{33}^{*}$ have significant effect on piezoelastic linear behavior of the active structure. Such sensitive factors should be carefully considered for the design of smart structures: for effective actuation, it is desirable to choose piezoelectric material with higher piezoelectric coupling constants.

It is important now to investigate the sensitivity of the present element with respect of mesh distortion. Here, we consider the mesh distortion proposed by Sze and Pan (1999) for an actuator case, as shown in Figure 8. The PVDF bimorph beam was subjected to the same boundary conditions as in the previous test but with a distorted mesh. In fact, it was distorted in increments of 5 mm starting from e = 0–20 mm and the comparison of static deflection of PVDF cantilevered beam with a distorted mesh to an undistorted mesh is provided in Figure 9.

An analytical solution is proposed by Tzou (1993) according to an Euler–Bernoulli beam theory. As may be seen, the present results are in good agreement with the analytical solution even with more distorted mesh. Moreover, the developed model is more robust compared to the mixed hexahedral elements of Sze and Pan (1999) and Sze and Ghali (1993). The worst result is achieved for the element H8D used in Sze and Pan (1999). Here, the displacement, the electric potential, and the dielectric displacements were integrated into the element formulation as an additional independent variable. If, in addition, the stress σ is introduced as an independent variable, the element is named H8DS. These two elements are very sensitive to mesh distortion as a result of shear locking. Using a selective scaling technique in order to reduce the shear locking effect, an excellent agreement is achieved between results obtained by H8DS* element (Sze and Ghali, 1993) and the analytical solution. In this work, the present proposed element does not show any shear locking, thus the results fit the analytical solution even for strong distorted meshes (Figure 10).

Shape control of adaptive composite plate

A simply supported cross-ply plate made of Graphite/Epoxy T300/976, with the internal sequence of layers [0/90/0] _S is initially subjected to a uniformly distributed load of $200 N / m^{2}$ as shown in Figure 11. Two piezoelectric layers, made of piezoceramic PZT (Lead Zirconate Titanate), bonded to the top and bottom surfaces. Starting from 0 V, an increase in the electrical potential up to 27 V is supplied to the outer active layers with the aim of achieving the voltage giving the shape that corresponds mostly to the initial undeformed shape. The structure has dimensions $a \times b = 254 \times 254 mm$ ; the thickness of each graphite/epoxy layer is 0.138 mm and of active layer is 0.254 mm. The material properties used for the calculation are shown in Table 3.

The piezoelectric excitation is achieved by supplying the same voltage to the oppositely polarized piezolayers, which results in bending moments uniformly distributed over the edges of the plate. An 8 × 8 finite element mesh is applied. This test is originally examined by Kioua and Mirza (2000), who used the conventional Ritz analysis based on the shallow-shell theory, and recently by Marinković and Zehn (2015), who used a nine-node degenerated shell element based on the Reissner–Mindlin kinematical assumptions. The obtained results illustrated in Figure 11 for normalized centerline deflection are in good agreement with those given by Kioua and Mirza (2000) and Marinković and Zehn (2015) (Figure 12; Table 4). Hence, the present finite element was able to compute the deformed shape quite accurately, even though the solution showed a slight fluctuation for a 27 V input. In fact, the FEM results show that the structure subjected to the voltage of 27 V is not exactly flat. However, the structure recovers the flat shape accordingly to Ritz analysis (Kioua and Mirza, 2000).

Table 4.

Results shape control of adaptive composite plate.

W/b $(\times 1 0^{- 3})$	Present model	Marinković and Zehn (2015)	Kioua and Mirza (2000)
Applied voltage 0 V	−1.1656	−1.0466	−1.0712
Applied voltage 15 V	−0.4930	−0.4816	−0.4696
Applied voltage 27 V	−0.0485	−0.0548	−0.0000

Simply supported composite cylindrical arch

In this example, a composite cylindrical arch made of Graphite/Epoxy T300/976, with a stacking sequence [45/–45/0]s, is considered in order to analyze layered shell structures. The top and bottom surfaces are bordered by piezoelectric PZT G1195 layers as given in Figure 13. The system is simply supported at the rectilinear edges and is modeled by $4 \times 20$ elements with 20 elements in circumference direction and 4 elements for the width. The radius of the middle surface is R = 100 mm and the axial length is b = 60 mm. Each composite ply is of 0.12 mm thick, and the thickness of each piezoelectric layer is 0.24 mm. The material parameters for the graphite epoxy and PZT layers according to Marinković et al. (2006) are given in Table 3.

The piezoactive layers work as actuator and are excited by a potential voltage of 100 V. The radial deflection w is measured along the centerline at b / 2. The results, normalized by the total thickness h, are depicted in Figure 14 dependent on the normalized hoop distance between $0^{°}$ and $180^{°}$ of the cylinder segment. Comparison between the present shell formulation and data from Marinković et al. (2006) for a full biquadratic (nine-noded) formulation of the degenerate shell element with linear approximation for the electric potential in thickness direction are carried out. As may be seen from Figure 15, the results of developed model correlated to the solution of Marinković et al. (2006). The good accordance of the results shows the reliable applicability of the present laminated four-node formulation for layered piezomechanical structures.

Clamped composite-laminated plate with piezolayers (nonlinear analysis)

The purpose of this test is to apply the proposed shell element to the geometrically nonlinear static analysis of a smart structure under mechanical and electrical loading. The example consists of a clamped rectangular plate made of orthotropic fiber reinforced material (Graphite/Epoxy) with two piezolayers, made of PZT, bonded to its upper and lower surfaces, as shown in Figure 15. The lamination sequence is [30/90/30/90]s. The plate dimensions are 200 × 195 mm. The thickness of each composite layer and active layer is 0.2 mm and 0.1 mm, respectively. The material properties are listed in Table 3. The structure is discretized by 3 × 200 rectangular shell elements.

The plate is initially subjected to mechanical loading. A concentrated force is applied on nodes A, B, and C, which has a magnitude of 12 N. It is noted that the magnitude of loads at node B is twice that at nodes A and C. Table 5 lists the displacements of the same nodes at which the forces act obtained by present model in comparison with three-node shell element elaborated by Marinkovic and Rama (2017) and ABAQUS standard element S3. As can be seen, results obtained by the developed model are in accordance with those obtained by the reference element (Marinkovic and Rama, 2017). Moreover, in linear analysis, deflection at nodes A, B, and C is higher than in nonlinear response. This is true because nonlinear stiffness matrix $K_{G}$ has influence on the response of the plate. In fact, stiffness of the plate increases and subsequently deflection of the plate decreases from linear deflection.

Table 5.

Mechanical excitation of the clamped composite plate—free edge displacements.

Mechanical excitation		Point A	Point B	Point C
		Linear analysis
Displacement w (mm)	Present model	−55.50	−41.29	−25.65
	Abaqus S3	−50.93	−36.75	−22.90
	Marinković and Rama (2017)	−50.94	−36.80	−22.90
		Geometrically nonlinear analysis
	Present model	−46.15	−35.84	−26.27
	Abaqus S3	−44.81	−33.19	−23.55
	Marinković and Rama (2017)	−44.80	−33.22	−23.58

Table 6 presents the linear and nonlinear deflection at three nodes (A, B, and C) of clamped plate under electrical loading. An electrical potential of 1000 V is applied over the electrodes of the piezolayers (see Figure 16). Due to the opposite polarization of the active layers, their activation produces internal bending moments uniformly distributed over the edges of the piezoelectric layers. As explained in the mechanical excitation case, the nonlinear displacements are lower than displacements from linear response.

Table 6.

Voltage excitation of the clamped composite plate—free edge displacements.

Voltage excitation		Point A	Point B	Point C
		Linear analysis
Displacement w (mm)	Present model	93.32	49.62	21.95
	Abaqus S3	91.78	48.00	20.64
	Marinković and Rama (2017)	91.84	48.00	20.64
		Geometrically nonlinear analysis
	Present model	77.97	45.56	22.70
	Abaqus S3	83.15	45.89	21.72
	Marinković and Rama (2017)	83.59	46.30	22.10

Adaptive composite plate with surface-bonded piezolayers

A simply supported cross-ply square made of S-glass/Epoxy, with the internal sequence of layers [45°/–45°/45°] is subjected to mechanical and electrical loading. Two piezoelectric layers, made of PXE-52, bonded to the top and bottom surfaces (see Figure 17). The side dimension is $a = 0.1 m$ ; the thickness of each S-glass/Epoxy layer is 0.0004 m and of active layer is 0.0002 m. The material properties of S-glass/Epoxy used for the calculation are listed in Table 3.

Two load cases are considered, and linear and geometrically nonlinear computations are performed. In the first case, the plate is exposed to a uniform distributed load L_mech. Starting from p₀ = 10 kN/m², an increase in the mechanical load up to $30 kN / m^{2}$ is subjected to the plate. Since bending of the shell is induced, the vertical displacement W_c is at center point. The results obtained by the present model and the finite element developed by Moita et al. (2002) for both linear and geometrically nonlinear are summarized in Table 7. It may be noted that the deflection predicted by the CLT model (Moita et al., 2002) is slightly less than that of the proposed modified FSDT model, due to the higher stiffness of the former model as shear deformation is not accounted for.

Table 7.

Central deflection W_c (mm) for different load cases.

Load level µ	Mechanical load L_mech = µ p₀				Electrical load L_electr = µ V₀
	Linear analysis		Nonlinear analysis		Linear analysis		Nonlinear analysis
	(a)	(b)	(a)	(b)	(a)	(b)	(a)	(b)
0.5	0.1091	0.1085	0.1087	0.1089	−0.1065	−0.1035	−0.1059	−0.1036
1.0	0.2183	0.2170	0.2146	0.2172	−0.2131	−0.2070	−0.2085	−0.2064
1.5	0.3275	0.3255	0.3155	0.3237	−0.3196	−0.3105	−0.3053	−0.3073
2.0	0.4367	0.4340	0.4102	0.4265	−0.4262	−0.4140	−0.3949	−0.4055
2.5	0.5459	0.5425	0.4984	0.5257	−0.5327	−0.5175	−0.4769	−0.5004
3.0	0.6551	0.6510	0.5800	0.6208	−0.6393	−0.6210	−0.5518	−0.5916

(a) Present model; (b) Moita et al. (2002).

In the second case, an actuation of piezolayers by a constant electric voltage of V₀ = 151.35 V is considered. Therefore, the induced internal stresses result in a bending moment which causes deflection of the plate. Table 7 shows the load-level versus deflections at center point of the structure (linear and geometrically nonlinear analysis) to indicate the development of the induced deformation with the increasing electric voltage L_electr. As observed from the results, the effect of nonlinearity is more pronounced for the highest applied voltage. In fact, the difference between linear and nonlinear deflection decreases as the electric voltage becomes little.

Conclusion

A novel geometrically nonlinear piezoelectric finite shell element has been derived to statically simulate plate and shell structures with integrated sensors/actuators. This element, based on a modified FSDT theory, is developed to include stiffness and the electromechanical coupling of the active layers. Several numerical examples for different piezolaminated geometries are conducted to highlight the efficiency of the proposed model compared to existing solutions available in literature. An excellent agreement among the results confirms the high accuracy of the current piezoelastic model. In the case of thin plates (“Simply supported composite cylindrical arch” and “Clamped composite laminated plate with piezolayers (nonlinear analysis)” sections), there are differences between the CLT and the improved FSDT models and the FSDT and the improved FSDT models. This large variation in the results between these models is due to inappropriate estimation of the shear deformation in the FSDT and CLT models. The difference observed is considerable when the load is increasing. The proposed formulation can be used to examine various piezoelastic behavior of shells involving nonlinear dynamics, control, and free vibration of smart structures.

Footnotes

Appendix

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through research groups program under grant number (R.G.P.1/70/40).

ORCID iDs

Hanen Mallek

Hanen Jrad

References

Balamurugan

Narayanan

(2008) A piezolaminated composite degenerated shell finite element for active control of structures with distributed piezosensors and actuators. Smart Materials and Structures 17(3): 035031.

Bathe

Dvorkin

(1985) A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering 21: 367–383.

Benjeddou

Deu

Letombe

(2002) Free vibrations of simply-supported piezoelectric adaptive plates: an exact sandwich formulation. Thin-Walled Structures 40: 573–593.

Bodaghi

Damanpack

Aghdam

(2012) Non-linear active control of FG beams in thermal environments subjected to blast loads with integrated FGP sensor/actuator layers. Composite Structures 94(12): 3612–3623.

Bodaghi

Damanpack

Aghdam

(2014) Geometrically non-linear transient thermo-elastic response of FG beams integrated with a pair of FG piezoelectric sensors. Composite Structures 107: 48–59.

Bohlooly

Mirzavand

(2015) Closed form solutions for buckling and postbuckling analysis of imperfect laminated composite plates with piezoelectric actuators. Composites Part B: Engineering 72: 21–29.

Buechter

Ramm

(1992) Shell theory versus degeneration—a comparison in large rotation finite element analysis. International Journal for Numerical Methods in Engineering 34: 39–59.

Cinefra

Carrera

(2013) Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures. International Journal for Numerical Methods in Engineering 93: 160–182.

Cinefra

Valvano

(2016) A variable kinematic doubly-curved mitc9 shell element for the analysis of laminated composites. Mechanics of Advanced Materials and Structures 23(11):1312–1325.

10.

Cinefra

Belouettar

Soave

et al . (2010) Variable kinematic models applied to free vibration analysis of functionally graded materials shells. European Journal of Mechanics—A/Solids 29: 1078–1087.

11.

Cinefra

Carrera

Della Croce

et al . (2012) Refined shell elements for the analysis of functionally graded structures. Composite Structures 94: 415–422.

12.

Correia

IFP

Soares

CMM

Soares

CAM

et al . (2002) Active control of axisymmetric shells with piezoelectric layers: a mixed laminated theory with a high order displacement field. Computers and Structures 80: 2256–2275.

13.

Daneshmehr

Shakeri

(2007) Three-dimensional elasticity solution of cross-ply shallow and non-shallow panels with piezoelectric sensors under dynamic load. Composite Structures 80: 429–439.

14.

Fereidoon

Danial

Hessameddin

(2016) Dynamic behavior of piezoelectric composite beams under moving loads. Journal of Composite Materials 50: 899–916.

15.

Foda

Almajed

ElMadany

(2010) Vibration suppression of composite laminated beams using distributed piezoelectric patches. Smart Materials and Structures 19: 115018.

16.

Frikha

Wali

Hajlaoui

et al . (2016) Dynamic response of functionally graded material shells with a discrete double directors shell element. Composite Structures 154: 385–395.

17.

Frikha

Zghal

Dammak

(2018a) Dynamic analysis of functionally graded carbon nanotubes-reinforced plate and shell structures using a double directors finite shell element. Aerospace Science and Technology 78: 488–451.

18.

Frikha

Zghal

Dammak

(2018b) Finite rotation three and four nodes shell elements for functionally graded carbon nanotubes-reinforced thin composite shells analysis. Computer Methods in Applied Mechanics and Engineering 329: 289–311.

19.

Gabbert

Duvigneau

Ringwelski

(2017) Noise control of vehicle drive systems. Facta Universitatis, Series: Mechanical Engineering 15(2): 183–200.

20.

Gabbert

Koppe

Seeger

et al . (2002) Modeling of smart composite shell structures. Journal of Theoretical and Applied Mechanics 40(3): 575–593.

21.

Gandhi

Thompson

(1992) Smart Materials and Structures. London: Chapman and Hall.

22.

Hajlaoui

Jarraya

El Bikri

et al . (2015) Buckling analysis of functionally graded materials structures with enhanced solid-shell elements and transverse shear correction. Composite Structures 132: 87–97.

23.

Sivashanker

et al . (2001) Active control of FGM plates with integrated piezoelectric sensors and actuators. International Journal of Solids and Structures 38: 1641–1655.

24.

Irons

(1976) The semiloof shell element. In: Ashwell

Gallagher

(eds) Finite Elements for Thin Shells and Curved Membranes, Vol. 11. New York: John Wiley, pp. 11197–11222.

25.

Joshi

Mukherjee

Schmauder

(2003) Numerical characterization of functionally graded active materials under electrical and thermal fields. Smart Materials and Structures 12(4): 571–579.

26.

Jrad

Mallek

Wali

et al . (2018a) Finite element formulation for active functionally graded thin-walled structures. Comptes Rendus Mecanique 346: 1159–1178.

27.

Jrad

Mars

Wali

et al . (2018b) Geometrically nonlinear analysis of elastoplastic behavior of functionally graded shells. Engineering with Computers. Epub ahead of print 28 August 2018. DOI: 10.1007/s00366-018-0633-3.

28.

Kidane

Shukla

(2008) Dynamic constitutive behavior of Ti/TiB FGM under thermo-mechanical loading. Journal of Materials Science 43: 2771–2777.

29.

Kim

K-D

Lomboy

Han

S-C

(2008) Geometrically non-linear analysis of functionally graded material (FGM) plates and shells using a four-node quasi-conforming shell element. Journal of Composite Materials 42: 485–511.

30.

Kioua

Mirza

(2000) Piezoelectric induced bending and twisting of laminated composite shallow shells. Smart Materials and Structures 9: 476–484.

31.

Kogl

Bucalem

(2005) Analysis of smart laminates using piezoelectric MITC plate and shell elements. Computers and Structures 83: 1153–1163.

32.

Kulikov

Plotnikova

(2013) A new approach to three-dimensional exact solutions for functionally graded piezoelectric laminated plates. Composite Structures 106: 33–46.

33.

Kulikov

Mamontov

Plotnikova

(2015) Coupled thermoelectroelastic stress analysis of piezoelectric shells. Composite Structures 124: 65–76.

34.

Kulkarni

Bajoria

(2007) Large deformation analysis of piezolaminated smart structures using higher-order shear deformation theory. Smart Materials and Structures 16: 1506–1516.

35.

Lam

(1999) Active control of composite plates with integrated piezoelectric sensors and actuators under various dynamic loading conditions. Smart Materials and Structures 8(2): 223–237.

36.

Lee

(2005) Layerwise laminate analysis of functionally graded piezoelectric bimorph beams. Journal of Intelligent Material Systems and Structures 16(4): 365–371.

37.

Lentzen

Klosowski

Schmidt

(2007) Geometrically nonlinear finite element simulation of smart piezolaminated plates and shells. Smart Materials and Structures 16: 2265–2274.

38.

Levinson

(1980) An accurate, simple theory of the statics and dynamics of elastic plates. Mechanics Research Communications 7: 343–350.

39.

Liew

Kitipornchai

(2004) Finite element method for the feedback control of FGM shells in the frequency domain via piezoelectric sensors and actuators. Computer Methods in Applied Mechanics and Engineering 193: 257–273.

40.

Liew

et al . (2001) Active control of FGM plates subjected to a temperature gradient: modelling via finite element method based on FSDT. International Journal for Numerical Methods in Engineering 52: 1253–1271.

41.

Mallek

Jrad

Wali

Dammak

(2018) Piezoelastic response of smart functionally graded structure with integrated piezoelectric layers using discrete double directors shell element. Composite Structures 210: 354–366.

42.

Marinković

Petrović

(2012) On efficiency of a single-layer shell element for composite laminated structures. Mechanical Engineering 10(2): 105–112.

43.

Marinković

Rama

(2017) Co-rotational shell element for numerical analysis of laminated piezoelectric composite structures. Composites Part B: Engineering 125: 144–156.

44.

Marinkovic

Zehn

(2015) Finite element formulation for active composite laminates. American Journal of Engineering and Applied Sciences 8(3): 328–335.

45.

Marinković

Koppe

Gabbert

(2006) Numerically efficient finite element formulation for modeling active composite laminates. Mechanics of Advanced Materials and Structures 13: 379–392.

46.

Marinković

Koppe

Gabbert

(2009) Aspects of modeling piezoelectric active thin-walled structures. Journal of Intelligent Material Systems and Structures 20(15): 1835–1844.

47.

Mars

Koubaa

Wali

et al . (2017) Numerical analysis of geometrically non-linear behavior of functionally graded shells. Latin American Journal of Solids and Structures 14(11): 1952–1978.

48.

Mellouli

Jrad

Wali

et al . (2018) Meshless implementation of arbitrary 3D-shell structures based on a modified first order shear deformation theory. Computers & Mathematics with Applications. Epub ahead of print 11 October 2018. DOI: 10.1016/j.camwa.2018.09.010.

49.

Mindlin

(1951) Influence of rotary inertia and shear on flexural motions of isotropic elastic plates. Journal of Applied Mechanics 18: 31–38.

50.

Moita

JMS

Soares

CMM

Soares

CAM

(2002) Geometrically non-linear analysis of composite structures with integrated piezoelectric sensors and actuators. Composite Structures 57(1–4): 253–261.

51.

Müller

Drašar

Schilz

et al . (2003) Functionally graded materials for sensor and energy applications. Materials Science and Engineering: A 362: 17–39.

52.

Murthy

MVV

(1981) An improved transverse shear deformation theory for laminated anisotropic plates. NASA technical paper, pp. 1–37. Hampton, VA: Langley Research Center.

53.

Nestorović

Marinković

Chandrashekar

et al . (2012) Implementation of a user defined piezoelectric shell element for analysis of active structures. Finite Elements in Analysis and Design 52: 11–22.

54.

Nestorović

Marinković

Shabadi

et al . (2014) User defined finite element for modeling and analysis of active piezoelectric shell structures. Meccanica 49(8): 1763–1774.

55.

Neto

Leal

(2012) A triangular finite element with drilling degrees of freedom for static and dynamic analysis of smart laminated structures. Computers and Structures 108(2): 61–74.

56.

Liew

(2002) Finite element modeling of active control of functionally graded shells in frequency domain via piezoelectric sensors and actuators. Computational Mechanics 28: 1–9.

57.

Pompea

Worch

Epple

et al . (2003) Functionally graded materials for biomedical applications. Materials Science and Engineering: A 362: 40–60.

58.

Rama

(2017) A 3-node piezoelectric shell element for linear and geometrically nonlinear dynamic analysis of smart structures. Facta Universitatis, Series: Mechanical Engineering 15(1): 31–44.

59.

Reddy

(1984) A simple higher-order theory for laminated composite plates. Journal of Applied Mechanics 51: 745–752.

60.

Ren

Azar

(2006) A model for shape control of cross-ply laminated shells using a piezoelectric actuator. Journal of Composite Materials 40: 1271–1285.

61.

Schulz

Klinkel

(2008) A geometrically nonlinear mixed finite element formulation for the simulation of piezoelectric shell structures. Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 6932: 69320G.

62.

Shi

(2007) A new simple third-order shear deformation theory of plates. International Journal of Solids and Structures 44: 4399–4417.

63.

Simo

Fox

Rifai

(1990) Resultant geometrically exact shell model. Part iii: computational aspects of the nonlinear theory. Computer Methods in Applied Mechanics and Engineering 79: 21–70.

64.

Sudhakar

Kamal

(2003) Finite element modeling of smart plates/shells using higher order shear deformation theory. Composite Structures 62: 41–50.

65.

Sze

Ghali

(1993) Hybrid hexahedral element for solids, plates, shells and beams by selective scaling. International Journal for Numerical Methods in Engineering 36: 1519–1540.

66.

Sze

Pan

(1999) Hybrid finite element models for piezoelectric materials. Journal of Sound and Vibration 226(3): 519–547.

67.

Tanov

Tabiei

(2000) A simple correction to the first order shear deformation shell finite element formulations. Finite Elements in Analysis and Design 35: 189–197.

68.

Torres

DAF

Paulo de Tarso

(2010) HSDT-layerwise analytical solution for rectangular piezoelectric laminated plates. Composite Structures 92(8): 1763–1774.

69.

Tzou

(1993) Piezoelectric Shells: Distributed Sensing and Control of Continua. Dordrecht: Kluwer Academic.

70.

Tzou

(1996) Analysis of piezoelectric structures with laminated piezoelectric triangle shell elements. AIAA Journal 34(1): 110–115.

71.

Wali

Hentati

Jaraya

et al . (2015) Free vibration analysis of FGM shell structures with a discrete double directors shell element. Composite Structures 125: 295–303.

72.

Chen

Shen

et al . (2002) A high order theory for functionally graded piezoelectric shells. International Journal of Solids and Structures 39: 5325–5344.

73.

Xiang

Shi

(2008) Static analysis for multi-layered piezoelectric cantilevers. International Journal of Solids and Structures 45: 113–128.

74.

Zghal

Frikha

Dammak

(2018) Non-linear bending analysis of nanocomposites reinforced by graphene-nanotubes with finite shell element and membrane enhancement. Engineering Structures 158: 95–109.

75.

Zhang

Meng

Chen

(2007) Stability, nonlinearity and reliability of electrostatically actuated MEMS devices. Sensors 7(5): 760–796.

76.

Zhang

Tabata

Tsuchiya

et al . (2011) Noise-induced chaos in the electrostatically actuated MEMS resonators. Physics Letters A 375: 2903–2910.