## Abstract

Nonlinearities in practical systems can arise in contacts between components, possibly from friction or impacts. However, it is also known that quadratic and cubic nonlinearity can occur in the stiffness of structural elements undergoing large amplitude vibration, without the need for local contacts. Nonlinearity due purely to large amplitude vibration can then result in significant energy being found in frequency bands other than those being driven by external forces. To analyse this phenomenon, a method is developed here in which the response of the structure in the frequency domain is divided into frequency bands, and the energy flow between the frequency bands is calculated. The frequency bands are assigned an energy variable to describe the mean response and the nonlinear coupling between bands is described in terms of weighted summations of the convolutions of linear modal transfer functions. This represents a nonlinear extension to an established linear theory known as statistical energy analysis (SEA). The nonlinear extension to SEA theory is presented for the case of a plate structure with quadratic and cubic nonlinearity.

## 1. Introduction

### (a) Energy scattering phenomenon

This paper addresses a particular type of nonlinear high-frequency vibration problem. The nonlinearity of interest produces an effect termed ‘energy scattering’ or an ‘energy cascade’ [1] which is characterized by a power input to one frequency band, through either impact or random forcing, resulting in a response that contains energy across a wider frequency range than the bandwidth of excitation. This paper contains a possible extension of an established energy balance technique used for linear vibration analysis [2] to the nonlinear energy cascade scenario.

It is known that the energy cascade phenomenon in structural elements such as plates arises from large amplitude vibration, although the resulting strain is small and can be modelled using Green’s strain [3]. The resulting Von Kármán equations of motion are capable of describing nonlinear effects [4] due to either: (i) large amplitude vibration in a uniform structure giving rise to a dominant cubic nonlinearity or (ii) curvature in the profile of the structure resulting in an additional quadratic nonlinearity, where the contribution of cubic and quadratic behaviour to the vibration response will depend strongly on the amplitude of excitation.

An example structure of type (ii) above is shown in figure 1*a*. It consists of a nearly flat plate, with a slightly wavy profile, as shown by the mesh plot. This structure is a musical instrument called a ‘thundersheet’ [5] and shows nonlinearity that is revealed in both impact and steady forcing test regimes. While this example in particular is a musical instrument, its dimensions (0.8 m×0.5 m×3 mm) mean that it resembles a typical engineering component of a built-up structure, for example an aircraft or satellite panel.

The thundersheet shown in figure 1*a* has been subjected to vibration tests by suspending the structure from two points (which approximates free edge boundary conditions). The point acceleration response of the structure has been measured in one corner, when excitation is applied by a shaker attached at the location shown in figure 1*a*. The data have been post-processed by performing a fast Fourier transform (FFT) of the time domain acceleration signal.

A random band-limited force with frequency content in the range 110–140 Hz has been considered. The frequency response is shown in figure 1*b*, where energy is found in frequency bands (region ii) other than those of the initial excitation (region i). Figure 1*b* can be interpreted in terms of the overall power balance within the system: there is an energy input to a given frequency range, energy flows between frequencies and finally energy is removed by structural damping. As this phenomenon arises so readily it perhaps implies that many slightly irregular structures will show a similar behaviour, and this is supported by the literature on the topic, as discussed in what follows.

The presence of nonlinearity giving rise to a cascade of energy into higher frequencies has been studied in detail by Camier, Touze and colleagues [6,7], for the case of structural plate elements. The various categories of behaviour exhibited as the vibration amplitude is varied have been explored in [8], while [6,9] reveal how the spectra of the energy response show self-similar scaling laws for both transient and steady-state responses. The self-scaling laws of [9] are compared favourably with experimental results for the case of harmonic excitation of isolated plates, at amplitudes that produce frequency content across a range of the order of several kilohertz.

While the self-similar scaling results of [8] are potentially applicable to the single-plate structures considered in this paper, no general formulation has been developed for different bandwidths of excitation, and consideration has been given to the effect of uncertainty in the geometry and properties of the system. Furthermore, the more general case of the evolution of energy cascades through a built-up structure is considered here, and this has not been investigated in the context of scaling laws.

Beyond the extensive consideration given by Touze and co-workers there exists literature on energy flow between frequencies within structures, in particular the transition to a state of equipartition across all frequencies (rather than just within a frequency band). A significant volume of work is built upon the Fermi–Pasta–Ulam problem [10], studying thermal equipartition in crystal lattices, modelled as a set of nonlinearly coupled oscillators. However, while the underlying equations of motion are similar to those arising for nonlinear structural vibration, differences exist such as the absence of damping, lack of external forcing and the use of special distributions of oscillator frequencies to encourage equipartition. Other studies have been performed into equipartition criteria, for instance by Magionesi & Carcaterra [11], where equipartition is found as a limit for large values of nonlinear parameters (i.e. high nonlinearity). While the state of equipartition is relevant, it is only one particular energy distribution of interest and not universal across all forcing bandwidths and amplitudes.

