Elastic Response Spectra Script

silvia · Post by **silvia** » Thu May 12, 2005 12:25 pm

All,
I have created a script that calculates elastic response spectra for a series of ground motions. This script will eventually be incorporated into the examples manual. If anyone is interested in this script I am placing it here. If you have any questions, comments or recommendations, I'd love to hear them. As usual with my input files, it consists of a series of files:

main: ExSpectraTZ.tcl, which calls the following files:
LibGMfiles.tcl
LibReadSMDFile.tcl
LibSDOFanalysis.tcl
LibUnits.tcl

the files in LibGMfiles.tcl are files I have, you should use files in your own bin and put them in a subdirectory called GMfiles

I am also attaching a script for post-processing in matlab:
spectra.m

please make sure you set up the proper directories for post-processing. Plus, change the ground-motions files to ground-motions you already have (do this in both tcl script and matlab post-processing).

here are the files:
ExSpectraTZ.tcl

# -------------------------------------------------------------------------
# ExSpectraTZ.tcl
#
# Silvia Mazzoni, May 2005
#
# this script is used to calcuate elastic response spectra for a series of
# ground motions over a specified range of dampings and natural periods.

file mkdir spettra

source LibUnits.tcl; # define units and constants
source LibGMfiles.tcl; # read names for ground-motion files
source LibSDOFanalysis.tcl; # define elastic dynamic time-history analysis

set GMfact 1.0; # Ground-motion scaling factor

set iTs "0.1 0.25 0.35 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0"
set iZs "0.01 0.02 0.04 0.05 0.06 0.08 0.1 0.15 0.20 0.25 "

# output periods
set OUTfile [open spettra/Periods.out w];
foreach T $iTs {
puts $OUTfile $T;
}; # end of period
close $OUTfile;

# output dampings
set OUTfile [open spettra/Dampings.out w];
foreach Z $iZs {
puts $OUTfile $Z;
}; # end
close $OUTfile;

foreach GroundFile $iGroundFile {
puts Spectrum_ANALYSIS........GroundMotion$GroundFile
set SdOUTfile [open spettra/SdSpett$GroundFile.out w];
set SvOUTfile [open spettra/SvSpett$GroundFile.out w];
set SaOUTfile [open spettra/SaSpett$GroundFile.out w];
foreach T $iTs {
set omega [expr 2*$PI/$T];
set iDi "$T "
set iVi "$T "
set iAi "$T "
foreach Z $iZs {
SDOFanalysis $omega $Z $GroundFile $GMfact; # elastic dynamic SDOF: Dmax

set Dcontent [read [open tmpEnvNodalDisp.out ]];
set Di [lindex $Dcontent 2];
puts Sd$Di
set iDi "$iDi $Di";

set Vcontent [read [open tmpEnvNodalVel.out ]];
set Vi [lindex $Vcontent 2];
set iVi "$iVi $Vi";

set Acontent [read [open tmpEnvNodalAccel.out ]];
set Ai [lindex $Acontent 2];
set iAi "$iAi $Ai";
}
puts $SdOUTfile $iDi; # output converged Dmax for this ground motion
puts $SvOUTfile $iVi; # output converged Vmax for this ground motion
puts $SaOUTfile $iAi; # output converged Amax for this ground motion
}; # end of period
close $SdOUTfile;
close $SvOUTfile;
close $SaOUTfile;
}; # end of $GroundFile

LibGMfiles.tcl

# GMFILES.tcl
# ------------------------------------------------------------------------------------------------------------------------- set up GROUND MOTIONS
# # 7 August 2001
# labels for Ground-Motion Files

#small magnitude, small distance
set iGroundFileSMSR {"IV79cal" "IV79chi" "IV79e01" "IV79e12" "IV79e13" "IV79qkp" "IV79wsm" "LV80kod" "LV80srm" "MH84agw" "MH84g02" "MH84g03" "MH84gmr" "PM73phn" "PS86psa" "WN87cas" "WN87cat" "WN87flo" "WN87w70" "WN87wat" }
#small magnitude, large distance
set iGroundFileSMLR {"BO42elc" "CO83c05" "CO83c08" "IV79cc4" "IV79cmp" "IV79dlt" "IV79nil" "IV79pls" "IV79vct" "LV80stp" "MH84cap" "MH84hch" "MH84sjb" "PS86h06" "PS86ino" "WN87bir" "WN87cts" "WN87har" "WN87sse" "WN87stc" }

