Force based model for Cyclic Pushover

[arpitak]

Hello,
I am a new Opensees user and I am analyzing an RC Cantilever column using cyclic pushover but I want to do this using force based formulation. I am using the loadcontrol integrator for this and I get a convergence error. Can somebody please check my code and suggest ways to solve this issue.

Thank you!
Arpita

# --------------------------------------------------------------------------------------------------
# Example 3. 2D Cantilever -- Build Model
# nonlinearBeamColumn element, uniaxial inelastic section
# Silvia Mazzoni & Frank McKenna, 2006
#
# ^Y
# |
# 2 __
# | |
# | |
# | |
# (1) LCol
# | |
# | |
# | |
# =1= _|_ -------->X
#

# SET UP ----------------------------------------------------------------------------
wipe; # clear memory of all past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm=#dimension, ndf=#dofs
set dataDir Cyclic4D_Force; # set up name for data directory
file mkdir $dataDir/; # create data directory
#set GMdir "../GMfiles"; # ground-motion file directory

# define UNITS ----------------------------------------------------------------------------
set in 1.; # define basic units -- output units
set kip 1.; # define basic units -- output units
set sec 1.; # define basic units -- output units
set LunitTXT "inch"; # define basic-unit text for output
set FunitTXT "kip"; # define basic-unit text for output
set TunitTXT "sec"; # define basic-unit text for output
set ft [expr 12.*$in]; # define engineering units
set ksi [expr $kip/pow($in,2)];
set psi [expr $ksi/1000.];
set lbf [expr $psi*$in*$in]; # pounds force
set pcf [expr $lbf/pow($ft,3)]; # pounds per cubic foot
set in2 [expr $in*$in]; # inch^2
set in4 [expr $in*$in*$in*$in]; # inch^4
set cm [expr $in/2.54]; # centimeter, needed for displacement input in MultipleSupport excitation
set PI [expr 2*asin(1.0)]; # define constants
set g [expr 32.2*$ft/pow($sec,2)]; # gravitational acceleration
set Ubig 1.e10; # a really large number
set Usmall [expr 1/$Ubig]; # a really small number

# define GEOMETRY -------------------------------------------------------------
set Weight [expr "1130.4*$kip"]; # superstructure weight
# define section geometry
set DCol [expr 5.*$ft]; # Column Diameter
#set HCol [expr 5.*$ft]; # Column Depth
#set BCol [expr 5.*$ft]; # Column Width
set LCol [expr 4*$DCol]; # column length
puts $LCol
puts $Weight
# calculated parameters
set PCol [expr $Weight]; # nodal dead-load weight per column
set Mass [expr $PCol/$g]; # nodal mass
# calculated geometry parameters
set ACol [expr 3.14*pow($DCol,2)/4]; # cross-sectional area circular column
#set ACol [expr $BCol*$HCol]; # cross-sectional area
set IzCol [expr 1./4.*3.14*pow($DCol/2,4)]; # Column moment of inertia circular column
#set IzCol [expr 1./12.*$BCol*pow($HCol,3)]; # Column moment of inertia

# nodal coordinates:
node 1 0 0; # node#, X, Y
node 2 0 $LCol

# Single point constraints -- Boundary Conditions
fix 1 1 1 1; # node DX DY RZ

# we need to set up parameters that are particular to the model.
set IDctrlNode 2; # node where displacement is read for displacement control
set IDctrlDOF 1; # degree of freedom of displacement read for displacement control
set iSupportNode "1"; # define support node, if needed.

# nodal masses:
mass 2 $Mass 1e-9 0.; # node#, Mx My Mz, Mass=Weight/g, neglect rotational inertia at nodes

# Define ELEMENTS & SECTIONS -------------------------------------------------------------
set ColSecTag 1; # assign a tag number to the column section
# define section geometry
set coverCol [expr 2.*$in]; # Column cover to reinforcing steel NA.
set numBarsCol 20; # number of longitudinal-reinforcement bars in column. (symmetric top & bot)
set barAreaCol [expr 2.826*$in2]; # area of longitudinal-reinforcement bars


