Nonlinear dynamics of a piezoelectric precision drive system used for robot joint

Abstract

To improve the robot joint’s control precision, a piezoelectric precision drive system is proposed. The proposed system has merits of fast response and compact size. Using continuous dynamic theory, the nonlinear coupled dynamic models and equations of the drive system are established. Applying Linz Ted-Poincaré Method, the approximate solutions of nonlinear solutions are solved. The magnitude-frequency characteristics are also analyzed. What’s more, the effects of nonlinear piezoelectric characteristic on the system’s dynamic performance are investigated yet. Results show that the jump phenomenon of vibration amplitude occurs in nonlinear resonance. Meanwhile, the nonlinear piezoelectric characteristic influences the dynamic performance of the drive system obviously. These results can be used to predict the dynamic load of the piezoelectric precision drive system.

Keywords

Nonlinear dynamic magnitude-frequency characteristic dynamic response piezoelectric drive system

1 Introduction

Piezoelectric smart material has many excellent merits, such as small size, fast response and no electromagnetic interference. As a result, it has been successfully utilized in the field such as precision actuating [1], accurate positioning [2] and energy harvesting [3].

In the area of robot, a large number of studies have been carried out. O. U. David et al. [4] proposed an enhanced tactile sensor with a piezoelectric bimorph. As the sensor can differentiate soft materials with similar mechanical characteristics, so the tactile sensor was designed for brain tumour resection. B. Ghosh et al. [5] used a bimorph piezoelectric actuator to provide the dexterous behavior during robotic assembly. J. Hadi et al. [6] designed a micro-robot in stick-slip mode, and the robot was driven by two piezoelectric actuators with harmonic signal. Y. Tajitsu [7] proposed a piezoelectric Poly-L-lactic acid fabric, and it could feel the human complex movement. Hence, the fabric was utilized for controlling humanoid robot. What’s more, a steerable miniature legged robot driven by a single piezoelectric bending unimorph actuator was designed and fabricated by A. Dharmawan et al. [8], and it could do the actions such as move forward, turn right and turn left. Y. Zou et al. [9] proposed an insect-scale flapping-wing robot with a wingspan of 35 mm. The robot was actuated by piezoelectric actuators and it can generate sufficient thrust to take off with a flapping amplitude approximately ±60° under the resonant wingbeat frequency of 100 Hz. Besides, J. Ma et al. [10] designed a micro-robot system used for abdominal surgery. In the robot system, the piezoelectric actuators were used to drive flexible hinge to accomplish the joint motions.

Thus it can be seen that piezoelectric smart material has been widely utilized in robot area. However, in robot joint aspect, relatively less studies of piezoelectric driving have been investigated now. Therefore, the authors [11] proposed a piezoelectric precision drive system used for robot joint.

In this paper, the nonlinear dynamic models and equations of the piezoelectric precision drive system used for the robot joint are established with continuous system dynamics theory. Using Linz Ted-Poincaré method, the approximate solutions of the nonlinear vibration equations are presented. The results can be used to reduce the vibration and predict the dynamic load of the piezoelectric precision drive system.

2 Operating principle and applications

Figure 1 shows the structure of the piezoelectric precision drive system. Compare with traditional piezoelectric drive, the proposed drive system is a new integrated system. It uses rolling contact of movable tooth to replace the frictional contact of traditional piezoelectric system. Thus the novel system has the advantage of low speed, large output torque, long operating life and no electromagnetic interference.

The driving source of the piezoelectric precision drive system consists of two piezoelectric actuators with a phase difference of 90 degrees. When the piezoelectric actuators are subject to cosine signals with a positive bias and 90 degrees phase difference, they generate an axial elongation. With z-shaped rod, the swaying rod is pushed toward the side of the springs. When the voltage of the exciting signals returns to zero, the swaying rod returns to the initial position under the elastic force of the springs. Then, the swaying rod swings back and forth in two directions performing a continuous harmonic wave on wave generator. Here, 30 movable teeth and 29 rigid cog teeth are used. With the effect of movable tooth drive, the rotation can achieve a transmission ratio of 30. Here, the movable teeth are steel balls and its diameter is 1 mm. When working, the movable teeth locate in the groove of the rotor. The rigid cog is formed by tooth profile equation, and the equation can be expressed by $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}X=b\cos \left({\theta}-\arcsin \left[a\sin \left(\left(i_{cp}-1\right){\theta}\right)/b\right]\right)+a\cos \left(i_{cp}{\theta}\right)\pm r_{p}\cos {\psi}\\[3.60004pt] Y=b\sin \left({\theta}-\arcsin \left[a\sin \left(\left(i_{cp}-1\right){\theta}\right)/b\right]\right)+a\sin \left(i_{cp}{\theta}\right)\pm r_{p}\sin {\psi}\end{array}\right. & & \displaystyle\end{eqnarray}$ (1) where θ is rotational angle of the rotor; 𝜓 is the angle between x axis and tangent of central trajectory of movable tooth; b is coefficient, b = r _s + r _p; r _s and r _p are radius of wave generator and movable tooth.

As mentioned above, the proposed piezoelectric precision drive system can achieve the functions of fast response, small speed and large torque, so it can be used in small robot, precision positioning system and lunar probe. In aerospace field, Zhao Chunsheng and his partners [12] applied their research results to the ‘Chang’e 3’ spacecraft. In this paper, the proposed system is just in experimental stage, and it can be used in robot and other fields after extensively testing.

The proposed piezoelectric precision drive system has characteristics of compact structure, large output torque, quick response long lifespan and no electromagnetic interference. Hence, it can work at extreme conditions such as magnetic environments and quiet occasions. This is the reason why we choose this piezoelectric precision drive system for robot use.

For the piezoelectric precision drive system, the free vibration and forced vibration of the driving system for the drive system was investigated [13,14]. The proposed piezoelectric drive system is a coupled system, and the coupled dynamics affect the overall behavior of the robot joint system. In Ref. [15], the coupled free vibration of the system was investagated. However, the nonlinear output characteristic of piezoelectric actuator influence the system’s dynamic performance, which will decrease the load-carrying capacity of the robot joint system. Moreover, the nonlinear coupled vibrations of the piezoelectric precision drive system have not been investagated yet. So, the study of the nonlinear coupled vibration for the proposed piezoelectric drive system is necessary.

3 Nonlinear coupled dynamic equations

Figure 2 shows dynamics model of the piezoelectric precision drive system, where the z-shaped rod and swaying rod constitute displacement amplification mechanism. What’s more, the z-shaped rod consists of OD part and DB part. The driving system is a coupled dynamic system. And the exciting forces of the displacement amplification mechanism come from piezoelectric actuator. Due to the existence of hysteresis effect, the output strain and stress of piezoelectric actuator are nonlinear variation along with the electric field intensity. Applying nonlinear piezoelectric theory [16], the nonlinear strain of the piezoelectric actuator can be written as $\begin{eqnarray}\displaystyle S_{3}=d_{33}E_{3}+\frac{1}{2}d_{333}E_{3}^{2} & & \displaystyle\end{eqnarray}$ (2) where S ₃, d ₃₃, E ₃ and d ₃₃₃ are the strain, strain constant, electric field strength and nonlinear piezoelectric strain constant of the piezoelectric actuator, respectively. As the piezoelectric actuator is good at resisting compression instead of tension, so the exciting signal of the system is sinusoidal signal with positive bias, so E ₃ = V _p−p[cos(ω_e t) +1]∕(2 l _p), here V _p−p is peak--peak voltage, ω_e is exciting frequency, and l _p is the thickness of per piezoelectric piece.

