gravity static analysis problem?

If you have a script you think might be useful to others post it
here. Hopefully we will be able to get the most useful of these incorporated in the manuals.

Moderators: silvia, selimgunay, Moderators

Post Reply
khili83
Posts: 4
Joined: Sat Aug 09, 2008 10:03 am
Location: tufts university

gravity static analysis problem?

Post by khili83 » Sat Aug 09, 2008 10:19 am

I am modeling a three story irregular building? it is strange that it fails under static force controled analysis for gravity analysis, but the model fulfills the dynamic analysis ?!!!!!!!!!!
*************************************

# --------------------------------------------------------------------------------------------------
# Example 8. 3D Steel W-section frame
# Silvia Mazzoni & Frank McKenna, 2006
# nonlinearBeamColumn element, inelastic fiber section
#

# SET UP ----------------------------------------------------------------------------
wipe; # clear memory of all past model definitions
model BasicBuilder -ndm 3 -ndf 6; # Define the model builder, ndm=#dimension, ndf=#dofs
set dataDir Data; # set up name of data directory
file mkdir $dataDir; # create data directory
set GMdir "GMfiles"; # ground-motion file directory
source LibUnits.tcl; # define units
source DisplayPlane.tcl; # procedure for displaying a plane in model
source DisplayModel3D.tcl; # procedure for displaying 3D perspectives of model
source Wsection.tcl; # procedure to define fiber W section

# ------ frame configuration
set NStory 3; # number of stories above ground level
set NBay 3; # number of bays in X direction
set NBayZ 3; # number of bays in Z direction
puts "Number of Stories in Y: $NStory Number of bays in X: $NBay Number of bays in Z: $NBayZ"
set NFrame [expr $NBayZ + 1]; # actually deal with frames in Z direction, as this is an easy extension of the 2d model

# define GEOMETRY -------------------------------------------------------------
# define structure-geometry paramters
set LCol [expr 10.4987*$ft]; # column height (parallel to Y axis)
set LBeam [expr 19.685*$ft]; # beam length (parallel to X axis)
set LGird [expr 19.685*$ft]; # girder length (parallel to Z axis)

# define NODAL COORDINATES
set Dlevel 10000; # numbering increment for new-level nodes
set Dframe 100; # numbering increment for new-frame nodes
for {set level 1} {$level <=[expr $NStory+1]} {incr level 1} {
if {$level==1} {
set a 1}
if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}
for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
set Z [expr ($frame-1)*$LGird]
set Y [expr ($level-1)*$LCol];
for {set pier $a} {$pier <= [expr $NBay+1]} {incr pier 1} {
set X [expr ($pier-1)*$LBeam];
set nodeID [expr $frame*$Dframe+$level*$Dlevel+$pier]
node $nodeID $X $Y $Z; # actually define node
}
}
}

# rigid diaphragm nodes
set RigidDiaphragm ON ; # options: ON, OFF. specify this before the analysis parameters are set the constraints are handled differently.
set iMasterNode ""
for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {

set Y [expr ($level-1)*$LCol];

# rigid-diaphragm nodes in center of each diaphram
set MasterNodeID [expr 9900+$level]

if {$level==2} {
set Xa [expr ($NBay*$LBeam)/2]; # mid-span coordinate for rigid diaphragm
set Za [expr ($NFrame-1)*$LGird/2]; }
if {$level==3} {
set Xa [expr (($NBay-1)*$LBeam)]; # mid-span coordinate for rigid diaphragm
set Za [expr ($NFrame-1)*$LGird/2]; }
if {$level==4} {
set Xa [expr (($NBay-1/2)*$LBeam)]; # mid-span coordinate for rigid diaphragm
set Za [expr ($NFrame-1)*$LGird/2]; }

node $MasterNodeID $Xa $Y $Za; # master nodes for rigid diaphragm

if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}
for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
for {set pier $a} {$pier <= [expr $NBay+1]} {incr pier 1} {

set nodeID [expr $frame*$Dframe+$level*$Dlevel+$pier]

fix $MasterNodeID 0 1 0 1 0 1; # constrain other dofs that don't belong to rigid diaphragm control
lappend iMasterNode $MasterNodeID
set perpDirn 2; # perpendicular to plane of rigid diaphragm
rigidDiaphragm $perpDirn $MasterNodeID $nodeID; # define Rigid Diaphram,


}
}

}


# determine support nodes where ground motions are input, for multiple-support excitation
set iSupportNode ""
for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
set level 1
for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} {
set nodeID [expr $level*$Dlevel+$frame*$Dframe+$pier]
lappend iSupportNode $nodeID
}
}

# BOUNDARY CONDITIONS
fixY 0.0 1 1 1 1 1 1; # pin all Y=0.0 nodes

