## Abstract

Real-time dynamic substructuring is a powerful testing method, which brings together analytical, numerical and experimental tools for the study of complex structures. It consists of replacing one part of the structure with a numerical model, which is connected to the remainder of the physical structure (the substructure) by a transfer system. In order to provide reliable results, this hybrid system must remain stable during the whole test. A primary mechanism for destabilization of these type of systems is the delays which are naturally present in the transfer system. In this paper, we apply the dynamic substructuring technique to a nonlinear system consisting of a pendulum attached to a mechanical oscillator. The oscillator is modelled numerically and the transfer system is an actuator. The system dynamics is governed by two coupled second-order neutral delay differential equations. We carry out local and global stability analyses of the system and identify the delay dependent stability boundaries for this type of system. We then perform a series of hybrid experimental tests for a pendulum–oscillator system. The results give excellent qualitative and quantitative agreement when compared to the analytical stability results.

## 1. Introduction

There are several methods for predicting the dynamic response of structures to external forces. For simple linear systems, analytical methods can usually be applied. More often, for large, complex and/or nonlinear systems some type of numerical method will be required. In both cases, the approach is one where a mathematical model is used to approximate the behaviour of a physical system. For complex systems, mathematical models may contain significant uncertainty and it may be more appropriate to build prototype systems and to perform tests in the laboratory (Williams & Blakeborough 2001). In cases when the complex system is very large—many civil and aerospace engineering structures fall into this category—such an approach may be impractical and quite often prohibitively expensive. In these cases, it is desirable to use testing techniques which can include the benefits of both theoretical and experimental methods.

One such approach is real-time dynamic substructuring (for recent literature see, for instance, Blakeborough *et al*. 2001; Horiuchi & Konno 2001; Nakashima 2001; Bonelli & Bursi 2004; Pinto *et al*. 2004; Wallace *et al*. 2005*a*). This method consists of combining numerical and physical substructures using a set of transfer systems—usually a set of actuators. Delay naturally arises in this type of test system because of the non-instantaneous nature of the transfer system. In fact, there can be a number of different delays from several different sources, such as data acquisition, computation, digital signal processing and the actuator itself, which all contribute to produce an overall delay. Despite this seeming complexity, modelling the overall delay as a single fixed quantity, *τ*, has been shown to give very good agreement when compared to hybrid test results (Horiuchi *et al*. 1999; Darby *et al*. 2001, 2002; Wallace *et al*. 2005*a*). This is due to the fact that the overall delay is nearly always dominated by the actuator delay component, which for a small range of frequencies remains at a similar value. When significant variations in frequency exist, adaptive techniques can be used to maintain system stability (e.g. Wallace *et al*. 2005*b*).

The assumption of a fixed delay allows the system to be modelled using delay differential equations (DDEs). These could be partial or ordinary, depending on the complexity of the system. DDEs depend not only on the current state of the system, but also on the history of the system over some previous time interval (the delay). Hence, the initial state space and the solution space of the delayed dynamical system are infinite dimensional. The theory of these type of equations lies in the area of functional differential equations (a detailed discussed is given by Wu (1996)). Some methods and results with applications to engineering problems can be found in Xu & Chung (2003) (see also Stépan (1989); Kuang (1993); Diekmann *et al*. (1995), and references therein).

The application of DDE modelling techniques to substructuring has been considered by Wallace *et al*. (2005*a*,*b*), who considered modelling linear mass–spring–damper (MSD) systems. For these system, the authors demonstrated the dependence of the stability of the system on delay and linked their work with previous energy analysis approaches to this problem (e.g. Horiuchi *et al*. 1999).

In this paper, we consider an example of a system which consists of a pendulum coupled to a mechanical oscillator or MSD. The MSD is linear, and therefore taken to be the numerical model while the pendulum, the physical substructure, was constructed in the laboratory. This oscillator–pendulum system has been studied by Gonzalez-Buelga *et al*. (submitted), who showed that the autoparametric resonance behaviour, known to occur in this type of system, can be modelled using a hybrid numerical–experimental simulation. The equations of motion governing the motion of the oscillator–pendulum system have the added complexity that they are *neutral* DDEs, unlike those studied by Wallace *et al*. (2005*a*,*b*).

