Global bifurcation and circuit simulation of a soft-impacting system

1. Introduction

The dynamic study of piecewise-smooth mechanical systems aims to elucidate the complex behaviors induced by non-smooth factors such as clearance, constraints, and friction, and related research has attracted significant attention [1-6]. For mechanical systems, due to factors such as lubricant film or thermal expansion and contraction, clearance must be intentionally provided in the design of certain equipment. Some equipment may experience clearance due to wear and tear caused by service or errors during manufacturing and assembly. Due to the existence of clearance, there will be repeated impact phenomena of contact, detachment, re-contact and re-detachment between components, which will adversely affect the operation efficiency. On the other hand, some mechanical equipment designed based on the principle of vibro-impact has also been developed rapidly, such as vibration conveyor [7], impact hammer [8], vibration molding machine [9, 10], pendulum [11], vibrating screen [12] and so on. Therefore, for vibro-impact systems with clearance, how to seek benefits, avoid harm, improve reliability, and reduce noise has both theoretical value and important engineering practical significance. The dynamics of the vibro-impact system with clearance are mainly revealed by theoretical analysis and numerical simulation. Yue [13] deduced the symmetry property of the Poincaré map and determined the normal form of the Neimark–Sacker-pitchfork bifurcation, and finally obtained the Neimark-Sacker-pitchfork bifurcation by numerical simulation. For a gear system with multiple parameters, Gou [14] studied its dynamic response through the use of improved cell mapping. Zhang [15] deeply studied the coexisting attractors based on the Poincaré map and Lyapunov exponents and first demonstrated that high-period attractors and chaotic attractors can coexist for some parameter values. Shi [16] used chaotic characterization tools, such as Poincaré map, bifurcation diagram, maximum Lyapunov exponent diagram and phase diagram to reveal the dynamics of gear system. By combining the GPU parallel computing technique with the Lyapunov exponents and Poincaré section, bifurcation phenomenon of a mechanical speed control system has been conducted in-depth research by Deng [17]. Tan [18] used numerical methods such as GPU and numerical continuation to analyze the dynamic response of an impact oscillator with elastic constraints. Sun [19] proposed a Hopf bifurcation stability criterion applicable to maglev systems, and investigated the influence law of guideway parameters on vehicle-guideway coupled vibration through theoretical analysis and numerical simulation. Luo's team proposed a two-parameter bifurcation analysis method based on two-parameter bifurcation concept to reveal the vibro-impact system with clearance, including elastic constraint system [20-22], rigid constraint system [23, 24], plastic impact system [25] and so on.

Numerical simulation is based on the classical fourth-order Runge–Kutta method for solving ordinary differential equations. It uses programming software such as MATLAB or Visual C++ 6.0 to solve motion differential equations and generate dynamic results that characterize chaotic characteristics, such as bifurcation diagrams, phase plane diagrams, time history diagrams, etc. This method not only requires high programming ability, but also high-precision numerical simulation is limited by the operation performance of the computer, and the operation simulation speed is slow, which belongs to static calculation. As the system parameters are changed, the program needs to be recalculated, which takes up a large amount of computer capacity, consumes a lot of time and has poor economy. In fact, equivalent circuit design based on system mathematical models not only enables the analysis of nonlinear vibration characteristics, but also has advantages such as fast speed, dynamic adjustment of operation process, and convenient modal conversion. Ref [26] introduced a new four-scroll chaotic system consisting of two autonomous Duffing oscillators with memristive nonlinearity by using electronic circuit. The results show that the dynamic responses obtained by the circuit and numerical simulation are completely consistent. However, most of the studies focus on Chua's circuit [27-29] and Duffing oscillator [30, 31], and the simulation design of equivalent electronic circuit for non-smooth vibro-impact system is rare. An insight is obtained to the global dynamics for a new electro-vibro-impact system in Ref [32]. Ref [33] analyzed a piecewise-linear oscillator that consists of a three-resonant circuit with a hysteresis element.

In this study, we considered a soft-impacting system. In Multisim 14.0 environment, the integral circuit, addition circuit and inverse circuit composed of integrated operational amplifier, resistor and capacitor are used to design the circuit schematic diagram which is completely equivalent to the mechanical system. The nonlinear module is realized through the saturation characteristics of the operational amplifier. The output results of the virtual oscilloscope are in good agreement with the results of numerical simulation. The circuit simulation not only has fast operation speed, but also can realize the dynamic adjustment of parameters, which gives a way to analyze vibration, impact, and noise in nonlinear research and engineering fields.

2. Mechanical model