### (b) Analysis techniques for nonlinear energy cascades

A range of possible analytical solution techniques are available to study the nonlinear energy cascade problem. The self-similar energy distributions already discussed in §1a are a strong candidate for isolated plate structures, but more general nonlinear analysis tools also exist. Established series methods include perturbation techniques as well as Volterra [12] and Weiner [13] series. Higher order terms in a Volterra series propagate the input forcing into frequency bands beyond those covered directly by the excitation force: this is exactly the scenario that results in energy cascades, and this approach is potentially useful for weak nonlinearity. However, the application of the Volterra series to a significantly nonlinear system may require a large number of terms in the series expansion, and convergence of the series is not guaranteed. By contrast, the Weiner series has better convergence properties (successive terms in the series are uncorrelated from all previous terms), but the approach is applicable only to white noise forcing, and hence it is not suitable for predicting energy cascades. Further, Weiner series are governed by kernels that must be recalculated for each forcing case, reducing the potential computational savings available through this method. Ideally, a technique where the descriptors of energy flow between frequencies are not dependent on the details of the forcing is desired.

Equally, a technique that does not depend on the detail properties and geometry of a structure has benefit. Experience with high-frequency linear vibration has shown that, at high frequencies, small physical variations in a structure at scales similar to the wavelength of vibration can significantly alter the linear mode shapes and hence the frequency response [2]. These variations can occur due to damage to a structure, but may also occur due to natural variability, for instance on production lines. The majority of the literature [6,8,10,14] presents theories of energy cascades for the response of a particular realization of a nonlinear structure. A technique that inherently handles the variation across different structures, possible by predicting the mean and variance across an ensemble of similar systems, would be valuable.

The mean response of an ensemble of structures is considered by Akolzin & Weaver [15] for the case of nonlinear stress–strain laws that when truncated to third order give rise to the same equations of motion as for Green’s strain approach. In [15], perturbation methods are used to evaluate the ensemble response for the early time evolution of energy flow in a diffuse field (uncorrelated signals in each mode). A system is initially driven by external forces to establish a diffuse field, and then released. The time evolution of the spectral response is considered in the presence of weak nonlinearity but in the absence of dissipation or further forcing . While a full damped, forced response is not considered, this is a step towards considering the variation in response across an ensemble of slightly random nonlinear structures.

For high-frequency linear vibration, an established technique, called statistical energy analysis (SEA) [2], is used to analyse the variation in behaviour across an ensemble of similar structures by providing predictions for the mean and variance of the energy of the system at a specified excitation frequency. As an example, a linear system consisting of a set of plates that are coupled by springs would be modelled with SEA using the following philosophy:

(i) The structure is divided into ‘subsystems’. In this case, each plate is considered to be a single subsystem.

(ii) Each subsystem is assigned an energy variable

*E*_{j}, describing the ensemble average total vibrational energy of that subsystem. For*N*subsystems, SEA is formulated as a set of power balance equations of the form 1.1where*ω*is the response frequency,*η*_{j}is the loss factor,*n*_{j}is the modal density, defined as the average number of modes per unit frequency,*η*_{jk}is the linear coupling loss factor (CLF) and*P*_{ex,j}is the external power input to the*j*th subsystem.(iii) Coupling loss factors

*η*_{jk}are calculated based on the plate properties and the strength of the coupling springs between the plates. Equation (1.1) represents a power balance at each frequency, between external power inputs, energy flows between subsystems and the loss of energy from the system through structural damping.(iv) Equation (1.1) is solved by simple matrix methods to yield the energies in the subsystems. Because of the assumptions involved in the detailed derivation of equation (1.1), the energies represent ensemble average values, where the ensemble consists of a set of structures with random properties.

As described above, the fundamental principle of SEA is one of a power balance between physical components of a structure. To extend this approach to the energy cascade scenario, a formulation of SEA is presented here where, rather than describing subsystems as separate plates in a built-up structure, each subsystem is a frequency band within a particular plate. Again each of these subsystems has an energy variable assigned, and nonlinear coupling loss factors (NCLFs) are derived to describe the energy flows between frequency bands. Importantly, these NCLFs are not a function of the forcing amplitude or bandwidth, depending only on the properties of the nonlinear structure itself. Section 2b provides the details of a benchmark model used to validate the results of the nonlinear SEA equations derived in §2c. Section 3 provides validation of the SEA model for the case of an isolated nonlinear system, and a nonlinear plate coupled to other linear subsystems.