The specific focus of this paper is to provide a rigorous mathematical treatment of linear stability of the pendulum–oscillator system, when modelled as a hybrid system with a fixed delay. In particular, we are interested in the stability boundaries at which a hybrid test loses stability via a Hopf bifurcation. We obtain the analytical expressions for stability boundary as a function of parameters of the system. These analytical expressions are then verified by comparison with a series of hybrid experimental tests, where the delay is artificially increased until the stability boundary is located. The results show excellent qualitative and quantitative agreement between analytical predictions and experimental data.

Further to this, using the results from the stability analysis, we demonstrate how for a particular value of mass ratio, critical delay values and instability frequencies can be predicted using the DDE model. The ability to predict these quantities for a particular hybrid testing configuration is significant from a practical point of view. The critical delay values give information on when a test will be destabilized and can be used to design an appropriate delay compensation scheme (e.g. Wallace *et al*. 2005*a*). More fundamentally, for systems of limited actuator capacity, critical delay values can be used to determine the viability of a hybrid test. Some of the possible applications which are directly relevant to the coupled pendulum–MSD systems include the modelling of human leg motion (e.g. Lafortune & Lake 1995; Coveney *et al*. 2001) and cable-stay bridge vibration Gatulli *et al*. (2005).

The paper is organized as follows: in §2, we introduce a system of DDEs describing the coupled pendulum–MSD system. Section 3 is devoted to the stability analysis of the neutral DDE. In §4, stability of periodic solutions is investigated by means of multiple scale analysis for the case of near resonant soft excitation. Section 5 contains numerical simulations, which confirm our analytical findings. The linearized stability analysis of the coupled nonlinear system of DDEs is presented in §6. Section 7 contains a comparison between experimental and analytical results. In §8, the effects of viscous damping are studied. The paper concludes with a summary of results in §9.

## 2. The equations of motion

We consider a mechanical system which consists of a mass *M* mounted on a linear spring, to which a pendulum of mass *m* is attached via a hinged massless rod of length *ℓ*. The angular deflection of the pendulum from the downward position is denoted by *θ*. We assume linear viscous damping *C* acting on the mass *M*. The system is shown schematically in figure 1*a*. The equations of motion for this system are(2.1)whereis the external excitation force applied in the *y*-direction, *C* and *K* are the damping and stiffness coefficients, respectively, and a dot indicates the derivative with respect to time *t*. We will refer to *y*(*t*), and as the position, velocity and acceleration of the MSD at time *t*. There are three important frequencies associated with this system: the natural frequency of the pendulum *ω*_{0}, where , the natural frequency of the mass–spring–damper *ω*_{1}, where and the frequency *Ω* of the external excitation force. The situation where the natural frequencies of the MSD and the pendulum are in 2 : 1 resonance has been widely studied in the literature (see Tondl *et al*. (2000) and references therein).

Since the transfer system produces an unavoidable delay in the system, the feedback force will be delayed. The time delay is denoted by *τ* and it is assumed to be constant and non-negative. To account for the delay in the displacement signal, the force in the system (2.1) has to be described by the delayed state of the numerical model of the MSD. Therefore, in the right-hand side of the first equation in (2.1) we replace *y*(*t*) and *θ*(*t*) by their delayed counterparts: *y*(*t*−*τ*) and *θ*(*t*−*τ*), respectively. The same has to be done for *all* terms in the second equation of the system (2.1). By inserting these expressions in system (2.1), we obtain the following:(2.2)A key objective of hybrid testing is to get (2.2), the equations of the hybrid system, to replicate the dynamics of (2.1) as closely as possible. This gives rise to two interlinked issues: stability and accuracy. The issue of accuracy is still an open problem. For recent discussions and attempts to quantify accuracy see Mosqueda (2003) and Wallace *et al*. (2005*b*) (and references therein). In this paper, we will concentrate our attention on the stability of equation (2.2). In particular, we will be looking at the stability of the so-called trivial solution in which the MSD moves (so *y*(*t*)≠0 in general) but the pendulum does not oscillate (so *θ*(*t*)=0 in general).

