1 Introduction

In recent years, piezoelectric elements have been commonly concerned because of their wide applications in sensors, actuators, transducers, and some micro intelligent structures^[1-2]. In fact, in the fields of aviation, aerospace, and automation, the working environment of equipment is harsh, and is always sensitive to the temperature change, which has a great effect on the implementation of control^[3-4]. Therefore, the thermo-piezoelectric-mechanical coupling effect should be considered in the process of accurate modeling for piezoelectric elements.

The finite element method (FEM) has the merits of high precision and excellent adaptation to structures and loads. It can be used to study the vibration response of a continuous medium conveniently^[5-9]. Huang and Shen^[10] obtained the vibration characteristics of functionally graded material plates under different thermo-piezoelectric-mechanical combined loads by using the FEM. They analyzed the effects of temperature, control voltage, and material parameters on the system vibration. Compared with quantitative research, the FEM has some shortcomings in the qualitative description. The harmonic balance method (HBM) does not depend on the time domain response. Therefore, it is convenient to carry out a qualitative study on the effects of different parameters. Another advantage of the HBM is that this technique is suitable for analyzing the steady-state response of strongly nonlinear systems.

The primary resonance of piezoelectric ceramic composites is a classical practical commonly concerned engineering issue^[11-13]. The nonlinear resonance problems of beams, plates, and shell structures with magneto-electric, electro-elastic, thermo-elastic, and thermo-magnetic couplings are studied, respectively^[14-23]. Based on Hamilton's principle and the Rayleigh-Ritz method, the primary resonance of a piezoelectric ceramic thin circular plate was studied in Ref. [24]. With the method of multiple scales, the primary resonance of the laminated disk was investigated. It was shown that the primary resonance had multiple coexisting solutions and jump phenomena when the excitation frequency varied from a low value to a high one. The realization of the primary resonance lay in the stability and initial conditions.

From the current research results, the works of thermo-piezoelectric-mechanical coupling based on the FEM mainly are about the transient response. The studies in the primary resonance analysis and the steady-state qualitative analysis of the piezoelectric rectangular thin plates are not enough.

Therefore, a piezoelectric rectangular thin plate model with the consideration of the effect of the temperature field is established, and the nonlinear governing equation of the piezoelectric rectangular thin plate is obtained by the energy method. Then, the HBM is applied to obtain the first-order approximate response. The effects of the temperature difference, the damping, the plate thickness, the excited charge, and the mode on the primary resonance response of the system are discussed. Finally, through the numerical response under the sweeping frequency excitation, a theoretical analysis for the nonlinear dynamics behaviors, e.g., the hardening nonlinearity, the jumped solution, and the coexistence of multiple solutions, are verified.

2 System modeling and governing equation of nonlinear vibration

Figure 1 shows the rectangular piezoelectric plate with the substrate material No. 45 steel. The length and the width of the plate are L_x and L_y, respectively. The thicknesses of the substrate material and the piezoelectric material are h_m and h_p, respectively. Both the upper and lower layers are vertically polarized along the z-axis. The piezoelectric rectangular thin plate is supposed to be symmetric with respect to the xOy-plane. It should be noted that the applied electric field is considered as an external excitation and it exists in the z-axis direction.

For the substrate material, the constitutive equations can be expressed as follows:

(1)

where (ε_x, ε_y, γ_xy) and (σ_x, σ_y, τ_xy) are the strain and stress vectors, respectively. c_ij^m (i, j = 1, 2, 6) are the stiffness coefficients of the substrate material. The above equations are for the orthotropic plate case. For the corresponding isotropic plate, we set

where E_m and v_m are Young's modulus and Poisson's ratio, respectively.

Considering the thermal effects, the constitutive equations for the piezoelectric material can be given as follows^[25]:

(2)

