Nonlinear column with harmonic load on the top

[diegobuitrago]

Hi everyone,
I have this nonlinear column with a harmonic load on the top and is working using a linear analysis. However, I am not sure why the nonlinear analysis is not working good. It seems that the system cross the elastic limit but after a couple of cycles come back to the linear range.

# column
# English System (ft, s, lb, lb-sec/ft)

wipe;
model basic -ndm 2 -ndf 3;
set data1Dir Data1; # set up name of data directory
file mkdir $data1Dir;

#INPUT
set HCol 144; # Column height [in]
set bcol 12; # Column section Width [in]
set hcol 12; # Column section Height [in]
set H [expr 2*$HCol]; # Height of the structure [in]
set grav 386.22; # [in/sec2]
set fc -3750; # Streght of concrete (Compression) No confinement [lb/in2]

set nuc 0.2; # Poisson concrete
set GammaConcrete 0.08680556; # [lb/in3] (150 lb/ft3)
set ColMatTagFlex 2; # assign a tag number to the column flexural behavior
set ColMatTagAxial 3; # assign a tag number to the column axial behavior
set ColSecTag 1;

# CALCULUS
set PI [expr 2*asin(1.0)];
set Ec 3417000; # Young's modulus of concrete [lb/in2]
set AgCol [expr $bcol*$hcol]; # Area ST Column[in2]
set IzCol [expr $bcol*pow($hcol,3)/12]; # Moment of Inertia [in4]

# COORDINATES
node 1 0.0 0.0;

node 2 0.0 $HCol;

fix 1 1 1 1;


set Qcol [expr $GammaConcrete*$AgCol]; # Column [lb/in]
set QcolT [expr $Qcol]; # Total load Column [lb/in]
set Wcol [expr $QcolT*$HCol]; # Total load Column [lb]

# MASS
set WeightNode [expr $Wcol]
set MassNode [expr (($Wcol)/$grav)]
mass 2 $MassNode 1e-9 0.0;


# COLUMN section
# calculated stiffness parameters
set EICol [expr $Ec*$IzCol]; # EI, for moment-curvature relationship
set EACol [expr $Ec*$AgCol]; # EA, for axial-force-strain relationship
set MyCol 934000; # yield moment [lb-in]
set PhiYCol 0.4354e-3; # yield curvature (1/in)
set EIColCrack [expr $MyCol/$PhiYCol]; # cracked section inertia (or 0.70*Ig*E)
set b 0.01 ; # strain-hardening ratio (ratio between post-yield tangent and initial elastic tangent)

uniaxialMaterial Steel01 $ColMatTagFlex $MyCol $EIColCrack $b; # bilinear behavior for flexure
uniaxialMaterial Elastic $ColMatTagAxial $EACol; # this is not used as a material, this is an axial-force-strain response
section Aggregator $ColSecTag $ColMatTagAxial P $ColMatTagFlex Mz; # combine axial and flexural behavior into one section (no P-M interaction here)
set ColTransfTag 1; # associate a tag to column transformation
geomTransf Linear $ColTransfTag ;

# element connectivity:
set numIntgrPts 5; # number of integration points for force-based element
element nonlinearBeamColumn 1 1 2 $numIntgrPts $ColSecTag $ColTransfTag; # self-explanatory when using variables





# Define RECORDERS -------------------------------------------------------------
recorder Node -file $data1Dir/Node2disp.out -time -node 2 -dof 1 2 3 disp; # displacements of support nodes
recorder Node -file $data1Dir/Node1react.out -time -node 1 -dof 1 2 3 reaction; # support reaction
recorder Element -file $data1Dir/FCol1.out -time -ele 1 globalForce; # element forces -- column
recorder Element -file $data1Dir/DCol.out -time -ele 1 deformation; # element deformations -- column
recorder EnvelopeElement -file $data1Dir/ENVFCol1.out -time -ele 1 globalForce; # element forces -- column
recorder Element -file $data1Dir/ForceColSec1.out -time -ele 1 section 1 force; # Column section forces, axial and moment, node i
recorder Element -file $data1Dir/DefoColSec1.out -time -ele 1 section 1 deformation; # section deformations, axial and curvature, node i
recorder Element -file $data1Dir/ForceColSec$numIntgrPts.out -time -ele 1 section $numIntgrPts force; # section forces, axial and moment, node j
recorder Element -file $data1Dir/DefoColSec$numIntgrPts.out -time -ele 1 section $numIntgrPts deformation; # section deformations, axial and curvature, node j