From Eq. (2), the outer force of the piezoelectric actuator under electric field strength E ₃ can be given $\begin{eqnarray}\displaystyle F_{p}=A_{p}c_{33}\left(d_{33}E_{3}+\frac{1}{2}d_{333}E_{3}^{2}\right) & & \displaystyle\end{eqnarray}$ (3) where A _p and c ₃₃ are the cross section area and elastic stiffness coefficient of the piezoelectric actuator.

Considering soft spring performance of the piezoelectric material, the relationship between the strain and the stress of the piezoelectric actuator is $\begin{eqnarray}\displaystyle {\sigma}_{p}=c_{33}\left({\varepsilon}_{y}-{\varepsilon}_{1}{\varepsilon}_{y}^{2}-{\varepsilon}_{1}{\varepsilon}_{y}^{3}\right) & & \displaystyle\end{eqnarray}$ (4) where ϵ₁ is the nonlinear constant; σ_p and ϵ_y are the axial stress and strain of piezoelectric actuator, ϵ_y = ∂v∕∂y.

The total internal force under electric field force on the piezoelectric actuator can be expressed by $\begin{eqnarray}\displaystyle F & = & \displaystyle A_{p}{\sigma}_{p}+\int _{0}^{y}\left(F_{p}/l_{np}\right)dy\nonumber\\ \displaystyle & = & \displaystyle A_{p}\left[c_{33}\left({\varepsilon}_{y}-{\varepsilon}_{1}{\varepsilon}_{y}^{2}-{\varepsilon}_{1}{\varepsilon}_{y}^{3}\right)\right]+\int _{0}^{y}\left[A_{p}c_{33}\left(d_{33}E_{3}+d_{333}E_{3}^{2}/2\right)/l_{np}\right]dy\nonumber\\ \displaystyle & = & \displaystyle A_{p}\left[c_{33}\left({\varepsilon}_{y}-{\varepsilon}_{1}{\varepsilon}_{y}^{2}-{\varepsilon}_{1}{\varepsilon}_{y}^{3}\right)\right]+{\varepsilon}_{2}\int _{0}^{y}F_{p0}dy+{\varepsilon}_{2}\int _{0}^{y}F_{p1}dy\cos \left({\omega}t\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +{\varepsilon}_{2}\int _{0}^{y}F_{p2}dy\cos \left(2{\omega}t\right)\end{eqnarray}$ (5) where ϵ₂ is the nonlinear constant, ϵ₂ = l _p∕l _np, as ϵ₁ and ϵ₂ are both small parameters, here letting ϵ = ϵ₁ = ϵ₂; l _np are the total length of piezoelectric actuator; F _p0, F _p1 and F _p2 are the coefficients $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}F_{p0}=\displaystyle \frac{c_{33}V_{p-p}\left(8l_{p}d_{33}+3d_{333}V_{p-p}\right)}{16l_{p}^{2}l_{np}}\\[9.60004pt] F_{p1}=\displaystyle \frac{c_{33}V_{p-p}\left(2l_{p}d_{33}+d_{333}V_{p-p}\right)}{4l_{p}^{2}l_{np}}\\[9.60004pt] F_{p2}=\displaystyle \frac{c_{33}d_{333}V_{p-p}^{2}}{16l_{p}^{2}l_{np}}\end{array}\right. & & \displaystyle\end{eqnarray}$ (6) As the piezoelectric actuator vibrates in its axis direction, so the nonlinear dynamic equation of the piezoelectric actuator can be written as $\begin{eqnarray}\displaystyle {\rho}_{p}A_{p}\frac{\partial ^{2}v}{\partial t^{2}}=A_{p}c_{33}\frac{\partial ^{2}v}{\partial y^{2}}\left(1-2{\varepsilon}\frac{\partial v}{\partial y}-3{\varepsilon}\left(\frac{\partial v}{\partial y}\right)^{2}\right)+{\varepsilon}\left(F_{p0}+F_{p1}+F_{p2}\right) & & \displaystyle\end{eqnarray}$ (7)

The design purpose of the drive system is to reduce the volume and external diameter under the precondition of satisfying the magnification factor requirement. After magnifying of displacements, the beam’s depth to length ratio is smaller than 5, so the Timoshenko beam model is used for OD beam of z-shaped rod, while Bernoulli model is used for the DB beam and swaying rod, as the same ratios for DB beam and swaying rod are larger than 5. So, the dynamic equations of displacement amplification mechanism can be written as $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\rho}S_{1}\frac{\partial ^{2}y_{1}}{\partial t^{2}}+EI\frac{\partial ^{4}y_{1}}{\partial x_{1}^{4}}-{\rho}I\frac{\partial ^{4}y_{1}}{\partial x_{1}^{2}\partial t^{2}}=f_{od}\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle EI\frac{\partial ^{4}y_{2}(x_{2},t)}{\partial x_{2}^{4}}+{\rho}\frac{\partial ^{2}y_{2}(x_{2},t)}{\partial t^{2}}=f_{db}\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle EI_{T}\frac{\partial ^{4}u(y,t)}{\partial ^{4}y}+{\rho}\frac{\partial ^{2}u(y,t)}{\partial t^{2}}=f_{s}\end{eqnarray}$ (10) where $\rho$ is density of displacement amplification mechanism; S ₁ is cross-sectional of z-shaped rod; I and I _T are inertia moment of z-shaped rod and swaying rod.

Meanwhile, the nonlinear exciting forces of displacement amplification mechanism can be given as $\begin{eqnarray}\displaystyle f_{od}=\frac{{\varepsilon}}{l_{np}}\int _{0}^{l_{np}}F\dot{{\phi}}\left(y\right)dy & & \displaystyle\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle f_{db}=\frac{{\varepsilon}l_{1}}{l_{np}l_{3}}\int _{0}^{l_{np}}F\dot{{\phi}}\left(y\right)dy & & \displaystyle\end{eqnarray}$ (12) $\begin{eqnarray}\displaystyle f_{s}=\frac{{\varepsilon}l_{1}}{l_{np}l_{3}}\int _{0}^{l_{np}}F\dot{{\phi}}\left(y\right)dy & & \displaystyle\end{eqnarray}$ (13) where 𝜙(y) is modal functions of piezoelectric actuator; variables l ₁ and l ₃ are length parameters (see Fig. 3).