A two-degree-of-freedom vibration system with clearance composed of symmetric elastic stops is depicted in Fig.1. Although the model simplifies by excluding complex factors such as oil film forces and nonlinear stiffness, it holds significant qualitative value: the revealed dynamic patterns provide a theoretical foundation for understanding real-world systems. For instance, even in the presence of oil film forces, stiffness transitions and impacts induced by clearances remain key mechanisms for vibration amplification. Table 1 presents the definitions and units of the variables. The two oscillators with mass $M_1$ and $M_2$ are connected by a pair of linear spring-damping elements ($K_1$, ${\bar{C}}_1$), and another pair of elements ($K_2$, ${\bar{C}}_2$) couple the oscillator $M_2$ and the foundation. The harmonic excitation force $P_i\sin (\mathrm{\Omega }T+\tau )$ acts on oscillator $M_i$, where $P_i$, $\mathrm{\Omega }$ and $\tau$ are the amplitude, excitation frequency and phase angle, respectively ($i=\mathrm{1,2}$). The coordinate system is established with the static equilibrium position of the system as the origin, and the displacements of the two oscillators are $X_1$ and $X_2$, respectively. There are elastic constraints $A_1$ and $A_2$ with stiffness $K_0$ at the distance $B$ between the two sides of the oscillator $M_1$. As the amplitude of the harmonic excitation is small, the system vibrates linearly, and the oscillator $M_1$ does not collide with the elastic constraints $A_1$ or $A_2$. As the displacement of the oscillator $M_1$ is greater than the designed clearance threshold $B$, i.e., $\left|X_1\right|>B$, the impact event occurs in the system. Elastic constraints cause the systems with linear characteristics to exhibit complex non-smooth features.

The motion differential equation of the piecewise linear system in Fig. 1 is written as:

The piecewise linear spring force $F\left(X_1\right)$ generated by the collision is given by:

In the present model, the contact interface is simplified as a constraint spring with stiffness $K_0$. Future work could draw inspiration from the approach to bolted joint interfaces described in Ref. [34], describing the local nonlinear stiffness (e.g., Hertzian contact stiffness) in the clearance contact region as an “equivalent virtual material layer” that varies with contact state. This would facilitate the development of dynamic models that more accurately reflect the micro-mechanical properties of the interface.

Fig. 1Two-degree-of-freedom periodic forced vibration system with elastic constraints

The dimensionless parameters, time and variables are technically selected as:

then, the dimensionless expression of Eq. (1) is rewritten as:

in which, the expression of the spring force $f\left(x_1\right)$ is:

3. Poincaré maps

For such harmonically excited systems featuring spring-damper cushioning mechanisms, the term “soft impact” is commonly used to characterize the collision behavior. Unlike rigid impact, rigid impact uses a discrete mapping of recovery coefficients to describe sudden changes in motion state, while soft impact considers the time effect of contact deformation during collision and typically uses piecewise functions to describe the non-smooth characteristics of the system. By varying the stiffness $K_0$ from zero to infinity within a specified interval, the system can transition from a piecewise linear model to a rigid impact model. The dynamic behaviors under these two limiting conditions have been thoroughly investigated in Refs. [35]. Various mathematical models corresponding to soft impact have been developed for different practical scenarios, such as the Kelvin-Voigt model in Refs. [36], the Hertz contact model in Refs. [37], and the piecewise linear model in Ref. [38].

Due to the existence of constraints, the system motion is separated into several subspaces, i.e., ${\mathrm{\Sigma }}_1$, ${\mathrm{\Sigma }}_2$ and ${\mathrm{\Sigma }}_3$, in which, ${\mathrm{\Sigma }}_1={\left\{\left(x_1,{\dot{x}}_1,x_2,{\dot{x}}_2,\theta \right)\right.}^T\in {\mathit{R}}^4\times \mathit{S}\left.\left|H_1\left(X\right)>0\right.\right\}$ indicates that the oscillator $M_1$ collides with the constraint $A_1$; ${\mathrm{\Sigma }}_2={\left\{\left(x_1,{\dot{x}}_1,x_2,{\dot{x}}_2,\theta \right)\right.}^T\in {\mathit{R}}^4\times \mathit{S}\left.\left|H_1\left(X\right)<0\cap \right.H_2\left(X\right)>0\right\}$ indicates that there is no collision between the oscillator $M_1$ and the elastic constraints ($A_1$and $A_2$) on both sides; ${\mathrm{\Sigma }}_3={\left\{\left(x_1,{\dot{x}}_1,x_2,{\dot{x}}_2,\theta \right)\right.}^T\in {\mathit{R}}^4\times \mathit{S}\left.\left|H_2\left(X\right)<\text{0}\right.\right\}$ indicates that the oscillator $M_1$ collides with the lower elastic constraint $A_2$, where, $H_1\left(X\right)=X_1-B$, $H_2\left(X\right)=X_2+B$.

The general solution of Eq. (4) is:

