Vibration control of a lead zirconate titanate structure considering controller–structure interactions

Abstract

This study focuses on integrating an active vibration controller into the finite element model of a piezoelectric laminated plate with the controller–structure interactions considered. A finite element model of a piezoelectric laminated plate is formulated using the third-order shear deformation theory. A state-space model is set up by performing a system identification technique. The state-space model is then used to design an optimal vibration controller. Considering that the finite element model is more appropriate than state-space model for dynamic simulation, the state-space model-based controller is integrated into the finite element model to capture the controller–structure interactions. The results obtained by applying vibration controller in state-space model are also presented to make a comparison. It is numerically demonstrated that the controller–structure interactions occur and cause performance degradation in case that the state-space model-based controller works with the finite element model. There is no prior guarantee that a state-space model-based controller satisfying the control requirements still works well in closed loop with the finite element model. The results of this study can be used to evaluate the controller performance for the piezoelectric smart structures during the preliminary design stage.

Introduction

The attempts at solving vibration control of lightweight, flexible, and lightly damped structures, such as satellites and other large-scale space structures, have stimulated extensive researches into smart structures over the past years. A smart structure is capable of sensing and reacting to its environment in a desired manner through the integration of sensors, actuators, and associated control systems. There are a number of smart materials available that can be integrated into the smart structures, such as electrorheological fluids, shape memory alloys, piezoelectric materials, and magnetostrictive materials. Among all these materials, piezoelectric materials, such as lead zirconate titanate (PZT) and polyvinylidene fluoride (PVDF), are widely used as sensors and actuators owing to their intrinsic piezoelectric-mechanical coupling effects, rapid response, and large force transmission. Consequently, the study of laminated composites with surface bonded or embedded piezoelectric materials has attracted many researchers recently because of their potential benefits in shape control, vibration suppression, and damage detection. The design of piezoelectric laminated smart structures relies upon multidiscipline knowledge, such as mechanics, mechanical engineering, and control theory.

For the purpose of vibration control of a piezoelectric laminated smart structure, it is first required to investigate its dynamic characteristics. The finite element (FE) method is capable of dealing with the piezoelectric smart structures.^1–3 It has been widely accepted that one of the pioneering works in FE modeling of piezoelectric smart structures is due to Allik and Hughes ⁴ who investigated the electricity–elasticity interactions using a solid element. Following this research, numerous FE models have been developed to analyze the dynamic behaviors of piezoelectric smart structures. Due to the geometric characteristics of piezoelectric laminated structures, one-dimensional models and two-dimensional models are often used, which include the classical laminated plate theory,^5–7 first-order shear deformation theory,^8,9 and third-order shear deformation theory (TSDT).^10,11 The classical laminated plate theory is based on the Kirchhoff–Love assumption in which the transverse shear deformation effects are neglected. The first-order shear deformation theory implies a linear displacement variation through the thickness and results in a constant transverse shear strain through the thickness, so it needs a shear correction factor which is very difficult to be determined especially for arbitrarily laminated composite structures with piezoelectric layers. The TSDT accommodates quadratic variation of transverse shear strains and eliminates the transverse shear stresses on the top and bottom of a laminated composite plate. Accordingly, no shear correction factor is required in the TSDT. Another model for thick laminated composite structures is based on a layerwise theory^12,13 which assumes a biquadratic model for the electric potential distribution along the thickness. These models have their own advantages and disadvantages in terms of accuracy, speed of convergence, and computational cost. Mackerle ¹⁴ and Benjeddou ¹⁵ have presented a detailed bibliographical review of the FE methods applied to the analysis and simulation of smart piezoelectric structures.

In general, a FE model must include many degrees of freedom (DOFs) in order to obtain an accurate prediction of the dynamical behavior. However, a model with a large number of DOFs is not really suitable for the design of a controller, as its evaluation requires considerable computational efforts. For this reason a model reduction technique, such as mode superposition technique, is applied to develop a reduced-order model, such as state-space model (SSM). Control algorithms are then designed based on SSM of the structure dynamics. A survey on various control algorithms employed for active vibration control of smart structures is presented by Alkhatib and Golnaraghi. ¹⁶ Classical constant gain velocity feedback (CGVF) control algorithm is widely employed for active vibration control of smart beams and plates.^17,18 Paper by Tjahyadi et al. ¹⁹ deals with the vibration control of a flexible cantilever beam utilizing adaptive resonant control method. Kircali et al. ²⁰ applied the modal truncation method for dynamic modeling of a beam and then designed a H∞ controller to suppress the first two flexural vibrations of the beam. Manning et al. ²¹ presented a smart structure vibration control scheme using system identification and pole placement technique. Yang et al. ²² presented a feedforward adaptive controller based on an adaptive filter for system dynamics identification of composite laminated smart structures with experimental verifications. Vasques and Dias Rodrigues ²³ presented a comparison between the CGVF control algorithm and optimal control strategies (linear quadratic regulator, i.e. LQR, and linear quadratic Gaussian controller, i.e., LQG) to investigate their effectiveness to suppress vibrations in piezoelectric smart structures. Unlike CGVF controller, the LQR controller with the assumption that the full state variables are available^24,25 does not require collocated sensor/actuator (S/A) pairs. The LQG controller which does not require the measurement of all the state variables is formed by employing an optimal state observer along with the LQR.

It can be seen that, in most of the current studies, if not all, the traditional method for design of the piezoelectric smart structures is divided into two stages: creating an open-loop FE model and performing a model reduction procedure, then setting up a closed-loop control system based on a simplified mathematical model or a SSM. Although the smart structure which includes sensors, actuators, and control algorithm is a highly integrated physical system, dynamic modeling of the smart structures and design of the controller are usually carried out separately without taking the controller–structure interactions into account. In practice, the controller operates in a closed loop with the actual structure and unwanted controller–structure interactions may occur due to unmodeled dynamics. The accuracy and efficiency of active vibration control models are highly dependent on the perfection of understanding the interactions between the structure and the controller.