#large magnitude, small distance
set iGroundFileLMSR {"IV40elc" "LD92yer" "LP89agw" "LP89cap" "LP89g03" "LP89g04" "LP89gmr" "LP89hch" "LP89hda" "LP89svl" "NR94cnp" "NR94far" "NR94fle" "NR94glp" "NR94hol" "NR94stc" "SF71pel" "SH87bra" "SH87icc" "SH87wsm" }
#large magnitude, large distance
set iGroundFileLMLR {"BM68elc" "LD92ind" "LD92psa" "LP89a2e" "LP89fms" "LP89hvr" "LP89sjw" "LP89slc" "NR94bad" "NR94cas" "NR94cen" "NR94jab" "NR94lh1" "NR94loa" "NR94lv2" "NR94php" "NR94pic" "NR94sor" "NR94sse" "NR94ver" }
set iGroundFileLMLR {"SH87wsm" "BM68elc" "LD92ind" "LD92psa" "LP89a2e" "LP89fms" "LP89hvr" "LP89sjw" "LP89slc" "NR94bad" "NR94cas" "NR94cen" "NR94jab" "NR94lh1" "NR94loa" "NR94lv2" "NR94php" "NR94pic" "NR94sor" "NR94sse" "NR94ver" }
# concatenate them all
set iGroundFile $iGroundFileSMSR$iGroundFileSMLR$iGroundFileLMSR$iGroundFileLMLR; #*****************

