help: time history analysis with concentrated plasticity

Forum for OpenSees users to post questions, comments, etc. on the use of the OpenSees interpreter, OpenSees.exe

Moderators: silvia, selimgunay, Moderators

Post Reply
setareh20
Posts: 26
Joined: Wed Mar 24, 2010 10:12 pm
Location: IIEES
Contact:

help: time history analysis with concentrated plasticity

Post by setareh20 »

Hi,
Thank you so much for your cooperation.
I want to do time history analysis with concentrated plasticity, but when I ran it for different records, it had the same error message for all of them and it did not work completely.
Please check my script to do time history analysis and guide me to run it completely without any error.
Sincerely yours
# --------------------------------------------------------------------------------------------------
# Example: 1-Story 1-bay 3D RC Moment Frame with Concentrated Plasticity- For Test Only
# Centerline Model with Concentrated Plastic Hinges at Beam-Column Joint (After FEMA-P695 Deterioration Model)
# Units: KN, m, sec

###################################################################################################
# Set Up & Source Definition
###################################################################################################
#wipe all; # clear memory of past model definitions
model BasicBuilder -ndm 3 -ndf 6; # Define the model builder, ndm = #dimension, ndf = #dofs
source DisplayModel3D.tcl; # procedure for displaying a 3D perspective of model
source DisplayPlane.tcl; # procedure for displaying a plane in a model
# source rotSpring3DbeamX.tcl; # procedure for defining a rotational spring (zero-length element)- X beams
# source rotSpring3DbeamZ.tcl; # procedure for defining a rotational spring (zero-length element)- Z beams
# source rotSpring3DCol.tcl; # procedure for defining a rotational spring (zero-length element)-
# #source rotLeaningCol.tcl; # procedure for defining a rotational spring (zero-length element) with very small stiffness

###################################################################################################
# Define Analysis Type
###################################################################################################
# Define type of analysis: "pushover" = pushover; "dynamic" = dynamic
set analysisType "dynamic";

if {$analysisType == "pushover"} {
set dataDir 1st-3d-Pushover-Output; # name of output folder
file mkdir $dataDir; # create output folder
}
if {$analysisType == "dynamic"} {
set dataDir 1st-3d-Dynamic-Output; # name of output folder
file mkdir $dataDir; # create output folder
}

###################################################################################################
# Define Building Geometry, Nodes, and Constraints ##################################################################################################
# define structure-geometry parameters
set NStories 1; # number of stories
set NBays 1; # number of frame bays (excludes bay for P-delta column)
set XBay 6.0; # bay width in meters- X direction
set ZBay 6.0; # bay width in meters-Z direction
set HStory1 3.5; # 1st story height in meters
set HStoryTyp 3.5; # story height of other stories in meters (same as story 1)
set HBuilding [expr $HStory1 + ($NStories-1)*$HStoryTyp]; # height of building
#set g 9.81;
# calculate locations of beam/column joints:
# set Floor1 0.0; # ground floor
# set Floor2 [expr $Floor1 + $HStory1];
# calculate joint offset distance for beam plastic hinges
# set phlat23 [expr 0.0]; # lateral dist from beam-col joint to loc of hinge on Floor 2
# calculate nodal masses -- lump floor masses at frame nodes
set g 9.81; # acceleration due to gravity
set FloorWeight 297; # weight of Floor 2 in KN
set totmass [expr $FloorWeight/$g];
set NODEMASS [expr ($FloorWeight/$g) / (4)]; # mass at each node on Floor 1 - Totally 4 nodes exist
set Negligible 1e-9; # a very smnumber to avoid problems with zero
# define nodes and assign masses to beam-column intersections of frame
# command: node nodeID xcoord ycoord zcoord -mass mass_dof1 mass_dof2 mass_dof3
# nodeID convention: "abx" where a = Axis x and b = Axis z and x=story #
node 141 0 0 0;
node 131 0 0 6;
node 241 6 0 0;
node 231 6 0 6;
#
node 142 0 3.5 0 #-mass $NODEMASS $Negligible $NODEMASS $Negligible $Negligible $Negligible;
node 132 0 3.5 6 #-mass $NODEMASS $Negligible $NODEMASS $Negligible $Negligible $Negligible;
node 242 6 3.5 0 #-mass $NODEMASS $Negligible $NODEMASS $Negligible $Negligible $Negligible;
node 232 6 3.5 6 #-mass $NODEMASS $Negligible $NODEMASS $Negligible $Negligible $Negligible;