## 2. Nonlinear statistical energy analysis

### (a) Governing equations

The equations of motion that govern a nonlinear multi-degree-of-freedom system are often expressed in terms of a set of oscillators with nonlinear couplings up to third order, which physically might represent nonlinearity arising from the large amplitude Green’s strain phenomenon [3], or material [15] or other nonlinearity. The derivation in [6] gives a detailed presentation of the origin of equation (2.1); however, it in general arises when a quartic and cubic energy expression is written in terms of a summation of linear oscillators. Each oscillator can be identified with a basis function that contributes to the total displacement of the plate; it is convenient to consider the basis functions to be the set of mode shapes that arise in the absence of any nonlinearity. The degree of freedom associated with oscillator *j* is then *a*_{j}, say, which corresponds to the displacement amplitude of the *j*th mass normalized mode shape. This approach results in a set of equations of motion of the form
2.1where *β*_{j}, *ω*_{j} and *F*_{j} are the natural frequency, damping factor and generalized forcing of the *j*th oscillator, respectively, and **C** and **D** are tensors that describe the strength of the second- and third-order coupling, respectively, between the oscillators.

### (b) Numerical benchmark

To validate the theory to be developed in §2c, an idealized benchmark has been developed consisting of simply supported linear flat plates (figure 2) coupled in series with linear springs. The nonlinearity is introduced by adding a set of randomly located springs with point couplings between the out-of-plane degrees of freedom of the plate and fixed ground locations. This model may be further randomized with the addition of point masses to the coupled plates. The coupled structure is then driven with a random point force. While this system is an arbitrary creation, it will be shown that it conveniently illustrates the energy cascade phenomenon described in §1a.

This idealized coupled plate model is governed by the nonlinearly coupled differential equations given in equation (2.1) which must be solved in the time domain, then post-processed to acquire a frequency domain description of the total energy in the system. A direct time domain simulation of these equations is computationally expensive for a system containing hundreds of linear modes, in particular the evaluation of the final term of equation (2.1) that describes the cubic nonlinearity. Using the coupled plate benchmark with nonlinear springs helps to limit the computation cost as time domain simulations now only require the calculation of the nonlinear forcing in each spring, rather than the full evaluation of the matrix products on the right-hand side of equation (2.1).

The springs have been chosen to produce a nonlinear force *F*_{n} of the form
2.2for a spring stiffness *k*_{nl}, offset *b*_{nl} and displacement *x*. In generalized modal coordinates this produces nonlinear terms that fully populate the tensors **C** and **D** given on the right-hand side of equation (2.1).

A further benefit of the random spring model is that it allows an ensemble of structures with different natural frequencies to be generated, by randomizing the spring locations and adding random point masses to the plates. The global strength of nonlinearity will remain approximately constant across the ensemble if the spring parameters *k*_{nl} and *b*_{nl} are kept constant and a sufficient number of springs are employed. Example cases from the benchmark simulation will now be presented to further illustrate the nonlinear energy cascade phenomenon.

For all the results presented in this paper, the following parameters are selected. Thirty nonlinear springs are used with the stiffnesses chosen as *k*_{nl}=10^{8} Nm^{−3} and offsets *b*_{nl}=1×10^{−3} Nm^{−2}. Where relevant, 10 linear coupling springs are used between each pair of plates and have a stiffness of *k*_{l}=10^{2} Nm^{−1}. The plates have identical material properties with flexural rigidity *D*_{r}=1.14 Nm, density *ρ*=8100 kg m^{−3} and thickness *h*=0.5 mm. The length along one edge is *L*_{1}=0.6 m, for all plates, and the second dimension *L*_{2} is 0.5, 0.6 and 0.7 m, respectively, for plates 1, 2 and 3. The *Q*-factor of the plates varies linearly with frequency *ω* and constant plate modal density *n*_{j}, following a distribution
2.3which has been chosen because a linearly increasing *Q*-factor with mode count matches that found in the experimental structure described in §1a. This distribution of the *Q*-factor gives a typical modal overlap factor (defined as *ωn*_{j}*η*, with loss factor *η*) of 2–3%, which, while low for typical SEA applications [16] (and hence requiring larger ensembles for convergence to a mean), keeps the damping at higher frequencies low and encourages nonlinear energy scattering into higher frequencies, highlighting the nonlinear cascade phenomenon. Figure 3 illustrates an example output of the benchmark code for case (*a*), in figure 2. A random forcing is applied either with frequency content across a narrow frequency band 100–110 Hz or from excitation from 0 to 100 Hz. The plate contains 100 modes. For clarity the output from a single realization of the structure is shown in figure 3*a*,*b* and then a full ensemble in grey in figure 3*c* with the ensemble average energy in black. Note that all decibel scales throughout this paper have a reference level of 1 J. All plots show the frequency content of the total vibrational energy (calculated as twice the time-averaged kinetic energy of the structure).