The displacement expressions can be written as $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}v\left(y,t\right)=\left(y\right)q_{p}\left(t\right)\\ y_{1}\left(x_{1},t\right)={\varphi}\left(x_{1}\right)q_{od}\left(t\right)\\ y_{2}\left(x_{2},t\right)={\eta}\left(x_{2}\right)q_{db}\left(t\right)\\ u\left(y,t\right)={\psi}\left(y\right)q_{s}\left(t\right)\end{array}\right. & & \displaystyle\end{eqnarray}$ (14) Substituting Eq. (14) into Eqs. (7)–(10), the nonlinear coupled dynamic equations of the motor can be expressed by $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\ddot{q}_{p}\left(t\right)+b_{11}q_{p}\left(t\right)+b_{12}{\varepsilon}q_{p}^{2}\left(t\right)+b_{13}{\varepsilon}q_{p}^{3}\left(t\right)={\varepsilon}F_{10}+{\varepsilon}F_{11}\cos {\omega}t+{\varepsilon}F_{12}\cos 2{\omega}t\\ \ddot{q}_{od}\left(t\right)+b_{21}q_{od}\left(t\right)={\varepsilon}F_{20}+{\varepsilon}F_{21}\cos {\omega}t+{\varepsilon}F_{22}\cos 2{\omega}t\\ \ddot{q}_{db}\left(t\right)+b_{31}q_{db}\left(t\right)={\varepsilon}F_{30}+{\varepsilon}F_{31}\cos {\omega}t+{\varepsilon}F_{32}\cos 2{\omega}t\\ \ddot{q}_{s}\left(t\right)+b_{41}q_{s}\left(t\right)={\varepsilon}F_{40}+{\varepsilon}F_{41}\cos {\omega}t+{\varepsilon}F_{42}\cos 2{\omega}t\end{array}\right. & & \displaystyle\end{eqnarray}$ (15) where 𝜙(y), 𝜑(x ₁), 𝜂(x ₁) and 𝜓(y) are modal functions of piezoelectric actuator, OD beam, DB beam and swaying rod. b _ij is coefficients of dynamic equations. Here $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\bar{{\phi}}=\displaystyle \frac{1}{l_{np}}\displaystyle \int _{0}^{l_{np}}\sum {\phi}_{i}\left(y\right)dy\\[12.0pt] \bar{{\varphi}}=\displaystyle \frac{1}{l_{2}}\displaystyle \int _{0}^{l_{2}}\sum {\varphi}_{i}\left(x_{1}\right)dx_{1}\\[12.0pt] \bar{{\eta}}=\displaystyle \frac{1}{l_{3}}\displaystyle \int _{0}^{l_{3}}\sum {\eta}_{i}\left(x_{2}\right)dx_{2}\\[12.0pt] \bar{{\psi}}=\displaystyle \frac{1}{l_{8}}\displaystyle \int _{0}^{l_{8}}\sum {\psi}_{i}\left(y\right)dy\end{array}\right.\end{eqnarray}$ (16) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}b_{11}=\displaystyle -\frac{c_{33}\bar{{\phi}^{\prime \prime }}}{{\rho}_{p}\bar{{\phi}}}\quad b_{12}=\frac{2c_{33}\bar{{\phi}^{\prime }}\bar{{\phi}^{\prime \prime }}}{{\rho}_{p}\bar{{\phi}}}\quad b_{13}=\frac{3c_{33}\bar{{\phi}^{\prime }}^{2}\bar{{\phi}^{\prime \prime }}}{{\rho}_{p}\bar{{\phi}}}\\[12.0pt] b_{21}=\displaystyle \frac{EI\bar{{\phi}}^{(4)}}{{\rho}\left(S_{1}\bar{{\phi}}-I\bar{{\phi}}^{\prime \prime }\right)}\quad b_{31}=\frac{EI\bar{{\eta}}^{(4)}}{{\rho}\bar{{\eta}}}\quad b_{41}=\frac{EI_{T}\bar{{\psi}}^{(4)}}{{\rho}\bar{{\psi}}}\end{array}\right.\end{eqnarray}$ (17) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}{\Lambda}_{od}=\displaystyle \frac{{\rho}_{p}\bar{{\phi}}}{{\rho}\left(S_{1}\bar{{\varphi}}-I\bar{{\varphi}^{\prime \prime }}\right)}\left(\displaystyle \int _{0}^{l_{np}}\dot{{\phi}}dy\right)\\[12.0pt] {\Lambda}_{db}=\displaystyle \frac{{\rho}_{p}\bar{{\phi}}l_{1}}{{\rho}\bar{{\eta}}l_{3}}\left(\displaystyle \int _{0}^{l_{np}}\dot{{\phi}}dy\right)\\[12.0pt] {\Lambda}_{s}=\displaystyle \frac{{\rho}_{p}\bar{{\phi}}l_{1}}{{\rho}\bar{{\psi}}l_{3}}\left(\displaystyle \int _{0}^{l_{np}}\dot{{\phi}}dy\right)\end{array}\right.\end{eqnarray}$ (18) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}F_{10}=\displaystyle \frac{F_{p0}}{{\rho}_{p}\bar{{\phi}}}\quad F_{11}=\frac{F_{p1}}{{\rho}_{p}\bar{{\phi}}}\quad F_{12}=\frac{F_{p2}}{{\rho}_{p}\bar{{\phi}}}\\[9.60004pt] F_{20}={\Lambda}_{od}F_{10}\quad F_{21}={\Lambda}_{od}F_{11}\quad F_{22}={\Lambda}_{od}F_{12}\\ F_{30}={\Lambda}_{db}F_{10}\quad F_{31}={\Lambda}_{db}F_{11}\quad F_{32}={\Lambda}_{db}F_{12}\\ F_{40}={\Lambda}_{s}F_{10}\quad F_{41}={\Lambda}_{s}F_{11}\quad F_{42}={\Lambda}_{s}F_{12}\end{array}\right.\end{eqnarray}$ (19) From Eq. (15), the matrix form nonlinear dynamic equations of the drive system can be given by $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}\left(\boldsymbol{t}\right)+\boldsymbol{B}_{\mathbf{1}}\boldsymbol{Q}\left(\boldsymbol{t}\right)+{\varepsilon}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}^{\mathbf{2}}\left(\boldsymbol{t}\right)+{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}^{\mathbf{3}}\left(\boldsymbol{t}\right)={\varepsilon}\boldsymbol{F}_{\mathbf{0}}+{\varepsilon}\boldsymbol{F}_{\mathbf{1}}\cos \boldsymbol{\it\omega}t+{\varepsilon}\boldsymbol{F}_{\mathbf{2}}\cos 2\boldsymbol{\it\omega}t & & \displaystyle\end{eqnarray}$ (20) where $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{Q}=\left[\begin{array}{@{}cccc@{}}q_{p} & q_{od} & q_{db} & q_{s}\\ \end{array}\right]^{T}\end{eqnarray}$ (21) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\boldsymbol{F}_{\mathbf{0}}=\left[\begin{array}{@{}cccc@{}}F_{10} & F_{20} & F_{30} & F_{40}\\ \end{array}\right]^{T}\\ \boldsymbol{F}_{\mathbf{1}}=\left[\begin{array}{@{}cccc@{}}F_{11} & F_{21} & F_{31} & F_{41}\\ \end{array}\right]^{T}\\ \boldsymbol{F}_{\mathbf{2}}=\left[\begin{array}{@{}cccc@{}}F_{12} & F_{22} & F_{32} & F_{42}\\ \end{array}\right]^{T}\end{array}\right.\end{eqnarray}$ (22) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\boldsymbol{B}_{\mathbf{1}}=\text{diag}\left(\left[\begin{array}{@{}cccc@{}}b_{11} & b_{21} & b_{31} & b_{41}\\ \end{array}\right]\right)\\ \boldsymbol{B}_{\mathbf{2}}=\text{diag}\left(\left[\begin{array}{@{}cccc@{}}b_{12} & b_{22} & b_{32} & b_{42}\\ \end{array}\right]\right)\\ \boldsymbol{B}_{\mathbf{3}}=\text{diag}\left(\left[\begin{array}{@{}cccc@{}}b_{13} & b_{23} & b_{33} & b_{43}\\ \end{array}\right]\right)\end{array}\right.\end{eqnarray}$ (23)

4 Solution of the nonlinear dynamic equations

4.1 Nonlinear free vibration