where, symbols ${\mathrm{\Phi }}_{ij}$, ${\breve{\mathrm{\Phi }}}_{ij}$, ${\widehat{\mathrm{\Phi }}}_{ij}$ are matrix elements corresponding to state paces ${\mathrm{\Sigma }}_1$, ${\mathrm{\Sigma }}_2$, ${\mathrm{\Sigma }}_3$, respectively. ${\eta }_j=\zeta {{\omega }_{nj}}^2$, ${\breve{\eta }}_j=\zeta {{\breve{\omega }}_{nj}}^2$, ${\widehat{\eta }}_j=\zeta {{\widehat{\omega }}_{nj}}^2$, ${\omega }_{dj}=\sqrt{{{\omega }_j}^2-{{\eta }_j}^2}$, ${\breve{\omega }}_{dj}=\sqrt{{{\breve{\omega }}_j}^2-{{\breve{\eta }}_j}^2}$, ${\widehat{\omega }}_{dj}=\sqrt{{{\widehat{\omega }}_j}^2-{{\widehat{\eta }}_j}^2}$. Symbols ${\omega }_j$, ${\breve{\omega }}_j$ and ${\widehat{\omega }}_j$ represent the undamped natural frequencies in three subspaces (${\mathrm{\Sigma }}_1$, ${\mathrm{\Sigma }}_2$, ${\mathrm{\Sigma }}_3$), respectively. Symbols $A_j$, ${\breve{A}}_j$, ${\widehat{A}}_j$ and $B_j$, ${\breve{B}}_j$, ${\widehat{B}}_j$ are amplitude constants in three subspaces (${\mathrm{\Sigma }}_1$, ${\mathrm{\Sigma }}_2$, ${\mathrm{\Sigma }}_3$), respectively. The integral constants $a_j$, ${\breve{a}}_j$, ${\widehat{a}}_j$ and $b_j$, ${\breve{b}}_j$, ${\widehat{b}}_j$ are determined by the initial conditions of the vibration system shown in Fig. 1. $x_{0i}={\left(x_{01},x_{02}\right)}^T={\mathrm{\Phi }}^T{\left(\frac{k_0\left(x_1-b\right)}{{\omega }_{n1}^2,0}\right)}^T$ and ${\widehat{x}}_{0i}={\left({\widehat{x}}_{01},{\widehat{x}}_{02}\right)}^T={\widehat{\mathrm{\Phi }}}^T{\left(\frac{-k_0\left(x_1-b\right)}{{\widehat{\omega }}_{n1}^2,0}\right)}^T$.

We introduce the symbol $n-p-q$ to describe the basic characteristics of the periodic motion, where $n=\mathrm{1,2},3,...$ denotes the number of periods of the excitation force, $p$ and $q$ denote the number of impacts associated with the oscillator $M_1$ and the elastic constraints $A_1$and $A_2$ on both sides, respectively ($p,q=\mathrm{1,2},3,...$). In particular, for $1-p-p$ response with $n=1$ as well as $p=q$, we call it symmetric single-period multi-impact motion.

Based on the meaning of the symbol $n-p-q$, three Poincaré maps are defined:

Fixed points corresponding to section ${\sigma }_n$ can determine the period number $n$. Fixed points corresponding to ${\sigma }_{A_1}$ and ${\sigma }_{A_2}$ sections can respectively determine the number of collisions $q$ and $p$ of the oscillator $M_1$ at the upper ($A_1$) / lower ($A_2$) limiting constraints on both sides. The impact Poincaré map is briefly represented as:

where, $X^{\left(i\right)}=({\dot{x}}_{1+}^{\left(i\right)},x_2^{\left(i\right)},{\dot{x}}_{2+}^{\left(i\right)},\hspace{0.33em}{\tau }^{\left(i\right)}{)}^T$, $X^{(i+1)}=({\dot{x}}_{1+}^{(i+1)},x_2^{(i+1)},{\dot{x}}_{2+}^{(i+1)},{\tau }^{(i+1)}{)}^T$, $X\in R^4$. The symbol $\mu$ is the system parameter, which can be any of the system parameters, i.e., $m$, $k$, $k_0$, $c$, $\zeta$, $f_{20}$, $\omega$ and $b$.

The response $n$-1-1 with $n$ -period and double-impact is likely to appear in the vibro-impact system with symmetric elastic constraints in Fig. 1. The characteristic of the response $n$-1-1 lies in the fact that the instantaneous time after the collision between the oscillator $M_1$ and the elastic stop $A_1$ is zero, then the instantaneous time before the next impact is exactly $t=2n\pi /\omega$. Periodicity conditions for symmetric $n$-1-1 response are given as follows:

4. 1-p-p periodic motion and chattering-impact characteristic

