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

Authors: H. R. Liu; J. Z. Li; D. X. Hou; R. R. Yin; J. S. Jiang
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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.


## 2. Model of Piecewise-Nonlinear Asymmetric Hysteresis

Considering the single DOF vertical vibration system, as shown in Figure1, x is the oscillating amplitude of the elastic body in the vertical direction and x˙ is the vibration velocity. m is the mass of the elastic body. k is the restoring stiffness. c is the damping of the system. The external excitation of the system is Fcos⁡ωt.Figure 1
Asymmetric hysteresis dynamical system.The dynamical equation of the system is(1)mx¨+cx˙+kx+fx,x˙=Fcos⁡ωt,where(2)fx,x˙=a1x+b1x˙+a2x3+b2x˙3x˙≤0,a1′x+b1′x˙+a2′x3+b2′x˙3x˙>0.f
(
x
,
x
˙
) 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.


## 3. 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:(3)τ=ωt.Then (1) takes the dimensionless time form of(4)mω2x′′+cωx′+kx+fx,x′=Fcos⁡τ,where(5)x′′=d2xdτ2,x′=dxdτ,fx,x′=a1x+b1ωx′+a2x3+b2ω3x′3x′≤0,a1′x+b1′ωx′+a2′x3+b2′ω3x′3x′>0.For the convenience to analyze the system, transfer the form of the equation into(6)fx,x′=fx+fx′,where(7)fx=f1x+f3x,f1x=a1xx′≤0,a1′xx′>0,f3x=a2x3x′≤0,a2′x3x′>0,fx′=ωf1x′+ω3f3x′,f1x′=b1x′x′≤0,b1′x′x′>0,f3x′=b2x′3x′≤0,b2′x′3x′>0.Then the following incremental equation is derived:(8)mω2x′′+cωx′+kx+f1x+f3x+ωf1x′+ω3f3x′=Fcos⁡τ.The second step of the IHB method is the incremental procedure of Newton-Raphson. A substitution is then made forx and ω such that(9)xτ=x0τ+Δxτ(10)ω=ω0+Δω,where x0(τ) is related to the response of a nearby state of excitation at ω0 and Δx(τ) is the difference in the response between the two states when the excitation frequency is changed by Δω.f
1
(
x
), f3(x), f1(x′), and f3(x′) can be represented approximately by using their first-order Taylor series expansion, respectively:(11)f1x=f1x0+f1′x0Δx,f3x=f3x0+f3′x0Δx,f1x′=f1x0′+f1′x0′Δx′,f3x′=f3x0′+f3′x0′Δx′.By substituting (9), (10), and (11) into (8), and eliminating the high-order items, the following linearized differential equation on Δx is obtained:(12)mω02Δx′′+cω0Δx′+kΔx+f1′x0+f3′x0Δx+ω0f1′x0′+ω03f3′x0′Δx′=R+ΔωP,(13)R=Fcos⁡τ-mω02x0′′+cω0x0′+kx0+f1x0+f3x0+ω0f1x0′+ω03f3x0′,(14)P=-2mω0x0′′-cx0′-f1x0′+3ω02f3x0′,where R is the correction term. R=0 when x0 and ω0 are the accuracy solution. P is the unbalanced force generated by Δω.The next step of the IHB method is the harmonic balance process. Thereforex0 and Δx0 are expanded to be Nth order Fourier series: (15)x0=a02+∑n=1Nancos⁡nτ+bnsin⁡nτ=CsAΔx0=Δa02+∑n=1NΔancos⁡nτ+Δbnsin⁡nτ=CsΔA,where(16)A=a0,a1,…,aN,b1,b2,…,bNTΔA=Δa0,Δa1,…,ΔaN,Δb1,Δb2,…,ΔbNTCs=12,cos⁡τ,…,cos⁡Nτ,sin⁡τ,…,sin⁡Nτ.The solution for system response is represented in form such that(17)x0=CsA,Δx0=CsΔA.The Galerkin procedure is then applied to handle (12) after substituting (17) into it: (18)∫02πmω02Δx′′+cω0Δx′+kΔx+f1′x0+f3′x0Δx+ω0f1′x0′+ω03f3′x0′Δx′δΔxdτ=∫02πR+ΔωPδΔxdτwhich is equivalent to the linearized equations with ΔA and Δω being variables: (19)TKCΔA=R+ΔωP,where(20)TKC=TKC11TKC12TKC21TKC22,R=RKC1RKC2,P=P1P2.The expressions ofTKC,RKCi,Pi(i=1,2;j=1,2) are in the Appendix.


## 4. 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 x(τ) can be given by the expression (9). Thirdly, the updated parameter TKC 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 thatm=1, k=1, c=0.01, F=2, ω=1, a1=0.2, a2=3, b1=0.02, b2=0.3, a1′=0.1, a2′=2.5, b1′=0.01, b2′=0.25, setting the number of terms of the Fourier series as N=3; even-order 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(21)xt=p1cos⁡ωt+q1sin⁡ωt+p2cos⁡2ωt+q2sin⁡2ωt+p3cos⁡3ωt+q3sin⁡3ωt.The lagging forcef(x,x˙) 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.Figure 2
Asymmetric hysteresis force.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.Figure 3
Amplitude frequency response characteristics.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 itemsp1,p2,p3,q1,q2,q3 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.Figure 4
Relationship betweenp1 and ω.Figure 5
Relationship betweenp2 and ω.Figure 6
Relationship betweenp3 and ω.Figure 7
Relationship betweenq1 and ω.Figure 8
Relationship betweenq2 and ω.Figure 9
Relationship betweenq3 and ω.Decreasing the damping, the system became a weak nonlinearity system. Specifyingb2=0.03,b2′=0.025, the bistable regions A and B of the bifurcation curve were analyzed, as shown in Figure 10.Figure 10
Bistable region of bifurcation.And when the damping changed asb2=0.003,b2′=0.0025, the characteristic curves were plotted in Figure 12. The bistable regions A, B, and C were shown.It can be seen from Figures10–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.Figure 11
Local amplification of bistable region.Figure 12
Bistable region of bifurcation.Figure 13
Local amplification of bistable region.


## 5. 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.



---

*Source: 296085-2015-10-26.xml*