# column hinges at bottom of Story 1 (base)

node 1418 0 0 0
node 1318 0 0 6
node 2418 6 0 0
node 2318 6 0 6
#----------------------------------------------------------------
# column hinges at top of Story 1 (floor 1)

node 1429 0 3.5 0
node 1329 0 3.5 6
node 2429 6 3.5 0
node 2329 6 3.5 6
#----------------------------------------------------------------
# beam hinges at Story 1 (floor 2)

node 1424 0 3.5 0
node 1427 0 3.5 0
node 2425 6 3.5 0
node 2427 6 3.5 0
node 1326 0 3.5 6
node 1324 0 3.5 6
node 2325 6 3.5 6
node 2326 6 3.5 6

#########################################################
# define rigid floor diaphragm
set RigidDiaphragm ON;
set Xmid 3
set Zmid 3
set RigidDiaphragm ON;
node 1612 $Xmid 3.5 $Zmid -mass $totmass $Negligible $totmass $Negligible 181.7 $Negligible; #Master node for disphragm at story 1
#-----------------------------------------------------------------
#Constraints for rigid diaphragm master nodes
fix 1612 0 1 0 1 0 1
#-----------------------------------------------------------------
#Define Rigid Diaphragm, dof 2 is normal to floor
set perDir 2
rigidDiaphragm $perDir 1612 142 242 232 132
#----------------------------------------------------------------
# assign boundary condidtions
#Fix all nodes at y=0;
fix 131 1 1 1 1 1 1;
fix 141 1 1 1 1 1 1;
fix 231 1 1 1 1 1 1;
fix 241 1 1 1 1 1 1;
###################################################################################################
# Define Section Properties and Elements ###################################################################################################
# define material properties
set Ec 23979157; # Concrete Young's modulus for f'c=23 Mpa

set nu 0.2;
set Gc [expr $Ec/2./[expr 1+$nu]];
set Jcol 0.0036;

# define 1st story column section C408T16
set Acol_1 0.16; # cross-sectional area, m^2
set Izcol_1 2.13e-3; # moment of inertia about z-z, m^4
set Iycol_1 2.13e-3; # moment of inertia about y-y, m^4

set Mycol_1 71.51; # yield moment, KN.m
#------------------------------------------------------------------------------
# define beam section B30x40 for Floor 1 - Z Direction

set Abeam1_zdir 0.12; # cross-sectional area (full section properties), m^2
set Izbeam1_zdir 1.6e-3; # moment of inertia about z local axis for z-dir beams (full section properties), m^4
set Iybeam1_zdir 9e-4; # moment of inertia about y local axis for z-dir beams (full section properties), m^4
set Jbeamz 0.00194;

set Mybeam1_zdir 78.261; # yield moment at plastic hinge location about z local axis , KN.m
#------------------------------------------------------------------------------
# define beam section B30x30 for Floor 1 - x Direction

set Abeam1_xdir 0.09; # cross-sectional area (full section properties), m^2
set Jbeamx 0.00114;
set Izbeam1_xdir 6.75e-4; # moment of inertia about x local axis for x-dir beams (full section properties), m^4
set Iybeam1_xdir 6.75e-4; # moment of inertia about x local axis for x-dir beams (full section properties), m^4
set Mybeam1_xdir 14.71;
#set Mybeam1_xdir 14.71; # yield moment at plastic hinge location about x local axis , KN.m
#------------------------------------------------------------------------------
# determine stiffness modifications to equate the stiffness of the spring-elastic element-spring subassembly to the stiffness of the actual frame member
# calculate modified section properties to account for spring stiffness being in series with the elastic element stiffness
set n 10; # stiffness multiplier for rotational spring

