the main program, which then calls a script posted at the bottom:
# -------------------------------------------------------------------------------------------------
# exampoFrame2Dpush.tcl: pushover analysis of frame in 2D, nonlinear beam-column elements
# by Silvia Mazzoni, 2005
#
# 7________(9)______8________(10)___________9 _
# | | | |
# | | | |
# | | | |
# |(4) |(5) |(6) LcolTop
# | | | |
# | | | |
# | | | |
# 4________(7)______5________(8)____________6 _
# | | | |
# | | | |
# | | | |
# | | | |
# | | | |
# |(1) |(2) |(3) LcolBot
# | | | |
# | | | |
# | | | |
# | | | |
# === 1 === 2 === 3 -
#
# |----LbeamLeft----|------LbeamRight-------|
#
#
#
# SET UP ----------------------------------------------------------------------------
set dataDir data/Frame2D; # this variable will be used with the recorders
file mkdir $dataDir; # create data directory
# UNITS-- define system of units used in the tcl script
set in 1.; # define basic units
set sec 1.; # define basic units
set kip 1.; # define basic units
set ft [expr 12.*$in]; # define engineering units
set ksi [expr $kip/pow($in,2)];
set psi [expr $ksi/1000.];
set in2 [expr $in*$in]; # inch^2
set in4 [expr $in*$in*$in*$in]; # inch^4
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
set cm [expr $in/2.54]; # SI centimeter unit
# set up OpenSees model ----------------------------------------------------------------------------
wipe; # clear memory of all past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder
# Define MODEL -----------------------------------------------------------------------
# geometry
set LcolBot [expr 25*$ft]; # bottom-column length
set LcolTop [expr 20*$ft]; # top-column length
set LbeamLeft [expr 25*$ft]; # left-beam length
set LbeamRight [expr 30*$ft]; # right-beam length
set wDL [expr 75.*$kip/$ft]; # distributed beam dead-load
# calculated parameters
set PnodeL [expr -$wDL*$LbeamLeft/2]; # force at left column
set PnodeM [expr -$wDL*($LbeamLeft+$LbeamRight)/2]; # force at middle columns
set PnodeR [expr -$wDL*($LbeamRight)/2]; # force at right columns
set MnodeL [expr -$wDL*pow($LbeamLeft,2)/12]; # moment at left column
set MnodeM [expr $wDL*(pow($LbeamLeft,2)-pow($LbeamRight,2))/12]; # moment at middle columns
set MnodeR [expr +$wDL*($LbeamRight)/2]; # moment at right columns
set MassL [expr -$PnodeL/$g]; # nodal mass at left columns
set MassM [expr -$PnodeM/$g]; # nodal mass at middle columns
set MassR [expr -$PnodeR/$g]; # nodal mass at right columns
set Weight [expr 3.*$wDL*($LbeamLeft+$LbeamRight)]; # total structural weight
set Mass [expr $Weight/$g]; # total superstructural mass
set Lcol [expr ($LcolBot+$LcolTop)/2]; # average column height for pushover drift calculation
# Define nodes
node 1 0 0
node 2 $LbeamLeft 0
node 3 [expr $LbeamLeft+$LbeamRight] 0
node 4 0 $LcolBot -mass $MassL 1e-9 0.
node 5 $LbeamLeft $LcolBot -mass $MassM 1e-9 0.
node 6 [expr $LbeamLeft+$LbeamRight] $LcolBot -mass $MassR 1e-9 0.
node 7 0 [expr $LcolBot+$LcolTop] -mass $MassL 1e-9 0.
node 8 $LbeamLeft [expr $LcolBot+$LcolTop] -mass $MassM 1e-9 0.
node 9 [expr $LbeamLeft+$LbeamRight] [expr $LcolBot+$LcolTop] -mass $MassR 1e-9 0.
set IDctrlNode 8; # ID for node for any lateral load/disp-ctrl analysis and for drift recorder
set iIDbaseNode "1 2 3"; # ID for support node for imposed support ground motion.
set IDdof 1; # DOF number for any lateral load/disp-ctrl analysis
set IDdofPerp 2; # perpendicural dof for drift recorder
set iIDpushNodeLevel1 "4 5 6"; # ID for the nodes where lateral loads are applied in reference load
set iIDpushNodeLevel2 "7 8 9"; # ID for the nodes where lateral loads are applied in reference load
# Single point constraints -- Boundary Conditions
fix 1 1 1 1
fix 2 1 1 1
fix 3 1 1 1
# MATERIAL properties -------------------------------------------------------------
# elastic-Material properties
set fc [expr -5.5*$ksi]; # CONCRETE Compressive Strength (+Tension, -Compression)
set Ec [expr 57*$ksi*sqrt(-$fc/$psi)]; # Concrete Elastic Modulus
set Es [expr 29000.*$ksi]; # modulus of steel
# confined concrete
set Ec [expr 57*$ksi*sqrt(-$fc/$psi)]; # 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 $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.003; # strain at maximum stress
set fc2U [expr 0.1*$fc]; # ultimate stress
set eps2U -0.006; # strain at ultimate stress
# concrete02 material properties:
set lambda 0.1 ; # ratio between unloading slope at $epscu and initial slope
set ftC [expr -$fc1C/10.]; # tensile strength +tension
set ftU [expr -$fc1U/10.]; # tensile strength +tension
set Ets [expr $Ec/10.]; # tension softening stiffness
# reinforcing steel
set Fy [expr 70.*$ksi]; # STEEL yield stress
set Es [expr 29000.*$ksi]; # modulus of steel
set epsY [expr $Fy/$Es]; # steel yield strain
set Fy1 [expr 95.*$ksi]; # steel stress post-yield
set epsY1 0.03; # steel strain post-yield
set Fu [expr 112.*$ksi]; # ultimate stress of steel
set epsU 0.08; # ultimate strain of steel
set Bs [expr ($Fy1-$Fy)/($epsY1-$epsY)/$Es]; # post-yield stiffness ratio of steel
set R0 10; # Steel02 -- control the transition from elastic to plastic branches
set cR1 0.925; # Steel02 -- control the transition from elastic to plastic branches
set cR2 0.15; # Steel02 -- control the transition from elastic to plastic branches
set a2 0.1; # Steel02 -- isotropic hardening parameter, associated with a1
set a1 [expr $a2*($Fy/$Es)]; # Steel02 -- isotropic hardening parameter, increase of compression yield envelope as proportion of yield strength after a plastic strain of $a2*($Fy/E0).
set a4 0.1; # Steel02 -- isotropic hardening parameter, associated with a3
set a3 [expr $a4*($Fy/$Es)]; # Steel02 -- isotropic hardening parameter, increase of tension yield envelope as proportion of yield strength after a plastic strain of $a4*($Fy/E0)
# associate tags to library of materials
set IDconcCore [set icount 1]; # material ID tag
set IDconcCover [incr icount 1]
set IDreinfInel [incr icount 1]
# set up library of materials
uniaxialMaterial Concrete02 $IDconcCore $fc1C $eps1C $fc2C $eps2C $lambda $ftC $Ets; # Core concrete (confined)
uniaxialMaterial Concrete02 $IDconcCover $fc1U $eps1U $fc2U $eps2U $lambda $ftU $Ets; # Cover concrete (unconfined)
uniaxialMaterial Steel02 $IDreinfInel $Fy $Es $Bs $R0 $cR1 $cR2 $a1 $a2 $a3 $a4
# SECTION properties -------------------------------------------------------------
source LibRCrectSection.tcl; # source-in a file with RC rectangular-section definition
# define COLUMN section
set IDcolSecInelBot [set icount 1]; # assign tag to column section -- bottom column
set IDcolSecInelTop [incr icount 1]; # assign tag to column section -- top column
# define BEAM sections
set IDbeamSecInel [incr icount 1]; # assign tag to beam section
# cover spacing is the same for all elements:
set coverH [expr 5.*$in]; # rectangular-RC-column cover - perpendicular to axis of bending (along H)
set coverB [expr 4.*$in]; # rectangular-RC-column side cover -- parallel to axis of bending (along B)
# define COLUMN section geometry -- bottom column
set HcolBot [expr 5.5*$ft]; # rectangular-column Depth -- parallel to local-y axis
set BcolBot [expr 5.5*$ft]; # rectangular-column Width -- parallel to local-z axis (local bending axis)
set numBarsBot 50; # number of longitudinal-reinforcement bars -- same numBars top and bottom
set barAreaBot [expr 1.0*$in*$in]; # longitudinal-reinforcement bar area
# calculated geometry parameters
set AsColBot [expr $numBarsBot*$barAreaBot]; # longitudinal-steel area
set IzCol0 [expr 1./12*$BcolBot*pow($HcolBot,3)]; # about-local-z Rect-column gross moment of inertial
set IyColBot [expr 1./12*$HcolBot*pow($BcolBot,3)]; # about-local-y Rect-column gross moment of inertial
set IzColBot [expr $IzCol0+2*$AsColBot*($HcolBot/2-$coverH)*($HcolBot/2-$coverH)]; #
set AgColBot [expr $HcolBot*$BcolBot]; # rectuangular-column cross-sectional area
# define COLUMN section geometry -- top column
set HcolTop [expr 5.*$ft]; # rectangular-column Depth -- parallel to local-y axis
set BcolTop [expr 5.*$ft]; # rectangular-column Width -- parallel to local-z axis (local bending axis)
set numBarsTop 50; # number of longitudinal-reinforcement bars -- same numBars top and bot
set barAreaTop [expr 1.0*$in*$in]; # longitudinal-reinforcement bar area
# calculated geometry parameters
set AsColTop [expr $numBarsTop*$barAreaTop]; # longitudinal-steel area
set IzCol0 [expr 1./12*$BcolTop*pow($HcolTop,3)]; # about-local-z Rect-column gross moment of inertial
set IyColTop [expr 1./12*$HcolTop*pow($BcolTop,3)]; # about-local-y Rect-column gross moment of inertial
set IzColTop [expr $IzCol0+2*$AsColTop*($HcolTop/2-$coverH)*($HcolTop/2-$coverH)]; #
set AgColTop [expr $HcolTop*$BcolTop]; # rectuangular-column cross-sectional area
# define BEAM section geometry
set Hbeam [expr 6*$ft]; # rectangular-beam Depth -- parallel to local-y axis
set Bbeam [expr 4*$ft]; # rectangular-beam Width -- parallel to local-z axis (local bending axis)
set numBars 50; # number of longitudinal-reinforcement bars -- same numBars top and bottom
set barArea [expr 1.0*$in*$in]; # longitudinal-reinforcement bar area
set AsBeam [expr $numBars*$barArea]; # longitudinal-steel area
# calculated geometry parameters
set IzBeam0 [expr 1./12*$Bbeam*pow($Hbeam,3)]; # about-local-z Rect-beam gross moment of inertial
set IyBeam [expr 1./12*$Hbeam*pow($Bbeam,3)]; # about-local-y Rect-beam gross moment of inertial
set AgBeam [expr $Hbeam*$Bbeam]; # rectuangular-beam cross-sectional area
set IzBeam [expr $IzBeam0+2*$AsBeam*($Hbeam/2-$coverH)*($Hbeam/2-$coverH)]; #
set nfCoreY 16; # number of fibers in the core patch in the y direction
set nfCoreZ 2; # number of fibers in the core patch in the z direction
set nfCoverY 16; # number of fibers in the cover patches with long sides in the y direction
set nfCoverZ 2; # number of fibers in the cover patches with long sides in the z direction
# define fiber section
set nfCoreY 16; # number of fibers in the core patch in the y direction
set nfCoreZ 2; # number of fibers in the core patch in the z direction
set nfCoverY 16; # number of fibers in the cover patches with long sides in the y direction
set nfCoverZ 2; # number of fibers in the cover patches with long sides in the z direction
RCrectSection $IDcolSecInelBot $HcolBot $BcolBot $coverH $coverB $IDconcCore $IDconcCover $IDreinfInel $numBarsBot $barAreaBot $nfCoreY $nfCoreZ $nfCoverY $nfCoverZ
RCrectSection $IDcolSecInelTop $HcolTop $BcolTop $coverH $coverB $IDconcCore $IDconcCover $IDreinfInel $numBarsTop $barAreaTop $nfCoreY $nfCoreZ $nfCoverY $nfCoverZ
RCrectSection $IDbeamSecInel $Hbeam $Bbeam $coverH $coverB $IDconcCore $IDconcCover $IDreinfInel $numBars $barArea $nfCoreY $nfCoreZ $nfCoverY $nfCoverZ
# Define ELEMENTS -------------------------------------------------------------
set np 5; # number of integration points for nonlinearBeamColumn element
set IDcolTrans 1; # associate a tag to column transformation
set IDbeamTrans 2; # associate a tag to column transformation
geomTransf PDelta $IDcolTrans ; # P-delta transformation: second-order effects
geomTransf Linear $IDbeamTrans ; # Linear transformation: no second-order effects
element nonlinearBeamColumn 1 1 4 $np $IDcolSecInelBot $IDcolTrans
element nonlinearBeamColumn 2 2 5 $np $IDcolSecInelBot $IDcolTrans
element nonlinearBeamColumn 3 3 6 $np $IDcolSecInelBot $IDcolTrans
element nonlinearBeamColumn 4 4 7 $np $IDcolSecInelTop $IDcolTrans
element nonlinearBeamColumn 5 5 8 $np $IDcolSecInelTop $IDcolTrans
element nonlinearBeamColumn 6 6 9 $np $IDcolSecInelTop $IDcolTrans
element nonlinearBeamColumn 7 4 5 $np $IDbeamSecInel $IDbeamTrans
element nonlinearBeamColumn 8 5 6 $np $IDbeamSecInel $IDbeamTrans
element nonlinearBeamColumn 9 7 8 $np $IDbeamSecInel $IDbeamTrans
element nonlinearBeamColumn 19 8 9 $np $IDbeamSecInel $IDbeamTrans
# ------------------------------------------------------------- recorders
recorder Node -file $dataDir/DNode.out -time -node 4 5 6 7 8 9 -dof 1 2 3 disp
recorder Node -file $dataDir/DBase.out -time -node 1 2 3 -dof 1 2 3 disp
recorder Node -file $dataDir/RBase.out -time -node 1 2 3 -dof 1 2 3 reaction
recorder Drift -file data/DrNodeLevel1.out -time -iNode 5 -jNode 2 -dof $IDdof -perpDirn $IDdofPerp
recorder Drift -file data/DrNodeLevel2.out -time -iNode 8 -jNode 5 -dof $IDdof -perpDirn $IDdofPerp
recorder Element -file $dataDir/DelSec1.out -time -ele 1 2 3 section 1 deformation
recorder Element -file $dataDir/FelSec1.out -time -ele 1 2 3 section 1 force
# --- ID tags for Loads ----------
set IDgravity 1
set IDpush 2
# -------------------------------------------------define and apply gravity load
pattern Plain $IDgravity Linear {
load 4 0.0 $PnodeL $MnodeL
load 5 0.0 $PnodeM $MnodeM
load 6 0.0 $PnodeR $MnodeR
load 7 0.0 $PnodeL $MnodeL
load 8 0.0 $PnodeM $MnodeM
load 9 0.0 $PnodeR $MnodeR
}
# Apply GRAVITY -------------------------------------------------
# perform eigenvalue analysis to determine fundamental period and initial stiffness
set lambda [eigen 1]; # perform eigenvalue analysis to determine fundamental mode
set omega [expr pow($lambda,0.5)]; # calculate fundamental frequency
set K [expr $lambda*$Mass]; # initial stiffness
puts "K=$K (kip/in)"; # this output assumes basic units: kip, inch
set Tperiod [expr 2*$PI/$omega]; # calculate fundamental period
constraints Transformation; # how it handles boundary conditions, enforce constraints through the transformation
numberer Plain; # renumber dof's to minimize band-width (optimization), if you want to
system BandGeneral; # construct the LinearSOE and LinearSolver objects to store and solve the system of equations in the analysis
set tol 1.e-6; # convergence tolerance
set maxNumIter 6; # maximum number of iterations that will be performed before "failure to converge" is returned
test NormDispIncr $tol $maxNumIter; #tests positive force convergence if the 2-norm of the x vector (the displacement increment) in the LinearSOE object is less than the specified tolerance
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
set dLambda1 [expr 0.1*$in]; # first load increment (pseudo-time step) in the next invocation of the analysis command
integrator LoadControl $dLambda1 ; # load-control integrator
analysis Static; # define what type of analysis is to be performed
set numIncr 10; # number of load steps
analyze $numIncr; # perform analysis based on previously-defined analysis objects
loadConst -time 0.0; # set pseudo time in the domain to zero
# define and perform static PUSHOVER ANALYSIS ------------------------------------------------------------------
set DxPush [expr 0.1*$in]; # Displacement increment for pushover analysis
set DmaxPush [expr 0.015*$Lcol]; # this is now defined with model!!! # maximum displamcement for pushover analysis -- should be like: 0.05*$Lcol
set Nsteps [expr int($DmaxPush/$DxPush)];
# Reference lateral load for pushover analysis
set HloadLevel1 [expr $Weight/(3*2+3)];
set HloadLevel2 [expr 2*$Weight/(3*2+3)];
# create load pattern for lateral loads
pattern Plain $IDpush Linear {
foreach IDpushNodeLevel1 $iIDpushNodeLevel1 IDpushNodeLevel2 $iIDpushNodeLevel2 {
load $IDpushNodeLevel1 $HloadLevel1 0.0 0.0
load $IDpushNodeLevel2 $HloadLevel2 0.0 0.0
}
}
# since this is a static analysis, just as gravity was, we don't need to re-define all analysis objects.
integrator DisplacementControl $IDctrlNode $IDdof $DxPush
analyze $Nsteps
puts donePUSHOVER!!!
# -----------------------------------------------------------------------
LibRCrectSection.tcl:
#---------------------------------------------------------------------
# RCrectSection.tcl:
# Define a procedure which generates a rectangular reinforced concrete section
# with one layer of steel at the top & bottom, along the local-y axis and a
# confined core.
# by: Silvia Mazzoni, 2005
# adapted from Michael H. Scott, 2003
#
# Formal arguments
# id - tag for the section that is generated by this procedure
# Hsec - depth of section, along local-y axis
# Bsec - width of section, along local-z axis
# cH - distance from section boundary to neutral axis of reinforcement
# cB - distance from section boundary to side of reinforcement
# coreID - material tag for the core patch
# coverID - material tag for the cover patches
# steelID - material tag for the reinforcing steel
# numBars - number of reinforcing bars around the section perimeter
# barArea - cross-sectional area of each reinforcing bar
# nfCoreY - number of fibers in the core patch in the y direction
# nfCoreZ - number of fibers in the core patch in the z direction
# nfCoverY - number of fibers in the cover patches with long sides in the y direction
# nfCoverZ - number of fibers in the cover patches with long sides in the z direction
#
# ___________ __ __
# | | | |cH
# | xxxxxxx | | -
# | | |
# | ^y | |
# | | | Hsec
# | z<---- | |
# | | |
# | | |
# | xxxxxxx | |
# |___________| _|_
#
# |----Bsec---|
#
# |--|cB
#
#
# y
# ^
# |
# ---------------------
# |\ cover /|
# | \---------------/ |
# |c| |c|
# |o| |o|
# z <--------|v| core |v| Hsec
# |e| |e|
# |r| |r|
# | /---------------\ |
# |/ cover \|
# ---------------------
# Bsec
#
#
#
# Notes
# The core concrete ends at the NA of the reinforcement
# The reinforcing bars are all the same size
# The center of the section is at (0,0) in the local axis system
#
proc RCrectSection {id Hsec Bsec cH cB coreID coverID steelID numBars barArea nfCoreY nfCoreZ nfCoverY nfCoverZ} {
# The distance from the section z-axis to the edge of the cover concrete
# in the positive y direction -- outer edge of cover concrete
set coverY [expr $Hsec/2.0]
# The distance from the section y-axis to the edge of the cover concrete
# in the positive z direction -- outer edge of cover concrete
set coverZ [expr $Bsec/2.0]
# Determine the corresponding values from the respective axes to the
# edge of the core concrete/inner edge of cover concrete
set coreY [expr $coverY-$cH]
set coreZ [expr $coverZ-$cB]
# Define the fiber section
section fiberSec $id {
# Define the core patch
patch quadr $coreID $nfCoreZ $nfCoreY -$coreY $coreZ -$coreY -$coreZ $coreY -$coreZ $coreY $coreZ
# Define the four cover patches
patch quadr $coverID 1 $nfCoverY -$coverY $coverZ -$coreY $coreZ $coreY $coreZ $coverY $coverZ
patch quadr $coverID 1 $nfCoverY -$coreY -$coreZ -$coverY -$coverZ $coverY -$coverZ $coreY -$coreZ
patch quadr $coverID $nfCoverZ 1 -$coverY $coverZ -$coverY -$coverZ -$coreY -$coreZ -$coreY $coreZ
patch quadr $coverID $nfCoverZ 1 $coreY $coreZ $coreY -$coreZ $coverY -$coverZ $coverY $coverZ
# define reinforcing layers along constant values of y (in the z direction)
layer straight $steelID $numBars $barArea -$coreY $coreZ -$coreY -$coreZ
layer straight $steelID $numBars $barArea $coreY $coreZ $coreY -$coreZ
}; # end of fibersection definition
}; # end of procedure