# calculated MODEL PARAMETERS, particular to this model
# Set up parameters that are particular to the model for displacement control
set IDctrlNode [expr int((($NStory+1)*$Dlevel+$NFrame*$Dframe)+1)]; # node where displacement is read for displacement control
set IDctrlDOF 1; # degree of freedom of displacement read for displacement control
set LBuilding [expr $NStory*$LCol]; # total building height

# Define SECTIONS -------------------------------------------------------------
set SectionType Elastic; # options: Elastic FiberSection

# define section tags:
set Col1SecTag 1
set Col2SecTag 2
set Col3SecTag 3
set Col4SecTag 4
set Col5SecTag 5
set Col6SecTag 6
set Col7SecTag 7

set Beam1SecTag 8
set Beam2SecTag 9
set Beam3SecTag 10
set Beam4SecTag 11
set Beam5SecTag 12
set Beam6SecTag 13
set Beam7SecTag 14


set ColSecTagFiber 15
set BeamSecTagFiber 16
set GirdSecTagFiber 17
set SecTagTorsion 18
set Ubig 10e10;
if {$SectionType == "Elastic"} {
# material properties:
set Es [expr 29000*$ksi]; # Steel Young's Modulus
set nu 0.3; # Poisson's ratio
set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus
set J $Ubig; # set large torsional stiffness
# column1 sections: Box 320T15
set Ag1Col [expr 28.365*pow($in,2)]; # cross-sectional area
set Iz1Col [expr 683.304*pow($in,4)]; # moment of Inertia
set Iy1Col [expr 683.304*pow($in,4)]; # moment of Inertia
# column2 sections: Box 300T15
set Ag2Col [expr 26.51*pow($in,2)]; # cross-sectional area
set Iz2Col [expr 557.7*pow($in,4)]; # moment of Inertia
set Iy2Col [expr 557.7*pow($in,4)]; # moment of Inertia
# column3 sections: Box 280T15
set Ag3Col [expr 24.65*pow($in,2)]; # cross-sectional area
set Iz3Col [expr 448.53*pow($in,4)]; # moment of Inertia
set Iy3Col [expr 448.53*pow($in,4)]; # moment of Inertia
# column4 sections: Box 260T15
set Ag4Col [expr 22.79*pow($in,2)]; # cross-sectional area
set Iz4Col [expr 354.64*pow($in,4)]; # moment of Inertia
set Iy4Col [expr 354.64*pow($in,4)]; # moment of Inertia
# column5 sections: Box 240T15
set Ag5Col [expr 20.93*pow($in,2)]; # cross-sectional area
set Iz5Col [expr 274.78*pow($in,4)]; # moment of Inertia
set Iy5Col [expr 274.78*pow($in,4)]; # moment of Inertia
# column6 sections: Box 220T15
set Ag6Col [expr 19.07*pow($in,2)]; # cross-sectional area
set Iz6Col [expr 208.09*pow($in,4)]; # moment of Inertia
set Iy6Col [expr 208.09*pow($in,4)]; # moment of Inertia
# column7 sections: Box 200T15
set Ag7Col [expr 17.21*pow($in,2)]; # cross-sectional area
set Iz7Col [expr 153.12*pow($in,4)]; # moment of Inertia
set Iy7Col [expr 153.12*pow($in,4)]; # moment of Inertia


# beam1 sections: HE320A
set Ag1Beam [expr 14.663*pow($in,2)]; # cross-sectional area
set Iz1Beam [expr 395.213*pow($in,4)]; # moment of Inertia
set Iy1Beam [expr 119.14*pow($in,4)]; # moment of Inertia
# beam2 sections: HE300A
set Ag2Beam [expr 13.7795*pow($in,2)]; # cross-sectional area
set Iz2Beam [expr 331.55*pow($in,4)]; # moment of Inertia
set Iy2Beam [expr 113.735*pow($in,4)]; # moment of Inertia

# beam3 sections: HE280A
set Ag3Beam [expr 12.09*pow($in,2)]; # cross-sectional area
set Iz3Beam [expr 253.71*pow($in,4)]; # moment of Inertia
set Iy3Beam [expr 88.03*pow($in,4)]; # moment of Inertia

# beam4 sections: HE260A
set Ag4Beam [expr 10.695*pow($in,2)]; # cross-sectional area
set Iz4Beam [expr 191.744*pow($in,4)]; # moment of Inertia
set Iy4Beam [expr 66.982*pow($in,4)]; # moment of Inertia

# beam5 sections: HE240A
set Ag5Beam [expr 9.36*pow($in,2)]; # cross-sectional area
set Iz5Beam [expr 129.38*pow($in,4)]; # moment of Inertia
set Iy5Beam [expr 49.9*pow($in,4)]; # moment of Inertia

# beam6 sections: HE220A
set Ag6Beam [expr 7.98*pow($in,2)]; # cross-sectional area
set Iz6Beam [expr 100.19*pow($in,4)]; # moment of Inertia
set Iy6Beam [expr 36.28*pow($in,4)]; # moment of Inertia

# beam7 sections: HE200A
set Ag7Beam [expr 6.84*pow($in,2)]; # cross-sectional area
set Iz7Beam [expr 70.73*pow($in,4)]; # moment of Inertia
set Iy7Beam [expr 25.66*pow($in,4)]; # moment of Inertia

section Elastic $Col1SecTag $Es $Ag1Col $Iz1Col $Iy1Col $Gs $J
section Elastic $Col2SecTag $Es $Ag2Col $Iz2Col $Iy2Col $Gs $J
section Elastic $Col3SecTag $Es $Ag3Col $Iz3Col $Iy3Col $Gs $J
section Elastic $Col4SecTag $Es $Ag4Col $Iz4Col $Iy4Col $Gs $J
section Elastic $Col5SecTag $Es $Ag5Col $Iz5Col $Iy5Col $Gs $J
section Elastic $Col6SecTag $Es $Ag6Col $Iz6Col $Iy6Col $Gs $J
section Elastic $Col7SecTag $Es $Ag7Col $Iz7Col $Iy7Col $Gs $J

section Elastic $Beam1SecTag $Es $Ag1Beam $Iz1Beam $Iy1Beam $Gs $J
section Elastic $Beam2SecTag $Es $Ag2Beam $Iz2Beam $Iy2Beam $Gs $J
section Elastic $Beam3SecTag $Es $Ag3Beam $Iz3Beam $Iy3Beam $Gs $J
section Elastic $Beam4SecTag $Es $Ag4Beam $Iz4Beam $Iy4Beam $Gs $J
section Elastic $Beam5SecTag $Es $Ag5Beam $Iz5Beam $Iy5Beam $Gs $J
section Elastic $Beam6SecTag $Es $Ag6Beam $Iz6Beam $Iy6Beam $Gs $J
section Elastic $Beam7SecTag $Es $Ag7Beam $Iz7Beam $Iy7Beam $Gs $J

set matIDhard 1; # material numbers for recorder (this stressstrain recorder will be blank, as this is an elastic section)

} elseif {$SectionType == "FiberSection"} {
# define MATERIAL properties
set Fy [expr 60.0*$ksi]
set Es [expr 29000*$ksi]; # Steel Young's Modulus
set nu 0.3;
set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus
set Hiso 0
set Hkin 1000
set matIDhard 1
uniaxialMaterial ElasticPP $matIDhard $Es [expr $Fy/$Es]
# ELEMENT properties
# Structural-Steel W-section properties
# column sections: box
set bs [expr 8.66142*$in]; # depth
set bt [expr 8.66142*$in]; # width
set tt [expr 3.81*$in]; # top thickness
set ts [expr 3.81*$in]; # side thickness
set nftl 16; # number of fibers along dw
set nftt 2; # number of fibers along tw
set nfsl 16; # number of fibers along bf
set nfst 4; # number of fibers along tf
boxsection $ColSecTagFiber $matIDhard $bs $bt $tt $ts $nfsl $nfst $nftl $nftt
# beam sections: HE200A
set d [expr 47.244*$in]; # depth
set bf [expr 50.8*$in]; # flange width
set tf [expr 2.032*$in]; # flange thickness
set tw [expr 1.397*$in]; # web thickness
set nfdw 16; # number of fibers along dw
set nftw 2; # number of fibers along tw
set nfbf 16; # number of fibers along bf
set nftf 4; # number of fibers along tf
Wsection $BeamSecTagFiber $matIDhard $d $bf $tf $tw $nfdw $nftw $nfbf $nftf
# girder sections: HE240A
set d [expr 56.896*$in]; # depth
set bf [expr 60.96*$in]; # flange width
set tf [expr 2.286*$in]; # flange thickness
set tw [expr 1.651*$in]; # web thickness
set nfdw 16; # number of fibers along dw
set nftw 2; # number of fibers along tw
set nfbf 16; # number of fibers along bf
set nftf 4; # number of fibers along tf
Wsection $GirdSecTagFiber $matIDhard $d $bf $tf $tw $nfdw $nftw $nfbf $nftf

# assign torsional Stiffness for 3D Model
uniaxialMaterial Elastic $SecTagTorsion $Ubig
section Aggregator $ColSecTag $SecTagTorsion T -section $ColSecTagFiber
section Aggregator $BeamSecTag $SecTagTorsion T -section $BeamSecTagFiber
section Aggregator $GirdSecTag $SecTagTorsion T -section $GirdSecTagFiber
} else {
puts "No section has been defined"
return -1
}