The basic idea of this paper is that the FE model is more accurate than SSM for modeling of physical systems. Consequently, integration of the controller into FE model not only removes the need to use a SSM in the closed-loop simulation but also leads to better solutions than ever before. Moreover, it simplifies and speeds up the design and simulation processes for vibration control of smart structures. Malgaca and Karagulle ²⁶ have shown the feasibility and efficiency of integration of controller into ANSYS to perform closed-loop simulations.

This paper presents a study on integrating an active vibration controller into FE code. The objective of the present study is to develop a general scheme for piezoelectric smart structure analysis and controller performance evaluation with the controller–structure interactions being taken into account. Based on TSDT, a FE model is formulated for active vibration control of laminated composite plates with integrated piezoelectric sensors and actuators. Nonconforming four-noded rectangular elements with potential DOFs are used to model the laminated composite plate. Considering the FE simulations are analogous to performing experimental investigations, a system identification technique known as observer/Kalman filter identification (OKID) approach together with eigensystem realization algorithm (ERA) method is employed to obtain a multi-input–multi-output (MIMO) SSM for the piezoelectric laminated smart structure from the inputs and outputs of the FE simulations. An optimal controller is then designed to control unwanted vibration response of the piezoelectric structure based on the SSM. At every time instant in the transient analysis, the actuator voltages are updated in accordance with the controller outputs. It is worth pointing out that by this way the controller outputs are eventually fed back to the FE model but not the SSM. Therefore, the controller performance can be evaluated in the context of a closed-loop FE environment which takes the controller–structure interactions into consideration. Numerical examples are presented to demonstrate the efficiency of the proposed scheme in simulating an actively controlled piezoelectric structure.

Mathematical modeling

The lamina constitutive equations

Figure 1 shows a laminated composite plate of total thickness $h$ with integrated piezoelectric sensors and actuators. The $x-y$ plane of the coordinate system $x-y-z$ is defined at the midplane. The $k\text{th}$ layer from the bottom has bottom surface at $z=z_{k-1}$ . The fiber reinforced lamina is assumed to be orthotropic and the fiber direction is indicated by ${\alpha }_k$ , which is the positive rotation angle of the principal axes from the arbitrary $x-y$ axes. Some of the layers can be piezoelectric materials with poling along $z$ .

Figure 1.

A piezoelectric laminated composite plate.

Considering the piezoelectric sensors and actuators, while retaining the orthotropic behavior of the fiber reinforced layer, the linear coupling between mechanical and electrical effects can be described by constitutive equations as follows ⁵

$$\begin{array}{l}\sigma =c\epsilon -e^{\text{T}}E\\ D=e\epsilon +{\mathrm{ϵ}}^{\epsilon }E\end{array}$$

(1)

where $\sigma$ and $\epsilon$ denote the components of stress tensor and strain tensor, respectively and $E$ and $D$ denote electric field and electric displacement, respectively. $c$ is the elasticity matrix, ${\mathrm{ϵ}}^{\epsilon }$ is the dielectric constant matrix, and $e$ is the piezoelectric constant matrix. Equation (1) is valid for both piezoelectric layer and fiber reinforced layer. In case of fiber reinforced layer, the material parameters $e$ and ${\mathrm{ϵ}}^{\epsilon }$ should be zero.

Approximations of displacements and electric potential

Using TSDT, ²⁷ the displacement components of the laminated composite plate take the following form

$$\begin{array}{l}u\left(x,y,z,t\right)=u_0\left(x,y,t\right)+z{\theta }_x\left(x,y,t\right)-\frac{4z^3}{3h^2}\left({\theta }_x+\frac{\partial w_0}{\partial x}\right)\\ v\left(x,y,z,t\right)=v_0\left(x,y,t\right)+z{\theta }_y\left(x,y,t\right)-\frac{4z^3}{3h^2}\left({\theta }_y+\frac{\partial w_0}{\partial y}\right)\\ w\left(x,y,z,t\right)=w_0\left(x,y,t\right)\end{array}$$

(2)

where $t$ denotes time; $u$ , $v$ , and $w$ are the displacement components in the piezoelectric composite plate space in the direction of $x$ , $y$ , and $z$ , respectively; $u_0$ , $v_0$ , and $w_0$ are the displacements of a generic point on the middle plane ( $z=0$ ). ${\theta }_x$ and ${\theta }_y$ are the rotations of the normal to the middle plane, about $y-$ and $x-$ axes, respectively.

We define

$$U={\left[\begin{array}{lll}u& v& w\end{array}\right]}^{\text{T}}$$

(3)

$$\overline{\mathbf{U}}={\left[\begin{array}{lllllll}{u}_0& v_0& w_0& \frac{\partial w_0}{\partial x}& \frac{\partial w_0}{\partial y}& {\theta }_x& {\theta }_y\end{array}\right]}^{\text{T}}$$

(4)

where $\overline{\mathbf{U}}$ is a generalized displacement vector. Then equation (2) can be expressed in matrix form as

$$U=Z\overline{\mathbf{U}}$$

(5)

where

$$Z=\left[\begin{array}{lllllll}1& 0& 0& -\frac{4z^3}{3h^2}& 0& z-\frac{4z^3}{3h^2}& 0\\ 0& 1& 0& 0& -\frac{4z^3}{3h^2}& 0& z-\frac{4z^3}{3h^2}\\ 0& 0& 1& 0& 0& 0& 0\end{array}\right]$$

