Explicit Integration Scheme

[suusie]

Dear Frank and users,

I want to use explicit integration method for the nonlinear response history analysis.
For the first step of this, I made a simple 1-bay-1story elastic frame and tried to analyze by explicit scheme (Central Difference).

However, the obtained results was very strange (contained ".#IND" in the output displacement history) compared to the results from the implicit scheme.

I have been wondering what is wrong with the .tcl for some time... The script is attached below. I appreciate if you can point out the mistake in it.

----------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3;

# define geometry
set m [expr 1000./9810]; # mass=kN/(mm/sec2)... Mass for vertical & horizontal directions
set Tn 0.7; # sec (Natural Period)
set pi [expr acos(-1.0)];

# define nodal coordinates:
node 1 0. 0. ;
node 2 5000. 0. ;
node 3 0. 3000. -mass [expr 0.5*$m] [expr 0.5*$m] [expr 0.5*$m];
node 4 5000. 3000. -mass [expr 0.5*$m] [expr 0.5*$m] [expr 0.5*$m];

# Single point constraints -- Boundary Conditions
fix 1 1 1 1;
fix 2 1 1 1;

# Columns and Beam Properties
set K [expr pow(2*$pi/$Tn,2)*$m]; # KN/mm
set E 200.0; # KN/mm2

# define geometric transformation:
set TransfTag 1;
geomTransf Linear $TransfTag ;

# Define Elements----------------------

element elasticBeamColumn 1 1 3 [expr 1e12] $E [expr $K*pow(3000., 3)/(24*$E)] $TransfTag; # Columns
element elasticBeamColumn 2 2 4 [expr 1e12] $E [expr $K*pow(3000., 3)/(24*$E)] $TransfTag; # Columns
element elasticBeamColumn 3 3 4 [expr 1e12] $E [expr 1e12*$K*pow(3000., 3)/(24*$E)] $TransfTag; # Beam

timeSeries Path 1 -dt 0.01 -filePath TAKY.th -factor [expr 0.50*9810.];

# set damping based on first eigen mode
set freq [expr [eigen 1]**0.5]
rayleigh [expr 2*0.02*$freq] 0. 0. 0.

pattern UniformExcitation 1 1 -accel 1;

# create the analysis (explicit scheme)
wipeAnalysis;
constraints Transformation
numberer RCM
system SparseGEN
algorithm Linear
integrator CentralDifference;
analysis Transient

analyze [expr 10*4096] 0.001; # apply 10*4096 steps for 0.001-sec time steps in analysis

puts "Done!"
wipe

[fmk]

could be that the response just blew up because you do not meet the stability requirements of CentralDifference for this model if using a time step of 0.01.

[suusie]

Thanks for your reply,

I just changed the analysis time step to 0.000001 but I still got very strange displacement response (extremely large)...

As you see the above script, the frame is essentially single-degree-of-freedom system, and the first modal period is set 0.7 second...
I think the stability requirements of CentralDifference was met in the above .tcl script if analysis time step is 0.000001..

