Static and Dynamic Analysis of Piezoelectric Laminated Composite Beams and Plates

Abstract

In this chapter, the mechanical behavior analysis of piezoelectric laminated composite beams and plates is influenced subjected to static, dynamic, and aerodynamic loads. Algorithm for dynamic, stability problem analysis and vibration control of laminated composite beams and plates with piezoelectric layers is presented. In addition, numerical calculations, considering the effect of factors on static, dynamic, and stability response of piezoelectric laminated composite beams and plates are also clearly presented. The content of this chapter can equip readers with the knowledge used to calculate the static, dynamic, and vibration control of composite beams, panels made of piezoelectric layers applied in the field different techniques.

1. Introduction

The content of this chapter is the inheritance and development of the research results of the authors and other authors by published scientific works on composite materials, piezoelectric and structural calculation by piezoelectric composite materials.

2. Electromechanical interaction of piezoelectric materials

2.1 Mechanical-electrical behavior relations

Let us consider a block of elastic material in an environment with an electric field of zero, the relationship between stress and strain is followed Hooke’s law, and written as follows [1, 2]:

$$\left\{\mathrm{\sigma }\right\}=\left[\mathrm{c}\right]\left\{\mathrm{\epsilon }\right\},$$E1

where {σ} is the mechanical stress vector, {ε} is the mechanical strain vector, and [c] is the material stiffness matrix of beam.

Mechanical-electrical relations in piezoelectric materials have an interactive relationship, strain {ε} will produce eε - polarization, where e is the voltage stress factor when there is no mechanical strain. The imposed electric field E produces the -eE stress in the piezoelectric material according to the reverse voltage effect. Therefore, we have a mathematical model that describes the mechanical-electrical interaction relationship in piezoelectric materials as follows [3, 4, 5, 6]:

$$\left\{\mathrm{\sigma }\right\}=\left[\mathrm{c}\right]\left\{\mathrm{\epsilon }\right\}-\left[\mathrm{e}\right]\left\{\mathrm{E}\right\},$$E2

$$\left\{\mathrm{D}\right\}={\left[\mathrm{e}\right]}^{\mathrm{T}}\left\{\mathrm{\epsilon }\right\}+\left[\mathrm{p}\right]\left\{\mathrm{E}\right\},$$E3

$$\text{or}\left\{\mathrm{D}\right\}={\left[\mathrm{d}\right]}^{\mathrm{T}}\left\{\mathrm{\sigma }\right\}+\left[\mathrm{p}\right]\left\{\mathrm{E}\right\},$$E4

where [e] is the piezoelectric stress coefficient matrix, [p] is the dielectric constant matrix, {E} is the vector of applied electric field (V/m), and {D} is the vector of electric displacement (C/m²).

For the linear problem and small strain, strain vector in the piezoelectric structures can be defined as follows:

$$\left\{\mathrm{\epsilon }\right\}=\left[\mathrm{s}\right]\left\{\mathrm{\sigma }\right\}+\left[\mathrm{d}\right]\left\{\mathrm{E}\right\},$$E5

in which [s] is the matrix of compliance coefficients (m²/N), [d] is the matrix of piezoelectric strain constants (m/V).

In the field of engineering, piezoelectric materials are used by two types. The first type, the piezoelectric layers or the piezoelectric patches act as actuators, called the piezoelectric actuators. In this case, the piezoelectric layers are strained when imposing an electric field on it. The second type, the piezoelectric layers or piezoelectric patches act as sensors, called piezoelectric sensors. In this case, the voltage is generated in piezoelectric layers when there is mechanical strain.

2.2 Piezoelectric actuators and sensors

2.2.1 Piezoelectric actuators

Eq. (5) can be written in the matrix form as follows [4, 6]:

$$\left\{\mathrm{\epsilon }\right\}=\left\{\begin{array}{c}{\mathrm{\epsilon }}_{11}\\ {\mathrm{\epsilon }}_{22}\\ {\mathrm{\epsilon }}_{33}\\ {\mathrm{\gamma }}_{23}\\ {\mathrm{\gamma }}_{13}\\ {\mathrm{\gamma }}_{12}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathrm{s}}_{11}& {\mathrm{s}}_{12}& {\mathrm{s}}_{13}& {\mathrm{s}}_{14}& {\mathrm{s}}_{15}& {\mathrm{s}}_{16}\\ {\mathrm{s}}_{21}& {\mathrm{s}}_{22}& {\mathrm{s}}_{23}& {\mathrm{s}}_{24}& {\mathrm{s}}_{25}& {\mathrm{s}}_{26}\\ {\mathrm{s}}_{31}& {\mathrm{s}}_{32}& {\mathrm{s}}_{33}& {\mathrm{s}}_{34}& {\mathrm{s}}_{35}& {\mathrm{s}}_{36}\\ {\mathrm{s}}_{41}& {\mathrm{s}}_{42}& {\mathrm{s}}_{43}& {\mathrm{s}}_{44}& {\mathrm{s}}_{45}& {\mathrm{s}}_{46}\\ {\mathrm{s}}_{51}& {\mathrm{s}}_{52}& {\mathrm{s}}_{53}& {\mathrm{s}}_{54}& {\mathrm{s}}_{55}& {\mathrm{s}}_{56}\\ {\mathrm{s}}_{61}& {\mathrm{s}}_{62}& {\mathrm{s}}_{63}& {\mathrm{s}}_{64}& {\mathrm{s}}_{65}& {\mathrm{s}}_{66}\end{array}\right]\left\{\begin{array}{c}{\mathrm{\sigma }}_{11}\\ {\mathrm{\sigma }}_{22}\\ {\mathrm{\sigma }}_{33}\\ {\mathrm{\tau }}_{23}\\ {\mathrm{\tau }}_{13}\\ {\mathrm{\tau }}_{12}\end{array}\right\}+\left[\begin{array}{ccc}{\mathrm{d}}_{11}& {\mathrm{d}}_{21}& {\mathrm{d}}_{31}\\ {\mathrm{d}}_{12}& {\mathrm{d}}_{22}& {\mathrm{d}}_{32}\\ {\mathrm{d}}_{13}& {\mathrm{d}}_{23}& {\mathrm{d}}_{33}\\ {\mathrm{d}}_{14}& {\mathrm{d}}_{24}& {\mathrm{d}}_{34}\\ {\mathrm{d}}_{15}& {\mathrm{d}}_{25}& {\mathrm{d}}_{35}\\ {\mathrm{d}}_{16}& {\mathrm{d}}_{26}& {\mathrm{d}}_{36}\end{array}\right]\left\{\begin{array}{c}{\mathrm{E}}_1\\ {\mathrm{E}}_2\\ {\mathrm{E}}_3\end{array}\right\},$$E6

Assuming that the device is pulled along the axis 3, and viewing the piezoelectric material as a transversely isotropic material, which is true for piezoelectric ceramics, many of the parameters in the above matrices will be either zero, or can be expressed through each other. In particular, the non-zero compliance coefficients are s₁₁, s₁₂, s₁₃, s₂₁, s₂₂, s₂₃, s₃₁, s₃₂, s₃₃, s₄₄, s₅₅, s₆₆, in which s₁₂ = s₂₁, s₁₃ = s₃₁, s₂₃ = s₃₂, s₄₄ = s₅₅, s₆₆ = 2(s₁₁ − s₁₂).

Finally, Eq. (6) becomes:

$$\left\{\mathrm{\epsilon }\right\}=\left\{\begin{array}{c}{\mathrm{\epsilon }}_{11}\\ {\mathrm{\epsilon }}_{22}\\ {\mathrm{\epsilon }}_{33}\\ {\mathrm{\gamma }}_{23}\\ {\mathrm{\gamma }}_{13}\\ {\mathrm{\gamma }}_{12}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathrm{s}}_{11}& {\mathrm{s}}_{12}& {\mathrm{s}}_{13}& 0& 0& 0\\ {\mathrm{s}}_{12}& {\mathrm{s}}_{22}& {\mathrm{s}}_{23}& 0& 0& 0\\ {\mathrm{s}}_{13}& {\mathrm{s}}_{23}& {\mathrm{s}}_{33}& 0& 0& 0\\ 0& 0& 0& {\mathrm{s}}_{44}& 0& 0\\ 0& 0& 0& 0& {\mathrm{s}}_{55}& 0\\ 0& 0& 0& 0& 0& {\mathrm{s}}_{66}\end{array}\right]\left\{\begin{array}{c}{\mathrm{\sigma }}_{11}\\ {\mathrm{\sigma }}_{22}\\ {\mathrm{\sigma }}_{33}\\ {\mathrm{\tau }}_{23}\\ {\mathrm{\tau }}_{13}\\ {\mathrm{\tau }}_{12}\end{array}\right\}+\left[\begin{array}{ccc}0& 0& {\mathrm{d}}_{31}\\ 0& 0& {\mathrm{d}}_{32}\\ 0& 0& {\mathrm{d}}_{33}\\ 0& {\mathrm{d}}_{24}& 0\\ {\mathrm{d}}_{15}& 0& 0\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\mathrm{E}}_1\\ {\mathrm{E}}_2\\ {\mathrm{E}}_3\end{array}\right\},$$E7