## 3. Stability analysis of the neutral DDE

First, we consider the case when the angle *θ* is small (*θ*≪1). In this case, the system (2.2) decouples and we concentrate on the first equation, which describes the vertical motion of the pendulum–MSD system. In the absence of external forcing (*k*=0), linearization of this equation around the trivial solution leads to the equation(3.1)Since the delayed system depends on the acceleration of the state variable, equation (3.1) is of a neutral type. We introduce the following non-dimensional variables:Under this rescaling and dropping the hats, equation (3.1) becomes(3.2)where dot means differentiation with respect to *t*. This equation has one trivial steady state . The corresponding characteristic equation is(3.3)When *τ=0*, one obtains , so the steady state *z*=0 is locally asymptotically stable. In the case |*p*|>1, this steady state is always unstable for any positive delay *τ*. (This already provides us with an important consideration when constructing the hybrid system, namely that for stability we must always have *M*>*m*.) Therefore, we assume |*p*|<1 in all our further analysis. The purely imaginary eigenvalues occur when *λ*=±i*ϖ* for *ϖ*≠0, so from (3.3),Separating into the real and imaginary parts, we arrive at the following:(3.4)Upon squaring and adding these equations, we obtainSolving for *ϖ* gives(3.5)We observe that there can be either zero, one (repeated) or two real positive roots, depending on the values of *p* and *ζ*. This relation provides us with the explicit expressions for the stability boundaries, i.e. a family of solutions for the delay time *τ* has the form(3.6)where *n*=0, 1, 2, … and Arctan corresponds to the principal value of arctan.

*The solution**of the system**(3.2)**is locally asymptotically stable for**and* |*p*|<1 *independent of the delay time τ*>0.

The proof of this lemma immediately follows from the expression (3.5) for *ϖ*_{±}. In terms of the original parameters, the conditions of lemma 3.1 state that *C*^{2}>2*MK*, *M*>*m*. Therefore, if the damping coefficient, *C*, is high enough, unconditional stability is guaranteed.

Next, consider the case , which is depicted in figure 2. The shaded area shows the stability region. There is a value of *p*, defined below, for which the system is stable for all *τ*. However, above this *p*-value, the stability of the trivial steady state strongly depends on the value of the delay *τ*. As *τ* increases from zero, the steady state undergoes a Hopf bifurcation, as a pair of complex conjugate roots of the characteristic equation crosses the imaginary axis. Under certain non-degeneracy conditions, this implies the birth of periodic solutions. We summarize our findings in lemma 3.2.

*Assume*. *The trivial solution of system**(3.2)**is locally asymptotically stable in the region*(3.7)*for all positive delay times τ. In the region**the trivial solution is locally asymptotically stable for values of delay satisfying**where n*=1, 2, ….

First, we note that when the condition (3.7) holds, the square rootin the definition of (3.5) is purely imaginary. Thus, the eigenvalues of the linearization matrix never cross the imaginary axis, and consequently, the trivial steady state is stable for all delay times *τ*≥0 in this region.

Next, we consider the situation when . In this case, the characteristic polynomial (3.3) has two imaginary solutions *λ*_{±}=i*ϖ*_{±} with *ϖ*_{+}>*ϖ*_{−}>0 defined in (3.5). In order to determine stability as *τ* varies, we need to find the sign of the derivative of Re(*λ*) at the points where *λ*(*τ*) is purely imaginary. The technique we shall use has been widely applied to various characteristic equations (e.g. Kuang 1993). Considering the function *λ*=*λ*(*τ*) in equation (3.3) and differentiating this equation with respect to *τ*, one obtainsFrom the last expression it follows thatand from (3.3) we haveHence,Substituting into the last expression, it is clear that the sign is positive for and negative for . This means that as *τ* increases and takes a value corresponding to *ϖ*_{+}, *λ* crosses from the left to the right-hand half of the complex plane. This implies the loss of stability of the trivial solution. Similarly, for *τ* corresponding to *ϖ*_{−}, the crossing is from right to left, and stability is regained. The whole process is graphically illustrated in figure 2. The proof of lemma is complete. ▪

