[gdbb89]

Hi All,

I'm working on a simple structural dynamics FEA program in MATLAB. For now I am using a thin rectangular plate as the object to be analyzed, and I'm meshing it with rectangular plate elements with three degrees of freedom per node (one translation, two rotation). I have successfully implemented a normal modes solution that works seamlessly, however I am now working on a transient base excitation solver and hitting a hard wall with the numerical integration. My solutions are blowing up except for the case of very small timesteps (~1e-7 sec), even when I am using only ~100 elements (~350 DOFs).

Questions:

1) How do I calculate the stable time increment for an MDOF system using a time-marching scheme, e.g. when using the Newmark-beta method to simulate a dynamic transient response?

2) If there is a better integration method for transient response in structural dynamics, in terms of greater stability and quicker computation (likely from allowing larger stable timesteps), what would you suggest?

Some background on my project: The test/control plate for my program is 6" x 4" x 0.125". I'm integrating the response in modal coordinates and only using contributions from the 10 lowest modes, with the highest having a frequency of 23,142 rad/s (3,863 Hz). The excitation is a 1 lb, 5 ms half-sine shock pulse on an interior node ((x, y) = (1.75 in, 1.75 in)) with the four corner nodes pinned. I am currently using the Newmark-beta method for integrating the response with beta = 1/6 and gamma = 1/2. These are all mostly arbitrary, but this is my starting point.

I've run a bit of a numerical stability case study using this plate, and I can easily see that the stable time increment is a function of (modal) damping. It follows the behavior of the function f(x) = x*e^(-x) or something very similar. At any rate, I cannot seem to find an answer on Google on stable time increments and don't have a formal academic background in numerical methods.

Any and all help is greatly appreciated!

[gdbb89]

I sorted this out. The main issue was my implementation of the scheme was completely wrong, which I have not had the time to figure out yet. I glossed over the fact that I was reading in many places that the Newmark-beta method is unconditionally stable for beta >= 1/4. Even when I used that value, I would still get the unstable behavior at sufficiently "large" timesteps. There is a completely different way of updating than what appears at face value from, for example, the Wikipedia page.

I was updating by taking {udd^(n+1)} = [M]^(-1) * [ [F] - [C]*{ud^(n)} - [K]*{u^(n)} ] as well, which is wrong. If anyone with the same problem finds this thread, pg. 777 of K.J. Bathe's "Finite Element Procedures" has the proper implementation laid out.