Nonlinear damping and quasi-linear modelling

1. Introduction

In a mechanical system, the damping force is a function of the system’s velocity. This function is nonlinear in a number of mechanical systems and similar nonlinear damping behaviour is also seen in many electrical, biological and other dynamic systems. Examples that will be discussed in more detail below include the damping in aircraft structures, flow through orifices, damping in nanoelectromechanical (NEM) systems and the vibration inside the cochlea.

In some applications, the stiffness of the system also changes to some extent with excitation amplitude, as well as the damping. This nonlinear stiffness can complicate the dynamic response, particularly if the nonlinearity in the stiffness is severe, and can give rise to jump and bifurcation phenomena. In many cases, however, it is the damping that is the dominant source of nonlinearity, and it is interesting and worthwhile to consider the behaviour of systems in which only the damping is nonlinear. We will see that such systems do not exhibit the kinds of jump phenomena seen in systems with nonlinear stiffness. In addition, their response for a given excitation can often be approximated by that of an equivalent linear system, whose parameters depend on the level of excitation. Such ‘quasi-linear’ models of a nonlinear system, which are implicit in many methods of analysis, including equivalent linearization and the describing function, have surprising generality and have been developed in several different fields. These quasi-linear models are very useful in understanding the dynamic response of a number of engineering and biological systems, despite their apparent simplicity.

The nonlinear damping force often increases with excitation level, so that the relative response of the system, compared with the excitation level, is reduced. This is a mechanism for reducing the range of the system’s response, compared with the range of the input, i.e. compressing its dynamic response, which we will see is particularly important in the dynamics of the cochlea.

If there is some reason for the linear damping of a system to become negative, so that the envelope of its linear response would otherwise increase exponentially, nonlinear damping can provide a mechanism by which the output level is stabilized to a fixed level. Rayleigh [1], for example, considers systems in which ‘vibration is maintained by wind (organ pipes, harmonium reeds, Aeolian harps, etc.), by heat (singing flames, Rijke’s tube, etc.), by friction (violin strings, finger glasses) and the slower vibration of clock pendulums and watch balance wheels’. He goes on to say that ‘we may form an idea of the state of things’ by including a damping term proportional to a higher power of velocity in the dynamic equation for a single-degree-of-freedom system. He then considers the specific case of cubic damping, so that in the notation to be used here

where x(t) is the displacement of the system, c₁ and c₃ are linear and cubic damping coefficients, respectively, and m and k are the system’s mass and stiffness, respectively. If c₁ is negative but c₃ is positive, then the system evolves into a steady oscillation, which we would now call a limit cycle oscillation, that Rayleigh showed would have a small third harmonic component, but would otherwise be approximately sinusoidal, at the system’s natural frequency, ω₀, equal to . The fundamental amplitude of the response, x(t), X, was also shown [1] to be given by the solution to the equation a result we will come back to later. More recently, the dynamics of such a ‘relaxation oscillator’ have often been described using the Van der Pol equation [2], whose interesting history is discussed by Ginoux & Letellier [3]: in which the natural frequency is assumed to be unity and the linear damping coefficient is assumed to be equal to −μ, where μ is a positive number. If, following [4,5], equation (1.1) is differentiated with respect to time, and it is assumed that is equal to y(t), then Van der Pol’s equation, (1.3), can be obtained as a special case of the cubic damping equation, (1.1), under the conditions that Although we will mostly focus on stable systems in this paper, with positive values of both linear and nonlinear damping, this connection with Van der Pol’s equations shows that even simple models of systems with nonlinear damping can give rise to a rich variety of behaviour. We will also focus on nonlinear dampers in which the force is a smooth function of velocity, as in equation (1.1), rather than discontinuous forms of nonlinear damping, such as Coulomb friction, although similar methods to those described below can still be used to analyse such nonlinearities.

The aim of this paper is first to review both the variety of mechanisms that give rise to nonlinear damping and their effects. Second, the analysis of systems with nonlinear damping is discussed, emphasizing the calculation of the dissipated power. This leads to a quasi-linear model, for both tonal and random excitations, in which the power dissipation in the nonlinear damper is matched by an equivalent linear damper. Finally, some practical applications are reviewed in which nonlinear damping plays an important role in determining their performance.

2. Examples of nonlinear damping

(a) Aerospace structures

The dynamic response of aerospace structures, such as wings, is important in determining the vibration levels in flight, and in the prediction of flutter [6]. Ground vibration testing is a common way of measuring the dynamic response of such structures and, although many of the results are commercially confidential, there are a number of studies that specifically address the changing damping with excitation level. Fearnow [7], for example, discussed tests on wing and fuselage sections of a Curtiss-Wright C-46D aircraft and showed a structural damping trend that increased nonlinearly with amplitude, as shown in figure 1 [7].

• Download figure
• Open in new tab
• Download PowerPoint

Fellowes et al. [8] discussed the results of a series of free vibration tests on an Airbus partial wing box and showed that the structural damping was also increased at higher amplitudes. Fellowes et al. also note that the reduced response at higher amplitudes, due to the nonlinear damping, can have significant benefits in reducing aircraft flight loads. Traditionally, the damping used in flight load calculations has been based on ground vibration tests, conducted at relatively low excitation amplitudes, and this damping will then be less than that experienced at the higher levels of excitation in flight. A calculation of the flight loads using these lower values of damping will then significantly over-predict the response in flight, which may be seen as providing a level of robustness in the predictions.

(b) Damping in a nanoelectromechanical system

Eichler et al. [9] measured the mechanical response of nanotubes and graphene sheets intended as components of NEM systems, such as atomic force microscopes. Eichler et al. showed that, for components built of both materials, the width of the resonance curves (Res. width), whose natural frequency was in the MHz region, increased with excitation amplitude (as indicated by the AC drive voltage V^AC), as shown in figure 2, for example. Nonlinear damping in these systems then makes it possible to tune the quality factor Q of these resonators.

• Download figure
• Open in new tab
• Download PowerPoint

(c) Acoustic nonlinearity of an orifice