set QdlCol1 [expr 96.565*$lbf/$ft]; # HE-section weight per length
set QBeam1 [expr 49.9*$lbf/$ft]; # HE-section weight per length
set QGird1 [expr 49.9*$lbf/$ft]; # HE-section weight per length
set QdlCol2 [expr 90.23*$lbf/$ft]; # HE-section weight per length
set QBeam2 [expr 46.9*$lbf/$ft]; # HE-section weight per length
set QGird2 [expr 46.9*$lbf/$ft]; # HE-section weight per length
set QdlCol3 [expr 83.9*$lbf/$ft]; # HE-section weight per length
set QBeam3 [expr 41.15*$lbf/$ft]; # HE-section weight per length
set QGird3 [expr 41.15*$lbf/$ft]; # HE-section weight per length
set QdlCol4 [expr 77.57*$lbf/$ft]; # HE-section weight per length
set QBeam4 [expr 36.41*$lbf/$ft]; # HE-section weight per length
set QGird4 [expr 36.41*$lbf/$ft]; # HE-section weight per length
set QdlCol5 [expr 71.24*$lbf/$ft]; # HE-section weight per length
set QBeam5 [expr 31.87*$lbf/$ft]; # HE-section weight per length
set QGird5 [expr 31.87*$lbf/$ft]; # HE-section weight per length
set QdlCol6 [expr 64.9*$lbf/$ft]; # HE-section weight per length
set QBeam6 [expr 27.18*$lbf/$ft]; # HE-section weight per length
set QGird6 [expr 27.18*$lbf/$ft]; # HE-section weight per length
set QdlCol7 [expr 58.57*$lbf/$ft]; # HE-section weight per length
set QBeam7 [expr 23.27*$lbf/$ft]; # HE-section weight per length
set QGird7 [expr 23.27*$lbf/$ft]; # HE-section weight per length


