Very large displacement result of timehistory analysis

[andrewfirst]

I'm doing nonlinear transient analysis of a simple one story-three bay RC frame. The lateral displacement result runs quickly to almost 100% storey height during first steps, which is very weird. I've checked the code and done the moment-curture analysis of column. The eigenvalue analysis is correct, so does the moment-curture analysis. This problem has tortured me for several days.

[silvia]

have you compared it to the examples manual files?

[andrewfirst]

Yes, I've checked it for many times and can't found what's wrong. Even if I scaled the input acceleration to very low, the displacement result is still too large. How weird!

[silvia]

does it scale proportionally?
if you scale the gm by 0.001 does the displacement scale accordingly?

[andrewfirst]

No. It seems that the scaling of ground motion does not influence the lateral displacement! Just in about 0.2s the lateral displacement runs to very high values (e.g. several times of storey height)

[silvia]

check your model, or test your procedure on a model from the examples manual.

[andrewfirst]

This is the full code. Anyone help me?
# Reinforced concrete multi-bay, two-story frame
# Distributed vertical load on girder
# Displacement-based beam-column elements
# Gravity load analysis followed by transient analysis
# Units: kg,m,sec ->> N,Pa

wipe

#--------------------------
# define GEOMETRY variables
#--------------------------
set pi [expr {2*asin(1.0)}]
set Bcol 500.e-3; # column section width
set Hcol 500.e-3; # column section height
set covcol 30.e-3; # column cover thickness = between outer surfaces of concrete and longi rebar
set Bbeam 250.e-3; # beam section width
set Hbeam 500.e-3; # beam section height
set covbeam 25.e-3; # beam cover thickness = between outer surfaces of concrete and longi rebar
set Bmbeam 250.e-3; # middle beam width
set Hmbeam 300.e-3; # middle beam height
# define REINFORCEMENT variables
set Dcolbar 25.e-3; # diameter of column longitudinal rebar
set Acolbar [expr {$pi * $Dcolbar * $Dcolbar / 4.0}]; # single bar area of column longi rebar
set Dbeambar 18.e-3; # diameter of beam longitudinal rebar
set Abeambar [expr {$pi * $Dbeambar * $Dbeambar / 4.0}]; # single bar area of beam longi rebar
# define DAMPING variables from $xDamp use stiffness proportional damping only
set xDamp 0.05; # modal damping ratio
# ------ set analysis variables
set gamma 0.5; # gamma value for newmark integration
set beta 0.25; # beta value for newmark integration
set DtAnalysis 0.005; # time-step Dt for lateral analysis
set DtGround 0.02; # time-step Dt for input grond motion
set TmaxGround [expr {832 * $DtGround}]; # maximum duration of ground-motion analysis

# ------------------------------
# Start of model generation
# ------------------------------


# Create ModelBuilder (with two-dimensions and 3 DOF/node)

model BasicBuilder -ndm 2 -ndf 3

# Create nodes

# tag x y
node 1 0.0 0.0
node 2 6000.e-3 0.0
node 3 8400.e-3 0.0
node 4 14400.e-3 0.0
node 5 0. 4200.e-3
node 6 6000.e-3 4200.e-3
node 7 8400.e-3 4200.e-3
node 8 14400.e-3 4200.e-3

# Define nodal mass
# tag X Y RZ
mass 5 11.8e3 11.8e3 0.
mass 6 10.8e3 10.8e3 0.
mass 7 10.8e3 10.8e3 0.
mass 8 11.8e3 11.8e3 0.

puts "Node coordinates and mass defined"

# Fix supports at base of columns
# node DX DY RZ
fix 1 1 1 1
fix 2 1 1 1
fix 3 1 1 1
fix 4 1 1 1
puts "Node fixed"