(6)

The infinitesimal strain components associated with the displacements are as follows

$$\epsilon =TU$$

(7)

where $T$ is the derivation operator defined as

$$T=\left[\begin{array}{lll}\frac{\partial }{\partial x}& 0& 0\\ 0& \frac{\partial }{\partial y}& 0\\ 0& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial z}& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}& 0\end{array}\right]$$

(8)

Supposing the $k\text{th}$ layer is a piezoelectric layer. Grounding the surface at $z=z_{k-1}$ , the transverse variation of the electric potential can be expressed as

$${\varphi }_3(x,y,z,t)=\varphi (x,y,t)\frac{z-z_{k-1}}{z_k-z_{k-1}}$$

(9)

where ${\varphi }_3$ is the electric potential which varies across the $z-$ direction and $\varphi$ is the electric potential of the top surface of the $k\text{th}$ layer.

FE formulation

The nonconforming four-noded rectangular element is used to discretize the piezoelectric laminated plate. Each node has seven mechanical DOFs and one electrical potential DOF for each piezoelectric layer. The efficiency of this type of element for modeling of piezoelectric laminated structures has been demonstrated by Reddy. ²⁷ The continuous mechanical displacements $u_0$ , $v_0$ , ${\theta }_x$ , ${\theta }_y$ and the electrical potential $\varphi$ can be approximated by using the FE shape functions as

$$u_0\left(x,y,t\right)={\displaystyle \sum _{i=1}^4{\psi }_i\left(x,y\right)u_i^e\left(t\right)}$$

(10)

$$v_0\left(x,y,t\right)={\displaystyle \sum _{i=1}^4{\psi }_i\left(x,y\right)v_i^e\left(t\right)}$$

(11)

$$\begin{array}{l}{w}_0(x,y,t)={\displaystyle \sum _{i=1}^4[g_{i1}(x,y)w_i^e(t)+g_{i2}(x,y)\frac{\partial w_i^e(t)}{\partial x}}\\ +g_{i3}(x,y)\frac{\partial w_i^e(t)}{\partial y}]\end{array}$$

(12)

$${\theta }_x(x,y,t)={\displaystyle \sum _{i=1}^4{\psi }_i(x,y){\theta }_{xi}^e(t)}$$

(13)

$${\theta }_y(x,y,t)={\displaystyle \sum _{i=1}^4{\psi }_i(x,y){\theta }_{yi}^e(t)}$$

(14)

$$\varphi (x,y,t)={\displaystyle \sum _{i=1}^4{\psi }_i(x,y){\varphi }_i^e}$$

(15)

where ${\psi }_i$ and $g_{ij}(j=1,2,3)$ are the shape functions which are given by Reddy ²⁷

For each element, the nodal electric potential vector is

$${\Phi }^e={\left[\begin{array}{llll}{\varphi }_1^e& {\varphi }_2^e& {\varphi }_3^e& {\varphi }_4^e\end{array}\right]}^{\text{T}}$$

(16)

and the nodal displacement vector is

$${\text{q}}^e={\left[\begin{array}{llll}{q}_1^{\text{T}}& q_2^{\text{T}}& q_3^{\text{T}}& q_4^{\text{T}}\end{array}\right]}^{\text{T}}$$

(17)

where ${\varphi }_i^e$ denotes the electric potential of the $i\text{th}$ node and $q_i$ is the displacement vector of the $i\text{th}$ node which can be written as

$$q_i={\left[\begin{array}{lllllll}{u}_i^e& v_i^e& w_i^e& \frac{\partial w_i^e}{\partial x}& \frac{\partial w_i^e}{\partial y}& {\theta }_{xi}^e& {\theta }_{yi}^e\end{array}\right]}^{\text{T}}$$

(18)

Substituting equations (10) to (14) into equation (5) yields

$$U=N_u{q}^e$$

(19)

where $N_u$ is the matrix of displacement shape functions.

Substituting equation (19) into equation (7) yields

$$\epsilon =T{N}_u{q}^e=B_u{q}^e$$

(20)

where $B_u$ is strain–displacement matrix. Based on equations (15) and (9), it is deduced that

$${\varphi }_3(x,y,z,t)=N_{\varphi }{\Phi }^e$$

(21)

where $N_{\varphi }$ is the shape function of the electric potential which can be expressed as

$$N_{\varphi }=\frac{z-z_{k-1}}{h_k}\left[\begin{array}{llll}{\psi }_1& {\psi }_2& {\psi }_3& {\psi }_4\end{array}\right]$$

(22)

The electric field vector is the negative gradient of the electric potential ${\varphi }_3$

$$E=-\nabla {\varphi }_3=-\nabla N_{\varphi }{\Phi }^e=-B_{\varphi }{\Phi }^e$$

(23)

where $B_{\varphi }$ is the electric field-potential matrix.

The Hamilton’s variational principle is adopted in the derivation of the elementary dynamic equation, which is as follows

$$\begin{array}{l}{M}_{uu}^e{\ddot{\mathbf{q}}}^e+K_{uu}^e{q}^e+K_{u\varphi }^e{\Phi }^e=F_u^e\\ K_{\varphi u}^e{q}^e+K_{\varphi \varphi }^e{\Phi }^e=F_{\varphi }^e\end{array}$$

(24)

where $M_{uu}^e$ is the elementary mass matrix, $K_{uu}^e$ is the mechanical stiffness matrix, $K_{u\varphi }^e$ and $K_{\varphi u}^e$ are the piezoelectric coupling matrices, $K_{\varphi \varphi }^e$ is the dielectric permittivity matrix, $F_u^e$ is the nodal external force vector, and $F_{\varphi }^e$ is the nodal externally applied charge vector

