Contact Impact Forces at Discontinuous 2-DOF Vibroimpact

Nonlinear problems are arising in many different domains of science and engineering. Often they are modeled using sets of ordinary differential equations with discontinuous right-hand side. For example they are the systems with mechanical impacts, stick-slip motion from friction, electronic switches, hybrid dynamics in control, and genetic networks [3]. Vibroimpact system is one example of such systems. Vibroimpact system is strongly nonlinear one; the set of its motion differential equations contains the discontinuous right-hand side. Many new phenomena unique to non-smooth systems are observed under variation of system parameters. Jumps and switches in a system state represent the grossest form of nonlinearity [4].

Recently the investigations of such systems are developed rapidly. The short literature survey was made in [2]. The investigations of stability in systems with impacts, periodic motions, bifurcations, singularities at vibroimpact dynamics were noted in this review. The works devoted to learning of discontinuous dynamical systems [4–16] were discussed shortly.

Ukrainian mathematician Professor, Academician A.N. Sharkovsky was at the outset of dynamical systems theory [17]. The powerful scientific school of nonlinear vibrations under Professor Yu.V.Mikhlin leadership is working in Kharkiv, in Ukraine [18].

Discontinuous bifurcations are the hard bifurcations. The hard bifurcations were the subject of Catastrophe theory. Catastrophe theory was introduced in the 1960s by the renowned Field Medal mathematician Rene Thom as a part the general theory of local singularities [19]. Since then it has found applications across many areas, including biology, economics, and chemical kinetics. By investigations the phenomena of bifurcation and chaos, Catastrophe theory proved to be fundamental to the understanding of qualitative dynamics. The famous books [20,21] are devoted to this topic. The theory was very fashionable at 70th years of 20th century. Then this fashion went away and terminology from catastrophe returned to singularities, discontinuous bifurcations and so on. But the terms "catastrophe, catastrophe point" are used now too [22–25]. The Catastrophe theory is used in works of Saratov (in Russia) scientific school [26]. This school is strong scientific school under Professor S.P. Kuznetsov leadership [27, 28]. There is another powerful scientific school in Nizhny Novgorod (in Russia). Professor L.P. Shilnikov was its leader during many years. Particularly Blue Sky Catastrophes are learnt by its scientists. Professor V.S.Afraimovich is the fosterling of this scientific school; Professor V.I.Nekorkin is its bright representative.

We investigate the dynamic behaviour of vibroimpact systems by parameter continuation method. Short review was made in [2]. The works [29, 30] were discussed in this survey. We simulate impact by nonlinear contact interaction force according to Hertz’s contact law. Such simulation gives us the possibility to find the motion law along the whole timebase including the impact phase, to determine the impact duration and to find the contact impact forces. The contact impact force is one of the most significant characteristics of vibroimpact system motion. Its determination is necessary for fulfillment of the strength calculations in engineering.

Its dynamic behaviour is determined by the system state and motion. Naturally that contact force has qualitative change there where the motion regime has the bifurcations. The dynamic behaviour of contact forces in strongly nonlinear discontinuous vibroimpact system have been learned insufficiently. We have seen little of such papers in scientific literature. At [31] the experimental investigations and the mathematical modeling of the impact force behaviour in a vibroimpact system are described. This is system where an impact pendulum is mounted on a cart that moves with a prescribed displacement. Some approach for determining the forces of impact interaction in vibroimpact systems is proposed at [32]. A short description of different impact force models is given at [33].

The aims of this paper are:

So far as this paper is the continuation of works [1, 2] the problem formulation is the same. We’ll repeat it shortly.

We analyze the dynamic behaviour of contact impact forces for discontinuous nonlinear vibroimpact system presuming it is a two-body two-degree-of-freedom one (Fig. 1).

Fig. 1

This vibroimpact system is formed by the main body and attached one, and the latter can play the role of percussive or non-percussive dynamic damper. Bodies are connected by linear elastic springs and dampers. The main body is under the effect of periodical external force:

F(t)=Pcos(ωt+φ0).$$\begin{array}{} \displaystyle F(t) = Pcos(\omega t + {\varphi _0}). \end{array}$$(1)

We consider impacts as low velocity elastic collinear collisions without friction. The contact surfaces are smooth curvilinear ones without roughness. Thus real surface geometry in contact zone may be approximated by "Herzian" geometry.

