Lack of convergence in linear analysis ?

[Jean Domínguez]

Hi. I´m currently doing an analysis using elastic section and nonlinearBeamColumn element with only one load case (gravity). I´m comparing my OpenSees results with other obtained from SAP2000 with exactly the same model. First, i did used linear analysis in 10 steps (Analyze 10 and Integrator 0.1 like most examples) and the support reactions didn´t match, next i repeated with 100 steps (analyze 100 and integrator 0.01) and the reactions become similar to those from SAP2000, i was smiled until my next attempt with 500 steps hoping a closer results with the surprising diverging reactions. There´s any problem with the element type ?, can i post my .tcl to get some help ?

[Jean Domínguez]

I did repeated the same model with elasticBeamColumn element and seems to be fine, the results match those from SAP2000 and they are not number of steps depending.
Theoretically the nonlinearbeamcolumn must be the same as elasticbeamcolumn element if the elastic section is used, can anyone help me with this ?

[fmk]

why don't you post your model .. they should give the same results.

[Jean Domínguez]

Ok, this is the model. with the proper change of comments (#) one can change the type of element.

# unidades Kgf-cm
# SET UP ----------------------------------------------------------------------------
wipe;
model basic -ndm 2 -ndf 3;
file mkdir "3 niveles"

# GEOMETRIA -------------------------------------------------------------
# Coordenadas nodales:
#node 1 0 0
node 2 0 285
node 3 0 570
node 4 0 855
node 5 200 0
node 6 200 285
node 7 200 570
node 8 200 855
node 9 800 0
node 10 800 285
node 11 800 570
node 12 800 855
node 13 1250 0
node 14 1250 285
node 15 1250 570
node 16 1250 855
node 17 1700 0
node 18 1700 285
node 19 1700 570
node 20 1700 855
node 21 2225 0
node 22 2225 285
node 23 2225 570
node 24 2225 855

# Restricciones

fix 5 1 1 1
fix 9 1 1 1
fix 13 1 1 1
fix 17 1 1 1
fix 21 1 1 1

# Constraints

# Masas nodales: OJO EL SOLVER SOLO PUEDE RESOLVER N-1 MODOS

mass 6 109.44 1.0 0.0
mass 7 109.44 1.0 0.0
mass 8 82.31 1.0 0.0


# Secciones de columnas
# E A Izz

#section Elastic 1 238752 900 67500; # C30X30
#section Elastic 2 238752 1500 312500; # C30X50
#section Elastic 6 238752 1500 112500; # C50X30

# Secciones de vigas

#section Elastic 3 238752 1200 160000; # V30X40
#section Elastic 4 238752 1650 415937.5; # V30X55
#section Elastic 5 238752 1500 78125; # VP25X60

# Elementos -------------------------------------------------------------

geomTransf Linear 1

# conectividad:

# Columnas ----------------------------------------------------------------

#element nonlinearBeamColumn 1 5 6 5 2 1;# ojo con la transformacion y masa
#element nonlinearBeamColumn 2 6 7 5 2 1
#element nonlinearBeamColumn 3 7 8 5 2 1

element elasticBeamColumn 1 5 6 1500 238752 312500 1
element elasticBeamColumn 2 6 7 1500 238752 312500 1
element elasticBeamColumn 3 7 8 1500 238752 312500 1

#element nonlinearBeamColumn 4 9 10 5 2 1
#element nonlinearBeamColumn 5 10 11 5 2 1
#element nonlinearBeamColumn 6 11 12 5 2 1

element elasticBeamColumn 4 9 10 1500 238752 312500 1
element elasticBeamColumn 5 10 11 1500 238752 312500 1
element elasticBeamColumn 6 11 12 1500 238752 312500 1


#element nonlinearBeamColumn 7 13 14 5 1 1
#element nonlinearBeamColumn 8 14 15 5 1 1
#element nonlinearBeamColumn 9 15 16 5 1 1

element elasticBeamColumn 7 13 14 900 238752 67500 1
element elasticBeamColumn 8 14 15 900 238752 67500 1
element elasticBeamColumn 9 15 16 900 238752 67500 1


#element nonlinearBeamColumn 10 17 18 5 6 1
#element nonlinearBeamColumn 11 18 19 5 6 1
#element nonlinearBeamColumn 12 19 20 5 6 1

element elasticBeamColumn 10 17 18 1500 238752 112500 1
element elasticBeamColumn 11 18 19 1500 238752 112500 1
element elasticBeamColumn 12 19 20 1500 238752 112500 1

#element nonlinearBeamColumn 13 21 22 5 6 1
#element nonlinearBeamColumn 14 22 23 5 6 1
#element nonlinearBeamColumn 15 23 24 5 6 1

element elasticBeamColumn 13 21 22 1500 238752 112500 1
element elasticBeamColumn 14 22 23 1500 238752 112500 1
element elasticBeamColumn 15 23 24 1500 238752 112500 1

# vigas -------------------------------------------------------------------

#element nonlinearBeamColumn 16 2 6 5 4 1
#element nonlinearBeamColumn 17 6 10 5 3 1
#element nonlinearBeamColumn 18 10 14 5 5 1
#element nonlinearBeamColumn 19 14 18 5 3 1
#element nonlinearBeamColumn 20 18 22 5 5 1

element elasticBeamColumn 16 2 6 1650 238752 415937.5 1
element elasticBeamColumn 17 6 10 1200 238752 160000 1
element elasticBeamColumn 18 10 14 1500 238752 78125 1
element elasticBeamColumn 19 14 18 1200 238752 160000 1
element elasticBeamColumn 20 18 22 1500 238752 78125 1

#element nonlinearBeamColumn 21 3 7 5 4 1
#element nonlinearBeamColumn 22 7 11 5 3 1
#element nonlinearBeamColumn 23 11 15 5 5 1
#element nonlinearBeamColumn 24 15 19 5 3 1
#element nonlinearBeamColumn 25 19 23 5 5 1

element elasticBeamColumn 21 3 7 1650 238752 415937.5 1
element elasticBeamColumn 22 7 11 1200 238752 160000 1
element elasticBeamColumn 23 11 15 1500 238752 78125 1
element elasticBeamColumn 24 15 19 1200 238752 160000 1
element elasticBeamColumn 25 19 23 1500 238752 78125 1

#element nonlinearBeamColumn 26 4 8 5 3 1
#element nonlinearBeamColumn 27 8 12 5 3 1
#element nonlinearBeamColumn 28 12 16 5 5 1
#element nonlinearBeamColumn 29 16 20 5 3 1
#element nonlinearBeamColumn 30 20 24 5 5 1

element elasticBeamColumn 26 4 8 1200 238752 160000 1
element elasticBeamColumn 27 8 12 1200 238752 160000 1
element elasticBeamColumn 28 12 16 1500 238752 78125 1
element elasticBeamColumn 29 16 20 1200 238752 160000 1
element elasticBeamColumn 30 20 24 1500 238752 78125 1

# RECORDERS -------------------------------------------------------------

recorder Node -file "3 niveles/3N82Cnode.out" -dof 1 2 3 disp; # displacements of free nodes
recorder Node -file "3 niveles/3N82Creact.out" -node 5 9 13 17 21 -dof 1 2 reaction; # support reaction
recorder Node -file "3 niveles/3N82Ceig1" -dof 1 2 3 "eigen 1"; # Modes 1 2 & 3
recorder Node -file "3 niveles/3N82Ceig2" -dof 1 2 3 "eigen 2";
recorder Node -file "3 niveles/3N82Ceig3" -dof 1 2 3 "eigen 3";

# CARGAS DE GRAVEDAD -------------------------------------------------------------

# Cargas distribuidas sobre vigas -------------------------------------------------------

pattern Plain 100 Linear {
eleLoad -ele 16 21 -type -beamUniform -38.962
eleLoad -ele 26 27 29 -type -beamUniform -29.276
eleLoad -ele 17 22 19 24 -type -beamUniform -36.197
eleLoad -ele 18 23 20 25 -type -beamUniform -36.002
eleLoad -ele 28 30 -type -beamUniform -29.081

# Cargas puntuales sobre vigas ---------------------------------------------------------

eleLoad -ele 16 21 -type beamPoint -3139.50 0

# Cargas puntuales sobre nodos ----------------------------------------------------------

load 6 0.0 -5358.557 0.0
load 7 0.0 -5358.557 0.0
load 8 0.0 -4563.412 0.0
load 10 0.0 -5284.728 0.0
load 11 0.0 -5284.728 0.0
load 12 0.0 -4563.412 0.0
load 14 0.0 -4857.228 0.0
load 15 0.0 -4857.228 0.0
load 16 0.0 -4135.912 0.0
load 18 0.0 -5284.728 0.0
load 19 0.0 -5284.728 0.0
load 20 0.0 -4563.412 0.0
load 22 0.0 -6153.781 0.0
load 23 0.0 -6153.781 0.0
load 24 0.0 -3882.881 0.0
}

set lambda [eigen 3]
set lambda1 [lindex $lambda 0]
set lambda2 [lindex $lambda 1]
set lambda3 [lindex $lambda 2]
set omega1 [expr pow($lambda1,0.5)]
set omega2 [expr pow($lambda2,0.5)]
set omega3 [expr pow($lambda3,0.5)]
set T1 [expr 2*3.14159/$omega1]
set T2 [expr 2*3.14159/$omega2]
set T3 [expr 2*3.14159/$omega3]

set filename "3 niveles/periodos3N82C.out"
set fileid [open $filename "w"]
puts $fileid {T1 T2 T3}
puts $fileid $T1
puts $fileid $T2
puts $fileid $T3
close $fileid

constraints Plain; # how it handles boundary conditions
numberer Plain; # renumber dof's to minimize band-width (optimization), if you want to
system BandGeneral; # how to store and solve the system of equations in the analysis
test NormDispIncr 1.0e-8 6; # determine if convergence has been achieved at the end of an iteration step
algorithm Linear; # use Newton's solution algorithm: updates tangent stiffness at every iteration
integrator LoadControl 0.1; # determine the next time step for an analysis, # apply gravity in 10 steps
analysis Static; # define type of analysis static or transient
analyze 10; # perform gravity analysis
#loadConst; # hold gravity constant and restart time
print "3 niveles/3N82C.out"

[Jean Domínguez]

The reactions from SAP2000 are:
node 5: fx:109.14 fy:75398.04
node 9: fx:-55.10 fy:68080.75
node 13: fx:412.43 fy:58492.12
node 17: fx:905.46 fy:65504.16
node 21: fx:1371.93 fy:41144.25

With elasticBeamColumn the reactions are right and the node displacement with a 1% error (normal?)
The results don´t change with number of step change.

With nonlinearBeamColumn the reactions diverge and change if a change in step is made.

Thanks FMK

[Jean Domínguez]

Ok, finally it converges...., but for nonlinearBeamColumn i need to spend more than 1000 steps while with elasticBeamColumn i get decent results with 1 step... is it normal?

How much amount error i must expect from this comparison ?

[fmk]

the problem is that the way eleLoads are applied .. you need to iterate to achieve convergence and when using the nonlinearBeamColumn (forceBeamColumn) you will need 1 more iteration .. you were actually not iterating because despite spec a convergence test with up tp 1000 iterations you were using a Linear solution algorithm.


here are results which shows the results for 3 basic beam column elements (elastic, dispBeamColumn and forceBEamColumn) .. you will notice the forceBeamColumn equires that 1 extra iteration:

CTestNormDispIncr::test() - iteration: 1 current Norm: 0.135728 (max: 1e-08, Norm deltaR: 1.3893e-10)
CTestNormDispIncr::test() - iteration: 2 current Norm: 4.8094e-15 (max: 1e-08, Norm deltaR: 1.89248e-10)

eleType: elastic reactions: 109.154 75397.4 -55.1077 68080.4 412.424 58491.8 905.457 65503.8 -1371.93 41144.1

CTestNormDispIncr::test() - iteration: 1 current Norm: 0.146498 (max: 1e-08, Norm deltaR: 224189)
CTestNormDispIncr::test() - iteration: 2 current Norm: 0.0227855 (max: 1e-08, Norm deltaR: 9.20576e-11)
CTestNormDispIncr::test() - iteration: 3 current Norm: 3.38659e-16 (max: 1e-08, Norm deltaR: 3.28663e-10)

eleType: forceBeamColumn reactions: 109.154 75397.4 -55.1077 68080.4 412.424 58491.8 905.457 65503.8 -1371.93 41144.1

CTestNormDispIncr::test() - iteration: 1 current Norm: 0.135728 (max: 1e-08, Norm deltaR: 1.6039e-10)
CTestNormDispIncr::test() - iteration: 2 current Norm: 1.01972e-14 (max: 1e-08, Norm deltaR: 1.44851e-10)

eleType: dispBeamColumn reactions: 109.154 75397.4 -55.1077 68080.4 412.424 58491.8 905.457 65503.8 -1371.93 41144.1


here is the script. note that instead of using comments you can do multiple analysis by putting them inside the foreach loop:

[code]
foreach eleType {elastic forceBeamColumn dispBeamColumn} {

wipe;
model basic -ndm 2 -ndf 3;

# GEOMETRIA -------------------------------------------------------------
# Coordenadas nodales:
#node 1 0 0
node 2 0 285
node 3 0 570
node 4 0 855
node 5 200 0
node 6 200 285
node 7 200 570
node 8 200 855
node 9 800 0
node 10 800 285
node 11 800 570
node 12 800 855
node 13 1250 0
node 14 1250 285
node 15 1250 570
node 16 1250 855
node 17 1700 0
node 18 1700 285
node 19 1700 570
node 20 1700 855
node 21 2225 0
node 22 2225 285
node 23 2225 570
node 24 2225 855

# Restricciones

fix 5 1 1 1
fix 9 1 1 1
fix 13 1 1 1
fix 17 1 1 1
fix 21 1 1 1

# Constraints

# Masas nodales: OJO EL SOLVER SOLO PUEDE RESOLVER N-1 MODOS

mass 6 109.44 1.0 0.0
mass 7 109.44 1.0 0.0
mass 8 82.31 1.0 0.0


# Secciones de columnas
# E A Izz

section Elastic 1 238752 900 67500; # C30X30
section Elastic 2 238752 1500 312500; # C30X50
section Elastic 6 238752 1500 112500; # C50X30

# Secciones de vigas

section Elastic 3 238752 1200 160000; # V30X40
section Elastic 4 238752 1650 415937.5; # V30X55
section Elastic 5 238752 1500 78125; # VP25X60

# Elementos -------------------------------------------------------------

geomTransf Linear 1

if {$eleType != "elastic"} {
element $eleType 1 5 6 5 2 1;# ojo con la transformacion y masa
element $eleType 2 6 7 5 2 1
element $eleType 3 7 8 5 2 1

element $eleType 4 9 10 5 2 1
element $eleType 5 10 11 5 2 1
element $eleType 6 11 12 5 2 1

element $eleType 7 13 14 5 1 1
element $eleType 8 14 15 5 1 1
element $eleType 9 15 16 5 1 1

element $eleType 10 17 18 5 6 1
element $eleType 11 18 19 5 6 1
element $eleType 12 19 20 5 6 1

element $eleType 13 21 22 5 6 1
element $eleType 14 22 23 5 6 1
element $eleType 15 23 24 5 6 1

element $eleType 16 2 6 5 4 1
element $eleType 17 6 10 5 3 1
element $eleType 18 10 14 5 5 1
element $eleType 19 14 18 5 3 1
element $eleType 20 18 22 5 5 1

element $eleType 21 3 7 5 4 1
element $eleType 22 7 11 5 3 1
element $eleType 23 11 15 5 5 1
element $eleType 24 15 19 5 3 1
element $eleType 25 19 23 5 5 1
element $eleType 26 4 8 5 3 1
element $eleType 27 8 12 5 3 1
element $eleType 28 12 16 5 5 1
element $eleType 29 16 20 5 3 1
element $eleType 30 20 24 5 5 1

} else {
element elasticBeamColumn 1 5 6 1500 238752 312500 1
element elasticBeamColumn 2 6 7 1500 238752 312500 1
element elasticBeamColumn 3 7 8 1500 238752 312500 1

element elasticBeamColumn 4 9 10 1500 238752 312500 1
element elasticBeamColumn 5 10 11 1500 238752 312500 1
element elasticBeamColumn 6 11 12 1500 238752 312500 1

element elasticBeamColumn 7 13 14 900 238752 67500 1
element elasticBeamColumn 8 14 15 900 238752 67500 1
element elasticBeamColumn 9 15 16 900 238752 67500 1

element elasticBeamColumn 10 17 18 1500 238752 112500 1
element elasticBeamColumn 11 18 19 1500 238752 112500 1
element elasticBeamColumn 12 19 20 1500 238752 112500 1

element elasticBeamColumn 13 21 22 1500 238752 112500 1
element elasticBeamColumn 14 22 23 1500 238752 112500 1
element elasticBeamColumn 15 23 24 1500 238752 112500 1

# vigas -------------------------------------------------------------------


element elasticBeamColumn 16 2 6 1650 238752 415937.5 1
element elasticBeamColumn 17 6 10 1200 238752 160000 1
element elasticBeamColumn 18 10 14 1500 238752 78125 1
element elasticBeamColumn 19 14 18 1200 238752 160000 1
element elasticBeamColumn 20 18 22 1500 238752 78125 1
element elasticBeamColumn 21 3 7 1650 238752 415937.5 1
element elasticBeamColumn 22 7 11 1200 238752 160000 1
element elasticBeamColumn 23 11 15 1500 238752 78125 1
element elasticBeamColumn 24 15 19 1200 238752 160000 1
element elasticBeamColumn 25 19 23 1500 238752 78125 1

element elasticBeamColumn 26 4 8 1200 238752 160000 1
element elasticBeamColumn 27 8 12 1200 238752 160000 1
element elasticBeamColumn 28 12 16 1500 238752 78125 1
element elasticBeamColumn 29 16 20 1200 238752 160000 1
element elasticBeamColumn 30 20 24 1500 238752 78125 1
}
# RECORDERS -------------------------------------------------------------

recorder Node -file reactions.out -node 5 9 13 17 21 -dof 1 2 reaction; # support reaction

# CARGAS DE GRAVEDAD -------------------------------------------------------------

# Cargas distribuidas sobre vigas -------------------------------------------------------

pattern Plain 100 Linear {
eleLoad -ele 16 21 -type -beamUniform -38.962
eleLoad -ele 26 27 29 -type -beamUniform -29.276
eleLoad -ele 17 22 19 24 -type -beamUniform -36.197
eleLoad -ele 18 23 20 25 -type -beamUniform -36.002
eleLoad -ele 28 30 -type -beamUniform -29.081

# Cargas puntuales sobre vigas ---------------------------------------------------------

eleLoad -ele 16 21 -type beamPoint -3139.50 0

# Cargas puntuales sobre nodos ----------------------------------------------------------
eleLoad -ele 26 27 29 -type -beamUniform -29.276
eleLoad -ele 17 22 19 24 -type -beamUniform -36.197
eleLoad -ele 18 23 20 25 -type -beamUniform -36.002
eleLoad -ele 28 30 -type -beamUniform -29.081

# Cargas puntuales sobre vigas ---------------------------------------------------------

eleLoad -ele 16 21 -type beamPoint -3139.50 0

# Cargas puntuales sobre nodos ----------------------------------------------------------

load 6 0.0 -5358.557 0.0
load 7 0.0 -5358.557 0.0
load 8 0.0 -4563.412 0.0
load 10 0.0 -5284.728 0.0
load 11 0.0 -5284.728 0.0
load 12 0.0 -4563.412 0.0
load 14 0.0 -4857.228 0.0
load 15 0.0 -4857.228 0.0
load 16 0.0 -4135.912 0.0
load 18 0.0 -5284.728 0.0
load 19 0.0 -5284.728 0.0
load 20 0.0 -4563.412 0.0
load 22 0.0 -6153.781 0.0
load 23 0.0 -6153.781 0.0
load 24 0.0 -3882.881 0.0
}
set lambda [eigen 3]
set lambda1 [lindex $lambda 0]
set lambda2 [lindex $lambda 1]
set lambda3 [lindex $lambda 2]
set omega1 [expr pow($lambda1,0.5)]
set omega2 [expr pow($lambda2,0.5)]
set omega3 [expr pow($lambda3,0.5)]
set T1 [expr 2*3.14159/$omega1]
set T2 [expr 2*3.14159/$omega2]
set T3 [expr 2*3.14159/$omega3]

set filename "periodos3N82C.out"
set fileid [open $filename "w"]
puts $fileid {T1 T2 T3}
puts $fileid $T1
puts $fileid $T2
puts $fileid $T3
close $fileid

constraints Plain;
numberer Plain;
system BandGeneral;
test NormDispIncr 1.0e-8 3 1;
algorithm Newton;
integrator LoadControl 0.1;
analysis Static;
analyze 1;


test NormDispIncr 1.0e-8 3 0;
analyze 9
wipe

set fileID [open reactions.out r]
set lines [split [read $fileID] \n]
puts "eleType: $eleType reactions: [lindex $lines 9]\n"
close $fileID

}



[/code]

[Jean Domínguez]

Thanks FMK, very useful reply. I did changed the algorithm to Newton and achieved good results, and the forearch commands is very nice too, i´m learning tcl yet.
Where my problem exactly was? only the type of algorithm? or something else?. I don´t understand very much the twice use of test and analyze in the corrected script, isn´t the same to use a unique analyze 10?.

[fmk]

the problem was the force beam column required the extra iteration in a step which the Linear algoritm did not provide.

the use of the two tests is to demonstrate this in the output .. the first test has a 1 at the end which tells the test to print some info each time test() is called. the second test does not have it. i do one analysis with this printing on and then 9 with it off. know you understand you can remove the first step and analyze and change he second analyze to do 10 analysis steps.

the optional 1 and sometimes 2 at the end of the tests are useful to see how the model is converging if you have problems.

[Jean Domínguez]

FMK, once again, thanks... very useful comments, i have another issues but i will try to solve it myself, if i can´t i will ask for help again .