Non-stationary random response of MDOF Duffing systems

This paper presents a method to solve non-stationary random responses of nonlinear multi-degrees-of-freedom (MDOF) Duffing systems subjected to evolutionary random excitations. Specific phase-lags between the excitations can also be taken into account. The power spectral density (PSD) of the input excitations is not confined to simple white noise or filtered white noise, in fact it can also take more complicated forms. The MDOF nonlinear random differential equations are iteratively solved by means of the Equivalent Linearization Method (ELM) combined with the Pseudo Excitation Method (PEM). This combined method is easy and efficient. Two examples are given in which this method is well justified by the Monte-Carlo numerical simulations. Although only a Duffing model is dealt with in this paper for computational simplicity, the proposed method is in fact quite general, e.g. it can also deal with nonlinear hysteretic structures that will be dealt with in a separate paper.


Introduction
The response of many real structural systems subjected to external random excitations can be described by nonlinear random differential equations.In general, approximation methods have to be used to solve such nonlinear equations.Only in very few cases can an exact solution for stationary response be found in closed form.Although the non-stationary random response analysis of MDOF nonlinear systems is much more difficult than such problems, some pioneering work has laid a good foundation in this field [1][2][3][4][5][6][7][8][9][10][11].
One of the most popular methods used in the study and engineering applications of nonlinear random vibration problems is equivalent linearization, particularly if the random excitations are stationary.For nonlinear systems subjected to non-stationary random excitations, however, the equivalent linear systems are time-dependant, its solution is still not so easy.In particular, if the excitations are not uniform, i.e. phase lags exist between the excitations, the solution will become more difficult.
It is well known that for an arbitrary nonlinear system subjected to stationary random excitations, an equivalent time-independent linear system can be found readily based on the least-quadratic criterion.If the nonlinear system is instead subjected to non-stationary random excitations, then at any time t, an equivalent linear system can be found similarly based on the least quadratic criterion.And based on this idea, a nonlinear SDOF (Single Degree of Freedom) system subjected to a non-stationary random excitation is solved iteratively [11], wherein the analytical solution of the time-dependent mean-square response for the linear time-independent system must be available.Therefore its application field is limited.
In this paper, the pseudo excitation method (PEM) [12][13][14][15][16] is extended to solve the non-stationary random responses of MDOF Duffing systems subjected to evolutionary random excitations.It is not necessary to know the analytical expressions of the time-dependent mean-square responses for the linear MDOF time-independent system, therefore the proposed method is not confined to structures subjected to white-noise or filtered white-noise excitations.In addition, phase-lags (or time-lags) are permitted between the excitations.
Two examples are given in which the Monte-Carlo simulation is used to justify the proposed method.

Linear systems subjected to evolutionary multi-phase random excitations
Consider a linear structure subjected to an evolutionary random excitation in which g(t) is a given slowly varying modulation function, while x(t) is a zero-mean-valued stationary random process of which the auto-PSD S xx (ω) is also known.The PSD curve of f (t) has constant shape, however its amplitude varies with g(t).In other words, a stationary random process x(t) is modulated by a deterministic function g(t).In order to compute the PSD matrices of two arbitrary response vectors {y(t)} and {z(t)}, the following pseudo excitation should be used [13,15] Under the action of this deterministic loading, the equations of motion Eq. (1) of the linear structure are in which [M ], [C] and [K] are time independent mass, damping and stiffness matrices, {p} is a given constant vector.For seismic problems, the structure is initially at rest, namely, {y 0 } = { ẏ0 } = {0} when t = 0. Thus the time-histories of two arbitrary responses of interest, denoted as y(ω, t) and z(ω, t), can be computed in terms of the Newmark, Wilson-θ schemes or the precise integration method [15] with ω regarded as a parameter.It has been verified that [13] the auto-and cross-PSD functions of these responses can be accurately computed by using the following equations If the interested responses are two vectors, denoted as {y(t)} and {z(t)}, the corresponding time history {y(ω, t)} and {z(ω, t)} should be solved from Eq. (3) before the following equations are used to give their accurate auto-and cross-PSD matrices [13] [S yy (ω)] = {y(ω, t)} * {y(ω, t)} T (7) It has also been verified that if the precise integration method [15,17] is used in solving Eq. ( 4), the efficiency would be much higher than using Newmark method.
If the random excitation f (t) in Eq. ( 2) is a vector, specific time-lags exist between all its components, then in which a j (j = 1, 2, . . ., n) are given non-negative real numbers, t 1 , t 2 , . . ., t n represent the time lags between the forces applying on different points of the structure.a k = 0 means no load applies on the k-th DOF.F (t) is a stationary random process, its auto-PSD S F F (ω) is known.Using Wiener-Khintchine formula and denoting herein [R 0 ] and {q 0 } are matrix and vector with all elements be unity.Assume that {y(t)} is an arbitrary response vector due to the excitation {f (t)}, then If {y k (t k )} and {y l t l )} are two kinds of response vectors, then their cross-correlation matrix would be Substituting Eqs (10), (11) and (12) into Eq.( 14) gives in which It is known from Eq. ( 15) that when t k = t l = t, the cross-PSD matrix of {y k } and {y l } is It is also known from Eq. ( 16), that I k (ω, t) is the transient response due to the deterministic excitation [G(t)][V ]{q 0 }e iωt .Therefore by using the following pseudo excitation . . .
the resulted response must be Therefore, according to Eq. ( 17), one obtains Let k = l, Eq. ( 20), gives the auto-PSD matrix of {y k }.
By neglecting the subscript k, Eq. ( 7) is obtained.For two responses {y(t)} and {z(t)}, just follow the above process, Eq. ( 8) can also be derived.