A few points of particular interest are worth highlighting. For the case of narrow band forcing, the frequency content of the response contains dominant bands at integer multiples of the forcing frequency band, so between 220–280 Hz and 330–420 Hz and so on. For case (*b*) in figure 2, the response contains energy across the full frequency range above the driving frequencies. It is worth highlighting that energy cascades are most apparent when the energy is scattered into frequencies which themselves do not have an external power input, as external inputs swamp the effect of the nonlinearity.

Further, it is apparent from the plots of the ensembles the degree of variability that can occur across a range of similar systems. This effect is well known for manufacturing imperfections in linear structures, when the variations between structures are of a similar size to the wavelength of vibration [2]. The effect is even more pronounced here as, in addition, the details of the nonlinear coupling strength between any set of frequencies is partly dependent on the details of the natural frequencies for a given ensemble member. This provides further support for an SEA-type theory as it provides a model of the mean response rather than the detailed response of an individual ensemble member.

### (c) Nonlinear statistical energy analysis equations

To describe the response of nonlinear systems showing energy cascades, this section will provide a formulation of a set of nonlinear power flow equations that model the coupling between frequency bands of a given structure. The equations will be derived in detail for the quadratic case only; however, the derivation is analogous for the cubic nonlinearity. Consider the equations of motion, equation (2.1), truncated to second-order terms,
2.4To derive an energy balance equation analogous to SEA, initially an assumption will be made that the nonlinear term on the right-hand side of equation (2.4) is an external force acting on a linear system. Linear vibration theory [17] may then be used to describe the response of the linear structure consisting of the oscillators *a*_{j} to a random external force *F*_{j}. For a linear structure with randomly distributed natural frequencies *ω*_{j}, undergoing random excitation, it is common [18] to approximate the zero mean, random process of each oscillator *a*_{j} as being statistically uncorrelated, defined as *E*[*a*_{i}*a*_{j}]=0, where *E* indicates the ensemble average across the different realizations of the random structure. Then for random forcing the spectral response of oscillator *j* can be written as
2.5and
2.6where *S*_{k} is the spectral response of the *k*th oscillator and *H*_{j} is the linear transfer function between force and displacement for the *j*th oscillator if all the nonlinear parameters *C* in equation (2.4) were zero. Solving for the mean-square response of the *j*th oscillator by integrating over *ω* and assuming that the forcing is flat near the resonant peak of each oscillator so a white noise approximation is valid for *S*_{FF} [17] gives
2.7To reduce the number of degrees of freedom, the oscillators are grouped into neighbouring frequency bands, each of which will be termed a subsystem. For instance, the first 10 oscillators could be placed in one band, the next 10 in the next band and so on. Note that there is no requirement that the bands all contain the same number of oscillators. The energy of the individual oscillators in a given band summed provides the energy of a frequency band, and hence equation (2.7) becomes for oscillators *j* in band *r*:
2.8where *Q*_{r} is the range of mode indices associated with subsystem *r* (and similarly for *Q*_{s} etc.). To complete the derivation, the spectrums *S*_{k},*S*_{l} on the right-hand side of equation (2.8) must be written in terms of the energies of the frequency bands. To achieve this, two further assumptions are required: (i) equipartition between the oscillators in an ensemble sense (i.e. averaged across an ensemble of structures the energy in nearby oscillators is approximately equal); (ii) the response is dominated by the excitation near the resonance, where the forcing is nearly flat (as already stated for the approximation of equation (2.5)).

Given these assumptions, the oscillator spectrum near resonance would be proportional to the linear transfer function squared, so for a band containing *n*_{e} oscillators,
2.9equation (2.8) then gives
2.10where
2.11and nonlinear coupling loss factors are given by the expression
2.12For a single subsystem, the *D*_{rs} contains only the system dissipation (case (*a*) in figure 2),
2.13To extend equation (2.10) to coupled subsystems, for instance the three-plate system (case (*b*) in figure 2), further assumptions are required regarding wave propagation through the nonlinear structures, and results will be drawn from linear SEA theory. This paper will only consider the theory for the simplest case of linear coupling between subsystems in different plates. While restrictive, as this does not allow for many practical applications (for instance, nonlinear impacts, or friction at joints), an extension of this work would be to consider how joint nonlinearity can transfer energy between frequency bands, in addition to nonlinearity internal to each plate.