In this section, the variable-step fourth-order Runge-Kutta integration method is employed to numerically simulate the vibro-impact system with clearance shown in Fig. 1. To eliminate the influence of transient processes and ensure the reliability of the analysis, the transient responses from the first 500 periods were discarded, and the steady-state response data from the subsequent 30 periods were retained for analyzing the dynamics of the system. The 1-$p$-$p$ response and chattering-impact characteristics are analyzed. The step size is $\mathrm{\Delta }t=$ 0.005 and the system parameters are fixed as $b=$ 0.25, $f_{20}=$ 0.0, $c=$ 0.5, $k=$ 0.5, $k_0=$ 0.95, $m=$ 0.5, $\zeta =$ 0.1. Figs. 2(a)-2(a2) intuitively show the dynamic response of the system with frequency in the range of $\omega \in \left[\mathrm{0,1.8}\right]$ by means of bifurcation diagrams ${\dot{x}}_1\left(\omega \right)$ and $x_{1min}$ as well as Lyapunov exponent ${\lambda }_{max}$. Specifically, Figs. 2(a)-2(c) are bifurcation diagrams $x_{1min}$ based on Poincaré mapping ${\sigma }_n$, where the number of branches in the curve $x_{1min}$ determines the number of periods $n$; Figs. 2(a1)-2(c1) are bifurcation diagrams based on Poincaré mapping ${\sigma }_{A1}$, where the number of branches in the curve ${\dot{x}}_1\left(\omega \right)$ determines the number of impacts $q$ between oscillator $M_1$ and upper constraint $A_1$; Figs. 2(a2)-2(c2) are the largest Lyapunov exponent corresponding to Figs. 2(a1)-2(c1). According to whether the Lyapunov exponent is greater than zero, the stability of the system can be effectively determined. Figs. 2(b)-2(c) are the detailed descriptions of local intervals in Fig. 2(a). Figs. 2(b1)-2(c1) are the detailed descriptions of Fig. 2(a1). Figs. 2(b2)-2(c2) are the detailed descriptions of Fig. 2(a2).

Symmetrical 1-1-1 motion is a typical dynamic response of the system with symmetric constraints, which usually occurs Pitchfork bifurcation. As shown in Fig. 2(c1), Pitchfork bifurcation (abbreviated as PF Bif in Fig. 2(c1)) causes the symmetric 1-1-1_S response to transit to asymmetric response 1-1-1_AS with the decrease of frequency $\omega$ (the subscript $S$ is the abbreviation for “symmetric”, and $AS$ is the abbreviation for “asymmetric”). At the Pitchfork bifurcation point ${\omega }_1=$ 1.47, the largest Lyapunov exponent is ${\lambda }_{{\omega }_1}=0$, as shown in Fig. 2(c2). Fig. 3(a) shows the phase plane portrait $(x_1,{\dot{x}}_1)$ of the symmetric 1-1-1_S response at $\omega =\text{1.55}$ before the Pitchfork bifurcation point, and Fig. 3(b) shows the phase plane portrait $(x_1,{\dot{x}}_1)$ of the asymmetric 1-1-1_AS response at $\omega =$ 1.43 after Pitchfork bifurcation point. The asymmetric response 1-1-1_AS presents two different forms of motion trajectories due to different initial values. The two are antisymmetric, such as the black and red motion trajectories in Fig. 3(b). As $\omega$ continues to decrease, the 1-1-1_AS response transits to 1-1-2 through the Grazing bifurcation (abbreviated as G Bif in Fig. 2(c1)), and then the system enters chaos through the period-doubling cascade of the 1-1-2 response.

As $\omega$ continues to decrease, 1-2-2 response degenerates from chaotic response and again embeds chaos through the period-doubling cascade. Repeatedly, until the number of impacts $p$ increases to a large enough number, the system exhibits chattering-impact characteristics, which produces noise pollution, causing loose parts and fatigue, and even lead to machine parts fracture. This paper conducts an in-depth study on the low-frequency dynamics of the system with clearance in Fig. 1, laying a theoretical foundation for high-precision condition monitoring. It is noteworthy that the novel intelligent sensing structure provides a new approach to directly obtain high signal-to-noise ratio vibration responses [39]. By sensing the dynamic load of the bearing body, it can effectively capture rich dynamic behaviors such as impacts and bifurcations induced by clearance nonlinearity.

For small clearance, as frequency $\omega$ decreases, the transition law from symmetric $1-1-1_S$ response to chattering response is $\omega$: 1-$\tilde{p}$-$\tilde{p}$ ←…← 1-$p$-$p$ ←$$chaos ←…← PD Bif ← 4-8-8 ← PD Bif ← 2-4-4 ← PD Bif ← 1-2-2 ← chaos ←… PD Bif ← 4-4-8 ← PD Bif ← 2-2-4 ← PD Bif ← 1-1-2 $\leftarrow$G Bif ← 1-1-1_AS ← PF Bif ← 1-1-1_Swhere, PF Bif indicates Pitchfork Bifurcation; G Bif indicates Grazing Bifurcation; PD Bif indicates Period-doubling Bifurcation.

Fig. 2Bifurcation diagrams x1min and x˙1(ω), and Lyapunov exponent diagram λmax, b= 0.25

a) Bifurcation diagram ${\mathrm{x}}_{1\mathrm{m}\mathrm{i}\mathrm{n}}$

b) Enlarged view 1 of Fig. 2(a)

c) Enlarged view 2 of Fig. 2(a)

