Closed loop finite element modeling of piezoelectric smart structures

The objective of this paper is to develop a general design and analysis scheme for actively controlled piezoelectric smart structures. The scheme involves dynamic modeling of a smart structure, designing control laws and closed-loop simulation in a finite element environment. Based on the structure responses determined by finite element method, a modern system identification technique known as Observer/Kalman filter Identification (OKID) technique is used to determine the system Markov parameters. The Eigensystem Realization Algorithm (ERA) is then employed to develop an explicit state space model of the equivalent linear system for control law design. The Linear Quadratic Gaussian (LQG) control law design technique is employed to design a control law. By using ANSYS parametric design language (APDL), the control law is incorporated into the ANSYS finite element model to perform closed loop simulations. Therefore, the control law performance can be evaluated in the context of a finite element environment. Finally, numerical examples have demonstrated the validity and efficiency of the proposed design scheme. Without any further modifications, the design scheme can be readily applied to other complex smart structures.


Introduction
Advanced applications, such as large-scale space structures, jet fighters and concept automobiles, have structural requirements including high strength, lightweight, high structural damping and low acoustic noise. Attempts at achieving these characteristics have recently stimulated extensive research into smart structures, which have been widely applied in the field of active vibration control of structures. Piezoelectric materials, such as lead zirconate titanate (PZT), have coupled mechanical and electrical properties. Bonding or embedding piezoelectric patches in a smart structure can act as sensors to monitor or as actuators to control the response of the structure [1].
To design piezoelectric smart structures for active vibration control, both structural dynamics and control theory need to be considered. The finite element method, a widely accepted and powerful tool for analyzing complex structures, is capable of dealing with the piezoelectric smart structures. Special finite elements have been developed to account for piezoelectric effect [2][3][4]. ANSYS, one of the most widely-used commercially available finite element analysis (FEA) codes, has the ability to model piezoelectric materials. Unfortunately, finite element models are always of high-order and can not be used directly in controller design whereas most control law design methods require an explicit mathematical model of the system to be controlled. To overcome this difficulty, the mode superposition technique is widely used to transform the coupled finite element equations of motion in the physical coordinates into a set of reduced uncoupled equations in the modal coordinates, and then the state space model of the system is developed [5][6][7]. Commonly, the MatLAB Control System Toolbox is used for the control law design and then closed loop simulation is performed based on the state space model [6,8]. As can be seen, in all these research works [5][6][7][8], dynamic modeling of the smart structures and design of the control systems are usually carried out separately without taking the interaction between the structure and the control system into account. Due to the increasing interest in the design of complex piezoelectric smart structures and the need for fast and simple implementation of piezoelectric control systems, it is necessary to develop a general design scheme of actively controlled piezoelectric smart structures.
In this paper, a modern system identification technique is employed to develop an explicit state space model for control law design from the output of a commercial finite element code ANSYS. System identification technique deals with the realization of a model from the experiment results for the purpose of control law design. This requires the construction of a minimum order model from the test data that characterizes the dynamics of the system at the selected control and measurement positions. A finite element simulation provides time histories of the variables of a piezoelectric smart structure system, but does not generate an explicit mathematical model of the system. Finite element simulations are analogous to performing experimental investigations where the only direct outputs are time responses. Therefore, the structural responses determined by finite element simulations are employed as the input to an identification technique known as Observer/Kalman filter Identification (OKID) [9,10] approach that identifies the Markov parameters of an asymptotically stable observer. The system Markov parameters are then determined recursively from the Markov parameters of the observer system, from which a state space realization is obtained using the Eigensystem Realization Algorithm (ERA) [9,11]. The LQG design method is employed to design a control law based on the state space model. ANSYS parametric design language (APDL) is then used to integrate the control law into the ANSYS finite element model to perform closed loop simulations [12].
In summary, under the interaction of structure and control disciplines, the objective of this paper is to develop a general analysis and design scheme of piezoelectric smart structures. The proposed scheme can make use of ANSYS and implement a user selected control algorithm. By using APDL, the control law is incorporated into the ANSYS finite element model of the smart structure to perform closed loop simulations. The proposed design scheme can further reduce cost, as well as shorten design cycle time for smart structures. Without any further modifications, it can be readily applied to large flexible space structures.

