1D Site Response Analysis-I need help

[koweshuns]

Hello,
I'm very new to OpenSees and i'm trying to model a linear undamped elastic soil on a rigid bedrock as described in Kramer, S.L. (1996) book (Geotechnical Earthquake Engineering. Prentice Hall, Upper Saddle River, NJ).

I used the code found on http://opensees.berkeley.edu/wiki/index ... _Analysis) with slight modifications to reflect the rigid bedrock and zero damping.
When i run the code i don't get the same results which is provided by the analytical solution and i guess i'm doing something wrong. Kindly find below my code,
thanks for your help.
.....................................................................................................................


wipe

#------------------------------------------------------------------------------------------

-
# 1. DEFINE ANALYSIS PARAMETERS
#------------------------------------------------------------------------------------------

-

#---SOIL GEOMETRY
# thicknesses of soil profile (m)
set soilThick 40.0

#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
set rho 1.7
# soil shear wave velocity (m/s)
set Vs 250.0
# soil shear modulus (kPa)
set G [expr $rho*$Vs*$Vs]
# poisson's ratio of soil
set nu 0.0
# soil elastic modulus (kPa)
set E [expr 2*$G*(1+$nu)]

#---GROUND MOTION PARAMETERS
# time step in ground motion record
set motionDT 0.005
# number of steps in ground motion record
set motionSteps 7990

#---WAVELENGTH PARAMETERS
# highest frequency desired to be well resolved (Hz)
set fMax 50.0
# wavelength of highest resolved frequency
set wave [expr $Vs/$fMax]
# number of elements per one wavelength
set nEle 20

#---ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25

#------------------------------------------------------------------------------------------

-
# 2. DEFINE MESH GEOMETRY
#------------------------------------------------------------------------------------------

-

# trial vertical element size
set hTrial [expr $wave/$nEle]

# trial number of elements
set nTrial [expr $soilThick/$hTrial]
# round up if not an integer number of elements
if { $nTrial > [expr floor($soilThick/$hTrial)] } {
set numEleY [expr int(floor($soilThick/$hTrial)+1)]
} else {
set numEleY [expr int($nTrial)]
}
puts "number of verticl elements: $numEleY"

# final vertical size of elements (m)
set sizeEleY [expr $soilThick/$numEleY]
puts "vertical size of elements: $sizeEleY"

# define horizontal size of elements (m)
set sizeEleX $sizeEleY
puts "horizontal size of elements: $sizeEleX"

# number of nodes in vertical direction
set numNodeY [expr 2*($numEleY + 1)]

#------------------------------------------------------------------------------------------

-
# 3. DEFINE NODES FOR SOIL ELEMENTS
#------------------------------------------------------------------------------------------

-

# soil nodes are created in 2 dimensions, with 2 translational dof
model BasicBuilder -ndm 2 -ndf 2

set yCoord 0.0
set count 0
# loop over nodes
for {set j 1} {$j <= $numNodeY} {incr j 2} {
node $j 0.0 [expr $yCoord + $count*$sizeEleY]
node [expr $j+1] $sizeEleX [expr $yCoord + $count*$sizeEleY]

set count [expr $count+1]
}
puts "Finished creating all soil nodes..."

#------------------------------------------------------------------------------------------

-
# 5. DEFINE BOUNDARY CONDITIONS AND EQUAL DOF
#------------------------------------------------------------------------------------------

-

# define fixity of base nodes
fix 1 0 1
fix 2 0 1


# define equal DOF for simple shear deformation of soil elements
for {set k 3} {$k <= $numNodeY} {incr k 2} {
equalDOF $k [expr $k+1] 1 2
}

# define equal DOF base soil nodes
equalDOF 1 2 1

puts "Finished creating all boundary conditions and equalDOF..."

#------------------------------------------------------------------------------------------

-
# 6. DEFINE SOIL MATERIALS
#------------------------------------------------------------------------------------------

nDMaterial ElasticIsotropic 1 $E $nu

puts "Finished creating all soil materials..."

#------------------------------------------------------------------------------------------

-
# 7. DEFINE SOIL ELEMENTS
#------------------------------------------------------------------------------------------

-

set wgtX 0.0
set wgtY [expr -9.81*$rho]

