Analysis of the power flow in nonlinear oscillators driven by random excitation using the first Wiener kernel

Random excitation of mechanical systems occurs in a wide variety of structures and, in some applications, calculation of the power dissipated by such a system will be of interest. In this paper, using the Wiener series, a general methodology is developed for calculating the power dissipated by a general nonlinear multi-degree-of freedom oscillatory system excited by random Gaussian base motion of any spectrum. The Wiener series method is most commonly applied to systems with white noise inputs, but can be extended to encompass a general non-white input. From the extended series a simple expression for the power dissipated can be derived in terms of the first term, or kernel, of the series and the spectrum of the input. Calculation of the first kernel can be performed either via numerical simulations or from experimental data and a useful property of the kernel, namely that the integral over its frequency domain representation is proportional to the oscillating mass, is derived. The resulting equations offer a simple conceptual analysis of the power flow in nonlinear randomly excited systems and hence assist the design of any system where power dissipation is a consideration. The results are validated both numerically and experimentally using a base-excited cantilever beam with a nonlinear restoring force produced by magnets. © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Due to the diverse range of vibration sources encountered in engineering structures, a variety of different forms of excitation drive vibration. Random inputs constitute a large proportion of realistic excitation, although, depending on the application, the frequency content can range from narrow to broadband. Of particular interest in this paper are structures under stationary Gaussian random excitation of any spectrum.
In a number of applications the power dissipated by a vibrating system will be of interest. In some instances, such as vibration energy harvesting, the aim will be to maximise power, whereas in others, such as those involving concerns over fatigue or heat generation, the aim will be to minimise it. In addition, for single-degree-of-freedom (SDOF) systems with linear damping the power dissipated is proportional to the mean square velocity, a useful measure of the response of an oscillator.
Over the past two decades there has been great interest in vibration energy harvesting [1], the conversion of ambient vibrational energy into electrical energy in order to power small-scale electronic devices. Consequently, a large body of literature has been produced investigating the power dissipation of both linear and nonlinear oscillatory systems under various base excitations, although the theory is often applied to specific energy harvesting systems, see for example [2][3][4][5][6][7].
Of particular interest here are general methods for calculating power dissipation of nonlinear systems from random excitation. In general, the most prevalent technique is to solve the Fokker-Planck equation which governs the probability density function of the response [4,[8][9][10]. Whilst robust, this must be achieved computationally [11] or via simplifying the resulting equations by making assumptions. As the number of degrees of freedom increase, these solutions become significantly more involved. A noteworthy result in the case of power dissipated under white noise excitation is derived partially in Refs. [8,10,[12][13][14], and more generally in Ref. [6] and shows that for a general multi-degree-of-freedom (MDOF) nonlinear system subject to white noise base excitation, the power dissipated is simply proportional to the total oscillating mass and the magnitude of the noise spectrum regardless of the specific details of the system. This result is extended in Ref. [15] which uses the Wiener series to show that for systems exhibiting detailed balance the power dissipated under white noise excitation will be greater than or equal to the power dissipated under non-white excitation where the peak of the spectrum is taken as the magnitude of the white excitation.
The Wiener series is a useful method for analysing an output from a nonlinear system with a Gaussian white noise input via an orthogonal series expansion of the random output [16]. It is a commonly used tool for nonlinear system identification, particularly for physiological systems [17,18]. Whilst generally associated with white noise inputs, it can also be extended to non-white inputs [19] and this form, herein called the extended Wiener series, is applied in this paper. A thorough description and explanation of Wiener theory and its applications can be found in Ref. [16].
The aim of this paper is to provide a general methodology for calculating the power dissipated by a general nonlinear oscillator under non-white excitation. In what follows an introduction to the extended Wiener series is presented in Section 2 followed by the derivation of the method for calculating power dissipation in Section 3. The theory is then validated numerically and experimentally in Sections 4 and 5 respectively before conclusions are made in Section 6.

Wiener series for non-white excitation
In this section, the extended Wiener series for non-white input excitation is introduced. The series is very similar to the Wiener series for white noise and as such, the notation of [16] is used. A nonlinear system with a Gaussian random input, x(t), will produce a random output signal, y(t), that can be described as a sum of functionals When compared to the Wiener series of [16], the G-functionals, G n [k n ; x(t)], have been replaced with a lower case g n [k n ; x(t)] to represent that these are extended Wiener functionals for non-white noise. Each g-functional is defined as a sum of Volterra functionals, K j(n) [x(t)], up to the order of the g-functional such that where the Volterra functional K j(n) [x(t)] takes the form The order of the Volterra functional is j and k j(n) ( 1 , … , j ) is called an extended Wiener kernel of order j and must be calculated.
The (n) term in the subscript of both the functional and the kernel denotes that they both belong to the Volterra series of the nth order g-functional. When j = n the nth order kernel k n(n) will be rewritten as k n and is called the leading order kernel.
The relationship between the extended Wiener kernels in a g-functional can be found by enforcing the orthogonality condition where E [X] represents taking the ensemble average of the random variable X and H p [x(t)] is any Volterra functional of order p, where the form of a Volterra functional is given in Eq. (3). This condition is required to create an applicable orthogonal series, like the Wiener series, that converges and where the contributions from each g-functional can be isolated in order to calculate them. When the orthogonality condition is applied, the form of the g-functionals is attained and the first three and the nth gfunctionals are found as (8) where R xx ( ) is the autocorrelation function of the input at a time lag , represents time lags that are integrated over and in the nth functional the sum now includes only alternate order Volterra functionals. Equation (8) shows that each extended Wiener kernel within a g-functional can be derived from the leading order kernel and the autocorrelation function of the input and so the bracketed term in the subscript of the kernel in Eq. (3) can be dropped.
The g-functionals here are of a very similar form to the G-functionals, which can be retrieved if the autocorrelation function of white noise is applied, R xx ( ) = S 0 ( ). In Ref. [19] the first four functionals are given in this form and the nth functional can be assumed to follow the same pattern. For the sake of completeness the form of the nth functional is proven in the Appendix.

Calculating the first kernel
One of the key characteristics that makes the Wiener series desirable is that the kernels can be determined simply using cross-correlation [16]. It is important therefore to assess whether the kernels from the extended series can also be simply calculated. The cross-correlation method described in Ref. [16] is a time domain method where the cross-correlation between input and output measurements either from experiments or simulations are calculated at various relative time lags in order to calculate the Wiener kernels. A simple rearrangement of the method can convert it into the frequency domain such that crosscorrelations become cross-spectra. The frequency domain approach is more applicable for the extended Wiener kernels since the time domain approach for white noise relies on the useful white noise property that the autocorrelation function is a delta function, whereas the frequency approach does not.
In this paper only the first extended Wiener kernel is of interest and so the method for its calculation is briefly outlined as an extension of the method used for white noise in this section. The approach for the higher order kernels will be the same, but is not shown here. It can be noted that the zeroth order kernel is simply the ensemble average of the output, k 0 = E [y(t)] and will be assumed to be zero in this paper.
To calculate the first kernel, the cross-correlation of the output with the input at a time lag is taken and yields Since the time lag of the input can be interpreted as a first order Volterra functional, the orthogonality property of the kernels, Eq. (4), can be applied such that Eq. (9) becomes In the case of a white noise input, the autocorrelation function is simply R xx ( − 1 ) = S 0 ( − 1 ) and the cross-correlation of the input and output is equal to S 0 k 1 ( ). However, with a non-white input, the cross-correlation is less simple and Eq. (10) becomes less useful for calculating k 1 ( ). Equation (10) can be transformed to the frequency domain by taking the Fourier transform of both sides of the equation resulting in where S yx ( ) is the cross-spectrum between output and input and K 1 ( ) and S xx ( ) are Fourier transforms of k 1 ( ) and R xx ( ) using the Fourier transform convention The first extended Wiener kernel can therefore be calculated from Eq. (11) with known input spectrum and the cross-spectrum, S yx , calculated from either experimental or numerical results.

Calculating power from the extended Wiener series
The theory from the previous section can be applied in order to provide a framework for calculating the power dissipated under random excitation with a general spectrum. In order to derive a theory that encompasses as many systems as possible, this section analyses a general nonlinear N-degree-of-freedom-system identical to that of [15] of the form where y is a vector of the N generalised coordinates describing the system, ( ,̇, t) is a vector of length N that represents a general nonlinear dissipative and restoring force,b(t) is a random base acceleration and f(y) is a forcing vector that relates the base motion to the force on the system. The rth term of this vector can be calculated as ) (14) where m j is the jth mass, n is a unit vector in the direction of the base excitation in a Cartesian coordinate system, u j is the location of the jth mass in the same Cartesian coordinate system and y r is the rth generalised coordinate. The power, P, input by the base excitation and therefore dissipated by this system can be calculated by summing the effective force on each mass multiplied by its velocity such that where ( 16) and An extended Wiener series can be created for the random variable z(t) with input (t) such that and the power from Eq. (15) can be calculated with this series as . (18) This equation has the same form as Eq. (9) with = 0 and consequently where R ( ) is the autocorrelation of the excitation at time lag and the excitation has been assumed to have zero mean. This equation becomes more informative when converted into the frequency domain to yield Eq. (20) can be used to calculate the power where K 1 ( ) can be calculated using Eq. (11) either from numerical or experimental data. In order to understand how power will be dissipated, the properties of the first kernel are investigated in the following sections.

Properties of the first extended Wiener kernel
In what follows, a useful property of the first kernel, namely that influences understanding of power dissipation will be derived. Since k 1 ( ) and K 1 ( ) are Fourier transform pairs, Eq. (21) can be rewritten as and can be derived by following a similar process to convolution whereby the base excitation is considered as a series of impulses of magnitudes given by (t). A small change in z(t) due to an impulse (T) at time T, termed z(t)| (T) to differentiate it from the change in z due to previous impulses, can be assessed both from physical reasoning and from the Wiener series.
Physically, the effect of an impulse will instantaneously only change the acceleration term, ( )̈, and have no effect on the wherė| (T) t represents the accelerations of the generalised coordinates only due to the excitation at time T. This equation can be described physically as a unit impulse generating a unit change in momentum.
Rearranging Eq. (23), a small change iṅfrom the excitation,̇| (T) , is therefore The small time-step, t, here is from a time just before the impulse, T − , to a time just after the impulse, T + , where T + − T − = t and the magnitude of the impulse is (T) t. The change in the velocity vector from excitation (T) at time T,̇| (T) can also be calculated from the extended Wiener series. From Eq. (16) the change in the velocity vector provides a change in the z(t) variable where the Wiener series of Eq. (17) can be rewritten in a variational form and the Δg n terms represent the change in nth g-functional due to the excitation at time T, (T). The zeroth order Δg-functional is simply since it does not depend on the excitation. The first order Δg-functional is Since the instantaneous response to a single impulse of magnitude (T) t is being assessed, the integral of Eq. (28) disappears and the first Δg-functional becomes meaning that immediately after the impulse the first Δg-functional contributes k 1 (0 + ) (T) t to the response where the argument 0 + in k 1 ( ) represents the time instantaneously after = 0. This is analogous to the impulse response of a linear system where immediately after an impulse, the response would be given by the value of the impulse response just after the impulse multiplied by the magnitude of the impulse. However, whilst in the linear case k 1 ( ) would completely define the response, in the nonlinear case the higher order functionals may also contribute to this value so must also be assessed.
The second order Δg-functional is where the symmetry of the second kernel, k 2 ( 1 , 2 ) = k 2 ( 2 , 1 ), has been used and the subtraction of the double integral between − t/2 and t/2 on the first line accounts for the counting of the − t∕2 ≤ 1 ≤ t∕2, − t∕2 ≤ 2 ≤ t∕2 region in both of the first two terms on the right hand side. The value of this integral clearly goes to zero as t → 0 since it is of order t 2 . For all higher order Δg-functionals, this region will also be negligible since its value will be of order t n where n is the order of the functional. For the same reason as the first functional, one argument of the second kernel is evaluated at 1 where the limit of the summation is different from that of Eq. (8) because only terms that involve at least one excitation term can affect z.
Combining Eqs. (24)- (26) and (31) yields To simplify further and isolate the k 1 (0 + ) term, the (T) t terms can be cancelled and the expected value taken (noting that depending on the system both f and M(y) can be random variables) giving It is clear that for any even n, the nth term will be zero since it will include an ensemble average of an odd number of zero mean Gaussian excitation terms. The limit of the sum over m therefore becomes (n − 1)/2. The ensemble average containing n − 2m − 1 terms can be rewritten as a sum of products of autocorrelation functions. There are (n − 2m − 1)!/((n − 2m − 1)/2)!2 (n−2m−1)/2 distinct ways of pairing the terms so there will be this many terms in the sum and the integral over each term will be the same due to the symmetry property of the kernels. The sum in Eq. (33) therefore becomes where the integral is the same for every value of m. Analysis of this term shows that it is equal to zero since the sum becomes where M = (n − 1)/2 and a binomial expansion around the point −1 has been used. Eq. (33) therefore becomes The value of the first kernel, k 1 ( ), around the = 0 point will now be explored. When = 0 − , instantaneously before any excitation, the response must be zero because the system is causal so k 1 (0 − ) = 0. Combining with Eq. (36) suggests that at = 0 the value of the kernel is halfway between k 1 (0 − ) and k 1 (0 + ), so k 1 ∕2 as required by Eq. (22). This is equivalent to saying at k 1 (0), only half of the impulse from the excitation has occurred so the response is only half the magnitude. Whilst the simplicity of this result may seem surprising for nonlinear systems, it follows as a consequence of the power input being entirely dependent on the first kernel. The power input from any impulse of excitation is instantaneous and since it can only affect the first kernel, the value of the kernel at zero time is unaffected by any nonlinearity and therefore higher order kernels.
Although the triple product term E in Eq. (36) seems physically unintuitive, it has been discussed in detail in Ref. [15] where it has been shown that for a system with a total oscillating mass, M Tot , In Ref. [15] Langley shows that if a system is constrained to reduce its number of degrees of freedom then E However, if there is no constraint on the system, the inequality of Eq. (37) becomes an equality so Since the time domain kernel, k 1 ( ), is real, the frequency domain kernel, K 1 ( ), is Hermitian. The imaginary part therefore does not contribute to the integral in Eqs. (20) or (21) so the equations can be modified such that the integrals are only performed over positive frequencies and the real part of the kernel giving The real part of the first kernel completely defines the power dissipated.
To summarise the results so far: the power dissipated by a general nonlinear oscillator, Eq. (13), under random excitation of a general spectrum, S ( ), can be calculated using Eq. (39). To do this, the first Wiener kernel must be calculated from either simulations or experimentally using Eq. (11), but crucially for a designer of a system desiring a preliminary estimate of power dissipation, the first Wiener kernel has the property of Eq. (40) where the triple product term is simply equal to the oscillating mass provided the system is unconstrained.
It should be noted that this approach extends the results for white noise excitation of [6] and [15] to general random excitation. The white noise results are retrieved by using S ( ) = S 0 where the factor of remains in the spectrum due to the Fourier transform convention taken here. This is substituted into Eq. (39) and combining with Eqs. (40) and (37) the power dissipated by white noise is

Numerical validation
The results of Eqs. (39) and (40) are validated via numerical simulations in this section. Time domain simulations were conducted where an ensemble of results is built up by exciting a nonlinear oscillator with a number of realisations of random noise of a chosen spectrum. The ode45 routine in MATLAB was used and the first extended Wiener kernel can be found using Eq. (11) where S yx ( ) becomes S z ( ) and S xx ( ) becomes S ( ) with being the negative of the base motion.

Single-degree-of-freedom system
A simple example of a nonlinear system containing cubic damping and stiffness nonlinearities with equation of motion mÿ +̇y 3 + ky + y 3 = −mb(t) showing that the simulations agree well with the theory. The form of the kernels in Fig. 1 contains information about the dynamics of each system. As the stiffness nonlinearity increases, the peak frequency of the kernel increases because the system becomes stiffer. Additionally, the time domain first kernel decays more rapidly and the frequency domain kernel becomes lower and wider as nonlinearity increases, meaning that the kernel appears more damped with nonlinearity. This effect is due to the higher order kernels having greater influence at higher nonlinearity meaning more energy is distributed from the first kernel to the higher order ones thus increasing the energy loss, or damping, of the first kernel.

Multi-degree-of-freedom system
The numerical simulations show strong agreement with the theory for a SDOF oscillator so attention is now turned to a MDOF scenario where two coupled masses are free to move in one direction. The system, displayed in Fig. 2, has equations of motion In this case the vector triple product E

Experimental validation
To supplement the compelling numerical validation of Eqs. (39) and (40) in Section 4, experimental validation has been undertaken and is discussed in this section. Calculation of the power dissipation using Eq. (39) is experimentally difficult, and so only the property of the first kernel from Eq. (40) is examined. A base-excited cantilever with a tip mass as shown in Fig. 4 is used to approximate a SDOF oscillator. Neodymium cylinder magnets have been placed on the tip mass and base arms to generate a nonlinear stiffening restoring force. Clearly a cantilever beam is not a SDOF system, although with a tip mass the second resonance is found at a significantly higher frequency than the first such that the second mode has little effect on the first extended Wiener kernel.
The apparatus requires a single input, the excitation time signal to the shaker and two outputs, the base acceleration and tip velocity relative to the base. A computer with a National Instruments data-logging card was used to generate and measure the relevant signals, although due to the system dynamics, the input voltage time signal to the shaker did not correspond closely to the base acceleration signal. This is unimportant, however, as only the base acceleration spectrum is required for the calculation of the first kernel.
The base acceleration was measured using an accelerometer and, due to their strong performance at opposite ends of the frequency spectrum, both an accelerometer and a laser vibrometer were used to measure the relative tip velocity. The accelerometer was used to calculate the relative tip acceleration which can be found by subtracting the base from the tip acceleration and then integrating over time to provide velocity. The velocity calculated from the accelerometer data was poor at low frequencies, but strong at high frequencies due to time integration, therefore the relative tip velocity was also measured using a laser vibrometer which has stronger performance at low frequencies, but is noisier at high frequencies. To provide the smoothest data, the zero to 50 Hz part of the first kernel from the laser vibrometer data was combined with the kernel from the accelerometer data from 50 Hz upwards.
Since the experimental system is unconstrained and can be approximately modelled by a SDOF equation of motion similar to Eq. (42), the triple product, E , is equal to the oscillating mass. The value of this mass will be difficult to calculate as some of the beam mass must be included. It is therefore preferable to modify the output of the extended Wiener series from what would be z = mẏ according to Eq. (16) to z =̇y whereẏ is the relative tip velocity. The result is that the integral over the  kernel in Eq. (40) is no longer equal to m/2, but 1/2. This is easier to validate since no estimate of the oscillating mass is required, but is still validating the important property of the first kernel since its magnitude has been divided by its mass.
It is useful to examine the first extended Wiener kernel for the system with varying degrees of nonlinearity. A convenient method for achieving this is to vary the magnitude of the excitation in order to affect the magnitude of the tip motion between the magnets. With a low amplitude base motion, the tip motion is also low meaning its range of movement within the nonlinear magnetic force is low. Across any small displacement, the nonlinear force could be reasonably approximated by a linear force so the system is close to linear. However, when the base and therefore tip motion are large, the tip moves through a large region of nonlinear magnetic force and so the nonlinearity of the response is increased. Figure 5 displays the first extended Wiener kernel in the time and frequency domains with three different excitation magnitudes where each kernel was calculated from ten 30 s realisations. Similarly to the numerical results of Fig. 1 where the stiffness parameter is increased, the kernels show the natural frequency and the damping and therefore nonlinearity of the first kernel increasing with excitation amplitude. The critical result to verify is that the integral over the real part of the frequency kernel behaves according to Eq. (40). The values of this integral are 0.51, 0.51 and 0.50 for the low, medium and high magnitude excitation respectively providing excellent resemblance to the theory. It should be noted that any higher order resonance effects of the cantilever have a negligible effect on the integral of Eq. (40).
The spectrum of the base excitation at each three levels is displayed in Fig. 6 where the input is clearly non-white. In fact, it is evident from the dip at the resonant frequency in each spectrum, that the tip mass is acting as a vibration absorber.

Conclusions
A methodology has been presented for calculating the power dissipated by nonlinear MDOF systems under general random base excitation. The Wiener series approach for a Gaussian random input with general spectrum has been applied and the power dissipated is shown to be dependent only on the first kernel and the spectrum of the input according to Eq. (39). Calculation of the first kernel can be made either via simulations or from experimental data and it is shown to have the property that the integral over the frequency domain kernel is proportional to the total oscillating mass, Eq. (40), for an unconstrained system.
The combination of Eqs. (39) and (40) provides a simple conceptual understanding of power flow in nonlinear randomly excited systems. The first extended Wiener kernel displays peaks around the resonances of the system therefore for applications where power input to the system is to be minimised, the resonances of the system should be designed to be at a frequency where the input spectrum is low. Conversely, for applications such as energy harvesting where power is to be maximised, the resonances of the system should be tuned to frequencies where the input spectrum is high. Whilst these conclusions are conceptually obvious, the derivation of Eq. (39) provides a rigorous mathematical framework for calculating power and Eq. (40) provides valuable information concerning how much power can be dissipated.
A critical assumption of the theory is that the base motion is independent of the oscillating system, or no feedback is provided from the system to the base. In reality this will not be true, but may be a valid assumption when the oscillating mass is significantly smaller than the base mass.
Numerical simulations of two illustrative systems have been used to investigate the theory and found to closely match the values predicted by Eqs. (39) and (40). Additionally, a simple experiment involving a base-excited cantilever beam using magnets to generate a nonlinear restoring force on a tip mass was conducted and also produced strong results when compared to the values predicted by theory.
which will be shown to be equal to zero by extending the method of [16].
Inspecting the ensemble average term in Eq. (A.1), it can be rewritten as where q = n − 2m, s q represents x(t − q ) and l p represents x(t − p ) and are Gaussian random variables. Using the Gaussian property (see for example [20]), the ensemble average can be expanded into a sum of products of ensemble averages of one or two variables. Since the input is assumed to be zero mean if p + n is odd then the ensemble average in Eq. (A.1) becomes zero and the orthogonality condition is proven trivially. If p + n is even then the ensemble average becomes a sum of products of autocorrelations of the input. For example Due to the symmetry property of the kernel in Eq. (A.1) various terms in the sum of products of autocorrelations will provide the same result when integrated over both the and spaces. For example, the rightmost two terms in Eq. (A.3) will provide the same result from the integral because 1 and 2 are interchangeable whereas the first term will provide a different result. In fact, for a given value of m the same integral will be produced by any product of autocorrelations that has the same number of autocorrelations with arguments involving the difference between a and term such as  Note that N r|m depends on m since q depends on m as q = n − 2m and if q is odd (then p is odd) and r must be odd and if q is even (then p is even) and r must be even. Additionally, since the maximum number of s and s that can be taken out and paired is the minimum of q or p, r must be less than or equal to the smaller of q and p so that To prove that the series is orthogonal, this equation must be equal to zero. For a given value of r, the sum over m will contain the integrals of the same value, having the same number of autocorrelations that have both a and a term in their argument. However, the upper limit of the sum over m must change according to the value of r selected since according to Eq. (A.5), r ≤ n − 2m so for a given r, m ≤ (n − r)∕2. Therefore to show that Eq. (A.6) is equal to zero and hence prove orthogonality, it must be shown that for any permissible value of r