Finite element model
Modeling of piezoelectric smart structures by the finite element method and its implementation for active vibration control problems has been presented [2][3][4]13,14]. The global matrix equations governing a smart structure system can be written as: where u denotes structural displacement, φ denotes electric potential, and a dot above a variable denotes a time derivative; F u denotes the structural load and F φ denotes the electric load; M uu and K uu are the structural mass and stiffness matrix, respectively, K uφ is the piezoelectric coupling matrix and K φφ is the dielectric stiffness matrix, C uu is the damping matrix, which is usually assumed to be a linear combination of the structural mass matrix M uu and the structural matrix K uu The constant α and β are the damping parameters [15]. It is notable to point out that proportional damping is assumed but not required for system identification.
As mentioned before, the commercial finite element code ANSYS can be utilized to model the piezoelectric smart structure with prescribed inputs to obtain a set of corresponding output time histories. The next step is to employ a system identification technique, using the time histories of outputs and inputs of the ANSYS finite element model, to determine an "equivalent linear system" for use as a control law design model.

Identification of the system Markov parameters by OKID approach
The Observer/Kalman filter Identification (OKID) [9,10] method was developed to compute the Markov parameters of a linear system, which are the same as its pulse response samples. The method is formulated entirely in the time domain, and is capable of handling general response data. One of the keys to the OKID algorithm is the introduction of an observer into the identification process. The first step of the process is the calculation of the observer Markov parameters. Then the system Markov parameters are determined recursively from the Markov parameters of the observer system. In the following we summarize it in brief.
Firstly, consider a general discrete multivariable linear system expressed in the state space format where x is an n × 1 state vector, u an m × 1 input or control vector, y a q × 1 output vector. Matrices A, B, C and D are respectively the state matrix, input matrix, output matrix, and direct influence matrix. The integer i is the sample indicator. The input-output description of the above system with zero initial conditions can be obtained from Eq. (3) as where the parameters Y τ = CA τ B 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 into the OKID technique. If (A, C) is an observable pair, then there exists an observer of the form The matrix M can be interpreted as an observer gain. Consider the special case where M is a deadbeat observer gain such that all eigenvalues of A + M C are zero, the estimated statex converges to the true state x(i) after at most n steps where n is the order of the system. Equation (5) then becomes The input-output description of the system in Eq. (6) is τ ] Y τ and D are the Markov parameters of an observer system. A particular feature of the deadbeat observer is that the Markov parametersȲ τ will become identically zero after a finite number of time steps. The standard recursive least-squares technique is used to solve Eq. (7) and then the observer Markov parameters are computed. There is an algebraic relationship between the Markov parameters of the observer system and those of the actual system A complete description for the computation of the Markov parameters is given in reference [9,16]. Once the Markov parameters of the observer system are identified, the actual system Markov parameters can be recovered according to Eq. (8). The actual system Markov parameters are then used to obtain a state space model of the system by a realization procedure called the Eigensystem Realization Algorithm (ERA) [9,11].

Minimum realization by the Eigensystem Realization Algorithm
The Eigensystem Realization Algorithm (ERA) is a system realization method that has been studied extensively [9,11]. It determines a state space model (A r , B r , C r , D r ) of a structure from the system Markov parameters. The algorithm begins by forming the l × l block Hankel matrix denoted by H(l, τ ) The order of the system is determined from the singular value decomposition of H(l, 0) where the matrix U and V are unitary matrices. Assume with Choose a number δ as a zero threshold such that then the matrix H(l, 0) is considered to have rank n. Take Σ n as a n × n diagonal matrix of the first n singular values in Σ. U n and V n correspond to the singular values in Σ n . Equation (10) then becomes Defining a q × lq matrix E T q and an m × lm matrix E T m made up of identity and null matrices of the form a discrete-time minimal order realization of the system can be written as The direct influence matrix D r can be identified by solving Eq. (7). The above state space model of the system dynamics is used to design a LQG based robust controller.