# ground-motion information
#PEER GROUND MOTIONS, obtained from Eric Anderson. 7 august 2001
#small magnitude, small distance
#Record ID Event Year Magnitude Station R (km) NEHRP Soil Mechanism fHP (Hz) fLP (Hz) PGA (g) PGV (cm/s) PGD (cm) Dsm (s) Rec. Length (s)
#IV79cal Imperial Valley 1979 6.5 Calipatria Fire Station 23.8 D strike-slip 0.10 40.0 0.078 13.3 6.2 23.3 39.5
#IV79chi Imperial Valley 1979 6.5 Chihuahua 28.7 D strike-slip 0.05 0.270 24.9 9.1 20.1 40.0
#IV79e01 Imperial Valley 1979 6.5 El Centro Array #1 15.5 D strike-slip 0.10 40.0 0.139 16.0 10.0 8.9 39.5
#IV79e12 Imperial Valley 1979 6.5 El Centro Array #12 18.2 D strike-slip 0.10 40.0 0.116 21.8 12.1 19.4 39.0
#IV79e13 Imperial Valley 1979 6.5 El Centro Array #13 21.9 D strike-slip 0.20 40.0 0.139 13.0 5.8 21.2 39.5
#IV79qkp Imperial Valley 1979 6.5 Cucapah 23.6 D strike-slip 0.05 0.309 36.3 10.4 15.7 40.0
#IV79wsm Imperial Valley 1979 6.5 Westmorland Fire Station 15.1 D strike-slip 0.10 40.0 0.110 21.9 10.0 25.2 40.0
#LV80kod Livermore 1980 5.8 San Ramon Fire Station 21.7 D strike-slip 0.20 15.0 0.040 4.0 1.2 14.2 21.0
#LV80srm Livermore 1980 5.8 San Ramon - Eastman Kodak 17.6 D strike-slip 0.20 20.0 0.076 6.1 1.7 27.3 40.0
#MH84agw Morgan Hill 1984 6.2 Agnews State Hospital 29.4 D strike-slip 0.20 13.0 0.032 5.5 2.1 40.3 59.9
#MH84g02 Morgan Hill 1984 6.2 Gilroy Array #2 15.1 D strike-slip 0.10 29.0 0.069 2.9 1.3 16.4 30.0
#MH84g03 Morgan Hill 1984 6.2 Gilroy Array #3 14.6 D strike-slip 0.10 37.0 0.194 11.2 2.4 16.0 40.0
#MH84gmr Morgan Hill 1984 6.2 Gilroy Array #7 14.0 D strike-slip 0.10 30.0 0.113 6.0 1.8 10.7 30.0
#PM73phn Point Mugu 1973 5.8 Port Hueneme 25.0 D reverse-slip 0.20 25.0 0.112 14.8 2.6 10.8 23.2
#PS86psa N. Palm Springs 1986 6.0 Palm Springs Airport 16.6 D strike-slip 0.20 60.0 0.187 12.2 2.1 15.6 30.0
#WN87cas Whittier Narrows 1987 6.0 Compton - Castlegate St # 16.9 D reverse 0.09 25.0 0.332 27.1 5.0 8.0 31.2
#WN87cat Whittier Narrows 1987 6.0 Carson - Catskill Ave # 28.1 D reverse 0.18 25.0 0.042 3.8 0.8 20.6 32.9
#WN87flo Whittier Narrows 1987 6.0 Brea - S Flower Av # 17.9 D reverse 0.16 25.0 0.115 7.1 1.2 9.4 27.6
#WN87w70 Whittier Narrows 1987 6.0 LA - W 70th St # 16.3 D reverse 0.20 25.0 0.151 8.7 1.5 11.2 31.9
#WN87wat Whittier Narrows 1987 6.0 Carson - Water St # 24.5 D reverse 0.20 25.0 0.104 9.0 1.9 15.2 29.7
#small magnitude, large distance
#Record ID Event Year Magnitude Station R (km) NEHRP Soil Mechanism fHP (Hz) fLP (Hz) PGA (g) PGV (cm/s) PGD (cm) Dsm (s) Rec. Length (s)
#BO42elc Borrego 1942 6.5 El Centro Array #9 49.0 D 0.10 15.0 0.068 3.9 1.4 29.5 40.0
#CO83c05 Coalinga 1983 6.4 Parkfield - Cholame 5W 47.3 D reverse-oblique 0.20 22.0 0.131 10.0 1.3 14.6 40.0
#CO83c08 Coalinga 1983 6.4 Parkfield - Cholame 8W 50.7 D reverse-oblique 0.20 23.0 0.098 8.6 1.5 15.0 32.0
#IV79cc4 Imperial Valley 1979 6.5 Coachella Canal #4 49.3 D strike-slip 0.20 40.0 0.128 15.6 3.0 10.0 28.5
#IV79cmp Imperial Valley 1979 6.5 Compuertas 32.6 D strike-slip 0.20 0.186 13.9 2.9 21.7 36.0
#IV79dlt Imperial Valley 1979 6.5 Delta 43.6 D strike-slip 0.05 0.238 26.0 12.1 51.1 99.9
#IV79nil Imperial Valley 1979 6.5 Niland Fire Station 35.9 D strike-slip 0.10 30.0 0.109 11.9 6.9 21.7 40.0
#IV79pls Imperial Valley 1979 6.5 Plaster City 31.7 D strike-slip 0.10 40.0 0.057 5.4 1.9 10.7 18.7
#IV79vct Imperial Valley 1979 6.5 Victoria 54.1 D strike-slip 0.20 0.167 8.3 1.1 17.1 40.0
#LV80stp Livermore 1980 5.8 Tracy - Sewage Treatment Plant 37.3 D strike-slip 0.08 15.0 0.073 7.6 1.8 20.2 33.0
#MH84cap Morgan Hill 1984 6.2 Capitola 38.1 D strike-slip 0.20 30.0 0.099 4.9 0.6 17.2 36.0
#MH84hch Morgan Hill 1984 6.2 Hollister City Hall 32.5 D strike-slip 0.20 19.0 0.071 7.4 1.6 21.4 28.3
#MH84sjb Morgan Hill 1984 6.2 San Juan Bautista 30.3 C strike-slip 0.10 21.0 0.036 4.4 1.5 19.0 28.0
#PS86h06 N. Palm Springs 1986 6.0 San Jacinto Valley Cemetery 39.6 D strike-slip 0.20 31.0 0.063 4.4 1.2 17.0 40.0
#PS86ino N. Palm Springs 1986 6.0 Indio 39.6 D strike-slip 0.10 35.0 0.064 6.6 2.2 18.6 30.0
#WN87bir Whittier Narrows 1987 6.0 Downey - Birchdale # 56.8 D reverse 0.15 25.0 0.299 37.8 5.0 3.8 28.6
#WN87cts Whittier Narrows 1987 6.0 LA - Century City CC South 31.3 D reverse 0.20 25.0 0.051 3.5 0.6 20.2 40.0
#WN87har Whittier Narrows 1987 6.0 LB - Harbor Admin FF 34.2 D reverse 0.25 25.0 0.071 7.3 0.9 24.9 40.0
#WN87sse Whittier Narrows 1987 6.0 Terminal Island - S Seaside # 35.7 D reverse 0.20 25.0 0.042 3.9 1.0 16.3 22.9
#WN87stc Whittier Narrows 1987 6.0 Northridge - Saticoy St # 39.8 D reverse 0.20 25.0 0.118 5.1 0.8 19.8 40.0
#large magnitude, small distance
#Record ID Event Year Magnitude Station R (km) NEHRP Soil Mechanism fHP (Hz) fLP (Hz) PGA (g) PGV (cm/s) PGD (cm) Dsm (s) Rec. Length (s)
#IV40elc Imperial Valley 1940 7.0 El Centro Array #9 12.0 D strike-slip 0.20 15.0 0.215 30.2 23.9 24.0 40.0
#LD92yer Landers 1992 7.3 Yermo Fire Station # 24.9 D strike-slip 0.07 23.0 0.245 51.5 43.8 17.6 44.0
#LP89agw Loma Prieta 1989 6.9 Agnews State Hospital 28.2 D reverse-oblique 0.20 30.0 0.172 26.0 12.6 18.4 40.0
#LP89cap Loma Prieta 1989 6.9 Capitola 14.5 D reverse-oblique 0.20 40.0 0.443 29.3 5.5 13.2 40.0
#LP89g03 Loma Prieta 1989 6.9 Gilroy Array #3 14.4 D reverse-oblique 0.10 40.0 0.367 44.7 19.3 11.4 39.9
#LP89g04 Loma Prieta 1989 6.9 Gilroy Array #4 16.1 D reverse-oblique 0.20 30.0 0.212 37.9 10.1 14.8 40.0
#LP89gmr Loma Prieta 1989 6.9 Gilroy Array #7 24.2 D reverse-oblique 0.20 40.0 0.226 16.4 2.5 11.5 40.0
#LP89hch Loma Prieta 1989 6.9 Hollister City Hall 28.2 D reverse-oblique 0.10 29.0 0.247 38.5 17.8 17.4 39.1
#LP89hda Loma Prieta 1989 6.9 Hollister Differential Array 25.8 D reverse-oblique 0.10 33.0 0.279 35.6 13.1 13.2 39.6
#LP89svl Loma Prieta 1989 6.9 Sunnyvale - Colton Ave. 28.8 D reverse-oblique 0.10 40.0 0.207 37.3 19.1 21.2 39.3
#NR94cnp Northridge 1994 6.7 Canoga Park - Topanga Can 15.8 D reverse-slip 0.05 30.0 0.420 60.8 20.2 10.4 25.0
#NR94far Northridge 1994 6.7 LA - N Faring Rd 23.9 D reverse-slip 0.13 30.0 0.273 15.8 3.3 8.8 30.0
#NR94fle Northridge 1994 6.7 LA - Fletcher Dr 29.5 D reverse-slip 0.15 30.0 0.240 26.2 3.6 11.8 30.0
#NR94glp Northridge 1994 6.7 Glendale - Las Palmas 25.4 D reverse-slip 0.10 30.0 0.206 7.4 1.8 11.5 30.0
#NR94hol Northridge 1994 6.7 LA - Holywood Stor FF # 25.5 D reverse-slip 0.20 23.0 0.231 18.3 4.8 12.0 40.0
#NR94stc Northridge 1994 6.7 Northridge - 17645 Saticoy St 13.3 D reverse-slip 0.10 30.0 0.368 28.9 8.4 15.7 30.0
#SF71pel San Fernando 1971 6.6 LA - Hollywood Stor Lot 21.2 D reverse-slip 0.20 35.0 0.174 14.9 6.3 11.2 28.0
#SH87bra Superstition Hills 1987 6.7 Brawley 18.2 D strike-slip 0.10 23.0 0.156 13.9 5.4 13.5 22.1
#SH87icc Superstition Hills 1987 6.7 El Centro Imp. Co. Cent 13.9 D strike-slip 0.10 40.0 0.358 46.4 17.5 16.1 40.0
#SH87wsm Superstition Hills 1987 6.7 Westmorland Fire Station 13.3 D strike-slip 0.10 35.0 0.172 23.5 13.0 19.6 40.0
#large magnitude, large distance
#Record ID Event Year Magnitude Station R (km) NEHRP Soil Mechanism fHP (Hz) fLP (Hz) PGA (g) PGV (cm/s) PGD (cm) Dsm (s) Rec. Length (s)
#BM68elc Borrego Mountain 1968 6.8 El Centro Array #9 46.0 D strike-slip 0.20 12.8 0.057 13.2 10.0 28.7 40.0
#LD92ind Landers 1992 7.3 Indio - Coachella Canal # 55.7 D strike-slip 0.10 23.0 0.109 15.2 9.7 37.9 60.0
#LD92psa Landers 1992 7.3 Palm Springs Airport # 37.5 D strike-slip 0.07 23.0 0.089 13.8 5.3 35.8 60.0
#LP89a2e Loma Prieta 1989 6.9 APEEL 2E Hayward Muir Sch 57.4 D reverse-oblique 0.20 30.0 0.171 13.7 3.9 12.8 40.0
#LP89fms Loma Prieta 1989 6.9 Freemont - Emerson Court 42.4 D reverse-oblique 0.10 32.0 0.141 12.9 8.4 17.9 39.7
#LP89hvr Loma Prieta 1989 6.9 Halls Valley 31.6 D reverse-oblique 0.20 22.0 0.134 15.4 3.3 16.2 40.0
#LP89sjw Loma Prieta 1989 6.9 Salinas - John & Work 32.6 D reverse-oblique 0.10 28.0 0.112 15.7 7.9 20.3 40.0
#LP89slc Loma Prieta 1989 6.9 Palo Alto - SLAC Lab 36.3 D reverse-oblique 0.20 33.0 0.194 37.5 10.0 12.5 39.6
#NR94bad Northridge 1994 6.7 Covina - W. Badillo 56.1 D reverse-slip 0.20 30.0 0.100 5.8 1.2 17.4 35.0
#NR94cas Northridge 1994 6.7 Compton - Castlegate St 49.6 D reverse-slip 0.20 30.0 0.136 7.1 2.2 23.4 39.8
#NR94cen Northridge 1994 6.7 LA - Centinela St 30.9 D reverse-slip 0.20 30.0 0.322 22.9 5.5 12.4 30.0
#NR94jab Northridge 1994 6.7 Bell Gardens - Jaboneria 46.6 D reverse-slip 0.13 30.0 0.068 7.6 2.5 20.1 35.0
#NR94lh1 Northridge 1994 6.7 Lake Hughes #1 # 36.3 D reverse-slip 0.12 23.0 0.087 9.4 3.7 13.9 32.0
#NR94loa Northridge 1994 6.7 Lawndale - Osage Ave 42.4 D reverse-slip 0.13 30.0 0.152 8.0 2.6 23.3 40.0
#NR94lv2 Northridge 1994 6.7 Leona Valley #2 # 37.7 D reverse-slip 0.20 23.0 0.063 7.2 1.6 12.5 32.0
#NR94php Northridge 1994 6.7 Palmdale - Hwy 14 & Palmdale # 43.6 D reverse-slip 0.20 46.0 0.067 16.9 8.0 18.2 60.0
#NR94pic Northridge 1994 6.7 LA - Pico & Sentous # 32.7 D reverse-slip 0.20 46.0 0.186 14.3 2.4 14.8 40.0
#NR94sor Northridge 1994 6.7 West Covina - S. Orange Ave 54.1 D reverse-slip 0.20 30.0 0.063 5.9 1.3 19.3 36.5
#NR94sse Northridge 1994 6.7 Terminal Island - S Seaside 60.0 D reverse-slip 0.13 30.0 0.194 12.1 2.3 13.4 35.0
#NR94ver Northridge 1994 6.7 LA - E Vernon Ave 39.3 D reverse-slip 0.10 30.0 0.153 10.1 1.8 15.9 30.0