When exciting voltage V _p−p is zero and small parameter ϵ is small enough, the system changes into nonlinear free vibration system. So, the coupled dynamic equation can be written as $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}\left(\boldsymbol{t}\right)+\boldsymbol{B}_{\mathbf{1}}\boldsymbol{Q}\left(\boldsymbol{t}\right)+{\varepsilon}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}^{\mathbf{2}}\left(\boldsymbol{t}\right)+{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}^{\mathbf{3}}\left(\boldsymbol{t}\right)=\mathbf{0} & & \displaystyle\end{eqnarray}$ (24)

Using Linz Ted-Poincaré method [17] to solve Eq. (20), letting the solution Q ( t , 𝜺) and natural frequency 𝝎 expand into power series $\begin{eqnarray}\displaystyle \boldsymbol{Q}\left(\boldsymbol{t},\boldsymbol{\it\varepsilon}\right) & = & \displaystyle \boldsymbol{Q}_{\mathbf{0}}\left(\boldsymbol{t}\right)+{\varepsilon}\boldsymbol{Q}_{\mathbf{1}}\left(\boldsymbol{t}\right)+{\varepsilon}^{2}\boldsymbol{Q}_{\mathbf{2}}\left(\boldsymbol{t}\right)+\cdots\end{eqnarray}$ (25) $\begin{eqnarray}\displaystyle \boldsymbol{\it\omega}^{2} & = & \displaystyle \boldsymbol{\it\omega}_{\mathbf{0}}^{2}\left(1+{\varepsilon}\boldsymbol{\it\sigma}_{\mathbf{1}}+{\varepsilon}^{2}\boldsymbol{\it\sigma}_{\mathbf{2}}+\cdots \right)\end{eqnarray}$ (26) where 𝝎 and 𝝎₀ are nonlinear and linear natural frequency, respectively.

Introducing a variable 𝜿 = 𝝎t, and substituting 𝜅 and Eqs. (25), (26) into Eq. (24), then the approximate linear equations can be obtained $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{0}}+\boldsymbol{Q}_{\mathbf{0}} & = & \displaystyle \mathbf{0}\end{eqnarray}$ (27) $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{1}}+\boldsymbol{Q}_{\mathbf{1}} & = & \displaystyle -\boldsymbol{\it\sigma}_{\mathbf{1}}\ddot{\boldsymbol{Q}}_{\mathbf{0}}-\frac{1}{\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{2}}-\frac{1}{\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{3}}\end{eqnarray}$ (28) $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{2}}+\boldsymbol{Q}_{\mathbf{2}} & = & \displaystyle -\boldsymbol{\it\sigma}_{\mathbf{2}}\ddot{\boldsymbol{Q}}_{\mathbf{0}}-\boldsymbol{\it\sigma}_{\mathbf{1}}\ddot{\boldsymbol{Q}}_{\mathbf{1}}-\frac{2}{\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}_{\mathbf{0}}\boldsymbol{Q}_{\mathbf{1}}-\frac{3}{\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{2}}\boldsymbol{Q}_{\mathbf{1}}\end{eqnarray}$ (29) Setting initial conditions as follows $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\boldsymbol{Q}_{\mathbf{0}}\left(0\right)=\boldsymbol{\it\delta}_{\mathbf{0}},\quad \dot{\boldsymbol{Q}}_{\mathbf{0}}\left(0\right)=\mathbf{0}\\ \boldsymbol{Q}_{\mathbf{1}}\left(0\right)=\mathbf{0},\quad \dot{\boldsymbol{Q}}_{\mathbf{1}}\left(0\right)=\mathbf{0}\\ \boldsymbol{Q}_{\mathbf{2}}\left(0\right)=\mathbf{0},\quad \dot{\boldsymbol{Q}}_{\mathbf{2}}\left(0\right)=\mathbf{0}\end{array}\right. & & \displaystyle\end{eqnarray}$ (30) where 𝜹₀ is initial displacement of the system.

From Eqs. (27) and (30), the zero order approximate solution can be obtained $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{0}}\left(\boldsymbol{t}\right)=\boldsymbol{\it\delta}_{\mathbf{0}}\cos \boldsymbol{\it\kappa} & & \displaystyle\end{eqnarray}$ (31)

Substituting Eq. (31) into Eq. (28), the first order approximate solution can be given as $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{1}}\left(\boldsymbol{t}\right)=-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{2\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{2}}+\left(\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{3}}\right)\cos \boldsymbol{\it\kappa}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{6\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{2}}\cos 2\boldsymbol{\it\kappa}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{3}}\cos 3\boldsymbol{\it\kappa} & & \displaystyle\end{eqnarray}$ (32)