# --------------------------------------------------------------------------------------------------------------------------------
# define ELEMENTS
# set up geometric transformations of element
# separate columns and beams, in case of P-Delta analysis for columns
set IDColTransf 1; # all columns
set IDBeamTransf 2; # all beams
set IDGirdTransf 3; # all girds
set ColTransfType PDelta; # options for columns: Linear PDelta Corotational
geomTransf $ColTransfType $IDColTransf 0 0 1; # orientation of column stiffness affects bidirectional response.
geomTransf PDelta $IDBeamTransf 0 0 1
geomTransf PDelta $IDGirdTransf 1 0 0

# Define Beam-Column Elements
set numIntgrPts 5; # number of Gauss integration points for nonlinear curvature distribution
# columns
set N0col [expr 10000-1]; # column element numbers
set level 0

for {set level 1} {$level <=$NStory} {incr level 1} {
if {$level==1} {
set a 1}
if {$level==2} {
set a 2}
if {$level==3} {
set a 3}

for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
for {set pier $a} {$pier <= [expr $NBay+1]} {incr pier 1} {
set elemID [expr $N0col +$level*$Dlevel + $frame*$Dframe+$pier]
set nodeI [expr $level*$Dlevel + $frame*$Dframe+$pier]
set nodeJ [expr ($level+1)*$Dlevel + $frame*$Dframe+$pier]
if {$frame==1 || $frame==4} {
if {$level==1} {
if {$pier==1 || $pier==4} {
set ColSecTag $Col5SecTag;
} elseif {$pier==2 || $pier==3} {
set ColSecTag $Col3SecTag;}
} elseif {$level==2} {
if {$pier==2 || $pier==3} {
set ColSecTag $Col4SecTag;
} elseif {$pier==4} {
set ColSecTag $Col6SecTag;}
} elseif {$level==3} {
if {$pier==3} {
set ColSecTag $Col5SecTag;
} elseif {$pier==4} {
set ColSecTag $Col6SecTag;}
}
} elseif {$frame==2 || $frame==3} {
if {$level==1} {
set ColSecTag $Col3SecTag;
} elseif {$level==2} {
set ColSecTag $Col4SecTag;
} elseif {$level==3} {
set ColSecTag $Col5SecTag;
}
}

element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $ColSecTag $IDColTransf; # columns
}
}
}

# beams -- side frames
set N0beam 1000000; # beam element numbers
set N0gird 2000000; # gird element numbers

for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {
for {set frame 1} {$frame <= $NFrame} {incr frame 1} {
if {$level==1} {
set a 1}
if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}

for {set bay $a} {$bay <= $NBay} {incr bay 1} {
set elemID [expr $N0beam +$level*$Dlevel + $frame*$Dframe+ $bay]
set nodeI [expr $level*$Dlevel + $frame*$Dframe+ $bay]
set nodeJ [expr $level*$Dlevel + $frame*$Dframe+ $bay+1]
if {$frame==1 || $frame==4} {
if {$level==2} {
set BeamSecTag $Beam5SecTag;
} elseif {$level==3} {
set BeamSecTag $Beam6SecTag;
} elseif {$level==4} {
set BeamSecTag $Beam7SecTag; }
} elseif {$frame==2 || $frame==3} {
if {$level==2} {
set BeamSecTag $Beam3SecTag;
} elseif {$level==3} {
set BeamSecTag $Beam4SecTag;
} elseif {$level==4} {
set BeamSecTag $Beam5SecTag;
}
}

element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $BeamSecTag $IDBeamTransf; # beams
}
}
}