where c_ij^p (i, j = 1, 2, 6) are the stiffness coefficients of the piezoelectric material. σ_i and τ_ij are the normal and shear stresses, respectively, while ε_i and γ_ij are the normal and shear strains, respectively. e₃₁ and ε₃₃ are the piezoelectric constant and the dielectric permittivity, respectively. λ is the thermo-mechanical coupling constant expressed as λ = (c₁₁^p + c₁₂^p)α^[26-27], where α is the coefficient of linear thermal expansion. d₃ is the thermal-piezoelectric coupling constant. E_z is the electric field component. D_z is the electric displacement component. η is the entropy. a_T is defined as a_T = c_E/T₀, where c_E is the heat capacity. θ is the change of temperature from the initial temperature T₀.

Based on the generalized Hamilton's principle, we can get the following variational equation:

(3)

where T, U, W_e, and Θ are the kinetic energy, the elastic potential energy, the electric energy, and the thermal energy of the system, respectively. W is the work done by the electric field force^[28]. For the plate, the expressions of T, U, W_e, Θ, and W are as follows:

(4)

where ρ_i (i = m, p) are the densities of the materials. V_i (i = m, p) are the volumes occupied by the materials. The subscripts m and p stand for the substrate layer and the piezoelectric layer, respectively. is the vertical velocity. S = (ε_x, ε_y, γ_xy)^T and Q = (σ_x, σ_y, τ_xy)^T are the strain and stress vectors, respectively. E = (0, 0, E_z)^T and D = (0, 0, D_z)^T are the electric field and the electric displacement, respectively. V(t) is the voltage. q(t) is the surface charge of the piezoelectric layer.

According to the von Karman theory of large deflection for plates^[29], the nonlinear geometric relation of the strain and the displacement can be expressed as follows:

(5)

where u, v, and w denote the displacements in the x-, y-, and z-directions, respectively.

The Rayleigh-Ritz method is the most effective method to solve boundary value problems. It can reduce the order of equations and simplify the calculation. To find the governing equations of this problem, the Rayleigh-Ritz method is used, i.e.,

(6)

where φ_i (i = u, v, w) are the vibration mode functions, and X(t) is the displacement function.

According to the classical thin plate theory, it is assumed that the vibration mode function of any point (x, y, z) in the thin plate is

(7)

For the laminated piezoelectric simply-supported rectangular thin plate, the vibration mode function should satisfy the following boundary conditions:

Thus, the vibration mode function of the middle surface is assumed as follows:

(8)

where both m and n are positive integers, representing different modes of the plate, and the coefficients A and B are determined by

(9)

in which Π is the potential energy of the system defined by

(10)

In the classical plate theory, the electric potential of the piezoelectric plate is assumed to be uniformly distributed. Therefore, the electric field boundary condition can be expressed as follows:

(11)

The expressions of the kinetic energy, the elastic potential energy, the electric energy, and the thermal energy obtained from Eq. (4) are as follows:

The detailed expressions of the above coefficient are shown in Appendix A.

Substituting the above formulas into Eq. (3), we have

(12)

The voltage can be solved from the second equation of Eq. (12), and it can be written as follows:

(13)

where

We assume that the damping is viscous and is expressed as . Then, with the consideration of the effects of damping and external harmonic electric excitation, we can write the surface charge of the piezoelectric layer as q(t) = q₀ cos(ωt). Substituting Eq. (13) into the first equation of Eq. (12) yields

(14)

where

Through a transformation X = X^∗ + Y, the constant F₄ in Eq. (14) is eliminated, where X^∗ is the solution of the equation

Then, we obtain the governing equation of the nonlinear vibration system as follows:

(15)

where

3 Static and steady state dynamics 3.1 Buckling analysis

When the effects of damping and electrical external excitation are not taken into account, Fig. 2(a) shows the static bifurcation diagram of the system with the temperature difference change. In Fig. 2(a), the ordinate represents the displacement, and the abscissa represents the temperature difference. The system presents a single stable equilibrium when the temperature difference is below 6 ℃. If the temperature difference exceeds 6 ℃, the system has two stable equilibriums and one unstable equilibrium. Figure 2(b) shows the effect of the temperature difference on the natural frequency. As shown in the figure, the natural frequency attains its minimum when the temperature difference is 6 ℃, which corresponds to the critical temperature of the buckling configuration.