From lemma 3.2, it follows that as *τ* increases, the trivial equilibrium undergoes successive changes in its stability. However, it is worth noting that for fixed parameter values there are only a *finite* number of stability switches, and eventually this equilibrium becomes unstable. We can now find the number of these stability switches and the maximal delay beyond which the stability will never be recovered. The problem of constructing sequences of stability switches was addressed for first-order neutral equations by Wei & Ruan (2002).

Let us introduce a sequencewhere *p*_{n}, *n*=1, …, solve the equation(3.8)Here, we have used the notation

*If p*_{1}≤*p*<1*, there is one stability switch at τ*_{+}(0). *The trivial equilibrium is stable for* 0≤*τ*≤(0) *and unstable for τ*>*τ*_{+}(0). *If p*_{k+1}≤*p*<*p*_{k}*for k*=1, …, *then there are exactly* (2*k*+1) *stability switches, and the trivial equilibrium is stable for**and unstable for*

For a given value of *p*, theorem 3.3 gives the regions of linear stability in terms of delay time *τ*. Unless *p* is smaller than the lower stability boundary , there always exists a sufficiently large delay *τ* after which the system will be unstable. The stable boundaries of the origin can be located by the branches *τ*_{+} and *τ*_{−} which satisfy *τ*_{+}(*n*)>*τ*_{−}(*n*) for *n*=1, 2, …. For example, when *τ* crosses *τ*_{−}(1), the unstable origin becomes stable. When *τ* is increased further to cross *τ*_{+}(1), the stability is lost again as a pair of eigenvalues cross the imaginary axis into the right half plane. Similarly, other stability regions can be found for *τ*_{−}(*n*)<*τ*<*τ*_{+}(*n*). Each time we cross the stability boundary into an unstable region, the delayed action of the pendulum on the MSD leads to a destabilization of numerical model.

In the points (*p*, *τ*)=(*p*_{n}, *τ*_{n}), the system undergoes a codimension-two Hopf bifurcation. There is a pair of complex conjugate eigenvalues crossing the imaginary axis from left to right, and there is another pair crossing from right to left. Therefore, at the above points, the system has two frequencies simultaneously present. The first three (*p*_{n}, *τ*_{n}) points are indicated in figure 2. Possible resonances in this case were studied by Campbell (1997).

## 4. Hopf bifurcation

In this section, the Hopf bifurcation will be investigated. Linearizing equation (2.2) with *k*≠0, and rescaling using and , we obtain, after dropping the hats:(4.1)which is simply equation (3.2) with the forcing included. In order to study the criticality of the bifurcation, we will use a multiple scales method. Taking a delay time , where *τ*_{c} is the critical delay time, we let *T*=*ϵt* be the new time scale and rescale equation (4.1) with , . Physically this implies that the excitation is soft and a resonant case will be considered. In particular, we set , where is a detuning parameter. Before embarking on series expansion analysis, it is worth noting that since we are crossing the stability boundary, the following relation holds:Looking for solutions of equation (4.1) in the formwe obtain at order where . The solution of this equation is(4.2)where *A*(*T*) is an as yet undetermined amplitude of oscillation. At order we have(4.3)with . Substituting *z*_{0} from (4.2) in the last expression and using complex notation, giveswhere and c.c. denotes the complex conjugate. To avoid secular terms, we set the bracket on the right-hand side equal to zero. This gives us an equation for the amplitude *A* in the form(4.4)This equation can be transformed intoWriting *A*=*u*+i*v* and separating real and imaginary parts, we obtain the following linear system of differential equations:(4.5)where we have introduced the notationSystem (4.5) can be solved to yield(4.6)Thus, one can observe that the amplitude *A*_{T} of the bifurcating periodic solution is proportional to the amplitude *k* of the forcing, and is also modulated on a longer time scale as shown by the above expressions. If *k*=0 (and consequently *σ*=0), then depending on the sign of *τ*_{1}, the *u*(*T*) and *v*(*T*) exponentially tend to either 0 or ±∞, thus indicating that the periodic solution is unstable. However, if *k*≠0, then there is a periodic solutionwhich is stable if(4.7)and unstable otherwise.

