Eigenvalue Solvers for Structural Dynamics

[Ronankeating]

hi community,

In general I would like to know the what's the general way of solving the moderate or large-scale eigenvalue or algorithms in structural dynmics.

The simple motion equation is as follow.
M*d²X(t)/dt²+C*dX(t)/dt+K*X(t) = F(t).
The bolded expressions are known before the equation solution , mass, damping, stiffness and loading matrix, respectively, so those are known values either time-dependent or constant.

Solution of above equation, including various combinations, e.g. undamped-damped system, constant force, time dependent force etc.., has been given in books, in most of them eigenvalue has been introduced or transformed to A*λ=X*λ form. Among the solution methods it's been mentioned that, power method, eigenvalue decomposition, householder reflectionm, sturm secuence, Wielandt deflation etc.. are most effective, they usually introduce an A matrix and trying to find the eigenvalues and eigenvectors of it.

My question is: Since I have M,C,K,F(t) four matrices how to convert them into A-matrix form, depending on the problem type, in order to commence the solution for the fore-mentioned methods?

Regards,

 

[AlephZero]

A full answer to this would probably fill more than one book, but here's a quick overview.

The right hand side forces F are not relevant for eigensolutions.

Often C is defined in a mathematically convenient way rather than as a physical model of something in the system. Two ways are to define it in terms of the undamped modes of the structure, or as $C = \alpha K + \beta M$ for some constants $\alpha$ and $\beta$. Both those approaches mean you only have to solve an eignenproblem of the form $(M - \omega^2K)x = 0$.

Up to about 1980, the most popular way to solve that was first to reduce the problem to a manageable size (say 500 degree of freedom or less) using methods like Guyan reduction. Since the reduced M was always positive definite, you could do a Cholesky factorization $M = LL^T$ and then get a symmetric A matrix like this:

$(LL^T - \omega^2K)x = 0$
$(L - \omega^2 KL^{-T})L^Tx = 0$
$(I - \omega^2 L^{-1}KL^{-T})L^Tx = 0$
i.e. solve the eigenproblem $(I - \omega^2 L^{-1}KL^{-T})y = 0$
and then find $x = L^{-T}y$

An apparently simpler version of the above, just solving
$(I - \omega^2 M^{-1}K)x = 0$
is not such a good idea, because $M^{-1}K$ is not symmetric and numerical rounding errors may give negative or even complex values of $\omega^2$.

Given more powerful computers after about 1980, an alternative was to find a subset of the modes (with the lowest frequencies) of the "generalized eigenproblem" $(M - \omega^2K)x = 0$, without reducing the model size. A straightforward method that was popular for a while was the "Wilson subspace iteration". This is quite easy to program, and usually works fine for finding say 50 or 100 modes of a model with up to 100,000 DOF.

Most modern commercial FE systems now use some verson of Lanczos iteration rather than Wilson subspace. This is significantly faster for larger models but isn't something you want to code up yourself. The basic math looks simple enough, but it only works practice with careful attention to the effects of numerical rounding.

There are some types of problem where the C matrix is a model of the physical behaviour of the system, not just a mathematically convenient approximation. In that case, you have to calculate the complex eigenvalues (representing the frequency and the corresponding damping) and the mode shapes are also complex (different parts of the structure may have different phase responses in the same vibration mode). This can be solved by converting it into a double-sized matrix, treating the displacements $x$ and velocities $\omega x$ as independent variables:
$$\left(\omega\begin{bmatrix} M & 0 \\ 0 & I \end{bmatrix}
+ \begin{bmatrix} C & K \\ -I & 0\end{bmatrix}\right)
\begin{bmatrix}\omega x \\ x\end{bmatrix} = 0$$
(Here $\omega$ is a complex number, the imaginary part reperesents the frequency and the real part the amount of damping).

In this formulation the matrix containg C and K does not have any "nice" properties like symmetry or positive definiteness, but both "direct" solution methods for small problems, and iterative methods for large ones exist.

 

[Ronankeating]

Thank you in advance,

As you mentioned in your post that, RHS of motion equation is not taken into consideration during solution of eigen-problem process. That means that loading type is completely disregarded, is that valid even if I use the Ritz, Lanczos vectors for solution?

Furthermore IMHO that eigen-problem is just mathematical solution of dynamic problems in engineering, so due to that it introduces some round-off errors, e.g. modal participation ratio, mass participation ratio etc...For example, solving the problem with subspace iteration and with simlutaneous iteration produces different eigen values, solution is volatile in itself, so how can that be reliable solution process, I think it's just serves as a preliminary phase for more complicated calculations e.g. for non-linear analysis. Could you shed some light on this, as well?

Instead of using eigen, why not using direct integration like Wilson theta, central difference or linear accelaration methods which I believe that they are the exact solution for time dependent loading. More interestingly, codes/regulations forces the engineers to use modal superposition which relies on eigen -values , on the other hand if engineer prefers to use the direct integration methods more limits and restrictions are introduced in the codes. Isn't that weird ?