3.2 Primary resonance

The HBM is adopted to carry out the primary analysis. The steady-state displacement response is assumed to be as follows^[31]:

(16)

Then, the velocity and the acceleration response can be expressed as follows:

(17)

Substitute Eqs. (16) and (17) into Eq. (15), and set the harmonic coefficients of sin(ωt) and cos(ωt) to be zero. Then, we have

(18)

where r represents the displacement in the steady state, corresponding to the maximum displacement response amplitude in the numerical frequency sweeping response, and a²(t)+b²(t) = r². Based on Eq. (18), the frequency-response equation can be obtained as follows:

(19)

4 Numerical analysis

The substrate material is No. 45 steel, and the piezoelectric layer is PZT-2. The length and the width of the plate are L_x = 300 mm and L_y = 200 mm, respectively. The thickness of No. 45 steel is h_m = 1 mm, and the thickness of PZT-2 is h_p = 1 mm. The specific parameters are shown in Table 1^[32-33].

Figure 3 shows the effect of the temperature difference on the response under different modes, where the system presents the hardening nonlinear stiffness. The solid curves represent the stable results, while the dashed curves represent the unstable results. The stability of the solution can be obtained according to the Routh-Hurwitz stability criterion^[34]. Three temperature differences are selected above the critical temperature of buckling configuration. With the increase in the temperature difference, the resonant amplitude value decreases gradually, while the corresponding jumping frequency increases. As seen in this figure, when the temperature difference is 40 ℃, the piezoelectric rectangular thin plate exhibits two saddle-node bifurcations at ω/(2π) = 21.60 Hz and ω/(2π) = 17.80 Hz, corresponding to two jumps in the system response. More specifically, when the frequency increases from ω/(2π) = 17.80 Hz, the amplitude of the response increases accordingly until reaching Point A, where ω/(2π) = 21.60 Hz. At this point, the response jumps abruptly to the stable branch of the smaller amplitude. Decreasing the frequency results in a jump at a different point, i.e., Point B, where ω/(2π) = 17.80 Hz, from the lower-amplitude stable branch to the higher-amplitude one. The frequency range between the two saddle-node points is also the region of multiple coexistence solutions. When the temperature difference increases to 80 ℃, there appears two saddle-node bifurcations at ω/(2π) = 25.01 Hz and ω/(2π) = 23.05 Hz, respectively. Thus, with the increase in the temperature difference, the frequency ranges of multiple coexistence solutions have a tendency to increase. For higher order modes, the system presents the hardening nonlinearity. In addition, with the increase in the temperature difference, the resonance amplitude value and the bandwidth decrease gradually.

Figure 4(a) represents the primary resonance response of the system under different modes. The parameters are as follows: the temperature difference is θ = 40 ℃, the plate thickness is h_m = 1 mm, the damping coefficient is c = 2.5, and the excited charge is q₀ = 0.001 C. With the increase in the modal parameter, the jumping frequency increases gradually, while the resonant amplitude value and the bandwidth decrease gradually. Figure 4(b) illustrates the damping effect on the primary resonance response. It is shown that, when the system damping increases, both the resonant amplitude value and the jumping frequency are suppressed. Figure 4(c) demonstrates that when the plate thickness h_m increases gradually, both the hardening nonlinearity and the region of the multiple coexistence solutions shrink enormously. Figure 4(d) shows the primary resonance responses of the system under different excited charge levels. When the excited charge level increases, the resonant amplitude, the jumping frequency, the region of the multiple coexistence solutions, and the operational bandwidth increase.