for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {

if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}

for {set bay $a} {$bay <=[expr $NBay+1]} {incr bay 1} {
for {set frame 1} {$frame <=[expr $NFrame-1]} {incr frame 1} {
set elemID [expr $N0gird + $level*$Dlevel +$frame*$Dframe+ $bay]
set nodeI [expr $level*$Dlevel + $frame*$Dframe+ $bay]
set nodeJ [expr $level*$Dlevel + ($frame+1)*$Dframe+ $bay]

if {$bay==1} {
if {$level==2} {
set BeamSecTag $Beam5SecTag;}
} elseif {$bay==2} {
if {$level==2} {
set BeamSecTag $Beam3SecTag;
} elseif {$level==3} {
set BeamSecTag $Beam4SecTag;}
} elseif {$bay==3} {
if {$level==2} {
set BeamSecTag $Beam3SecTag;
} elseif {$level==3} {
set BeamSecTag $Beam4SecTag;
} elseif {$level==4} {
set BeamSecTag $Beam5SecTag;}
} elseif {$bay==4} {
if {$level==2} {
set BeamSecTag $Beam5SecTag;
} elseif {$level==3} {
set BeamSecTag $Beam6SecTag;
} elseif {$level==4} {
set BeamSecTag $Beam7SecTag;
}
}


element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $BeamSecTag $IDGirdTransf; # Girds
}
}
}

# --------------------------------------------------------------------------------------------------------------------------------
# Define GRAVITY LOADS, weight and masses
# calculate dead load of frame, assume this to be an internal frame (do LL in a similar manner)
# calculate distributed weight along the beam length

# assign masses to the nodes that the columns are connected to
# each connection takes the mass of 1/2 of each element framing into it (mass=weight/$g)
set iFloorWeight ""
set WeightTotal 0.0
set sumWiHi 0.0; # sum of storey weight times height, for lateral-load distribution

for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {
set FloorWeight 0.0;

if {$level == [expr $NStory+1]} {
set ColWeightFact 1; # one column in top story
} else {
set ColWeightFact 1; # two columns elsewhere
}


for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {

if {$frame == 1 || $frame == $NFrame} {
set GirdWeightFact 1; # 1x1/2girder on exterior frames
} else {
set GirdWeightFact 2; # 2x1/2girder on interior frames
}


if {$level==2} {
set a 1}
if {$level == 3} {
set a 2}
if {$level==4} {
set a 3}

for {set pier $a} {$pier <= [expr $NBay+1]} {incr pier 1} {

if {$pier == $a || $pier == [expr $NBay+1]} {
set BeamWeightFact 1; # one beam at exterior nodes
} else {
set BeamWeightFact 2; # two beams elewhere
}



if {$frame==1 || $frame==4} {
if {$level==2} {
if {$pier==1 || $pier==4} {
set QdlCol $QdlCol5;
} elseif {$pier==2 || $pier==3} {
set QdlCol $QdlCol3;}
} elseif {$level==3} {
if {$pier==2 || $pier==3} {
set QdlCol $QdlCol4;
} elseif {$pier==4} {
set QdlCol $QdlCol6;}
} elseif {$level==4} {
if {$pier==3} {
set QdlCol $QdlCol5;
} elseif {$pier==4} {
set QdlCol $QdlCol6;}
}
} elseif {$frame==2 || $frame==3} {
if {$level==2} {
set QdlCol $QdlCol3;
} elseif {$level==3} {
set QdlCol $QdlCol4;
} elseif {$level==4} {
set QdlCol $QdlCol5;
}
}
if {$frame==1 || $frame==4} {
if {$level==2} {
set QGird $QGird5;
} elseif {$level==3} {
set QGird $QGird6;
} elseif {$level==4} {
set QGird $QGird7; }
} elseif {$frame==2 || $frame==3} {
if {$level==2} {
set QGird $QGird3;
} elseif {$level==3} {
set QGird $QGird4;
} elseif {$level==4} {
set QGird $QGird5;
}
}

if {$pier==1} {
if {$level==2} {
set QBeam $QBeam5;}
} elseif {$pier==2} {
if {$level==2} {
set QBeam $QBeam3;
} elseif {$level==3} {
set QBeam $QBeam4;}
} elseif {$pier==3} {
if {$level==2} {
set QBeam $QBeam3;
} elseif {$level==3} {
set QBeam $QBeam4;
} elseif {$level==4} {
set QBeam $QBeam5;}
} elseif {$pier==4} {
if {$level==2} {
set QBeam $QBeam5;
} elseif {$level==3} {
set QBeam $QBeam6;
} elseif {$level==4} {
set QBeam $QBeam7;
}
}

set GammaConcrete [expr 150*$pcf]; # Reinforced-Concrete floor slabs
set Tslab [expr 6*$in]; # 6-inch slab
set Lslab [expr 2*$LGird/2]; # slab extends a distance of $LGird/2 in/out of plane
set DLfactor 2.0; # scale dead load up a little
set Qslab [expr $GammaConcrete*$Tslab*$Lslab*$DLfactor];
set QdlBeam [expr $Qslab + $QBeam]; # dead load distributed along beam (one-way slab)
set QdlGird $QGird; # dead load distributed along girder

set WeightCol [expr $QdlCol*$LCol]; # total Column weight
set WeightBeam [expr $QdlBeam*$LBeam]; # total Beam weight
set WeightGird [expr $QdlGird*$LGird]; # total Beam weight



set WeightNode [expr $ColWeightFact*$WeightCol + $BeamWeightFact*$WeightBeam/2 + $GirdWeightFact*$WeightGird/2]
set MassNode [expr $WeightNode/$g];
set nodeID [expr $level*$Dlevel+$frame*$Dframe+$pier]
mass $nodeID [expr 1*$MassNode] 10e-5 [expr 1*$MassNode] 10e-5 10e-5 0; # define mass
set FloorWeight [expr $FloorWeight+$WeightNode];
}
}
lappend iFloorWeight $FloorWeight
set WeightTotal [expr $WeightTotal+ $FloorWeight]
set sumWiHi [expr $sumWiHi+$FloorWeight*($level-1)*$LCol]; # sum of storey weight times height, for lateral-load distribution
}