Actually similar question was posted previously but it was not answered at that time, so I couldn't have any clue to find a good way.
(http://opensees.berkeley.edu/community/ ... od#p104733)

I am using OpenSees Version 2.5.0 (rev 6248) 32-Bit.

I am still wondering what's wrong with my .tcl... Any further suggestion regarding this problem is very appreciated!

[Jhno]

Regardless of any method .. I see you make mistake when defining the integer.
Maybe just changing that would help. I had weird result coming from little error such a these in the past.
set m [expr 1000./9810] -> set m [expr 1000./9810.] ; since m seems to enter in alot of your caclulation.

element elasticBeamColumn 1 1 3 [expr 1e12] $E [expr $K*pow(3000., 3)/(24*$E)] $TransfTag -> same
Easy to check

[Jhno]

So I ran your code and made some minor mod in the way I told you about (integer number). Also there is a error in the script when you define the damping and also with the analyze line.
I also used m and N and placed random value. Other values will work. I get these value at the end.
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 -1.11348e-038 -1.11348e-038
0 0 -5.01468e-024 -5.01468e-024
0 0 -9.13525e-024 -9.13525e-024
0 0 -1.87886e-023 -1.87886e-023
0 0 -2.61329e-023 -2.61329e-023
0 0 -4.06249e-023 -4.06249e-023
0 0 -4.9522e-023 -4.9522e-023
0 0 -7.30007e-023 -7.30007e-023
0 0 -4.30502e-023 -4.30502e-023
0 0 -5.81547e-022 -5.81547e-022
0 0 6.27717e-021 6.27717e-021
0 0 -8.65058e-020 -8.65058e-020
0 0 1.16792e-018 1.16792e-018
0 0 -1.57995e-017 -1.57995e-017
0 0 2.13708e-016 2.13708e-016
0 0 -2.8907e-015 -2.8907e-015
0 0 3.91009e-014 3.91009e-014
0 0 -5.28895e-013 -5.28895e-013
0 0 7.15407e-012 7.15407e-012
0 0 -9.67689e-011 -9.67689e-011
0 0 1.30894e-009 1.30894e-009
0 0 -1.77053e-008 -1.77053e-008
0 0 2.39489e-007 2.39489e-007
0 0 -3.23943e-006 -3.23943e-006
0 0 4.38179e-005 4.38179e-005
0 0 -0.0005927 -0.0005927
0 0 0.00801711 0.00801711
0 0 -0.108443 -0.108443
0 0 1.46684 1.46684
0 0 -19.8412 -19.8412


Which, if you take a look..are very very large values!
model BasicBuilder -ndm 2 -ndf 3;

# define geometry
set m [expr 1000./9.81]; # mass=kN/(mm/sec2)... Mass for vertical & horizontal directions
set Tn 0.7; # sec (Natural Period)
set pi [expr acos(-1.0)];

# define nodal coordinates:
node 1 0. 0. ;
node 2 5.0 0. ;
node 3 0. 3. -mass [expr 0.5*$m] [expr 0.5*$m] [expr 0.5*$m];
node 4 5. 3. -mass [expr 0.5*$m] [expr 0.5*$m] [expr 0.5*$m];

# Single point constraints -- Boundary Conditions
fix 1 1 1 1;
fix 2 1 1 1;

# Columns and Beam Properties
set K [expr pow(2*$pi/$Tn,2)*$m]; # KN/mm
set E 2.e9; # N/m2

# define geometric transformation:
set TransfTag 1;
geomTransf Linear $TransfTag ;

# Define Elements----------------------

element elasticBeamColumn 1 1 3 [expr 1.e12] $E 1.e10 $TransfTag; # Columns
element elasticBeamColumn 2 2 4 [expr 1.e12] $E 1.e10 $TransfTag; # Columns
element elasticBeamColumn 3 3 4 [expr 1.e12] $E 1.e10 $TransfTag; # Beam

timeSeries Path 1 -dt 0.01 -filePath TAKY.th -factor [expr 0.50*9.810*0.001];

# set damping based on first eigen mode
set freq [expr [eigen 1]*0.5]
rayleigh [expr 2*0.02*$freq] 0. 0. 0.

pattern UniformExcitation 1 1 -accel 1;
recorder Node -file Noeud_2.out -nodeRange 1 4 -dof 1 disp
# create the analysis (explicit scheme)
wipeAnalysis;
constraints Transformation
numberer RCM
system SparseGEN
algorithm Linear
integrator CentralDifference;
analysis Transient

analyze [expr int(10.*4.096)] 0.001; # apply 10*4096 steps for 0.01-sec time steps in analysis

puts "Done!"
#wipe - I like to type this one myself... cuz if you do in script, you cannot check your node and other thing..the model's gone!

Of course the data I used are random values but the model does run and give result. Now you have to check the values and make sure they make sense.

[suusie]

Thanks for your comments,

I think your script gives me the same results of what I got from the script I wrote above (Mon May 21, 2012 5:16 am) - the nodal displacement goes blew up.

The point is that when I replace the analysis commands of explicit scheme:

constraints Transformation
numberer RCM
system SparseGEN
algorithm Linear
integrator CentralDifference;
analysis Transient

to the implicit scheme:

constraints Transformation;
numberer RCM;
system UmfPack;
test EnergyIncr 1.0e-10 100;
algorithm Newton;
integrator Newmark 0.5 0.25 ;
analysis Transient;

then the analysis by implicit scheme gives result of what I expect to get.
I want to get the results from the explicit scheme that is same as what I can get from the implicit scheme.

Again, my analysis time step is small enough and I changed the other values several times to see if they worked but those weren't working well.

If there is anything that can help this problem, I would appreciate!

[Jhno]

Well, I did use the following script and obtained the same results with both methods (they are with # so you can switch from one to the other).
I used a sin function to generate the acceleration (I went from 0 to 1 second with dt : 0.00001)
The Central difference is more sensible to any change in values so you must be carefull.
Clearly the data I use aren't relevant.. but the methods in itselfs seems to be working just fine.

wipe
model BasicBuilder -ndm 2 -ndf 3;


# define geometry
set m 825000000.

set pi 3.14159

# define nodal coordinates:
node 1 0. 0. ;
node 2 5.0 0. ;
node 3 0. 3.
node 4 5. 3.
mass 3 $m $m $m
mass 4 $m $m $m
# Single point constraints -- Boundary Conditions
fix 1 1 1 1;
fix 2 1 1 1;
equalDOF 3 4 1 ;# I added this one... but I doubt it is "necessary"

# Columns and Beam Properties

set E 200.

# define geometric transformation:

geomTransf Linear 1

# Define Elements----------------------

element elasticBeamColumn 1 1 3 1. $E 1.e1 1; # Columns
element elasticBeamColumn 2 2 4 1. $E 1.e1 1; # Columns
element elasticBeamColumn 3 3 4 1. $E 1.e1 1; # Beam


set w [expr sqrt([eigen fullGenLapack 1])] ;#my class on those are long done so from memorie...maybe not the good expression
set T [expr 2.*$pi/$w]
puts "Period : $T"


timeSeries Path 1 -dt 0.01 -filePath TAKY.th ;

# set damping based on first eigen mode
set freq $w
rayleigh [expr 2*0.02*$freq] 0. 0. 0.

pattern UniformExcitation 1 1 -accel 1;
recorder Node -file Noeud_2.out -nodeRange 1 4 -dof 1 disp


constraints Transformation;
numberer RCM;
system FullGeneral
test EnergyIncr 1.0e-10 100;
algorithm Newton;
integrator Newmark 0.5 [expr 1./6.] ;
analysis Transient

# create the analysis (explicit scheme)
#wipeAnalysis;
#constraints Transformation
#numberer RCM
#system FullGeneral
#algorithm Linear
#integrator CentralDifference;
#analysis Transient

analyze [expr int(10000.*4.096)] 0.00001; # apply 10*4096 steps for 0.01-sec time steps in analysis

puts "Done!"
#wipe

[suusie]

Thanks for your comment,

The script you wrote above maybe helpful, but I just want to know which part of my script is wrong (the script was written above - on Wed Jun 08, 2016 2:03 pm).
I am not interested in if OpenSees' explicit methods in itself work fine or not.

The input I used "TAKY.th" was the earthquake acceleration history.
As I wrote above, I could get reasonable results if I use implicit methods. Also I tried small analysis time step so stability criteria was met.

Can anyone help to solve this problem ?
Thanks in advance.

[Jhno]

The idea was to show the method does work in good circumstances.
Now, since your model fails, I see three explanation :
1- the model is ill defined
2-the data you used are not consistant
3-the resolution has issues (example : dt too large in the context).
Your problems looks correct and if you are confident with your data, I would kindly suggest to reduce the time step.
I did use the same exact script as your in your first post with one of my acc time history and used a smaller dt than yours and it worked (trial and error until it works).
I did not check the validity of the result thought.
Maybe I am wrong and someone will point you the real problem but it is my opinion. Hope it helps, good luck.