$$M_{uu}^e={\displaystyle {\int }_{v_e}\rho N_u^{\text{T}}{N}_u\text{d}v}$$

(25)

$$K_{uu}^e={\displaystyle {\int }_{v_e}{B}_u^{\text{T}}c}{B}_u\text{d}v$$

(26)

$$K_{u\varphi }^e={\displaystyle {\int }_{v_e}{B}_u^{\text{T}}{e}^{\text{T}}}{B}_{\varphi }\text{d}v$$

(27)

$$K_{\varphi \varphi }^e=-{\displaystyle {\int }_{v_e}{B}_{\varphi }^{\text{T}}{\mathrm{ϵ}}^{\epsilon }}{B}_{\varphi }\text{d}v$$

(28)

$$K_{\varphi u}^e={\left[K_{u\varphi }^e\right]}^{\text{T}}$$

(29)

$$F_u^e={\displaystyle {\int }_{v_e}{N}_u^{\text{T}}{P}_b\text{d}v}+{\displaystyle {\int }_{s_1}{N}_u^{\text{T}}{P}_s\text{d}s}$$

(30)

$$F_{\varphi }^e={\displaystyle {\int }_{s_2}-N_{\varphi }^{\text{T}}{q}_s\text{d}s}$$

(31)

in which $P_b$ is the body force and $P_s$ is the surface force, respectively. For the element without piezoelectric layers, $K_{u\varphi }^e$ , $K_{\varphi \varphi }^e$ , and $F_{\varphi }^e$ are to be omitted so that equation (24) is valid for a conventional laminated composite element. Supposing that the electric potential vector is composed of a sensor component ${\Phi }^s$ and an actuator component ${\Phi }^a$ , the global system governing equations of motion derived from equation (24) can be written as

$$M_{uu}\ddot{\mathbf{q}}+C_{uu}\dot{\mathbf{q}}+K_{uu}q+K_{u\varphi }^a{\Phi }^a+K_{u\varphi }^s{\Phi }^s=F_u$$

(32)

$$K_{\varphi u}^{a}q+K_{\varphi \varphi }^a{\Phi }^a=F_{\varphi }$$

(33)

$$K_{\varphi u}^{s}q+K_{\varphi \varphi }^s{\Phi }^s=\text{0}$$

(34)

in which $C_{uu}$ is defined as a proportional damping

$$C_{uu}=\alpha M_{uu}+\beta K_{uu}$$

(35)

where $\alpha$ is the mass proportional Rayleigh damping coefficient and $\beta$ is the stiffness proportional Rayleigh, respectively. Using equation (34), the sensory electric potentials can be written as

$${\Phi }^s=-{\left(K_{\varphi \varphi }^s\right)}^{-1}{K}_{\varphi u}^{s}q$$

(36)

Considering equation (36), equation (32) can be written as

$$M_{uu}\ddot{\mathbf{q}}+C_{uu}\dot{\mathbf{q}}+\left(K_{uu}-K_{u\varphi }^s{\left(K_{\varphi \varphi }^s\right)}^{-1}{K}_{\varphi u}^s\right)q=F_u-K_{u\varphi }^a{\Phi }^a$$

(37)

where $F_u$ is the mechanical force and $-K_{u\varphi }^a{\Phi }^a$ is the electric force exerted by the actuators, respectively. Equation (37) can be used to calculate the response under external excitations, whereas equation (36) can be used to calculate the sensor voltages which provide feedback to actuators for vibration control application.

System realization and model reduction

Identification of system Markov parameters

System identification techniques are widely employed for developing a mathematical model by given experimental data. The electromechanical coupling FE model of piezoelectric smart structures provides time histories of the output signals of piezoelectric sensors; therefore, FE simulations using equation (37) are analogous to performing experimental investigations where the only direct outputs are time histories. Consequently, a system identification technique is a logical choice to set up a SSM for control law design. In this study, a system identification technique known as OKID^28,29 technique is employed. The OKID method is formulated entirely in the time domain and is capable of handling general response data. In the following we will summarize it in brief.

Supposing the FE model of a piezoelectric smart structure described by equation (37) can also be represented by the set of first-order difference equations of the form

$$\begin{array}{l}x(i+1)=Ax(i)+Bu(i)\\ y(i)=Cx(i)+Du(i)\end{array}$$

(38)

where $x$ is an $n\times 1$ state vector, $u$ is an $m\times 1$ input or control vector, and $y$ is a $q\times 1$ output vector. As to the piezoelectric smart structure mentioned above, the inputs are the externally exerted mechanical force and the applied voltages over the actuators and the outputs are the sensor voltages. Matrices $A$ , $B$ , $C$ , and $D$ are the state matrix, input matrix, output matrix, and direct influence matrix, respectively. The integer $i$ is the sample indicator. The input–output description of the above system with zero initial conditions can be obtained from equation (38) recursively as

$$y\left(i\right)={\displaystyle \sum _{\tau =0}^{i-1}{Y}_{\tau }u\left(i-\tau -1\right)}+Du\left(i\right)$$

(39)

where the parameters

$$Y_{\tau }=C{A}^{\tau }B$$

(40)

together with $D$ are known as the Markov parameters of the system, which are also the system pulse response samples. To reduce the number of Markov parameters needed to adequately model a system, an observer is introduced. If $(A,C)$ is an observable pair, then there exists an observer of the form

$$\begin{array}{l}\stackrel{⌢}{x}(i+1)=A\stackrel{⌢}{x}(i)+Bu(i)-G[y(i)-\stackrel{⌢}{y}(i)]\\ =(A+GC)\stackrel{⌢}{x}(i)+(B+GD)u(i)-Gy(i)\\ \widehat{\mathbf{y}}(i)=C\stackrel{⌢}{x}(i)+Du(i)\end{array}$$

