Dynamic Stability Analysis

1.2. Nonlinear Stability Analysis

The equation of motion of a geometrically nonlinear structural model is given by

Mx + r(x, x) = f(t). (45)

For neighboring trajectories, the tangential equation of motion may be utilized to describe temporal evolution of the difference y(t)

My + Cy + Ky = 0. (46)

To analyze the dynamic stability behaviour of nonlinear systems an integration of Equation (45) is necessary until stochastic stationarity is reached. In each time step, the tangential stiffness matrix K has to be determined. With this kind of analysis a criterion for sample stability is developed. In order to speed up explicit time integration, this equation can be projected into a subspace of dimension m as spanned by the eigenvectors of the undamped system corresponding to the m smallest natural frequencies (Bucher, 2001). These eigenvectors are the solutions to

^K(xstat) — &>2M^ Ф = 0; i = 1.. .m

In this equation, xstat is chosen to be the displacement solution of Equation (45) under static loading conditions. The mode shapes are assumed to be mass normalized. A transformation x = Ф v and a multiplication of Equation (45) with Фт represents a projection of the differential equation of motion for the reference solution into the subspace of dimension m as spanned by the eigenvectors:

v + Ф Tr(x, X) = Фтf. (48)

The integration of Equation (48) by the central difference method (Bathe, 1996) requires a minimal time step.

The time integration in the subspace and the computing of the restoring forces on the full system causes the following problem: If the start displacement or velocity vector of the time integration is not zero, for example by static loading, the projection of this vectors into the subspace is an optimiziation problem caused by the higher number of variables in the full space. By using a least square approach:

v = Ф-1х; Ф-1 = (Фтф Фт (49)

this projection is optimal approximated, but not suitable for a subspace spanned by a small number of eigenvectors. A possibility to handle this, is to start the time integration in the subspace with a displacement and velocity vector equal to zero. The start vectors have to be saved in the full system and the restoring force vector has to be computed by addition of the start and the time integration vectors:

Г(x, X) = r(Xstart + Фv, Xstart + ФV);

v(t = 0) = v (t = 0) = 0. (50)

In the investigated cases the start vector xstart is the static displacement vector, the start velocities are assumed to be vanished.

To analyze the stability behaviour of the reference solution x0(t), the long-term behavior of the neighboring motion (Equation (46)) is investigated. To reduce the dimension of the equation system, this equation can be projected into the same or a smaller subspace as Equation (48). Transformed into the state space description we obtain:

From this equation, the Lyapunov exponent X can be determined by a limiting process:


A.(x0, s)= lim – log ||0(xo, t)s\ (52)

t УЖ t

in which s is an arbitrary unit vector. In Equation (52), 0(x0, t) is the transition matrix from time 0 to t associated with Equation (51). Based on the multiplicative ergodic theorem (e. g. Arnold and Imkeller, 1994) the Lyapunov exponent can also be calculated as an expected value:

In the current investigation, the norm ||0(x0, t)s|| is expressed in terms of

0(x0, t)s|| < ||0(x0, t)|| ■ ||s|| = ||0(x0, t)||.

Finally, this result is used in calculating the Lyapunov exponent according to Equation (52) by using a matrix norm equal to the eigenvalue pmax of 0(x0, t) with the maximum absolute value. The time domain t has to be taken large enough that the Lyapunov exponent convergences to a stationary value. For the statistical estimation of the convergence of the Lyapunov exponent, Equation (53) is suitable.