## Abstract

In a recent paper we examined a model of an arch bridge with viscous damping subjected to a sinusoidally varying central load. We showed how this yields a useful archetypal oscillator which can be used to study the transition from smooth to discontinuous dynamics as a parameter, *α*, tends to zero. Decreasing this smoothness parameter (a non-dimensional measure of the span of the arch) changes the smooth load–deflection curve associated with snap-buckling into a discontinuous sawtooth. The smooth snap-buckling curve is not amenable to closed-form theoretical analysis, so we here introduce a piecewise linearization that correctly fits the sawtooth in the limit at *α*=0. Using a Hamiltonian formulation of this linearization, we derive an analytical expression for the unperturbed homoclinic orbit, and make a Melnikov analysis to detect the homoclinic tangling under the perturbation of damping and driving. Finally, a semi-analytical method is used to examine the full nonlinear dynamics of the perturbed piecewise linear system. A chaotic attractor located at *α*=0.2 compares extremely well with that exhibited by the original arch model: the topological structures are the same, and Lyapunov exponents (and dimensions) are in good agreement.

## 1. Introduction

A simple shallow arch model, proposed by Thompson & Hunt (1973), comprises a linked pair of inclined springs (capable of resisting both tension and compression) which are pinned to rigid supports. This model is widely used in physics and engineering to illustrate and model the snap-through buckling phenomenon (e.g. Pippard 1990; Savi & Pacheco 2003). Initial studies were given via energy methods for the pitch fork bifurcation and the nonlinear response. In our recent paper (Cao *et al.* 2006) we have extended this model to define a new oscillator named the SD oscillator, which allows one to study the transition from smooth to discontinuous systems. We have also conducted preliminary studies of its dynamics: a more comprehensive treatment will be given in a future paper.

An experimental study carried out by Tufillaro *et al.* (1992) shows the complex nonlinear dynamical behaviour of such a system. Tufillaro also modelled a string oscillating in its fundamental mode, consisting of a single mass fastened to the central axis by a pair of linearly elastic springs (e.g. Tufillaro 1989, 1990).

In our present arch (as in the Tufillaro model), it is the large changes in geometry that provide a highly nonlinear resultant force on the central oscillating mass, despite the fact that the springs (strings) themselves are assumed to have a linear load–deflection relationship (in both tension and compression). It is this nonlinearity that generates complex dynamical behaviour with period doubling and chaos, etc., when the arch is subjected to periodic excitation. The equation of motion can indeed be simplified to a Duffing equation (see Nayfeh & Mook 1979; Wiggins 1990; Thompson & Stewart 2002), which is also known to exhibit a wealth of chaotic phenomena (e.g. Pikovsky 1992; Anna & Vered 1997; Aguirre & Sanjuán 2002).

Our present work introduces a piecewise linearization, which reproduces well the full original behaviour of the arch, while allowing an analytical study of the nonlinear responses (see appendix A). This approach enables us to investigate the nonlinear dynamics theoretically for both the unperturbed (undamped and undriven) and perturbed systems.

This paper is organized as follows. In §2 the equations of motion are formulated for both the smooth and discontinuous oscillators, and the Hamiltonian functions are derived: these are then used to classify the forms of the trajectories. We see that the equilibrium point (0, 0) in the limit (discontinuous) case is a saddle-like point which creates a homoclinic-like orbit composed of this point and two circles. In §3, the piecewise linearized system is derived to approximate the original system. Analytical expressions are developed for the homoclinic orbit and families of periodic orbits surrounding the separate centre points inside and outside the homoclinic orbit. This approach is valid for both the smooth and the limit discontinuous systems. In §4, the Melnikov method is applied to determine the distance between the stable and unstable manifolds (see Melnikov 1963). The comparison between the dynamics of the original oscillator and its piecewise linear approximation shows a good degree of correlation.

## 2. The SD oscillator and preliminary dynamic analysis

Consider a simple elastic arch rigidly supported at each end (as shown in figure 1*a*). Assuming that one is only interested in the first (symmetric) vibrational mode, this system can be simplified to a nonlinear oscillator with a lumped mass, *m*, linked by a pair of linearly elastic springs of stiffness *k* and equilibrium length *L* which are pinned to rigid supports (see figure 1*b*).