# calculate modified moment of inertia for elastic elements
set Izcol_1mod [expr 0.39*$Izcol_1*($n+1.0)/$n]; # modified moment of inertia for columns in Story 1 about z-axis
set Iycol_1mod [expr 0.39*$Iycol_1*($n+1.0)/$n]; # modified moment of inertia for columns in Story 1 about z-axis
#
set Izbeam1_zdirmod [expr 0.655*$Izbeam1_zdir*($n+1.0)/$n]; # modified moment of inertia for beams in Floor 1 about z axis- z dir
set Iybeam1_zdirmod [expr 0.655*$Iybeam1_zdir*($n+1.0)/$n]; # modified moment of inertia for beams in Floor 1 about y axis- z dir

set Izbeam1_xdirmod [expr 0.80*$Izbeam1_xdir*($n+1.0)/$n]; # modified moment of inertia for beams in Floor 1 about z axis -x dir
set Iybeam1_xdirmod [expr 0.80*$Iybeam1_xdir*($n+1.0)/$n]; # modified moment of inertia for beams in Floor 1 about y axis-x dir
#
# calculate modified rotational stiffness for plastic hinge springs
#
set Ksz_col_1 [expr $n*6.0*$Ec*$Izcol_1mod/$HStory1]; # rotational stiffness of Story 1 column springs about z local
set Ksy_col_1 [expr $n*6.0*$Ec*$Iycol_1mod/$HStory1]; # rotational stiffness of Story 1 column springs about y local

set Ksz_beam1_z [expr $n*6.0*$Ec*$Izbeam1_zdirmod/6]; # rotational stiffness of Floor 1 beam springs about z local- z dir
set Ksz_beam1_x [expr $n*6.0*$Ec*$Izbeam1_xdirmod/6]; # rotational stiffness of Floor 1 beam springs about z local -x dir
#---------------------------------------------------------------------------------------------
# set up geometric transformations of element
set IDColTransf 1; #all columns
set IDBeamTransfz 2; # all beams in z direction
set IDBeamTransfx 3; # all beams in x direction
set ColTransfType PDelta

geomTransf $ColTransfType $IDColTransf 0 0 -1
geomTransf Linear $IDBeamTransfz 1 0 0
geomTransf Linear $IDBeamTransfx 0 0 1
# source rotSpring3DbeamZ.tcl;
#---------------------------------------------------------------------------------------------
# define elastic column elements using "element" command
# command: element elasticBeamColumn $eleID $iNode $jNode $A $E $G $J $Iy $Iz $transfID
# eleID convention: "9abcdef" where 9 = col, abc = node i ID#, def = node j ID#
# Columns Story 1
element elasticBeamColumn 9141142 1418 1429 $Acol_1 $Ec $Gc $Jcol $Iycol_1mod $Izcol_1mod $IDColTransf; # column 141 to 142
element elasticBeamColumn 9131132 1318 1329 $Acol_1 $Ec $Gc $Jcol $Iycol_1mod $Izcol_1mod $IDColTransf; # column 131 to 132
element elasticBeamColumn 9231232 2318 2329 $Acol_1 $Ec $Gc $Jcol $Iycol_1mod $Izcol_1mod $IDColTransf; # column 231 to 232
element elasticBeamColumn 9241242 2418 2429 $Acol_1 $Ec $Gc $Jcol $Iycol_1mod $Izcol_1mod $IDColTransf; # column 241 to 242
#----------------------------------------------------------------------------------------------------------------------------------
# define elastic beam elements
# command: element elasticBeamColumn $eleID $iNode $jNode $A $E $G $J $Iy $Iz $transfID
# eleID convention: "8abcdef" where 8 = beam, abc = node i ID#, def = node j ID#

# Beams in Z dir
element elasticBeamColumn 8132142 1326 1427 $Abeam1_zdir $Ec $Gc $Jbeamz $Iybeam1_zdirmod $Izbeam1_zdirmod $IDBeamTransfz; # beam 132 to 142
element elasticBeamColumn 8232242 2326 2427 $Abeam1_zdir $Ec $Gc $Jbeamz $Iybeam1_zdirmod $Izbeam1_zdirmod $IDBeamTransfz; # beam 232 to 242