set MassTotal [expr $WeightTotal/$g]; # total mass

# --------------------------------------------------------------------------------------------------------------------------------
# LATERAL-LOAD distribution for static pushover analysis
# calculate distribution of lateral load based on mass/weight distributions along building height
# Fj = WjHj/sum(WiHi) * Weight at each floor j
set iFj ""
for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { ;
set FloorWeight [lindex $iFloorWeight [expr $level-1-1]];
set FloorHeight [expr ($level-1)*$LCol];
lappend iFj [expr $FloorWeight*$FloorHeight/$sumWiHi*$WeightTotal]; # per floor
}
set iNodePush $iMasterNode; # nodes for pushover/cyclic, vectorized
set iFPush $iFj; # lateral load for pushover, vectorized

# Define RECORDERS -------------------------------------------------------------
set SupportNodeFirst [lindex $iSupportNode 0]; # ID: first support node
set SupportNodeLast [lindex $iSupportNode 15]; # ID: last support node
set Masterid01 [expr 3*$Dframe+1*$Dlevel+3];
set Masterid11 [expr 3*$Dframe+2*$Dlevel+3];
set Masterid12 [expr 3*$Dframe+2*$Dlevel+3];
set Masterid22 [expr 3*$Dframe+3*$Dlevel+3];
set Masterid23 [expr 3*$Dframe+3*$Dlevel+3];
set Masterid33 [expr 3*$Dframe+4*$Dlevel+3];
set FreecorID1 [expr ($NFrame)*$Dframe+($NStory-1)*$Dlevel+($NBay+1)]; # ID: free node
set FreecorID2 [expr ($NFrame)*$Dframe+($NStory)*$Dlevel+($NBay+1)]; # ID: free node
set FreecorID3 [expr ($NFrame)*$Dframe+($NStory+1)*$Dlevel+($NBay+1)]; # ID: free node


recorder Drift -file $dataDir/Drnode01d1.out -time -iNode $Masterid01 -jNode $Masterid11 -dof 1 -perpDirn 2; # lateral drift
recorder Drift -file $dataDir/Drnode01d3.out -time -iNode $Masterid01 -jNode $Masterid11 -dof 3 -perpDirn 2; # lateral drift
recorder Drift -file $dataDir/Drnode12d1.out -time -iNode $Masterid12 -jNode $Masterid22 -dof 1 -perpDirn 2; # lateral drift
recorder Drift -file $dataDir/Drnode12d3.out -time -iNode $Masterid12 -jNode $Masterid22 -dof 3 -perpDirn 2; # lateral drift
recorder Drift -file $dataDir/Drnode23d1.out -time -iNode $Masterid23 -jNode $Masterid33 -dof 1 -perpDirn 2; # lateral drift
recorder Drift -file $dataDir/Drnode23d3.out -time -iNode $Masterid23 -jNode $Masterid33 -dof 3 -perpDirn 2; # lateral drift

recorder Node -file $dataDir/DFrees1.out -time -node $Masterid11 -dof 1 2 3 disp;
recorder Node -file $dataDir/DFrees2.out -time -node $Masterid22 -dof 1 2 3 disp;
recorder Node -file $dataDir/DFrees3.out -time -node $Masterid33 -dof 1 2 3 disp;
recorder Node -file $dataDir/Dcorner1.out -time -node $FreecorID1 -dof 1 2 3 disp;
recorder Node -file $dataDir/Dcorner2.out -time -node $FreecorID2 -dof 1 2 3 disp;
recorder Node -file $dataDir/Dcorner3.out -time -node $FreecorID3 -dof 1 2 3 disp;