## 5. Numerical simulations

In this section, we perform a numerical simulation of the system (4.1). The equation was discretized using an explicit finite difference scheme. The damping term was approximated by central differences to improve numerical stability.

For all simulations, we set *p*=0.2. In terms of the original parameters of the system this means that the mass of the pendulum is five times smaller than the mass of the MSD. At the same time, the scaled damping coefficient *ζ* is taken to be *ζ*=0.05. It readily follows that , and one can expect a succession of stability switches for increasing values of *τ*.

The numerical timestep for the simulations is kept at Δ*t*=0.01. To illustrate the behaviour of solutions in stable and unstable regions, as discussed in §3, we will move along the *τ*-axis and observe changes which occur as the stability boundaries are crossed.

We begin our numerical analysis by considering first the case when there is no external force, i.e. *k*=0 in equation (4.1). In this situation, results from §3 indicate that the trivial steady state undergoes stability switches. Figure 3 shows temporal dynamics of the non-dimensionalized displacement *z*(*t*) governed by equation (3.2) for several values of time delay *τ*. As *τ* is increased from zero (*τ*=1), the solution remains asymptotically stable and quickly decays to the trivial steady state (figure 3*a*). Since oscillatory decay is observed, this indicates that the eigenvalues of characteristic equation have non-zero imaginary part. Figure 3*b* shows the situation when the first switch from stability to instability has occurred (*τ*=5). The characteristic eigenvalues have crossed the imaginary axis, and hence the trivial state has become unstable. Stability is regained in figure 3*c* for *τ*=7.5. It can be observed in this picture that there exist different frequencies of oscillations. When *τ* is larger still, stability is lost again (figure 3*d*, *τ*=11).

As the system recovers its stability as shown in figure 4*a*, one can clearly see the appearance of *beats*. Even though the eigenvalues of the characteristic equation still have negative real part (as can be seen in the way the oscillations decay), it now takes a much longer time for the system to eventually settle onto the trivial equilibrium. Figure 4*b* and figure 4*c* bear a close resemblance to earlier pictures illustrating the dynamics of the substructure in the unstable regime. Figure 4*c* indicates the beats increase in their amplitude while the decay slows down.

It is interesting to see the behaviour when the curves in figure 2 cross each other. These points correspond to the case when relation (3.8) holds. For *ζ*=0.05, the first point *p*_{1} can be found at (*p*_{1}, *τ*_{1})≈(0.44537, 7.1859). As expected, the system undergoes a codimension-two Hopf bifurcation and exhibits quasi-periodic oscillations as shown in figure 5*a*. Finally, figure 5*b* shows the dependence of the Hopf frequency on the delay time *τ*. At the points of discontinuity, there are two frequencies of oscillations.

Next, the case when the external force is present (*k*≠0) will be studied. In figure 6, we show the solution *z*(*t*) together with its Fourier spectrum for delay time *τ*=1.55, amplitude of perturbation *k*=0.01, *p*=0.75166. For this value of *p*, the stability switch occurs at *τ*_{c}≈1.604 with the corresponding frequency *ϖ*_{+}=2. We have taken perturbation frequency *Ω*=2.1 to be close to this frequency in order to study possible resonance. For the above-mentioned parameter values, condition (4.7) holds. One can observe that after some transient period, the solution settles onto a stable limit cycle. The Fourier spectrum shown in figure 6 clearly possesses a sharp peak indicating periodic behaviour of the solution. The numerical value of the corresponding oscillation frequency, as found from the Fourier spectrum is *ω*≈2.05. This value lies within 3% of the perturbation frequency *Ω*=2.1.

In the previous case, the value of the time delay was chosen to lie inside the region of asymptotic stability. In order to illustrate the unstable behaviour of the periodic orbit, we take the delay time *τ*=1.608 just outside the stability boundary. For this delay time, condition (4.7) is violated, and the periodic solution, even though it exists, is unstable. The other parameter values are as in figure 6. The temporal dynamics of the solution *z*(*t*) is plotted in figure 7 together with the Fourier spectrum. From both pictures one can observe that the periodic solution is modulated on a longer time scale. The Fourier spectrum now contains both the oscillation frequency and a much smaller frequency of the modulation.