# MATERIAL parameters -------------------------------------------------------------------
set IDconcU 1; # material ID tag -- unconfined cover concrete
set IDconcC 2; # material ID tag -- unconfined cover concrete
set IDreinf 3; # material ID tag -- reinforcement
# nominal concrete compressive strength
set fc [expr -4.0*$ksi]; # CONCRETE Compressive Strength, ksi (+Tension, -Compression)
set Ec [expr 57*$ksi*sqrt(-$fc/$psi)]; # Concrete Elastic Modulus
# unconfined concrete
set fc1U $fc; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U -0.003; # strain at maximum strength of unconfined concrete
set fc2U [expr 0.2*$fc1U]; # ultimate stress
set eps2U -0.01; # strain at ultimate stress
set lambda 0.1; # ratio between unloading slope at $eps2 and initial slope $Ec
#confined concrete Concrete04
set ROs 0.00775; #ratio of the volume of transverse confining steel to the volume of confined concrete core = 4Asp/ds = 0.00775 (as given in the book)
set ke 0.95; #typical value for circular section
set fyh [expr 66.*$ksi]; # yield strength of transverse reinforcement
set fl [expr -0.5*$ke*$ROs*$fyh]; #Effective lateral confining stress
set fc1C [expr $fc*(2.254*sqrt(1+7.94*$fl/$fc)-2*$fl/$fc-1.254)]; #compressive stress of confined concrete
#set Ecc [expr 57*sqrt(-$fc1C*1000)]; # floating point values defining initial stiffness Ec = 57000*sqrt|fcc|in psi
set epssu 0.1; #steel strain at max tensile stress-Park Reinforcing Model- Complex Strain Hardening
set eps1C [expr $eps1U*(1+5*($fc1C/$fc-1))]; #Strain at peak compressive stress
#set eps2C -0.02; #Strain at peak compressive stress
set eps2C [expr (-0.004+(1.4*$ROs*$fyh*$epssu/$fc1C))]; #Ultimate compressive strain the confined core-Priestly equation
set X [expr $eps1U/$eps1C]; #x=epsc/eps1C
set Esec [expr $fc1C/$eps1C];
set r [expr $Ec/($Ec-$Esec)];
set fc2C [expr $fc1C*$X*$r/($r-1+$X)]; # ultimate stress of the confined core
#set beta 0.1; #floating point value defining the exponential curve parameter to define the residual stress (as a factor of $ft) at $etu
# tensile-strength properties
set ftU [expr -0.14*$fc1U]; # tensile strength +tension
set Ets [expr $ftU/0.002]; # tension softening stiffness
# -----------
set Fy [expr 66.*$ksi]; # STEEL yield stress
set Es [expr 29000.*$ksi]; # modulus of steel
set Bs 0.01; # strain-hardening ratio
set R0 18; # control the transition from elastic to plastic branches
set cR1 0.925; # control the transition from elastic to plastic branches
set cR2 0.15; # control the transition from elastic to plastic branches
uniaxialMaterial Concrete02 $IDconcU $fc1U $eps1U $fc2U $eps2U $lambda $ftU $Ets; # build cover concrete (unconfined)
uniaxialMaterial Concrete04 $IDconcC $fc1C $eps1C $eps2C $Ec; #<$ftC $epssu>; #<$beta>; # build core concrete (confined)
uniaxialMaterial Steel02 $IDreinf $Fy $Es $Bs $R0 $cR1 $cR2; # build reinforcement material

set ri 0.0; # inner radius of the section, only for hollow sections
set ro [expr $DCol/2]; # overall (outer) radius of the section
set nfCoreR 20; # number of radial divisions in the core (number of "rings")
set nfCoreT 40; # number of theta divisions in the core (number of "wedges")
set nfCoverR 2; # number of radial divisions in the cover
set nfCoverT 40; # number of theta divisions in the cover