recorder Node -file $dataDir/aFrees1.out -time -node $Masterid11 -dof 1 2 3 accel;
recorder Node -file $dataDir/aFrees2.out -time -node $Masterid22 -dof 1 2 3 accel;
recorder Node -file $dataDir/aFrees3.out -time -node $Masterid33 -dof 1 2 3 accel;

recorder Node -file $dataDir/RBase.out -time -nodeRange $SupportNodeFirst $SupportNodeLast -dof 1 2 3 reaction;
# Define DISPLAY -------------------------------------------------------------
DisplayModel3D DeformedShape ; # options: DeformedShape NodeNumbers ModeShape

# GRAVITY -------------------------------------------------------------
# define GRAVITY load applied to beams and columns -- eleLoad applies loads in local coordinate axis
# apply GRAVITY-- # apply gravity load, set it constant and reset time to zero, load pattern has already been defined
pattern Plain 101 Linear {
for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
for {set level 1} {$level <=$NStory} {incr level 1} {
if {$level==1} {
set a 1}
if {$level==2} {
set a 2}
if {$level==3} {
set a 3}

for {set pier $a} {$pier <= [expr $NBay+1]} {incr pier 1} {


if {$frame==1 || $frame==4} {
if {$level==1} {
if {$pier==1 || $pier==4} {
set QdlCol $QdlCol5;
} elseif {$pier==2 || $pier==3} {
set QdlCol $QdlCol3;}
} elseif {$level==2} {
if {$pier==2 || $pier==3} {
set QdlCol $QdlCol4;
} elseif {$pier==4} {
set QdlCol $QdlCol6;}
} elseif {$level==3} {
if {$pier==3} {
set QdlCol $QdlCol5;
} elseif {$pier==4} {
set QdlCol $QdlCol6;}
}
} elseif {$frame==2 || $frame==3} {
if {$level==1} {
set QdlCol $QdlCol3;
} elseif {$level==2} {
set QdlCol $QdlCol4;
} elseif {$level==3} {
set QdlCol $QdlCol5;
}
}

set elemID [expr $N0col+ $level*$Dlevel +$frame*$Dframe+$pier]
eleLoad -ele $elemID -type -beamUniform 0. 0. -$QdlCol; # COLUMNS }
}
}

for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} {
for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {
if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}

for {set bay $a} {$bay <= $NBay} {incr bay 1} {
set elemID [expr $N0beam + $level*$Dlevel +$frame*$Dframe+ $bay]
if {$pier==1} {
if {$level==2} {
set QBeam $QBeam5;}
} elseif {$pier==2} {
if {$level==2} {
set QBeam $QBeam3;
} elseif {$level==3} {
set QBeam $QBeam4;}
} elseif {$pier==3} {
if {$level==2} {
set QBeam $QBeam3;
} elseif {$level==3} {
set QBeam $QBeam4;
} elseif {$level==4} {
set QBeam $QBeam5;}
} elseif {$pier==4} {
if {$level==2} {
set QBeam $QBeam5;
} elseif {$level==3} {
set QBeam $QBeam6;
} elseif {$level==4} {
set QBeam $QBeam7;
}
}
eleLoad -ele $elemID -type -beamUniform -$QdlBeam 0.; # BEAMS

}
}
}
for {set frame 1} {$frame <=[expr $NFrame-1]} {incr frame 1} {
for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} {
if {$level==1} {
set a 1}
if {$level==2} {
set a 1}
if {$level==3} {
set a 2}
if {$level==4} {
set a 3}

for {set bay $a} {$bay <= $NBay+1} {incr bay 1} {

if {$frame==1 || $frame==4} {
if {$level==2} {
set QGird $QGird5;
} elseif {$level==3} {
set QGird $QGird6;
} elseif {$level==4} {
set QGird $QGird7; }
} elseif {$frame==2 || $frame==3} {
if {$level==2} {
set QGird $QGird3;
} elseif {$level==3} {
set QGird $QGird4;
} elseif {$level==4} {
set QGird $QGird5;
}
}

set elemID [expr $N0gird + $level*$Dlevel +$frame*$Dframe+ $bay]
eleLoad -ele $elemID -type -beamUniform -$QdlGird 0.; # GIRDS
}
}
}

}

