Asymmetric and chaotic responses of dry friction oscillators with different static and kinetic coefficients of friction

The responses of a simple harmonically excited dry friction oscillator are analysed in the case when the coefficients of static and kinetic coefficients of friction are different. One- and two-parameter bifurcation curves are determined at suitable parameters by continuation method and the largest Lyapunov exponents of the obtained solutions are estimated. It is shown that chaotic solutions can occur in broad parameter domains—even at realistic friction parameters—that are tightly enclosed by well-defined two-parameter bifurcation curves. The performed analysis also reveals that chaotic trajectories are bifurcating from special asymmetric solutions. To check the robustness of the qualitative results, characteristic bifurcation branches of two slightly modified oscillators are also determined: one with a higher harmonic in the excitation, and another one where Coulomb friction is exchanged by a corresponding LuGre friction model. The qualitative agreement of the diagrams supports the validity of the results.


Introduction
The effects of dry friction on the relative motion of contacting bodies are difficult to predict. On the one hand, the underlying physical mechanisms are very complex. On the other hand, various mathematical difficulties can appear during the analysis of frictional systems, that require special tools to solve [1]. Moreover, friction models applied to characterize the contact of rigid bodies can lead to indeterminate or self-inconsistent solutions [2,3].
This model is often examined in theoretical contributions (see e.g., [1,[16][17][18]), but there are several engineering applications, too-from wiper blades and cotton picking machines to ocean wave energy converters [19][20][21] where it can serve as a proper description of the dynamics.
The dimensionless equation of motion assumes the form where t ¼ s ffiffiffiffiffiffiffiffi ffi k=m p , x ¼ zk=F 0 , _ ðÞ denotes d/dt, and . The kinetic friction is characterized by S ¼ lmg=F 0 , while the static friction parameter is denoted by S 0 ¼ l 0 mg=F 0 [ S.
Coordinate x is measured from the untensioned state of the spring and the friction characteristic is described by the function where l and l 0 denote the kinetic and static coefficient of friction, respectively, and l 0 =l ¼ S 0 =S. Den Hartog [10], one of the pioneers in the field of the theory of vibrations, obtained the exact, nonsticking solution of the steady state motion of this oscillator, assuming that the solutions are symmetric. Shaw [11] made a big leap towards the understanding of frictional systems by considering different static and kinetic coefficients of friction and complementing Den Hartog's work with a linear stability analysis. Pratt and Williams [12] and Hong and Liu [13,14] extended the aforementioned studies by the numerical examination of solutions which have multiple sticking segments during a period. Natsiavas [15] and Deimling [22] discussed the stability of solutions-the latter study is based on the theory of differential inclusions (set-valued differential equations). The combined effects of different friction coefficients and viscous damping were discussed in [23], where the stability properties of various solution types were determined with the corresponding bifurcation diagrams, as well.
Our study provides an extension of these results by a more detailed study and identification of bifurcations and the examination of chaotic responses. Systems with Coulomb friction belong to the class of piecewise smooth systems that are usually analysed based on the Filippov theory [1,18]. However, the extension of Filippov's method to the case of different static and dynamic friction coefficients is not straightforward. This is why special methods are necessary for the analysis of this case, as it was shown in [24], where the so-called Switch Model was introduced for this purpose. In the present contribution, we use a properly modified variant of the Switch Model.
Several authors reported the appearance of numerically found asymmetric solutions (with different maximal displacements in the positive and negative directions), e.g., [11,23,24], especially at negative values of the viscous damping factor. The conditions of symmetric, non-sticking solutions at pure dry friction damping were derived in [25]. These results led to the analytical formulation of a family of asymmetric responses [26] at subharmonic resonances, where the ratio of the excitation frequency and the natural frequency is X ¼ 1=2, 1/4, 1/6,.... The possibility of chaotic responses in general frictional systems was examined in [27], based on the ideas laid down in [28]. It was pointed out in [29] that the nonreversible nature of dry friction can lead to chaotic vibrations.
Regarding system (1), it was shown in [30] that the parameter domain of asymmetric periodic solutions opens up if the static and dynamic coefficients of friction are different and the occurrence of chaotic and transient chaotic behaviour was pointed out at X ¼ 1=2 and S 0 =SJ8. Such a big ratio between the coefficients of friction is rather extreme in models of physical systems, thus, this finding seemed to be irrelevant. However, the recently published experiments of Marino and Cicirello [31] showed that the theoretically predicted asymmetric solutions can be detected in reality. Moreover, certain data sets could be interpreted as transient chaotic responses. These results motivated the present study. Fig. 1 Harmonically excited frictional oscillator. F 0 denotes the forcing amplitude, x 0 is the angular frequency of excitation, while m and k are the mass and spring stiffness, respectively. The displacement coordinate is z and s denotes the real, physical time Our goal is the extension of the bifurcation analysis reported in [30] and the determination of the connection between asymmetric and chaotic solutions. It is worth to note that in real physical systems the friction parameters vary stochastically during motion due to various mechanical and physical interactions. Thus, a potential sensible and up-to-date technique to the analysis of these systems is the application of statistical tools [32,33]. We opted another approach: to check the robustness of the qualitative results obtained for the harmonically excited dry-friction oscillator, characteristic bifurcation curves of two slightly modified oscillators were also determined. One of the oscillators exhibits a higher harmonic in the excitation, while Coulomb friction is exchanged by a corresponding LuGre friction model [34] in the other one. The obtained diagrams show good qualitative agreement with the corresponding diagrams of the original model. The LuGre model is more realistic than the Coulomb model (e.g., it can describe the socalled pre-sliding phenomenon) and it is often used in high-precision control applications [35]. Notably, changes in the friction force were found to follow the changes in relative velocity with some delay, especially in case of lubricated interfaces [36]. As a consequence of this frictional lag or frictional memory, the friction force is larger for increasing velocities than for decreasing ones. The LuGre model exhibits this type of hysteresis in the friction force-velocity relationship, but it is worth to mention that modified versions of the model can incorporate more refined lag/memory effects, too [37].
The present study is based on the reformulation of the system's description in Sect. 2 for the application of the continuation method. One and two-parameter bifurcation curves are determined in Sects. 3 and 4, respectively, supplemented by the estimation of the largest Lyapunov exponent in broad parameter domains. Section 5 deals with the bifurcation analysis of the aforementioned, slightly modified oscillators. Finally, Sect. 6 is devoted to the discussion of the results.

Reformulation of the problem
In order to reveal the bifurcations leading to the emerge of asymmetric solutions among the symmetric ones, the continuation method [38] was applied. The continuation method is capable of the quick, reliable calculation of bifurcation diagrams. One of the most widely used continuation software is AUTO [39]. Unfortunately, AUTO is not able to deal with nonsmooth systems. However, a recently implemented Matlab toolbox COCO (COntinuation COre, [40]) provides the possibility for the analysis of these kinds of systems, too.
To be able to analyse a piecewise smooth system with COCO, one must have an initial guess about the solution's structure. It means that the temporal order of the subsequent smooth segments in the periodic solution must be known, together with the switching conditions-formulated as zeros of smooth functions. Physically, the oscillator can be in one of the following three states: sliding with positive speed, sliding with negative speed or sticking, as it is formulated in Eqs. (1) and (2). Mechanically, Eq. (2) is interpreted as follows: if the velocity is non-zero (i.e., the body slips), the kinetic friction force always opposes the velocity. Thus, different differential equations describe the motion at positive and negative velocities. This is why the zero-velocity subset of the phase-space is called switching manifold.
The switching conditions between the sticking and sliding cases are given by an equation and a nonsmooth inequality: Thus, the reformulation of these conditions as zeros of smooth functions is not a trivial task. To examine similar systems, Leine and Nijmeijer use a special method in their book [24], by introducing the concept of the so-called Switch Model to be able to deal velocity-related and force-related switches separately.
To solve this problem, we used a similar approach: we partitioned the three-dimensional (x, v, t) phasespace into domains in such a way that the sign change of a smooth function (indicating that a solution segment crosses a boundary) uniquely determines the type of the solution in the next segment. Taking this requirement into account the phase-space was divided into 7 domains-6 different sliding states and a sticking state-according to the values of v and x. The corresponding domains are given in Table 1 and illustrated in Fig. 2.
The domains are separated by three switching surfaces that are also depicted in Fig. 2. Two of them separate the states with respect to the displacement xthese are the stick-slip boundaries Fþ and FÀ, respectively. The third switching surface is the v ¼ 0 plane. Thus, the corresponding event functions can be formulated as If the actual state of the system is known, the triggering of an event function unambigously indicates the next state that the system enters. For example, if the oscillator's state is x þ v þ (see Table 1), two possible events are possible: There are three possible switches from the state x 0 v þ : The remaining four sliding cases ( can be analysed similarly. If the oscillator is in the sticking state x 0 v 0 , the collision with the boundary FÀ means that the resultant of forces acting on the block becomes positive. Thus, the next state will be x À v þ . Similarly, crossing Fþ leads to x þ v À : These two switches can happen only if the static coefficient of friction S 0 is greater than the kinetic one. In the opposite case, sliding in the opposite direction occurs. This possibility will be excluded from the analysis due to its non-physical nature. Once the states and event functions are properly defined, the bifurcations of the oscillator can be explored by COCO. However, this process is not straightforward, since besides the conventional bifurcations (e.g., period-doubling), several other so-called discontinuity-induced or non-smooth bifurcations (e.g., grazing-sliding, crossing-sliding, switching-sliding or adding-sliding) can occur. These bifurcations are accompanied by the change of the type and order of segments (the signature of the solution) during a period. Thus, one has to make an educated guess about the possible order of states after each non-smooth bifurcation. Besides the bifurcations, other events can also take place that are related to the order of switches among the possible states shown in Table 1, but do not correspond to physically meaningful changes in the periodic solutions. For example, segments can appear or disappear when a sliding solution curve touches one of the stick-slip boundary surfaces Fþ or FÀ. This event is denoted by T in Fig. 4. The disappearance of segments could be predicted by monitoring the durations of solution segments. To detect the appearance of new segments inside an existing one, the validity of solutions had to be checked during each run. Our experiences show that bifurcations often lead to a completely different solution signature. While the local qualitative change of the trajectory can be known in these cases as well, the signature of the whole periodic orbit is required for the continuation, that is often difficult or practically impossible to predict. The For the initiation of the calculation, the analysis by COCO was accompanied by the solution of the initial value problem, i.e., numerical simulation. Clearly, only stable solutions can be found with this approach, but once a solution type is found by COCO, it can be continued even after the loss of stability. Besides the calculation of the numerical solution, several files must be initialized for COCO-these are also done automatically by our Matlab code. In addition to that, utility functions were implemented for the proper continuation from a previously obtained solution point, with the possibility to change the signature of the solution.
As it will be shown, four different periodic solution types can be detected in large domains of the examined parameter space. Each of these typical periods consists of 6 segments, thus, the number of segments will be denoted only if it is different from 6. The four typical solutions will be denoted by 2s, 0s, 2a and 1a, where 2, 1 or 0 shows the number of sticking segments in a period, while s or a tells whether the solution is symmetric or asymmetric.
The bifurcation points or curves will be denoted accordingly. For instance, SN 1a and PD 1a denote the saddle-node and period-doubling bifurcation of a 1stick asymmetric solution, respectively. Sliding bifurcations lead to the change of the number of sticking segments. In these cases only the symmetry property will be shown as CS s , SS a or AS a .
Since the symmetry plays an important role in the classification of solutions, the measure of asymmetry was introduced as follows: This quantity was also monitored during each continuation run.
3 One-parameter bifurcation analysis

Bifurcation parameter: S
To reveal the source of the asymmetric solutions, the dimensionless kinetic coefficient of friction S was chosen to be the first bifurcation parameter, while the other parameters were fixed at X ¼ 0:5 and S 0 ¼ 0:4, just as in [30]. An initial solution guess was generated by numerical simulation at these parameters in Matlab-see Fig. 3.
The following solution types were detected: • The initial guess is symmetric and has the following structure at S ¼ 0:4: Fig. 4) loses stability via a saddle-node (SN 2s ) bifurcation at S % 0:33358.
• Below this parameter, the initialization process found only a pair of stable solutions (denoted by 1a). These solutions are asymmetric with one sticking segment per period. One of them has the following structure: There are two sticking segments, both at x [ 0.
Thus, this solution is asymmetric. Consequently, its mirror image exists, too. The solution's structure breaks down at S % 0:05578, where the 11 th solution segment (x 0 v À ) just touches the stick-slip boundary F À . This event-denoted by T in the figures-is related to the change of the order and number of segments.
• Nevertheless, no stable periodic solution was found in the range S ¼ ½0:043479; 0:055787 by COCO. Interestingly, the numerical initial problem solver detected a seemingly periodic solution, but it turned out to be unstable according to COCO. This experience indicates that the found solution is only slightly unstable, such that the stabilizing effect of the numerical algorithm could overcome its instability. Since no attracting periodic solution was found, the existence of multi-periodic solutions or chaos is expected in this parameter domain. • In the range S ¼ ½0:035620; 0:043479, an asymmetric solution (1a,9seg) with one stick per period was detected. This solution emerges at S ¼ 0:043479 via a switching sliding (SS a ) bifurcation. It can be divided into 9 segments: solution could not be continued past the value S ¼ 0:035620 due to the divergence of the numerical method (see point N in Fig. 4). No stable periodic solution was found at even smaller values of S-here the occurrence of a chaotic domain is expected, again. The more detailed analysis of the parameter domains S ¼ ½0; 0:035620 and S ¼ ½0:043479; 0:055787 is postponed to Sect. 4.

Bifurcation parameter:
As the previous calculation showed, asymmetric solutions (denoted by 1a in the figures) appeared at S % 0:33358 at the discontinuity of the bifurcation diagrams. In the present section, we search for the origin of these asymmetric solutions. For this purpose, continuation was initiated in the vicinity of their appearance, at S ¼ 0:32 and S 0 ¼ 0:4, varying the • The continuation is initiated at X ¼ 0:5 with the numerically detected solution type This asymmetric solution (1a) is stable in the range X 2 ½0:49182; 0:54191. At the lower end, a new solution (2a) with two sticking segments appears via a degenerate crossing-sliding (CS a ) bifurcation, when the second segment (x þ v þ ) disappears as the trajectory touches the v ¼ 0 plane and the Fþ boundary simultaneously. This bifurcation is asymmetric since no such event occurs at FÀ. The stability of solution 1a is lost via a saddlenode/fold (SN 1a ) bifurcation at X ¼ 0:54191. Finally, another crossing-sliding bifurcation (CS s ) takes place at X ¼ 0:53529, where the solution becomes symmetric and both sticking segments touch the v ¼ 0 plane and the Fþ or FÀ boundary simultaneously.
• Solution 2a, that is born at X ¼ 0:49182 via crossing-sliding, has the following structure: This solution branch undergoes a saddle-node (SN 2s ) bifurcation at X ¼ 0:48874. Here the solution is just symmetric, thus, this bifurcation is the crossing point between two coexisting asymmetric solutions that are each other's mirror images.
• A symmetric solution (2s) also appears at the frequency X ¼ 0:48874 with the same segment structure as the asymmetric solution 2a. Solution 2s is stable below X ¼ 0:48874, but unstable in the domain X 2 ½0:48874; 0:53529 • At X ¼ 0:53529, both sticking segments become sliding ones via crossing-sliding (CS s ) and the symmetric, non-sticking stable solution 0s appears with the structure This is the same parameter where the the asymmetric solution 1a branches out.
According to these results, the abrupt appearance of asymmetric solutions is related to the fold bifurcation (SN 1a ) of solution 1a that emerges from the symmetric solution via a crossing-sliding bifurcation CS s . At the critical parameter, where the symmetric non-sticking solution disappears, the trajectory just touches the stick-slip boundaries on both the positive and the negative side. However, the appearing symmetric sticking solution is unstable, while a pair of asymmetric sticking solutions appears.
The bifurcation structure depicted in Fig. 5 remains qualitatively unchanged at larger values of the kinetic friction parameter S, but the domain of asymmetric solutions shrinks and disappears at S % 0:38. However, significant changes occur if S is decreased: the connection between the solutions 1a and 2a breakes, giving rise to multi-periodic and chaotic solutions (see Fig. 6).
Similar bifurcation structures can be observed at higher values (S 0 ¼ 0:5; . . .0:9) of the static coefficient of friction as well, but the frequency domain of asymmetric solutions extends towards higher values of X.
However, the reduction of parameter S 0 leads to the separation of the domains of the asymmetric solutions 1a and 2a, as it is shown in Fig. 7 for S 0 ¼ 0:3 and S ¼ 0:22. If S becomes larger, the two loops shrink and disappear at S ¼ S 0 ¼ 0:3.
If the static coefficient of friction is reduced further (S 0 ¼ 0:2) , the loop in the diagram corresponding to solution 1a persists at X ¼ 0:5, but solution 2a disappears. Still, multi-periodic or chaotic responses can be detected in the interval X 2 ð0:18; 0:38Þ.

Two-parameter bifurcation analysis and search for chaos
As it was mentioned earlier, the occurrence of multiperiodic or chaotic solutions is expected in certain parameter domains. To localize these domains, twoparameter bifurcation analysis was performed on the S À X parameter plane at various values of the static friction parameter S 0 . Certain bifurcation curves could be continued automatically by COCO, but it did not work in all the cases due to convergence problems. These latter curves were contructed by several one-parameter runs. The results are shown in Figs. 8, 9, 10 and 11. In addition to the bifurcation analysis with COCO, the largest Lyapunov exponents of solutions were also estimated on grids of size 100 Â 100 to 150 Â 150 on  the depicted X À S parameter planes. During each run, the evolution of two nearby orbits was traced and the variation of their phase-space distance d was monitored in time steps Dt. In the most simple case, the largest Lyapunov exponent can be estimated as However, in case of the examined frictional system, it may happen that one or the other solution is just in a sticking phase at t ¼ ði þ 1ÞDt for some index i. In this case, a later time instant t ¼ ði þ jÞDt was chosen such that both solutions slide, and the distance d iþj was taken into account considering 1 jDt ln jd iþj j jd i j in the sum. This approach provides adequate results, as it was pointed out in [41]. Moreover, a similar procedure was already tested on the same system in [30] and the results agreed well with the Lyapunov exponents calculated by a more specific method, published in [42].
The initial conditions ðx 0 ; v 0 Þ for the two solutions were ð1 þ 0:9S 0 ; 0Þ and ð1 þ 0:9S 0 þ 10 À12 ; 0Þ, respectively. The value of Lyapunov exponent was estimated in each step at most until 600 exitation periods. The simulation finished earlier if the distance d reduced below 10 À18 or the relative error between subsequently calculated estimations dropped below 10 À3 . Since we are interested in chaotic responses, only parameter domains with positive Lyapunov exponents are shown in Figs. 8, 9, 10 and 11.

Check of the results with modified models
The oscillator model examined in the previous sections is quite simple, so its validity is limited if the dynamic behaviour of engineering structures with frictional contact is to be analysed. Thus, before drawing conclusions, two similar but more complex models are also examined, to check how sensitive the obtained results are to small changes of the model.

Two-frequency excitation in the Coulomb model
In realistic scenarios, the periodic excitation signals are not pure sinusoidal, but contain higher harmonics, too. In order to describe the effect of harmonics by a single parameter, we restrict ourselves to the case where parameter h characterizes the amplitude of the higher harmonic. The remaining notations are the same as in the case of Eq. (1), and the phase-space can be partitioned similarly, exchanging cosðXtÞ by cosðXtÞ þ h cosð3XtÞ in Table 1. As Fig. 12 illustrates, positive and negative values of parameter h lead to sharper peaks or more rounded peaks in the excitation signals, respectively.
Extensive numerical experiments showed that the main qualitative features of the bifurcation diagrams remain unchanged at reasonably small values of the higher-frequency perturbation h. For example, the dependence of asymmetry W on the dynamic friction parameter S is depicted in Fig. 13 for X ¼ 0:5 and h ¼ 0:05. Comparing this diagram with Fig. 4, one can check that the same bifurcations could be detected as in the case of the pure harmonic excitation.
Certainly, the additional excitation term leads to quantitative changes in the bifurcation structure. As Fig. 14 illustrates, both the crossing-sliding and the period-doubling bifurcations tend to appear at smaller values of the kinetic friction coefficient S at X ¼ 0:5 if the amplitude of the additional term h is increased. Consequently, the parameter domain of asymmetric solutions is extended at negative values of h.
The main structure of the two-parameter bifurcation diagrams is also preserved at small absolute values of parameter h, as it is shown in Fig. 15 (h ¼ 0:05) and Fig. 16 (h ¼ À0:05). However, comparing these diagrams with Fig. 9 that depicts the scenario at h ¼ 0, one can observe that the domain of positive Lyapunov exponents is extended towards higher values of the excitation frequency X as a result of the additional, higher-frequency excitation term.
It is worth to note here that certain domains with positive Lyapunov exponent appear in areas where one or more stable asymmetric solution types exist. Since the existence of an asymmetric solution means that the mirror image of that solution exists, too, this result implies that the detected irregular solutions are chaotic transients.
A remarkable difference between the cases h ¼ 0:05 and h ¼ À0:05 is that the parameter domain of coexisting solutions 1a and 0s is significantly larger in the latter case (see the area between the bifurcation curves SN 1a and CS s in the figures). Moreover, the CS a crossing-sliding bifurcation curve, where 1-stick solutions turn to 2-stick asymmetric ones, shrinks to a small closed curve at h ¼ À0:05, with a stable (solid line in Fig. 16) and an unstable (dashed) branch. 2stick asymmetric (2a) solutions branch out from the 2stick symmetric one along the bifurcation curve SN 2s . Note, that at h ¼ 0:05 these asymmetric solutions appear above the SN 2s curve, as well, due to fold bifurcations.
Perhaps the most important effect of the additional forcing term is that the relative position of curves SN 2s , SN 1a and CS s vary such that the cases h [ 0 and h\0 correspond to higher and lower static friction coefficients S 0 at h ¼ 0, respectively (compare Figs. 15 and 16 with Figs. 10 and 8).

LuGre friction model
Certain experimentally observed phenomena like presliding, rate-dependence and hysteresis of force w.r.t. velocity cannot be described by the Coulomb model. To capture these characteristics of friction, the so-called rate and state friction models incorporate at least one internal state variable that varies dynamically in time. The friction force depends both on this state variable and on the relative velocity of contacting surfaces. One of the most popular rate and state models is the LuGre model [34] that has been applied to a wide variety of systems, including high-precision applications [35]. Moreover, this model often serves as a starting point to develop more advanced friction models.
By exchanging the Coulomb model with the LuGre model in Eq. (1), one obtains the following set of equations: where z is the internal variable and It can be shown that in case of steady-state sliding (i.e., when _ z ¼ 0), the friction force can be given as F f ¼ gðvÞsgnðvÞ þ cv. Consequently, parameters S and S 0 correspond to the kinetic and static coefficients of friction in the Coulomb model. Still, the comparison of systems (1) and (7) is certainly not straightforward due to the additional parameters of the LuGre model. In order to obtain an acceptable level of similarity between these systems, a small Stribeck velocity v S and a large internal stiffness r 0 is required. After some numerical experiments, we used a ¼ 1 and v S ¼ 0:05 to ensure the quick transition between the kinetic and static friction parameters. Since viscous damping was neglected in the previous sections, we chose c ¼ 0. The stiffness and damping associated to the internal variable were set to r 0 ¼ 1000 and r 1 ¼ 0, respectively.
Our goal was to reconstruct some of the characteristic bifurcation diagrams of the Coulomb friction model. The application of the continuation method at X ¼ 0:5 and S 0 ¼ 0:4 led to the results depicted in Fig. 17.
Comparing this diagram with the corresponding figures belonging to the Coulomb-friction oscillators (see Figs. 4 and 13), one can see that the scenario is basically the same in all the three cases: the symmetric solution loses stability at a certain value of S, giving rise to asymmetric solutions. The previously found discontinuity between the branches of the symmetric and asymmetric solutions is replaced by a meandering curve in case of the LuGre model. Unfortunately, the numerical treatment of this stiff system is rather difficult, this is why the identification of the branch point and the further continuation of branches was not successful.
Another series of continuation runs resulted in Fig. 18.
Just as before, the structure of the bifurcation diagram is similar to the corresponding figure (Fig. 5) obtained by the analysis of the Coulomb friction oscillator: the symmetric solution is exchanged by an asymmetric one in an interval near X ¼ 0:5. Now the asymmetrical branch became even more tortuous and additional loops appeared in the domain X 2 ð0:44; 0:47Þ. As a consequence, the numerical algorithm was unable to find the links between the depicted solution branches, hindering the construction of two-parameter diagrams. Note that the similar pattern in Fig. 5 is formed via non-smooth bifurcations that certainly could not occur in this case.
The calculation of the Lyapunov exponent revealed positive values at very small kinetic friction parameters S-in partial agreement with the previous results. However, the calculated exponents were scattered seemingly randomly close to zero in large domains. The relative inaccuracy of the numerical method stems from the fact that while we were able to use analytical formulas between stops in case of the Coulomb friction oscillators, the numerical simulation of the LuGre model required the application of a stiff differential equation solver. The unavoidable numerical damping led to the disappearance of weakly chaotic solutions.

Conclusions
Themaingoalofthepresentcontributionwastoestablish the connection between asymmetric and chaotic solutions of simple dry-friction oscillators: one of them is excited purely harmonically, while an additional excitation term h cosð3XtÞ appears in the other case. In order to find the typical solution types, broad parameter domains were explored by continuation method. To facilitate the application of the method, the phasespaces of the examined systems were divided into seven domains that are separated by three switching surfaces such that the crossing of a switching surface unambiguously determines the state of the solution in the next segment. This reformulation of the problem makes easier to treat the case of different static and kinetic coefficients of friction. The application of the continuation method was supplemented by the determination of the Lyapunov exponent.
The results show that chaotic solutions are typically enclosed by the bifurcation curves SN 2s (where the 2stick symmetric solution loses stability, giving rise to 2-stick asymmetric solutions) and SN 1a (where the 1stick asymmetric solution loses stability). Thus, the chaotic responses are bifurcating from asymmetric solutions, as it was expected in [30]. Since asymmetric solutions were detected at X ¼ 0:5 for the case S ¼ S 0 and h ¼ 0, it was assumed in [30] that this is the parameter where it is worth to search for chaotic solutions. However, our results clearly show that the choice X ¼ 0:5 was rather unlucky and more pronounced chaotic regimes can be found at lower values of the frequency X, or at nonzero values of h.
It is remarkable that positive Lyapunov exponents were found in certain parameter domains where stable asymmetric solutions exist. The explanation of this fact is that each asymmetric cycle has a counterpart with opposite displacements and different basin of attraction. Thus, certain solutions can exhibit chaotic transients in the neighbourhood of the corresponding basin boundary. As the parameters are changed, the asymmetry-related transient chaotic solutions may turn to chaotic ones.
According to the results, quite simple connections can be found between the friction coefficients and the location of the chaotic domain: • Irrespectively of the static coefficient, the domain of asymmetric solutions-and the enclosed chaotic domain-typically shrinks if the kinetic coefficient of friction S is increased. • The increase of the static coefficient of friction implies that the domain of atypical solutions is extended towards larger frequencies.
It can also be seen in the figures that S\0:2 is needed for chaotic responses, at almost all values of the static friction parameter S 0 . However, a narrow parameter regime was found at S 0 ¼ 0:3 and h ¼ 0 where the Lyapunov exponent is positive even at S % 0:265. The corresponding ratio S 0 =S % 1:13 is quite realistic, thus, the occurrence of chaotic responses cannot be excluded in similar models of real mechanical systems. It is worth to note here that the case S 0 ¼ 0:3 seems to be an exception, since the further reduction of S 0 does not lead to the increase of the chaotic domain at h ¼ 0: only small islands of positive Lyapunov exponents were detected at S 0 ¼ 0:2. Because the considered Coulomb friction model is rather simple, a third oscillator model was also analysed where the friction force was calculated according to the LuGre model. Two characteristic bifurcation diagrams of the harmonically excited Coulomb friction oscillator were qualitiatively reproduced with this model, as well. This result confirms our opinion that the traces of certain theoretically derived solution types can be observed in real systems of contacting bodies as well, as it was shown by the experiments in [31]. Since quite simple qualitative rules could be established between the frequency domain of irregular solutions and the coefficients of friction, these rules may be applicable to real mechanical systems, too.