Abstract
Nonlinear mechanical vibrations are widely found in engineering structures. Traditional linear analysis methods are failing to meet the growing requirements of analysis accuracy. Because of the complexity of nonlinear system, there is no general method for nonlinear system analysis. Based on the describing function method, a polynomial inversion method may be used in nonlinear system identification. Using this method, the measured DF can be inverted into the form of force which has a clearer physical meaning and is easier to use for parameter fit. A MDOF system with two types of nonlinearities at the same location is studied to demonstrate this method.
1. Introduction
In the past, most attention was paid to linear systems, and theories and methods have been developed for linear system modeling, simulation and prediction. However, there exist lots of nonlinearities in an actual engineering system and some of them may have a significant effect in system dynamic characteristics. The study of nonlinearities is indispensable and inevitable. However, one should be forced to admit that there are no general analysis methods that can be applied to all nonlinear system in all instances [1].
During the past decades, lots of work has been done in identification of system nonlinearities and they can be roughly divided into two categories, time domain methods and frequency domain methods. Time domain identification of nonlinear structural models relies on processing time series while the data processing in frequency domain methods can take form of Fourier spectra, frequency response and transmissibility functions, or power spectral densities [2]. The describing function, DF, based on the definitions of the frequency response function and the firstorder harmonic balance method, was used in analyzing nonlinear structures by Tanrikulu O., et al. [3]. A further research was made by Özer and Özgüven [4] in calculating the describing function’s approximate solutions by using the ShermanMorrison Method. Aykan and et al [5] proposed a method to transform the describing function into the restoring force to make it more visualized.
Finding the nonlinear forces corresponding to the measured describing functions is called describing function inversion [6]. In this paper, nonlinear system identification using describing function is introduced briefly. A polynomial inversion method of inverting the measured describing function is proposed. At last, an illustrative example is given to demonstrate this method.
2. Theory
2.1. Describing function method
The use of nonlinear describing function matrix to represent the nonlinear forces in a MDOF system was proposed and developed by Özer et al. [4]. They established a set of feasible methods for nonlinear system identification, including early detection, nonlinear localization and parameter estimation. This method was confirmed to have an acceptable application to a MDOF system with several nonlinearities at different locations in a case study [7]. As this method was discussed in detail in ref. [7], a brief introduction is given below for readers’ comprehension.
The equation of motion for a nonlinear MDOF system under excitation can be written as:
where $\left[M\right]$, $\left[C\right]$, and $\left[K\right]$ are mass, viscous damping and stiffness matrices of the linear system, respectively. $\left\{f\right\}$ and $\left\{x\right\}$ stand for the external forces and the responses of the system. $\left\{nlf\left(x,\dot{x}\right)\right\}$ represents the nonlinear internal forces in the system, and it is a function of the displacement and/or velocity responses, depending on the type of nonlinearity present in the system.
A harmonic excitation with circular frequency $\omega $ may be written:
Using first order harmonic balance method, the response is assumed to consist of only the first fundamental harmonic, disregarding higher order harmonics. Since the response is harmonic as the assume, $\left\{x\right\}=\left\{X\right\}{e}^{j\omega t}$, the internal nonlinear forces can be harmonic too:
where $[\u2206\left(x,\dot{x}\right)]$ is the ‘nonlinearity matrix’ that is dependent on response amplitude. The “nonlinearity matrix” can be formed by using describing function $v$ representations of internal nonlinear forces as follows:
The relationship between describing functions and nonlinear forces will be introduced in Section 2.2. From the equations above:
Then the nonlinear frequency response function (FRF) matrix $\left[{H}^{N}\right]$ can be calculated as:
The FRF of underlying linear system, $\left[{H}^{L}\right]$, can be written as:
Notice that the nonlinear FRF may be different at different excitation levels due to the nonlinearity matrix. The nonlinearity matrix can formally be obtained by:
Eq. (9) is the core of this method. Once we get the nonlinearity matrix, we have the ability to localize and identify the nonlinearities in a MODF system. Although there still exist some challenges, such as the calculation of $\left[{H}^{L}\right]$ and $\left[{H}^{N}\right]$, they were discussed and partially solved in previous studies mentioned above. Readers who are interested can have a more comprehensive understanding in ref. [4, 7].
2.2. Polynomial inversion method
By the use of DF method, the nonlinear force describing function of the response’s amplitude can be obtained. Generally, they are represented in the form of two related series, the response amplitude, $\left\{dfx\right\}$, and corresponding DF values, $\left\{df\right\}$. The describing function, which can be used for determining the type of nonlinearity as well as for parametric identification, has a quite complex form and it will be worse if there exist multiple nonlinearities in one location of the nonlinear system. This not only makes it necessary for a wealth of engineering experience but also greatly increases the difficulty for identification. An inversion method that works for general describing functions will effectively solve this problem. In the rest of this section, a polynomial inversion method will be given in detail.
The describing function for general dynamic system nonlinearity can be written as:
where the response $x$ is supposed to be:
Eq. (10) indicates that the describing function is a complex function, and a complex function can be divided into two parts, real part and imaginary part:
The real part ${v}_{p}$ corresponds to stiffness nonlinearities while the imaginary part $j{\bullet v}_{q}$ corresponds to damping nonlinearities [8]. For each part, the derivation for the inverse is given as follows:
An arbitrary nonlinearity can be expressed as a sum of n nonlinearities, ${N}_{1}$, ${N}_{2}$,..., ${N}_{n}$:
The DF for composite nonlinearities is given by:
Eq. (14) and Eq. (15) indicate that the DF and the inverse DF transformations satisfy the superposition principle. By polynomial fitting, an arbitrary function can be approximately represented in the form of a sum of polynomials. Theoretically, the higher the polynomial order we use, the more accurate result we will get.
Because the DF yields zero when the corresponding force is even, we represent the inputoutput characteristic for a nonlinearity as an odd nthorder polynomial:
The $m$th item of DF can be computed by integrating over the interval $0<\tau <\pi /2$, as:
$=\left\{\begin{array}{c}{c}_{m}\frac{m\left(m2\right)\left(m4\right)\cdots \left(3\right)}{\left(m+1\right)\left(m1\right)\left(m3\right)\cdots \left(4\right)}{A}^{m1},m>1\text{i}\text{s}\text{}\text{a}\text{n}\text{}\text{o}\text{d}\text{d}\text{}\text{i}\text{n}\text{t}\text{e}\text{g}\text{e}\text{r},\\ {c}_{m}\frac{4}{\pi}\frac{m\left(m2\right)\left(m4\right)\cdots \left(2\right)}{\left(m+1\right)\left(m1\right)\left(m3\right)\cdots \left(3\right)}{A}^{m1},m>0\text{i}\text{s}\text{}\text{a}\text{n}\text{}\text{e}\text{v}\text{e}\text{n}\text{}\text{i}\text{n}\text{t}\text{e}\text{g}\text{e}\text{r}.\end{array}\right.$
Thus, the DF is written as:
The polynomial inversion method is composed by the following major steps: 1) Get the approximate describing function by the polynomial fitting of DF values calculated; 2) Transform the polynomial coefficients by Eq. (17) respectively and get the approximate nonlinear force function (NLFa); 3) Plot the NLFa, make a guess of the type of nonlinearity and give an initial value of its parameters from the graphic information (taking the saturation for example, $NLF=nlf\left(k,M\right)$, $k$ is the gain of linear part and $M$ is the saturation value).
In order to improve the accuracy of identification, the $fm$ in search function of MATLAB is used. Calculate the DF from the NLF we get in Step 3. Comparing this DF with the original DF we can finally get the parameters’ values with accuracy through iterations. An illustrative example is given in the next section.
3. An illustrative example
A four DOFs system with two nonlinear elements in the same DOF is studied here.
Fig. 14 DOFs nonlinear system used in case
System parameters:
$M=$ [1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0; 0, 0, 0, 5];
$C=$ [10, –5, 0, 0; –5, 10, –5, 0; 0, –5, 10, –5; 0, 0, –5, 10];
$K=$ [1000, –500, 0, 0; –500, 1000, –500, 0; 0, –500, 1000, –500; 0, 0, –500, 1000].
The nonlinear elements between m1 and ground are ‘cubic stiffness’, ${nlf}_{1}$ and ‘gap with end stiffness’ ${nlf}_{2}$:
The sinusoidal excitation is applied at DOF1. The system is linear under small excitation, while the nonlinearities appear when the excitation gets big enough.
The simulation tests are performed using harmonic force with frequency 2 Hz, Fig. 2. The FRFs obtained from 0.1 N amplitude harmonic force are taken as linear FRFs of system. For the nonlinear FRFs the force amplitude is stepped up to10 N. The simulation is performed with functions in a nonlinear MATLAB toolbox [9].
Fig. 2Linear and Nonlinear FRF H11
Fig. 3Sum of NLN of different DOFs
Fig. 4Comparison of calculated and fitted DF
Fig. 5Comparison of theoretical and fitted nlf
From Fig. 3 it is found that the nonlinearities only exist between DOF1 and ground. By using measured FRF the describing function of system nonlinearities can be calculated by Eq. (9) and is shown in Fig. 4 (black). By observing, a turning point can obviously be found at the displacement 0.008 [m]. Polynomial fitting in interval [0, 0.008] and [0.008, $+\infty $] is made with degrees 3 and 5, respectively. The polynomial fitted DF curve is obtained and shown in Fig. 4 (red). For the convenience of display, all values of displacement are multiplied by 1e+3. The polynomial coefficients, ${[c}_{n},{c}_{n1},\cdots ,{c}_{1}]$, for each part are:
Dividing ${coef}_{df1}$ and ${coef}_{df2}$ by the corresponding transform coefficients from Eq. (17), the polynomial coefficients for nonlinear force expressions are found. The comparison of fitted (red) and theoretical (black) nonlinear force functions is shown in Fig. 5:
By observing Fig. 5, when $0<x<0.008$, the expression of force should be in the form of $nlf\left(x\right)={a}_{3}\times {x}^{3}$ and the corresponding theoretical describing function, with parameter ${a}_{3}$, is $df\left(A\right)={a}_{3}\times {0.75\times A}^{2}$. Setting the difference between DF with ${a}_{3}$ and DF calculated as the FUN in $fm$ in search (FUN, ${a}_{3}$) of MATLAB and giving ${a}_{3}$ an initial value 1.706$e5$, we can get the value of ${a}_{3}$ with accuracy, ${a}_{3}=$ 169687. The change of parameter ${a}_{3}$ with the number of iterations is shown in Fig. 6. As it’s a pure cubic stiffness in the first interval, the nonlinearities could be a combination of a cubic stiffness and a nonlinearity with gap. Compare the values of fitted $nlf$ in Fig. 5 and $169687\times {x}^{3}$$\left(x\ge 0.008\right)$. The guess for the nonlinear force is $nlf\left(x\right)=169687\times {x}^{3}+k\left(x0.008\right)$. Similar to parameter ${a}_{3}$, we do the parameter fitting using fminsearch. Give an initial value $1e3\bullet 2.0955$ and parameter $k$ stays at the value $k=$ 228.0690 after iterations.
Fig. 6Iteration of parameter a3
Fig. 7Iteration of parameter k
4. Conclusions
Identification of nonlinear system is considered by the inversion of describing functions. The method in this paper reduces the difficulty but improves the accuracy of nonlinear system identification. Thanks to the broad applicability of polynomial fitting, the polynomial inversion method should work out well for most situations in which the nonlinearities are independent of frequency. The result of the illustrative example shows that this method is capable to identify the nonlinear systems with acceptable accuracy.
References