LQG design
In practice, the modeling of the smart structure is seldom perfect and the system parameters like density, elastic stiffness and damping usually vary during the course of operation of the structure due to corrosion, damage, etc. Moreover, the system might be subjected to varying dynamic loads which may excite multiple modes of vibration [7]. A LQG based robust controller which is less sensitive to system parameter variations, and also works efficiently when multiple modes of vibration are present in the response, is an excellent choice for controlling the vibrations of such systems. In this section, the robust controller design using the LQG is discussed. Considering the process and measurement noise w(i) and v(i), the state space equation for the piezoelectric smart structure can be written as In the LQG algorithm, the process noise w(i) and the measurement noise v(i) are both assumed to be stationary, zero-mean, Gaussian white, and to have covariance matrices satisfying where E{·} denotes the expected value, δ i,j denotes the Kronecker delta. In addition, w(i) and v(i) are assumed to be uncorrelated. The controller consists of a state feedback module and a state estimation module. The former corresponds to the following control law in which K denotes the feedback gain matrix. The feedback gain matrix is obtained so as to minimize the quadratic cost function of the form where Q is a symmetric positive semi-definite matrix and R is a symmetric positive-definite matrix. The optimal feedback gain matrix that minimizes the above cost function is given by where the positive definite matrix P 1 is the solution to the following algebraic Riccati equation In Eq. (18) it is assumed that all of the states are available for feedback. However, in practice only the system outputs are available for feedback. To estimate the states of the system from sensor outputs, a state estimator based on the Kalman filter is designed wherex(i) denotes the estimated state and the observer gain matrix is obtained form in which the positive P 2 satisfies the following algebraic Riccati equation Note that the state vector in Eq. (18) should be replaced by the estimated statex(i). One of the merits of the LQG technique lies in the separation principle, which states that the design of the feedback and the estimator can be carried out separately. A block diagram of the LQG controller is shown in Fig. 1. The selection of Q and R is vital in the control design process. Q and R are the free parameters of design and stipulate the relative importance of the control result and the control effort. A large Q puts higher demand on control result, and a large R puts more limit on control effort [17]. The optimal values for Q and R are typically obtained by trial-and-error methods. In this paper, Q is computed from C r with Weighting matrix R can be set as ρI with ρ as a scalar design parameter. Therefore, the challenging task of choosing Q and R reduces to choosing one parameter ρ. Consequently, the control performance such as settling time and the maximum value of the actuator voltage can be tuned through changing the value of ρ. This procedure is used in the illustrative example presented next.

Closed loop simulation in finite element environment
In this section, by using ANSYS parametric design language (APDL), the ANSYS finite element model of a smart structure is modified to accept control laws and perform closed-loop simulations. Therefore, the control law can be evaluated in the context of finite element environment. The block diagram of the analysis is shown in Fig. 2. F e is the vibration generating force. The instantaneous value of the vibration generating force is defined at each time steps. v s is the sensor signal at a sensor location. For instance, for a cantilever plate, vertical velocity can serve as sensor signal. v a is the actuator voltage and it is determined by the LQG controller using F e and v s .
In this study, to construct an ANSYS finite element model for a piezoelectric smart structure, SOLID45 elements are used for the metal part and SOLID5 elements are used for the piezoelectric part of the structure. A macro has been written to perform the transient analysis, in which the control law is accepted. In the transient analysis, the instantaneous value of the sensor signal is calculated at a time step by using the classical Newmark integration method. Then the feedback voltage, which will be used as the input to the piezoelectric actuator during the next time step, is determined by the LQG controller. 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 finite element simulation. If the control law performance is not adequate, the control law can be redesigned and evaluated again until the desired performance is obtained.