The equation of motion can be written in the following form:(2.1)where *X* is the mass displacement and *l* is the half distance between the rigid supports. System (2.1) is made dimensionless by letting , and ,(2.2)This system can be smooth or discontinuous depending on the value of the smoothness parameter *α* (as shown in figure 2). Note, though, that the case of *α*=0 has become physically unrealistic because the half-span, *l*, has become zero: we discuss this limit more fully in a forthcoming paper.

When *α*=0, equation (2.2) can be written in the following form:(2.3)Equation (2.3) is discontinuous and it is derived from the smooth system (2.2) as *α*→0.

Now suppose system (2.2) has viscous damping and an external harmonic excitation of amplitude *F*_{0} and frequency *Ω*. This leads to the forced dissipative oscillator of the following form:(2.4)Again equation (2.4) can be written in a dimensionless form by letting , and ,(2.5)where . The equation (2.5) has been named as the SD oscillator in Cao *et al.* (2006).

The equilibria of the system (2.5) with zero damping and no external force can be easily obtained and they areFor *α*>0, at the equilibrium point (0, 0) eigenvalues of the system can be obtained by examining the Jacobian matrix of the unperturbed system and they are , which means that this equilibrium is a saddle point.

If *α*=0, the singularity (0, 0) of the discontinuous system is the limit of the saddle point of the smooth system as *α*→0. The hyperbolicity is broken as the eigenvalues do not have finite limit as *α*→0. Therefore, the singularity point (0, 0) for the discontinuous system will be called a saddle-like singularity (SLS).

In the same way, the eigenvalues of the other two equilibrium points () can be obtained as for both the smooth and discontinuous cases.

To obtain the Hamiltonian function, equation (2.5) has to be rewritten in the form of two first-order differential equations. Then introducing the notation *y* for the velocity *x*′, one can get the Hamiltonian function as follows:(2.6)All types of trajectories are listed in table 1 for different values of the Hamiltonian function , while the corresponding phase portraits are shown in figure 3.

With the help of the Hamiltonian (2.6), the trajectories can be classified and analysed. For both continuous and discontinuous cases, the phase portraits of systems (2.2) and (2.3) are plotted for different values of the Hamiltonian . For instance, for the smooth nonlinearity, *α*=0.2, the dynamic behaviour is similar to that of the double-well Duffing oscillator (Guckenheimer & Holmes 1983; as shown in figure 3*a*). For the discontinuous case, *α*=0, the behaviour is shown in figure 3*b*: the orbits for *E*>0 comprise two large segments of circles with their centres located at (−1, 0) and (1, 0) connected at *x*=0. The case *E*<0 is represented by two families of circles.

If *α*=0, the unperturbed system (2.2) is reduced to system (2.3), whose Jacobian at (0, 0) does not exist. This case is the most interesting as the discontinuous system (2.3) exhibits two circles centred at (±1, 0), connected at the singularity point (0, 0), and forming a special singular homoclinic-like orbit. This point (0, 0) is a new type of isolated singularity which has neither eigenvalue nor eigenvector. The two circles excluding the point (0, 0) are not the manifolds of the singularity, but the flow along these circles approaches the point (0, 0) at a rate of as *x*→0 and it will be trapped by the singularity.

## 3. Unperturbed piecewise linear system

The Hamiltonian (2.6) is not readily analysed (see appendix A), so to allow further theoretical investigation of trajectories (particularly under harmonic excitation), it is useful to make a piecewise linear approximation to the load–deflection curve.

The nonlinear restoring force (shown in figure 4*a* by the dashed curve) can be expanded into a Taylor series in the neighbourhood of *x*=0 as(3.1)where *P*_{n}(*x*) is the Legendre polynomial defined byThe first few values areIt is worth noticing that truncating the Taylor expansion (equation (3.1)) after the cubic term and substituting it into equation (2.5) leads to a Duffing oscillator with the restoring force, , which is marked as a dotted curve in figure 4*a*. This truncation is valid for *α*>0; it fails for *α*=0.

To make our piecewise linear approximation of the function *F*(*x*), we introduce the trilinear function shown in figure 4*a* by solid lines. The middle section of this trilinear curve, i.e. , is obtained by taking the first term in the Taylor expansion equation (3.1). The other two sections of the trilinear function for are straight lines which are tangential to the function *F*(*x*) at the equilibrium points . To fit these lines, the derivative of the function *F*(*x*) is calculated aswhich is evaluated at as(3.2)Thus, we obtain the trilinear function for *F*(*x*) as(3.3)where , and *x*_{0} is the slope changing point for the trilinear function obtained asThis piecewise linearized force, , is shown by the solid curve in figure 4*a*. The corresponding potential energies are shown in figure 4*b*: the dashed, dotted and the solid curves are for the nonlinear, Duffing and piecewise linear approaches, respectively. In figure 4*a*,*b* we can see a satisfactory agreement between the piecewise linear approach and the original nonlinear problem. The cubic (Duffing) approximation fails badly in both the restoring force and potential energy plots.