As already described in §1b, linear SEA assigns a single energy variable to each subsystem (where each subsystem in figure 2*b* would be one of the plates), and the response is evaluated at each frequency individually. The linear coupling loss factors between subsystems (which will be introduced rigorously, shortly) would then populate the off-diagonal elements of matrix **D** in an equation of the form of equation (2.10), although without the nonlinear tensor **G**. For the case given in figure 2*b* in the absence of nonlinearity, **D** would be a three-by-three matrix, with one energy variable for each plate. Note the difference from the nonlinear model of equation (2.10) where a single plate has several energy variables, one associated with each frequency band.

To combine the linear SEA model (describing coupling between plates) and the nonlinear SEA model (describing nonlinear coupling between frequency bands in a single plate), it is necessary to describe each plate (linear or otherwise), with a set of energy variables associated with frequency bands. The energies *E*_{s} in equation (2.10) then include values for each frequency band in each plate, regardless of whether or not the plate is nonlinear.

The results for linear coupling loss factors are only valid for coupling between linear diffuse wave fields (i.e. one where there are many propagating waves in different directions) [2] in the neighbouring subsystems. To apply this to the nonlinear structure of figure 2*b*, it will be assumed that linear waves propagate through the structures from the linear couplings. The nonlinearity then acts within a plate and scatters energy into waves at other frequencies, but while scattered waves arrive at the linear coupling springs at frequencies in addition to those of the driving spectrum, those waves also propagate linearly.

Consider first the case of three linearly coupled plates, but with energy variables associated with frequency bands, rather than the whole plate as normal for linear SEA. In the absence of nonlinearity and if identical frequency bands are defined for each plate, then only couplings between bands with the same centre frequency would be expected, i.e. the first band in plate 1 would be coupled to only the first band in plate 2, and the first band in plate 2 would have coupling only to the first band in plates 1 and 3. This would yield a power balance of the form of equation (2.10), with *G*_{rsn}=0 and a matrix *D*_{rs} of size 3×3 m for *m* frequency bands in each plate.

If we maintain the assumption that all nonlinearity is internal to an individual plate, i.e. all coupling springs between plates are linear, then the addition of nonlinear springs to ground in a single plate would mean the elements of matrix *G*_{rsn} that correspond to coupling between frequency bands within the linear plate would all be zero. Elements of *G*_{rsn} corresponding to a nonlinear plate in the combined plate structure could be obtained from equation (2.10), as for the case of a single nonlinear plate.

For a set of coupled linear plates, the hybrid method [19] for linear vibration analysis provides analytical results for the linear coupling loss factors between bands of equal centre frequency in each structure, i.e. the first band in system 1 is coupled only to the first band in system 2.

These linear coupling loss factors *η*_{rs} are calculated in terms of modal density of the *r*th subsystem *n*_{r} from [19] as
2.14where **D**_{dir} is the direct field dynamic stiffness matrix, i.e. the dynamic stiffness of the subsystem with no reflections, **D**_{d} is the dynamic stiffness associated with the linear coupling springs and **D**_{tot} is the dynamic stiffness matrix including the spring dynamics and is of the form . A detailed understanding of the derivation of the parameters can be found in [19], but the results necessary to generate the coupling loss factors are given in equations (2.15)–(2.17). For a flat plate of uniform thickness *h*, density *ρ* and plate flexural rigidity *D*_{r}, the direct field dynamic stiffness at a point coupling is given by
2.15For the case of two plates coupled by *r* linear springs, the matrices **D**_{dir} and **D**_{d} are of size 2*r*,
2.16aand
2.16band for linear spring stiffness *k*_{l},
2.17For a three-plate system, this procedure can be repeated to obtain the linear coupling loss factors between each pair of plates in isolation. For a system consisting of linear plates and springs, i.e. where the tensor *G* in equation (2.10) is populated only with zero entries, the matrix **D**_{rs} then consists of a combination of dissipation losses and energy transferred between plates through linear springs:
2.18where *η*_{rs}=0 for any pair of bands at different frequencies, and *ω*_{c} is the frequency of the centre of the band.

This result can now be extended to include the cubic nonlinearity by analogy to provide the complete version of equation (2.10): 2.19The result in the above equation, however, is not sufficient, as, in deriving the driving spectrum of equation (2.5), it has been assumed that all the nonlinear terms on the right-hand side act as power sources for the response, whereas they are actually a function of the response. To complete the nonlinear SEA derivation, it is necessary to enforce a power balance across the whole system, requiring that the total dissipation in the linear terms is equal to the power input from external forces.