The acoustic flow through a number of different orifices, when excited by sound at different amplitudes, was measured by Ingard & Ising [10]. The flow showed a complicated pattern at very high amplitudes, when the flow becomes supersonic. At very low excitation amplitudes, the acoustic resistance of such an orifice was found to be constant, since the flow was dominated by the viscous loss in the fluid, but then begins to rise in proportion to the acoustic particle velocity at higher amplitudes, owing to the generation of turbulence, as shown in figure 3. The pressure across the orifice is then proportional to the square of the velocity, which is an example of the quadratic dependence of drag force on velocity [11].

• Download figure
• Open in new tab
• Download PowerPoint

(d) Nonlinear suspension and isolation systems

The characteristics of the suspension system are important in determining the ride and handling qualities of vehicles [12]. Part of the suspension system is the shock absorber, which provides both stiffness and damping characteristics. The properties of the shock absorber are often described by a ‘characteristic plot’, of force against velocity. An example of such a plot is shown in figure 4, which shows an almost bilinear damping response, with a ‘rebound’ damping that is significantly greater than the ‘compression’ damping [13]. Surace et al. [13] emphasize that the hysteretic loops in this characteristic plot depend on the frequency of excitation and that a more complete representation would plot the force against damper displacement as well as velocity. In this application, equivalent linearization has been found not to provide an adequate model to predict the vehicle dynamics, however, since it is the details of the waveform which are important, in determining the contact force, for example.

• Download figure
• Open in new tab
• Download PowerPoint

Nonlinear damping has also been incorporated into vibration isolation systems to improve the trade-off between the high-frequency isolation performance and that at resonance, as described in more detail below.

(e) The cochlear amplifier

The mechanical damping in the examples considered above is entirely passive, so that the nonlinear damping force opposes motion and can be considered as a form of negative feedback. Increased damping with level thus requires an expanding nonlinear function, so that the damping force is proportional to the cube of the velocity, for example.

In the cochlea, the mechanical response is amplified by a number of local positive feedback loops (e.g. [14,15]). These loops are formed by the 12 000 or so outer hair cells, whose force response is proportional to the velocity in the system, with a phase inversion, so that they provide negative damping. These cells also include a saturating nonlinearity, owing to their mechanoelectrical conduction channels, as shown in figure 5 [16]. At low excitation levels, the response of this function is almost linear and the negative damping provided by the outer hair cells enhances its response by about 40 dB. As the sound pressure, and thus the excitation of the hair cells, increases, the nonlinearity begins to saturate, lowering the loop gain of the positive feedback loop, reducing the active enhancement of the mechanical response and thus increasing the effective damping. This provides a form of automatic gain control within the cochlea, as discussed in more detail in §5d.

• Download figure
• Open in new tab
• Download PowerPoint

3. Response at single frequencies

(a) Theoretical approaches

The response of a nonlinear system to single-frequency excitation can, in general, be complicated and involve harmonic and sub-harmonic generation and even chaos [17,18]. If the nonlinearity is only in the system’s damping, however, the response is guaranteed to be at the driving frequency and its harmonics, as will be discussed in more detail below. The response at higher harmonics also tends to be smaller than the response at the fundamental driving frequency, and so, for many purposes, only this fundamental response needs to be calculated. It is also important to note that if a system is forced sinusoidally, but the velocity response contains a number of harmonics, the power dissipated in the nonlinear damper, which is equal to the time average product of the force multiplied by the velocity, depends only on the fundamental component of the velocity response, because of the orthogonality of the Fourier series.

Several methods have been developed for calculating the response at the fundamental frequency for such systems, including equivalent linearization, the describing function method and harmonic balance. Although if only the fundamental component is retained in the harmonic balance method, these solutions are equivalent [19,20].

In the following sections, the harmonic balance method will be used to derive the response at the fundamental frequency, initially for a single-degree-of-freedom system having a nonlinear damping described by a power law, and then for an example of a multi-degree-of-freedom (MDOF) system.

(b) Analysis of a single-degree-of-freedom system

The dynamic equation for a single-degree-of-freedom system with a general nonlinear damping can be written as

where m,c₁ and k are the mass, linear damping and stiffness, respectively; x(t) is the response of the system and f(t) is the input force excitation. If the nonlinear damping force is assumed to be proportional to the nth power of the velocity, where n is assumed to be an integer here, we have For odd powers of n, the nonlinear damping in equation (3.2) can be replaced by and for even powers of n, equation (3.2) becomes where the sign function is equal to 1 when the velocity is greater than zero and is equal to −1 if is less than zero. We will investigate the response of the system subject to harmonic (single-frequency) excitation for odd and even power of nonlinear damping. The force excitation is assumed to be sinusoidal and of the form where F is the amplitude of the sinusoidal excitation and θ is the phase shift necessary to ensure that the fundamental component of the response, x(t), has no phase shift.

In the harmonic balance method, the response of a nonlinear system to a sinusoidal input is analysed in terms of its Fourier series, but only a certain number of the output harmonics are retained. Approximating the response at the fundamental frequency, ω, as

and substituting equations (3.3a,b)–(3.5) into equation (3.1), yields and Trigonometric formulae can then be used in equation (3.6a), where The Fourier series for in equation (3.6b) is the same as that for a square wave, having the fundamental component of 4/π cos (ωt). Substituting this into equation (3.6b) leads to an odd power of , so we can again use equation (3.7) to calculate the fundamental component with (n−1) now replacing n.

Ignoring the higher order harmonics, the response at the fundamental frequency can be obtained from

where A_n is a constant coefficient, equal to and Partitioning the sine and cosine terms in equation (3.8), and adding the squares of these two terms, yields For a given power of nonlinear damping, the response of the nonlinear system, X, is given by the solution to equation (3.10). Comparing equation (3.10) with the modulus squared response of a linear system with total damping (c₁+c_eqe), where c_eqe is an equivalent linear damper for the nonlinearity, shows that c_eqe can be written as Defining as the amplitude of the velocity, which is equal to ωX, where leads to the following expression for the equivalent damping in terms of the velocity amplitude: For n=2, for example, the constant coefficient is A₂=3/π, and the equivalent damping becomes where is always positive, and for n=3, the constant coefficient is A₃=3/4, and the equivalent damping becomes which is the same expression derived in [21] and is equal to (3/4)c₃ω²X².