Kerschen G., et al. Review: past, present and future of nonlinear system identification in structural dynamics. Mechanical Systems and Signal Processing, Vol. 20, Issue 3, 2006, p. 505593.

Noël J. P., Kerschen G. Nonlinear system identification in structural dynamics: 10 more years of progress. Mechanical Systems and Signal Processing, Vol. 83, 2017, p. 235.

Tanrikulu O., et al. Forced harmonic response analysis of nonlinear structures using describing functions. AIAA Journal, Vol. 31, Issue 7, 1993, p. 13131320.

Özer M. B., Özgüven H. N., Royston T. J. Identification of structural nonlinearities using describing functions and ShermanMorrison method. Mechanical Systems and Signal Processing, Vol. 23, Issue 1, 2009, p. 3044.

Aykan M., Özgüven H. N. Parametric identification of nonlinearity in structural systems using describing function inversion. Mechanical Systems and Signal Processing, Vol. 40, Issue 1, 2013, p. 356376.

Jalali H., Bonab B. T., Ahmadian H. Identification of weakly nonlinear systems using describing function inversion. Experimental Mechanics, Vol. 51, Issue 5, 2011, p. 739747.

Zhongge Zhao, Chuanri Li, Kjell Ahlin, Huan Du Nonlinear system identification with the use of describing functions – a case study. Vibroengineering Procedia, Vol. 8, 2016, p. 33‑38.

Gelb A., Vander Velde W. E. MultipleInput Describing Functions and Nonlinear System Design. McGrawHill, 1968.

Ahlin K. Toolbox for simulation and parameter identification of nonlinear mechanical systems. Proceedings of IMAC 25, 2007.