Dynamic analysis diverges

[ja_abell]

Dear OpenSEES community,

I am performing a linear transient analysis on a big 3D earth-dam model.
This is a 3D model composed only of bbarBrick elements with ElasticIsotropic3D materials (different material for each brick).
I am using an acceleration record of 'El Centro' by assigning a pattern UniformExcitation load in the transverse direction.

I am assigning damping with the rayleigh command. The alpha and beta parameters are obtained by:

alpha = xi*w1
beta = xi/w1

Where xi is the percent critical damping at mode 1, and w1 is the angular frequency of mode 1. (I've also tried setElementRayleighDampingFactors, which is what I want to use in the first place to model strain dependent damping in the future).

My analysis options are:

system SparseSPD
constraints Transformation
test NormDispIncr 1.0e-3 40 1
algorithm Linear
numberer RCM
analysis Transient

I've tried many different combinations, but my analysis always diverges after a few seconds of analysis.

The analysis is performed inside an external for loop as in

for {set step 0} {$step < $NSteps} {incr step} {

...code

#Perform one analysis step
analyze 1 $DtAnalysis

....more code
}

This is so I can retrieve the element strains at each time step (i dont want big output files, I just need the element maximum strains), and perform some checking and estimate analysis time (usually about ~20hrs on my PC for a full analysis).

[fmk]

what integration scheme are you using?

[ja_abell]

I've tried Newmark 0.5 0.25, HHT with values 1 0.9 0.8 0.7 (it is 1+alpha right?), and other integrators such as AlphaOS integrator.
All diverge quickly.

Thank you for your reply.
José Abell

[fmk]

try using Newton for the algorithm and put a 1 at end of convergence test to have a look at the norms of disp incr and unbalance at each step (the Linear algorithm ignores the convergence test)

use the rayleigh command for damping. it does work with the integrator but not for long.

alpha=1 in HHT corresponds to Newmark .. you want alpha between 0.5 and 1.0 to introduce some numerical dissipation.

[ja_abell]

Thank you for the help. Indeed i needed to iterate more and the linear algorithm wasn't providing accurate enough solutions.

Is it normal that within a given step I have to do more than 50 iterations to get convergence?

Sometimes even that is not enough.

Should I just use trial and error to specify the maxunbalance and number of iterations?

Thank you once more.
José Abell

[ja_abell]

An update. I took a longer analysis. The integration still diverges though slower. I reach unbalance tolerances of 1e-16, but with each iteration it becomes harder to converge until the algorithm fails altogether.

Something is introducing instability and I don't know what.

What else can I test?

[fmk]

it is not normal for an elastic problem to reuire more than 40 iterations .. 2 for a linear problem should be about the max .. have you tied the regular brick instead of bbar?

[ja_abell]

By changing to stdBrick element, i got a better convergence but still there were some problems at the end. After a while the program exceeded the maximum of 100 iterations.

I isolated the problem to a few elements where the instability starts and then spreads to the rest of the model until theres no convergence at all. There seems to be nothing wrong there, but it is in the zone where the maximum displacements/strains occur.

I'm going to try with HHT and some added numerical damping. Should another solver work better? I have many (~4000) brick elements and active nodes (~6000 about 18000 Dofs).

Does setElementRayleighDampingFactors work?

Thanks for all the advice.
Jose A.

[fmk]

you can use the region command to increase the amount of damping in that area where the problem starts.

[ja_abell]

I solved the problem. Mi final settings are

system UmfPack
constraints Plain
test NormDispIncr 3.0e-5 100 1
algorithm Newton
numberer RCM
integrator HHT 0.95 (Newmark also works now).
analysis Transient

The change from SparseSPD to UmfPack proved vital.

I changed the bbarBrick elements to stdBrick, and the ElasticIsotropic3D material to ElasticIsotropic. Now it takes two solver iterations for each time step. I'm assigning setElementRayleighDampingFactors.

The solutions 'look' nice, so i think I'll stick to them.

Thanks for all the help, I really appreciate. I'll work up a demo of what I'm doing so that others can benefit from it.

José A.

[fmk]

it's interesting that you probably needed pivoting to get the solution .. SPD problems in theory should not need it .. if you can you should try to see if you can get ProfileSPD to work on the problem (it might not due to mem requirements if large 3d model) .. this is just to make sure the problem is SPD (as SparseSPD ignores negative values on the diagonal) .. if it is not SPD you might have some bricks with incorrectly numbered nodes.

[ja_abell]

I tried ProfileSPD, the model uses about 1GB ram. The program is not complaining, so I guess it's not a numbering problem. I've checked this a couple of times, but I'm not feeling too confident after all.

The time it's taking to solve the problem with ProfileSPD is about 4 or 5 times the time UmfPack takes, so I'll revert to the latter.

While using UmfPack I noticed that some instability occurs for a real acceleration record when using a fine time step, a coarser time step shows no instability. Also, this instability disappears if I fix the X direction and only consider accelerations in the Y direction, using the small time step.

It would be nice to compile a 'recomendations' document for each function in the OpenSEES workframe from the experience in the forums.

[fmk]

it's something people can add to the wiki document for each command.

[ja_abell]

The problem is still occurring, just not as fast as before so I missed it. It occurs further along the integration.

I'm out of ideas now.

[ja_abell]

There is one freak symptom. If I make the integration step shorter, the instability gets worse!

The lower the step size the earlier it occurs.