(41)

The matrix $G$ can be interpreted as an observer gain. Consider the special case where $G$ is a deadbeat observer gain such that all eigenvalues of $A+GC$ are zero, the estimated state $\widehat{\mathbf{x}}$ converges to the true state $x(i)$ after at most $n$ steps where $n$ is the order of the system. Then equation (41) becomes

$$\begin{array}{l}x(i+1)=(A+GC)x(i)+(B+GD)u(i)-Gy(i)\\ y(i)=Cx(i)+Du(i)\end{array}\text{ for }i\ge n$$

(42)

The input–output description of the system in equation (42) is

$$y\left(i\right)={\displaystyle \sum _{\tau =0}^{n-1}{\overline{\mathbf{Y}}}_{\tau }{\left[\begin{array}{ll}u\left(i-\tau -1\right)& y\left(i-\tau -1\right)\end{array}\right]}^{\text{T}}}+Du(i)\text{ for }i\ge n$$

(43)

where

$$\begin{array}{l}{\overline{\mathbf{Y}}}_{\tau }=\left[C{\left(A+GC\right)}^{\tau }\left(B+GD\right)-C{\left(A+GC\right)}^{\tau }G\right]\\ =\left[\begin{array}{ll}{\overline{\mathbf{Y}}}_{\tau }^{(1)}& {\overline{\mathbf{Y}}}_{\tau }^{(2)}\end{array}\right]\end{array}$$

${\overline{\mathbf{Y}}}_{\tau }$ and $D$ are the Markov parameters of an observer system. A particular feature of the deadbeat observer is that the Markov parameters ${\overline{\mathbf{Y}}}_{\tau }$ will become identically zero after a finite number of time steps. For simplicity we define

$$\begin{array}{l}\alpha =\left[\begin{array}{llll}{\overline{\mathbf{Y}}}_0^{(1)}& {\overline{\mathbf{Y}}}_1^{(1)}& \cdots & {\overline{\mathbf{Y}}}_{n-1}^{(1)}\end{array}\right]\\ \beta =\left[\begin{array}{llll}{\overline{\mathbf{Y}}}_0^{(2)}& {\overline{\mathbf{Y}}}_2^{(2)}& \cdots & {\overline{\mathbf{Y}}}_{n-1}^{(2)}\end{array}\right]\end{array}$$

$$\begin{array}{c}\underset{¯}{u}(i-n)=\left[\begin{array}{cccc}{\mathbf{u}}^{\text{T}}(i-1)& {\mathbf{u}}^{\text{T}}(i-2)& \cdots & {\mathbf{u}}^{\text{T}}(i-n)\end{array}\right]\\ \underset{¯}{y}(i-n)=\left[\begin{array}{cccc}{\mathbf{y}}^{\text{T}}(i-1)& {\mathbf{y}}^{\text{T}}(i-2)& \cdots & {\mathbf{y}}^{\text{T}}(i-n)\end{array}\right]\end{array}$$

$$\gamma =[\begin{array}{lll}\alpha & \beta & D\end{array}]$$

$$\Gamma (i-1)={\left[\begin{array}{lll}{\underset{¯}{u}}^{\text{T}}(i-n)& {\underset{¯}{y}}^{\text{T}}(i-n)& {\underset{¯}{u}}^{\text{T}}(i)\end{array}\right]}^{\text{T}}$$

Equation (43) becomes

$$y(i)=\gamma \Gamma (i-1)$$

(44)

Note that equation (44) is in linear form thus the unknown observer parameter vector $\gamma$ can be solved recursively. Let $\stackrel{⌢}{\gamma }(i)$ denotes the estimated parameter vector at time step $i$ , then the standard recursive least-squares solution to equation (44) is simply

$$\begin{array}{l}\stackrel{⌢}{\gamma }\left(i\right)=\stackrel{⌢}{\gamma }\left(i-1\right)+\Delta \stackrel{⌢}{\gamma }\left(i-1\right)\\ \Delta \stackrel{⌢}{\gamma }\left(i-1\right)={\left[\frac{\Theta \left(i-2\right)\Gamma \left(i-1\right){\left[y\left(i\right)-\stackrel{⌢}{\gamma }\left(i-1\right)\Gamma \left(i-1\right)\right]}^{\text{T}}}{1+{\Gamma }^{\text{T}}\left(i-1\right)\Theta \left(i-2\right)\Gamma \left(i-1\right)}\right]}^{\text{T}}\\ \Theta \left(i-1\right)=\Theta \left(i-2\right)-\frac{\Theta \left(i-2\right)\Gamma \left(i-1\right){\Gamma }^{\text{T}}\left(i-1\right)\Theta \left(i-2\right)}{1+{\Gamma }^{\text{T}}\left(i-1\right)\Theta \left(i-2\right)\Gamma \left(i-1\right)}\end{array}$$

(45)

with an arbitrary initial guess $\stackrel{⌢}{\gamma }(0)$ , and $\Theta (-1)$ is any positive definite matrix. The observer parameter vector $\gamma$ can then be used to recover the actual system Markov parameters by

$$Y_{\tau }=C{A}^{\tau }B={\overline{\mathbf{Y}}}_{\tau }^{(1)}+{\displaystyle \sum _{i=0}^{\tau -1}{\overline{\mathbf{Y}}}_i^{(2)}{Y}_{\tau -i-1}}+{\overline{\mathbf{Y}}}_{\tau }^{(2)}D$$

(46)

Juang et al.^28,29 presented a complete description for the computation of the Markov parameters.