# -------------------------------------------------------------------------------------------------------------------------end of GROUND MOTIONS

LibReadSMDFile.tcl

# READSDMFILE.TCL
# ------------------------------------------------------------------------------------------------------------read gm input format
#
# Written: MHS
# Date: July 2000
#
# A procedure which parses a ground motion record from the PEER
# strong motion database by finding dt in the record header, then
# echoing data values to the output file.
#
# Formal arguments
# inFilename -- file which contains PEER strong motion record
# outFilename -- file to be written in format G3 can read
# dt -- time step determined from file header
#
# Assumptions
# The header in the PEER record is, e.g., formatted as follows:
# PACIFIC ENGINEERING AND ANALYSIS STRONG-MOTION DATA
# IMPERIAL VALLEY 10/15/79 2319, EL CENTRO ARRAY 6, 230
# ACCELERATION TIME HISTORY IN UNITS OF G
# NPTS= 3930, DT= .00500 SEC

proc ReadSMDFile {inFilename outFilename dt} {

# Pass dt by reference
upvar $dt DT

# Open the input file and catch the error if it can't be read
if [catch {open $inFilename r} inFileID] {
puts stderr "Cannot open $inFilename for reading"
} else {
# Open output file for writing
set outFileID [open $outFilename w]

# Flag indicating dt is found and that ground motion
# values should be read -- ASSUMES dt is on last line
# of header!!!
set flag 0

# Look at each line in the file
foreach line [split [read $inFileID] \n] {

if {[llength $line] == 0} {
# Blank line --> do nothing
continue
} elseif {$flag == 1} {
# Echo ground motion values to output file
puts $outFileID $line
} else {
# Search header lines for dt
foreach word [split $line] {
# Read in the time step
if {$flag == 1} {
set DT $word
break
}
# Find the desired token and set the flag
if {[string match $word "DT="] == 1} {
set flag 1
}
}
}
}

# Close the output file
close $outFileID

# Close the input file
close $inFileID
}
}