puts goGravity
# Gravity-analysis parameters -- load-controlled static analysis
set Tol 1.0e-8; # convergence tolerance for test
variable constraintsTypeGravity Plain; # default;
if { [info exists RigidDiaphragm] == 1} {
if {$RigidDiaphragm=="ON"} {
variable constraintsTypeGravity Lagrange; # large model: try Transformation
}; # if rigid diaphragm is on
}; # if rigid diaphragm exists
constraints $constraintsTypeGravity ; # how it handles boundary conditions
numberer RCM; # 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 (large model: try UmfPack)
test EnergyIncr $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 NstepGravity 10; # apply gravity in 10 steps
set DGravity [expr 1./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
analysis Static; # define type of analysis static or transient
analyze $NstepGravity; # apply gravity
# ------------------------------------------------- maintain constant gravity loads and reset time to zero
loadConst -time 0.0

# -------------------------------------------------------------
puts "Model Built"

# --------------------------------------------------------------------------------------------------
# Example 8. 3D Uniform Earthquake Excitation
# Silvia Mazzoni & Frank McKenna, 2006
# execute this file after you have built the model, and after you apply gravity
#

# source in procedures
source ReadSMDfile.tcl; # procedure for reading GM file and converting it to proper format

# Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set iGMdirection "1 3"; # ground-motion direction
set iGMfile "G06090 G06000"; # ground-motion filenames
set iGMfact "4.32 4.32"; # ground-motion scaling factor

# Define DISPLAY -------------------------------------------------------------
# the deformed shape is defined in the build file
recorder plot $dataDir/Drnode23d1.out Displ-X 800 5 400 400 -columns 1 2 ; # a window to plot the nodal displacements versus time
recorder plot $dataDir/Drnode23d3.out Displ-Z 800 410 400 400 -columns 1 2 ; # a window to plot the nodal displacements versus time

# set up ground-motion-analysis parameters
set DtAnalysis [expr .005*$sec]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 12*$sec]; # maximum duration of ground-motion analysis -- should be 50*$sec

# ----------- set up analysis parameters
source LibAnalysisDynamicParameters.tcl; # constraintsHandler,DOFnumberer,system-ofequations,convergenceTest,solutionAlgorithm,integrator

# ------------ define & apply damping
# RAYLEIGH damping parameters, Where to put M/K-prop damping, switches (http://opensees.berkeley.edu/OpenSees/m ... l/1099.htm)
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
set xDamp 0.05; # damping ratio
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 1.0;
set KinitSwitch 0.0;
set nEigenI 1;
# mode 1
set nEigenJ 3; # mode 3
set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr $nEigenI-1]]; # eigenvalue mode i
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j
set omegaI [expr pow($lambdaI,0.5)];
set omegaJ [expr pow($lambdaJ,0.5)];
set alphaM [expr $MpropSwitch*$xDamp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]; # M-prop. damping; D = alphaM*M
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # current-K; +beatKcurr*KCurrent
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # initial-K; +beatKinit*Kini
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping

# --------------------------------- perform Dynamic Ground-Motion Analysis
# the following commands are unique to the Uniform Earthquake excitation
set IDloadTag 400; # for uniformSupport excitation
set dt $DtAnalysis;
# Uniform EXCITATION: acceleration input
foreach GMdirection $iGMdirection GMfile $iGMfile GMfact $iGMfact {
incr IDloadTag;
set inFile $GMdir/$GMfile.at2
set outFile $GMdir/$GMfile.g3; # set variable holding new filename (PEER files have .at2/dt2 extension)
ReadSMDFile $inFile $outFile dt; # call procedure to convert the ground-motion file
set GMfatt [expr $g*$GMfact]; # data in input file is in g Unifts -- ACCELERATION TH
set AccelSeries "Series -dt $dt -filePath $outFile -factor $GMfatt"; # time series information
pattern UniformExcitation $IDloadTag $GMdirection -accel $AccelSeries ; # create Unifform excitation
}
set Nsteps [expr int($TmaxAnalysis/$DtAnalysis)];
set ok [analyze $Nsteps $DtAnalysis]; # actually perform analysis; returns ok=0 if analysis was successful

if {$ok != 0} { # analysis was not successful.
# --------------------------------------------------------------------------------------------------
# change some analysis parameters to achieve convergence
# performance is slower inside this loop
# Time-controlled analysis
set ok 0;
set controlTime [getTime];
while {$controlTime < $TmaxAnalysis && $ok == 0} {
set controlTime [getTime]
set ok [analyze 1 $DtAnalysis]
if {$ok != 0} {
puts "Trying Newton with Initial Tangent .."
test NormDispIncr $Tol 1000 0
algorithm Newton -initial
set ok [analyze 1 $DtAnalysis]
test $testTypeDynamic $TolDynamic $maxNumIterDynamic 0
algorithm $algorithmTypeDynamic
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 $DtAnalysis]
algorithm $algorithmTypeDynamic
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch .8
set ok [analyze 1 $DtAnalysis]
algorithm $algorithmTypeDynamic
}
}
}; # end if ok !0

puts "Ground Motion Done. End Time: [getTime]"

silvia
Posts: 3912
Joined: Tue Jan 11, 2005 7:44 am
Location: Degenkolb Engineers
Contact:

Post by silvia » Thu Oct 09, 2008 6:30 am

check to see what's going on during the gravity analysis.
Silvia Mazzoni, PhD
Structural Consultant
Degenkolb Engineers
235 Montgomery Street, Suite 500
San Francisco, CA. 94104

Post Reply