Considering a single term in the summation of the quadratic products of energy in equation (2.19), for instance *E*_{r}*E*_{s}, this represents a forcing term on the system as a result of the presence of energy in bands *r* and *s*. It would be reasonable to assume that the energy supplied to band *n* from this term is provided in some fraction from the bands that generate the forcing *r*,*s*. To ensure a power balance, equation (2.19) is completed by removing terms *αE*_{r}*E*_{s} and (1−*α*)*E*_{r}*E*_{s} from the two bands *r* and *s* for some value of the parameter *α*.

In the absence of a compelling reason for bias, a default value of *α*=0.5 is adopted, and it is shown in §3a that the solution is insensitive to the precise value of this parameter. A similar approach can be taken for the cubic term, but removing energy in thirds from each band that contributes to the nonlinear term *E*_{s}*E*_{m}*E*_{n}. The following equation presents the results of this energy correction:
2.20If all the nonlinear terms in the above equation are summed across the equations *r*, then they contribute zero net power to the system, as required for a power balance between external power inputs and system dissipation.

### (d) Nonlinear coupling loss factors

Key to the nonlinear SEA equations given in equation (2.20) is the strength of the forcing on a given frequency band due to the nonlinearity generated by pairs, or triplets, of energies in other bands. The strength of the coupling is governed by the nonlinear coupling loss factors, described by tensors *G* and *H* in equation (2.20).

The properties of both the quadratic and cubic coupling loss factors are analogous and so only the quadratic will be discussed here as it is easier to visualize the behaviour; however, all the arguments given can be extended to the cubic case.

Physically, the coupling between frequency bands can be interpreted by considering each part of equation (2.12) in turn. Initially, consider the coefficients *C*, describing in some sense the raw strength of nonlinearity in the system. For the benchmark, these coefficients depend primarily on the spring stiffnesses chosen in equation (2.2). For a real structure with irregular geometry, the coefficients *C* would depend on the details of the irregularities [6], while for a system with material nonlinearity, it would depend on the details of the nonlinear stress–strain law [15]. The coefficients can be interpreted as simply increasing the strength of coupling across the frequency range; the same would be achieved by increasing the amplitude of the forcing—it results in a gross increase in the strength of nonlinearity.

The coupling loss factor also contains a part that depends on the transfer functions that would occur in the linear system if all the nonlinearity strength terms *C* were zero. This product consists of the convolution of pairs of these linear transfer functions, which will produce a frequency-dependent function that has peaks at the sum and difference frequencies of the linear transfer functions that form the input to the convolution. The result of the convolution is scaled by the term *H*_{j}, which is the linear transfer function corresponding to a mode in the driven band. It can therefore be stated that the nonlinear coupling loss factor will have peaks, at frequency combinations where the sum or difference frequencies of driving oscillators in bands fall near another underlying linear resonance.

Figure 4 shows the surface generated by evaluating the integrals for a single term in the summation of equation (2.12). Evaluation is for a single mode with frequency 230 Hz, in a system with 100 modes. Peaks are observed at frequencies where the sum or difference frequencies of the driving oscillators fall near the frequency of oscillator *n*. It is of note that these surfaces and equations resemble those of the second-order Volterra kernels [12,13]; however, they are derived from a different set of underlying assumptions and arise from an internal power balance rather than a series expansion.

Performing the summations in equation (2.12) across all the oscillators in the frequency bands, the nonlinear coupling loss factors required to solve the nonlinear SEA equations are obtained. An example is again shown in figure 4, which has the same general pattern as described for the individual oscillator but averaged across frequency bands.

When evaluating the coupling loss factors described by equation (2.12), a summation across the lowest frequency band with the lower edge at *ω*=0 results in a large value for the coupling loss factor, when comparing the first band and any other that is found to introduce numerical ill-conditioning. To resolve these numerical difficulties, the first band is made with a frequency range that contains only a small number of oscillators, and then is decoupled from the rest of the system by setting all terms in the nonlinear coupling surfaces that contain the first band to zero. The amount of energy contained in these first few oscillators is small compared with the whole structure, so any error due to the absence of this portion of the power input in the global energy balance is acceptable.

For all results presented in §3, it should be recognized that the first band has been decoupled (in a nonlinear fashion) from all other frequency bands, so will slightly overestimate the energy compared with any benchmark estimates.

## 3. Results and validation

### (a) Solution methodology

This section will present example solutions to the nonlinear SEA equations described in §2b along with comparisons with the benchmark method described in §2a. Examples will be given of different excitation regimes and amplitudes to highlight the range of cases which approach can be applied to.

The solution of equation (2.20) for a given case can be carried out using a suitable algebraic minimization technique, which constrains the energy to positive values, to ensure a physically meaningful solution. However, as equation (2.20) contains both quadratic and cubic terms in the energy, it allows for the possibility of multiple real solutions, not all of which may be stable.