where E₁, E₂, and E₃ are electric fields in the 1, 2, and 3 directions, respectively.

2.2.2 Piezoelectric sensors

The induction charge equation of piezoelectric sensor layers is derived from Eq. (4) can be written in the matrix form as [4, 6, 7]:

$$\left\{\mathrm{D}\right\}=\left\{\begin{array}{c}{\mathrm{D}}_1\\ {\mathrm{D}}_2\\ {\mathrm{D}}_3\end{array}\right\}=\left[\begin{array}{cccccc}{\mathrm{d}}_{11}& {\mathrm{d}}_{12}& {\mathrm{d}}_{13}& {\mathrm{d}}_{14}& {\mathrm{d}}_{15}& {\mathrm{d}}_{16}\\ {\mathrm{d}}_{21}& {\mathrm{d}}_{22}& {\mathrm{d}}_{23}& {\mathrm{d}}_{24}& {\mathrm{d}}_{25}& {\mathrm{d}}_{26}\\ {\mathrm{d}}_{31}& {\mathrm{d}}_{32}& {\mathrm{d}}_{33}& {\mathrm{d}}_{34}& {\mathrm{d}}_{35}& {\mathrm{d}}_{36}\end{array}\right]\left\{\begin{array}{c}{\mathrm{\sigma }}_{11}\\ {\mathrm{\sigma }}_{22}\\ {\mathrm{\sigma }}_{33}\\ {\mathrm{\tau }}_{23}\\ {\mathrm{\tau }}_{13}\\ {\mathrm{\tau }}_{12}\end{array}\right\}+\left[\begin{array}{ccc}{\mathrm{p}}_{11}& {\mathrm{p}}_{12}& {\mathrm{p}}_{13}\\ {\mathrm{p}}_{21}& {\mathrm{p}}_{22}& {\mathrm{p}}_{23}\\ {\mathrm{p}}_{31}& {\mathrm{p}}_{32}& {\mathrm{p}}_{33}\end{array}\right]\left\{\begin{array}{c}{\mathrm{E}}_1\\ {\mathrm{E}}_2\\ {\mathrm{E}}_3\end{array}\right\},$$E8

The non-zero piezoelectric strain constants are d₃₁, d₃₂, d₁₅, d₂₄, and d₃₃, in which d₃₁ = d₃₂, d₁₅ = d₂₄. And the non-zero dielectric coefficients are p₁₁, p₂₂, and p₃₃, where p₁₁ = p₂₂. Eq. (8) becomes:

$$\left\{\mathrm{D}\right\}=\left\{\begin{array}{c}{\mathrm{D}}_1\\ {\mathrm{D}}_2\\ {\mathrm{D}}_3\end{array}\right\}=\left[\begin{array}{cccccc}0& 0& 0& 0& {\mathrm{d}}_{15}& 0\\ 0& 0& 0& {\mathrm{d}}_{24}& 0& 0\\ {\mathrm{d}}_{31}& {\mathrm{d}}_{32}& {\mathrm{d}}_{33}& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\mathrm{\sigma }}_{11}\\ {\mathrm{\sigma }}_{22}\\ {\mathrm{\sigma }}_{33}\\ {\mathrm{\tau }}_{23}\\ {\mathrm{\tau }}_{13}\\ {\mathrm{\tau }}_{12}\end{array}\right\}+\left[\begin{array}{ccc}{\mathrm{p}}_{11}& 0& 0\\ 0& {\mathrm{p}}_{22}& 0\\ 0& 0& {\mathrm{p}}_{33}\end{array}\right]\left\{\begin{array}{c}{\mathrm{E}}_1\\ {\mathrm{E}}_2\\ {\mathrm{E}}_3\end{array}\right\},$$E9

where D₁, D₂, D₃, p₁₁, p₂₂, and p₃₃ are the displacement charge, dielectric constant in the 1, 2, and 3 directions, respectively.

Normally, the voltage is transmitted through the thickness of the actuator layers.

3. Static and dynamic analysis of laminated composite beams with piezoelectric layers

3.1 Displacement and strain

Based on the first-order shear deformation theory (FSDT), the displacement field at any point of the beam is defined as [1, 2]:

$$\begin{array}{l}\mathrm{u}\left(\mathrm{x},\mathrm{z}\right)={\mathrm{u}}_0\left(\mathrm{x}\right)+z{\mathrm{\theta }}_y\left(\mathrm{x}\right),\\ \mathrm{w}\left(\mathrm{x},\mathrm{z}\right)={\mathrm{w}}_0\left(\mathrm{x}\right),\end{array}$$E10

where u, w denotes the displacements of a point (x, z) in the beam; u₀, w₀ are the displacements of a point at the beam neutral axis, and θ_y is the rotation of the transverse normal about the y axis. The bending and shear strains associated with the displacement field in Eq. (10) are defined as:

$$\left\{\mathrm{\epsilon }\right\}=\left\{\begin{array}{c}{\mathrm{\epsilon }}_x\\ {\mathrm{\gamma }}_{\mathit{xz}}\end{array}\right\}=\left\{\begin{array}{c}\frac{\mathrm{du}}{\mathrm{dx}}\\ \frac{\mathrm{du}}{\mathrm{dz}}+\frac{\mathrm{dw}}{\mathrm{dx}}\end{array}\right\}=\left\{\begin{array}{c}\frac{{\mathrm{du}}_0}{\mathrm{dx}}+\mathrm{z}\frac{{\mathrm{d\theta }}_z}{\mathrm{dx}}\\ {\mathrm{\theta }}_y+\frac{\mathrm{d}{\mathrm{w}}_0}{\mathrm{dx}}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{dx}}& 0& \mathrm{z}\frac{\mathrm{d}}{\mathrm{dx}}\\ 0& \frac{\mathrm{d}}{\mathrm{dx}}& 1\end{array}\right]\left\{\begin{array}{c}{\mathrm{u}}_0\\ {\mathrm{w}}_0\\ {\mathrm{\theta }}_z\end{array}\right\},$$E11

in which ε_x, γ_xz are the normal strain, and shear strain, respectively.

Using finite element method, we consider 2-node bending elements with 3 degrees of freedom per node (Figure 1).

Figure 1.

Two noded beam element.

The displacements of the beam neutral axis are expressed in local coordinate system in the form:

$$\left\{{\mathrm{d}}_0\right\}=\left\{\begin{array}{c}{\mathrm{u}}_0\\ {\mathrm{v}}_0\\ {\mathrm{\theta }}_z\end{array}\right\}=\left\{\begin{array}{c}\left[{\mathrm{N}}^{\mathrm{u}}\right]\left\{{\mathrm{q}}^{\mathrm{u}}\right\}\\ \left[{\mathrm{N}}^{\mathrm{v}}\right]\left\{{\mathrm{q}}^{\mathrm{v}}\right\}\\ \left[{\mathrm{N}}^{{\mathrm{\theta }}_z}\right]\left\{{\mathrm{q}}^{{\mathrm{\theta }}_z}\right\}\end{array}\right\}=\left[{\mathrm{N}}^{\mathrm{M}}\right]{\left\{{\mathrm{q}}_{\mathrm{b}}\right\}}_{\mathrm{e}},$$E12

where {q_b}_e is the vector of vector of nodal displacements of element, [N^M] is the matrix mechanical shape functions:

$${\left\{{\mathrm{q}}_{\mathrm{b}}\right\}}_{\mathrm{e}}={\left\{{\mathrm{q}}_1{\mathrm{q}}_2{\mathrm{q}}_3{\mathrm{q}}_4{\mathrm{q}}_5{\mathrm{q}}_6\right\}}^{\mathrm{T}},$$E13

$$\underset{3\times 6}{\underbrace{\left[{\mathrm{N}}^{\mathrm{M}}\right]}}=\left[\begin{array}{ccc}\left[{\mathrm{N}}^{\mathrm{u}}\right]& 0& 0\\ 0& \left[{\mathrm{N}}^{\mathrm{v}}\right]& 0\\ 0& 0& \left[{\mathrm{N}}^{{\mathrm{\theta }}_{\mathrm{z}}}\right]\end{array}\right],$$E14

in which [N^u], [N^v], [N^θz] are, in this order, the row vectors of longitudinal, transverse along y, and rotation about z shape functions.

Substituting Eq. (12) into Eq. (11), we obtain:

$$\underset{2\times 1}{\underbrace{\left\{\mathrm{\epsilon }\right\}}}=\underset{2\times 6}{\underbrace{\left[{\mathrm{B}}_{\mathrm{b}}\right]}}\underset{6\times 1}{\underbrace{{\left\{{\mathrm{q}}_{\mathrm{b}}\right\}}_{\mathrm{e}}}},$$E15

$$\text{where}\left[{\mathrm{B}}_{\mathrm{b}}\right]=\left[\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{dx}}& 0& z\frac{\mathrm{d}}{\mathrm{dx}}\\ 0& \frac{\mathrm{d}}{\mathrm{dx}}& 1\end{array}\right]\left[{\mathrm{N}}^{\mathrm{M}}\right].$$E16

The electric potential is constant over the element surface:

$${\varphi }_k=\sum _{i=1}^n{N}_i{\varphi }_i,$$E17

where n is the element node number.

A voltage ϕ is applied across an actuator of layer thickness t_p generates an electric field vector {E}, such that [4, 8, 9, 10]:

$$\left\{E_k\right\}=-\nabla {\varphi }_k=\left\{\begin{array}{ccc}0& 0& E_k^z\end{array}\right\},$$E18

in which

$$E_k^z=-\frac{{\varphi }_k}{t_{\mathit{pk}}}=\left[B_{\varphi }\right]\left\{\varphi \right\}={\left[\begin{array}{cccccc}0& 0& \frac{1}{t_{p1}}& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{1}{t_{p2}}\end{array}\right]}^T\left\{\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right\},$$E19

where t_pk is the thickness of the k^th piezoelectric layer.

Substituting Eq. (19) into Eq. (18), the electric field vector {E} can also be defined in terms of nodal variables as:

$$\left\{\mathrm{E}\right\}=-\left[{\mathrm{B}}_{\mathrm{\varphi }}\right]{\left\{\mathrm{\varphi }\right\}}_{\mathrm{e}},$$E20

Using Eqs. (15), and (20), the linear piezoelectric constitutive equations coupling the elastic and electric fields will be completely determined by Eqs. (2) and (3).

3.2 Finite element equations

Using Hamilton’s principle, we have [11, 12, 13]:

$$\underset{t_1}{\overset{t_2}{\int }}\left(T^e-U^e-W^e\right)\mathit{dt}=0,$$E21

where T^e, U^e are the kinetic and potential energy, respectively and W^e is the work done by external forces. They are determined by:

$$T^e=\frac{1}{2}\underset{V_e}{\int }\rho {\left\{\dot{q}\right\}}_e^T{\left\{q\right\}}_e\mathit{dV},$$E22

$$U^e=\frac{1}{2}\underset{V_e}{\int }{\left\{\epsilon \right\}}_e^T{\left\{\sigma \right\}}_e\mathit{dV},$$E23

$$W^e=\underset{V_e}{\int }{\left\{q\right\}}_e^T{\left\{f_b\right\}}_e\mathit{dV}+\underset{S_e}{\int }{\left\{q\right\}}_e^T{\left\{f_s\right\}}_e\mathit{dS}+{\left\{q\right\}}_e^T{\left\{f_c\right\}}_e,$$E24

in which ${\left\{f_b\right\}}_e,{\left\{f_s\right\}}_e,{\left\{f_c\right\}}_e$ are the body, surface, and concentrated forces acting on the element, respectively. V_e and S_e are elemental volume and area.

Substituting Eqs. (15), (2), (20), (22), (23), and (24) into Eq. (21), one obtains:

$$\left[M_{\mathit{bb}}^e\right]{\left\{\ddot{q}\right\}}_e+\left[K_{\mathit{bb}}^e\right]{\left\{q\right\}}_e+\left[K_{\mathit{b\varphi }}^e\right]{\left\{\varphi \right\}}_e={\left\{f\right\}}_e,$$E25

$$\left[K_{\mathit{\varphi b}}^e\right]{\left\{q\right\}}_e-\left[K_{\mathit{\varphi \varphi }}^e\right]{\left\{\varphi \right\}}_e={\left\{Q\right\}}_e,$$E26

where

$$\text{Element mass matrix}:\left[M_{\mathit{bb}}^e\right]=\underset{V_e}{\int }\rho {\left[N^M\right]}^T\left[N^M\right]\mathit{dV},$$E27

$$\text{Element mechanical stiffness matrix}:\left[K_{\mathit{bb}}^e\right]=\underset{S_e}{\int }{\left[B_b\right]}^T\left[H\right]\left[B_b\right]\mathit{dS},$$E28

Element mechanical-electrical coupling stiffness matrix:

$$\left[K_{\mathit{b\varphi }}^e\right]=\underset{S_e}{\int }{\left[B_b\right]}^T\left[\overline{e}\right]\left[B_{\varphi }\right]\mathit{dS},$$E29

Element electrical-mechanical coupling stiffness matrix:

$$\left[K_{\mathit{\varphi b}}^e\right]={\left[K_{\mathit{b\varphi }}^e\right]}^T,$$E30

Element piezoelectric permittivity matrix:

$$\left[K_{\mathit{\varphi \varphi }}^e\right]=-\underset{S_e}{\int }{\left[B_{\varphi }\right]}^T\left[\overline{p}\right]\left[B_{\varphi }\right]\mathit{dS},$$E31

$$\text{where}\left[H\right]=\left[\begin{array}{cc}{c}_{11}& 0\\ 0& c_{22}\end{array}\right],\left[\overline{e}\right]=\left[\begin{array}{cc}{e}_{11}& e_{12}\\ e_{21}& e_{22}\end{array}\right],\left[\overline{p}\right]=\left[\begin{array}{cc}{t}_{p1}{p}_{11}& 0\\ 0& t_{p2}{p}_{22}\end{array}\right],$$E32

{f}_e, {Q}_e are the applied external load and charge, respectively.

3.2.1 Static analysis

In the case of beams subjected to static loads, zero acceleration, from Eqs. (25) and (26), we obtain the static equations of the beam as follows:

$$\left[K_{\mathit{bb}}^e\right]{\left\{q\right\}}_e+\left[K_{\mathit{b\varphi }}^e\right]{\left\{\varphi \right\}}_e={\left\{f\right\}}_e,$$E33

$$\left[K_{\mathit{\varphi b}}^e\right]{\left\{q\right\}}_e-\left[K_{\mathit{\varphi \varphi }}^e\right]{\left\{\varphi \right\}}_e={\left\{Q\right\}}_e,$$E34

Assembling the element equations yields general static equation:

$$\left[K_{\mathit{bb}}\right]\left\{q\right\}+\left[K_{\mathit{b\varphi }}\right]\left\{\varphi \right\}=\left\{f\right\},$$E35

$$\left[K_{\mathit{\varphi b}}\right]\left\{q\right\}-\left[K_{\mathit{\varphi \varphi }}\right]\left\{\varphi \right\}=\left\{Q\right\}.$$E36

where [K_bb], [K_ϕϕ] are the overall mechanical stiffness and piezoelectric permittivity matrices respectively; [K_bϕ] and [K_ϕb] are the overall mechanical - electrical and electrical - mechanical coupling stiffness matrices, respectively, and {q}, {ϕ} are respectively the overall mechanical displacement, and electric potential vector.

Substituting Eq. (36) into Eq. (35) yields:

$$\left(\left[K_{\mathit{bb}}\right]+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\right)\left\{q\right\}=\left\{f\right\}+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left\{Q\right\},$$E37

Substituting {q} from Eq. (37) into Eq. (36), we obtain the vector {ϕ}.

3.2.2 Dynamic analysis

From Eqs. (25) and (26), assembling the element equations yields general dynamic equation of motion:

$$\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left[K_{\mathit{bb}}\right]\left\{q\right\}+\left[K_{\mathit{b\varphi }}\right]\left\{\varphi \right\}=\left\{f\right\},$$E38

$$\left[K_{\mathit{\varphi b}}\right]\left\{q\right\}-\left[K_{\mathit{\varphi \varphi }}\right]\left\{\varphi \right\}=\left\{Q\right\},$$E39

Substituting Eq. (39) into Eq. (38), we obtain:

$$\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left(\left[K_{\mathit{bb}}\right]+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\right)\left\{q\right\}=\left\{f\right\}+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left\{Q\right\},$$E40

3.2.3 Free vibration analysis

For free vibrations, from Eq. (40), the governing equation is:

$$\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left(\left[K_{\mathit{bb}}\right]+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\right)\left\{q\right\}=\left\{0\right\}.$$E41

The beam vibrations induce charges and electric potentials in sensor layers. Therefore, the control system allows current to flow and feeds back to the actuators. In this case, if we apply no external charge Q to a sensor, from Eq. (39), we will have:

$${\left[K_{\mathit{\varphi \varphi }}\right]}_s^{-1}{\left[K_{\mathit{\varphi b}}\right]}_s{\left\{q\right\}}_s={\left\{\varphi \right\}}_s.$$E42

