Simple model, pushover analysis problem

[Ricsey]

Dear All!

I have a lot of experiment about using opensees, because I'm using opensees for my thesis work modeling bridges with nonlinear RC columns for modal pushover analysis. I want to find my errors and I defined a very simple model to check what am I modeling wrong.
The model is:
- HEA200 steel column, 3.0 m height
- steel01 material model (S235, fy = 235 MPa, b = 0.01)
- linear transformation
- 1 concentrated force at the topmost node for pushover analysis
This simple model doesn't converge after some steps, and when i look in the output the results shows that the force - deformation curve is only in the LINEAR state.
I would be really appreciate if someone could check where did I go wrong, it would be very helpful.

MY MODEL:
#units kN, m, sec, tonna
wipe
model BasicBuilder -ndm 3 -ndf 6

########################
#Initial datas
set columnLength 3.
set division 10.
set elementLength [expr $columnLength / $division]
set elementNumbers $division

set g 9.81
#set gammaSteel 76.982
set rhoSteel 7.85
set fy 235e3
set Es 210e6
set slope 0.01
set EBig [expr $Es * 1.e4]

set h 0.19
set b 0.20
set tw 0.006
set tf 0.010
set areaHEA [expr 2*($b*$tf) + ($h-2*$tf)*$tw]

set nodalMass [expr $areaHEA*$rhoSteel*$elementLength]
set nodalMassHalf [expr $areaHEA*$rhoSteel/2*$elementLength]
######################################


set currentLength 0.0
for {set i 0} {$i <= $division} {incr i} {
if {$i == 0} {
node 1 0. 0. 0.
puts "node 1 0. 0. 0."
} else {
set currentLength [expr $currentLength + $elementLength]
set nodeNumber [expr $i + 1]
node $nodeNumber 0. 0. $currentLength
puts "node $nodeNumber 0. 0. $currentLength"
}
}
set controlNode $i
fix 1 1 1 1 1 1 1


uniaxialMaterial Steel01 1 $fy $Es $slope
section Fiber 2 {
set flangeTopY1 [expr -$b/2]
set flangeTopZ1 [expr $h/2 - $tf]
set flangeTopY2 [expr $b/2]
set flangeTopZ2 [expr $h/2]
set flangeBottomY1 [expr -$b/2]
set flangeBottomZ1 [expr -$h/2 + $tf]
set flangeBottomY2 [expr $b/2]
set flangeBottomZ2 [expr -$h/2]
set webY1 [expr -$tw/2]
set webZ1 [expr -$h/2 + $tf]
set webY2 [expr $tw/2]
set webZ2 [expr $h/2 - $tf]

patch rect 1 10 4 $flangeTopY1 $flangeTopZ1 $flangeTopY2 $flangeTopZ2
patch rect 1 10 4 $flangeBottomY1 $flangeBottomZ1 $flangeBottomY2 $flangeBottomZ2
patch rect 1 4 10 $webY1 $webZ1 $webY2 $webZ2
}


uniaxialMaterial Elastic 10 $EBig
section Aggregator 1 10 T -section 2

################
#Defining masses
################

for {set nodeNumber 1} {$nodeNumber <= [expr $division+1]} {incr nodeNumber} {
if {$nodeNumber == 1 || $nodeNumber == $division} {
mass $nodeNumber $nodalMassHalf $nodalMassHalf $nodalMassHalf 1.e-9 1.e-9 1.e-9
puts "mass $nodeNumber $nodalMassHalf $nodalMassHalf $nodalMassHalf 1.e-9 1.e-9 1.e-9"
} else {
mass $nodeNumber $nodalMass $nodalMass $nodalMass 1.e-9 1.e-9 1.e-9
puts "mass $nodeNumber $nodalMass $nodalMass $nodalMass 1.e-9 1.e-9 1.e-9"
}
}

#################
#defining elements
#################
geomTransf PDelta 100 1 0 0
set integrationPoints 5
for {set elementNumber 1} {$elementNumber <= $division} {incr elementNumber} {
set nodeTagStart $elementNumber
set nodeTagEnd [expr $elementNumber + 1]
#element forceBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag <-mass $massDens> <-iter $maxIters $tol> <-integration $intType>
element dispBeamColumn $elementNumber $nodeTagStart $nodeTagEnd $integrationPoints 1 100
puts "element dispBeamColumn $elementNumber $nodeTagStart $nodeTagEnd $integrationPoints 1 100"
}