The initial point of x coordinate is chosen in the main body mass center at the moment when all springs are not deformed. The initial distance between bodies at this moment is D. The structure of the system is experiencing transformation during the movement. The reason is its dynamic states modification forced by the impact contacts between elements.

Motion equations of the system have got the form:

x1..=−2ξ1ω1x˙1−ω12x1−2ξ2ω2χ(x˙1−x˙2)−−ω22χ(x1−x2+D)+1m1[F(t)−Fcon(x1−x2)],x2..=−2ξ2ω2(x2.−x1.)−ω22(x2−x1−D)+1m2Fcon(x1−x2),$$\begin{array}{l} \displaystyle \begin{array}{*{20}{l}} {\mathop {{x_1}}\limits^{..} = }&{ - 2{\xi _1}{\omega _1}{{\dot x}_1} - \omega _1^2{x_1} - 2{\xi _2}{\omega _2}\chi ({{\dot x}_1} - {{\dot x}_2}) - }\\ {}&{\quad \quad - \omega _2^2\chi ({x_1} - {x_2} + D) + \frac{1}{{{m_1}}}[F(t) - {F_{con}}({x_1} - {x_2})],}\\ {\mathop {{x_2}}\limits^{..} = }&{ - 2{\xi _2}{\omega _2}(\mathop {{x_2}}\limits^. - \mathop {{x_1}}\limits^. ) - \omega _2^2({x_2} - {x_1} - D) + \frac{1}{{{m_2}}}{F_{con}}({x_1} - {x_2}),} \end{array} \end{array}$$(2)

where ω1=k1m1$\begin{array}{} \displaystyle {\omega _1} = \sqrt {\frac{{{k_1}}}{{{m_1}}}} \end{array}$,ω2=k2m2$\begin{array}{} \displaystyle \omega_2=\sqrt{\dfrac{k_2}{m_2}} \end{array}$; ξ1=c12m1ω1$\begin{array}{} \displaystyle \xi_1=\dfrac{c_1}{2m_1\omega_1} \end{array}$, ξ2=c22m2ω2$\begin{array}{} \displaystyle \xi_2=\dfrac{c_2}{2m_2\omega_2} \end{array}$; χ=m2m1$\begin{array}{} \displaystyle \chi=\dfrac{m_2}{m_1} \end{array}$; ω₁,ω₂ – partial oscillation frequencies,F_con(x₁ − x₂) - contact interaction force, it is simulating the impact and working only during the impact.

Initial conditions are:

x1(0)=0,x2(0)=D,x˙1(0)=0,x˙2(0)=0.$$\begin{array}{} \displaystyle {x_1}(0) = 0,{x_2}(0) = D,{\dot x_1}(0) = 0,{\dot x_2}(0) = 0. \end{array}$$(3)

We considered in detail the impact simulation manner in [34, 35]. To simulate the impact here we use the Hertz’s contact interaction force based on quasistatic Hertz’s theory [36, 37]:

Fcon(x1−x2)=K[(x1−x2)H(x1−x2)]3/2,K=43q(δ1+δ2)A+B,δ1=1−ν12E1π,δ2=1−ν22E2π,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{F_{con}}({x_1} - {x_2})} \hfill & = \hfill & {K{{[({x_1} - {x_2})H({x_1} - {x_2})]}^{3/2}},} \hfill \\ {K = \frac{4}{3}\frac{q}{{({\delta _1} + {\delta _2})\sqrt {A + B} }},{\delta _1}} \hfill & = \hfill & {\frac{{1 - {\nu _1}^2}}{{{E_1}\pi }},{\delta _2} = \frac{{1 - {\nu _2}^2}}{{{E_2}\pi }},} \hfill \\ \end{array} \end{array}$$(4)

where (x₁ − x₂) is the relative bodies rapprochement due the local deformation in contact zone, H(x₁ − x₂) is the Heaviside step function, v_i and E_i are respectively Poisson’s ratios and Young’s modulus for both bodies, A, B and q are the geometry characteristics of contact zone. We consider these surfaces as spherical ones, then A=B=12R1+12R2$\begin{array}{} \displaystyle A=B=\dfrac{1}{2R_1}+\dfrac{1}{2R_2} \end{array}$, where R₁,R₂ are the contact surfaces radiuses. Only local deformations in contact zone are taken into account by the Hertz’s theory.

