Time History Forced Response in Nonlinear Mechanical Systems

A formulation of a digital filter method for computing the forced response of a linear MDOF mechanical system is proposed. It is shown how aliasing error effects can be avoided at the expense of a bias error. The bias error is however completely known and it is system independent, as it only depends on the sampling frequency used. The mechanical system is described by its modal parameters, poles and residues. The method is extended to include non-linear elements. A toolbox in MATLAB has been created where nonlinear elements with and without memory can be treated, as well as system described by coupled non-linear equations.


Introduction
There exist several methods to obtain the time history response of a mechanical system to an arbitrary excitation. A few of the time domain methods are: the Duhamel integral, the state transition matrix method, and the Runge-Kutta method with variations. The Duhamel integral, also known as the convolution integral, is a wellestablished method in the literature, its main drawback is the computational burden associated with its nonrecursive nature. However, the integral can be used to derive a recursive algorithm, a digital filter, which leads to a faster solution than the naive Duhamel integral method. In the paper, the digital filter method is introduced on linear systems. The method is then extended to nonlinear systems with localized nonlinearities.

Establishing a Recursive Solution
Working on discrete sampled data necessitates transformation of a model in the continuous time domain into the discrete time domain. The basis for a family of methods is the convolution integral: where x(t) is the input signal, h(t) is the system impulse response and y(t) is the output signal. The convolution integral may be used for any linear, time invariant system. We start by studying the simple analog system: where s is the Laplace variable. The system in equation (2) has the impulse response: We now calculate y(nT+T), where T is the sampling interval: From equation (4) we find that we may calculate y(nT+T) with a recursion formula using only y(nT) and the input signal in the time interval [nT, nT+T]. The system influence shows in the terms exp(-aT) and exp(au) in the last integral in equation (4). The way we use samples of x(t) to evaluate the last integral in equation (4) defines the design method, see figure 1. A recursion formula may be thought of as a digital filter.
As an example, the step invariant method uses the constant value x(nT) in the whole interval, which, using equation (4), gives:  A difference equation, such as equation (5), may be regarded as a digital filter. A general expression would be:

Application to Linear Mechanical Systems
There are various ways to represent a mechanical system: through a Finite Element Model (FEA) model, a lumped MCK model or an analytical model. All of them can provide quantitative information about the modal parameters, i.e. poles and residues of the system. Here we note that these parameters can as well be obtained by performing an experimental modal analysis on the mechanical system. See figure 2.

Fig. 2.
A linear mechanical system may be characterized by its residues and poles. These may come from different models of the system, or from experiments. The residues and poles then define the corresponding filters.
According to the modal superposition theorem for a MDOF system, consisting of N modes, the frequency response function H(s) connecting one input degree-offreedom with one output degree-of-freedom can be expressed using partial fraction expansion in terms of the residues R r and poles s r as follows): For the impulse invariant transform, equation (4) gives: For a complex conjugate pair with residue R r and pole s r , as in equation (7), we will get: Equation (9)  To get the complete response, the contributions from all wanted modes are just added.

Error Calculations
As in all sampled systems, we will have aliasing. It doesn't help us if our input signal x(t) has aliasing protection, as the system response, h(t), in general has not. As a result, the spectrum of our response will have contributions from the original spectrum Y(f) centred 03002-p.2 around multiples of the sampling frequency f s . The aliased spectrum Y a (f) will be: The error resulting from aliasing is system dependent. If the system frequency response has large amplitude at the sampling frequency, the error in the low frequency region may be considerable, due to the fold back. This error when we use the impulse invariant method is well known and has led to the recommendation not to use impulse invariant technique when handling response to SDOF systems, for instance when calculating shock response spectra. In an international standard, ISO 18431-4, the ramp invariant method is recommended instead. It is not so well known, however, why other methods lead to better result. The handling of the input x(t) in the last integral in equation (4) and figure 1 may be seen as a convolution between the input samples and a convolution kernel. In the case of ramp invariant, the convolution kernel w(τ) is: which is a triangle with base from -T to +T. As we have a convolution in the time domain, this corresponds to a multiplication in the frequency domain. The multiplier is the Fourier transform of the convolution kernel. For the ramp invariant the amplitude of the Fourier transform is: This function is zero at the sampling frequency and quite flat in the frequency region around, which helps to reduce the aliasing error considerably. The price to pay is a bias error in the low frequency domain, but it is known and starts at zero for low frequencies. The aliasing error is dependent on the system and the chosen sampling frequency, while the bias error (for a certain frequency) is only depending on the sampling frequency. This makes the ramp invariant version the preferred one in many applications, such as shock response spectrum calculations. Another benefit is that (as opposed to step invariant) there is no phase error for the ramp invariant. This is very important when we study nonlinear systems.