LibSDOFanalysis.tcl

proc SDOFanalysis {omega zdamp GFile Gfact} {
source LibUnits.tcl
set DtAnalysis [expr 0.005*$sec]; # time-step Dt for lateral analysis
set TmaxGround [expr 50*$sec]; # maximum duration of ground-motion analysis
set Nsteps [expr int($TmaxGround/$DtAnalysis)];
set M 1;
set K [expr $M*pow($omega,2)];

wipe
model basic -ndm 1 -ndf 1
node 1 0.0; fix 1 1
node 2 0.0; mass 2 $M
uniaxialMaterial Elastic 1 $K
element zeroLength 1 1 2 -mat 1 -dir 1
set Bdamp [expr 2.0*$zdamp/$omega]; # Rayleigh damping factor

recorder Node -file tmpNodalDisp.out -node 2 -dof 1 disp ; # not needed
recorder EnvelopeNode -file tmpEnvNodalDisp.out -node 2 -dof 1 disp; # negMax, posMax, absMax
recorder EnvelopeNode -file tmpEnvNodalVel.out -node 2 -dof 1 vel; # negMax, posMax, absMax
recorder EnvelopeNode -file tmpEnvNodalAccel.out -node 2 -dof 1 accel; # negMax, posMax, absMax

# read a PEER strong motion database file, extracts dt from the header and converts the file
# to the format G3 expects for uniform ground motions
source LibReadSMDFile.tcl; # read in procedure difinition
set outFile GMfiles/$GFile.g3; # define output files
ReadSMDFile GMfiles/$GFile.th $outFile dt; # call procedure to convert the ground-motion file
set GMfatt [expr $g*$Gfact]; # data in input file is factor of g
set Gaccel "Series -dt $dt -filePath $outFile -factor $GMfatt"; # time series information
pattern UniformExcitation 1 1 -accel $Gaccel; # create uniform excitation
constraints Plain
algorithm Linear
system ProfileSPD
numberer Plain
set gamma 0.5; # gamma value for newmark integration
set beta 0.25; # beta value for newmark integration
integrator Newmark $gamma $beta
rayleigh $Bdamp 0 0 0

analysis Transient
analyze $Nsteps $DtAnalysis
wipe
}; # END of Analysis procedure