There are works where proposals of different ways for making Hertz formula more precise are given [38–40]. For instance, at [38] different contact-impact force models for both spherical and cylindrical shape surfaces collisions in multi-body systems were reviewed. Various types of friction force models based on the Coulomb law were also listed and discussed. In [41] the authors examine limitations of the Hertz’s theory use for different individual cases. For example, the contact surfaces aren’t Herzian; impact velocities are large and the plastic deformations occur; there is conformal contact between bodies; the impacts aren’t pure elastic and energy dissipation is necessary to be taken into account. In [41] there is a set of contact forces models presented that were making the Hertz’s theory wider and more precisely defined. They complete contact force expression by additional terms that take into account damping and energy dissipation during the impact. In [42] the authors consider the evolution from the Hertzian contact model to non-Hertzian conditions for fast dynamic simulation.

Nevertheless Hertz’s theory is widely used for analysis of vibroimpact system dynamics now too. Just impact simulation by nonlinear contact interaction force allows to find the motion law at all timebase including impact phase, to define impact duration and contact forces values.

Let us note that less rough impact models use for vibroimpact system (for example the wave theory) causes considerable difficulties due to repeated impacts.

We don’t introduce the non-dimensional variables into the motion Eq. (2) having a hope the dimensional characteristics will clearly demonstrate the dynamical behaviour of specific vibroimpact system model with defined numerical parameters.

Numerical analysis of dynamic behaviour of contact impact forces was fulfilled by parameter continuation method in conjunction with shooting and Newton-Raphson methods [1]. Periodic motion stability or instability is determined by matrix monodromy eigenvalues that is by Floquet multipliers’ values. The periodical solution is becoming unstable one if even though one Floquet multiplier leaves the unit circle in complex plane that is its modulus becoming more than unit. Such multiplier value characterizes the bifurcation kind of this bifurcation point.

We have described the theoretical basis for analysis of two-body two-degree-of-freedom system in [1]. Numerical system parameters are given further in Table 6.

3.1

Contact forces behaviour depending on excitation amplitude

Global view of contact forces behaviour is presented at Fig. 2 in wide range of excitation amplitude.

Fig. 2

The points O, B, C, D, K, L, N - are the bifurcation points. At these points the motion regime changes its stability (or instability), new regimes - (5,2), (2,3), (4,6)

We call nT -periodic regime with k impacts per cycle as (n,k)-regime in [43].

, quasiperiodic - are rising at these points. First of all one can see that the values of contact forces at new regimes are larger (sometimes considerably larger) than ones at the main (1,1)-regime. Contact forces at these regimes are depicted at Fig. 3. Its phase portraits are given to the right of this plot. The excitation force is plotted by green colour at all graphs. The time origin is any point in steady-state regime. The numerical coefficient before excitation amplitude is introduced in order to have the possibility to show it at the same plot. We see the contact force increasing more than an order of magnitude at (2,3)-regime. The contact forces for impacts occuring during 9 s from 0 to 9 s at quasiperiodic regime (see Fig.3,(b)) are depicted at Fig. 4.

Fig. 3

Fig. 4

Point O at Fig. 2 is the point of discontinuous bifurcation. It is phenomenon unique for non-smooth nonlinear system whose equations have discontinuous right-hand side. The system motion was impactless under small values of excitation amplitude. The contact force is absent under impactless motion. At point O stable impactless motion became unstable (1,1)-regime with one impact per cycle. Other regimes - stable and unstable branches of (5,2)-regime - are arising here. The set-valued Floquet multipliers cross the unit circle by jump that is their moduli becoming more than unit by jump. Fig. 5 shows well these jumps for multipliers µ₁ and µ₂ . Table 1 shows these jumps by numbers.

Fig. 5

Table 1

P	Re(μ1)	Im(μ1)	|μ1|	Re(μ2)	Im(μ2)	|μ1|
98.48	-0.25	0.83	0.87	-0.25	0.83	0.87
98.98	-0.25	0.83	0.87	-0.25	0.83	0.87
99.48	3.57	0	3.57	1.12	0	1.12

Let us note by the way that Floquet multipliers during the impactless motion under small excitation amplitude values have constant values (Fig. 5).

At bifurcation point C the transition of multipliers’ moduli via unit is smooth one (Table 2). Here two conjugate multipliers leave the unit circle fluently. It is the point of quasiperiodic motion arising whose phase trajectories take the form of torus.

Table 2

P	Re(μ1)	Im(μ1)	|μ1|	Re(μ2)	Im(μ2)	|μ1|
392.5	0.197	0.977	0.997	0.197	-0.977	0.997
392.0	0.197	0.981	1.001	0.197	-0.981	1.001
391.0	0.197	0.983	1.003	0.197	-0.983	1.003