Forced Response in Nonlinear Systems
As a simple case, we consider a non-linear MDOF system for which both the applied force and the non-linearity are at the same degree of freedom as shown in Figure 3.
More specifically, the single non-linear component is located between the input force DOF and ground.

Fig. 3. A non-linear component is introduced to a linear system.
The input force and the non-linearity are applied at the same degree-of-freedom.
The equation of motion for the system of interest in discrete time domain is governed by the following equation: with the usual notations and where g(n) represents a single non-linear component which is a function of x. Equation (13) in the frequency domain (z-plane) is conveniently written as: r r r r r Equation (16) is identified as a recursive equation which can be written down for each mode r (i.e. for each digital filter). Using the superposition theorem then, the total response is the sum of each of the individual digital filters. Taking a step further, Equation (16) can be re-CSNDD 2012 03002-p.3 written in a compact fashion for each mode r and at time n as follows: where it can be seen that the term Q r embodies the terms in the second line of equation (16) and also the term . In this respect we note that Q r is completely determined by past parameter values together with the current applied force f(n). Now, equation (17) can be recognized as being a non-linear equation in x with g being some non-linear function of x. Equation (17) is computed for each mode r and then summed up for all the modes; after which the resulting equation is solved to yield the desired total response at time step n. The procedure is summarized in figure 4.   Fig. 4. Illustration of the forced response simulation method. The underlying linear system is described by its residues and poles and the nonlinearities are modeled as additional external forces.

Extensions
The idea described here has been extended in many ways, and a toolbox with simulation functions has been created in MATLAB. The functions have been extended to incorporate many non-linearities, coupled between different degrees of freedom for the linear system. This means that we have to solve a non-linear system of equations for each time step.
Very efficient routines to accomplish that have been implemented into the scripts. In order to make the execution fast, the non-linearities have to be defined within the script.
As an example, the function nonlindxmck has the following syntax: x = nonlindxmck (F, fs, M, C, K, FD, nlD, modes, resps) The input is the force vector F with a sampling frequency of fs. The linear system is in this case defined by the matrices M, C and K. The force is applied to degree of freedom FD, and the non-linearities are applied to degrees of freedom, DOFs, defined by nlD. Modes is a vector containing the modes that we want in the calculation. Finally, resps is a vector defining in which degrees of freedom we want the response(s).
As an example, let the system be a simple four degree of freedom linear system. We apply one cubic nonlinearity from one DOF to ground and we have another cubic nonlinearity between two other DOFs. We then apply a random force to the remaining DOF and calculate all four responses. See figure 5. To calculate all the responses in 100,000 samples takes 11 seconds on a laptop computer.

Summary
In this paper, a formulation of a digital filter method for computing the forced response of a linear MDOF mechanical system has been proposed. It is shown how aliasing (and phase) error effects can be avoided at the expense of a bias error. The bias error is however completely known and it is system independent, as it only depends on the sampling frequency used. The mechanical system is described by its modal parameters, poles and residues. The method is extended to introduce non-linear elements. A toolbox in MATLAB has been created where nonlinear elements with and without memory can be treated, as well as system described by coupled nonlinear equations. The functions are very effective and the computation time is minimized.