LibUnits.tcl

# ------------------------------------------------------------------------------------------------------------------------------------------------------------- unit definition
# UNITS.tcl -- define system of units used in the tcl script
# OpenSees output will be in the basic units
# by Silvia Mazzoni, 2001
#

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

# -------------------------------------------------------------------------------------------------------------------------end of unit definition

spectra.m

clear
close all
cd C:\Users\AAprojects\OpenSees\Examples\spectrum\spettra

in =1; sec = 1; kip = 1; % basic units
ft =12*in; ksi = kip/in^2; psi = ksi/1000; g=32.2*ft/sec^2; pi = 2*asin(1);

% gm=50 is funky (high Sd) at T>2sec: LP89svl Loma Prieta 1989 6.9 Sunnyvale - Colton Ave. 28.8 D reverse-oblique 0.10 40.0 0.207 37.3 19.1 21.2 39.3
%small magnitude, small distance
fill1 = ['IV79cal'; 'IV79chi'; 'IV79e01'; 'IV79e12'; 'IV79e13'; 'IV79qkp'; 'IV79wsm'; 'LV80kod'; 'LV80srm'; 'MH84agw'; 'MH84g02'; 'MH84g03'; 'MH84gmr'; 'PM73phn'; 'PS86psa'; 'WN87cas'; 'WN87cat'; 'WN87flo'; 'WN87w70'; 'WN87wat'];
%small magnitude, large distance
fill2 = ['BO42elc'; 'CO83c05'; 'CO83c08'; 'IV79cc4'; 'IV79cmp'; 'IV79dlt'; 'IV79nil'; 'IV79pls'; 'IV79vct'; 'LV80stp'; 'MH84cap'; 'MH84hch'; 'MH84sjb'; 'PS86h06'; 'PS86ino'; 'WN87bir'; 'WN87cts'; 'WN87har'; 'WN87sse'; 'WN87stc'];
%large magnitude, small distance
fill3 = ['IV40elc'; 'LD92yer'; 'LP89agw'; 'LP89cap'; 'LP89g03'; 'LP89g04'; 'LP89gmr'; 'LP89hch'; 'LP89hda'; 'LP89svl'; 'NR94cnp'; 'NR94far'; 'NR94fle'; 'NR94glp'; 'NR94hol'; 'NR94stc'; 'SF71pel'; 'SH87bra'; 'SH87icc'; 'SH87wsm'];
%large magnitude, large distance
fill4 = ['BM68elc'; 'LD92ind'; 'LD92psa'; 'LP89a2e'; 'LP89fms'; 'LP89hvr'; 'LP89sjw'; 'LP89slc'; 'NR94bad'; 'NR94cas'; 'NR94cen'; 'NR94jab'; 'NR94lh1'; 'NR94loa'; 'NR94lv2'; 'NR94php'; 'NR94pic'; 'NR94sor'; 'NR94sse'; 'NR94ver'];
fill = [fill1
fill2
fill3
fill4];
iSMSR=1:20;
iSMLR=21:40;
iLMSR=41:60;
iLMLR=61:80;

%fill = ['IV79cal']; % temporarily plot only one ground motion
[Ngm,junk]=size(fill);