Figure 5 shows the frequency sweep diagrams for the temperature differences of 40 ℃ and 60 ℃. The solid bold curves are the steady-state displacement response obtained from the HBM, and the solid plain curves represent the numerical responses obtained by the Runge-Kutta method. The parameters are given in Table 1. The green curves represent the forward frequency sweep responses, while the blue curves represent the backward frequency sweep responses. When the temperature difference is 40 ℃, the jump-down frequency appears at 21.60 Hz, and the jump-up frequency appears at 17.80 Hz. Thus, the region of multiple coexistence solutions appears from 17.80 Hz to 21.60 Hz. In the processes of forward frequency sweeping and backward frequency sweeping, the large responses are generated with the jump emergence, and the amplitudes of the forward sweeping and backward sweeping match well with the analytical results from the HBM. When the temperature difference increases to 60 ℃, the region of multiple coexistence solutions moves to the region from 20.50 Hz to 23.31 Hz, which means that the bandwidth becomes narrower.

Figure 6 represents the frequency sweep diagrams with the damping coefficients c = 3.5 and c = 4.5. It can be seen that the damping coefficient c = 3.5 leads to higher peaks and broader bandwidth. When the damping coefficient increases gradually, both the resonance amplitude and the bandwidth decreases gradually.

Figure 7 shows the frequency sweep diagrams of the plate thicknesses h_m = 0.5 mm and h_m = 2 mm. It is interesting to find out that the thickness change largely affects the jumping frequency, the amplitude, and the bandwidth. When the plate thickness is h_m = 0.5 mm, the jump-down frequency appears at 23.69 Hz, and the jump-up frequency appears at 17.3 Hz. The frequency bandwidth of multiple coexisting solutions appears from 17.3 Hz to 23.69 Hz. When the plate thickness increases to h_m = 2.0 mm, the frequency bandwidth of multiple coexisting solutions appears from 18.32 Hz to 20.33 Hz. Thus, a thinner plate is better for the actuation in the low frequency range and the operation with a wider bandwidth.

Figure 8 represents the frequency sweep diagrams for q₀ = 0.000 5 C and q₀ = 0.001 5 C. A nonlinear system with distinct hardening and jump phenomena is observed in the sinusoidal sweep analysis. Stronger excited charges are beneficial for eliciting greater responses and generating a wider region of co-existing solutions. Both jump-up and jump-down frequencies have a tendency to increase when the excited charge increases. However, this tendency is only noticeable for the jump-down frequency that climbs from 18.40 Hz to 23.86 Hz when the excited charge is boosted up from 0.000 5 C to 0.001 5 C. For the jump-up point, where the system jumps from the low-energy branch to the high-energy branch, the increase is trivial (17.30 Hz–18.07 Hz in the simulation). The jump-up frequency, unlike its jump-down counterpart, is not stable either. Take Fig. 8(b) for example. When the excitation is set as q₀ = 0.001 5 C with a frequency of 23.86 Hz, the displacement response amplitude is 2.08 mm for the high-energy branch and 0.02 mm for the low-energy branch. Whether it is a high-energy response or a low-energy state physical response is determined by the initial conditions.

5 Conclusions

In this paper, according to the von Karman theory of large deflection, the governing equation of a rectangular thin plate with thermo-piezoelectric-mechanical coupling is established. The HBM is used to study the first-order primary resonance response of the system. The effects of the parameters, e.g., the temperature difference, the damping coefficient, the plate thickness, and the excited charge, on the system response are analyzed by the numerical method. The conclusions are summarized as follows:

(ⅰ) The critical temperature difference of the thermal buckling is obtained through the static bifurcation analysis. When the temperature difference increases, transcritical pitchfork bifurcation occurs, and the number of equilibriums increases from one to three.

(ⅱ) The governing equation of the rectangular thin plate has a hardening nonlinearity, which obviously affects the primary resonance response of the structure. When the excited charge level increases, both the resonance peak and the bandwidth increase gradually. However, when the damping coefficient increases, both the resonance amplitude and the bandwidth decrease gradually.

(ⅲ) With the increase in the temperature difference, the jumping frequency increases, while the resonant amplitude decreases gradually. With the increase in the plate thickness, both the hardening nonlinearity and the regions of multiple coexisting solutions shrink enormously.

(ⅳ) All the above conclusions are based on fundamental modes. For high order modes, the system exhibits hardening nonlinearity. In addition, the resonance amplitude value and the bandwidth decrease gradually because the deformation energy mainly distributes within the low frequency range.

Appendix A

where