# Beams in x dir
element elasticBeamColumn 8142242 1424 2425 $Abeam1_xdir $Ec $Gc $Jbeamx $Iybeam1_xdirmod $Izbeam1_xdirmod $IDBeamTransfx; # beam 142 to 242
element elasticBeamColumn 8132232 1324 2325 $Abeam1_xdir $Ec $Gc $Jbeamx $Iybeam1_xdirmod $Izbeam1_xdirmod $IDBeamTransfx; # beam 132 to 232

#----------------------------------------------------------------
###################################################################################################
# Define Rotational Springs for Plastic Hinges
###################################################################################################
# define rotational spring properties and create spring elements using "rotSpring3DCol" procedure
# input values for Story 1 column springs
set McMy 1.13; # ratio of capping moment to yield moment, Mc / My
set LS 68.12; # basic strength deterioration (a very large # = no cyclic deterioration)
set LK 68.12; # unloading stiffness deterioration (a very large # = no cyclic deterioration)
set LA 68.12; # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
set LD 100.0; # post-capping strength deterioration (a very large # = no deterioration)
set cS 1.0; # exponent for basic strength deterioration (c = 1.0 for no deterioration)
set cK 1.0; # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
set cA 1.0; # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
set cD 1.0; # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)
set th_pP 0.0304; # plastic rot capacity for pos loading
set th_pN 0.0304; # plastic rot capacity for neg loading
set th_pcP 0.0399; # post-capping rot capacity for pos loading
set th_pcN 0.0399; # post-capping rot capacity for neg loading
set ResP 0.1; # residual strength ratio for pos loading
set ResN 0.1; # residual strength ratio for neg loading
set th_uP 0.0775; # ultimate rot capacity for pos loading
set th_uN 0.0775; # ultimate rot capacity for neg loading
set DP 1.0; # rate of cyclic deterioration for pos loading
set DN 1.0; # rate of cyclic deterioration for neg loading
set a_mem [expr ($n+1.0)*($Mycol_1*($McMy-1.0)) / ($Ksy_col_1*$th_pP)]; # strain hardening ratio of spring
set b [expr ($a_mem)/(1.0+$n*(1.0-$a_mem))]; # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: Eqn B.5 is incorrect)
# define column springs
# Spring ID: "8abcdef", where 8 = col spring about y, abc = Master Node #, defg=seconf node #
# command: rotSpring3DCol id ndR ndC K asPos asNeg MyPos MyNeg LS LK LA LD cS cK cA cD th_p+ th_p- th_pc+ th_pc- Res+ Res- th_u+ th_u- D+ D-
###############################################
# col springs @ bottom of Story 1 (at base)
#uniaxialMaterial Elastic 999 23979157
################################################

uniaxialMaterial Bilin 999 $Ksy_col_1 $b $b $Mycol_1 [expr -$Mycol_1] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
element zeroLength 81311318 131 1318 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0

# Constrain the translational DOF with a multi-point constraint
# retained constrained DOF_1 DOF_2 ... DOF_n
# #
equalDOF 131 1318 1 2 3 5

element zeroLength 81411418 141 1418 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0

# Constrain the translational DOF with a multi-point constraint
# retained constrained DOF_1 DOF_2 ... DOF_n
# #
equalDOF 141 1418 1 2 3 5

element zeroLength 82412418 241 2418 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0

# Constrain the translational DOF with a multi-point constraint
# retained constrained DOF_1 DOF_2 ... DOF_n
# #
equalDOF 241 2418 1 2 3 5

element zeroLength 82312318 231 2318 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0

# Constrain the translational DOF with a multi-point constraint
# retained constrained DOF_1 DOF_2 ... DOF_n
# #
equalDOF 231 2318 1 2 3 5
#######################################################################################top#####################
element zeroLength 81321329 1329 132 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0
equalDOF 1329 132 1 2 3 5


element zeroLength 81421429 1429 142 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0
equalDOF 1429 142 1 2 3 5

element zeroLength 82422429 2429 242 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0
equalDOF 2429 242 1 2 3 5


