explicit methods for dynamical analysis

[Alexzhou]

Hi,

Recently I developed a user element and some material models with OpenSee. And the codes work well in static analysis, dynamical analysis with Newmark method. However, when I try to use explicit methods for dynamical analysis, such as CentralDifference, HHT, ect. , my analysis quickly blows out no matter how small the time step is.

I can hardly find any examples on how to use these explicit methods in the Manual and example scripts. Can anyone point out some sample codes just like quickmain or example tcl for this topic? Any hint will be appreciated.

P.S. what's the difference between Newmark1 and Newmark, HHT1 and HHT, and so on? They doesn't appear in the manul, and the code itself doesn't say much.

Thanks.

Alex Zhou

[Alexzhou]

Any comments?

[fmk]

Sorry about the delay .. i have been on holiday.

CentralDifference is explicit, HHT implicit.

The usage in a script for the explicit method is no different than the implicit; the only difference being that the time step must be carefully chosen to avoid stability problems and depending
on your model, this may require a very very small dt. with no damping: dt < Tn/PI (Tn smallest period), with damping present it is even smaller.


Newmark and HHT should be used instead of Newmark1 and HHT1 for nonlinear problems. The
difference is in the initial prediction of the accel at t+dt, in Newmark, HHT it is the accel at last time step; in Newmark1 and HHT1 it is 0.0 (for linear problms the 0.0 results in 1 less step to get correct disp at t+dt)

[Alexzhou]

Hi, Frank

You reponse make things more clear now. Thanks.

I did do some tests on the central difference(CD) methods, however no luck.
I checked the CD code, and still not very clear how the RHS of the SOE was formed since it is different than what I expected. I think I will dig deeper into the code later.

Do you have any woked example C++/script file for CD test, which may help me to figure out things I didn't got.

Alex

[fmk]

this is a modified version of example6.1 .. script does both Newmark & CentralDifference(NoDamping) .. will look into CentralDifference itself some more .. nore that for CentralDifference the algorithm is Linear

proc run {integratorType} {
wipe
model basic -ndm 2 -ndf 2

nDMaterial ElasticIsotropic 1 1000 0.25 6.75
set Quad quad
set eleArgs "1 PlaneStrain2D 1"

set nx 10; # NOTE: nx MUST BE EVEN FOR THIS EXAMPLE
set ny 4
set bn [expr $nx + 1 ]
set l1 [expr $nx/2 + 1 ]
set l2 [expr $l1 + $ny*($nx+1) ]

# now create the nodes and elements using the block2D command
block2D $nx $ny 1 1 $Quad $eleArgs {
1 0 0
2 40 0
3 40 10
4 0 10
}
# Single point constraints
# node u1 u2
fix 1 1 1
fix $bn 0 1

# Gravity loads
pattern Plain 1 Linear {
load $l1 0.0 -1.0
load $l2 0.0 -1.0
}

integrator LoadControl 1.0 1 1.0 10.0
test EnergyIncr 1.0e-12 10 0
algorithm Newton
numberer RCM
constraints Plain
analysis Static
# Perform the analysis
analyze 10

recorder Node -file Node$integratorType.out -time -dT 0.5 -node $l1 -dof 2 disp
recorder plot Node$integratorType.out CenterNodeDisp$integratorType 625 10 625 450 -columns 1 2 -dT 0.5


# ---------------------------------------
# Create and Perform the dynamic analysis
# ---------------------------------------

# Remove the static analysis & reset the time to 0.0
wipeAnalysis
setTime 0.0

# Now remove the loads and let the beam vibrate
remove loadPattern 1

# Create the transient analysis

if {$integratorType == "Newmark"} {
test EnergyIncr 1.0e-12 10 0
algorithm Newton
numberer RCM
constraints Plain
integrator Newmark 0.5 0.25
system ProfileSPD
integrator Newmark 0.5 0.25
analysis Transient
analyze 500 0.5
} else {
algorithm Linear
numberer RCM
constraints Plain
integrator CentralDifferenceNoDamping
system Diagonal
analysis Transient
analyze 10000 0.025
}
}

run Newmark
run CentralDifference

[Alexzhou]

I got the CenterDifferenceNoDamping method right immediately after I use the diagonal SOE. And better when I change my element routine to give lumped diagonal mass matrix.

After I did some debugging, I found the integration orders I used is just not high enough for correct consistent mass matrix. Now it works ok with the CenterDiffence method.

The script file you gave is of great help.

And I found the CenterDifference code are O.K. except at one place:
CentralDifference::update(const Vector &X)
{
....
...

*Utm1 = *U;
*U = X;

*Udot = X;
*Udot -= *Utm1;
*Udot *= c2;

Udotdot->addVector(0.0, *Utm1, 1.0);
Udotdot->addVector(1.0, X, 1.0); //X=*U here;
Udotdot->addVector(1.0, *U, -2.0);
*Udotdot *= c3;
...
}
It doesn't seems like the way CenterDifference calculate the acceleration, which is as in the head file:
// An = 1/dT^2 (Un+1 - 2Un + Un-1)
Anyway, this seems doesn't affect the solution since the velocity and acceleration vector will be overwritten in the newStep() routine by the "garbage quantities needed for vel & accel at nodes", which I finally found that they were calculated in purpose to account for the two parts of RHS of SOE: [(1/dT*dT)M](2Un-Un-1) + [(2/dT)C] Un-1; And I ignored these at the the first sight, my mistakes.