Substituting Eq. (30) and (31) into Eq. (29), the second order approximate equation can be written as $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{2}}+\boldsymbol{Q}_{\mathbf{2}}=\mathbf{Z}_{\mathbf{0}}+\left(\boldsymbol{\it\sigma}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}+\mathbf{Z}_{\mathbf{1}}\right)\cos \boldsymbol{\it\kappa}+\mathbf{Z}_{\mathbf{2}}\cos 2\boldsymbol{\it\kappa}+\mathbf{Z}_{\mathbf{3}}\cos 3\boldsymbol{\it\kappa}+\mathbf{Z}_{\mathbf{4}}\cos 4\boldsymbol{\it\kappa}+\mathbf{Z}_{\mathbf{5}}\cos 5\boldsymbol{\it\kappa} & & \displaystyle\end{eqnarray}$ (33) where $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\mathbf{Z}_{\mathbf{0}}=\displaystyle -\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}-\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{32\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}+\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{4\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\\[12.0pt] \mathbf{Z}_{\mathbf{1}}=\displaystyle \frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{128\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{2\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}\\[12.0pt] \mathbf{Z}_{\mathbf{2}}=\displaystyle \frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{4\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{4\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\\[12.0pt] \mathbf{Z}_{\mathbf{3}}=\displaystyle \frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{16\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{6\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{4\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}\\[12.0pt] \mathbf{Z}_{\mathbf{4}}=-\displaystyle \frac{5\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{32\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}\\[12.0pt] \mathbf{Z}_{\mathbf{5}}=-\displaystyle \frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{128\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\end{array}\right. & & \displaystyle\end{eqnarray}$ (34)

The second order approximate solution can be obtained by solving Eq. (33) $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{2}}\left(\boldsymbol{t}\right)=\mathbf{Z}_{\mathbf{0}}-\mathbf{Z}_{\mathbf{6}}\cos \boldsymbol{\it\kappa}-\frac{1}{3}\mathbf{Z}_{\mathbf{2}}\cos 2\boldsymbol{\it\kappa}-\frac{1}{8}\mathbf{Z}_{\mathbf{3}}\cos 3\boldsymbol{\it\kappa}-\frac{1}{15}\mathbf{Z}_{\mathbf{4}}\cos 4\boldsymbol{\it\kappa}-\frac{1}{24}\mathbf{Z}_{\mathbf{5}}\cos 5\boldsymbol{\it\kappa} & & \displaystyle\end{eqnarray}$ (35) where $\begin{eqnarray}\displaystyle \mathbf{Z}_{\mathbf{6}}=\frac{67\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{144\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{6\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}-\frac{5\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{6\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{23\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{1024\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}} & & \displaystyle\end{eqnarray}$ (36)

So, from Eqs. (31), (32) and (35), the approximate solution of the nonlinear free vibration can be expressed by $\begin{eqnarray}\displaystyle \boldsymbol{Q}\left(\boldsymbol{t},\boldsymbol{\it\varepsilon}\right) & = & \displaystyle \boldsymbol{\it\delta}_{\mathbf{0}}\cos \boldsymbol{\it\kappa}+\frac{{\varepsilon}}{\boldsymbol{\it\omega}_{\boldsymbol{ 0}}^{2}}\left[-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{2}\boldsymbol{B}_{\mathbf{2}}+\left(\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{3}\boldsymbol{B}_{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32}\boldsymbol{B}_{\mathbf{3}}\right)\cos \boldsymbol{\it\kappa}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{6}\boldsymbol{B}_{\mathbf{2}}\cos 2\boldsymbol{\it\kappa}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32}\boldsymbol{B}_{\mathbf{3}}\cos 3\boldsymbol{\it\kappa}\right]\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\frac{{\varepsilon}^{2}}{\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\left[\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{4}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{3}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}-\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{32}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}+\left(\frac{5\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{6}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}+\frac{23\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{1024}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{67\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{144}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{6}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}\right)\right.\cos \boldsymbol{\it\kappa}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\left(\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{9}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{12}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{12}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\right)\cos 2\boldsymbol{\it\kappa}+\left(\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{48}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{32}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}-\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{128}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\right)\cos 3\boldsymbol{\it\kappa}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -∼\left.\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{96}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}\cos 4\boldsymbol{\it\kappa}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{1024}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\cos 5\boldsymbol{\it\kappa}\right]\end{eqnarray}$ (37)

To avoid secular terms, the coefficients for cos 𝜿 in first order and second order approximate equations must be zero. Hence, coefficients σ₁ and σ₂ can be written as $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\boldsymbol{\it\sigma}_{\mathbf{1}}=3\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\boldsymbol{B}_{\mathbf{3}}/\left(4\boldsymbol{\it\omega}_{\mathbf{0}}^{2}\right)\\[4.79993pt] \boldsymbol{\it\sigma}_{\mathbf{2}}=\displaystyle \frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{128\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{2\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}\end{array}\right. & & \displaystyle\end{eqnarray}$ (38) Therefore, the relationship between nonlinear vibration frequencies and initial displacement can be written as $\begin{eqnarray}\displaystyle \boldsymbol{\it\omega}^{2}=\boldsymbol{\it\omega}_{\mathbf{0}}^{2}\left[1+{\varepsilon}\left(\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{4\boldsymbol{\it\omega}_{\mathbf{0}}^{2}}\boldsymbol{B}_{\mathbf{3}}\right)+{\varepsilon}^{2}\left(\frac{3\boldsymbol{\it\delta}_{\mathbf{0}}^{4}}{128\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}-\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{2\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}+\frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{3\boldsymbol{\it\omega}_{\mathbf{0}}^{4}}\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}\right)\right] & & \displaystyle\end{eqnarray}$ (39)

4.2 Nonlinear forced vibration

Eq. (20) is the forced vibration equation of the system. Here, when exciting frequency is close to natural frequency, introducing frequency-variance ${\varepsilon}\boldsymbol{\it\sigma}=\boldsymbol{\it\omega}^{2}-\boldsymbol{\it\omega}_{\mathbf{0}}^{2}$ . Substituting ϵ𝝈 and Eq. (25) into Eq. (20), the approximate linear equations can be obtained $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{0}}+\boldsymbol{\it\omega}^{2}\boldsymbol{Q}_{\mathbf{0}} & = & \displaystyle \mathbf{0}\end{eqnarray}$ (40) $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{1}}+\boldsymbol{\it\omega}^{2}\boldsymbol{Q}_{\mathbf{1}} & = & \displaystyle \boldsymbol{\it\sigma}Q_{\mathbf{0}}-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{2}}-\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{3}}+\boldsymbol{F}_{\mathbf{0}}+\boldsymbol{F}_{\mathbf{1}}\cos \boldsymbol{\it\omega}t+\boldsymbol{F}_{\mathbf{2}}\cos 2\boldsymbol{\it\omega}t\end{eqnarray}$ (41) $\begin{eqnarray}\displaystyle \ddot{\boldsymbol{Q}}_{\mathbf{2}}+\boldsymbol{\it\omega}^{2}\boldsymbol{Q}_{\mathbf{2}} & = & \displaystyle \boldsymbol{\it\sigma}Q_{\mathbf{1}}-2\boldsymbol{B}_{\mathbf{2}}\boldsymbol{Q}_{\mathbf{0}}\boldsymbol{Q}_{\mathbf{1}}-3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{Q}_{\mathbf{0}}^{\mathbf{2}}\boldsymbol{Q}_{\mathbf{1}}\end{eqnarray}$ (42) Assuming the solution of Eq. (40) is $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{0}}\left(\boldsymbol{t}\right)=\mathbf{M}_{\mathbf{0}}\cos \boldsymbol{\it\omega}t+\mathbf{N}_{\mathbf{0}}\sin \boldsymbol{\it\omega}t & & \displaystyle\end{eqnarray}$ (43)

Substituting Eq. ((43)) into Eq. ((41)), yields $\begin{eqnarray}\displaystyle &&\ddot{\boldsymbol{Q}}_{\mathbf{1}}+\boldsymbol{\it\omega}^{2}\boldsymbol{Q}_{\mathbf{1}} = \displaystyle -\frac{1}{2}\boldsymbol{B}_{\mathbf{2}}\left(\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}+\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)+\boldsymbol{F}_{\boldsymbol{ 0}}+\left[-\frac{3}{4}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{M}_{\mathbf{0}}\left(\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}+\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)+\boldsymbol{\it\sigma}M_{\boldsymbol{ 0}}+\boldsymbol{F}_{\mathbf{1}}\right]\cos \boldsymbol{\it\omega}t\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad+∼\left[-\frac{3}{4}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{N}_{\mathbf{0}}\left(\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}+\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)+\boldsymbol{\it\sigma}N_{\boldsymbol{ 0}}\right]\sin \boldsymbol{\it\omega}t+\left[-\frac{1}{2}\boldsymbol{B}_{\mathbf{2}}\left(\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}-\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)+\boldsymbol{F}_{\boldsymbol{ 2}}\right]\cos 2\boldsymbol{\it\omega}t\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad-∼\boldsymbol{B}_{\mathbf{2}}\boldsymbol{M}_{\mathbf{0}}\boldsymbol{N}_{\mathbf{0}}\sin 2\boldsymbol{\it\omega}t-\frac{1}{4}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{M}_{\mathbf{0}}\left(\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}-3\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)\cos 3\boldsymbol{\it\omega}t-\frac{1}{4}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{N}_{\mathbf{0}}\left(3\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}-\boldsymbol{N}_{\boldsymbol{ 0}}^{\mathbf{2}}\right)\sin 3\boldsymbol{\it\omega}t\end{eqnarray}$ (44)