element zeroLength 82322329 2329 232 -mat 999 999 -dir 5 6 -orient 0 1 0 1 0 0
equalDOF 2329 232 1 2 3 5


#rotSpring3DCol 8232 2329 232 $Ksy_col_1 $b $b $Mycol_1 [expr -$Mycol_1] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
#-------------------------------
# create region for frame column springs
# command: region $regionID -ele $ele_1_ID $ele_2_ID...
region 1 -ele 81311318 81411418 82412418 82312318 81321329 81421429 82422429 82322329 ;
############################################################################################################
# define rotational spring properties for 1st story beams and create spring elements using "rotSpring3DbeamX" procedure

# input values for Story 1 beams prings-----### x ###### direction ---
set McMy 1.13; # ratio of capping moment to yield moment, Mc / My
set LS 83.67; # basic strength deterioration (a very large # = no cyclic deterioration)
set LK 83.67; # unloading stiffness deterioration (a very large # = no cyclic deterioration)
set LA 83.67; # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
set LD 100.0; # post-capping strength deterioration (a very large # = no deterioration)
set cS 1.0; # exponent for basic strength deterioration (c = 1.0 for no deterioration)
set cK 1.0; # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
set cA 1.0; # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
set cD 1.0; # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)
set th_pP 0.0336; # plastic rot capacity for pos loading
set th_pN 0.0336; # plastic rot capacity for neg loading
set th_pcP 0.0519; # post-capping rot capacity for pos loading
set th_pcN 0.0519; # post-capping rot capacity for neg loading
set ResP 0.1; # residual strength ratio for pos loading
set ResN 0.1; # residual strength ratio for neg loading
set th_uP 0.0969; # ultimate rot capacity for pos loading
set th_uN 0.0969; # ultimate rot capacity for neg loading
set DP 1.0; # rate of cyclic deterioration for pos loading
set DN 1.0; # rate of cyclic deterioration for neg loading
set a_mem [expr ($n+1.0)*($Mybeam1_xdir*($McMy-1.0)) / ($Ksz_beam1_x*$th_pP)]; # strain hardening ratio of spring
set b [expr ($a_mem)/(1.0+$n*(1.0-$a_mem))]; # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: Eqn B.5 is incorrect)

# define beam springs @ Story 1 (Floor2)
# Spring ID: "6abcdefg", where 6 = beam spring for x-dir beams, abc = master #, defg = second node #

uniaxialMaterial Bilin 888 $Ksz_beam1_x $b $b $Mybeam1_xdir [expr -$Mybeam1_xdir] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;

element zeroLength 61321324 132 1324 -mat 888 -dir 6 -orient 1 0 0 0 1 0
equalDOF 132 1324 1 2 3 4 5

element zeroLength 62322325 2325 232 -mat 888 -dir 6 -orient 1 0 0 0 1 0
equalDOF 232 2325 1 2 3 4 5
element zeroLength 62422425 2425 242 -mat 888 -dir 6 -orient 1 0 0 0 1 0
equalDOF 242 2425 1 2 3 4 5

element zeroLength 61421424 142 1424 -mat 888 -dir 6 -orient 1 0 0 0 1 0
equalDOF 142 1424 1 2 3 4 5

############################################################################################################
# define rotational spring properties for 1st story beams and create spring elements using "rotSpring3DbeamZ" procedure

# input values for Story 1 beams prings-----### z ###### direction ---
set McMy 1.13; # ratio of capping moment to yield moment, Mc / My
set LS 92.19; # basic strength deterioration (a very large # = no cyclic deterioration)
set LK 92.19; # unloading stiffness deterioration (a very large # = no cyclic deterioration)
set LA 92.19; # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
set LD 100.0; # post-capping strength deterioration (a very large # = no deterioration)
set cS 1.0; # exponent for basic strength deterioration (c = 1.0 for no deterioration)
set cK 1.0; # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
set cA 1.0; # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
set cD 1.0; # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)