a1) Bifurcation diagram ${\dot{\mathrm{x}}}_1\left(\mathrm{\omega }\right)$

b1) Enlarged view 1 of Fig. 2(a1)

c1) Enlarged view 2 of Fig. 2(a1)

a2) Diagram ${\mathrm{\lambda }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$

b2) Enlarged view 1 of Fig. 2(a2)

c2) Enlarged view 2 of Fig. 2(a2)

Fig. 3Phase plane portraits (x1,x˙1) of 1-p-p responses (p= 1-8), b= 0.25

a) 1-1-1_S, $\omega =$ 1.55

b) 1-1-1_AS, $\omega =$ 1.43

c) 1-1-2, $\omega =$ 1.34

d) 1-2-2_S, $\omega =$ 0.85

e) 1-2-2_AS, $\omega =$ 0.75

f) 1-3-3_S, $\omega =$ 0.54

g) 1-3-3_AS, $\omega =$ 0.515

h) 1-4-4_S, $\omega =$ 0.41

i) 1-4-4_AS, $\omega =$ 0.4

j) 1-5-5_S, $\omega =$ 0.335

k) 1-5-5_AS, $\omega =$ 0.32

l) 1-6-6_S, $\omega =$ 0.285

m) 1-6-6_AS, $\omega =$ 0.27

n) 1-7-7_S, $\omega =$ 0.24

o) 1-8-8_S, $\omega =$ 0.22

Figs. 3(a)-3(t) shows the phase plane portraits of 1-$p$-$p$ ($p=$1-8) responses at different excitation frequencies. For 1-$p$-$p$ response characterized by $p\ge 2$, there are also two antisymmetric trajectories due to different initial values, which are represented in black and red trajectories, respectively (see Fig. 3(g), Fig. 3(j), Fig. 3(m), Fig. 3(p), Fig. 3(r)).

The symmetric and asymmetric 1-$p$-$p$ periodic motions discussed above are all verified by the equivalent electronic circuit designed in Section 4.

This study reveals the nonlinear dynamics induced by clearances under constant excitation. Future research will build upon the approach outlined in Ref. [40], extending into two key directions: transient process analysis and fault-coupled modeling.

5. Equivalent circuit model and simulation

5.1. Linear module design

In the previous section, the basic feature of 1-$p$-$p$ response and chattering characteristics of the system with elastic constraints are analyzed through the use of bifurcation diagrams and phase portraits based on Poincaré mapping, as well as Lyapunov exponent. This method not only requires high programming ability, but also high-precision numerical simulation is limited by the operation performance of the computer, and the operation simulation speed is slow, which belongs to static calculation. As the system parameters are changed, the program needs to be recalculated, which takes up a large amount of computer capacity, consumes a lot of time and has poor economy. In this section, the equivalent electronic circuit model of the system shown in Fig. 1 is constructed using NI Multisim 14.0 software. Owing to its extensive component library and highly flexible parameter configuration capabilities, this software is particularly suitable for modeling vibration system circuits that involve nonlinear characteristics. The simulation method based on the equivalent circuit not only enables effective analysis of nonlinear dynamic behavior but also offers advantages such as high computational efficiency and the ability to dynamically adjust parameters.

In order to design the fully equivalent electronic circuit corresponding to the Eq. (1), the variable transformation $v_1=X_1$, $v_2={\dot{v}}_1$, $v_3=X_2$, $v_4={\dot{v}}_3$ is introduced, and the Eq. (1) can be transformed into the following state equations by reducing the order:

10

$\left\{\begin{array}{l}{v}_2={\dot{v}}_1,\\ P_1\sin (\mathrm{\Omega }T+\tau )-M_1{\dot{v}}_2-{\overline{C}}_1(v_2-v_4)-K_1(v_1-v_3)-F(v_1)=0,\\ v_4={\dot{v}}_3,\\ P_2\sin (\mathrm{\Omega }T+\tau )-M_2{\dot{v}}_4-{\overline{C}}_2{v}_4-{\overline{C}}_1(v_4-v_2)-K_2{v}_3-K_1(v_3-v_1)=0,\end{array}\right.$

in which, the piecewise linear spring force $F\left(v_1\right)$ can be rewritten as:

11

$F\left(v_1\right)=\left\{\begin{array}{ll}{K}_0\left(v_1-B\right),& v_1>B,\\ 0,& \left|v_1\right|<B,\\ K_0(v_1+B),& v_1<-B.\end{array}\right.$

In order to reduce the noise, the TL074CN integrated operational amplifier is selected here, and the integration circuit, addition circuit and amplification circuit composed of resistors, capacitors and other components are combined to realize Eq. (10). The schematic diagram of the designed electronic circuit is shown in Fig. 4, the structure of which is described as follows. Firstly, four integral operational amplifier circuits, i.e., $\text{U1A}$, $\text{U1C}$, $\text{U2C}$ and $\text{U2D}$, are the main components of the circuit, which are expressed as integral treatment of four parameters, i.e., $v_1$, $v_2$, $v_3$ and $v_4$, respectively; secondly, $\text{U2A}$, $\text{U1B}$, $\text{U2B}$ and $\text{U1D}$ are amplification circuits, which represent the reverse amplification of the four variables $v_1$, $v_2$, $v_3$ and $v_4$, respectively; thirdly, $\text{U2C}$ and$\text{U2D}$ are integral summation circuits.

Fig. 4Linear circuit design diagram

a) Integrating circuit

b) Summing amplifier

