i am seting up a 3D model with nonlinear fiber section and dispBeamColumn elements.
it`s one storey and i made it rigidDiaphraghm. but when i run eigen analysis, 2 first my eigenvalues are negative and furthermore i can`t set up rayleigh damping matrix?
what is wrong with my model?
model have 4 frame in each direction, it`s noticeable that Node 1 is mass center and assumed as master node
here is the code:
#****************************************************************SET UP**********************************************************************************
wipe
model BasicBuilder -ndm 3 -mdf 6
set dataDir Outputdata;
file mkdir $dataDir;
set GMdir "Groundmotion/";
source Units.tcl;
source DisplayPlane.tcl;
source DisplayModel3D.tcl;
source Boxsection.tcl;
source Beamsection.tcl;
#-------------------------------------------------------------Define geometry-----------------------------------------------------------------------
set LCol [expr 330*$cm];
set LBeamX [expr 600*$cm];
set LBeamZ [expr 500*$cm];
#----------calculated MODEL PARAMETERS, particular to this model
#----------Set up parameters that are particular to the model for displacement control
set IDctrlNode 1; # node where displacement is read for displacement control
set IDctrlDOF 1; # degree of freedom of displacement read for displacement control
#----------frame configuration
set NStory 1; # number of stories above ground level
set NBayX 1; # number of bays in X direction
set NBayZ 1; # number of bays in Z direction
puts "Number of Stories in Y: $NStory; Number of bays in X: $NBayX; Number of bays in Z: $NBayZ"
#-----------define NODAL COORDINATES
#-----------calculate locations of beam/column intersections:
set X1 0.;
set X2 [expr $X1+$LBeamX]
set X3 [expr $X2+$LBeamX]
set X4 [expr $X3+$LBeamX]
set Z1 0.;
set Z2 [expr $Z1+$LBeamZ]
set Z3 [expr $Z2+$LBeamZ]
set Z4 [expr $Z3+$LBeamZ]
set Y1 0
set Y2 $LCol
set HBuilding [expr $Y2]; # total building height
#-------------------------------------------------------------------Node-----------------------------------------------------------------------------
#----------Base
#----------frame 1
node 111 $X1 $Y1 $Z1
node 112 $X2 $Y1 $Z1
node 113 $X3 $Y1 $Z1
node 114 $X4 $Y1 $Z1
#----------frame 2
node 121 $X1 $Y1 $Z2
node 122 $X2 $Y1 $Z2
node 123 $X3 $Y1 $Z2
node 124 $X4 $Y1 $Z2
#----------frame 3
node 131 $X1 $Y1 $Z3
node 132 $X2 $Y1 $Z3
node 133 $X3 $Y1 $Z3
node 134 $X4 $Y1 $Z3
#----------frame 4
node 141 $X1 $Y1 $Z4
node 142 $X2 $Y1 $Z4
node 143 $X3 $Y1 $Z4
node 144 $X4 $Y1 $Z4
#----------first storey
#----------frame 1
node 211 $X1 $Y2 $Z1
node 212 $X2 $Y2 $Z1
node 213 $X3 $Y2 $Z1
node 214 $X4 $Y2 $Z1
#----------frame 2
node 221 $X1 $Y2 $Z2
node 222 $X2 $Y2 $Z2
node 223 $X3 $Y2 $Z2
node 224 $X4 $Y2 $Z2
#----------frame 3
node 231 $X1 $Y2 $Z3
node 232 $X2 $Y2 $Z3
node 233 $X3 $Y2 $Z3
node 234 $X4 $Y2 $Z3
#----------frame 4
node 241 $X1 $Y2 $Z4
node 242 $X2 $Y2 $Z4
node 243 $X3 $Y2 $Z4
node 244 $X4 $Y2 $Z4
#-----------------------------------------------------------------------*****----------------------------------------------------------------
#--------------------------------------------------------------define Rigid Floor Diaphragm------------------------------------------------------------
set RigidDiaphragm ON; # options: ON, OFF
set Xa [expr ($X1+$X4)/2]; # mid-span coordinate for rigid diaphragm
set Za [expr ($Z4+$Z1)/2];
#-------rigid-diaphragm node in center of each diaphram
node 1 $Xa $Y2 $Za; # master nodes for rigid diaphragm
#-------Constraints for rigid diaphragm master node
fix 1 0 1 0 1 0 1 # 0 means constrainted
# ------define Rigid Diaphram, dof 2 is normal to floor
set perpDirn 2; # direction Y is prependecular to ZX plane
rigidDiaphragm $perpDirn 1 211 212 213 214 221 222 223 224 231 232 233 234 241 242 243 244;
#-------determine support nodes where ground motions are input, for multiple-support excitation
set iSupportNode "111 112 113 114 121 122 123 124 131 132 133 134 141 142 143 144"
# BOUNDARY CONDITIONS
fixY 0.0 1 1 1 1 1 1; # restrained all Y=0.0 nodes
#----------------------------------------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------Define SECTIONS -------------------------------------------------------------
set SectionType FiberSection; # options: Elastic;FiberSection
#------------define section tags:
set ColSecTag 1
set BeamSecTag 2
set ColSecTagFiber 3
set BeamSecTagFiber 4
set SecTagTorsion 70
set matIDhard 1;
#-----------define MATERIAL properties
set Fy [expr 2400.*$kgBcm2]
set Es [expr 2.1*pow(10,6)*$kgBcm2]; # Steel Young's Modulus
set nu 0.3;
set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus
set b 0.02
set matIDsteel01 1
uniaxialMaterial Steel01 $matIDsteel01 $Fy $Es $b
#-------------ELEMENT properties
#-------------Structural-Steel I-section properties
# Beam sections: PL200X8+2PL200X12
set d [expr 20*$cm+(2*(1.2*$cm))]; # depth
set bf [expr 20*$cm]; # flange width
set tf [expr 1.2*$cm]; # flange thickness
set tw [expr 0.8*$cm]; # 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
Beamsection $BeamSecTagFiber $matIDsteel01 $d $bf $tf $tw $nfdw $nftw $nfbf $nftf
#-------------Structural-Steel Box-section properties
# Column sections:BOX220X220X12X12
set a [expr 22*$cm]
set t [expr 1.2*$cm]
set nfdw 16;
set nftw 2;
set nfbf 16;
set nftf 4;
Boxsection $ColSecTagFiber $matIDsteel01 $a $t $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
#--------------------------------------------------------------------------******--------------------------------------------------------------------------------
set ABeam 64*$cm2
set ACol 99.84*$cm2
set QCol [expr 78.3*($kg/$m)]; # BOXsection weight per length
set QBeam [expr 52*($kg/$m)]; # Isection weight per length
#-----------------------------------------------------------------------define ELEMENTS -----------------------------------------------------------------------------
#---------set up geometric transformations of element
#---------separate columns and beams, in case of P-Delta analysis for columns
#---------in 3D model, assign vector vecxz
set IDColTransf 1; # all columns
set IDBeamTransf 2; # all beams
set ColTransfType PDelta ; # options, Linear PDelta Corotational
geomTransf $ColTransfType $IDColTransf 1 0 0 ; # only columns can have PDelta effects (gravity effects)
geomTransf Linear $IDBeamTransf 0 1 0
#----------Define Beam-Column Elements
set np 5; # number of Gauss integration points for nonlinear curvature distribution
#----------Columns
element dispBeamColumn 11 111 211 $np $ColSecTag $IDColTransf;
element dispBeamColumn 12 112 212 $np $ColSecTag $IDColTransf;
element dispBeamColumn 13 113 213 $np $ColSecTag $IDColTransf;
element dispBeamColumn 14 114 214 $np $ColSecTag $IDColTransf;
element dispBeamColumn 21 121 221 $np $ColSecTag $IDColTransf;
element dispBeamColumn 22 122 222 $np $ColSecTag $IDColTransf;
element dispBeamColumn 23 123 223 $np $ColSecTag $IDColTransf;
element dispBeamColumn 24 124 224 $np $ColSecTag $IDColTransf;
element dispBeamColumn 31 131 231 $np $ColSecTag $IDColTransf;
element dispBeamColumn 32 132 232 $np $ColSecTag $IDColTransf;
element dispBeamColumn 33 133 233 $np $ColSecTag $IDColTransf;
element dispBeamColumn 34 134 234 $np $ColSecTag $IDColTransf;
element dispBeamColumn 41 141 241 $np $ColSecTag $IDColTransf;
element dispBeamColumn 42 142 242 $np $ColSecTag $IDColTransf;
element dispBeamColumn 43 143 243 $np $ColSecTag $IDColTransf;
element dispBeamColumn 44 144 244 $np $ColSecTag $IDColTransf;
#----------Beams
element dispBeamColumn 1112 211 212 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1213 212 213 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1314 213 214 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2122 221 222 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2223 222 223 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2324 223 224 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3132 231 232 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3233 232 233 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3334 233 234 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 4142 241 242 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 4243 242 243 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 4344 243 244 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1121 211 221 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2131 221 231 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3141 231 241 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1222 212 222 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2232 222 232 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3242 232 242 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1323 213 223 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2333 223 233 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3343 233 243 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 1424 214 224 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 2434 224 234 $np $BeamSecTag $IDBeamTransf;
element dispBeamColumn 3444 234 244 $np $BeamSecTag $IDBeamTransf;
#-----------------------------------------------------------------------------******--------------------------------------------------------------
#--------------------------------------------------------------Define GRAVITY LOADS, weight and mass------------------------------------------------
#------------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
set QDead [expr 650*$kgBm2];
set Qlive [expr 200*$kgBm2]
set LBeam1X [expr $LBeamX/2]
set LBeam1Z [expr $LBeamZ/2]
set DLfactor 1.0; # scale dead load up a little
set QslabZ [expr ($QDead+0.2*$Qlive)*$LBeam1X*$DLfactor];
set QslabX [expr ($QDead+0.2*$Qlive)*$LBeam1Z*$DLfactor];
set QdlBeamX [expr $QslabX + $QBeam];
set QdlBeamZ [expr $QslabZ + $QBeam]; # dead load distributed along beam (one-way slab)
set MassCol [expr $QCol*$LCol]; # total Column weight
puts "hello world"
set MassBeamX [expr $QdlBeamX*$LBeamX]; # total Beam weight
set MassBeamZ [expr $QdlBeamZ*$LBeamZ];
#------------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 M1 [expr ($MassCol/2 +$MassBeamX/2)]; # 211 214 241 244
set M2 [expr ($MassCol/2 +$MassBeamX/2+$MassBeamZ/2)]; #221 231 212 213 241 243 224 234
set M3 [expr ($MassCol/2 +$MassBeamX+$MassBeamZ)]; # 222 232 233 223
set b 1e-9;
mass 211 $M1 $b $M1 $b $b $b;
mass 214 $M1 $b $M1 $b $b $b;
mass 241 $M1 $b $M1 $b $b $b;
mass 244 $M1 $b $M1 $b $b $b;
mass 221 $M2 $b $M2 $b $b $b;
mass 231 $M2 $b $M2 $b $b $b;
mass 212 $M2 $b $M2 $b $b $b;
mass 213 $M2 $b $M2 $b $b $b;
mass 241 $M2 $b $M2 $b $b $b;
mass 243 $M2 $b $M2 $b $b $b;
mass 224 $M2 $b $M2 $b $b $b;
mass 234 $M2 $b $M2 $b $b $b;
mass 222 $M3 $b $M3 $b $b $b;
mass 232 $M3 $b $M3 $b $b $b;
mass 233 $M3 $b $M3 $b $b $b;
mass 223 $M3 $b $M3 $b $b $b;
set MassTotal [expr (16*$MassCol)/2 + 12*$MassBeamX + 12*$MassBeamZ];
set WeightTotal [expr ((16*$MassCol)/2 + 12*$MassBeamX + 12*$MassBeamZ)*$g];
#------------------------------------------------------------------------------******---------------------------------------------------------
#----------------------------------------------------------------------------Define DISPLAY -------------------------------------------------------------
set xPixels 1200; # height of graphical window in pixels
set yPixels 800; # height of graphical window in pixels
set xLoc1 10; # horizontal location of graphical window (0=upper left-most corner)
set yLoc1 10; # vertical location of graphical window (0=upper left-most corner)
set dAmp 2; # scaling factor for viewing deformed shape, it depends on the dimensions of the model
DisplayModel3D NodeNumbers $dAmp $xLoc1 $yLoc1 $xPixels $yPixels
#-------------------------------------------------------------------------------******-------------------------------------------------------------
#----------------------------------------------------------------------------Define RECORDERS -------------------------------------------------------------
recorder Node -file $dataDir/DFree.out -time -node 1 -dof 1 disp; # displacements of free node
recorder Node -file $dataDir/edgePoints.out -time -node 211 214 241 244 -dof 1 disp; # displacements of edge nodes
recorder Node -file $dataDir/RBase.out -time -node 111 -dof 1 2 3 reaction; # support reaction
recorder Node -file $dataDir/eigen1.out -time -node 1 -dof 1 2 3 4 5 6 eigen 1;
recorder Node -file $dataDir/eigen2.out -time -node 1 -dof 1 2 3 4 5 6 eigen 2;
recorder Node -file $dataDir/eigen3.out -time -node 1 -dof 1 2 3 4 5 6 eigen 3;
recorder Node -file $dataDir/eigen4.out -time -node 1 -dof 1 2 3 4 5 6 eigen 4;
recorder Node -file $dataDir/eigen5.out -time -node 1 -dof 1 2 3 4 5 6 eigen 5;
recorder Node -file $dataDir/eigen6.out -time -node 1 -dof 1 2 3 4 5 6 eigen 6;
recorder Element -file $dataDir/Fel1.out -time -ele 11 localForce; # element forces in local coordinates
set yFiber [expr 0.]; # fiber location for stress-strain recorder, local coords
set zFiber [expr 0.]; # fiber location for stress-strain recorder, local coords
recorder Element -file $dataDir/SSreinfEle1sec1.out -time -ele 22 32 33 23 section $np fiber $yFiber $zFiber stressStrain; # steel fiber stress-strain, node i
eigen frequency 3
#-------------------------------------------------------------------------------******-----------------------------------------------------------------------------
#---------------------------------------------------------------------------define GRAVITY -------------------------------------------------------------
#-----------------GRAVITY LOADS
#-----------------define gravity load applied to beams and columns -- eleLoad applies loads in local coordinate axis
pattern Plain 1 Linear {
#------Column
eleLoad -ele 11 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 12 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 13 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 14 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 21 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 22 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 23 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 24 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 31 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 32 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 33 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 34 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 41 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 42 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 43 -type -beamUniform 0. 0. -[expr $QCol*$g];
eleLoad -ele 44 -type -beamUniform 0. 0. -[expr $QCol*$g];
#-----BeamX
eleLoad -ele 1112 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 1213 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 1314 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 2122 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 2223 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 2324 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 3132 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 3233 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 3334 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 4142 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 4243 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 4344 -type -beamUniform -[expr $QdlBeamX*$g] 0.;
eleLoad -ele 1121 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 2131 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 3141 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 1222 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 2232 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 3242 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 1323 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 2333 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 3343 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 1424 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 2434 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
eleLoad -ele 3444 -type -beamUniform -[expr $QdlBeamZ*$g] 0.;
}
#----------------------------------------------------------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 2 ; # determine if convergence has been achieved at the end of an iteration step
algorithm Linear ; # 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
set Tol 1.0e-6;
puts "Model Built"
negative eigenvalue?
Moderators: silvia, selimgunay, Moderators
Re: negative eigenvalue?
I run source eigen.tcl in this file
# ---------------------------------------------------------------Transient analysis----------------------------------------------------------------
# one storey symmetric model. Uniform Earthquake Excitation
# execute this file after you have built the model, and after you apply gravity
#
#--------source in procedures
source model.tcl
source ReadSMDfile.tcl; # procedure for reading GM file and converting it to proper format
puts "stage 2"
#---------Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set GMdirection 1; # ground-motion excitation in X direction
set GMfile "taft" ; # ground-motion filenames
set GMfact 1.5; # ground-motion scaling factor
#---------------------------------------------------------------------Define DISPLAY -------------------------------------------------------------
set xPixels 1200; # height of graphical window in pixels
set yPixels 800; # height of graphical window in pixels
set xLoc1 10; # horizontal location of graphical window (0=upper left-most corner)
set yLoc1 10; # vertical location of graphical window (0=upper left-most corner)
set ViewScale 10; # scaling factor for viewing deformed shape, it depends on the dimensions of the model
DisplayModel3D DeformedShape $ViewScale $xLoc1 $yLoc1 $xPixels $yPixels
recorder plot $dataDir/DFree.out Disp1DOF$GMdirection 1200 10 400 400 -columns 1 2; # a window to plot the nodal displacements versus time
puts "stage 3"
#------------set up ground-motion-analysis parameters
set DtAnalysis [expr 0.005*$sec]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 20. *$sec]; # maximum duration of ground-motion analysis -- should be 50*$sec
#------------set up analysis parameters
source LibAnalysisDynamicParameters.tcl;
source eigen.tcl
# --------------------------------------------------------------- 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.02; # 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]];
set lambdaI [lindex $lambdaN [expr $nEigenI-1]];
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]];
puts "lambdaI is $lambdaI"
puts "lamdaJ=$lambdaJ"
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
#---------Uniform EXCITATION: acceleration input
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";
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
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]"
# ---------------------------------------------------------------Transient analysis----------------------------------------------------------------
# one storey symmetric model. Uniform Earthquake Excitation
# execute this file after you have built the model, and after you apply gravity
#
#--------source in procedures
source model.tcl
source ReadSMDfile.tcl; # procedure for reading GM file and converting it to proper format
puts "stage 2"
#---------Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set GMdirection 1; # ground-motion excitation in X direction
set GMfile "taft" ; # ground-motion filenames
set GMfact 1.5; # ground-motion scaling factor
#---------------------------------------------------------------------Define DISPLAY -------------------------------------------------------------
set xPixels 1200; # height of graphical window in pixels
set yPixels 800; # height of graphical window in pixels
set xLoc1 10; # horizontal location of graphical window (0=upper left-most corner)
set yLoc1 10; # vertical location of graphical window (0=upper left-most corner)
set ViewScale 10; # scaling factor for viewing deformed shape, it depends on the dimensions of the model
DisplayModel3D DeformedShape $ViewScale $xLoc1 $yLoc1 $xPixels $yPixels
recorder plot $dataDir/DFree.out Disp1DOF$GMdirection 1200 10 400 400 -columns 1 2; # a window to plot the nodal displacements versus time
puts "stage 3"
#------------set up ground-motion-analysis parameters
set DtAnalysis [expr 0.005*$sec]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 20. *$sec]; # maximum duration of ground-motion analysis -- should be 50*$sec
#------------set up analysis parameters
source LibAnalysisDynamicParameters.tcl;
source eigen.tcl
# --------------------------------------------------------------- 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.02; # 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]];
set lambdaI [lindex $lambdaN [expr $nEigenI-1]];
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]];
puts "lambdaI is $lambdaI"
puts "lamdaJ=$lambdaJ"
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
#---------Uniform EXCITATION: acceleration input
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";
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
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]"