# loop over elements
for {set j 1} {$j <= $numEleY} {incr j 1} {

set nI [expr 2*$j - 1]
set nJ [expr $nI + 1]
set nK [expr $nI + 3]
set nL [expr $nI + 2]


element quad $j $nI $nJ $nK $nL 1.0 "PlaneStrain" 1 0.0 $rho $wgtX $wgtY
}

puts "Finished creating all soil elements..."

#------------------------------------------------------------------------------------------

-
# 8. CREATE GRAVITY RECORDERS
#------------------------------------------------------------------------------------------

-

# record nodal displacments, velocities, and accelerations at each time step
recorder Node -file Gdisplacement.out -time -nodeRange 1 $numNodeY -dof 1 2 disp
recorder Node -file Gvelocity.out -time -nodeRange 1 $numNodeY -dof 1 2 vel
recorder Node -file Gacceleration.out -time -nodeRange 1 $numNodeY -dof 1 2 accel


#------------------------------------------------------------------------------------------

-
# 9. APPLY GRAVITY LOADING
#------------------------------------------------------------------------------------------

-

constraints Transformation
test NormDispIncr 1e-5 30
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
analysis Transient

analyze 10 5.0e2

puts "Finished with elastic gravity analysis..."



#------------------------------------------------------------------------------------------

-
# 10. CREATE POST-GRAVITY RECORDERS
#------------------------------------------------------------------------------------------

-

# reset time and analysis
setTime 0.0
wipeAnalysis

# record nodal displacments, velocities, and accelerations at each time step
recorder Node -file displacement.out -time -dT $motionDT -nodeRange 1 $numNodeY -dof 1 2

disp
recorder Node -file velocity.out -time -dT $motionDT -nodeRange 1 $numNodeY -dof 1 2

vel
recorder Node -file acceleration.out -time -dT $motionDT -nodeRange 1 $numNodeY -dof 1 2

accel


puts "Finished creating all recorders..."

#------------------------------------------------------------------------------------------

-
# 11. DYNAMIC ANALYSIS
#------------------------------------------------------------------------------------------

# define acceleration time history file
set accelerationFile GilroyNo1EW.out

# timeseries object for applied force history

set mseries "Series -dt $motionDT -filePath $accelerationFile -factor 9.81"

pattern UniformExcitation 10 1 -accel $mseries

puts "Dynamic loading created..."

#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# trial analysis time step
set kTrial [expr $sizeEleX/(pow($Vs,0.5))]
# define time step and number of steps for analysis
if { $motionDT <= $kTrial } {
set nSteps $motionSteps
set dT $motionDT
} else {
set nSteps [expr floor($duration/$kTrial)+1]
set dT [expr $duration/$nSteps]
}
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"

# analysis objects
constraints Transformation
test NormDispIncr 1e-3 15
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh 0.0 0.0 0.0 0.0
analysis Transient

analyze $nSteps $dT

puts "Finished with dynamic analysis..."

wipe

[mcganncr]

Are you comparing your analysis to example 7.1 from Kramer's text? If so, the profile you are specifying in the input file is not the same, and that is most likely why you can't get your results to match up. Check your profile thickness, 10 ft does not equal 40 m. The shear wave velocity and mass density also look like they might not match, but I'm converting units in my head, so they could be fine. I don't see any other real problems in your input file, though it could also be an issue of relative vs. absolute acceleration, since you are using the uniform excitation command.

Hope that helps,

Chris McGann
University of Washington

[koweshuns]

Hello Chris,
Thanks for your reply and you are right about the units as they don't match up.I am not comparing the example from Kramer's to that of the code but rather trying to do was to modify the code(which represents a elastic bedrock) to represent a rigid bedrock. In doing that i removed the dashpot in the code and just applied the motion(accelerations) to the base of the model.

The result from this approach didn't match with the analytical solution( i wrote a matlab code using the frequency approach described in Kramer's text).
Thanks

[mcganncr]

OK, I see. Have you tried running these with a sine curve instead of the Gilroy motion? The simplicity of the solution for that case might help you pinpoint why the two analyses aren't matching up. It also could be helpful to compare both sets of results to something that you know is correct, e.g., the examples in Kramer's textbook.

Chris McGann
University of Washington