Minimal state-space realization

The ERA method ²⁸ uses the system Markov parameters to construct Hankel matrices and then obtain the minimal SSM by singular value decomposition. The algorithm begins by forming the $l\times l$ block Hankel matrix denoted by $H(l,\tau )$

$$H(l,\tau )=\left[\begin{array}{llll}{Y}_{\tau }& Y_{\tau +1}& \cdots & Y_{\tau +l-1}\\ Y_{\tau +1}& Y_{\tau +2}& \cdots & Y_{\tau +l}\\ ⋮& ⋮& & ⋮\\ Y_{\tau +l-1}& Y_{\tau +l}& \cdots & Y_{\tau +2l-2}\end{array}\right]$$

(47)

The order of the system is determined from the singular value decomposition of $H(l,0)$

$$H(l,0)=V_1\Sigma V_2^{\text{T}}$$

(48)

where the matrix $V_1$ and $V_2$ are unitary matrices, $\Sigma$ is an $n\times n$ diagonal matrix of positive singular values, and $n$ is the order of the system. Defining a $q\times lq$ matrix $E_q^{\text{T}}$ and a $m\times lm$ matrix $E_m^{\text{T}}$ made up of identity and null matrices of the form

$$E_q^{\text{T}}=\left[\begin{array}{ll}{I}_q& {\text{0}}_{q\times (l-1)q}\end{array}\right],E_m^{\text{T}}=\left[\begin{array}{ll}{I}_m& {\text{0}}_{m\times (l-1)m}\end{array}\right]$$

(49)

a discrete-time minimal order realization of the system can be written as

$$A_r={\Sigma }^{-1/2}{V}_1^{\text{T}}H(l,1)V_2{\Sigma }^{-1/2}$$

(50)

$$B_r={\Sigma }^{1/2}{V}_2^{\text{T}}{E}_m$$

(51)

$$C_r=E_q^{\text{T}}{V}_1{\Sigma }^{1/2}$$

(52)

The direct influence matrix $D_r$ can be identified by solving equation (43). The above SSM of the system dynamics is used to design a LQG-based robust controller.

It is well known that the rank of the Hankel matrix is the order of the system. For noise-free data, the minimum order realization can be easily obtained by keeping only the nonzero Hankel singular values. However, for real data or data contaminated by noise, the Hankel matrix tends to be full rank, thus making the problem of determining a minimal SSM nontrivial. A criterion deciding which singular values are sufficiently small is very problem dependent. But in all investigated cases there always exists a significant drop in singular values for an increasing position index that signals the order of the system.

Vibration control in FE environment

An optimal controller

In optimal control the feedback control system is designed to minimize the cost function or performance index. A discrete-time minimal realization of the piezoelectric smart structure is represented by $\left(A_r,B_r,C_r,D_r\right)$ . LQR controller^24,25 is employed to determine the optimal state feedback gain matrix by minimizing a cost function of the form

$$J=E\left\{{\displaystyle \sum _{i=1}^{\infty }\left[x^{\text{T}}\left(i\right)Qx\left(i\right)+u^{\text{T}}\left(i\right)Ru\left(i\right)\right]}\right\}$$

(53)

where $Q$ matrix defines the weights on the states while $R$ matrix defines the weights on the control input. The input voltages applied on the actuators are written as

$$u\left(i\right)=-Lx\left(i\right)$$

(54)

in which $L$ is the state feedback gain matrix and is given by

$$L={\left(R+B_r^{\text{T}}P{B}_r\right)}^{-1}{B}_r^{\text{T}}P{A}_r$$

(55)

where the positive definite matrix $P$ is the solution to the following algebraic Riccati equation

$$A_r^{\text{T}}P{A}_r-P-A_r^{\text{T}}P{B}_r{\left(R+B_r^{\text{T}}P{B}_r\right)}^{-1}{B}_r^{\text{T}}P{A}_r+Q=\text{0}$$

(56)

In the state feedback design, it is assumed that all of the state variables are available for feedback. However, this is not always the case and only the system outputs are available for feedback. Thus, one should design a suitable state estimator for the purpose of estimating the states of the system. Among several proposed observers, a state estimator based on the Kalman filter is selected in this study. The dynamic behavior of the Kalman filter is already given by equation (41) and the observer gain matrix $G$ can be obtained by solving equation (43). Note that in practice, the state vector in equation (54) should be replaced by the estimated state $\widehat{\mathbf{x}}(i)$ .

The selection of $Q$ and $R$ is vital in the control design process. A large $Q$ puts higher demand on control result, and a large $R$ puts more limits on control effort. An efficient way of choosing $Q$ is to compute it from $C_r$ with

$$Q=C_r^{\text{T}}{C}_r$$

(57)

$R$ matrix is assumed to be

$$R=\lambda I$$

(58)

where $\lambda$ is a design parameter which is adjusted to meet the control performance.

Control-oriented FE simulation

The overall methodology for integration of a LQR controller into a FE code will be illustrated briefly. The process begins with performing an open-loop FE simulation with prescribed actuator inputs to obtain a set of corresponding output time histories. Band-limited white noise is preferred to be used as the input signal since it satisfies the condition of persistent excitation, as required by a reliable identification. The next step is to employ OKID method and ERA method, using the time histories of inputs and outputs obtained previously, to determine a discrete linear state-space system for use as a control law design model. Next, the LQR controller gain matrix $L$ is computed along with the Kalman filter gain matrix $G$ . Finally, the FE code is to be modified to accept the control law and perform closed-loop simulations. Therefore, the active vibration controller can be integrated into the FE model.