3.2

Contact forces behaviour depending on excitation frequency

Global view of contact forces behaviour is presented at Fig. 6 in wide range of excitation frequency.

Fig. 6

One can see several changes of motion stability and several bifurcation points. They are the points A, B, C, D, E, K, L, and M. Other regimes are arising at these points: (3,1), (4,2), qusiperiodic, chatter regimes. The contact force at chatter regimes (nT -periodic with multiple impacts per cycle) is depicted by dark green circles. We see that contact force is increasing at such regimes. The points A and M are the discontinuous bifurcation points. At these points the impactless motion is changed at the motion with impacts. At point M we observe phenomenon unique to discontinuous system - discontinuous fold bifurcation. The discontinuous fold bifurcation connects a stable branch to an unstable branch. Here set-valued Floquet multiplier makes huge jump along the positive real axis. Its motion along positive real axis is demonstrated by Table 3.

Table 3

ωi, rads.s-1	8.03	8.04	8.05	8.06	8.07	8.10	8.16
Re(μ1)	-0.595	-0.593	151.4	90.1	65.9	37.4	19.7
IM(μ1)	0.654	0.652	0	0	0	0	0
|μ1|	0.8828	0.8829	151.4	90.1	65.9	37.4	19.7

3.3

Contact forces behaviour depending on stiffness characteristics

Naturally the stiffness characteristics of vibroimpact system exert influence at contact forces values. We have found the contact forces dependence on connecting spring stiffness k₂ (see Eq. (1) and Fig. 1), bodies’ materials, and contact surfaces geometry. The contact force graphs during one impact under different stiffless values are depicted at Fig. 7.

Fig. 7

One can see that the impact force weakens when spring stiffness increases. Under small spring stiffness motion regime becomes nT -periodic with multiple impacts - chatter. It is shown at Fig. 8 under k₂=1050 N·m^-1.

Fig. 8

In order to see the influence of bodies’ materials we supposed that material of both bodies was identical that is their Yung’s moduli were equal E₁=E₂. Then we have obtained results that are shown at Table 4.

Table 4

Material	Steel	Copper	Aluminium	Rubber
Young’s modulus, N m-2	2.10	1.11	0.69	0.00008
Maximum Fcon,N·105	1.71	1.33	1.10	0.0309

Let us note the interesting result: the rubber elasticity modulus is less at 26000 times than steel elasticity modulus. The contact force is less at 55 times and main body vibration amplitude is less only at 1.22 time!

Contact Hertz’s theory supposes that contact surfaces are curvilinear "Hertzian" ones [36]. The local geometry of contact zone characterises by constants A, B, q that are inside the formula for Hertz contact force ( Eq.(4)). We consider the contact surfaces are spherical ones with radiuses R₁ and R₂. Let these radiuses are equal ones, then we have the results presented at Table 5.

Table 5

R1=R2, m	A=B, m-1	Fcon,N·105
1	1	1.49
2	0.5	1.71
5	0.2	2.05
10	0.1	2.36

It is interesting that hereat main body vibration amplitude isn’t change practically!

Thus we can get the decreasing (or increasing) of contact force by changing the system stiffness parameters.

3.4

Contact forces behaviour under change of impact kind

There is the vibroimpact system classification by different aspects [44]. One of them is impact kind characteristic - rigid or soft impact. Some principal difference between rigid and soft contact were formulated in [45]. The main sign is its duration. Just impact duration dictates the way of its simulation. If impact duration is large then impact isn’t instantaneous. Its simulation by boundary conditions with Newton’s restitution coefficient using based on stereomechanic theory isn’t possible [34,35]. The stiffness of vibroimpact system elements causes the impact softness. The soft impacts take place in engineering very often. The authors write in [46]: "Soft impacts occur in many practical mechanical systems where there is some "cushioning" at the impacting surfaces - meant for reducing the noise and chatter. It can be visualized as a mass impacting not with a hard wall, but with a spring-damper support in front of a wall. The existence of the spring-damper type cushion introduces some special features in the system dynamics".

The investigation of systems which constituted by softly impacting beams and rods of non-negligible mass is fulfilled in [47]. The impact is simulated by force which depends linearly on bodies penetrations one into another. Some virtual spring imitates this force. The values of bodies’ displacements and penetrations characterize the impact softness. The authors investigate the influence of virtual spring stiffness at system dynamic behaviour and vibratory motion character.