In order to determine the values of capacitors ($C_1$, $C_2$, $C_3$ and $C_4$) and resistors ($R_1-R_{12}$, $R_{100}-R_{107}$ and $R_{50}$) of the circuit in Fig. 4, we introduce Kirchhoff's law and active electronic circuit, then, Eq. (10) can be transformed into the following circuit state space equations:

12

$\left\{\begin{array}{l}{R}_1{C}_1{\dot{v}}_1=v_2,\\ {\dot{v}}_2+\frac{1}{R_2{C}_2}{v}_2-\frac{1}{R_3{C}_2}{v}_4+\frac{1}{R_4{C}_2}{v}_1-\frac{1}{R_5{C}_2}{v}_3+\frac{1}{R_{50}{C}_2}F\left(v_1\right)-\frac{V_{M1}}{R_6{C}_2}\sin (\mathrm{\Omega }T+\tau )=0,\\ R_7{C}_3{\dot{v}}_3=v_4,\\ {\dot{v}}_4-\frac{1}{R_8{C}_4}{v}_2+\frac{1}{R_9{C}_4}{v}_4-\frac{1}{R_{10}{C}_4}{v}_1+\frac{1}{R_{11}{C}_4}{v}_3-\frac{V_{M2}}{R_{12}{C}_4}\sin (\mathrm{\Omega }T+\tau )=0.\end{array}\right.$

Transform Eq. (12) into the following second-order differential equations:

13

$\left\{\begin{array}{l}{R}_1{C}_1{\ddot{v}}_1+\frac{R_1{C}_1}{R_2{C}_2}{\dot{v}}_1-\frac{R_7{C}_3}{R_3{C}_2}{\dot{v}}_3+\frac{1}{R_4{C}_2}{v}_1-\frac{1}{R_5{C}_2}{v}_3+\frac{1}{R_{50}{C}_2}F\left(v_1\right)\\ -\frac{V_{M1}}{R_6{C}_2}\sin (\mathrm{\Omega }T+\tau )=0,\\ R_7{C}_3{\ddot{v}}_3-\frac{R_1{C}_1}{R_8{C}_4}{\dot{v}}_1+\frac{R_7{C}_3}{R_9{C}_4}{\dot{v}}_3-\frac{1}{R_{10}{C}_4}{v}_1+\frac{1}{R_{11}{C}_4}{v}_3-\frac{V_{M2}}{R_{12}{C}_4}\sin (\mathrm{\Omega }T+\tau )=0.\end{array}\right.$

In order to make Eq. (1) and Eq. (10) completely equivalent, the following relationship should be satisfied:

14

$\left\{\begin{array}{c}\frac{R_1{C}_1}{R_2{C}_2}=\frac{R_7{C}_3}{R_3{C}_2}=\frac{R_1{C}_1}{R_8{C}_4}={\overline{C}}_1,\\ R_4{C}_2=R_5{C}_2=R_{10}{C}_4=\frac{1}{K_1},\\ \frac{1}{M_1}=R_1{C}_1, \frac{1}{M_2}=R_7{C}_3, P_1=\frac{V_{M1}}{R_6{C}_2},\\ R_9{C}_4\ll R_3{C}_2, R_{11}\ll R_{10}.\end{array}\right.$

The values of the resistors and capacitors in Fig. 4 are determined as follows. First, based on the dimensionless system parameters given in Section 4, i.e., $b=\text{0.25}$, $f_{20}=\text{0.0}$, $c=\text{0.5}$, $k=\text{0.5}$, $k_0=\text{0.95}$, $m=\text{0.5}$, $\zeta =\text{0.1}$, and by introducing state variables to transform Eq. (1) into state-space form, the initial component values are determined as $R_1=R_6=R_7=R_{12}=$ 100 kΩ and $C_1=C_3=$ 1 nF, $C_2=C_4=$100 nF. Next, according to the proportional relationships of the coefficients in Eq. (19), the remaining resistor values are calculated as $R_2=R_3=R_8=$ 22.37 kΩ, $R_4=R_5=R_{10}=$ 20 kΩ, $R_9=$ 11.19 kΩ, $R_{11}=$ 10 kΩ, $R_{50}=$ 20.88 kΩ. Finally, for the operational amplifiers $\text{U2A}$, $\text{U1B}$, $\text{U2B}$ and $\text{U1D}$, which are used solely for signal inversion, the peripheral resistors are uniformly set as $R_{101}=R_{102}=R_{103}=R_{104}=R_{105}=R_{106}=R_{107}=$ 100 kΩ.

According to the proportional relationship of coefficients $b=(B\times K_1)/(P_1+P_2)$, $P_2=V_{M2}/(R_{12}\times C_4)$, $K_1=1/(R_4\times C_2)$, it can be concluded that the voltage amplitude $V_{M2}$ of the sine signal is related to the clearance threshold $b$. Therefore, we can obtain:

15

$V_{M2}=\frac{B{R}_6}{b{R}_4}=\frac{0.064\times 100}{0.25\times 20}=1.28 \mathrm{V}.$

The correlation between the input frequency $f$ and the excitation frequency $\omega$ is:

16

$\omega =\mathrm{\Omega }\sqrt{\frac{M_1}{K_1}}=\pi f\sqrt{R_1{C}_1{R}_4{C}_2},$

i.e., $f=357.419\omega$.

5.2. Nonlinear module design

The nonlinear factor of the system is derived from the collision of the oscillator $M_1$ at the elastic constraints $A_1$ or $A_2$, which is described by three linear functions, see Eq. (17):

17

$f\left(x\right)=\mu x-0.5\mu (\left|x+b\right|-\left|x-b\right|)=\left\{\begin{array}{c}\mu (x-b), x>b,\\ 0, \left|x\right|<b,\\ \mu (x+b), x<-b,\end{array}\right.$

where, $\mu =k_0/(1-k_0)$.

Therefore, the linear saturation characteristics of the operational amplifier can be utilized to design the nonlinear circuit represented by Eq. (17). As shown in Fig. 5, the nonlinear circuit is composed of two circuit modules with the same magnification, opposite phase and different saturation parameters, which are linearly superimposed with the addition circuit. $\text{U3B}$ and $\text{U3C}$ are two in-phase operational amplifier units, and the output characteristic curves are shown in Fig. 6(a); $\text{U3A}$ is a reverse operational amplifier unit, and the output characteristic curve is shown in Fig. 6(b); $\text{U3D}$ is an additive operational amplifier unit. By utilizing the different saturation parameter intervals in Figs. 6(a)-6(b), Eq. (17) can be achieved through linear superposition of the additive operational amplifier unit $\text{U3D}$, as shown in Fig. 6(c).

Fig. 5Non-linear circuit design diagram

Fig. 6Piecewise linear function implementation diagram

a) Phase amplification

b) Reverse amplification

c) Superimposed output

The power supply voltage of the integrated operational amplifier TL074CN is selected as ±15 V, and the saturation voltage is selected as $E_{sat}=$ ±13 V. Then, we can obtain $R_{200}=$ 5 kΩ, $R_{201}=$ 100 kΩ, $R_{202}=$ 1 kΩ, $R_{203}=$ 10.53 kΩ, $R_{204}=$9 kΩ, $R_{205}=$1 kΩ, in which resistors $R_{204}$ and $R_{205}$ are used for voltage division. It can be calculated that the amplification factor $A_{+}$ of the forward amplifier circuit is $A_{+}=(R_{201}\times R_{203})/(R_{200}\times R_{202})=$ 210.6, and the output $v_{11}$ is:

18

$v_{11}=\left\{\begin{array}{c}\frac{R_{205}}{R_{204}+R_{205}}{E}_{sat},\\ \frac{R_{201}}{R_{200}}\frac{R_{203}}{R_{202}}\frac{R_{205}}{R_{204}+R_{205}}{E}_{sat}\\ -\frac{R_{205}}{R_{204}+R_{205}}{E}_{sat},\end{array}\right.,$

that is:

19

$v_{11}=\left\{\begin{array}{ll}\frac{1}{10}{E}_{sat}, & v_1>\frac{E_{sat}}{A_{+}},\\ 21.06v_1,& \left|v_1\right|\le \frac{E_{sat}}{A_{+}},\\ -\frac{1}{10}{E}_{sat},& v_1<-\frac{E_{sat}}{A_{+}}.\end{array}\right.$

In order to make the inverter $\text{U3A}$ achieve saturation, the inverse magnification is $A\_=R_{207}/R_{206}=$ 21.06, from which it can be inferred that $R_{206}=$ 1 kΩ, $R_{207}=$ 21.06 kΩ. Due to the fact that the saturated input voltage of $\text{U3A}$ is much higher than the actual input voltage of $v_1$, only the majority of the output is taken, and the inverted output $v_{12}$ is $v_{12}=-A\_\times v_1=-21.06v_1$, and $-E_{sat}/A\_\le v_1\le E_{sat}/A\_$. Finally, the sum of $v_{11}$ and $v_{12}$ is completed through the adder $U3D$, and the inverted output $f\left(v_1\right)$ can be obtained as follows ($R_{208}=R_{209}=R_{210}$):

20

$f\left(v_1\right)=\left\{\begin{array}{ll}21.06v_1-\frac{E_{sat}}{10},& v_1>\frac{E_{sat}}{A_{+}},\\ 0,& \left|v_1\right|\le \frac{E_{sat}}{A_{+}},\\ 21.06v_1+\frac{E_{sat}}{10},& v_1<-\frac{E_{sat}}{A_{+}}.\end{array}\right.$

Based on the design of the nonlinear equivalent circuit in Eq. (20), the output characteristics of the piecewise linear function $f\left(v_1\right)$ simulated in the MultiSim14.0 environment are shown in Fig. 7, where the abscissa is $v_1$ and the ordinate is $f\left(v_1\right)$.

Fig. 7Output characteristics of nonlinear circuit

5.3. MultiSim14.0 simulation

Fig. 8 is the circuit phase plane portraits simulated in MultiSim14.0 environment based on the equivalent circuit. The input voltage is $V_{M2}=$1.28 V and the input frequency is $f=357.419\omega$. We can observe that the simulation results obtained in MultiSim14.0 environment are in good agreement with the numerical simulation results in Fig. 3.

Fig. 8Circuit phase plane portraits, b= 0.25

a) 1-1-1_S, $f=$ 554 Hz

b) 1-1-1_AS, $f=$ 511 Hz

c) 1-2-2_S, $f=$ 304 Hz

d) 1-2-2_AS, $f=$ 268 Hz