## 6. Stability analysis of the full system

Now, we return to system (2.2), and rewrite it (without loss of generality) by omitting the delay in the second equation and setting *F*_{ext}=0 to obtain(6.1)It is more convenient to rewrite these equations as a first-order system. In order to do so, we introduce the new variables

With these variables, system (6.1) becomes(6.2)

The equilibria for this system are and , *n*=1, 2, …. Linearization near the trivial steady state ** x**=0 giveswhere matrices

*A*and

*B*are given byThe characteristic polynomial becomes(6.3)It is clear that stability is determined by the roots of the second multiplier in equation (6.3).

Similarly, for the steady state the matrices *A* and *B* are given byThe corresponding characteristic polynomial is now(6.4)One can now observe that for *n* odd, which corresponds to an upright position of the pendulum, there are always eigenvalues with positive real part. This implies that such steady states are always unstable for any value of delay *τ*, including *τ*=0. For the case of *n* even, the roots are determined by the same equation as in the case *n*=0 considered above. The second bracket in equations (6.3) and (6.4) is exactly equation (3.3) multiplied by *M* in its dimensional form.

## 7. Experimental results

In the previous sections, we studied the stability of the model using analytical and numerical tools. In order to confirm our findings, we need to perform some experimental tests, using the real-time dynamic substructuring technique. As discussed in §1, the pendulum is the physical substructure, the MSD is modelled numerically in the computer. The numerical model is used to calculate the displacement at the interface due to some external excitation. The displacement is applied to the substructure (pendulum) in real-time using an electro-mechanical actuator (the transfer system). The force acting on the physical substructure is measured via a load cell and fed back to the numerical model. This feedback force is used to calculate the displacement at the interface for the next time step. This process is repeated until the end of the test. To implement real-time tasks a dSpace DS1104 RD Controller Board is used; MATLAB/Simulink is employed to build the numerical model. The dSpace module ControlDesk is used for online analysis and control. The transfer system is an electrically driven ball-screw actuator with an in-line mounted synchronous servo-motor controlled by a servo-drive which applies a displacement to the pendulum pivot point in the vertical direction. The values of the system parameters are given in table 1. It is worth noting that, since the MSD is represented by a numerical model, *M*, *K* and *C* can be changed easily from one test to another to observe different situations. Further details of the experimental set up can be found in Gonzalez-Buelga *et al*. (submitted).

As we have used the assumption that the actuator delay, which dominates other sources of delay, is a constant we show the transfer system response plots in figure 8*a*, the step response and experimental frequency sweep test in figure 8*b*. From figure 8*a*, the approximate transfer function (written via Laplace transform) is found to be or a fixed delay of approximately 0.018 s.

During the experiment, the instability due to delay appears as a new frequency in the numerical model displacement. We start with a delay of 0.018 s (the natural delay time of the transfer system) and the delay will be artificially increased by 0.001 s increments until the system becomes unstable, and the new instability frequency appears. The delay time is deliberately increased by holding the signals going to the actuator so that the stability boundary can be found and a comparison made with the analytical results.

Figure 9*a* shows the theoretical stability boundary in the original non-rescaled system (3.1) with the parameters from table 1. It also shows the values of the stability switches. Figure 9*b* illustrates the way in which we will explore the stability switches experimentally. In particular, we fix the value of *p*=0.225 as determined by the parameters, and start simulations in the stable region for a very small delay time *τ*. Then, we follow path 1 (figure 9*b*) until the stability boundary is reached. Once the boundary is crossed, the response of the system grows exponentially making it impossible to continue experimental testing (the actuator cannot achieve the required displacements). Consequently, it is impossible to increase the delay further in the experimental system, in order to reach the next stability area which starts at *τ*_{2}. In order to find the stability boundaries, we start with the value of delay time *τ* inside the next stability region, i.e. *τ*_{2}<*τ*<*τ*_{3}. We then decrease or increase the delay time accordingly, so that either path 2 or path 3 in figure 9*b* are followed until the instability frequency appears. In this way, the experimental values of stability switches *τ*_{2} and *τ*_{3} will be found.