What influence will impact kind change exert at contact force behaviour?

The clear criterion of impact rigidness or softness is absent. The typical trait of impact softness is its duration. Is it instantaneous or not? Let us examine the value - the coefficient of the relative impact duration kcon=TconT⋅100%$\begin{array}{} \displaystyle % MathType!Translator!2!1!LaTeX.tdl!LaTeX 2.09 and later! {k_{con}} = \frac{{{T_{con}}}}{T} \cdot 100\% % MathType!End!2!1! \end{array}$ at (1,1)-regime. Here T is period of external loading ( Eq. (1)), T = 2π·ω⁻¹, T_con is the time of impact that is the time of contact between bodies. The impact duration is T_con = 0.782 ms and the coefficient of relative impact duration is k_con = 0.09 % for the motion with rigid impact at (1,1)-regime. Rigid impact is almost instantaneous because its duration is very small. The graph of contact impact force in dependence on time in every impact has the form of "stick" (Fig. 3, Fig. 7). How can we change the vibroimpact system parameters for decrease the coefficient of relative impact duration k_con i.e. for decrease the impact duration T_con ?

We use the modified steepest descent method – gradient projection method with correction of residual with constraints [48,49]. We formulate the minimax problem like that: to find such parameters of vibroimpact system which will provide the biggest value of relative impact duration coefficient k_con that is the largest value of impact duration T_con .We seek for the objective function kcon=TconT⋅100%$\begin{array}{} \displaystyle {k_{con}} = \frac{{{T_{con}}}}{T} \cdot 100\% \end{array}$ by deciding the problem about steady-state exiting vibrations of vibroimpact system under the concrete system parameters. The calculation of objective function gradient takes into account the constraints laid on the parameters. The obtained value determines the next optimization step. We investigated such parameters influence at impact duration: the attached body’s mass m₂; the joining spring’s stiffness k₂ (it comes into the partial frequency expression Eq. (3)); Young’s modulus of both system bodies E₁, E₂; the radiuses of contact surfaces for both bodies R₁, R₂ in assumption that these surfaces are spherical (these radiuses comes into expressions of coefficients A and B in formula for Hertz’s force Eq. (2)). The other parameters weren’t changed. We have obtained such results. The impact duration was enlarged considerably under influence of these parameters changing [50, 51]. Table 6 shows system parameters under rigid and soft impacts and the results of such optimisation.

Table 6

Parameters of vibroimpact system

Bodies’ characteristic	Rigid impact		Soft impact
	Main body	Attached body	Main body	Attached body
Mass mi, kg	1000	100	1000	310
Partial vibration frequency ωi,rad·s-1	6.283	5.646	6.283	3.606
Young’s modulus Ei, N·m-2	2.1 · 1011	2.1 · 1011	2.44 · 105	2.1 · 1011
Contact surface radius Ri, m	2	2	1	0.5
Coefficients A,B, m-1	A = 0.5 B = 0.5		A = 1.5 B = 1.5
Poisson’s ratio vi	0.3	0.3	0.3	0.3
Damper coefficient ξi	0.036	0.036	0.036	0.036
Initial distance between bodies D, m	0.05		0.05
Impact duration Tcon, s	7.82 · 10-4		0.19
Coefficient of impact duration kcon, %	0.09		20.9

The impact becomes the soft one. The graph of contact impact force in dependence on time in every impact has the form of "bell" ( Fig. 9).

Fig. 9

At Fig. 9 one can see the graphs of contact forces dependence on time for rigid and soft impacts at (1,1)-regime. We see well the "stick" for rigid impact and the "bell" for soft one. We see also that contact force at rigid impact is more than 2 orders. The excitation force is plotted by green colour. The time origin is any point in steady-state regime. The numerical coefficient before excitation amplitude is introduced in order to have the possibility to show it at the same plot.

The contact force behaviour in dependence on excitation amplitude is shown at Fig. 10.

Fig. 10

At points C, D, B the unstable (1,1)-regime turns into stable one. The bifurcation point C is the point of period doubling. The (2,2)-regime is arising at this point. At point A the impactless motion turns into impact regime. The contact force behaviour in dependence on excitation frequency is shown at Fig. 11.

Fig. 11

The stable motion takes turns into unstable one at points C, D, L, M, N. The impact regime turns into impactless at points B, E, K. The (2,2)-regime is arising at this point L. The discontinuous bifurcations are absent under soft impact. The contact forces are considerably less than ones under rigid impact.