To ensure a dynamically stable solution to these equations, which would be expected for steady vibration of a real structural component, instead the following equation will be solved:
2.21where for compactness the new nonlinear tensors ** J** and

**are the outcome of combining the three realizations of tensor**

*K***associated with quadratic nonlinearity and four realizations of tensors**

*G***associated with cubic nonlinearity in equation (2.20). The addition of the term d**

*H**E*

_{r}/d

*t*allows the equations to be solved using a standard time domain solver, for instance the Matlab ode45 algorithm, but once the solution stabilizes to a fixed energy, the time-varying terms will be zeros and give the same solution to the algebraic equation (2.20). An arbitrary initial condition can be used and the presence of the linear term proportional to the energy acts as a damping term and ensures that the solution settles to a constant energy in each frequency band. A typical solution time is several seconds to solve for a given case, once the nonlinear coupling loss factors have been evaluated. This compares with hours for a full Monte Carlo simulation of the time domain benchmark, for every individual forcing amplitude or bandwidth considered.

### (b) Single subsystems

The nonlinear SEA approach summarized in equation (3.1) has been benchmarked against full time domain simulations of the sets of coupled oscillators given in equation (2.1). Four example cases are presented in figure 5, with the single nonlinear plate divided into 12 frequency bands of width 25 Hz (and, as discussed in §2d, a narrow first band of width 5 Hz), each band containing approximately 10 linear modes. Cases (*a*) and (*b*) in figure 5 are for constant forcing, *S*_{ff}, for each of the first four bands, with values of 0.00035 m^{2} s^{−3} and 0.001 m^{2} s^{−3}, respectively. An ensemble of 150 members has been used to calculate the mean and two-sigma error bars (calculated as the mean±twice the standard deviation across the ensemble divided by the square root of the ensemble size) for the ensemble energy from the benchmark simulation. For a linear system this forcing energy would be expected in the first three bands only, but the nonlinear effect is to spread the energy over a range three to four times larger than the initial excitation and the nonlinear SEA technique correctly predicts the energy to within 3–4 dB over the full frequency range. This is within the error bars for cases (*a*) and (*b*), up to above approximately 200 Hz.

Cases (*c*) and (*d*) in figure 5 are the response for forcing only the fourth frequency band at approximately 90 Hz, with values for *S*_{ff} for this band of 0.0002 m^{2} s^{−3} and 0.0005 m^{2} s^{−3}, respectively. Again the response shows energy in bands outside that of the initial forcing, this time dominated by the bands at two and three times the forcing frequency, approximately 175 and 250 Hz. The prediction of energy in the bands at resonance (i.e. the directly driven bands, and those at integer multiples of the driving frequency) again is accurate to within a few decibels; however, two distinct errors are present: (i) the SEA methodology has a tendency to smear the energy across neighbouring frequency bands, for example in case (*c*) the benchmark predicts frequency content in the seventh band, while the SEA method makes a prediction that includes energy in both the seventh and eighth bands; (ii) outside of the resonant bands the prediction is less accurate. Errors outside of strongly driven bands are to be expected as the assumption of equation (2.9) required that each frequency band contains resonant oscillators. A resonant assumption would incorrectly predict the energy if the oscillators were in fact only weakly driven. Noticeably, the nonlinear SEA results fall within the error bars of the bench mark code for the bands at resonance.

At the end of the derivation in §2c, it was necessary to enforce a power balance between frequency bands, by artificially removing a fraction of power according to a parameter *α*. It was stated that a value of *α*=0.5 would be used to remove equal energy from each frequency band involved in any given energy flow. To illustrate that this assumption appears to be valid, case (*d*) in figure 5 is rerun using values for the parameter *α* of 0.3, 0.4 and 0.5. The results for different values of *α* are shown in figure 6 and they are all similar to within 1 decibel.

### (c) Coupled structures

The general plate system consisting of coupled substructures will be considered as an extension of the nonlinear SEA theory and is illustrated diagrammatically in figure 2*b*. The three plates are simply supported at their boundaries, having the same material properties but different sizes, as defined in §2b, giving rise to slightly different model densities in each plate. They are coupled through their main structure by linear springs, and any of these plates can then be made nonlinear with nonlinear springs to ground, as for the case of the isolated subsystem. The structure properties are again perturbed with randomized masses across the three plates. Each plate is divided into six frequency bands, of width 30 Hz (and an additional narrow first band of width 5 Hz). The forcing is in the first four frequency bands at level *S*_{ff}=0.00013 m^{2} s^{−3}, for case (i) in figure 2*b* and into the first four bands at *S*_{ff}=0.00025 m^{2} s^{−3} for case (ii).