If c₁ were negative, the system would begin to go unstable at the natural frequency, ω₀, but would be stabilized when the power dissipated in the equivalent damper was equal to the power generated by the negative linear damper. The system would thus settle into a limit cycle oscillation with an amplitude X such that (3/4)c₃ω²X² was equal to −c₁, as in equation (1.2) in the Introduction.

In the case of the forced response of a system, the equivalent linear damper is dependent on the velocity response of the system, which is itself dependent on the equivalent linear damper. Figure 6, for example, shows the dependence of the equivalent linear damper, c_eqe, on the velocity response, , in the case of cubic damping (equation (3.15)), with the numerical values given by the specific case considered in the next section. Also shown in this figure is the inverse of the dependence of the response on the equivalent linear damper, when driven at resonance, as in equation (3.16). The two equations are, clearly, both satisfied when the two curves cross. In the case of cubic damping, the crossing point can be calculated analytically, as below, but for higher orders of nonlinear damping, numerical methods need to be used to calculate the solution.

• Download figure
• Open in new tab
• Download PowerPoint

Since the response of the system at the fundamental frequency, is given by the solution to the polynomial equation (3.10), where only one of the possible solutions is real and thus physically realistic, the equivalent damping can, in principle, be calculated directly from this equation. For low powers of nonlinear damping, for example for (n=3), a closed-form analytical solution for the response at resonance can be obtained. The damping force at resonance can be written as

Substituting equation (3.15) into equation (3.16), yields a cubic polynomial in If then the amplitude of the velocity response is given by F/c₁ and changes linearly with the input amplitude. If, however, , the amplitude of the velocity response is given by so that the level curve, of against , has the slope of 1/3 dB/dB, meaning that the velocity response increases by 1/3 dB for every dB increase in the amplitude of the excitation. The transition in the slope of the level curve, from 1 dB/dB at low levels to 1/3 dB/dB at higher levels, occurs when is equal to (4c₁/3c₃)^1/2.

