Analysis of Dynamical Characteristic of Piecewise-Nonlinear Asymmetric Hysteretic System Based on Incremental Harmonic Balance Method

Considering a sort of elastic mass with asymmetric hysteresis characteristic which widespread existing in engineering field, a piecewise-nonlinear dynamical equation, which contains asymmetric hysteretic loop, is established. By using the method of Incremental Harmonic Balance (IHB), the analytic linearized algebraic equation of the system is obtained. On the basis of this algebraic equation, the coefficients of the algebraic expression are figured out by the incremental procedure and the iterative process of the regulated variable. Through the emulation, the amplitude frequency response curve and the relation between the value of the harmonic component and the external excitation are researched; the bistable regions of the bifurcation diagram of the system under the variation of the excitation amplitude are studied. The above results can be used to guide the research of the asymmetric hysteretic system with polynomial expression.


Introduction
Hysteresis phenomena exist widely in nature, such as magnetic hysteresis, internal damping of material, and damping generated by dry friction [1].Elastic body will become inelastic and presents hysteresis characteristics for the plastic deformation under overloading.The restoring force and displacement of elastic body form hysteresis loop [2].
Generally the hysteresis restoring force is symmetric during the process of load-on and load-off, and the corresponding hysteresis models are established according to the symmetry, for instance, bilinear hysteretic nonlinear system, dry friction ideal model [3], Davidenkov model [4], Bouc-Wen equation [5], and polynomial model [6].These symmetric systems have been researched deeply.
Under certain condition, the deformation feature of the elastic body change because of its different vibration direction and the hysteresis restoring force become asymmetric.Wang researched the asymmetric hysteresis dynamic characteristic of symmetric magnetohydrodynamic damper caused by the asymmetry of system parameters under different direction of motion [7,8].But the researches on the asymmetric hysteresis model of polynomial expression have not been made profoundly.
Usually, the average method is used to solve the piecewise-nonlinear dynamical equation.However, the global characteristics of piecewise-nonlinear system are always analyzed by using the IHB method whose accuracy is better than the average method.The IHB method was extracted by Lau and Cheung [9]; it has been widely applied to analyze several kinds of nonlinear equation.Lau and Zhang solved the piecewise-nonlinear equation with linear rigidity and got the amplitude frequency curves and the relationship between the harmonic constant and the external excitation frequency [10].Xu [11,12] analyzed the piecewise-nonlinear equation with linear damping by applying the IHB method and compared the simulation results between the IHB method and the Runge-Kutta method.The multiple DOF equation with Coulomb damping and quadratic damping was established by Zhou and Zhang [13] and solved by using the IHB method.
And the simulation under different number of harmonic items certificated that the more the number of harmonic items, the more the accuracy of analytic solution.
Hou et al. [14] analyzed the modeling and simulation of dynamically rolling force under the asymmetric hysteretic deformation of rolling mill; the asymmetric hysteretic rolling force was facilitated to multinomial asymmetric hysteretic system.
Considering a type of asymmetric hysteresis model of polynomial expression with cubic stiffness and cubic damping, the piecewise-nonlinear dynamical system was established in this paper.The analytic linear algebraic equation was obtained after solving the dynamical equation of the system by applying the IHB method.After the alternative incremental process and repeated iterative process of the equation, each coefficient of the linear algebraic equation can be computed.Then the amplitude frequency response characteristics and the relationship between the parameters of the harmonic items and the external frequency under different vibration direction were researched.The bistable regions of the bifurcation diagram of the system under the variation of the excitation amplitude were studied.These works can be used to guide the research of the multinomial asymmetric hysteretic system existing widely in engineering field.

Model of Piecewise-Nonlinear Asymmetric Hysteresis
Considering the single DOF vertical vibration system, as shown in Figure 1,  is the oscillating amplitude of the elastic body in the vertical direction and ẋ is the vibration velocity. is the mass of the elastic body. is the restoring stiffness. is the damping of the system.The external excitation of the system is  cos().
The dynamical equation of the system is where (, ẋ ) represents the asymmetric hysteresis force of the vertical vibration equation.The different vibration directions of the elastic body result in the different coefficients of the nonlinear model, and then the hysteresis force is generated.So (1), which contains asymmetric hysteresis model, is a vertical vibration dynamical equation with piecewise-nonlinear characteristic.

Incremental Harmonic Balance Method
The first step in the IHB method is to change the variables to transform (1) from the time domain to a dimensionless domain with the substitution: Then (1) takes the dimensionless time form of where For the convenience to analyze the system, transfer the form of the equation into where Then the following incremental equation is derived: The second step of the IHB method is the incremental procedure of Newton-Raphson.A substitution is then made for  and  such that where  0 () is related to the response of a nearby state of excitation at  0 and Δ() is the difference in the response between the two states when the excitation frequency is changed by Δ. 1 (),  3 (),  1 (  ), and  3 (  ) can be represented approximately by using their first-order Taylor series expansion, respectively: By substituting ( 9), (10), and ( 11) into (8), and eliminating the high-order items, the following linearized differential equation on Δ is obtained: where  is the correction term. = 0 when  0 and  0 are the accuracy solution. is the unbalanced force generated by Δ.The next step of the IHB method is the harmonic balance process.Therefore  0 and Δ 0 are expanded to be th order Fourier series: where The solution for system response is represented in form such that The Galerkin procedure is then applied to handle (12) after substituting (17) into it: which is equivalent to the linearized equations with ΔA and Δ being variables: where The expressions of T KC , R KC , P  ( = 1, 2;  = 1, 2) are in the Appendix.

