A unified Gram-Schmidt-Ritz formulation for vibration and active flutter control analysis of honeycomb sandwich plate with general elastic support

2.1. The equivalent model of honeycomb sandwich plate

Generally, there are three equivalent models for mechanical analysis of the honeycomb sandwich structures, the sandwich theory, cellular equivalent plate theory and equivalent plate theory. The validity of the equivalent plate theory is verified by Paik tests [1], which gives the elasticity modulus and thickness of plate according to the equal principle of stiffness and density, ultimately making the honeycomb sandwich plate equivalent to an isotropic plate. Supposing the sheet satisfies Kirchhoff hypothesis, the equation can be obtained according to the equivalent bending stiffness:

where, $E_{eq}$, $t_{eq}$ are the equivalent elastic modulus and thickness of the plate, $\nu$ is Poisson ratio. Using the equivalent principle between the equivalent plate and honeycomb sandwich plate, the equation can be written as:

The equivalent plate parameters are finally obtained based on Reissner theory [2]:

where, ${\rho }_c$ is the average mass density of honeycomb core, ${\rho }_c=\rho \left(1.54t/l\right)$.

A FEM model is set up to solve the honeycomb sandwich plate problem adopting the sandwich theory (ANSYS Shell 91), where the core is equivalent to the orthotropic material, and the face sheet is treated as the isotropic material. The material parameters of core layer are determined by means of $Y$ model using the energy method:

where, $E_{hci}$ ($i=x,y,z$) is the elastic modulus of honeycomb core, $G_{hcxy}$, $G_{hcxz}$ and $G_{hcyz}$ are the shear modulus of honeycomb core, ${\rho }_c$ is the density of the honeycomb core. The FEM result using the sandwich theory is regarded as the reference value.

2.2. Constitutive relation of equivalent sandwich plate and piezoelectric material

The equivalent plate structure with the piezoelectric layer in a supersonic airflow, is shown in Fig. 2. The length, width and thickness of the plate are $a$, $b$ and $t_{eq}$ respectively, the airflow is along the x axis direction. The artificial springs are distributed on the edges of plate to model the elastic support conditions. The different boundaries can be obtained by accordingly setting the artificial spring stiffness to be extremely large, small or arbitrary value.

Fig. 2Equivalent plate with the piezoelectric layer under general elastic support

The piezoelectric actuator is pasted on the upper surface of plate, the position coordinates are $\left(x_1,y_1\right)$ and $\left(x_2,y_2\right)$. The plane displacements $u$, $v$ and transverse displacement $w$ of the plate are expressed as:

where, $z$ is the coordinate of any point along the $z$-axis direction. The strain-displacement relation is written as:

where, ${\epsilon }_x$, ${\epsilon }_y$ and ${\gamma }_{xy}$ are the strains in $x$, $y$ direction and shear strain in $x$-$y$ plane, respectively. The constitutive relation of the plate structure is:

where, ${\sigma }_x$, ${\sigma }_y$ and ${\tau }_{xy}$ are the stress in $x$, $y$ direction and shear stress in $x$-$y$ plane, respectively. $Q_{11}=Q_{22}=E_{eq}/(1-{\upsilon }^2)$, $Q_{12}=\upsilon E_{eq}/(1-{\upsilon }^2)$, $Q_{66}=G_{eq}$. The constitutive relation of piezoelectric layer is expressed as:

where, ${\sigma }_x^p$, ${\sigma }_y^p$ and ${\tau }_{xy}^p$ are the stress of piezoelectric material in the $x$, $y$ direction and shear stress in $x$-$y$ plane, respectively. $e_{31}$, ${\epsilon }_{33}$ and $D_3$ are the piezoelectric constants, dielectric constant and electric displacement, respectively. $E_3$ is the electric field strength in the $z$ direction:

where, $h_p$ is the thickness of the piezoelectric layer, $V_0\left(t\right)$ is the applied voltage.

The total kinetic energy of the plate system with the piezoelectric layer is calculated as:

where, $V$ and $V_p$ is the volume of the plate and piezoelectric layer, ${\rho }_p$ is the density of piezoelectric material.