e) 1-3-3_S, $f=$ 193 Hz

f) 1-3-3_AS, $f=$184 Hz

g) 1-4-4_S, $f=$ 146 Hz

h) 1-4-4_AS, $f=$ 142 Hz

i) 1-5-5_S, $f=$ 120 Hz

j) 1-5-5_AS, $f=$ 114 Hz

k) 1-6-6_S, $f=$ 102 Hz

l) 1-6-6_AS, $f=$ 96 Hz

m) 1-7-7_S, $f=$ 86 Hz

n) 1-8-8_S, $f=$78 Hz

In addition, symmetric and asymmetric 1-$p$-$p$ responses, as well as periodic motions with large number of impacts, such as 1-7-7 and 1-8-8 responses, can be clearly captured. The test results demonstrate the advantages of the circuit simulation approach: on identical hardware, its computational time is only one-third of that required by conventional numerical solution methods, while entirely eliminating the cumbersome steps of modifying code and recompiling after parameter adjustments.

6. Conclusions

A combined approach of the fourth-order Runge-Kutta method and equivalent circuit simulation is utilized to explore the mechanisms of the 1-$p$-$p$ response and chattering-impact phenomena in a two-degree-of-freedom vibro-impact system with elastic constraints.

1) In the case of small clearance, the system undergoes the following changes as $\omega$ changes from high to low: (i) symmetric 1-$p$-$p$ response occurs Pitchfork bifurcation to produce asymmetric 1-$p$-$p$ response, which has two antisymmetric trajectories due to different initial values; (ii) Asymmetric 1-$p$-$p$ response is embedded into chaos through periodic-doubling bifurcation sequence; (iii) The system degenerates from chaotic to $1-(p+1)-(p+1)$ response; (4) as $p$ increases sufficiently, the system is characterized by chattering-impact response. However, 1-$p$-$p$ response may undergo Grazing bifurcation, resulting in $1-(p+1)-p$ response.