Figure 10 shows experimental records while following path 1. In figure 10*a*, the measured solution *x* is plotted against its analytical counterpart *y*^{*} as determined by the solution of equation (3.1)—note that in this section we denote analytically computed values with a ( )^{*}. The results show that the numerical model is stable until the delay *τ* reaches the value of 0.093 s, figure 10*c*. At this point an instability frequency of *ω*=6.3 Hz appears, as shown in figure 10*d*. The experimental eigentriple at the first stability boundary is, therefore, as compared with analytical values of .

Figure 11 illustrates the recorded signals when following path 2, when the delay is decreasing. The numerical model is stable until the delay *τ* reaches 0.162 s, figure 11*c*. The corresponding value of the instability frequency is *ω*=5.2 Hz, as shown in figure 11*d*. Hence, the experimental eigentriple for this crossing of stability boundary is as compared with analytical values of . The experimental results when following path 3 are depicted in figure 12. Now the numerical model remains stable until the delay time *τ* reaches 0.255 s, and the solution becomes unstable, figure 12*c*. In this case, the instability frequency is *ω*=6.25 Hz, as can be found from the Fourier transform of the solution at this point, as shown in figure 12*d*. Therefore, this last boundary point in our analysis is characterized by the experimental eigentriple as compared with analytical values of .

Figure 13 shows an excellent agreement between theoretical model and experimental results for stability boundaries. The experimental points follow the stability curve and thereby confirm the existence of stability/instability transitions in the hybrid system. When the stability boundary is crossed, the solution very quickly becomes unstable and develops higher amplitude oscillations and irregular motions. This makes finding experimental values for the stability boundary an extremely sensitive problem and explains the slight deviation of experimental values from the analytical prediction as illustrated in figure 13*a*. In figure 13*b*, we show experimental stability border points corresponding to different values of *p*.

We summarize analytical and experimental values of the instability frequency in table 2.

## 8. Viscous damping

In this section, we consider the effect of including viscous damping in the pendulum model. This can have a significant influence on the velocity feedback of the pendulum–MSD system (see Gonzalez-Buelga *et al*. (submitted) for a detailed discussion). When including viscous damping, the vertical equation of motion for the numerical model in the absence of forcing takes the form(8.1)where *C*_{1} is the damping of the numerical model and *C*_{2} is the viscous damping of the pendulum. Using the same scaling as in §3 and introducing the new variables,equation (8.1) transforms into(8.2)The characteristic equation is now(8.3)Transition to instability occurs when , *ϖ*>0. Substituting this into the characteristic equation and separating real and imaginary parts gives the following system:(8.4)Squaring and adding these two equations, we obtain an equation for *ϖ*, which can be solved to give(8.5)Employing the same methods as the ones used in §3, one can prove the following theorem concerning the stability of the trivial solution.

*Let*. *Then for* |*p*|<1*, the trivial solution is asymptotically stable for any delay time τ*>0. *If, however,**, then the trivial solution of system**(3.2)**is locally asymptotically stable in the region*(8.6)*for all positive delay times τ. In the region**the trivial solution is locally asymptotically stable for values of delay satisfying**where n*=1, 2, ….

It is worth noting that in the limit *ζ*_{2}→0 the results of this lemma coincide with those of lemma 3.2. Furthermore, as it is clear from relation (8.6), these stability results are only valid if one imposes an additional but experimentally reasonable requirement *ζ*_{2}<*ζ*_{1}. Physically, this means that the viscous damping *C*_{2} of the pendulum is smaller than that of the mass–spring–damper.

We illustrate in figure 14*a* the changes to the stability boundary as the viscous damping is increased from 0 until it reaches the value of the damping of MSD. As the viscosity increases, the stability area gets smaller and stability boundary shifts to the left. Eventually, as can be seen in figure 14*b*, after touching zero the stability boundary disintegrates and splits into separate stability regions. As any vibration will eventually die down in these regions, they are called *amplitude death* regions or *death islands* (see, for instance, Xu & Chung 2003). These regions are experimentally important for controlling the stability of the system. Experimental results regarding the stability of the system with viscous damping can be found in Gonzalez-Buelga *et al*. (submitted). We note also that in figure 13*b*, the experimental points are all to the left-hand side of the theoretically predicted boundary. This indicates that (as one would expect) there is a small amount of damping in the physical pendulum which is shifting these points slightly, as explained by figure 14*a*.