# Define materials for nonlinear columns
# CONCRETE
# Confined concrete:
set fc -24.e6; # CONCRETE Compressive Strength, ksi (+Tension, -Compression)
set Ec 30000.e6; # Concrete Elastic Modulus
set fc1C [expr {1.26394 * $fc}]; # CONFINED concrete (mander model), maximum stress
set eps1C [expr {2. * $fc1C / $Ec}]; # strain at maximum stress
set fc2C [expr {1.0 * $fc}]; # ultimate stress
set eps2C [expr {5. * $eps1C}]; # strain at ultimate stress
# Unconfined concrete:
set fc1U $fc; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U -0.002; # strain at maximum stress
set fc2U [expr {0.1 * $fc}]; # ultimate stress
set eps2U -0.006; # strain at ultimate stress
# Concrete02 variables:
set lambda 0.1 ; # ratio between unloading slope at $epscu and initial slope
set ftC [expr {-$fc1C / 10.}]; # tensile strength +tension confined
set ftU [expr {-$fc1U / 10.}]; # tensile strength +tension unconfined
set Ets [expr {$Ec / 10.}]; # tension softening stiffness
# reinforcing steel
set Fy 335.e6; # STEEL yield stress
set Es 200000.e6; # modulus of steel
set epsY [expr {$Fy / $Es}]; # steel yield strain
set Fu 469.e6; # ultimate stress of steel
set epsU 0.1; # ultimate strain of steel
set E2 [expr {($Fu-$Fy)/($epsU-$epsY)}]; # post-yield tangent stiffness
set Bs [expr {$E2 / $Es}]; # post-yield stiffness ratio of steel
# set up parameters for column section and element definition
set IDcore 1; # ID tag for core concrete
set IDcover 2; # ID tag for cover concrete
set IDsteel 3; # ID tag for steel
uniaxialMaterial Concrete02 $IDcore $fc1C $eps1C $fc2C $eps2C $lambda $ftC $Ets; # CORE CONCRETE
uniaxialMaterial Concrete02 $IDcover $fc1U $eps1U $fc2U $eps2U $lambda $ftU $Ets; # COVER CONCRETE
uniaxialMaterial Steel01 $IDsteel $Fy $Es $Bs; # REINFORCING STEEL

puts "Material defined"


# Define cross-section for columns

# column section

set coreiy [expr {$Hcol/2. - $covcol - $Dcolbar/2.}]; # concrete core point i y
set coreiz [expr {$Bcol/2. - $covcol - $Dcolbar/2.}]; # concrete core point i z

# core, bottom, top, left, right
section Fiber 1 {

# mat nfIJ nfJK yI zI yK zK
patch rect $IDcore 15 1 $coreiy $coreiz [expr {-$coreiy}] [expr {-$coreiz}]

patch rect $IDcover 10 1 [expr {$Hcol/2.}] [expr {$Bcol/2.}] [expr {-$Hcol/2.}] $coreiz

patch rect $IDcover 10 1 [expr {$Hcol/2.}] [expr {-$coreiz}] [expr {-$Hcol/2.}] [expr {-$Bcol/2.}]

patch rect $IDcover 2 1 [expr {$Hcol/2.}] $coreiz $coreiy [expr {-$coreiz}]

patch rect $IDcover 2 1 [expr {-$coreiy}] $coreiz [expr {-$Hcol/2.}] [expr {-$coreiz}]

# mat nBars area yStart zStart yEnd zEnd
layer straight $IDsteel 3 $Acolbar $coreiy $coreiz $coreiy [expr {-$coreiz}]
layer straight $IDsteel 2 $Acolbar 0. $coreiz 0. [expr {-$coreiz}]
layer straight $IDsteel 3 $Acolbar [expr {-$coreiy}] $coreiz [expr {-$coreiy}] [expr {-$coreiz}]
}
puts "column section defined"

# exterior Beam section
set coreiy [expr {-($Hbeam/2. - $covbeam - $Dbeambar/2.)}]; # concrete core point i y
set coreiz [expr {$Bbeam/2. - $covbeam - $Dbeambar/2.}]; # concrete core point i z