The total potential energy of the plate system with the piezoelectric layer is calculated as:

where, $\mathbf{\epsilon }$ and $\mathbf{\sigma }$ are the strain vector and stress vector of plate structure, ${\mathbf{\epsilon }}^p$ and ${\mathbf{\sigma }}^p$ are the strain vector and stress vector of piezoelectric layer.

2.3. The simulation of general elastic support boundary

The artificial springs are used to describe the connection between the plate and the support members. The springs are evenly distributed on the edges of the plate, and the elastic supports adopt the tension spring and torsion spring, where the corresponding stiffness are $k^w$ and $k^{\vartheta }$, shown in Fig. 3. The key of the remaining problem is to find the admissible functions to expand the displacement components of the plate structure, which are introduced in Part 2.4.

Since the artificial springs have no mass, only the potential energy is calculated as:

where, $k_{0\left(a\right)}^w$, $k_{a\left(a\right)}^w$, $k_{0\left(b\right)}^w$, $k_{b\left(b\right)}^w$ and $k_{0\left(a\right)}^{\vartheta }$, $k_{a\left(a\right)}^{\vartheta }$, $k_{0\left(b\right)}^{\vartheta }$, $k_{b\left(b\right)}^{\vartheta }$ are the stiffness of the tension spring and torsion spring at $x=$ 0 and $x=a$, $y=$ 0 and $y=b$, respectively.

Fig. 3Elastic support plate with a set of artificial springs

2.4. Method of solution

The Rayleigh-Ritz method is used to obtain the frequency and mode functions. The transverse displacement can be expressed as:

where, $\omega$ is the circular frequency, $W(x,y)$ is mode shape which can be determined by Gram-Schmidt orthogonal polynomials, i.e.:

where, $a_{mn}$ is the unknown coefficient, ${\varphi }_m\left(x\right)$ and ${\phi }_n\left(y\right)$ are the orthogonal polynomials in the $x$ and $y$ directions, respectively. $m_t$ and $n_t$ are the truncation number of polynomials.

Substituting Eqs. (5-7) into Eq. (12), only the potential energy of plate structure is considered:

where, $\xi =x/a$, $\eta =y/b$.

Substituting Eqs. (14, 15) into Eqs. (11, 13, 16), the Rayleigh quotient to corresponding factors is written as:

The minimizing condition with energy function yields the generalized eigenvalue problem:

where, $\mathbf{K}$ and $\mathbf{M}$ are $m_t\times n_t$ matrices, $\mathbf{a}={\left[a_{11},a_{12},\cdots ,a_{m_t{n}_t}\right]}^T$. The matrices ${\mathbf{K}}_{\epsilon }$, ${\mathbf{K}}_s$ and $\mathbf{M}$ are determined as:

where, ${\mathbf{K}}_{\epsilon }$ is the stiffness matrix of vibratory stress, ${\mathbf{K}}_s$ is the stiffness matrix of the artificial springs, $m=\mathrm{1,2},\dots ,m_t$, $n=\mathrm{1,2},\dots ,n_t$, $i=\mathrm{1,2},\dots ,m_t$, $j=\mathrm{1,2},\dots ,n_t$.

The construction of Gram-Schmidt orthogonal polynomials is as follows, giving the first term of orthogonal polynomial cluster ${\varphi }_0^{}\left(\xi \right)$, which satisfies the certain boundary condition. According to the steps of orthogonal polynomial cluster configuration, the remaining items can be calculated as:

where, $B_k={\int }_0^1\xi {\varphi }_{{}_{k-1}}^{}(\xi {)}^{2}d\xi /{\int }_0^1{\varphi }_{{}_{k-1}}^{}(\xi {)}^{2}d\xi$, $C_k={\int }_0^1\xi {\varphi }_{k-1}^{}\left(\xi \right){\varphi }_{k-2}^{}\left(\xi \right)d\xi /{\int }_0^1{\varphi }_{{}_{k-2}}^{}(\xi {)}^{2}d\xi .$

Each polynomial is normalized as:

If the first term satisfies the geometry boundary conditions, it is easy to verify that other terms do. The relationship between any two orthogonal polynomials is:

The first item of polynomials with different boundary conditions is shown in Table 1. The letter, F represents for freedom, C for clamped, S for simply supported.

2.5. The governing equations

Using the mode function $W\left(x,y\right)$ in Part 2.4 to derive the governing equation based on the assumed mode method, the displacement of the plate is written as the following form:

where, $q_{ij}\left(t\right)$ is the generalized coordinate, $W_{ij}$ is the corresponding mode function obtained from the Rayleigh-Ritz approach, $m$ and $n$ are the modal order, and:

Substituting Eq. (26) into Eqs. (11-12), the expressions of potential energy and kinetic energy can be obtained:

where:

The potential energy of the plate system with the piezoelectric layer is expressed as:

where:

The aerodynamic pressure caused by the supersonic flow can be described by the quasi-steady-order piston theory, which has a good accuracy in the Mach number $\sqrt{2}<M_{\infty }<5$, the expression can be written [17, 20]:

where, $\zeta ={\rho }_{\infty }{U}_{\infty }^2/\sqrt{M_{\infty }^2-1}$, $\mu ={{\rho }_{\infty }{U}_{\infty }(M_{\infty }^2-2)/(M_{\infty }^2-1)}^{3/2}$ are the aerodynamic stiffness and aerodynamic damping coefficient, respectively. ${\rho }_{\infty }$, $U_{\infty }$ and $M_{\infty }$ are the free flow density, flow velocity and Mach number.

The virtual work can be expressed as:

where, $F$ is the external force, $\mathrm{\Delta }p$ is the aerodynamic load, $A$ is the plate area, and $c_0$ is the damping coefficient, ${\mathbf{F}}_{\mathrm{\Delta }p2}$, ${\mathbf{F}}_{\mathrm{\Delta }p1}$ are the damping and aerodynamic stiffness matrix, respectively:

The governing equations of the plate system is established based on Hamilton principle:

where, $\delta T$, $\delta U$ and $\delta W$ are the virtual work of the kinetic energy, strain energy and external loads. Substituting Eqs. (27-28) and Eq. (30) into Eq. (31), the final motion equation is:

where, $\mathbf{M}$, ${\mathbf{K}}_1$ and ${\mathbf{K}}_p$ are modal mass matrix, stiffness matrix and electro-mechanical coupling matrix, respectively.

3.1. Flutter control based on the displacement feedback method

The piezoelectric material can be used to design sensors and actuators [17]. In this paper, a voltage is applied on the piezoelectric material attached to the surface of the structures, which is designed to the piezoelectric actuator. The centre position displacement of plate is designed as an input, and the voltage value is calculated by the control strategy. The control voltage of the actuator based on the center position displacement can be written as:

where, $K_d$ is a piezoelectric actuator feedback control gain, $w(x_m,y_m,t)$ is the center position displacement of the plate structure.

Substituting Eq. (26) into Eq. (33):

Substituting Eq. (34) into Eq. (32), the control equation can be obtained:

where, $\mathbf{K}=K_d{\mathbf{K}}_p^T{\mathbf{W}}^T(x_m,y_m)$.

3.2. Flutter control based on the LQR method

The LQR method is adopted to control the flutter problem of plate structure, the governing equation is converted into the form of state space:

where $\mathbf{x}$ is the state vector. The LQR algorithm is regarded as a feedback control, making the objective function to obtain the minimum value, and the objective function is usually proportional to the system response. The objective function is defined as:

where, $\mathbf{Q}$ and $\mathbf{R}$ are positive definite and semi-definite matrix, which represent a measure of the effect of input and output on the control result. The feedback voltage can be determined:

where, ${\mathbf{G}}_{Lc}={\mathbf{R}}^{-1}{\mathbf{B}}^T\mathbf{P}$ is a feedback gain, $\mathbf{P}$ can be calculated from the Riccati equation:

For the LQR control algorithm, the Lyapunov function of the system is defined as:

First derivative of time on Eq. (40), the equation can be written as:

It is noted that the Lyapunov function of the system is negative, the system is asymptotically stable, and the critical flutter point can disappear.