The fact that *x*_{0}→0 for *α* tending to zero indicates that the piecewise linear approximation gives a precisely correct fit for the case *α*=0. This precise fit (illustrated in figure 2) means that we can use the piecewise linear model to study the change in dynamics as we decrease *α* and sweep from the continuous to the discontinuous system.

### (a) Hamiltonian function

The piecewise linearized system, obtained for the unperturbed system using the restoring force approximation equation (equation 3.3), can be written as a system of two first-order differential equations as(3.4)It should be noted that *α*>0 is being assumed in the following discussion. If *α*=0, it implies *x*_{0}=0. It is reasonable to write the first part of the second equation of system (3.4) as *y*′=0 for the case of *α*=0 when *x*=0.

The Hamiltonian function of equation (3.4) can be obtained and written in the following form:(3.5)The phase portrait plotted via the Hamiltonian function (3.5) for the piecewise linear system (3.4) for *α*=0.2 is shown in figure 5*a*.

### (b) Homoclinic orbit

The homoclinic orbit of the system (3.4) obtained by taking intersects the *x*-axis at *x*_{c}, so that , i.e.This homoclinic orbit can be divided by vertical lines into two elliptical segments, marked as and , and four line segments, marked as , , and , meeting at the point (0, 0) as *τ*→+∞ (as shown in figure 5*b*).

The elliptical segments of the homoclinic orbit for have been obtained in analytical form as a parametric function of the time variable *τ*,(3.6)whereThe segments of the homoclinic orbit in the regime have been found similarly as(3.7)The limiting tangents to these line segments are(3.8)which tend to infinity as *α*→0.

In the limit case of *α*=0, it can be seen that *R*=1 from (3.3) and *ω*_{2}=1 from (3.3). Equations (3.6) and (3.7) lead to(3.9)where This is the solution of the special homoclinic-like orbit of the discontinuous system, which is shown in figure 3*b* marked as *E*=0. This analytical solution can also be obtained by considering both (3.6) and (3.7) for *α*→0. This trajectory is the limit case of the homoclinic orbit of the piecewise linear system (3.4).

### (c) Periodic orbits

In this subsection, periodic solutions are presented for the piecewise linear system for *α*>0. As for *α*=0, the corresponding solutions can be obtained by considering . At the equilibrium points , the Hamiltonian constant 2*E* is equal to . If , the two families of periodic orbits exist inside the homoclinic orbit surrounding the centre equilibrium points. The parts of the orbits determined for can be obtained and written as(3.10)whereIn the same way, the parts of the orbits have been determined for ,(3.11)For *E*>0, a family of periodic orbits exist outside the homoclinic orbit surrounding all the singularities. For , these parts of the orbits have been obtained as(3.12)whereFurthermore, for , these orbits can be obtained and written as(3.13)

## 4. Perturbed system

In this section we focus our attention on the piecewise linear system obtained above, but adding driving and damping to give us what we call the perturbed system. The Melnikov method will be used to detect the homoclinic tangency (the Melnikov boundary) beyond which (possibly localized) chaos will be generated. Subsequently, a semi-analytical method will be used to investigate the nonlinear dynamics of the perturbed system.

### (a) Melnikov analysis

We consider the perturbation of our piecewise linearized system (3.4) by viscous damping and an external harmonic forcing in the form(4.1)where is the piecewise linearized force given in (3.3).

Then we introduce following notation in system (4.1):For the homoclinic orbit, we get(4.2)where the operator is defined as(4.3)for any and .

The corresponding Melnikov function of equation (4.1) (see Melnikov 1963; Guckenheimer & Holmes 1983; Chen & Leung 1998), for instance, is given by(4.4)The calculation for this integral leads to the following Melnikov function (see details in appendix B):(4.5)whereIt can be seen that(4.6)has simple zero for *τ*_{1} if and only if the following inequality holds:(4.7)The Melnikovian detected chaotic boundary for system (4.1) obtained for *α*=0.2 is shown in figure 6*a*, and for *α*=0 is shown in figure 6*b*, for the different values of parameter and 0.060. In the area above the detected boundaries, chaotic motions will be generated, i.e. in the area of