To avoid secular terms, letting the coefficients for cos 𝝎t and sin 𝝎t in Eq. (44) are zero, it can be obtained $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}3{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{M}_{\mathbf{0}}^{\mathbf{3}}+3{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{M}_{\mathbf{0}}\boldsymbol{N}_{\mathbf{0}}^{\mathbf{2}}-4\boldsymbol{M}_{\mathbf{0}}\left(\boldsymbol{\it\omega}^{2}-\boldsymbol{\it\omega}_{\mathbf{0}}^{2}\right)-4{\varepsilon}\boldsymbol{F}_{\mathbf{1}}=\mathbf{0}\\[6.0pt] 3{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{M}_{\mathbf{0}}^{\mathbf{2}}\boldsymbol{N}_{\mathbf{0}}+3{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{N}_{\mathbf{0}}^{\mathbf{3}}-4\boldsymbol{N}_{\mathbf{0}}\left(\boldsymbol{\it\omega}^{2}-\boldsymbol{\it\omega}_{\mathbf{0}}^{2}\right)=\mathbf{0}\end{array}\right. & & \displaystyle\end{eqnarray}$ (45)

From Eq. (45), the amplitude-frequency relationship can be expressed by $\begin{eqnarray}\displaystyle \boldsymbol{\it\omega}^{2}=\boldsymbol{\it\omega}_{\mathbf{0}}^{2}+\frac{3}{4}{\varepsilon}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-\frac{1}{\boldsymbol{\it\delta}_{\mathbf{0}}}{\varepsilon}\boldsymbol{F}_{\mathbf{1}} & & \displaystyle\end{eqnarray}$ (46)

From initial condition and Eq. (45), it can be obtained M ₀ = 𝜹₀, N ₀ = 0. So the zero order approximate solution of nonlinear force vibration can be written as $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{0}}\left(\boldsymbol{t}\right)=\boldsymbol{\it\delta}_{\mathbf{0}}\cos \boldsymbol{\it\omega}t & & \displaystyle\end{eqnarray}$ (47)

Substituting Eq. (47) into Eq. (41), the first order approximate solution of nonlinear force vibration can be expressed by $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{1}}\left(\boldsymbol{t}\right) & = & \displaystyle \frac{2\boldsymbol{F}_{\mathbf{0}}-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{2\boldsymbol{\it\omega}^{2}}+\left(\frac{\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-3\boldsymbol{F}_{\boldsymbol{ 0}}+\boldsymbol{F}_{\mathbf{2}}}{3\boldsymbol{\it\omega}^{2}}-\frac{\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}^{2}}\right)\cos \boldsymbol{\it\omega}t\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -∼\frac{2\boldsymbol{F}_{\mathbf{2}}-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{6\boldsymbol{\it\omega}^{2}}\cos 2\boldsymbol{\it\omega}t+\frac{\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}^{2}}\cos 3\boldsymbol{\it\omega}t\end{eqnarray}$ (48)

By the same token, the second order approximate solution and the expression of 𝝈 can be expressed by $\begin{eqnarray}\displaystyle \boldsymbol{Q}_{\mathbf{2}}\left(\boldsymbol{t}\right) = \displaystyle \frac{1}{\boldsymbol{\it\omega}^{2}}\left(\boldsymbol{C}_{\mathbf{0}}+\boldsymbol{C}_{\mathbf{6}}\cos \boldsymbol{\it\omega}t-\frac{1}{3}\boldsymbol{C}_{\mathbf{2}}\cos 2\boldsymbol{\it\omega}t-\frac{1}{8}\boldsymbol{C}_{\mathbf{3}}\cos 3\boldsymbol{\it\omega}t-\frac{1}{15}\boldsymbol{C}_{\mathbf{4}}\cos 4\boldsymbol{\it\omega}t-\frac{1}{24}\boldsymbol{C}_{\mathbf{5}}\cos 5\boldsymbol{\it\omega}t\right)\text{}&&\end{eqnarray}$ (49) $\begin{eqnarray}\displaystyle \boldsymbol{\it\sigma} = \displaystyle \frac{\boldsymbol{\it\delta}_{\mathbf{0}}^{3}\left(224\boldsymbol{B}_{\boldsymbol{ 2}}^{\mathbf{2}}-144\boldsymbol{B}_{\boldsymbol{ 2}}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}+9\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\right)+48\boldsymbol{\it\delta}_{\boldsymbol{ 0}}\boldsymbol{F}_{\mathbf{0}}\left(9\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}-8\boldsymbol{B}_{\mathbf{2}}\right)-16\boldsymbol{\it\delta}_{\mathbf{0}}\boldsymbol{F}_{\mathbf{2}}\left(4\boldsymbol{B}_{\mathbf{2}}+9\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}\right)}{2\left(32\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}-96\boldsymbol{F}_{\mathbf{0}}+32\boldsymbol{F}_{\mathbf{2}}\right)}\text{}\quad&&\end{eqnarray}$ (50) where $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\boldsymbol{C}_{\mathbf{0}}=\displaystyle \frac{1}{96\boldsymbol{\it\omega}^{2}}\left[-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\left(48\boldsymbol{\it\sigma}-75\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}+32\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}+12\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}\right)+48\boldsymbol{F}_{\mathbf{0}}\left(2\boldsymbol{\it\sigma}+2\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}-3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\right)\right.\\[6.0pt] \quad \quad \,\,+\left.8\boldsymbol{\it\delta}_{\mathbf{0}}\boldsymbol{F}_{\mathbf{2}}\left(4\boldsymbol{B}_{\mathbf{2}}+3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}\right)\right]\\[6.0pt] \boldsymbol{C}_{\mathbf{2}}=\displaystyle \frac{1}{96\boldsymbol{\it\omega}^{2}}\left[\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\left(45\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-16\boldsymbol{\it\sigma}-32\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}\right)+48\boldsymbol{\it\delta}_{\mathbf{0}}\boldsymbol{F}_{\mathbf{0}}\left(2\boldsymbol{B}_{\mathbf{2}}-3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}\right)\right.\\[6.0pt]\quad\quad\,\,+\left.16\boldsymbol{F}_{\mathbf{2}}\left(2\boldsymbol{\it\sigma}-2\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}+3\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\right)\right]\\[9.60004pt] \boldsymbol{C}_{\mathbf{3}}=\displaystyle \frac{1}{384\boldsymbol{\it\omega}^{2}}\left[32\boldsymbol{\it\delta}_{\mathbf{0}}\boldsymbol{F}_{\mathbf{0}}\left(4\boldsymbol{B}_{\mathbf{2}}+9\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}\right)-96\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\boldsymbol{F}_{\mathbf{2}}\right.\\[6.0pt]\quad\quad\,\,\left.-\boldsymbol{\it\delta}_{\mathbf{0}}^{3}\left(64\boldsymbol{B}_{\mathbf{2}}^{\mathbf{2}}+96\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}+9\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-12\boldsymbol{\it\sigma}B_{\mathbf{3}}\right)\right]\\[9.60004pt] \boldsymbol{C}_{\mathbf{4}}=-\displaystyle \frac{1}{32\boldsymbol{\it\omega}^{2}}\left(5\boldsymbol{B}_{\mathbf{2}}\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{4}-8\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}\boldsymbol{F}_{\mathbf{2}}\right)\\[9.60004pt] \boldsymbol{C}_{\mathbf{5}}=\displaystyle \frac{3\boldsymbol{B}_{\mathbf{3}}^{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{5}}{128\boldsymbol{\it\omega}^{2}}\\[9.60004pt] \boldsymbol{C}_{\mathbf{6}}=\boldsymbol{C}_{\mathbf{0}}+\displaystyle \frac{1}{3}\boldsymbol{C}_{\mathbf{2}}+\frac{1}{8}\boldsymbol{C}_{\mathbf{3}}+\frac{1}{15}\boldsymbol{C}_{\mathbf{4}}+\frac{1}{24}\boldsymbol{C}_{\mathbf{5}}.\end{array}\right. & & \displaystyle\end{eqnarray}$ (51) Hence, from Eqs. (47)–(49), the approximate solution of the nonlinear forced vibration can be written as $\begin{eqnarray}\displaystyle \boldsymbol{Q}\left(\boldsymbol{t},\boldsymbol{\it\varepsilon}\right) & = & \displaystyle \boldsymbol{\it\delta}_{\mathbf{0}}\cos \boldsymbol{\it\omega}t+{\varepsilon}\left[\left(\frac{\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}-3\boldsymbol{F}_{\boldsymbol{ 0}}+\boldsymbol{F}_{\mathbf{2}}}{3\boldsymbol{\it\omega}^{2}}-\frac{\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}^{2}}\right)\cos \boldsymbol{\it\omega}t-\frac{2\boldsymbol{F}_{\mathbf{2}}-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{6\boldsymbol{\it\omega}^{2}}\cos 2\boldsymbol{\it\omega}t\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\left.\frac{\boldsymbol{B}_{\mathbf{3}}\boldsymbol{\it\delta}_{\mathbf{0}}^{3}}{32\boldsymbol{\it\omega}^{2}}\cos 3\boldsymbol{\it\omega}t+\frac{2\boldsymbol{F}_{\mathbf{0}}-\boldsymbol{B}_{\mathbf{2}}\boldsymbol{\it\delta}_{\mathbf{0}}^{2}}{2\boldsymbol{\it\omega}^{2}}\right]+\frac{{\varepsilon}^{2}}{\boldsymbol{\it\omega}^{2}}\left(\boldsymbol{C}_{\mathbf{0}}+\boldsymbol{C}_{\mathbf{6}}\cos \boldsymbol{\it\omega}t-\frac{1}{3}\boldsymbol{C}_{\mathbf{2}}\cos 2\boldsymbol{\it\omega}t\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -∼\left.\frac{1}{8}\boldsymbol{C}_{\mathbf{3}}\cos 3\boldsymbol{\it\omega}t-\frac{1}{15}\boldsymbol{C}_{\mathbf{4}}\cos 4\boldsymbol{\it\omega}t-\frac{1}{24}\boldsymbol{C}_{\mathbf{5}}\cos 5\boldsymbol{\it\omega}t\right)\end{eqnarray}$ (52)