Duffing systems subjected to evolutionary multi-phase random excitations
The nonlinear equations of motion can be written as in which {g} represents the nonlinear property of the structure.For the MDOF Duffing system as shown in Fig. 1, {g} has the following form {x} is the story drift vector, {µ} is the vector which characterizes the intensity of system non-linearity, {f (t)} is a zero-mean-valued evolutionary random process.
Clearly, in the above steps, the analytical expressions of the time-dependent mean-square response for the equivalent linear MDOF time-independent system need not be known, therefore this method applies not only to structures subjected to simple white-noise or filtered white-noise excitations, but also to structures subjected to more complicated random excitations.It can also be seen from Eq. ( 9) or Eq. ( 11) that the proposed method permits phase-lags (or time-lags) existing between the excitations.

Sample time history curves
A uniformly modulated evolutionary random process with the form of f (t) in Eq. ( 2) can be simulated by the following expression [16,18,19] )∆ω, k = 1, 2, . . ., N ω l and ω u are the lower and upper bounds of the frequency region, ϕ k is a random real number uniformly distributed in the region [0, 2π].

Numerical simulation for nonlinear systems
Provided the equations of motion of a nonlinear dynamic system is in which G( Ẏ , Y ) is the nonlinear term.The simulation process for the nonlinear system is as follows: (1) Substitute an excitation according to Eq. ( 29) into the right hand side of Eq. (32), compute the displacement vectors {Y } at a series of time points by means of a time history method; For the displacement values corresponding to each element in the vector {Y } at a series of discrete time points, the Fast Fourier Transformation (FFT) scheme is used to obtain the PSD values at a series of frequencies which (2) will then be used to compute the PSD vector {S Y Y (ω)} of the displacement vector at these frequencies.
(3) Integrate the PSD vector {S Y Y (ω)} in the effective frequency region to obtain the corresponding response variances {σ 2 Y }. (4) The above steps (1)∼(3) are repeated for a large amount of excitation samples to generate corresponding {σ 2 Y } samples, and the statistics.For nonlinear equations, iterations must be executed at each time step.The statistics of the responses at any arbitrary time will then be computed by means of the assemble average of a great deal of response samples.
It is noted that the time step size must be small enough to ensure the integration precision [18,19].
in which T is the time-lag between f 1 (t) and f 2 (t).The PSD functions of f 1 (t) and f 2 (t) are equal to each other, and both will simultaneously take one of the following forms: The modulation function has the form When the excitation PSD is a white noise expressed by Eq. ( 34) and for T = 0.0 or T = 1.0, without or with timelags between the excitations, the displacement standard deviations computed by PEM are shown in Figs 2(a)-(b).Such standard deviations are also computed using Monte-Carlo simulation method, with 4000 samples, ω l = 0.0, ω u = 100.0,N = 2000; and the results are also shown in Figs 2(a)-(b).Clearly, these two methods give quite close results.When T = 0.0, the maximum difference is less than 3.0%, and when T = 1.0 the errors are mostly within 5.0%, the maximum is 7.0%.As the PEM is an accurate method for linear systems, the errors come from the linearization of the nonlinear system.
When using Eq.(34) as the PSD of both excitations, Fig. 3 gives the comparison of the two methods for T = 0.4, 1.0 and 1.6.It is shown that for all cases, the two methods give quite close results; it is also shown that the non-stationary responses of the nonlinear system rely heavily on the time lag T .The PSD form of Eq. ( 35) is somewhat similar to a wind-gust spectrum.Under this excitation and with T = 1.6, the responses due to the two methods still meet quite well, as shown in Fig. 4.
The standard deviation curves of the non-stationary displacement responses are computed using PEM and Monte-Carlo method with 1000 samples, ω l = 0.0, ω u = 50.0,N = 2000, respectively.And the results are both shown in Fig. 5. Once again, It is found that the results from PEM agree well with those from Monte-Carlo simulation.The differences are all within 6%.The first eigenvalue of the structure is 19.11 (for µ i = 0) and 23.79 (for µ i = 0.15) respectively, the difference is 11%, which also signifies a rather strong nonlinear characteristics.Figure 5 also shows that the maximum difference, about 6%, takes place in the strong excitation region.This phenomenon shows that stronger excitations will lead to stronger non-linearity, and so more severe errors will take place due to the equivalent linearization.It should be noted that in essence the proposed method does not apply to problems with strong non-linearity.

Conclusions
The pseudo excitation method [12][13][14][15] is extended to the analysis of non-stationary random responses of nonlinear structures.Duffing models are used in the given examples for computational simplicity.In fact, the proposed method is quite general, for instance it can also deal with structures with nonlinear hysteretic damping [16].This method not only has good precision, but also has high computational efficiency, especially for structures with many degrees of freedom.Its implementation is also rather convenient.

Fig. 4 .
Fig. 4. Displacement standard deviation for non white spectrum and time lag T = 1.6 s.