and ${\left\{Q\right\}}_s={\left[K_{\mathit{\varphi b}}\right]}_s{\left\{q\right\}}_s$ is the induced charge due to strain.

The operation of the amplified control loop implies, the actuating voltage is determined by the following relationship [1, 10, 14]:

$${\left\{\varphi \right\}}_a=G_d{\left\{\varphi \right\}}_s+G_v{\left\{\dot{\varphi }\right\}}_s,$$E43

where G_d and G_v are the feedback control gains for displacement and velocity.

Substituting Eq. (43) into Eq. (39), the charge in the actuator due to actuator strain in response to the beam vibration modified by control system feedback is:

$${\left[K_{\mathit{\varphi b}}\right]}_a{\left\{q\right\}}_a-{\left[K_{\mathit{\varphi \varphi }}\right]}_a\left(G_d{\left\{\varphi \right\}}_s+G_v{\left\{\dot{\varphi }\right\}}_s\right)={\left\{Q\right\}}_a.$$E44

Substituting (42) into (44) leads to:

$${\left\{Q\right\}}_a={\left[K_{\mathit{\varphi b}}\right]}_a{\left\{q\right\}}_a-G_d{\left[K_{\mathit{\varphi \varphi }}\right]}_a{\left[K_{\mathit{\varphi \varphi }}\right]}_s^{-1}{\left[K_{\mathit{\varphi b}}\right]}_s{\left\{q\right\}}_s-G_v{\left[K_{\mathit{\varphi \varphi }}\right]}_a{\left[K_{\mathit{\varphi \varphi }}\right]}_s^{-1}{\left[K_{\mathit{\varphi b}}\right]}_s{\left\{\dot{q}\right\}}_s.$$E45

Substituting Eq. (45) into (40), we obtain:

$$\begin{array}{l}\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left(\left[K_{\mathit{bb}}\right]+\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\right)\left\{q\right\}=\left\{f\right\}+\\ +\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left(\begin{array}{l}{\left[K_{\mathit{\varphi b}}\right]}_a{\left\{q\right\}}_a-G_v{\left[K_{\mathit{\varphi \varphi }}\right]}_a{\left[K_{\mathit{\varphi \varphi }}\right]}_s^{-1}{\left[K_{\mathit{\varphi b}}\right]}_s{\left\{\dot{q}\right\}}_s-\\ -G_d{\left[K_{\mathit{\varphi \varphi }}\right]}_a{\left[K_{\mathit{\varphi \varphi }}\right]}_s^{-1}{\left[K_{\mathit{\varphi b}}\right]}_s{\left\{q\right\}}_s\end{array}\right),\end{array}$$E46

in which {q}_s ≡ {q}_a ≡ {q} is the beam displacement vector, [K_ϕϕ]_a = [K_ϕϕ]_s = [K_ϕϕ] is the piezoelectric permittivity matrix, and [K_ϕb]_a = [K_ϕb]_s = [K_ϕb] is the mechanical-electrical coupling stiffness matrix.

Therefore, Eq. (46) becomes:

$$\begin{array}{l}\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left[K_{\mathit{bb}}\right]\left\{q\right\}+G_v\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi \varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\left\{\dot{q}\right\}+\\ +G_d\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi \varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\left\{q\right\}=\left\{f\right\}.\end{array}$$E47

In the case of considering the structural damping, the equation of motion of the beam is:

$$\left[M_{\mathit{bb}}\right]\left\{\ddot{q}\right\}+\left(\left[C_A\right]+\left[C_R\right]\right)\left\{\dot{q}\right\}+K^{\ast }\left\{q\right\}=\left\{f\right\},$$E48

where $\left[C_A\right]=G_v\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi \varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]$ is the active damping matrix, $\left[K^{\ast }\right]=\left(\left[K_{\mathit{bb}}\right]+G_d\left[K_{\mathit{b\varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi \varphi }}\right]{\left[K_{\mathit{\varphi \varphi }}\right]}^{-1}\left[K_{\mathit{\varphi b}}\right]\right)$ is the total of mechanical stiffness matrix and piezoelectric, $\left[C_R\right]={\alpha }_R\left[M_{\mathit{bb}}\right]+{\beta }_R\left[K_{\mathit{bb}}\right]$ is the overall structural damping matrix, α_R, and β_R are respectively the Rayleigh damping coefficients, which are generally determined by the first and second natural frequencies (ω₁, ω₂) and ratio of damping ξ, {f} is the overall mechanical force vector.

Eq. (48) can be solved by the direct integration Newmark’s method.

3.3 Numerical analysis

An example for free vibration of laminated beam affected by piezoelectric layers is presented here. The beam is made of four layers symmetrically (0°/90°/90°/0°) of epoxy-T300/976 graphite material with 2.5 mm thickness per layer, and with one layer piezo ceramic materials bonded to the top and bottom surfaces, 2.0 mm thickness per layer as shown in Figure 2 is considered (a = 0.254 m, b = 0.0254 m). The material properties of the piezo ceramic layers and graphite-epoxy are shown in Table 1.

Figure 2.

Piezoelectric composite cantilever beam.

Properties	PZT G1195 N	T300/976
E11 [N/cm^2]	0.63 × 10^6	1.50 × 10^6
E22 = E33 [N/cm^2]	0.63 × 10^6	0.09 × 10^6
ν12 = ν13 = ν23	0.3	0.3
G12 = G13 [N/cm^2]	0.242 × 10^6	0.071 × 10^6
G23 [N/cm^2]	0.242 × 10^6	0.025 × 10^6
ρ [kg/m^3]	7600	1600
d31 = d32 (m/V)	254 × 10^−12	—
p11 = p22 (F/m)	15.3 × 10^−9	—
p33 (F/m)	15.0 × 10^−9	—

Table 1.

Relevant mechanical properties of respective materials.

The direct integration Newmark’s method is used with parameters α_R = 0.5, β_R = 0.25; integral time step Δt = 0.005 s with total time calculated t = 15 s.

Figures 3 and 4 illustrate the vertical displacement w at the free end of the beam for two cases:

1. Case 1: With structural damping, and without piezoelectric damping (G_v = 0, G_d = 0).

2. Case 2: With structural damping, with piezoelectric damping (G_v = 0.5, G_d = 30).

Figure 3.

Vertical displacement response (G_v = 0, G_d = 0 − Case 1).

Figure 4.

Vertical displacement response (G_v = 0.5, G_d = 30 − Case 2).

4. Dynamic analysis of laminated piezoelectric composite plates

4.1 The electromechanical behavioral relations in the plate

Consider laminated composite plates with general coordinate system (x, y, z), in which the x, y plane coincides with the neutral plane of the plate. The top and bottom surfaces of the plate are bonded to the piezoelectric patches or piezoelectric layers (actuator and sensor). The plate under the load acting on its neutral plane has any temporal variation rule (Figure 5).

Figure 5.

Piezoelectric composite plate and coordinate system of the plate (a), and lamina details (b).

Hypothesis: The piezoelectric composite plate corresponds with Reissner-Mindlin theory. The material layers are arranged symmetrically through the neutral plane of the plate, ideally adhesive with each other.

4.1.1 Strain - displacement relations

Based on the first-order shear deformation theory, the displacement fields at any point in the plate are [7, 8]:

$$\begin{array}{l}u\left(x,y,z,t\right)=u_0\left(x,y,t\right)+z{\theta }_y\left(x,y,t\right),\\ v\left(x,y,z,t\right)=v_0\left(x,y,t\right)-z{\theta }_x\left(x,y,t\right),\\ w\left(x,y,z,t\right)=w_0\left(x,y,t\right),\end{array}$$E49

where u, v and w are the displacements of a general point (x, y, z) in the laminate along x, y and z directions, respectively. u₀, v₀, w₀, θ_x and θ_y are the displacements and rotations of a midplane transverse normal about the y-and x-axes respectively.

The components of the strain vector corresponding to the displacement field (49) are defined as:

For the linear strain:

$$\begin{array}{l}{\epsilon }_x=\frac{\partial u}{\partial x}=\frac{\partial u_0}{\partial x}+z\frac{\partial {\theta }_y}{\partial x},{\epsilon }_y=\frac{\partial v}{\partial y}=\frac{\partial v_0}{\partial y}-z\frac{\partial {\theta }_x}{\partial y},\\ {\gamma }_{\mathrm{xy}}=\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)+\frac{\mathrm{\partial w}}{\partial x}\cdot \frac{\mathrm{\partial w}}{\partial y}=\left(\frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\right)+z\left(\frac{\partial {\theta }_y}{\partial x}-\frac{\partial {\theta }_x}{\partial y}\right),\\ {\gamma }_{\mathrm{xz}}=\frac{\mathrm{\partial u}}{\partial z}+\frac{\mathrm{\partial w}}{\partial x}=\frac{\partial {\mathrm{w}}_0}{\partial x}+{\theta }_y,{\gamma }_{\mathrm{yz}}=\frac{\mathrm{\partial v}}{\partial z}+\frac{\mathrm{\partial w}}{\partial y}=\frac{\partial {\mathrm{w}}_0}{\partial y}-{\theta }_x,\end{array}$$E50