5 Results and discussion

5.1 Analysis of nonlinear free vibration

The equations in this paper are utilized for the nonlinear vibration analysis of the piezoelectric precision drive system whose parameters are shown in Table 1. The linear and nonlinear natural frequencies are shown in Table 2, where Δω is difference between linear and nonlinear frequencies. The influence of system’s parameters on its natural frequencies are selected to be investigated for each parameter. The results obtained are shown in Figs 3–11.

In Table 2, the nonlinear frequencies are less than linear frequencies at ϵ = 0.2, the system exhibits characteristic of soft spring. As the increase of frequency order, the frequency difference Δω and chaning rate Δω∕ω₀ both decrease, hence, the nonlinear characteristic weakens. What’s more, the largest changing rate between linear and nonlinear frequencies occurs at the first order, the value of which is 11.82%. The reason for the higher changing rate is that in Linz Ted-Poincaré Method, the nonlinear frequencies are the function of small parameter ϵ. So when small parameter ϵ changes, the nonlinear frequencies change yet. In this paper, ϵ = 0.2 is selected to calculate the nonlinear frequencies, whereas 0.2 is a quite big value for the small parameter. Hence, the higher small parameter generates the larger frequency changing rate. Another reason for the large frequency change is that the experiment results [16] of piezoelectric material show that the piezoelectric material has strong nonlinear properties.

From Figs 3 and 4, with the increase of elastic stiffness coefficient of piezoelectric actuator c ₃₃, the nonlinear frequencies decrease, while Δω increases, the system’s nonlinear characteristic strengthens. However, as total length of piezoelectric actuator l _np grows, the nonlinear frequencies increase, Δω decreases. The reason for above phenomenon is that the nonlinear coefficients B ₂ and B ₃ increase with the grow of c ₃₃, while the ones decrease with the grow of l _np. So the changing of c ₃₃ strengthens the nonlinear characteristic of system, while l _np weakens the nonlinear characteristic. From another aspect, nonlinear forces increase with c ₃₃, while they decrease with l _np. These are the reasons for the nonlinear characteristic strengthens with c ₃₃ and weakens with l _np. So, in order to decrease the influence on nonlinear characteristic on dynamic performance of the system, parameters c ₃₃ and l _np should be designed reasonably.

From Figs 5–9, as length parameters l ₁ and l ₂ grow, nonlinear frequencies decrease first and then increase, but Δω increases first and then decreases. While as length l ₃ grow, nonlinear frequencies decrease marginally, and Δω increases slowly. Therefore, l ₃ has little influence on nonlinear characteristic. The law of nonlinear frequencies change with length l ₅ is similar with that of l ₁ and l ₂. However, as length l ₈ increases, nonlinear frequencies increase first and then decrease, while Δω increases slowly. So the influence of parameter l ₈ on nonlinear vibration is not obvious. The length parameters l ₁, l ₂, l ₃, l ₅ and l ₈ do not change the nonlinear forces directly, but they influence the exciting forces of each element. In addition, the effect of these length parameters on nonlinear frequencies mainly through changing the linear frequencies. On the other hand, exciting forces are influenced by nonlinear forces, and then influence the nonlinear frequencies indirectly.

From Figs 10 and 11, with the increase of sectional width of z-shaped rod l _B and sectional diameter of swaying rod d, nonlinear frequencies decrease, Δω increases. Both l _B and d have an important influence on nonlinear vibration of the system. What’s more, as the frequency order increases, the changing rate between nonlinear and nonlinear frequencies decrease. These show that the nonlinear characteristic weakens with the increase of frequency order. And the nonlinear effect has an important influence of lower order frequencies. Hence, during the piezoelectric precision drive system operates in the robot, lower order frequencies should be paid much attention, and parameters l _B and d should be chosen properly.

In a word, changes of parameters c ₃₃, l _B and d strengthens the nonlinear vibration of the system, while changes of l _np weakens the nonlinear vibration. Besides, parameters l ₃ and l ₈ have litter influence on nonlinear dynamics.

5.2 Analysis of nonlinear forced vibration

Eq. (46) is used for analyzing the amplitude-frequency relationship of piezoelectric actuator. Figure 12 shows the changes of amplitude--frequency curves of piezoelectric actuator along with parameters.

From Fig. 12, the nonlinear resonance of piezoelectric actuator occurs at far from natural frequencies. When exciting frequency takes a certain value, jump phenomenon of vibration amplitude occurs in nonlinear resonance. Besides, the amplitude-frequency curves deviate from its skeleton-line, and the deviation degree of which is in proportion to the value of F ₁. Hence, the amplitude-frequency characteristic of the piezoelectric actuator is mainly determined by the natural characteristic of the vibration system. And the jump of amplitude is a peculiar phenomenon in nonlinear dynamic system. In practical application, the special issue should be prohibited.

As small parameter ϵ increases, bending degree of the skeleton-line increases. For a certain exciting frequency, the maximum amplitude decreases. Therefore, ϵ has a great influence on both nonlinear frequencies and vibration amplitudes. However, as peak--peak voltage V _p−p grows, the skeleton-line remains unchanged, while the scope of amplitude--frequency curves deviate from its skeleton-line becomes larger. The reason is that the amplitude--frequency curves is affected by F ₁, while the skeleton-line is not. Hence, the increase of V _p−p changes the value of F ₁, and in turn affects the nonlinear resonance curves.

Piezoelectric strain constant d ₃₃ has litter influence on amplitude--frequency curves. With the increase of d ₃₃, resonance curves change slightly. However, as with V _p−p, the influence of d ₃₃ on amplitude--frequency is by altering F ₁. Unlike the other conditions, with the increase of elastic stiffness coefficient of piezoelectric actuator c ₃₃, the amplitude-- frequency curves and its skeleton-line move along the positive direction, while their curve shape does not change. So the influence of c ₃₃ on amplitude--frequency is by altering natural frequencies.