proc numfilter {strvar} {
set n [string length $strvar]
for {set i 0} {$i < $n} {incr i} {
set actnum [string index $strvar $i]
if {$actnum in {0 1 2 3 4 5 6 7 8 9 .}} {
append newvar $actnum
} else {
break
}
}

return $newvar
}

set NODES [getNodeTags]
set MASS {}
set GRAVITYLOAD {}
foreach node $NODES {
set actualMass [numfilter [nodeMass $node 1]]
lappend MASS $actualMass
lappend GRAVITYLOAD [expr $actualMass * $g]
}

pattern Plain 1 Linear {
foreach node $NODES gravityLoad $GRAVITYLOAD {
load $node 0.0 0.0 -$gravityLoad 0.0 0.0 0.0
}
}

constraints Plain
numberer RCM
system BandGeneral
test EnergyIncr 1.e-8 5
algorithm Newton
set NStepsGravity 10
set DGravity [expr 1./$NStepsGravity]
integrator LoadControl $DGravity
analysis Static
analyze $NStepsGravity
loadConst -time 0.0


##################
#Pushover analysis
##################

pattern Plain 100 Linear {
load $controlNode 1.0 0.0 0.0 0.0 0.0 0.0
}

recorder Node -file "baseShearForce.txt" -node 1 -dof 1 reaction
recorder Node -file "topDisplacement.txt" -node $controlNode -dof 1 disp

set displacementIncrement [expr $columnLength/1000]
set maximumDisplacement 1.0
integrator DisplacementControl $controlNode 1 $displacementIncrement
set currentDisplacement [nodeDisp $controlNode 1]
set ok 0
set count 0
while {$ok == 0 && abs($currentDisplacement) < abs($maximumDisplacement)} {
set ok [analyze 1]
set currentDisp [nodeDisp $controlNode 1]
if {$ok != 0} {
puts "****************analyze failed at $currentDisp"
puts "****************Trying Newton with INITIAL stiffness ..****************"
test NormDispIncr 1.0e-7 500
algorithm ModifiedNewton -initial
set ok [analyze 1]
if {$ok == 0} {puts "that worked .. back to regular newton"}
test NormDispIncr 1.0e-8 10
algorithm Newton
}
if {$ok != 0} {
puts "****************analyze failed at $currentDisp"
puts "****************Trying Broyden ..****************"
algorithm Broyden 8
set ok [analyze 1 ]
if {$ok == 0} {puts "****************that worked .. back to regular newton"}
algorithm Newton
}
if {$ok != 0} {
puts "****************analyze failed at $currentDisp"
puts "****************Trying NewtonWithLineSearch ..****************"
algorithm NewtonLineSearch 0.8
set ok [analyze 1]
if {$ok == 0} {puts "****************that worked .. back to regular newton"}
algorithm Newton
}
set maximumStep [expr int($maximumDisplacement/$displacementIncrement)]
set currentStepPercentage [expr $currentDisplacement/$maximumDisplacement*100]
puts [format "$count/$maximumStep step %0.2f percent, %0.6f OK" $currentStepPercentage $currentDisp]
set count [expr $count + 1]
}
wipe

[Ricsey]

I really need an answer pls!

[mhscott]

Try with linear geometric transformation instead of P-Delta. You could be losing geometric stiffness before material yielding

[Ricsey]

Thank you for your reply.
I made it right to get the correct results. I had to decrease the displacement increment if the solution is not converge.
This is an interesting thing. If I want to model the P-Delta effect, I HAD to set the transformation to PDelta. In this case I can lose the geometric stiffness so I cant get the right results. Am I thinking good?

[fmk]

if you want to model PDelta, use PDelta. If you don't, use Linear. Setting the transformation to PDelta to get PDelta is NOT an '"interesting thing", it should be an "obvious thing"! Decreasing the step size when things don't converge is what other black box programs do automatically and should be in everyones arsenal of things to try is the solution does not converge.