or in the vector form:

$$\begin{array}{l}\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{\mathit{xy}}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_x^o\\ {\epsilon }_y^o\\ {\gamma }_{\mathit{xy}}^o\end{array}\right\}+z\left\{\begin{array}{c}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{\mathit{xy}}\end{array}\right\}=\left[\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]\left\{\begin{array}{c}{u}_0\\ v_0\end{array}\right\}+z\left[\begin{array}{cc}-\frac{\partial }{\partial y}& 0\\ 0& -\frac{\partial }{\partial x}\\ -\frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]\left\{\begin{array}{c}{\theta }_x\\ {\theta }_y\end{array}\right\}=\\ =\left[D_{\epsilon }\right]\left\{\begin{array}{c}{u}_0\\ v_0\end{array}\right\}+\left[D_{\kappa }\right]\left\{\begin{array}{c}{\theta }_x\\ {\theta }_y\end{array}\right\}=\left\{{\epsilon }_0\right\}+z\left\{\kappa \right\}=\left\{{\epsilon }_b^L\right\},\end{array}$$E51

$$\left\{\begin{array}{c}{\gamma }_{\mathit{xz}}\\ {\gamma }_{\mathit{yz}}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\partial }{\partial x}& 0& 1\\ \frac{\partial }{\partial y}& -1& 0\end{array}\right]\left\{\begin{array}{c}{w}^o\\ {\theta }_x\\ {\theta }_y\end{array}\right\}=\left[\begin{array}{cc}\left\{{}^{w}D\right\}& -\left[I_s\right]\end{array}\right]\left\{\begin{array}{c}{w}_0\\ {\theta }_x\\ {\theta }_y\end{array}\right\}=\left\{{\epsilon }_s\right\}.$$E52

and for the nonlinear strain:

$$\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{\mathit{xy}}\end{array}\right\}=\left\{{\epsilon }_b^L\right\}+\left\{{\epsilon }^N\right\}=\left\{{\epsilon }_b^N\right\},$$E53

$$\left\{\begin{array}{c}{\gamma }_{\mathit{xz}}\\ {\gamma }_{\mathit{yz}}\end{array}\right\}=\left\{{\epsilon }_s\right\},$$E54

where $\left\{{\epsilon }^N\right\}=\frac{1}{2}\left[\begin{array}{cc}\frac{\partial w_0}{\partial x}& 0\\ 0& \frac{\partial w_0}{\partial y}\\ \frac{\partial w_0}{\partial y}& \frac{\partial w_0}{\partial x}\end{array}\right]\left\{\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right\}{w}_0$ is the non-linear strain vector, $\left\{{\epsilon }_b^L\right\}$ is the linear strain vector, {ε_s} is the shear strain vector.

4.1.2 Stress-strain relations

The equation system describing the stress-strain relations and mechanical-electrical quantities is respectively written as [8, 14]:

$$\begin{array}{l}\left\{{\sigma }_b\right\}=\left[Q\right]\left\{{\epsilon }_b^N\right\}-\left[e\right]\left\{E\right\},\\ \left\{{\tau }_b\right\}=\left[Q_s\right]\left\{{\epsilon }_s\right\},\end{array}$$E55

$$\left\{D\right\}=\left[e\right]\left\{{\epsilon }_b^N\right\}+\left[p\right]\left\{E\right\},$$E56

where $\left\{{\sigma }_b\right\}={\left\{\begin{array}{ccc}{\sigma }_x& {\sigma }_y& {\tau }_{\mathit{xy}}\end{array}\right\}}^T$ is the plane stress vector, $\left\{{\tau }_b\right\}={\left\{\begin{array}{cc}{\tau }_{\mathit{yz}}& {\tau }_{\mathit{xz}}\end{array}\right\}}^T$ is the shear stress vector, [Q] is the ply in-plane stiffness coefficient matrix in the structural coordinate system, [Q_s] is the ply out-of-plane shear stiffness coefficient matrix in the structural coordinate system. Notice that {τ_b} is free from piezoelectric effects.

The in-plane force vector at the state pre-buckling:

$$\left\{N^0\right\}={\left\{\begin{array}{ccc}{N}_x^0& N_y^0& N_{\mathit{xy}}^0\end{array}\right\}}^T=\sum _{k=1}^n\underset{h_{k-1}}{\overset{h_k}{\int }}{\left\{\begin{array}{c}{\sigma }_x^0\\ {\sigma }_y^0\\ {\tau }_{\mathit{xy}}^0\end{array}\right\}}_k\mathit{dz}.$$E57

4.1.3 Total potential energy

The total potential energy of the system is given by:

$$\mathrm{\Pi }=\frac{1}{2}\underset{V_p}{\int }{\left\{{\epsilon }_b^N\right\}}^T\left\{{\sigma }_b\right\}\mathit{dV}+\frac{1}{2}\underset{V_p}{\int }{\left\{{\epsilon }_s\right\}}^T\left\{{\tau }_b\right\}\mathit{dV}-\frac{1}{2}\underset{V_p}{\int }{\left\{E\right\}}^T\left\{D\right\}\mathit{dV}-W,$$E58

where W is the energy of external forces, V_p is the entire domain including composite and piezoelectric materials.

Introducing [A], [B], [D], [A_s], and vectors {N_p}, {M_p} as [8]:

$$\begin{array}{l}\left(\left[A\right],\left[B\right],\left[D\right]\right)=\underset{-h/2}{\overset{h/2}{\int }}\left(1,z,z^2\right)\left[Q\right]\mathit{dz},\\ \left[A_s\right]=\underset{-h/2}{\overset{h/2}{\int }}\left[Q_s\right]\mathit{dz},\left(\left\{N_p\right\},\left\{M_p\right\}\right)=\underset{-h/2}{\overset{h/2}{\int }}\left(1,z\right)\left[e\right]\left\{E\right\}\mathit{dz},\end{array}$$E59

where h is the total laminated thickness and combining with (5), (6) the total potential energy equation (8) can be written

$$\begin{array}{l}\mathrm{\Pi }=\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{\left\{{\epsilon }_0\right\}}^T\left[A\right]\left\{{\epsilon }_0\right\}d\mathrm{\Omega }+\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{\left\{\kappa \right\}}^T\left[D\right]\left\{\kappa \right\}d\mathrm{\Omega }+\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{\left\{{\epsilon }_s\right\}}^T\left[A_s\right]\left\{{\epsilon }_s\right\}d\mathrm{\Omega }+\\ +\underset{\mathrm{\Omega }}{\int }{\left\{{\epsilon }^N\right\}}^T\left(\left[A\right]\left\{{\epsilon }_0\right\}-\left[N_p\right]\right)d\mathrm{\Omega }-\underset{\mathrm{\Omega }}{\int }{\left\{{\epsilon }_0\right\}}^T\left[N_p\right]d\mathrm{\Omega }-\underset{\mathrm{\Omega }}{\int }{\left\{\kappa \right\}}^T\left[M_p\right]d\mathrm{\Omega }-W,\end{array}$$E60

where Ω is the plane xy domain of the plate.

4.2 Dynamic stability analysis of laminated composite plate with piezoelectric layers

4.2.1 Finite element models

Nine-node Lagrangian finite elements are used with the displacement and strain fields represented by Eqs. (49), (53), and (54). In the developed models, there is one electric potential degree of freedom for each piezoelectric layer to represent the piezoelectric behavior and thus the vector of electrical degrees of freedom is [6, 14]:

$$\left\{{\varphi }^e\right\}={\left\{\begin{array}{ccccc}.& .& {\varphi }_j^e& .& .\end{array}\right\}}^T,j=1,\dots ,{\mathit{NPL}}^e,$$E61

in which NPL^e is the number of piezoelectric layers in a given element.

The vector of degrees of freedom for the element {q^e} is:

$$\left\{q^e\right\}={\left\{\begin{array}{ccccc}\left\{q_1^e\right\}& \left\{q_2^e\right\}& \dots & \left\{q_9^e\right\}& {\varphi }^e\end{array}\right\}}^T,$$E62

where $\left\{q_i^e\right\}=\left\{\begin{array}{ccccc}{u}_i& v_i& w_i& {\theta }_{x_i}& {\theta }_{y_i}\end{array}\right\}$ is the mechanical displacement vector for node i.

4.2.2 Dynamic equations

The dynamic equations of piezoelectric composite plate can be derived by using Hamilton’s principle, accordingly, the vibration equation of the membrane (without damping) with in-plane loads is:

$$\left[M_{\mathit{ss}}\right]\left\{{\ddot{q}}_{\mathit{ss}}\right\}+\left[K_{\mathit{ss}}\right]\left\{q_{\mathit{ss}}\right\}=\left\{F\left(t\right)\right\}.$$E63

The equation of bending vibrations with out-of-plane loads is:

$$\left[\begin{array}{cc}\left[M_{\mathit{bb}}\right]& \left[0\right]\\ \left[0\right]& \left[0\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\ddot{q}}_{\mathit{bb}}\right\}\\ \left\{\ddot{\varphi }\right\}\end{array}\right\}+\left[\begin{array}{cc}\left[C_R\right]& \left[0\right]\\ \left[0\right]& \left[0\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\dot{q}}_{\mathit{bb}}\right\}\\ \left\{\dot{\varphi }\right\}\end{array}\right\}+\left[\begin{array}{cc}\left[K_{\mathit{bb}}\right]+\left[K_G\right]& \left[K_{\mathit{b\varphi }}\right]\\ \left[K_{\mathit{\varphi b}}\right]& -\left[K_{\mathit{\varphi \varphi }}\right]\end{array}\right]\left\{\begin{array}{c}\left\{q_{\mathit{bb}}\right\}\\ \left\{\varphi \right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{R\right\}\\ \left\{Q_{\mathit{el}}\right\}\end{array}\right\},$$E64

where [M_ss], [K_ss] are the overall mass, membrane elastic stiffness matrix respectively, and $\left\{q_{\mathit{ss}}\right\},\left\{{\dot{q}}_{\mathit{ss}}\right\},\left\{{\ddot{q}}_{\mathit{ss}}\right\}$ are respectively the membrane displacement, velocity, acceleration vector. [M_bb], [K_bb] and $\left\{q_{\mathit{bb}}\right\},\left\{{\dot{q}}_{\mathit{bb}}\right\},\left\{{\ddot{q}}_{\mathit{bb}}\right\}$ are the overall mass, bending elastic stiffness matrix and the bending displacement, velocity, acceleration vector; [K_G] is the overall geometric stiffness matrix; ([K_G] is a function of external in-plane loads); {F(t)} is the in-plane load vector, {R} is the normal load vector, {Q_el} is the vector containing the nodal charges and in-balance charges.

The element coefficient matrices are:

$$\left[K_G^e\right]=\left[K_{\mathit{Gx}}^e\right]+\left[K_{\mathit{Gy}}^e\right]+\left[K_{\mathit{Gxy}}^e\right],$$E65

$$\text{where}\begin{array}{l}\left[K_{\mathit{Gx}}^e\right]=\underset{A_e}{\int }{N}_x^0\left[N_x^{\text{'}}\right]{\left[N_x^{\text{'}}\right]}^T{\mathit{dA}}_e,\\ \left[K_{\mathit{Gy}}^e\right]=\underset{A_e}{\int }{N}_y^0\left[N_y^{\text{'}}\right]{\left[N_y^{\text{'}}\right]}^T{\mathit{dA}}_e,\\ \left[K_{\mathit{Gxy}}^e\right]=\underset{A_e}{\int }{N}_{\mathit{xy}}^0\left[N_x^{\text{'}}\right]{\left[N_y^{\text{'}}\right]}^T{\mathit{dA}}_e,\end{array}$$E66

$$\text{in which}\left[N_x^{\text{'}}\right]=\frac{\partial }{\partial x}\left[N\left(x,y\right)\right],\left[N_y^{\text{'}}\right]=\frac{\partial }{\partial y}\left[N\left(x,y\right)\right],$$E67

$$\frac{\mathrm{\partial w}}{\mathrm{\partial x}}=\left[\frac{\mathrm{\partial N}}{\mathrm{\partial x}}\right]\left\{q_{\mathit{bb}}^e\right\}=\left[N_x^{\text{'}}\right]\left\{q_{\mathit{bb}}^e\right\},\frac{\mathrm{\partial w}}{\partial y}=\left[\frac{\mathrm{\partial N}}{\partial y}\right]\left\{q_{\mathit{bb}}^e\right\}=\left[N_y^{\text{'}}\right]\left\{q_{\mathit{bb}}^e\right\}$$E68

$$\left[K_G\right]=\sum _{\mathit{ne}}\left[K_G^e\right]$$E69

4.2.3 Dynamic stability analysis

When the plate is subjected to in-plane loads only ({R} = {0}), the in-plane stresses can lead to buckling, from Eqs. (63) and (64) the governing differential equations of motion of the damped system may be written as:

$$\begin{array}{l}\left[M_{\mathit{ss}}\right]\left\{{\ddot{q}}_{\mathit{ss}}\right\}+\left[K_{\mathit{ss}}\right]\left\{q_{\mathit{ss}}\right\}=\left\{F\left(t\right)\right\},\\ \left[M_{\mathit{bb}}\right]\left\{{\ddot{q}}_{\mathit{bb}}\right\}+\left[C_R\right]\left\{{\dot{q}}_{\mathit{bb}}\right\}+\left(\left[K_{\mathit{bb}}\right]+\left[K_G\right]\right)\left\{q_{\mathit{bb}}\right\}+\left[K_{\mathit{b\varphi }}\right]\left\{\varphi \right\}=\left\{0\right\},\\ \left[K_{\mathit{\varphi b}}\right]\left\{q_{\mathit{bb}}\right\}-\left[K_{\mathit{\varphi \varphi }}\right]\left\{\varphi \right\}=\left\{Q_{\mathit{el}}\right\}.\end{array}$$E70

Eq. (70) is rewritten as:

$$\begin{array}{l}\left[M_{\mathit{ss}}\right]\left\{{\ddot{q}}_{\mathit{ss}}\right\}+\left[K_{\mathit{ss}}\right]\left\{q_{\mathit{ss}}\right\}=\left\{F\left(t\right)\right\},\\ \left[M_{\mathit{bb}}\right]\left\{{\ddot{q}}_{\mathit{bb}}\right\}+\left(\left[C_A\right]+\left[C_R\right]\right)\left\{{\dot{q}}_{\mathit{bb}}\right\}+\left(\left[K^{\ast }\right]+\left[K_G\right]\right)\left\{q_{\mathit{bb}}\right\}=\left\{0\right\}.\end{array}$$E71

The overall geometric stiffness matrix [K_G] is defined as follows:

• In the case of only tensile or compression plates (w = 0): Solving Eq. (71) helps us to present unknown displacement vector {q_ss}, and then stress vector:

$$\left\{{\sigma }_{\mathit{ss}}\right\}=\left[A_s\right]\left[B_s\right]\left\{q_{\mathit{ss}}\right\},$$E72

where [A_s] and [B_s] are the stiffness coefficient matrix and strain-displacement matrix of the plane problem.

• In the case of bending plate (w ≠ 0), the stress vector is:

$$\begin{array}{l}\left\{{\sigma }_{\mathit{sb}}\right\}=\left\{{\sigma }_{\mathit{ss}}\right\}+\left\{{\sigma }_{\mathit{bb}}\right\},\\ \left\{{\sigma }_{\mathit{bb}}\right\}=\left[A_b\right]\left[B_b\right]\left\{q_{\mathit{bb}}\right\},\end{array}$$E73

where [A_b] and [B_s] are the stiffness coefficient matrix and strain-displacement matrix of the plane bending problem.

Stability criteria [14]:

• In the case of plate subjected to periodic in-plane loads and without damping, the elastic stability problems become simple only by solving the linear equations to determine the eigenvalues.

• In case of the plate under any in-plane dynamic load and with damping, the elastic stability problems become very complex. This iterative method can be proved effectively and the following dynamic stability criteria are used:

  • Plate is considered to be stable if the maximum bending deflection is three times smaller than the plate’s thickness:Eq. (71)has the solution (w_i)_max satisfying the condition$0\le {\left|{\mathrm{w}}_{\mathrm{i}}\right|}_{\max }<3h$, where w_i is the deflection of the plate at node number i.

  • Plate is called to be in critical status if the maximum bending deflection of the plate is three times equal to the plate’s thickness.Eq. (71)has the solution (w_i)_max satisfying the condition${\left|{\mathrm{w}}_{\mathrm{i}}\right|}_{\max }=3h$.

  • Plate is called to be at buckling if the maximum deflection of the plate is three times larger than the plate’s thickness:Eq. (71)has the solution (w_i)_max satisfying the condition${\left|{\mathrm{w}}_{\mathrm{i}}\right|}_{\max }>3h$.

The identification of critical forces is carried out by the iterative method.

4.2.4 Iterative algorithm

1. Step 1. Defining the matrices, the external load vector and errors of load iterations.

2. Step 2. Solving Eq. (71) to present unknown displacement vector, {q_ss} and the stress vector is defined by (72), updating the geometric stiffness matrix [K_G].