set th_pP 0.0354; # plastic rot capacity for pos loading
set th_pN 0.0354; # plastic rot capacity for neg loading
set th_pcP 0.05191; # post-capping rot capacity for pos loading
set th_pcN 0.05191; # post-capping rot capacity for neg loading
set ResP 0.1; # residual strength ratio for pos loading
set ResN 0.1; # residual strength ratio for neg loading
set th_uP 0.097; # ultimate rot capacity for pos loading
set th_uN 0.097; # ultimate rot capacity for neg loading
set DP 1.0; # rate of cyclic deterioration for pos loading
set DN 1.0; # rate of cyclic deterioration for neg loading
set a_mem [expr ($n+1.0)*($Mybeam1_zdir*($McMy-1.0)) / ($Ksz_beam1_z*$th_pP)]; # strain hardening ratio of spring
set b [expr ($a_mem)/(1.0+$n*(1.0-$a_mem))]; # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: Eqn B.5 is incorrect)

# define beam springs @ Story 1 (Floor2)
# Spring ID: "7abcdefg", where 6 = beam spring for z-dir beams, abc = master #, defg = second node #

uniaxialMaterial Bilin 777 $Ksz_beam1_z $b $b $Mybeam1_zdir [expr -$Mybeam1_zdir] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
element zeroLength 61321326 132 1326 -mat 777 -dir 6 -orient 0 0 -1 0 1 0
equalDOF 132 1326 1 2 3 5 6

element zeroLength 61421427 1427 142 -mat 777 -dir 6 -orient 0 0 -1 0 1 0
equalDOF 1427 142 1 2 3 5 6

element zeroLength 62422427 2427 242 -mat 777 -dir 6 -orient 0 0 -1 0 1 0
equalDOF 2427 242 1 2 3 5 6

element zeroLength 62322326 232 2326 -mat 777 -dir 6 -orient 0 0 -1 0 1 0
equalDOF 232 2326 1 2 3 5 6

#---------------------------------------------------------------------------------------------
# create region for beam springs
region 2 -ele 61321324 62322325 62422425 61421424 61321326 61421427 62422427 62322326;

recorder Element -file $dataDir/Fgravs.out -ele 9241242 force;
#############################################################################
# Gravity Loads & Gravity Analysis
############################################################################
# apply gravity loads
#command: pattern PatternType $PatternID TimeSeriesType
pattern Plain 101 Constant {

# point loads on frame column nodes
# Floor 1 loads
set DL_Load [expr -78.25];

load 132 0.0 $DL_Load 0.0 0.0 0.0 0.0;
load 232 0.0 $DL_Load 0.0 0.0 0.0 0.0;
load 142 0.0 $DL_Load 0.0 0.0 0.0 0.0;
load 242 0.0 $DL_Load 0.0 0.0 0.0 0.0;
}
# Gravity-analysis: load-controlled static analysis
set Tol 1.0e-5;
# convergence tolerance for test
constraints Penalty 10e9 10e9;
numberer RCM;
system BandGeneral;
test EnergyIncr $Tol 600;
algorithm Newton;
set NstepGravity 10;
set DGravity [expr 1./$NstepGravity];
integrator LoadControl $DGravity;
analysis Static;
analyze $NstepGravity;

puts "Model Built"