[seocy]

Dear user,

I tried to compare implicit and explicit methods for a transient analysis. For some reasons, I have stability problems with explicit method. No matter how small time step I chosed, it always blew up in a very few seconds or fails even at the first time step.

For example, I tried the example 4 (simple 2D RC portal frame with transient analysis) and modified a beam column element to dispBeamColumn elements, used linear algorithm either with centralDifferenceNoDamping, centralDifference, or NewmarkExplicit method.
Here I post the modified example script. Would you tell me what I did wrong?

Code: Select all

# This is a modified LibAnalysisDynamicParameters.tcl
constraints Transformation; 
numberer RCM 
system BandGeneral 
# or system Diagonal 
algorithm  Linear        
integrator NewmarkExplicit 0.5 
# or integrator CentralDifference
# or integrator CentralDifferenceNoDamping 
analysis Transient 

one of the error message I encountered is
" ... algorithm failed at time 0.001"
after some warning message.
I used a time step 0.001. So, this did not even run at the first step. I wipeAnalysis after gravity analysis.

I am using OS1.7.5 on Windows.
I appreciate your help in advance.

[seocy]

Anybody tried explicit method for nonlinear transient analysis?
For implicit method, it is working fine but by switching to explicit method, the SOE seems to have zero diagonal terms. Why does this thing happen?

[mdcarpio]

Hi fmk,

I followed the input file you posted with explicit integrators. I followed it to construct my model with Central Difference but I still cannot get opensees to run it. The problem is the linear solver. Below is my input file, could please take a look at it?
Any help would be appreciated, Thanks
#####################################

# 3 m 4
# [+]-------------------------[+] -> dof1 (disp)
# o rigid beam o
# ¦ ¦ o : pin
# ¦ ¦ [+] : node
# ¦ ¦ @ : rotational spring
# 2 ¦ 5 ¦
# [+] dof2 (rot) [+] dof3 (rot)
# 1 @ 6 @
# ---[+]--- ---[+]---
# /////// ///////

wipe all;

# variable definition
set dt 0.01; # time step for analysis
set g 386.4; # acceleration of gravity
set m [expr 20/$g]; # mass, 20 kips /g

# model Command
model BasicBuilder -ndm 2 -ndf 3; # define the model builder, ndm = #dimension, ndf = #dofs
file mkdir Results; # create output folder

node 1 0 0;
node 2 0 0;
node 3 0 42 -mass [expr $m/2] 0 0;
node 4 42 42 -mass [expr $m/2] 0 0;
node 5 42 0;
node 6 42 0;

fix 1 1 1 1;
fix 6 1 1 1;

equalDOF 3 4 1; # diaphragm (roof)
equalDOF 1 2 1 2; # constrain displacement of rotational spring (left)
equalDOF 6 5 1 2; # constrain displacement of rotational spring (right)

geomTransf Linear 1; # Linear transformation

uniaxialMaterial Steel01 1 6000 26087 0.01; # bilinear model
uniaxialMaterial Steel01 2 60 26087 0.01; # bilinear model

element elasticBeamColumn 1 2 3 1000 29000 11.3 1; # left column
element elasticBeamColumn 2 5 4 1000 29000 11.3 1; # right column
element zeroLength 3 1 2 -mat 1 -dir 6; # rotational spring (left)
element zeroLength 4 6 5 -mat 2 -dir 6; # rotational spring (right)

recorder Node -file Results/Displacement.out -node 3 -dof 1 disp; # displacement

set duration 31.18; # duration of ground motion in seconds
set groundmotion "Path -filePath EL_Centro.txt -dt 0.02 -factor $g"

pattern UniformExcitation 1 1 -accel $groundmotion
numberer RCM # renumber dof's to minimize band-width (optimization)
constraints Transformation
system Diagonal
integrator CentralDifferenceNoDamping
algorithm Linear
analysis Transient

# perform the transient analysis
set no_steps [expr round($duration/$dt)+1];
analyze $no_steps $dt

wipe all;

[mdcarpio]

Problem solved, my model had massless DOF.

When working with Central Difference, all DOF must have mass for the scheme to work.

I wonder if you could overcome this by condensing out the massless DOF? Any thoughts?

Maikol

[sheikhakbari]

Hello fmk

I have a 3-span bridge model. As I rotate the ground motion or change the ground motion time history, I have convergence problem in some cases. I decided to try Explicit Newmark Method because I couldn't get the convergence no matter how much I decreased the time step. By knowing the number of DOF of my structure and running modal analysis, I calculated the minimum period of my structure . But when I decrease the time step to the amount which satisfies the stability issue of the method, the running never ends ! I even tried different system commands but no success. Can you please give me any suggestion to get the convergence?

[fmk]

the only convergence issue you can have when using an explicit method resides in the iterations needed by the force based element .. do you have these elements in the model?

[sheikhakbari]

Yes,I do.All of the columns are modeled as nonlinear beam column element.I tried to get the convergence by changing many factors like number of integration points,integration types,maximum number of integration and the tolerance.However,I got not success.

[fmk]

you have to beef up the sections or switch to displacement based elements (they need no convergence) however you need a few elements per column as you probably will not want linear curvature assumption along elements length if just 1 element

[Gholamreza]

So informative, Thanks