DFileTXT = 'SdSpett';
VFileTXT = 'SvSpett';
AFileTXT = 'SaSpett';
outTXT = '.out';
load Periods.out;
[NT,junk]=size(Periods);
iiT=1:NT;
load Dampings.out;
[NZ,junk]=size(Dampings);
iiZ=1:NZ;

for igm = 1:Ngm
DFilex = [DFileTXT fill(igm,: ) outTXT];
DVett=load(DFilex);
VFilex = [VFileTXT fill(igm,: ) outTXT];
VVett=load(VFilex);
AFilex = [AFileTXT fill(igm,: ) outTXT];
AVett=load(AFilex);

for iT=1:NT
omega = 2*pi/Periods(iT);
T(igm,iT) =DVett(iT,1);
for iZ=1:NZ
Sd(igm,iZ,iT)=DVett(iT,iZ+1);
Sv(igm,iZ,iT)=VVett(iT,iZ+1);
Sa(igm,iZ,iT)=AVett(iT,iZ+1);
% Sv(igm,iZ,iT)=omega*Sd(igm,iZ,iT);
% Sa(igm,iZ,iT)=omega*omega*Sd(igm,iZ,iT);
end
end
end

% Mean spectra
Sdmean(iiZ,iiT)=mean(Sd(:,iiZ,iiT));
Svmean(iiZ,iiT)=mean(Sv(:,iiZ,iiT));
Samean(iiZ,iiT)=mean(Sa(:,iiZ,iiT));

% make Ngm figures with spectra
for igm=1:Ngm
figure
for iZ=1:NZ
SD(iiT)=Sd(igm,iZ,iiT);
SV(iiT)=Sv(igm,iZ,iiT);
SA(iiT)=Sa(igm,iZ,iiT);
subplot(2,2,1); hold on; grid on;
plot(Periods,SD);
subplot(2,2,2); hold on; grid on;
plot(Periods,SV);
subplot(2,2,3); hold on; grid on;
plot(Periods,SA/g);
subplot(2,2,4); hold on; grid on;
plot(SD,SA/g);
subplot(2,2,1); title('Response Spectra');
xlabel('Period(sec)'); ylabel('Spectral Displacement (in)')
subplot(2,2,2); title('Response Spectra');
xlabel('Period(sec)'); ylabel('Spectral Velocity (in/sec)')
subplot(2,2,3); title('Response Spectra');
xlabel('Period(sec)'); ylabel('Spectral Acceleration (g)')
subplot(2,2,4); title('Response Spectra')
xlabel('Spectral Displacement (in)'); ylabel('Spectral Acceleration (g)')
end
end

% make 1figures with MEAN spectra
figure
for iZ=1:NZ
SD(iiT)=Sdmean(iZ,iiT);
SV(iiT)=Svmean(iZ,iiT);
SA(iiT)=Samean(iZ,iiT);
subplot(2,2,1); hold on; grid on;
plot(Periods,SD);
subplot(2,2,2); hold on; grid on;
plot(Periods,SV);
subplot(2,2,3); hold on; grid on;
plot(Periods,SA/g);
subplot(2,2,4); hold on; grid on;
plot(SD,SA/g);
subplot(2,2,1); title('MEAN Response Spectra');
xlabel('Period(sec)'); ylabel('MEAN Spectral Displacement (in)')
subplot(2,2,2); title('MEAN Response Spectra');
xlabel('Period(sec)'); ylabel('MEAN Spectral Velocity (in/sec)')
subplot(2,2,3); title('MEAN Response Spectra');
xlabel('Period(sec)'); ylabel('MEAN Spectral Acceleration (g)')
subplot(2,2,4); title('MEAN Response Spectra')
xlabel('MEAN Spectral Displacement (in)'); ylabel('MEAN Spectral Acceleration (g)')
end

Yoshi · Post by **Yoshi** » Sun May 15, 2005 10:54 pm

Hi Silvia,
Thank you for your useful scripts. I recommend that these scripts will be incorpolated into the examples manual. Here I report one bug and one request.

silvia wrote:All,
LibSDOFanalysis.tcl
set Bdamp [expr 2.0*$zdamp/$omega]; # Rayleigh damping factor

Bdamp is set to [expr 2.0*$zdamp*$omega], isn't it?

And I have one request on the acceleration spectra. The calculated spectra seems to be the spectra with relative acceleration value. But we usually use the spectra of absolute acceleration. Is it possible to generate the absolute acceleration spectra?

silvia · Post by **silvia** » Wed May 18, 2005 10:41 am

you are right about the wrong value of beta, i just moved the K-prop factor over to the M-prop factor without chaning it.
thank you

about absolute accelearation vs. relative, you can always add the base acceleration, I'll look into it.

Yoshi · Post by **Yoshi** » Thu May 19, 2005 2:25 am