3. Step 3. Solving Eq. (71) to present unknown bending displacement vector {q_bb}, and then testing stability conditions.

4. − If for all $\left|{\mathrm{w}}_i\right|=0$: increase load, recalculate from step 2;

5. − If at least one value $\left|{\mathrm{w}}_i\right|\ne 0$:

6. + In case: $0<{\left|{\mathrm{w}}_i\right|}_{m\mathrm{ax}}<3h$: Define stress vector by Eq. (73), update the geometric stiffness matrix [K_G]. Increase load, recalculate from step 2;

7. + In case: $0\le \frac{\left|{\left|{\mathrm{w}}_i\right|}_{\max }-3h\right|}{{\left|{\mathrm{w}}_i\right|}_{\max }}\le {\epsilon }_D$: Critical load p = p_cr. End.

4.2.5 Numerical analysis

Stability analysis of piezoelectric composite plate with dimensions a × b × h, where a = 0.25 m, b = 0.30 m, h = 0.002 m. Piezoelectric composite plate is composed of three layers, in which two layers of piezoelectric PZT-5A at its top and bottom are considered, each layer thickness h_p = 0.00075 m; the middle layer material is Graphite/Epoxy material, with thickness h₁ = 0.0005 m. The material properties for graphite/epoxy and PZT-5A are shown in Section 5.1 above. One short edge of the plate is clamped, the other three edges are free. The in-plane half-sine load is evenly distributed on the short edge of the plate: p(t) = p₀sin(2πft), where p₀ is the amplitude of load, f = 1/T = 1/0.01 = 100 Hz (0 ≤ t ≤ T/2 = 0.005 s) is the excitation frequency, voltage applied V = 50 V. The iterative error of the load ε_D = 0.02% is chosen.

Consider two cases: with damping (ξ = 0.05, G_v = 0.5, G_d = 15) and without damping (ξ = 0.0, G_v = 0.0, G_d = 15). The response of vertical displacement at the plate centroid over the plate thickness for the two cases is shown in Figure 6.

Figure 6.

Vertical displacement response at the plate centroid over the plate thickness.

The results show that the critical load of the plate with damping is larger than that without damping. In the two cases above, the critical load rises by 6.8%.

Analyze the stability of the plate with damping when a voltage of −200, −150, −100, −50, 0, 50, 100, 150 and 200 V is applied to the actuator layer of the piezoelectric composite plate.

Figure 7 shows the time history of the vertical displacement at the plate centroid over the plate thickness when a voltage of 0, 50, 100, 150 and 200 V is applied. The relation between critical load and voltages is shown in Figure 8.

Figure 7.

Vertical displacement response at the plate centroid over the plate thickness.

Figure 8.

Critical load-voltage relation.

The results show that the voltage applied to the piezoelectric layers affects the stability of the plate. As the voltage increases, the critical load of the plate also increases.

When the amplitude of the load changes from 0.25p_cr to 1.5p_cr (where p_cr is the amplitude of the critical load), a voltage of 50 V is applied to the actuator layer of the plate.

The results show the time history response of the vertical displacement at the plate centroid over the plate thickness as seen in Figure 9.

Figure 9.

Time history of the vertical displacement at the plate centroid over the plate thickness when p₀ = 0.25p_cr, 0.5p_cr, 0.75p_cr, 1.0p_cr, 1.25p_cr, and 1.5p_cr.

4.3 Dynamic analysis of piezoelectric stiffened composite plates subjected to airflow

Consider isoparametric piezoelectric laminated stiffened plate with the general coordinate system (x, y, z), in which the x, y plane coincides with the neutral plane of the plate. The top surface and lower surface of the plate are bonded to the piezoelectric patches (actuator and sensor). The plate subjected to the airflow load acting (Figure 10).

Figure 10.

Smart stiffened plate subjected to airflow. (a) Smart stiffened plate and coordinate system and (b) Lamina details.

The dynamic equations of a finite smart composite plate are written as follows:

$${\left[M\right]}_e{\left\{\ddot{u}\right\}}_e+{\left[C_A\right]}_e{\left\{\dot{u}\right\}}_e+\left({\left[K_{\mathit{bb}}^{\ln }\right]}_e+{\left[K_{\mathit{bb}}^{\mathrm{n}l}\right]}_e+{\left[K_A\right]}_e\right){\left\{u\right\}}_e={\left\{f\right\}}_e^m,$$E74

where ${\left[K_{\mathit{bb}}^{\ln }\right]}_e=\underset{V_e}{\int }{\left[B_b^{\ln }\right]}^T\left[Q\right]\left[B_b^{\ln }\right]\mathit{dV},$ and ${\left[K_{\mathit{bb}}^{\mathit{nl}}\right]}_e=\underset{V_e}{\int }{\left[B_b^{\mathit{nl}}\right]}^T\left[Q\right]\left[B_b^{\mathit{nl}}\right]\mathit{dV}$ are the element linear mechanical stiffness and nonlinear mechanical stiffness respectively, ${\left\{f\right\}}_e^m$ is element external mechanical force vector.

4.3.1 Formulation of Stiffener:

4.3.1.1 Formulation of x-Stiffener

$$\begin{array}{l}{U}_{\mathit{xs}}\left(x,z\right)=u_0\left(x\right)+z{\theta }_{\mathit{xs}}\left(x\right),\\ W_{\mathit{xs}}\left(x,z\right)=w_{\mathit{xs}}\left(x\right).\end{array}$$E75

where x-axis is taken along the stiffener centerline and the z-axis is its upward normal. The plate and stiffener element shown in Figure 11.

Figure 11.

Modeling of plate and stiffener element.

If we consider that the x-stiffener is attached to the lower side of the plate, conditions of displacement compatibility along their line of connection can be written as:

$${\left.u_p\right|}_{z=-t_p/2}={\left.u_{\mathit{xs}}\right|}_{z=t_{\mathit{xs}}/2},{\left.{\theta }_{\mathit{xp}}\right|}_{z=-t_p/2}={\left.{\theta }_{\mathit{xs}}\right|}_{z=t_{\mathit{xs}}/2},{\left.w_p\right|}_{z=-t_p/2}={\left.w_{\mathit{xs}}\right|}_{z=t_{\mathit{xs}}/2},$$E76

where t_p is the plate thickness and t_xs is the x-stiffener depth.

The element stiffness and mass matrices are defined as follows [2, 15]:

$${\left[K_{\mathit{xs}}\right]}_e=\underset{l_e}{\int }{\left[B_{\mathit{xs}}\right]}^T\left[D_{\mathit{xs}}\right]\left[B_{\mathit{xs}}\right]\mathit{dx},$$E77

$${\left[M_{\mathit{xs}}\right]}_e=\underset{A_e}{\int }\left[P\left({\left[N_{u^0}\right]}^T\left[N_{u^0}\right]+{\left[N_w\right]}^T\left[N_w\right]\right)+I_y\left({\left[N_{{\theta }_x}\right]}^T\left[N_{{\theta }_x}\right]\right)\right]\mathit{dA},$$E78

with [B_xs] is the strain-displacement relations matrix, [D_xs] is the stress-strain relations matrix and l_e is the element length, $\left[N_{u^0}\right],\left[N_w\right]\mathit{\text{and}}\left[N_{{\theta }_x}\right]$ are the shape function matrices relating the primary variables u₀, w, x, in terms of nodal unknowns, I_y is the area moment of inertia related to the y-axis and $P={\sum }_{k=1}^n\underset{h_{k-1}}{\overset{h_k}{\int }}{\rho }_k\mathit{dz},$ with ρ_k is density of k^th layer.

4.3.1.2 Formulation of y-Stiffener

The same as for x-stiffener, the element stiffness and mass matrices of the y-stiffener are defined as follows:

$${\left[K_{\mathit{ys}}\right]}_e=\underset{l_e}{\int }{\left[B_{\mathit{ys}}\right]}^T\left[D_{\mathit{ys}}\right]\left[B_{\mathit{ys}}\right]\mathit{dy},$$E79

$${\left[M_{\mathit{ys}}\right]}_e=\underset{A_e}{\int }\left[P\left({\left[N_{u^0}\right]}^T\left[N_{u^0}\right]+{\left[N_w\right]}^T\left[N_w\right]\right)+I_x\left({\left[N_{{\theta }_y}\right]}^T\left[N_{{\theta }_y}\right]\right)\right]\mathit{dA},$$E80

4.3.2 Modeling the effect of aerodynamic pressure and motion equations of the smart composite plate-stiffeners element

Based on the first order theory, the aerodynamic pressure l_h and moment m_θ, can be described as [15, 16, 17]:

$$\begin{array}{l}{l}_w=\frac{1}{2}{\rho }_a{\left(Ucos\alpha \right)}^{2}B\left[{\mathit{kH}}_1^{\ast }\frac{\dot{\mathrm{w}}}{Ucos\alpha }+{\mathit{kH}}_2^{\ast }\frac{B\dot{\theta }}{Ucos\alpha }+k^2{H}_3^{\ast }\theta \right]+\frac{1}{2}{C}_p{\rho }_{\mathrm{a}}{\left(Usin\alpha \right)}^2,\\ m_{\theta }=\frac{1}{2}{\rho }_a{\left(Ucos\alpha \right)}^2{B}^2\left[{\mathit{kA}}_1^{\ast }\frac{\dot{\mathrm{w}}}{Ucos\alpha }+{\mathit{kA}}_2^{\ast }\frac{B\dot{\theta }}{Ucos\alpha }+k^2{A}_3^{\ast }\theta \right],\end{array}$$E81

where $k=\mathit{b\omega }/U$ is defined as the reduced frequency, ω is the circular frequency of oscillation of the airfoil, U is the wind velocity, B is the half-chord length of the airfoil or half-width of the plate, ρ_a is the air density and α is the angle of attack.

The functions $A_i^{\ast }\left(K\right),H_i^{\ast }\left(K\right)$are defined as follows:

$$\begin{array}{l}{H}_1^{\ast }\left(K\right)=-\frac{\pi }{k}F\left(k\right),H_2^{\ast }\left(K\right)=-\frac{\pi }{4k}\left[1+F\left(k\right)+\frac{2G\left(k\right)}{k}\right],\\ H_3^{\ast }\left(K\right)=-\frac{\pi }{2k^2}\left[F\left(k\right)-\frac{\mathit{kG}\left(k\right)}{2}\right],A_1^{\ast }\left(K\right)=\frac{\pi }{4k}F\left(k\right),\\ A_2^{\ast }\left(K\right)=-\frac{\pi }{16k}\left[1-F\left(k\right)-\frac{2G\left(k\right)}{k}\right],A_3^{\ast }\left(K\right)=\frac{\pi }{8k^2}\left[\frac{k^2}{8}+F\left(k\right)-\frac{\mathit{kG}\left(k\right)}{2}\right],\end{array}$$E82

where F(k) and G(k) are defined as:

$$\begin{array}{l}F\left(k\right)=\frac{0.500502k^3+0.512607k^2+0.2104k+0.021573}{k^3+1.035378k^2+0.251293k+0.021508},\\ G\left(k\right)=-\frac{0.000146k^3+0.122397k^2+0.327214k+0.001995}{k^3+2.481481k^2+0.93453k+0.089318}.\end{array}$$E83

Using finite element method, aerodynamic force vector can be described as:

$${\left\{f\right\}}_e^{\mathit{air}}=-{\left[K^{\mathit{air}}\right]}_e{\left\{u\right\}}_e-{\left[C^{\mathit{air}}\right]}_e{\left\{\dot{u}\right\}}_e+{\left\{f\right\}}_e^n,$$E84

with ${\left[K^{\mathit{air}}\right]}_e,{\left[C^{\mathit{air}}\right]}_e$ and ${\left\{f\right\}}_e^n$ are the aerodynamic stiffness, damping matrices and lift force vector, respectively

$$\left[K_e^{\mathit{air}}\right]={\rho }_a{\left(Ucos\alpha \right)}^2{\mathit{Bk}}^2\underset{A_e}{\int }\left[H_3^{\ast }\left(k\right){\left[N_w\right]}^T\left[N_{\mathit{\theta x}}\right]+{\mathit{BA}}_3^{\ast }\left(k\right){\left[\frac{\partial N_{\mathit{\theta y}}}{\partial x}\right]}^T\left[N_{\mathit{\theta x}}\right]\right]\mathit{dA},$$E85

$$\left[C_e^{\mathit{air}}\right]={\rho }_a\left(\mathit{Uc}\mathrm{os}\alpha \right)\mathit{Bk}\left[\begin{array}{l}\underset{A_e}{\int }\left(H_1^{\ast }\left(k\right){\left[N_{\mathrm{w}}\right]}^T\left[N_{\mathrm{w}}\right]+{\mathit{BH}}_2^{\ast }\left(k\right){\left[N_{\mathrm{w}}\right]}^T\left[N_{\theta \mathrm{x}}\right]\mathit{dA}\right)\\ +\underset{A_e}{\int }\left({\mathit{BA}}_1^{\ast }\left(k\right){\left[\frac{\partial N_{\mathit{\theta y}}}{\partial x}\right]}^T\left[N_{\mathrm{w}}\right]+B^2{A}_2^{\ast }\left(k\right){\left[\frac{\partial N_{\mathit{\theta y}}}{\partial x}\right]}^T\left[N_{\theta \mathrm{x}}\right]\mathit{dA}\right)\end{array}\right],$$E86

$${\left\{f\right\}}_e^n=C_p{\rho }_{\mathrm{a}}{\left(Usin\alpha \right)}^2\underset{A_e}{\int }{\left[N_w\right]}^T\mathit{dA},$$E87

where A_e is the element area, [N_w], [N_θ] are the shape functions.

From Eqs. (74) and (84), the governing equations of motion of the smart composite plate-stiffeners element subjected to an aerodynamic force without damping can be derived as:

$${\left[M^{\ast }\right]}_e{\left\{\ddot{u}\right\}}_e+{\left[C_A\right]}_e{\left\{\dot{u}\right\}}_e+\left({\left[K^{\ast }\right]}_e+{\left[K_A\right]}_e+{\left[K^{\mathit{air}}\right]}_e\right){\left\{u\right\}}_e={\left\{f^{\ast }\right\}}_e^m,$$E88

where ${\left[M^{\ast }\right]}_e={\left[M\right]}_e+{\left[M_{\mathit{xs}}\right]}_e+{\left[M_{\mathit{ys}}\right]}_e$, ${\left[K^{\ast }\right]}_e={\left[K_{\mathit{bb}}^{\ln }\right]}_e+{\left[K_{\mathit{bb}}^{\mathrm{n}l}\right]}_e+{\left[K_{\mathit{xs}}\right]}_e+{\left[K_{\mathit{ys}}\right]}_e,$${\left\{f^{\ast }\right\}}_e^m={\left\{f\right\}}_e^m+{\left\{f\right\}}_e^n$.

4.3.3 Governing differential equations for total system

Finally, the elemental equations of motion are assembled to obtain the open-loop global equation of motion of the overall stiffened composite plate with the PZT patches as follows:

$$\left[M^{\ast }\right]\left\{\ddot{u}\right\}+\left(\left[C_R\right]+\left[C_A\right]\right)\left\{\dot{u}\right\}+\left(\left[K^{\ast }\right]+\left[K_A\right]+\left[K^{\mathit{air}}\right]\right)\left\{u\right\}={\left\{f^{\ast }\right\}}^m,$$E89

where $\left[C_R\right]={\alpha }_R\left[M_{\mathit{bb}}\right]+{\beta }_R\left(\left[K_{\mathit{bb}}^{\ln }\right]+\left[K_{\mathit{bb}}^{\mathrm{n}l}\right]\right)$.

The solution of nonlinear Eq. (89) is carried out by using Newmark direct and Newton-Raphson iteration method.

4.3.4 Numerical applications

A rectangle cantilever laminated composite plate is assumed to be [0°/90°]_s with total thickness 4 mm, length of 600 mm and width of 400 mm with three stiffeners along each direction x and y. The geometrical dimension of the stiffener is 5 mm of high and 10 mm of width. The plate and stiffeners are made of graphite/epoxy with mechanical properties: E₁₁ = 181 GPa, E₂₂ = E₃₃ = 10.3 GPa, E₁₂ = 7.17 GPa, ν₁₂ = 0.35, ν₂₃ = ν₃₂ = 0.38, ρ = 1600 kg·m⁻³. Material properties for piezoelectric layer made of PZT-5A are: d₃₁ = d₃₂ = −171 × 10⁻¹² m/V, d₃₃ = 374 × 10⁻¹² m/V, d₁₅ = d₂₄ = −584 × 10⁻¹² m/V, G₁₂ = 7.17 GPa, G₂₃ = 2.87 GPa, G₃₂ = 7.17 GPa, ν_PZT = 0.3, ρ_PZT = 7600 kg·m⁻³ and thickness t_PZT = 0.15876 mm, ξ = 0.05, G_v = 0.5, G_d = 15. The effects of the excitation frequency and location of the actuators are presented through a parametric study to examine the vibration shape of the composite plate activated by the surface bonded piezoelectric actuators. The iterative error of the load ε_D = 0.02% is chosen. The piezoelectric stiffened composite plate is subjected to the airflow in the positive x direction as shown in Figure 10a.

Dynamic response of the piezoelectric stiffened composite plate is shown in Figure 12.

Figure 12.

History of the plate at a critical airflow velocity U_cr = 30.5 m/s. (a) Displacement response and (b) Piezoelectric voltage response.