Two particular cases of this general system will be highlighted in this section: (i) plate 1 nonlinear with plates 2 and 3 linear (figure 7); (ii) plate 2 nonlinear with plates 1 and 3 linear (figure 8). These two cases have been chosen to highlight particular characteristics of the propagation of energy through a built-up structure, owing to nonlinearity internalized to particular subsystems.

Figure 7 shows the response for cases (i) above. For case (i), there is nonlinearity in the driven plate, and a subsequent coupling of system energy into the higher frequencies, as for the isolated system in §3b. The linear spring couplings between the subsystems then result in transfer of energy to plates 2 and 3 at progressively lower energy levels, as would be predicted by linear theory. The absence of nonlinearity in these plates results in no further scattering. Figure 8 presents results for case (ii)—no nonlinearity is found in the driven plate; instead, energy is coupled in a linear fashion to plate 2, where the nonlinearity then scatters energy into the higher frequencies. The distribution in plate 2 above the driving frequencies is then coupled linearly onto the third plate and back to the first. Nonlinearity in the second plate therefore results in the situation where the plate with the greatest energy at higher frequencies is perhaps counterintuitively (from a linear perspective) not the plate being directly driven. It is not clear why the SEA results for case (i) do not provide a smooth curve with increasing frequency, tending to jump above and below the error bounds on the mean, while for case (ii), it fits well particularly for plates 1 and 2.

Of note is that the third plate energy distribution is similar in both cases (i) and (ii); however, as the amount of energy coupled into higher frequencies is dependent on the energy levels in the second plate for case (ii), less energy is seen above 100 Hz in plate 1 for case (ii) compared with case (i). For case (i), the nonlinearity was operating on the energy in the directly driven plate so there was more energy available to be transferred to the higher frequencies, and then spread through the structure. The energies outside the driven band in the first plate are orders of magnitude different between cases (i) and (ii), but for plates 2 and 3 there is less difference. It is notable that the quality of the fit gets worse for the plates furthest from the drive point (i.e. plate 3), and worse than you would expect from a purely linear analysis using SEA. This perhaps implies that the assumptions made regarding the interaction of the nonlinearity with the linear couplings in §2c are not completely valid.

From a design perspective, this introduces a challenge if the aim is to control the flow of energy through a structure. Depending on the objective, if the aim was to isolate the third plate from the higher frequency response, it does not matter which of the plates exhibits the nonlinearity, the energy follows a path encompassing the nonlinearity and the high-frequency content occurs regardless of the exact location of nonlinearity, i.e. the whole superstructure must be protected from nonlinearity to prevent this behaviour. However, if it is the higher frequency response of the directly driven plate that is of interest there is merit in ensuring the structural element represented by plate 1 remains linear, as nonlinearity in other plates has less influence as the energy levels the nonlinearity is operating on are lower.

## 4. Conclusion and further work

This paper has presented an efficient energy method for predicting the response of individual nonlinear structural elements showing energy cascades as well as the response of more complex coupled structures consisting of a mixture of linear and nonlinear subsystems coupled with linear springs. A time domain benchmark consisting of an idealized model of a set of coupled plates that illustrates the physics of interest has been designed and used to provide a validation of the SEA methodology. The benchmark and SEA methodologies have been shown to demonstrate the same qualitative behaviour for energy scattering in both the isolated and coupled plate structures. Detailed comparisons of the results show good agreement in many cases, although some loss of accuracy is notable for the nonlinear SEA method, in particular for the coupled nonlinear structures. This perhaps implies that the assumption that the linear and nonlinear elements can be considered in isolation may require further consideration.

The nonlinear SEA theory is open to further development, in particular prediction of the variance of the response across an ensemble of systems. A further extension is to consider the effect of nonlinear joints. Linear coupling loss factors can be derived from the consideration of propagating waves across joints, and it may be possible to formulate an equivalent result with joints that can generate waves at frequencies different from those of the incoming waves.

There is also the potential to consider the transient response of these structures to an impact using SEA methods analogous to linear SEA.

## Authors' contributions

Both authors contributed to and helped draft the manuscript and developed the theory in §2c. The simulations were performed by G.M.S.

## Data accessibility

The Matlab scripts used to generate figures 3–8 are available on request: please contact G. M. Spelman (gms41{at}cam.ac.uk).

## Competing interests

We declare we have no competing interests.

## Funding

The authors would like to acknowledge the ‘Engineering Nonlinearity’ programme grant, funded by EPSRC, that has funded the research content of this paper.

## Footnotes

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

- Accepted June 16, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.