Since the discriminant of equation (3.17) is negative, this equation has only one real root and the other solutions are complex conjugate roots (http://en.wikipedia.org/wiki/Cubic_function). The analytical solution for can then be obtained, in a more general form, from

where and hence c_eqe can be calculated by substituting the value of into equation (3.15).

(c) Numerical simulation for a single-degree-of-freedom system

In practice, for higher order forms of nonlinear damping, deriving an analytical expression for the response becomes more difficult, but direct time domain simulations can still be used to calculate the result. The analytic result for the fundamental component of the response of the system with a cubic damper, derived above, is compared here with the results of a numerical simulation for a particular case. We consider a system with parameters, m=1 kg, c₁=0.1 Ns m⁻¹, k=1 N m⁻¹ and c₃=0.1 Ns³ m⁻³. Figure 7a,b shows the normalized amplitude and the phase of the fundamental component of the velocity, plotted as a function of the excitation frequency, ω, calculated when the input amplitude is F=0.1, 1 and 5 N, using both the harmonic balance method described in §3 and time domain simulation using Matlab ode45. Figure 7c shows the variation of the equivalent linear damping ratio, ζ_eqe=c₀/(2mω₀), as a function of the excitation frequency, ω, for this example at different excitation levels. Since the velocity response is small for excitation frequencies well below or well above the resonance frequency, the equivalent linear damping is, according to equation (3.15), also small at these frequencies, only rising as the excitation frequency approaches resonance and the velocity response becomes more significant. This is because, when the system is driven off resonance, the response is mainly controlled by either the linear stiffness, at low frequencies, or the linear mass, at high frequencies, and the damping plays little part in determining the response. Near the resonance, the analytic and numerical results are very similar, although the response is controlled by the nonlinear damper. The response at the resonance frequency is much greater than at three times the resonance frequency, so that the harmonics at this frequency and the higher harmonics are filtered out by the dynamics of the system to leave a largely sinusoidal response. Figure 7d shows the normalized mean square harmonic distortion (MSD) obtained by analysing the time domain simulations, which is defined as

where is the amplitude of the nth harmonic of the velocity response. This is equal to the square of the ‘total harmonic distortion’, which is normally defined as the ratio of the RMS quantities. The MSD is chosen here for comparison with the ratio of the non-coherent to coherent power, calculated for the case of random excitation in §4c. The MSD is only about 2% for the largest excitation force in these simulations at resonance, so that all the waveforms obtained from these time domain simulations are close to being sinusoidal.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 8a shows the relationship between the level of the velocity response and the level of the input force at resonance, ω=1 rad s⁻¹, for a linear system, with c₃=0, and for the nonlinear system defined above. Nonlinear damping has the effect of limiting the displacement response at high excitation amplitudes, giving rise to the change of slope, from 1 dB/dB at low levels to 1/3 dB/dB at higher levels, predicted in §3b. The equivalent linear damping ratio for the cubic damper at resonance is also plotted in figure 8b, for different levels of input force excitation, which increases with the input amplitude.

• Download figure
• Open in new tab
• Download PowerPoint

(d) Generalization to multi-degrees of freedom

The method of harmonic balance can be extended to the general MDOF system. The dynamic equation, including the nonlinear damping, can be written as

where q(t) is the displacement vector and f(t) is the forcing vector (both of size N×1 for an N degree-of-freedom model), M, K and C are the linear mass, stiffness and damping matrices, respectively, and are the nonlinear damping forces.

Backbone curves are plots of undamped natural frequency against excitation amplitude, and their shape can indicate the presence of jump phenomena or bifurcation in the forced response [22]. The backbone curves are derived by solving the differential equations for the undamped forced system. The backbone curves are straight vertical lines for linear systems, but in the presence of stiffness nonlinearities, such as hardening (or softening) nonlinearities, the backbone curves bend to the right (or left). Assuming that all the nonlinearities are in the damping matrix, however, so that the mass and stiffness matrices are constant with amplitude, the backbone curves would remain as straight lines, and so the system is guaranteed to have no jump or bifurcation behaviour. This may be a significant advantage of implementing nonlinear damping as in some applications the jump phenomenon is undesirable, since it may be required to have a response that is independent of the previous excitation of the system, for example. To demonstrate the harmonic balance method for MDOF systems, we have assumed a special case where the nonlinear damping matrix is diagonal, so that h_j is a function of only . This assumption is adopted in order to simply illustrate the method, but is not necessary, since the method can still be applied in the more general case of non-diagonal damping. If an equivalent linear damping matrix is now assumed to replace the nonlinear function, then the equivalent linear system is

Assuming nth power law damping and that the responses are at the fundamental driving frequency ω, the elements of the matrix C_eqe can be obtained from the harmonic balance method as described above for each nonlinear term where A_n is calculated from equation (3.9) and q_j is the amplitude of the displacement for the jth degree of freedom. For harmonic excitation, the velocity vector of the equivalent linear system is given by For MDOF systems, even for low powers of nonlinear damping, it is difficult to obtain a closed-form solution for the equivalent damping, owing to the coupling between the degrees of freedom, and hence iterative methods are required. An iterative method that can be used to obtain the equivalent linear damping for the nonlinear damping first involves calculating the response of the linear system, including linear damping C and considering C_eqe=0. The velocity amplitudes are then calculated from equation (3.25), and these amplitudes are substituted into equation (3.24) to obtain the elements of an equivalent damping matrix C_eqe. The updated damping matrix is then used to generate the next response amplitudes. This procedure continues until convergence is achieved.

(e) Numerical simulation for a multi-degree-of-freedom system

Consider the two d.f. system shown in figure 9. A harmonic excitation acts on the mass m₁ with the forcing amplitude of F and excitation frequency of ω. The displacements of the two masses are denoted by q₁ and q₂, respectively. A linear spring k₁ and a linear damper c₁₁ are attached to the mass m₁, whereas a linear spring k₂ and a nonlinear damper connects the two masses m₁ and m₂. The nonlinear damping force between the two masses is assumed to be .

• Download figure
• Open in new tab
• Download PowerPoint

The dynamic equations in terms of the two coordinates q₁ and q₂ can be written in matrix form as

Since the nonlinear damping depends on the relative velocity between the two masses, it is convenient to introduce the transformation Therefore, the dynamic equations can be rewritten as So that The equivalent damping matrix can be obtained using equation (3.24), Therefore, the equivalent system matrix becomes where . The equations can be written in the normalized form, by dividing the first equation by m₁ and the second equation by m₂ to give where the parameters are To obtain the frequency response, an iterative method is used at every individual excitation frequency. The amplitude of the velocities can be obtained from equation (3.32) at every iteration and then be used to construct the equivalent damping matrix, which is used to recalculate the response from equation (3.32). The iteration continues until convergence is achieved. For numerical simulation, the following system parameters are used: so that the natural frequencies of the system are 0.618 and 1.618 rad s⁻¹. The amplitudes of the resulting mobilities and are plotted in figure 10. A very good agreement is achieved between the analytical approach, using the harmonic balance method (red solid line), and the numerical approach, using time domain simulation (blue dashed line).

• Download figure
• Open in new tab
• Download PowerPoint

The equivalent damping ratio, ζ_2eqe, at every iteration is plotted in figure 11a, for excitation at the two undamped frequencies, which indicates a good convergence after about three iterations in this case. In addition, the equivalent damping ratio c_2eqe, after convergence, is plotted in figure 11b as a function of the individual excitation frequencies ω. The equivalent damping ratio when the system is driven at the natural frequencies is higher than at other frequencies, since the nonlinear response is greater at resonances.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 12a shows the response level curve for the first mass, in comparison with the linear case. At high excitation amplitudes, the nonlinear response is again less than the linear response. The total equivalent damping is also obtained numerically at every excitation amplitude, as shown in figure 12b, and again increases with excitation amplitude, as expected, but in a slightly different way from that for the single-degree-of-freedom system, shown in figure 8b.

• Download figure
• Open in new tab
• Download PowerPoint

4. Response to random excitation

(a) Theoretical approaches

In this section, the response of systems with nonlinear damping subject to random excitation will be considered. In practice, many engineering systems are subjected to loadings that are random in nature and various methods to predict the response of nonlinear systems subject to random excitation have been developed [23,24]. A probabilistic description of the random process can be obtained from the statistical moments [25]. When the input is a Gaussian process and the system is linear, it is well known that the response also has a Gaussian probability density function (PDF). This enables closed-form solutions for the statistics of the response to be calculated. However, in many engineering systems, owing to the nonlinear structural behaviour, the response to a Gaussian input is non-Gaussian. For weakly nonlinear systems, however, the response is still approximately Gaussian and the method of ‘equivalent linearization’ or ‘stochastic linearization’ can be used to approximate the response [23,24]. This method can only predict the response moments up to second order with reasonable accuracy, however. To take into account the non-Gaussian characteristic of the nonlinear response, the alternative method, of equivalent nonlinearization, has been developed [25,26]. Equivalent nonlinearization is based on replacing the nonlinear system with a simpler equivalent nonlinear model, which depends on the energy level of the system. However, this method is restricted to single-degree-of-freedom systems, as discussed in [24]. A different approach to random vibration problems is the use of Fokker–Planck–Kolmogorov (FPK) theory, to obtain the state transition PDF. For white noise excitations and for particular nonlinearities in the stiffness and damping, the method provides a direct approach to obtain the exact response [27]. However, for general nonlinear systems, numerical approaches are still required to solve the FPK equation [28].

Similar problems have been treated in a different way for communication systems, using the invariance property of separable random process [29]. Nuttall [29] used the separable class of random processes to prove their invariance under nonlinear transformations [29,30] and showed that the cross-correlation function between the input and output of a single-valued nonlinearity is proportional to the auto-correlation of the input. This invariance property allows the nonlinear function to be replaced with a level-dependent linear gain. A number of different approaches to analysing a single-degree-of-freedom system are compared in the next section and the links with the cross-correlation approach are emphasized [31].

(b) Analysis of a single-degree-of-freedom system

The equation of motion of a single-degree-of-freedom system with nonlinear damping that depends only on the velocity has the general form

where the forcing f(t) is now considered to be a stationary random process. In order to perform a quasi-linear analysis, the nonlinear damping that appears in equation (4.1) must be replaced with a linear damping term that is optimal in some sense. There are three possible approaches to determining the appropriate linear damping coefficient which all lead to the same mathematical result: (i) error minimization, (ii) series truncation, and (iii) power balance. It is useful to consider each of these approaches in turn, since taken together they provide a physical insight into the nature of the approximation process. With the error minimization approach, the nonlinear damping function is approximated by a linear function and the error involved in the approximation of the nonlinear damping force is written as The coefficient c_eqe is now found by minimizing the mean squared error, which yields where E[] represents the ensemble average and ε is the error. Therefore, since the equivalent linear damping can be obtained by substituting this and equation (4.2) into equation (4.3) and solving for c_eqe to give With the series truncation approach, the nonlinear damping force is written in the form where the functions are a set of orthonormal polynomials (the nth function being of order n) with weighting function a_n. Thus, by definition, and where is the PDF of the response and δ_nm is the Kronecker delta function, which is equal to 1 when n=m, and otherwise is zero. If the random response is stationary and the velocity has zero mean, it can readily be deduced that the first two orthonormal polynomials are and With this approach, the linearization of the system is achieved by retaining only the first two terms in equation (4.6). Assuming that has zero mean, it is readily shown from equations (4.5)–(4.8) that so that the orthogonal series approach agrees with the error minimization approach. This is clearly to be expected, since from equations (4.2), (4.5) and (4.6) the mean squared error can be written in the form so that c_eqe=a₁ is the optimal linear approximation. Equation (4.5) is often presented without reference to the underlying orthonormal series, as in [32], where the method is referred to as equivalent linearization. As an aside, it can be noted that if the nonlinear damping is a function of the velocity, rather than an instantaneous, memory-less function, then equation (4.5) would be replaced by a functional series which is analogous to the Wiener series [33], and equation (4.7) would be replaced by a cross-correlation approach analogous to the Lee–Schetzen algorithm [30].

The third and final approach to obtaining the quasi-linear damping coefficient is to note that the power dissipation rates under linear and nonlinear damping are, respectively,

and If the constant c_eqe is chosen to equate the ensemble average of the two powers, E[P_L]=E[P_NL], then the result is which is in agreement with equation (4.4). Thus, the mathematical optimization technique leading to equation (4.14) ensures that a quantity of great physical importance, the dissipated power, is conserved. In the expression for c_eqe, and can be interpreted as the cross-correlation between the nonlinear damping force and the velocity and the auto-correlation function of the squared velocity, both at zero lag, in agreement with Nuttall’s approach mentioned above.

If the nonlinearity were to be in the stiffness, rather than the damping, then the conserved quantity, in the equivalent to equation (4.14), would be E[h(x)x], which does not have such physical significance, thus suggesting that the quasi-linear approach can be expected to be more effective for nonlinear damping than for nonlinear stiffness.

One key feature of the linearization approach represented by equation (4.14) is that the optimal linearization constant c_eqe depends on the statistics of the response: the PDF is needed to evaluate the expectations that appear in the equation. A common approximation is to assume that the velocity is Gaussian, in which case it can be shown that equation (4.14) becomes

which is called stochastic linearization [23,24]. The resulting expression for c_eqe can be expressed as a function of the standard deviation of the velocity.

The foregoing presentation of the quasi-linear technique could also have been employed for the harmonic case presented in §3b. In that case, equation (4.5) would be replaced by a Fourier series expansion of the nonlinear damping, and the ensemble average E[] would be replaced by a time average. The assumption that the response is Gaussian is replaced by the assumption that the response contains only the first harmonic of the excitation frequency, and the harmonic balance approach ensures that (under the harmonic response assumption) the time average of the power dissipated over a cycle of the excitation is preserved.

An analytical expression for c_eqe can be obtained in some cases by substituting for into equation (4.15). With cubic damping, for example, the following analytical expression for c_eqe can be derived for white noise excitation having a velocity variance of as

which can be solved using integration by parts, as shown in appendix A, to give . The ‘total equivalent’ damping consists of the linear damping, c₁, and the ‘equivalent linear damping’, c_eqe, In a single-degree-of-freedom system, can also be expressed as a function of c_eqe, as described in appendix A, to give a closed-form expression for c_eqe, as where S_ff is the uniform power spectral density of the random, white noise, excitation force.

Figure 13b shows the increase of the equivalent linear damping, calculated from equation (4.18), as the power spectral density increases, for the numerical example considered in §4c. At high levels of random excitation, the equivalent linear damping can be obtained by assuming that c₁ is negligible in equation (4.18), so that c_eqe is equal to (3c₃S_ff/2m)^1/2. The power spectral density of the velocity at resonance, , can also be obtained from

Substituting the equivalent linear damping, equation (4.18), into equation (4.19), yields For purely cubic nonlinear damping, (c₁=0), the auto-power spectrum of the velocity at resonance becomes 2 m/3c₃, which is independent of the power spectral density of the random excitation, S_ff. Figure 13a shows the amplitude level of the auto-power spectrum of the velocity response at resonance, for the numerical example above when plotted as against . At low level of excitation, the slope is 1 dB/dB; however, at high levels of excitation the slope changes until reaches the constant value of 2 m/3c₃. The overall mean square velocity, given by , is provided using the results in appendix A as c_eqe/3c₃, so that for high levels of forcing, when c_eqe is given by (3c₃S_ff/2 m)^1/2, then becomes equal to (S_ff/6c₃m)^1/2 and its level increases at 1/2 dB/dB compared with that of the force power spectral density.

• Download figure
• Open in new tab
• Download PowerPoint

More generally, an iterative technique can [24] be used to obtain the power spectrum of the response, and hence c_eqe. For the first iteration (j=1), it is assumed that c_eqe(1)=0, and the auto-power spectrum of the velocity is obtained from

and the cross-power spectrum of the velocity and force can be obtained from where, in general, and c_eqe(j) is the equivalent linear damping at the jth iteration. The variance of the velocity can then be obtained from The next iteration of the equivalent linear damping c_eqe(2) for the nonlinear function can then be calculated from where it is assumed that is Gaussian at every iteration. This value is then used to obtain the next power spectrum and the variance of the velocity. This iteration continues until the convergence is achieved.

(c) Numerical simulation of a single-degree-of-freedom system

As an example of the iterative technique, consider a system with cubic damping having the same parameters used in §3c, m=1 kg, c₁=0.1 Ns m⁻¹, c₃=0.1 Ns³ m⁻³, k=1 N m⁻¹ and an intensity of input S_ff=1 N². The values of and have been calculated at each iteration, as described above, and the results used to calculate the equivalent damping and the variance of the velocity, as shown in figure 14. The final value for the equivalent damping and the variance are found to be c_eqe=0.344 Ns m⁻¹ and , which are in close agreement with the theoretical values obtained from equations (4.18) and (4.24). The iterative method converges to the closed-form solution in this particular case of cubic damping, but if the power of damping is high, then the analytical solution is difficult to obtain although the iterative method can still be used.

• Download figure
• Open in new tab
• Download PowerPoint

The amplitude and phase of the normalized H₁ frequency response estimate, [34], is plotted as a function of the frequency in figure 15a,b, for three power spectral densities of the input excitation S_ff=0.1, 1 and 5 N², using both the analytical method described in §4 and time domain simulation using Matlab ode45. Simulations were performed with a Gaussian random signal of 2×10⁷ points, having a sampling frequency of 100 Hz, divided into 128 blocks for averaging, and using a Hanning window with 50% overlap between each segment. The analytical results are in good agreement with time domain simulations in figure 15a,b. The equivalent damping is calculated for different power spectral densities, using equation (4.18), which agrees well with that obtained from the simulation, and is plotted in figure 15c against frequency, for comparison with figure 7c, even though it is independent of frequency in this case. Figure 15d shows the ratio of the non-coherent to the coherent power (NCP) calculated as [35]

where γ²(ω) is the coherence function [35], calculated from the cross- and power spectral densities obtained by analysing the results of the time domain simulation as This NCP ratio for the velocity response is relatively independent of frequency for S_ff=5 N² and has a value of about 5% at the resonant frequency, corresponding to a coherence, γ²(ω), which is also relatively independent of frequency, and has a value of about 0.95. The NCP can be compared with the normalized MSD plotted in figure 7d for tonal excitation, where the MSD changed significantly with excitation frequency in that case, as did the equivalent linear damping.

• Download figure
• Open in new tab
• Download PowerPoint

(d) Generalization to multi-degrees of freedom

For an MDOF system,

the method of stochastic linearization can be used to obtain the response of the system to random excitation having the power spectral density matrix of S_ff(ω)=E[f(ω)f^H(ω)].

The equivalent linear system equation can be written as

The elements of the matrix C_eqe can be obtained from The spectral density matrix of the output, S_qq(ω)=E[q(ω)q^H(ω)], is given by where ()^H denotes the Hermitian, complex conjugate transpose and The cross-correlation of the responses can be calculated from Consider the previous two-degrees-of-freedom system shown in figure 9. We again assume a special case of nonlinear diagonal damping, in which the jth term of the nonlinear damping is a function of the jth velocity. The force, which acts on the mass m₁, is now considered to be random white noise. A cubic nonlinear damping is considered between the two masses. Therefore, The equivalent damping matrix can be obtained using equation (4.23), The total equivalent damping matrix C_e becomes where . To obtain the equivalent damping, we need to calculate By definition where S_uu is the spectral density of u(t)=f₁(t)/m₁ and G₁₁ is the term in the first row and column of G(ω) in equation (4.32). Similarly, where G₂₁ is in the second row and first column of G(ω). Substituting G₁₁ and G₂₁ into the variance yields where Π(ω) and Λ(ω) are, respectively, the numerator and denominator polynomials of ω²|G₂₁(ω)|². An analytical closed-form solution for the variance and is provided in appendix B. Since the two variances are dependent on the equivalent damping c_2eqe, an iterative approach is used to solve the equations. For the first iteration, the initial value of c_2eqe=0 is considered and initial estimates of and are obtained. Then these values are used to update the equivalent damping c_2eqe. The iteration continues until convergence is achieved. For numerical simulation, the following system parameters are used: The equivalent damping c_2eqe and the two variances and are plotted at each iteration in figure 16. The results converge after about six iterations. The first iteration shows the values for the linear system, when cubic damping c₃=0.

• Download figure
• Open in new tab
• Download PowerPoint

The auto-spectrum and cross-spectrum of velocities from the iterative analytical method after convergence are plotted in figure 17 with the red solid line, and the results are compared with time domain simulations of this system using Matlab ode45, denoted with the blue dashed line. The analytical and numerical results are in good agreement.

• Download figure
• Open in new tab
• Download PowerPoint

5. Practical applications of nonlinear damping

(a) Microspeakers

The small loudspeakers that are required to reproduce speech and, increasingly, music in mobile phones are called microspeakers, and their construction is illustrated in figure 18a. In general, the nonlinearity of larger sized loudspeakers is due to the number of mechanisms, particularly the nonlinear stiffness of the diaphragm surround and the nonlinear coupling between the coil and the magnet [36]. In microspeakers, however, these sources of nonlinearity are not so important and the dominant nonlinearity is due to the nonlinear damping, as shown by the measurements of Klippel, reproduced in figure 18b [36]. This nonlinearity in the damping will reduce the response of the microspeaker at higher drive levels and so limit its throw and hence its distortion. Microspeakers generally have small holes in the construction, as indicated by the air flow through the leak in figure 18a, and the mechanism of nonlinear damping is probably due to the nonlinear flow through these holes, as discussed in §2c.

• Download figure
• Open in new tab
• Download PowerPoint

In a number of applications, the practical effect of nonlinear damping can be understood by considering the relationship between the level of the output, in dB, as a function of the level of the input, in dB. An idealized version of this level curve for a microspeaker is shown by the solid line in figure 19a, where the dashed line is the 1 dB/dB response of an entirely linear system. If the microspeaker is driven at resonance and it is assumed to have both linear and cubic damping, then the output level is less than that of the linear system above some limiting value of the driving signal, beyond which the response level only rises by 1/3 dB for every dB increase in the drive level, as discussed in §3, and so will be significantly restricted if the microspeaker is driven hard. For tonal excitation well away from the resonance frequency, the response is again dominated by either the stiffness or mass, both of which are linear, so that the level curves always have a slope of 1 dB/dB in this and the other cases discussed below.

• Download figure
• Open in new tab
• Download PowerPoint

(b) Vibration isolation

The transmission of vibration from a source structure to a receiving structure is often controlled with vibration isolators. In its simplest form, such an isolator consists of a spring and damper in parallel. If the damper is linear, there is a well-known trade-off between choosing a high damping coefficient, to control the response at resonance, and choosing a low damping coefficient, to reduce the vibration transmission well above resonance [37]. A number of methods of overcoming this problem have been suggested, including the use of cubic nonlinear dampers [38]. When driven one frequency at a time, cubic damping will produce a high equivalent linear damping value at resonance, when the response level is high, but a low equivalent linear damping value at excitation frequencies well above resonance, where the response level is low, as required. Less advantage is gained for broadband excitation, however, since, as seen in §4, the equivalent linear damper will then have a single value, mainly governed by the resonant response, resulting in poor high-frequency isolation. Experimental demonstrations of such systems have been developed using an electromagnetic actuator and nonlinear velocity feedback [39].

The level curve for such a device, when driven by a single frequency at resonance, will look similar to that in figure 19a, if the damper has both linear and cubic components, so that the transmitted vibration is significantly reduced for high-level excitation in this case.

(c) Energy harvesting

Electrical energy may be harvested from ambient motion using an inertial system, in which the damping is provided by an electromechanical transducer attached to an electrical load. The scale of such systems can vary from miniaturized devices for low-power wireless sensors [40], to devices designed to harvest power from human motion [41], to large-scale devices powered by ocean waves or large structures [42,43]. If driven at resonance, the harvested power is, in principle, maximized if the electromechanical damping is very small, so that the motion of the inertial mass is very high [44], but the maximum extent of this motion, the maximum throw, is generally limited by practical construction constraints. Sufficient damping then needs to be introduced so that, at the highest level of excitation, the response is within the maximum throw of the device.

In many applications, however, the excitation of such an energy harvester is not stationary, and the maximum throw is only observed at the highest levels of excitation, which are not experienced for the majority of the time. Under these conditions, the use of a nonlinear electromechanical damper can increase the response of the harvester, and hence its power output, when it is driven at excitation levels less than those which would give the maximum throw [21]. Figure 19b, for example, shows an idealized version of the level of power that can be harvested from such a device (solid curve) compared with that harvested using a linear damper (dashed line), where both devices are designed to have the same maximum throw at the maximum excitation level. At this level of excitation, the linear and nonlinear devices produce almost the same power output, since the power dissipated in the nonlinear device is the same as that in the equivalent linear damper. For levels of excitation below this maximum, the power from the nonlinear harvester, which is assumed to have cubic electromechanical damping, falls by 2/3 dB/dB [21] and can produce more output power than the linear device.

At very low levels of excitation, the dynamics of such a harvester is, in practice, dominated by parasitic mechanical sources of damping. If this parasitic damping is linear, the slope of the level curve returns to 1 dB/dB once the excitation level falls below a certain limit, as illustrated in figure 19b. If, however, if there are elements with Coulomb friction in the parasitic damping, the motion of the harvester and hence the power output drops to zero at very low levels of excitation [45]. Despite these limitations, it is clear that the use of nonlinear electromechanical damping has the potential to increase the harvested power over a range of excitation levels.

(d) The mechanical response of the cochlea

The cochlea is a coiled structure in the inner ear that converts sound into neural signals. It has remarkable sensitivity and selectivity and also has a huge dynamic range, of about 120 dB, compared with the operating dynamic range of the 3000 or so inner hair cells that convert the internal motion of the cochlea into neural signals, which is only about 30 dB. There are several mechanisms that provide the required compression in the dynamic range of hearing, the most important of which is probably the nonlinearity in the cochlear amplifier [46–48]. As noted above, the cochlear amplifier is driven by the 12 000 or so outer hair cells acting as local positive feedback loops.

Figure 19c shows an idealized level curve in this case, between the input sound pressure level and the resulting level of vibration at a point inside the cochlea. Similar level curves have been observed experimentally, using laser measurements of cochlea motion [48,49]. At low sound pressure levels, the cochlear amplifier can provide about 40 dB of enhancement to the vibration in the cochlea, as the positive feedback loop operates with a high gain over an almost linear part of the operating curve of the hair cell, which is shown in figure 5.

When the sound pressure level reaches about 30 dB, however, the nonlinearity on the operating curve starts to have an effect and the gain of the positive feedback loop is reduced, causing the enhancements in the response to be less than 40 dB. For sound pressure levels above about 90 dB, the motion inside the cochlea is much larger than the range of the operating curve shown in figure 5, and the cochlear amplifier is then completely saturated and can provide no enhancements to the motion. For sound pressure levels above 90 dB, the response is thus almost linear, as indicated by the dashed, 1 dB/dB, line in figure 19c. For sound pressure levels between about 30 and 90 dB, the slope of the level curve is about 1/3 dB/dB. A 60 dB dynamic range of sound pressure level is thus converted to a 20 dB dynamic range of cochlear motion, which can be faithfully converted into neural signals by the inner hair cells.

It is the outer hair cells that are mainly damaged by loud sounds or the ageing process, and this damage is one of the main causes of hearing loss. In the absence of these outer hair cells, the cochlear amplifier no longer functions and the level curve in figure 19c will revert to the dashed, linear, line. Not only will the low-level amplification then be lost, which could be compensated for with a high-gain hearing aid, but the compression effect of the cochlear amplifier will also disappear, making a higher level sound feel uncomfortably loud if such a high-gain hearing aid is used.

Nonlinear damping thus plays an important role in the normal functioning of our sense of hearing and has been used for some years to explain the dynamic behaviour of the cochlea [5,50,51]. In fact, the analysis of cochlear mechanics is more complicated than suggested above, since the dynamic response at each point along the cochlea is coupled with all the other points, owing to the fluid motion in the cochlea fluid chambers. A full nonlinear analysis has to take account of this interaction, for example, by iteratively taking a first approximation to the nonlinear local behaviour, linearly coupling each of these via the fluids and then obtaining a better approximation to the nonlinear local responses at each individual position along the cochlea. This iterative quasi-linear approach was originally developed for tonal excitation and for combinations of tones [52,53] but has also been used for broadband random excitation [54,55], using cross-correlation methods with a theory that is similar to those discussed in §5, which in this field is called the equivalent nonlinear or EQ-NL theorem [56]. Similar iterative methods have been used for other distributed nonlinear systems, such as the response of underwater structures and cables that are damped to different extents along their lengths due to the local flow velocity, which itself depends on the motion [57].

6. Conclusion

Nonlinear damping arises due to a variety of mechanisms in engineering and biological systems and commonly acts to limit the resonant response at higher excitation amplitudes. It thus compresses the dynamic range of the response compared with that of a linear system. Even though the stiffness of some physical systems, as well as the damping, can change with amplitude, it is interesting to consider the effect of nonlinear damping alone, since this can provide some clear physical insights and is useful in a number of practical applications. It is shown that, if the nonlinearity of a system is confined to its damping, the backbone curves of such a system are vertical straight lines and so the system will not exhibit any jump or bifurcation phenomena. This significantly simplifies their analysis. The analysis of systems in which nonlinear damping is proportional to the nth power of velocity is considered, initially for sinusoidal and then for random excitation. In both cases, the system is analysed in terms of an equivalent linear damper, whose value changes with excitation level. Although the focus of this paper is on systems with positive values of both linear and nonlinear damping, which are thus stable, it is noted that the amplitude response of an unstable system, with negative linear damping, can be limited by positive nonlinear damping, as described by Rayleigh [1]. In fact, amplitude stabilization of such instabilities occurs in the human cochlea, where the resulting limit cycle oscillations are known as spontaneous otoacoustic emissions [15].

For sinusoidal excitation, it is seen that such an equivalent linear damper can be obtained by retaining the fundamental term in a harmonic balance analysis, which is equivalent to the describing function approach used in control engineering. Closed-form solutions for the equivalent linear damper can be obtained, in terms of the driving level, for a single-degree-of-freedom system will low-order nonlinear damping, but iterative methods need to be used for higher order nonlinearities or in MDOF systems. For random excitation, it is shown that an error minimization technique, a series truncation approach and an equivalent power method all give the same expression for the equivalent linear damper, which is equal to the cross-correlation between the damping force and the velocity, divided by the auto-correlation of the velocity, both at zero lag. Closed-form solutions for the equivalent linear damper can again be obtained for single-degree-of-freedom systems with low-order nonlinear dampers, but iterative methods must be used in more complicated cases, in which case the assumption of a Gaussian response appears to be a reasonable one. The quasi-linear method using equivalent linear damping has been shown to accurately predict the responses from time domain simulations of a single-degree-of-freedom system with cubic damping at different levels of forcing, for both sinusoidal (figure 7) and random excitation (figure 15). Also shown in these figures are the low levels of mean square distortion in the responses.

The use of nonlinear damping in microspeakers, vibration isolation systems and vibration energy harvesters illustrates its practical application in extending the dynamic range of these devices using this mechanism. The nonlinear damping in the mammalian cochlea has a rather different origin, since it is due to the saturation in the positive feedback provided by the outer hair cells that amplify the mechanical motion, but once again this acts to compress the dynamic range of the response and provides a crucial aspect of our hearing.

The equivalent linear damper provides a very convenient quasi-linear model of systems with nonlinear damping. The idea also has strong theoretical support if the equivalent linear damper is arranged to dissipate the same power as the nonlinear damper, at the specified excitation level. Although similar linearization techniques can be used for systems with nonlinear stiffness, the conserved quantity, which can be written as E[h(x)x] in this case, where h(x) is the nonlinear stiffness, does not have such a clear physical significance as the power, , which is the conserved quantity for systems with nonlinear dampers. The equivalent linearization method is thus particularly well suited to the analysis of nonlinear damping [23,24,30].

Although the full behaviour of a system with nonlinear damping, such as the generation of harmonics or the non-Gaussianity of the random response, will not be completely captured by such a quasi-linear approach, it does provide a very useful engineering tool for their first-order analysis, since it reproduces the power dissipated at different points in the system. Clear physical insights can thus be brought to systems with nonlinear damping: as quasi-linear systems with equivalent linear dampers whose value depends on the local velocity response.

Authors' contributions

S.J.E. drafted the original paper, suggested the examples and contributed to the theoretical development. M.G.T. was most involved in developing the majority of the theoretical formulations and conducted the numerical simulations. R.S.L. provided physical insight into the method and suggested the theoretical approach for random data.

Competing interests

We declare we have no competing interests.

Funding

The authors are grateful for the support provided through the EPSRC programme grant ‘Engineering Nonlinearity’ (EP/K003836/1) and the first grant scheme (EP/ K005456/1).

Footnotes

One contribution of 11 to a theme issue ‘A field guide to nonlinearity in structural dynamics’.