The major modification to the FE code is the updating of the actuator voltages at each time instant in the transient analysis to account for the feedback control. First, the sensor signals are calculated at a given time instant based on equations (36) and (37) and the system states are then estimated using equation (41). And then, the feedback voltages are determined based on equation (54). This process will continue for the predefined time duration of the computed response. Through the results of transient analysis, the control law can be evaluated in the FE environment. If the controller performance is not satisfactory, the controller can be redesigned and evaluated again until the desired performance is obtained.

Case studies and discussion

Validation of numerical code

A simple benchmark problem is solved to validate the numerical code. The results obtained using the developed code are compared with the published results. In this case, a bimorph cantilever piezoelectric plate, which consists of two uniaxially bonded PVDF layers with opposite polarities (see Figure 2) is considered. The length and width of the laminate are 100 and 5 mm, respectively, and the thickness of each lamina is 0.5 mm. The material properties of the PVDF are shown in Table 1. The cantilever plate is modeled with 10 identical plate elements, and each element has the dimension of 10 mm × 5 mm × 1 mm. The cantilever boundary conditions are simulated by specifying zero displacement DOF for the x, y, and z directions at the fixed end of the plate. With a unit voltage applied across the thickness, the deflection of the bimorph plate is calculated by the present FE model and compared with those proposed by Tzou and Ye. ³⁰ To begin with, static analysis is performed to calculate the static deflection of the bimorph plate. For static analysis, equation (37) takes the form

$$\left(K_{uu}-K_{u\varphi }^s{\left(K_{\varphi \varphi }^s\right)}^{-1}{K}_{\varphi u}^s\right)q=F_u-K_{u\varphi }^a{\Phi }^a$$

(59)

Figure 2.

A cantilevered bimorph PVDF plate.

Table 1.

Material properties.

	PVDF	PZT	T300/976
E11 (GPa)	2.0	63.0	150.0
E22 (GPa)	2.0	63.0	9.0
E33 (GPa)	1.0	54.0	9.0
G12 (GPa)	0.775	24.6	7.1
G23 (GPa)	0.400	24.6	2.5
G13 (GPa)	0.775	24.6	7.1
${\nu }_{12}$	0.29	0.3	0.3
e31 (C/m2)	0.046	−5.4	–
e32 (C/m2)	0.046	−5.4	–
e33 (C/m2)	0.000	15.8	–
e15 (C/m2)	–	12.3	–
${\mathrm{ϵ}}_{11}^{\epsilon }$ (F/m)	1.062 × 10−10	1.53 × 10−8	–
${\mathrm{ϵ}}_{22}^{\epsilon }$ (F/m)	1.062 × 10−10	1.53 × 10−8	–
${\mathrm{ϵ}}_{33}^{\epsilon }$ (F/m)	1.062 × 10−10	1.50 × 10−8	–
ρ (kg/m3)	1800	7600.0	1600.0

PVDF: polyvinylidene fluoride; PZT: lead zirconate titanate.

The results obtained by numerical solution of equation (59) are compared with the theoretical and other FE predictions as shown in Figure 3. It can be seen that the present results are much closer to the theoretical solutions than those from the solid FE results of Tzou and Ye, ³⁰ since the solid FE model is relatively rigid and less efficient for thin plates due to the higher stiffness coefficients in the thickness direction. Furthermore, transient analysis is performed to verify the consistency between the steady-state value of the transient response and the static displacement. A time step $\Delta t=0.002\text{s}$ is used for the Newmark method, and Rayleigh type damping is considered with coefficients $\alpha =0.0\text{rad/s}$ and $\beta =0.003\text{s}$ . The tip deflection of the bimorph plate is shown in Figure 4. It can be observed that the steady-state value of transient response converges to the static displacement as time approaches infinity. Therefore, the effectiveness of the proposed FE model is demonstrated numerically.

Figure 3.

Static deflection of the bimorph plate for 1 V across the thickness. FE: finite element.

Figure 4.

Transient response of the tip displacement with constant voltage applied.

Vibration control in FE environment

Having validated the FE code, we begin our discussion on controller design and active vibration control of a piezoelectric smart structure in FE environment. A simply supported square laminated composite plate with four collocated pairs of PZT S/As bonded on the top and bottom surfaces of the plate, as shown in Figure 5, is considered. The plate with side dimension $400\text{mm}\times 400\text{mm}$ is constructed of four layers of T300/976 unidirectional graphite/epoxy composites with a stacking sequence of [−30/30/30/−30]. The thicknesses of each substrate layer and PZT layer are $0.2$ and $0.1\text{mm}$ , respectively. The material properties of the T300/976 composites and PZT are given in Table 1. In this analysis, the whole structure is evenly discretized into 64 identical rectangular elements, and each S/A pair covers four elements. It is supposed that the PZT layer contains metallic electrode on its top and bottom surfaces and thus they form equipotential surfaces. The electric potential DOFs belonging to the same electrode surface are condensed to one electric potential DOF by applying Guyan’s reduction method. ³¹ This electric potential DOF is the potential which is applied across the actuator or the voltage output of the sensor.

Figure 5.

A simply supported laminated plate with four pairs of S/A.

White noise signals, band limited in a frequency range varying from 2 to 200 Hz, which cover six low-frequency modes of the plate, are applied to the four PZT actuators as input data, and the corresponding voltage output data of the four PZT sensors are obtained by performing open-loop FE analysis. The sampling rate of the input–output time is set to 400 Hz for system identification. The prementioned identification technique, i.e. OKID method, is employed to determine the system Markov parameters from the input–output data. For brevity, only the system Markov parameters observed at sensor 4 are shown in Figure 6. It can be seen from Figure 6 that the dynamic characteristics of the structure can be described by the obtained system Markov parameters. The ERA method is then used to find a state-space realization from the system Markov parameters. The order of the SSM is determined upon inspecting the singular values of the corresponding Hankel matrix. Figure 7 shows the singular values of the Hankel matrix. One can clearly observe that there are 12 first-order modes in the data, as evidenced by the sharp drop after the 12th singular value, and the fact that the rest of the singular values are numerical zeros. Hence, a MIMO SSM with the system order 12 is obtained from singular value decomposition of the Hankel matrix. The identified model can be utilized in the controller design procedure.