P.S.: For now I can stick to that Rayleigh damping which is expressed as a proportional to stiff and mass matrix for mathematical convenience, doesn't do any harm in my idea.

Regards,

 

[AlephZero]

The basic fact of life is that dynamics is much more more computationally intensive than statics, especially in the time domain. Therefore if you are not prepared to tolerate very long run times, you need a "smaller" model for dynamics than for statics.

Modal reduction is one way to reduce the model size, but as you said it can introduce serious errors if the structure has "stiff" load paths that are only associated with high frequency modes that were excluded from the model. It's a very good idea to check that all the applied load is "going somewhere" in the structural response, and not just being ignored.

On the other hand, the displacement velocities and accelerations usually are dominated by the low frequency modes, which can lead to the situation where the response "looks right" but the corresponding stresses from a modal reduced model are seriously in error.

Static reduction methods like Guyan retain the load paths in the full model, but need too much human expertise to select a good set of variables to capture the full model's mode shapes and frequencies accurately.

One solution to this is to use a combination of modal and static reduction, for example the Craig-Bampton method (a.k.a. component mode synthesis).

Direct integration of the full model is not "exact" and not necessarily even accurate, if the time step size is not correctly chosen. Leaving the choice to the user can be good recipe for getting wrong answers quicker by using bogger time steps!

I don't understand your comment

...For example, solving the problem with subspace iteration and with simlutaneous iteration produces different eigen values, solution is volatile in itself, so how can that be reliable solution process,

The stiffness and mass matrices of a given structural model has only one set of modes. If different eigensolution methods are giving significantly different results, something is wrong with the software IMO.

I'm not a civil engineer, so I can't comment on the codes of practice or regulations for dynamic analysis in CE.

 

[Ronankeating]

Most of your posts, really meets whith what books say about that topic,

With little effort I think that all of the reduction concepts can be learned, for now Guyan works well improving the reduction technique is optimization phase for me.

There is one thing I didn't get or I didn't understand clearly, usually it's been assumed that at each 0.005 sec time-step seismic excitation can be assumed as linearly varying function. Even if I stick to that and doesn't allow the user to meddle up with it does it still stand as a "not exact" method?

Because I'm somehow tended to believe that if I use direct integration with literally enough small time step, that will produce exact results.

I don't understand your comment
...For example, solving the problem with subspace iteration and with simlutaneous iteration produces different eigen values, solution is volatile in itself, so how can that be reliable solution process,
The stiffness and mass matrices of a given structural model has only one set of modes. If different eigensolution methods are giving significantly different results, something is wrong with the software IMO.

Judging from your reponse you got it clearly. This proves that my algorithm is malfunctioning.

 

[AlephZero]

There is one thing I didn't get or I didn't understand clearly, usually it's been assumed that at each 0.005 sec time-step seismic excitation can be assumed as linearly varying function. Even if I stick to that and doesn't allow the user to meddle up with it does it still stand as a "not exact" method?

The time marching output should converge to the exact solution to the equation of motion as the timesteps are reduced.

A typical time marching algorithm samples the value of the applied force at one time point in each step. It can't be more "accurate" than assuming the forcing function is a piecewise linear function between those sampling points, because it doesn't have any more data to work with.

Time marching algorithms introduce errors in the frequencies of the response in each mode, and in the damping levels, compared with the equations of motion. These errors are small if the timestep is small compared with the period of the mode, but for your example of a 5 ms step, the output can not represent any vibration with a higher frequency than two step lengths (this is the same idea as the Nyquist frequency in digital signal processing) which corresponds to a frequency of 100 Hz.

The reduced model probably has many modes with frequences higher than that, which are of no particular interest in themselves, but they represent the "stiff" load paths in the structure. (Hence throwing them all away with modal reduction can be a bad idea!). The time marching algorithm effectively reduces the frequencies of all those modes to something of the order of 100 Hz (the precise details are obviously algorithm dependent).

Another way to look at this is to consider the speed that forces can propagate through the structure. That is simple to answer for an explicit algorithm: a force can only propagate from one point its nearest neighbours in one time step. This is part of the reason why explicit algorithms are conditionally stable - if dynamic effects can't propagate through the structure "one step at a time" at something approximating to the speed of sound in the material, the results won't make any sense.

Implicit unconditionally stable methods have the opposite defect, that dynamics effects can propagate everywhere in the model in one timestep, even though that is physically impossible.

If you want to look at this experimentally, set up a single degree of freedom model of a mass on a spring with no damping. With your step of 5ms, make models with frequencies of say 1, 3, 10, 30, 100, 300 Hz, give them some initial conditions to create free vibration (no applied right hand size force) and compare the frequency and damping of the output with the correct values. Then repeat including some damping in the model (say 10% of critical damping) and see what damping levels you get in the output. The results for 1Hz and 3Hz should be pretty good, but the others will probably not be! The amount of damping (if any) in the hgih frequency modes will probably be an artefact of the time marching algorithm itself, and almost indepenent of the damping matrix in the equations of motion.