# core, left, right, bottom, top cover
section Fiber 2 {

# mat nfIJ nfJK yI zI yK zK
patch rect $IDcore 1 15 $coreiy $coreiz [expr {-$coreiy}] [expr {-$coreiz}]

patch rect $IDcover 1 10 [expr {-$Hbeam/2.}] [expr {$Bbeam/2.}] [expr {$Hbeam/2.}] $coreiz

patch rect $IDcover 1 10 [expr {-$Hbeam/2.}] [expr {-$coreiz}] [expr {$Hbeam/2.}] [expr {-$Bbeam/2.}]

patch rect $IDcover 1 2 [expr {-$Hbeam/2.}] $coreiz $coreiy [expr {-$coreiz}]

patch rect $IDcover 1 2 [expr {-$coreiy}] $coreiz [expr {$Hbeam/2.}] [expr {-$coreiz}]

# mat nBars area yStart zStart yEnd zEnd
layer straight $IDsteel 3 $Abeambar $coreiy $coreiz $coreiy [expr {-$coreiz}]
layer straight $IDsteel 3 $Abeambar [expr {-$coreiy}] $coreiz [expr {-$coreiy}] [expr {-$coreiz}]
}
puts "Exterior beam section defined"

# middle Beam section
set coreiy [expr {-($Hmbeam/2. - $covbeam - $Dbeambar/2.)}]; # concrete core point i y
set coreiz [expr {$Bmbeam/2. - $covbeam - $Dbeambar/2.}]; # concrete core point i z

# core, left, right, bottom, top cover
section Fiber 3 {

# mat nfIJ nfJK yI zI yK zK
patch rect $IDcore 1 10 $coreiy $coreiz [expr {-$coreiy}] [expr {-$coreiz}]

patch rect $IDcover 1 10 [expr {-$Hmbeam/2.}] [expr {$Bmbeam/2.}] [expr {$Hmbeam/2.}] $coreiz

patch rect $IDcover 1 10 [expr {-$Hmbeam/2.}] [expr {-$coreiz}] [expr {$Hmbeam/2.}] [expr {-$Bmbeam/2.}]

patch rect $IDcover 1 2 [expr {-$Hmbeam/2.}] $coreiz $coreiy [expr {-$coreiz}]

patch rect $IDcover 1 2 [expr {-$coreiy}] $coreiz [expr {$Hmbeam/2.}] [expr {-$coreiz}]

# mat nBars area yStart zStart yEnd zEnd
layer straight $IDsteel 3 $Abeambar $coreiy $coreiz $coreiy [expr {-$coreiz}]
layer straight $IDsteel 3 $Abeambar [expr {-$coreiy}] $coreiz [expr {-$coreiy}] [expr {-$coreiz}]
}
puts "Middle beam section defined"

# Number of integration points
set np 5

# Geometric transformation column
geomTransf Linear 1

# Define elements
# tag iNode jNode intgpts section transf
element dispBeamColumn 1 1 5 $np 1 1
element dispBeamColumn 2 2 6 $np 1 1
element dispBeamColumn 3 3 7 $np 1 1
element dispBeamColumn 4 4 8 $np 1 1

# define beam
# Geometric transformation beam
geomTransf Linear 2

element dispBeamColumn 5 5 6 $np 2 2
element dispBeamColumn 6 6 7 $np 3 2
element dispBeamColumn 7 7 8 $np 2 2

puts "Element defined"

set omega2 [eigen 1]
puts "Before gravity analysis Tn = [expr {2 * $pi / pow($omega2,0.5)}]"

# set recorders
recorder Node -file Node5.out -time -node 5 -dof 1 3 disp
recorder Element -file Element1.out -time -ele 1 force
recorder Element -file ele5steel.out -time -ele 5 section 1 fiber 216.e-3 91.e-3 $IDsteel stressStrain

# Define gravity loads
# Constant gravity load

# Create a Plain load pattern with a Linear TimeSeries
pattern Plain 1 Linear {

# Create nodal loads at nodes
load 5 0.0 -75250. 0.0
load 6 0.0 -61750. 0.0
load 7 0.0 -61750. 0.0
load 8 0.0 -75250. 0.0

# create unifor-distributed load on beams
# tag1 tag2 Wz
eleLoad -ele 5 7 -type -beamUniform -13.6e3
eleLoad -ele 6 -type -beamUniform -3.1e3
}
puts "Gravity load exerted"