# Define the fiber section
section fiberSec $ColSecTag {
set rc [expr $ro-$coverCol]; # Core radius
patch circ $IDconcC $nfCoreT $nfCoreR 0 0 $ri $rc 0 360; # Define the core patch
patch circ $IDconcU $nfCoverT $nfCoverR 0 0 $rc $ro 0 360; # Define the cover patch
set theta [expr 360.0/$numBarsCol]; # Determine angle increment between bars
layer circ $IDreinf $numBarsCol $barAreaCol 0 0 $rc $theta 360; # Define the reinforcing layer
}; # end of fibersection definition

# define geometric transformation: performs a linear geometric transformation of beam stiffness and resisting force from the basic system to the global-coordinate system
set ColTransfTag 1; # associate a tag to column transformation
set ColTransfType Linear ; # options, Linear PDelta Corotational
geomTransf $ColTransfType $ColTransfTag ;


# element connectivity:
set numIntgrPts 2; # 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 $dataDir/DFree.out -time -node 2 -dof 1 2 3 disp; # displacements of free nodes
recorder Node -file $dataDir/DBase.out -time -node 1 -dof 1 2 3 disp; # displacements of support nodes
recorder Node -file $dataDir/RBase.out -time -node 1 -dof 1 2 3 reaction; # support reaction
recorder Drift -file $dataDir/Drift.out -time -iNode 1 -jNode 2 -dof 1 -perpDirn 2 ; # lateral drift
recorder Element -file $dataDir/FCol.out -time -ele 2 globalForce; # element forces -- column
recorder Element -file $dataDir/ForceColSec1.out -time -ele 1 section 1 force; # Column section forces, axial and moment, node i
recorder Element -file $dataDir/DefoColSec1.out -time -ele 1 section 1 deformation; # section deformations, axial and curvature, node i
recorder Element -file $dataDir/ForceColSec$numIntgrPts.out -time -ele 1 section $numIntgrPts force; # section forces, axial and moment, node j
recorder Element -file $dataDir/DefoColSec$numIntgrPts.out -time -ele 1 section 1 deformation; # section deformations, axial and curvature, node j
recorder Element -xml $dataDir/PlasticRotation.out -time -ele 1 plasticRotation; # section deformations, axial and curvature, node j

# define GRAVITY -------------------------------------------------------------
pattern Plain 1 Linear {
load 2 0 -$PCol 0
}


# create analysis & perform analysis
constraints Transformation
numberer RCM
system BandGeneral
test NormDispIncr 1.0e-6 6 2
algorithm Newton
integrator LoadControl 0.1
analysis Static
analyze 10
# look at what happened to node 2
print node 2

# set gravity loads constant & reset time in domain
loadConst -time 0.0

# create load pattern for lateral pushover load
set Hload [expr $Weight]; # define the lateral load as a proportion of the weight so that the pseudo time equals the lateral-load coefficient when using linear load pattern
set iPushNode "2"; # define nodes where lateral load is applied in static lateral analysis
pattern Plain 200 Linear {; # define load pattern -- generalized
foreach PushNode $iPushNode {
load $PushNode $Hload 0.0 0.0
}
}

# do some cyclic analysis
foreach {numSteps stepSize} {10 0.1 10 -0.1 10 -0.1 10 0.1 10 0.1} {
integrator LoadControl $stepSize
analyze $numSteps
set time [getTime]
set disp [nodeDisp 2 1]
puts “Time: $time Displacement $disp”
}

[selimgunay]

Generally loadcontrol would not work after yield unless you use very small increments to increase the force and if there is a negative stiffness after yield, it will not work for sure as the forces need to be increased while they would always increase with loadcontrol

[arpitak]

Hi Selimgunay,

Thank you for your answer. So can I just use nonlinear beamColumn element to specify force based formulation and keep the integrator as displacement control in the solver as used in example 4-Portal frame cyclic pushover?

Thanks

[selimgunay]

Yes that is correct.