recorder Node -file $data1Dir/DFree.out -time -node 2 -dof 1 2 3 disp; # displacements of free nodes
recorder Node -file $data1Dir/DBase.out -time -node 1 -dof 1 2 3 disp; # displacements of support nodes
recorder Node -file $data1Dir/RBase.out -time -node 1 -dof 1 2 3 reaction; # support reaction
recorder Drift -file $data1Dir/Drift.out -time -iNode 1 -jNode 2 -dof 1 -perpDirn 2 ; # lateral drift
recorder Element -file $data1Dir/FCol.out -time -ele 2 globalForce; # element forces -- column
recorder Element -file $data1Dir/ForceColSec1.out -time -ele 1 section 1 force; # Column section forces, axial and moment, node i
recorder Element -file $data1Dir/DefoColSec1.out -time -ele 1 section 1 deformation; # section deformations, axial and curvature, node i
recorder Element -file $data1Dir/ForceColSec$numIntgrPts.out -time -ele 1 section $numIntgrPts force; # section forces, axial and moment, node j
recorder Element -file $data1Dir/DefoColSec$numIntgrPts.out -time -ele 1 section $numIntgrPts deformation; # section deformations, axial and curvature, node j


# define GRAVITY -------------------------------------------------------------
pattern Plain 1 Linear {
load 2 0. -[expr $Wcol] 0.; # node#, FX FY MZ -- superstructure-weight
}
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
set Tol 1.0e-8; # determine if convergence has been achieved at the end of an iteration step
test NormDispIncr $Tol 6 ; # determine if convergence has been achieved at the end of an iteration step
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
set NpGrav 10; # aplico la gravedad en 10 pasos
set DGrav [expr 1./$NpGrav]; # incrementos para la aplicación de carga
integrator LoadControl $DGrav; # determine the next time step for an analysis, # apply gravity in 10 steps
analysis Static; # define type of analysis static or transient
analyze $NpGrav; # Indicate the number of steps to carry out the analysis
loadConst -time 0.0; # hold gravity constant and restart time


# define HARMONIC FORCE
set T 0.18; # T= Period of the harmonic force
set Fy 30000. ; # magnitude
set SineSeries "Trig 0. 15. $T";

pattern Plain 2 $SineSeries {
load 2 $Fy 0.0 0.0;
}

# apply Rayleigh DAMPING from $xDamp
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial

set xDamp 0.05; # 5% damping ratio
set lambda [eigen 1]; # eigenvalue mode 1
set omega [expr pow($lambda,0.5)];
set alphaM 0.; # M-prop. damping; D = alphaM*M
set betaKcurr 0.; # K-proportional damping; +beatKcurr*KCurrent
set betaKcomm [expr 2.*$xDamp/($omega)]; # K-prop. damping parameter; +betaKcomm*KlastCommitt
set betaKinit 0.; # initial-stiffness proportional damping +beatKinit*Kini
# define damping
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping


# create the analysis
wipeAnalysis; # clear previously-define analysis parameters
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
algorithm Linear # use Linear algorithm for linear analysis
integrator Newmark 0.5 0.25 ; # determine the next time step for an analysis
analysis Transient; # define type of analysis: time-dependent
analyze 15000 0.001; # apply 15 seconds step 0.001 15000*0.001= 15 in analysis

puts "Done!"