# ------------------------------
# End of model generation
# ------------------------------
# Start of analysis generation for gravity analysis
# -------------------------------------------------
# Create the convergence test, the norm of the residual with a tolerance of

# 1e-8 and a max number of iterations of 10

test NormDispIncr 1.0e-8 10 0

# Create the solution algorithm, a Newton-Raphson algorithm
algorithm Newton

# Create the integration scheme, the LoadControl scheme using steps of 0.1

integrator LoadControl 0.1 1 0.1 0.1

# Create the system of equation, a SPD using a profile storage scheme
system BandGeneral

# Create the DOF numberer, the reverse Cuthill-McKee algorithm
numberer Plain

# Create the constraint handler, the transformation method
constraints Plain

# Create the analysis object
analysis Static

puts "Static analysis generated"
# ------------------------------------------------
# End of analysis generation for gravity analysis
# -------------------------------------------------

# ------------------------------
# Perform gravity load analysis
# ------------------------------
# initialize the model, done to set initial tangent

initialize


# perform the gravity load analysis, requires 10 steps to reach the load level

analyze 10

# set gravity loads to be const and set pseudo time to be 0.0
loadConst -time 0.0
puts "End of Gravity analysis"

# ------------------------------
# Define transient analysis
# ------------------------------
# create load pattern
set maxacc [expr {0.4 * 9.8}]; # max ground acceleration
set accelSeries "Series -dt $DtGround -filePath eq.txt -factor $maxacc";
# tag direction
pattern UniformExcitation 2 1 -accel $accelSeries

# eigen value
rayleigh 0. 0. 0. [expr {2. * $xDamp / pow($omega2,0.5)}]
puts "damping defined"

# create the analysis
wipeAnalysis
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-12 10
algorithm Newton
integrator Newmark $gamma $beta
analysis Transient
puts "Transient analysis defined"

set ok [analyze 3328 $DtAnalysis]

# if the analysis fails at a step we try an initial stiffness iteration for that step
# and then go back to the current stiffness for the rest

if {$ok != 0} {
set tCurrent [getTime]
set ok 0
while {$ok == 0 && $tCurrent < $TmaxGround} {

set ok [analyze 1 $DtAnalysis]

# if the analysis fails try initial tangent iteration
if {$ok != 0} {
puts "regular newton failed... lets try an initial stiffness for this step"
test NormDispIncr 1.0e-12 100 1
algorithm ModifiedNewton -initial
set ok [analyze 1 $DtAnalysis]
if {$ok == 0} {puts "that worked...back to regular newton"}
test NormDispIncr 1.0e-12 10
algorithm Newton
}

set tCurrent [getTime]

}
}

if {$ok == 0} {
puts "Transient analysis completed successfully";
} else {
puts "Transient analysis completed failed";
}

[silvia]

what does the gm file look like?
many files are in inches/sec^2, not g, it depends on the database

[andrewfirst]

The grounmotion file is processed already with PGA=1 and no units. So in the Pattern UniformExicition.. I set the 'factor' to 0.4g.

I've tried using the 'nonlinearbeamcolumn' and removed the uniform load on girders. Then the computation failed at about 0.3s because of 'Can't obtain equilibrium within several elements.

Interesting it is that when I set the beam to elastic section, the computation runs well.

[silvia]

test your section moment curvature analysis

[andrewfirst]

All the sections has been run for moment curveture analysis. No problem.

[silvia]

don't know what to say, then.

[andrewfirst]

Dear Silvia:
I've changed all the element into nonlinearbeamcolumn and defined uniform girder load as many concentrated load.
The TH analysis can obtain reasonable results but can only run for several seconds. Then it display many messages about "can not get compatible forces and deformations for elemtnt...".
I used several numerical methods such as modified Newton, Broyden and so on. But it still couldn't work.
Would you please give me some advice about the set for values of some parameters, such as Tol, maxIternum? Thanks!

[silvia]

your sections have reached really high deformation demands.

[Trirat_to]

Hi, Andrewfirst

I have a problem like you. The OpenSees got a very large displacement result from timehistory analysis (very small EQ) .

How did you solve your problem?