############################################################################
# Eigenvalue Analysis
############################################################################
set pi [expr 2.0*asin(1.0)]; # Definition of pi
set nEigenI 1; # mode i = 1
set nEigenJ 3; # mode j = 3
set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr 0]]; # eigenvalue mode i = 1
set lambdaJ [lindex $lambdaN [expr $nEigenJ-2]]; # eigenvalue mode j = 2
set lambdak [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j = 2
set w1 [expr pow($lambdaI,0.5)]; # w1 (1st mode circular frequency)
set w2 [expr pow($lambdaJ,0.5)]; # w2 (2nd mode circular frequency)
set w3 [expr pow($lambdak,0.5)]; # w2 (2nd mode circular frequency)set w2 [expr pow($lambdaJ,0.5)]; # w2 (2nd mode circular frequency)
set T1 [expr 2.0*$pi/$w1]; # 1st mode period of the structure
set T2 [expr 2.0*$pi/$w2]; # 2nd mode period of the structure
set T3 [expr 2.0*$pi/$w3]; # 2nd mode period of the structure
puts "T1 = $T1 s"; # display the first mode period in the command window
puts "T2 = $T2 s";
puts "T3 = $T3 s";
############################################################################
############################################################################
# Recorders
############################################################################
# record drift histories
# drift recorder command: recorder Drift -file $filename -iNode $NodeI_ID -jNode $NodeJ_ID -dof $dof -perpDirn $Record.drift.perpendicular.to.this.direction
recorder Drift -file $dataDir/Drift-132.out -iNode 131 -jNode 132 -dof 3 -perpDirn 2;
recorder Drift -file $dataDir/Drift-142.out -iNode 141 -jNode 142 -dof 3 -perpDirn 2;

# record base shear reactions
recorder Node -file $dataDir/Vbase.out -node 1318 1418 2318 2418 -dof 3 reaction;

# record story 1 column forces in global coordinates
recorder Element -file $dataDir/Fcol111.out -ele 9131132 force;
recorder Element -file $dataDir/Fcol121.out -ele 9141142 force;
recorder Element -file $dataDir/Fcol131.out -ele 9231232 force;
recorder Element -file $dataDir/Fcol141.out -ele 9241242 force;
# record response history of all frame column springs (one file for moment, one for rotation)
recorder Element -file $dataDir/MRFcol-Mom-Hist.out -region 1 force;
recorder Element -file $dataDir/MRFcol-Rot-Hist.out -region 1 deformation;

# record response history of all frame beam springs (one file for moment, one for rotation)
recorder Element -file $dataDir/MRFbeam-Mom-Hist.out -region 2 force;
recorder Element -file $dataDir/MRFbeam-Rot-Hist.out -region 2 deformation;

#######################################################################################
# #
# Analysis Section #
# #
#######################################################################################


puts "Running time history..."

recorder display "Mode Shape 1" 10 510 200 200 -wipe
prp $HStory1 $HStory1 1;
vup 0 1 0;
vpn 0 0 0;
viewWindow -950 950 -950 950
display -1 5 20

recorder display "Mode Shape 2" 210 510 200 200 -wipe
prp $HStory1 $HStory1 1;
vup 0 1 0;
vpn 0 0 1;
viewWindow -650 650 -650 650
display -2 5 20

DisplayModel3D DeformedShape ; # options: DeformedShape NodeNumbers ModeSh


set accelSeriesN "Path -filePath tabasFN.txt -dt 0.02 -factor $g"

set accelSeriesP "Path -filePath tabasFP.txt -dt 0.02 -factor $g"

# Define the ground motion excitation using Tabas fault parallel and fault normal records

# tag dir accel series args

pattern UniformExcitation 2 1 -accel $accelSeriesN

pattern UniformExcitation 3 2 -accel $accelSeriesP


# Record DOF 1 and 2 displacements at nodes 9, 14, and 19

#recorder Node RigidFrame3D.out disp -time -node 1612 -dof 1 2


# Source in commands to display the structure

source RigidFrame3Ddisplay.tcl



# DOF numberer

numberer RCM

# Constraint handler

# aSP aMP

#constraints Lagrange 1.0e6 1.0e6

constraints Penalty 10e9 10e9;

# Convergence test

# tol maxIter printFlag

test EnergyIncr 1.0e-8 20 1

# Solution algorithm

algorithm Newton
# System of equations solver

#system SparseGeneral -piv

system BandGeneral


# Transient integrator

# gamma beta

integrator Newmark 0.5 0.25

# Transient analysis

# dt T

analysis Transient

# Perform the analysis

#analyze 1000 0.01

# set some variables
set tFinal [expr 1000 * 0.01]
set tCurrent [getTime]
set ok 0

# Perform the transient analysis
while {$ok == 0 && $tCurrent < $tFinal} {

set ok [analyze 1 .01]

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

set tCurrent [getTime]
}

# Print a message to indicate if analysis succesfull or not
if {$ok == 0} {
puts "Transient analysis completed SUCCESSFULLY";
} else {
puts "Transient analysis completed FAILED";
}
KARIMIYAN
Post Reply