Illustrative example
In this section, a cantilevered plate is considered. The structure consists of a host aluminum plate (700 mm × 150 mm × 1.2 mm), APC-850 piezoelectric patches (70 mm × 15 mm × 0.5 mm) and velocity sensors. The adhesive layers are neglected. The piezoelectric patches work as actuators are perfectly bonded on the upper and lower surfaces of the plate structure at the same location. It is assumed that the same magnitude but the opposite direction of electric field is applied to the upper and lower piezoelectric patches so as to create a pure bending moment for vibration suppression. Velocity sensors are bonded on top of the actuators. The configuration of the structure is shown in Fig. 3. As can be seen in Fig. 3, there are three actuator pairs marked with actuator 1, actuator 2 and actuator 3 separately. Also, there are three velocity sensors marked with sensor 1, sensor 2 and sensor 3 separately. It should be pointed out that in ANSYS finite element model the elastic and inertial properties of the sensors are not taken into account. The transverse velocities at sensor locations are calculated at each time step and then serve as sensor signals. The material properties of APC850 and aluminum are listed below. APC850: The finite element model for the cantilevered plate with piezoelectric pitches is shown in Fig. 4. A total of 17248 three-dimensional solid elements are used to mesh the structure. Among those elements, 1008 SOLID5 elements are utilized to model piezoelectric patches, and 16240 SOLID45 elements are employed to model host structure. The time step is chosen as 0.0025.
White noise is selected as the input signal since it satisfies the condition of persistent excitation, as required by a reliable identification [18]. Finite element simulation using ANSYS is performed to obtain a set of corresponding output time histories. The input-output time history is discretized at a sampling rate of 400 Hz for system identification. The aforementioned system identification technique is then employed, using the time histories of inputs and outputs, to set up a state space model for use as a control law design model. Assume an order for the system to be identified, denoted by n. In this study, only the first 6 modes are taken into account, thus n is chosen as 12. The identified model is of order 12 and is given in the Appendix. Frequency response functions can be obtained from the state space model. For brevity, only two of them are shown in Fig. 5. The results from the ANSYS finite element model are also shown in Fig. 5 for the purpose of comparison. Excellent agreement is observed between the two results except for the frequency range above 300 rad/sec. For most structure systems under practical loading conditions, only the first few modes play an important part in the vibration response of the system. Therefore, it can be concluded that the state space model developed by identification technique characterizes the dynamics of the system and can be used for control law design.
For the LQG controller, the free parameters ρ, V 1 and V 2 selected as ρ = 0.001, V 1 = 10 −3 I n and V 2 = 10 −6 I q . As mentioned previously, these parameters are determined by trail and error methods to most effectively control the  structure. At the same time, these parameters are limited for the sake of the breakdown voltage of the piezoelectric materials. The maximum voltage per thickness of the piezoelectric material is taken as 300 V mm −1 . The feedback gain matrix K is given in the Appendix. Once the LQG controller design is completed, APDL is used to integrate control law into the ANSYS finite element model to perform closed loop simulations. Finally, the control law is evaluated in the finite element simulation. If the control law performance is not adequate, the control law can be redesigned and evaluated again until the desired performance is obtained. The external load is assumed to be applied at the free end of the plate and the dynamic responses at the free end, i.e. the tip displacements, are given to evaluate the control law. For case 1, a transverse impulsive loading of 0.1 N is applied at the free end of the plate for 0.1 s. Transient responses of the plate tip displacement are shown in Fig. 6. Results when the control is off are also shown for comparison. Actuator voltages are shown in Fig. 7. As can be seen, the controller successfully damps the vibration of the plate. The actuator voltages are lower than the breakdown voltage of piezoelectric materials. For case 2, banded limited white noise in the frequency range 0-200 rad/sec is applied at the free end of the plate. The plate tip displacement and actuator voltages are shown in Figs 8 and 9 respectively. The controller is also effective in controlling the random vibration. Still, the actuator voltages are lower than the breakdown voltage of piezoelectric materials.

Conclusions
Under the interaction of structure and control disciplines, a general scheme of analyzing and designing piezoelectric smart structures is successfully developed in this paper. The scheme involves using a modern system identification technique to develop an equivalent linear model from the outputs of the commercial finite element code ANSYS. Standard control law design techniques are then used to design control laws. By using ANSYS parametric design language (APDL), the resulting control laws are then incorporated into transient analysis for evaluation. Results of a numerical study applying this methodology to vibration control of a cantilevered plate with piezoelectric actuators are presented.
Numerical results indicate that the equivalent linear models developed by employing a system identification technique can represent the input-output relationship of a finite simulation very well. The proposed scheme can integrate a user selected control algorithm into the ANSYS finite element model to perform closed loop simulation. The control law can be evaluated in the context of a finite element environment. The technique presented in this paper may be also applied to other different vibration control problems, based on passive or hybrid techniques.
The system matrices are listed below. For ease of presentation, the matrices are subdivided. The feedback gain matrix K is given below.