2) Based on the mechanical model, an equivalent circuit model is built through the integral circuit, inverse circuit and addition circuit designed by common electronic components such as resistor, capacitor and integrated operational amplifier. The nonlinear circuit is designed by using the amplification and saturation characteristics of the operational amplifier. The numerical simulation output results obtained by fourth-order Runge-Kutta integration method are completely consistent with the circuit simulation output results. Compared with numerical simulation, especially high-precision numerical simulation algorithms, circuit simulation is not only fast, efficient, but also economical. It has low requirements on computer performance and can realize real-time adjustment of parameters. The dynamic characteristics of the system under different excitation frequencies and clearance thresholds can be obtained by changing the input frequency and voltage amplitude. The equivalent circuit simulation can be further extended to the dynamics analysis of mechanical systems with multiple-clearances and multiple-degrees-of-freedom, providing a novel and feasible reference idea and technical approach for the dynamic analysis of various mechanical devices involving nonlinear factors.

3) In the follow-up arrangement of the research content, we can engrave the hardware PCB circuit according to the designed electronic circuit schematic diagram, and realize the experimental verification of the hardware circuit for the vibration system with clearances. In the design and manufacturing of PCBs, the following issues need to be considered to ensure the stability and reliability of the circuit. For example, fully considering the clearance between components to prevent short circuits or series problems caused by too close spacing between components during the welding process; It is necessary to shorten the length of the wires as much as possible, avoid excessive through holes, and increase the difficulty of PCB soldering. In the field of mechanical engineering, Men et al. [41] successfully adapted simulated vibration data to real-world scenarios via style transfer technology to address the small-sample problem in bearing fault diagnosis, and their research paradigm is highly consistent with this work.