We can explain figure 14 physically. The critical stability boundary first touches the *τ* axis (*p*=0) when *ζ*_{1}=*ζ*_{2} and *τ*=*π*. From equation (8.2), we can see that this corresponds to the case when the contribution to the damping, through the term , is exactly balanced by an equal (*ζ*_{1}=*ζ*_{2}) and opposite (out of phase, since *τ*=*π*) contribution, through the term . In this case, and the resulting solution is neutrally stable since . As *ζ*_{2} then increases, there is a finite range of *τ* for *p*=0 when the delayed damping due to the pendulum exceeds that due to the numerical model. In that range, the trivial solution is then unstable.

## 9. Conclusions

In this paper, the stability of a real-time dynamic substructuring model of a coupled oscillator–pendulum system has been investigated. We have modelled the hybrid testing of this system using two-coupled second-order DDEs, with one of them being neutral. This approach includes the assumption that the delay, *τ*, is a fixed value, but allows us to make significant advances in the mathematical analysis of the complex dynamics of the hybrid system. In particular, we started our analysis by considering the situation which occurs for small angles. This results in a neutral DDE for the vertical component of motion. For this equation, we identified regions of stability and instability, parameterized by the delay, *τ*. We established not only the points of stability switching but also showed that there is a finite number of them for fixed system parameters and proposed a scheme for calculating their number. After the presence of a Hopf bifurcation had been established, conditions for the stability of the bifurcating periodic solutions were established using the method of multiple scales. This has also provided analytical expressions for the amplitude and frequency of the periodic orbit.

The numerical simulations in §5 confirm the theoretical findings from §§3 and 4. Several phenomena are shown such as convergence to the steady state in the stable regions, fast unbounded growth of solutions in the unstable regions, quasi-periodic oscillations and the existence of the beats. The Hopf frequency as a function of the delay time shows the presence of two frequencies of oscillation at the codimension-two Hopf bifurcation points. In the case of non-zero forcing, our numerical results illustrate two possibilities: the solutions asymptotically approach a stable limit cycle or they are modulated on a long time scale by a growing harmonic.

In §6, we return to the original system and study its linearized stability. This indicates how the linear stability analysis fits into a more general framework. The practical significance of this is that for systems with many more degrees of freedom, this type of formulation could be extended to identify critical delay and frequency values by numerically approximating the eigenvalues.

The analytical results were then compared to results from a hybrid pendulum–oscillator experiment. The stability boundary of the trivial solution was found to give excellent qualitative and quantitative agreement with the analytical prediction, as were the instability frequencies. When viscous damping is included in the system model, the stability boundaries have been shown to distort. For very large damping, the stability area transforms into disjoint death islands separated by instability regions.

Overall, the study in this paper has focused on the analysis of the pendulum–oscillator system when modelled as a hybrid system with a fixed delay. The close agreement of these results with hybrid-experimental data demonstrate the utility of DDE modelling in predicting the stability boundaries and instability frequencies in hybrid testing models. This information is essential for assessing the viability of particular substructuring test configurations, and in designing suitable delay compensators. Future plans to extend this work include the modelling of large-scale engineering structures, with many degrees of freedom and multiple transfer systems. Practical applications which this will be applied to include cable–deck interaction of cable-stay bridges, and lag damper units for helicopter rotors.

## Acknowledgements

The authors would like to acknowledge the support of the EPSRC: Y.K. is supported by EPSRC grant (GR/72020/01), K.B. is supported by EPSRC grant (GR/S31662/01), A.G.B. is supported by EPSRC grant (GR/S49780) and D.J.W. via an Advanced Research Fellowship. S.J.H. is grateful to Centre de Recerca Matemàtica, Bellaterra, Catalunya for providing facilities to complete part of this project.

## Footnotes

- Received July 6, 2005.
- Accepted November 25, 2005.

- © 2006 The Royal Society