### (b) Semi-analytical analysis

In this section a semi-analytical approach (see Pavlovskaia & Wiercigroch 2004) is presented for the solutions of equation (4.1) to investigate the nonlinear behaviour of the piecewise linearized system in the chaotic regions detected by the Melnikov curves computed in §4*a*. The Lyapunov exponents and the corresponding dimension of the attractors are calculated as in Takens (1981), Grassberger & Procaccia (1983), Lai & Lerner (1998), Stefanski & Kapitaniak (2000).

The phase space of the system (4.1) is divided by the boundaries into three subspaces , and , which correspond to three different regimes of the system. The equations of motion are linear for each regime; therefore, the global solution can be constructed by joining the local explicit solutions for each regime at the points of discontinuity when the trajectory crosses the boundaries . The set of initial values defines in which regime the system will operate.

In the first subspace, where , the solution for the initial conditions is written as(4.8)where , andwithIn the second and third subspaces, where and , the solutions for the initial conditions are(4.9)whereThe constants and can be determined from initial conditions asWhen the conditions corresponding to the current regime fail, the system starts to operate in the next regime, and the final displacements and velocity for the preceding regime define the initial conditions for the next one. To construct the global solution, extra care is required in determining values of the crossing time *τ*_{i} since the response can be very sensitive to any inaccuracy on the discontinuity boundaries. This is achieved by solving numerically the appropriate nonlinear algebraic equations (describing the boundary crossings) with the required precision.

This semi-analytical method has been used to obtain the bifurcation diagrams shown in figures 7*a* and 9*a*, and the chaotic attractors shown in figures 8*a* and 9*b*. Figure 7 demonstrates the comparison between the bifurcation diagrams calculated for the original system (2.5) and its piecewise linear approximation (4.1) and shows corresponding views in figure 7*a*,*b* for *ω*=1.05, *α*=0.2 and *ξ*=0.015. As can be seen in figure 7*b*, no chaotic motion is observed below the boundary predicted by the Melnikov method which is marked by the dashed line in figure 7*a*.

The structure of a typical chaotic attractor is demonstrated in figure 8. Attractors displayed in figure 8 are calculated for the original smooth system (2.5) (figure 8*a*) and the piecewise linear approximation (4.1) (figure 8*b*) for the same parameters *f*_{0}=0.83, *α*=0.2, *ω*=1.05 and *ξ*=0.015. As can be seen in figure 8, a good degree of correlation is obtained both in the attractors' topological structure and in the values of the Lyapunov exponents and the fractal dimensions (the latter are shown in the figure legend). A pertinent historical observation is that the topology of the present chaotic attractor has intriguing similarities to the early attractor identified by Thompson & Ghaffari (1983), and also Thompson & Stewart (2002).

In the limit case *α*=0, the original system is piecewise linear and there is no difference between the original system and its approximation. Thus, the same semi-analytical method can be used to obtain the system responses and also for constructing the bifurcation diagrams. A bifurcation diagram under controlled variation of the amplitude of the external excitation is shown in figure 9*a* for *α*=0, *ω*=1.05 and *ξ*=0.015. Again the boundary predicted by the Melnikov method is marked by the dashed line, and there is no chaotic motion below it. As can be seen in figure 9*b*, where a chaotic attractor is presented for the case *α*=0 at *f*_{0}=0.83, the structure of the attractor differs significantly from those obtained at *α*≠0.

## 5. Conclusions

An archetypical oscillator introduced in Cao *et al.* (2006) in which the strength and smoothness of the nonlinearity can be controlled by the parameter *α* has been studied. A piecewise linear approximation of this oscillator has been developed and investigated. This piecewise linear system becomes coincident with the original oscillator at the limit of *α*=0 (where the smooth restoring force has become sawtoothed) and has the advantage of being amenable to closed-form analysis for all *α*. Using a Hamiltonian formulation, we have made a closed form Melnikov analysis to estimate the homoclinic tangency at which (at least localized) chaotic motions will occur. Chaotic attractors are presented for both the original oscillator and the piecewise linear approximation at *α*=0.2. They compare extremely well: the topological structures are the same, and Lyapunov exponents (and dimensions) are in good agreement.

## Acknowledgments

The first two authors acknowledge the financial support from EPSRC under grant no. GR/R85556, and the Scottish Enterprise.

## Footnotes

One contribution of 14 to a Theme Issue ‘Experimental chaos II’.

- © 2007 The Royal Society