In a word, small parameter ϵ influences the trend and outline of the amplitude--frequency curves, peak--peak voltage V _p−p and piezoelectric strain constant d ₃₃ affect the outline of the amplitude--frequency curves, while elastic stiffness coefficient c ₃₃ influences the position of the amplitude--frequency curves. Hence, small parameter ϵ has a great influence of the amplitude--frequency characteristic. However, V _p−p, d ₃₃ and c ₃₃ only influence part of the amplitude--frequency characteristic.

Letting exciting frequency close to natural frequency 240914 rad/s. The nonlinear dynamic response near to resonance is presented by solving Eq. (52), and the results are shown in Fig. 13.

From Fig. 13, as the input signal of piezoelectric actuator is non-zero positive bias sinusoidal signal, so the equilibrium positions of the dynamic responses for each part are both above zero line. And the average amplitude of DB beam is largest. The nonlinear and linear vibration responses have the same phase, while the nonlinear amplitudes is slightly smaller than that of linear one.

The maximum error between nonlinear and linear amplitude of piezoelectric actuator is 10.4%, while the ones for OD beam, DB beam and swaying rod are 9.9%, 8.5% and 3.2%, respectively. So, the influence of nonlinear effect on amplitude is weakened orderly as piezoelectric actuator, OD beam, DB beam and swaying rod.

In short, nonlinear effect mainly affects the amplitudes of vibration response curves. The effect of nonlinear effect on piezoelectric actuator is largest. The reason is that the nonlinear effect comes from piezoelectric actuator directly. Whereas the nonlinear effect of OD beam, DB beam and swaying rod derives from nonlinear output forces of piezoelectric actuator.

In order to verify the correctness of the analytical solution, fourth order Runge–Kutta method is used to solve the numerical solution of nonlinear dynamic Eq. (20). Here, the comparision of the two type solution is presented in Fig. 14.

From Fig. 14, in the resonance response of piezoelectric actuator, the amplitude of numerical solution is 3.3% larger than analytical solution. While the phase difference increases over time. At 0.1 ms, the maximum phase difference is 0.17π. In the resonance response of OD beam, DB beam and swaying rod, numerical solutions are accordant with the analytical ones.

In general, the consistence between the numerical solutions and analytical solutions shows the effectiveness of the derivation process of nonlinear dynamics in this paper. Runge–Kutta method and Linz Ted-Poincaré Method are two kind of solving nonlinear dynamic equations. Their solving procedures are different. But their solving precisions are similar. In this paper, Linz Ted-Poincaré Method is introduced to solve the nonlinear equations of the piezoelectric drive system. Runge–Kutta method is just as an auxiliary validation method to verify the correctness of the theoretical derivation in this paper.

5.3 FEM simulation verification

To verify the correctness of the theoretical derivation in this paper, Finit Element Method (FEM) simulation of natural frequencies is conducted. Figure 15 gives FEM simulation figures of meshing model and the first mode cloud picture. Table 3 shows comparison of simulation and theoretical frequencies. Here, the simulation frequencies are compared with linear and nonlinear frequencies. Results show the simulation frequencies are close to the nonlinear frequencies. The errors between simulation frequencies and nonlinear frequencies are less than 1%. The FEM simulation proves that the nonlinear dynamic model is more close to reality conditions.

6 Conclusions

In this paper, the nonlinear coupled dynamic models and equations of the piezoelectric precision drive system used for robot joint are established. The approximate solutions of nonlinear free vibrations and forced vibrations are investigated. The results show:

The elastic stiffness coefficient c ₃₃, sectional width of z-shaped rod l _B and diameter of swaying rod d have significant influences on the nonlinear vibration of the drive system. The maximum error between nonlinear and linear amplitude of the system is up to 10.4%, thus, it shows that the nonlinear piezoelectric characteristic influences the dynamic performance of the system obviously. The consistence between the numerical solutions and analytical solutions shows the effectiveness of the derivation process of nonlinear dynamic solutions. The results obtained in this paper can be used to predict the dynamic load of the piezoelectric precision drive system used for robot joint.

Footnotes

Acknowledgements

This project is supported by National Natural Science Foundation of China (51605423), Research Project of State Key Laboratory of Mechanical System and Vibration (MSV201808) and Youth Science and Technology Innovation Project of Jiangsu University of Science and Technology.

Nomenclature

References

Pan

and Zhang

, Piezotronics and piezo-phototronics - From single nanodevices to array of devices and then to integrated functional system, Nano Today 8(6) (2013), 619–642.

Lin

and Li

, A Neural-repetitive control approach for high-performance motion control of piezo-actuated systems, Arab J Sci Eng 39(5) (2014), 4131–4140.

Wen

and Zi

, Design of a new piezoelectric energy harvester based on compound two-stage force amplification frame, IEEE Sensors Journal 18(10) (2018), 3989–4000.

David

O.U.

Johan

and Ralf

, Improved tactile resonance sensor for robotic assisted surgery, Mechanical Systems and Signal Processing 99 (2018), 600–610.

Ghosh

Jain

Majumder

, Experimental characterizations of bimorph piezoelectric actuator for robotic assembly, Journal of Intelligent Material Systems and Structures 28(15) (2017), 2095–2109.

Hadi

Hassan

and Gholamreza

, Study of a piezo-electric actuated vibratory micro-robot in stick-slip mode and investigating the design parameters, Nonlinear Dynamics 89(3) (2017), 1927–1948.

Tajitsu

, Piezoelectric Poly-L-lactic acid fabric and its application to control of humanoid robot, Ferroelectrics 515(1) (2017), 44–58.

Dharmawan

Hariri

Foong

Steerable miniature legged robot driven by a single piezoelectric bending unimorph actuator, in: Proceedings - IEEE International Conference on Robotics and Automation, Singapore, 2017, pp. 6008–6013.

Zou

Zhang

W.P.

X.J.

, The design and micro fabrication of a sub 100 mg insect-scale flapping-wing robot, Micro and Nano Letters 12(5) (2017), 297–300.

10.

Y.H.

, Study on the joint driver of micro robot based on piezoelectric element for celiac operation, Robot 25(4) (2003), 335–338.

11.

, Design and mechanical analysis of a mesh-type rotary harmonic piezoelectric motor, International Journal of Applied Electromagnetics and Mechanics 56(3) (2018), 461–478.

12.

and Xu

, Non-linear primary resonance of driving system for a harmonic piezoelectric motor, Journal of Vibration, Measurement and Diagnosis 36(3) (2016), 575–579.

13.

Xing

and Xu

, Coupled vibration of driving sections for an electromechanical integrated harmonic piezoelectric system, AIP Advances 4(3) (2014), 031320.

14.

and Xu

, Vibration characteristics of driving part for a novel harmonic piezodrive system under piezoelectric excitation, ICIC Express Letters 9(9) (2015), 2583–2589.

15.

and Li

, Coupled dynamics for an electromechanical integrated harmonic piezodrive system, Arabian Journal for Science and Engineering 39(12) (2014), 9137–9159.

16.

Tan

and Tong

, A one-dimensional model for non-linear behavior of piezoelectric composite materials, Composite Structures 58(4) (2002), 551–561.

17.

Liu

and Chen

, Nonlinear Vibrations, Beijing, Higher Education Press, 2001.