Figure 6.

System Markov parameters observed at sensor 4: (a) stimulated by actuator 1, (b) stimulated by actuator 2, (a) stimulated by actuator 3, and (a) stimulated by actuator 4.

Figure 7.

Singular values of the Hankel matrix.

After several attempts we set the parameter $\lambda$ as $\lambda =0.1$ to damp the vibration effectively. However, it must be noted that this parameter should be limited for the sake of the breakdown voltage of the piezoelectric materials. The maximum voltage per thickness of the piezoelectric material is taken as $2000\text{V}/\text{mm}$ . Once the LQR controller design is completed, the FE code is then modified to integrate control law into transient analysis procedure. Newmark direct integration method is used to perform closed-loop simulations with a time step $\Delta t=0.001\text{s}$ . For the closed-loop simulations considered in the following numerical examples, the same controller developed previously is integrated into the SSM and FE model, respectively, to make a comparison.

Two cases are considered in vibration control of the simply supported square laminated plate. For case 1, the plate is initially exposed to a uniform distributed load of $50\text{N}/{\text{m}}^2$ . The load is removed by setting the plate into vibration. Figure 8 illustrates the uncontrolled and controlled time response of the transverse displacement of node 41, which can be localized through Figure 5. Time histories of sensor voltages and actuator voltages of the FE model-based closed-loop simulation are presented in Figures 9 and 10. For brevity, only two of them are shown. The sensor voltages and actuator voltages vary as the plate vibrates, and the vibrational period is the same as the period of the laminated plate. From Figures 8 to ¹⁰ , one can infer that the controller successfully damps the transverse vibration of the laminated plate while the actuator voltages do not exceed the breakdown voltage. However, the performances of the same controller observed in SSM and in FE model are not the same. It is observed that the controller is more efficient in SSM environment than in FE environment.

Figure 8.

Transverse displacement of node 41 under initial disturbance. FE: finite element; SSM: state-space model.

Figure 9.

Actuator voltages for case 1: (a) actuator 1 and (b) actuator 2.

Figure 10.

Sensor voltages for case 1: (a) sensor 1 and (b) sensor 2.

For case 2, a distributed transverse harmonic force $5\sin (2\pi ft)\text{N}/{\text{m}}^2$ with frequency $f=26\text{Hz}$ is applied on the plate. The forcing frequency is close to the first natural frequency of the system. The transverse displacement of node 41 is shown in Figure 11. Beat vibration phenomenon occurs in the open-loop vibration due to the small difference between the natural frequency and the forcing frequency. Figures 12 and 13 illustrate, respectively, the actuator voltages and the sensor voltages. As shown in the figures, the controller is also effective in controlling the harmonic vibration. And we can observe that the vibration frequency of the actuator voltages and the sensor voltages is the same as the forcing frequency. However again, the closed-loop responses show that the performances of the same controller observed in SSM and in FE model are not the same.

Figure 11.

Transverse displacement of node 41 under harmonic load. FE: finite element; SSM: state-space model.

Figure 12.

Actuator voltages for case 2: (a) actuator 1 and (b) actuator 2.

Figure 13.

Sensor voltages for case 2: (a) sensor 1 and (b) sensor 2.

Generally, the SSM contains only the low-frequency modes. The high-frequency modes which are not included in the SSM are termed as unmodeled dynamics. Compared with SSM, it is generally considered that the FE model contains the unmodeled dynamics. Therefore, when the SSM-based controller operates in closed loop with FE model, the unmodeled dynamics are to be excited. In this way, the controller–structure interactions can then be captured numerically. Usually, we feed back the SSM-based controller outputs to the SSM to evaluate the controller performance. However, one can infer from Figures 8 and 11 that the performances of the same controller observed in SSM and in FE model are not the same due to unmodeled dynamics. That is to say, it is not adequate to evaluate the controller performance in SSM environment. Furthermore, it cannot be concluded that a SSM-based controller satisfying the control requirements still works well when it is connected to the real plant.

Conclusions

An active vibration controller is integrated into the FE model of a piezoelectric laminated plate to conduct vibration control simulations, with a focus on taking the controller–structure interactions into consideration. By using TSDT, an accurate open-loop FE model is developed and a system identification technique known as OKID approach together with ERA method is employed to obtain a reduced-order SSM. An optimal feedback controller is designed based on SSM and then integrated into the FE code by updating the actuator voltages at each time instant in the transient analysis. By this way the controller outputs are eventually fed back to the FE model, but not the SSM as in some traditional active vibration control methods. Because the FE model represents the real system much better than a SSM, the proposed FE model-based closed-loop simulations can be used to capture the controller structure interaction. Numerical examples are presented to demonstrate the efficiency of the proposed scheme in the simulation of an actively controlled piezoelectric structure. It is numerically demonstrated that there is no prior guarantee that a SSM-based controller satisfying the control requirements still works well in closed loop with the FE model due to unmodeled dynamics. It is strongly recommended to evaluate the performance of an active vibration controller by using FE model-based closed-loop simulations but not SSM-based closed-loop simulations. Future research calls for investigating unwanted controller–structure interactions occurring due to unmodeled dynamics, such as control spillover, based on the proposed method.