silvia wrote: about absolute accelearation vs. relative, you can always add the base acceleration, I'll look into it.

It seems not easy to calculate the absolute accel. spectra efficiently.
1) The EnvelopNodalRecorder cannot record the time corresponding to the recorded value. Therefore the absolute value cannot be calculated after the analysis. Of cource the max of the relative accel. response is not the max of the absolute accel. response...
2) Since all time history can be recorded by a NodalRecorder, we can calculate the max value of abs. accel.. But if we calculate it every time after one combination of (T, h), the computational cost becomes very high...
3) In the present OpenSees, all responses are stored in the relative value if the UniformExcitation is used. If only the abs. accel. spectra is needed, the MultipleSupport can be used because all responses are stored in the absolute value (I comfirmed that the abs. accel. spectra was obtained although many warnings were displayed.).

In terms of the acceleration, the absolute value is meaningful. As far as I know, there is no method to get the abs. accel. response directly in not only the interpreter but the framework. Could you ask Frank whether it is possible to make the method or the recorder?

It must be useful.

silvia · Post by **silvia** » Thu May 19, 2005 6:30 am

all very good points. i'll look into it!

hachem · Post by **hachem** » Thu May 19, 2005 8:54 am

I can't see why we can't have a recorder that will return the time of occurence of the maximum. This will be extremely useful for other things too.

silvia wrote:all very good points. i'll look into it!

silvia · Post by **silvia** » Fri May 20, 2005 6:40 am

still,
the maximum relative acceleration does not have to coincide with the maximum absolute acceleration. we just need another argument for the envelopeNode recorder with maximum absolute.
this raises another question in dynamics, for me. are the disp, velo and accel spectra all calculated from the response independently, or is one calculated from the response and then the others are factors of that value, multiplying or dividing by 2pi and T??

kryan · Post by **kryan** » Wed Oct 12, 2005 10:52 pm

I would like to make another plug for having the ability to record absolute (or total) acceleration as an option for the Node or EnvelopeNode recorder. Relative acceleration is pretty much useless. I just assumed that the recorder output absolute acceleration and found out, much to my dismay, that it did not. If anyone has come up with a clever solution for how to get the absolute acceleration by manipulation in .tcl script I would like to see it.

Regarding the dynamics questions that were raised, pseudo-acceleration and pseudo-velocity spectra can be computed directly from the relative displacement spectra by the usual relations (A = wn*V = wn^2*D). Relative velocity, relative acceleration, and absolute acceleration spectra must be computed independently. Although in an undamped system, pseudo-acceleration equals absolute acceleration.

silvia · Post by **silvia** » Thu Oct 13, 2005 6:33 pm

wait a minute!!! I am sure the output is absolute acceleration!

Goulet · Post by **Goulet** » Wed Feb 28, 2007 2:02 pm

Dear Silvia,

I just found this thread and I do not understand your Beta term for the Rayleigh damping. It WAS correct to have Beta=2*Xi/Omega, no? I refer to user Yoshi's comment, May 15 2005.

Cheers!

silvia · Post by **silvia** » Wed Feb 28, 2007 3:55 pm

look at Chopra's book and let me know what to go with.

fatooshka · Post by **fatooshka** » Fri Nov 30, 2007 4:32 am

Hi,how can i do response spectrum analysis in OS after finding response spectra?

zhenxinli1999 · Post by **zhenxinli1999** » Fri Apr 25, 2008 8:19 pm

Dear Silvia
About the damps I found large difference according to K-prop factor and M-prop factor. Which one should I use?
When T>1s acceleration values uing M-prop factor are larger when uing the other softwares, also larger than the elastic response spectrum. How to deal with this problem?

Thank you!

silvia · Post by **silvia** » Sat Apr 26, 2008 6:09 am

it depends on how you are trying to model it, it is up to you.
see what the other programs use.
when T is large, it means that K is low. so K-proportional damping will be low.

zhenxinli1999 · Post by **zhenxinli1999** » Sat Apr 26, 2008 6:27 pm

Thank you for reply
I use your useful scripts to calculate elastic response spectra for a series of ground motions.
The ground motion I used is got from elastic design spectrum. With your scripts elastic response spectrum for this ground motion is not compatible with original spectrum. Results from other software is good. I am confused......
According to K-prop factor and M-prop factor the results(accelaration, velocity or displacement) are greatly different. Why?

[quote="when T is large, it means that K is low. so K-proportional damping will be low.[/quote]

In your scripts only one damp parameter is used. How to deal with the above problem?
In general how to define damping to calculate elastic response spectrum for a ground motion?

The OpenSees Community

Elastic Response Spectra Script

Elastic Response Spectra Script

Re: Elastic Response Spectra Script

Recorder

Absolute acceleration recorder

Rayleigh Beta