Dynamical Analysis and Simulation
The IHB is a type of nonlinear derivation method, which contains semianalytic process and seminumerical process.The analytic expression of the system is obtained by using the analysis, but the coefficients of the expression are determined by the numerical method [15].
Expression ( 19) is a linear matrix equation, and the unknown quantity is one more than the number of the equations.So one of the unknown quantities must be selected as the incremental quantity for manipulation.The unknown quantity Δ was taken as the manipulation incrementation.First, the initial value of Δ was preset.Second, the other unknown quantity ΔA was obtained by solving (19); then () can be given by the expression (9).Thirdly, the updated parameter T KC and R will be refreshed, and then the other unknown quantity ΔA will be refreshed by ( 19).This procedure was the iterative process.After carrying out this procedure uninterruptedly until the error vector R is less than the error constant set up previously, the amplitude of vibration A corresponding to  was obtained, which can be used to get the solution of the system.When this iteration procedure is finished the manipulation incrementation will be added for Δ, and the iteration procedure will be carried out again for the vibration amplitude A corresponding to the refreshed parameter .The solution of the nonlinear vibration equation ( 1) is the application of the incremental process and iterative process alternatively and repeatedly.The parameters of the dynamical system are that  = 1,  = 1,  = 0.01,  = 2,  = 1,  1 = 0.2,  2 = 3,  1 = 0.02,  2 = 0.3,   1 = 0.1,   2 = 2.5,   1 = 0.01,   2 = 0.25, setting the number of terms of the Fourier series as  = 3; evenorder harmonics are not included in the Fourier series from the nature of the nonlinear terms of the equation.Therefore the solution of the system is The lagging force (, ẋ ) of the polynomial asymmetric hysteretic model of system (4), as shown in Figure 2, can be obtained from expression (4).The lagging force appears asymmetric because of the parameters of the polynomial hysteresis loop under different vibration direction.
According to the incrementation process and iterative process of the manipulator Δ based on the analytic linear matrix equation, the relation curves between the amplitude of vibration and the external excitation of the system were plotted in Figure 3.The amplitude of vibration and the jumping phenomena of the system under different vibration direction exhibit different disciplinarians, and hysteretic phenomena exist between the two kinds of amplitude frequency characteristics curves.When the excitation frequency varies around the natural frequency of the system, the amplitude of vibration will enlarge and the hysteretic phenomena between the amplitude frequency characteristics curves will become apparent.
The coefficients of the harmonic items, whose variation determines the property of the system and relationship between each harmonic component and the system frequency, are significant in solving the equation by using the IHB method.The relationship curves between the coefficients of each of the harmonic items  1 ,  2 ,  3 ,  1 ,  2 ,  3 and the excitation frequency were shown in Figures 4-9.It can be seen that the low order harmonic item operates master effect and determines the shape of the amplitude frequency response characteristics.The high-order harmonic items have impact on the partial feature of the amplitude frequency response characteristics.The visualized regularity of the system can be studied by researching the relationship curves between each harmonic component and the system frequency.
Decreasing the damping, the system became a weak nonlinearity system.Specifying  2 = 0.03,   2 = 0.025, the bistable regions A and B of the bifurcation curve were analyzed, as shown in Figure 10.
And when the damping changed as  2 = 0.003,   2 = 0.0025, the characteristic curves were plotted in Figure 12.The bistable regions A, B, and C were shown.
It can be seen from Figures 10-13 that the pattern of the bifurcation characteristic of the system under the variation   of the external excitation in the different vibration direction exhibits hysteresis phenomena, and the decreasing of the damping will deepen these phenomena.

Conclusions
Asymmetric hysteresis characteristic is a sort of peculiar piecewise-nonlinear characteristic which is distinguished from traditional symmetric hysteresis feature.This kind of characteristic exists widely in engineering.In this paper, an elastic body having asymmetric hysteresis characteristic was taken into account, employing a type of asymmetric hysteresis polynomial model with cubic stiffness and cubic damping, and the piecewise-nonlinear dynamical system containing the asymmetric hysteretic loop was established.Based on the IHB method, the analytic linear algebraic equation was obtained by solving the dynamical equation of the system.After the incremental process and iterative process of the equation, each coefficient of the linear algebraic equation can be calculated.The amplitude frequency response characteristics and relation curves between the coefficients of the harmonic items and the frequency of the external excitation during the different vibration directions were simulated.The bistable regions of the bifurcation diagram of the system under the variation of the excitation amplitude were studied.The work enriched the research of the asymmetric hysteresis model, and the above results can be used to guide the research of the multinomial asymmetric hysteretic system existing widely in engineering field.