Code: Select all
# define velocity time history file
set filePath 1cycle.vel
puts "Input Motion File: $filePath"
# time step in ground motion record
set motionDT 0.02
# number of steps in ground motion record
set motionSteps 2000
set duration [expr $motionSteps * $motionDT]
puts "duration: $duration sec"
# delta time for analysis
set dT 0.001
# number of steps in analysis
set nSteps [expr int($duration / $dT)]
#CONSTRUCT RAYLEIGH DAMPING MATRIX (Match EQL in weak motion)
set pi 3.141592654
# natural frequency of soil column
set fn $fn
# critical damping ratio
set damp 0.02
# lower angular frequency
set omega1 [expr 2 * $pi * $fn]
puts "Lower Frequency = $omega1"
# upper angular frequency (5.0 is obtained from iteration)
set omega2 [expr 2 * $pi * (5.0*$fn)]
puts "Higher Frequency = $omega2"
# Rayleigh damping coefficients
set a0 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set a1 [expr 2*$damp/($omega1 + $omega2)]
puts "damping coefficients: a_0 = $a0; a_1 = $a1"
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"
#===ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25
#=============================================================================#
# 7. GRAVITY ANALYSIS #
#=============================================================================#
# Set OPENSEES analysis parameter
constraints Transformation
test NormDispIncr 1e-6 1000 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
analysis Transient
# Set material to behave elastically
for {set i 1} {$i <= $numLayers} {incr i 1} {
updateMaterialStage -material $i -stage 0
}
# RUN calculation ELASTIC gravity analysis
set startT [clock seconds]
analyze 10 5.0e2
puts "Finished with elastic gravity analysis..."
# Update material to consider elastoplastic behavior
for {set i 1} {$i <= $numLayers} {incr i 1} {
updateMaterialStage -material $i -stage 1
}
# Set FirstCall flag to 0 before plastic & dynamic analysis (PM4Sand)
set elmID 1
#loop over soil layer
for {set i 1} {$i <= $numLayers} {incr i 1} {
#loop over element within soil layer
for {set j 1} {$j <= $nElemY($i)} {incr j 1} {
setParameter -value 0 -ele $elmID FirstCall $i
puts "setParameter -value 0 -ele $elmID FirstCall $i"
set elmID [expr $elmID + 1]
}
}
# RUN calculation PLASTIC gravity analysis
analyze 100 5.0e-2
puts "Finished with plastic gravity analysis..."
#=============================================================================#
# 8. UPDATE ELEMENT PERMEABILITY VALUES FOR POST-GRAVITY ANALYSIS #
#=============================================================================#
set ctr 10000.0
# loop over elements to define parameter IDs
for {set i 1} {$i<=$nElemT} {incr i 1} {
parameter [expr int($ctr+1.0)] element $i vPerm
parameter [expr int($ctr+2.0)] element $i hPerm
puts "parameter [expr int($ctr+1.0)] element $i vPerm"
puts "parameter [expr int($ctr+2.0)] element $i hPerm"
set ctr [expr $ctr+2.0]
}
# Update Permeability value for each soil element
set ctr 10000.0
set elmID 1
#loop over soil layer
for {set i 1} {$i <= $numLayers} {incr i 1} {
#loop over element within soil layer
for {set j 1} {$j <= $nElemY($i)} {incr j 1} {
updateParameter [expr int($ctr+1.0)] $vPerm($i)
updateParameter [expr int($ctr+2.0)] $hPerm($i)
puts "updateParameter [expr int($ctr+1.0)] $vPerm($i)"
puts "updateParameter [expr int($ctr+2.0)] $hPerm($i)"
set ctr [expr $ctr+2.0]
set elmID [expr $elmID + 1]
}
}
puts "Finished updating permeabilities for dynamic analysis..."
#=============================================================================#
# 9. CREATE POST GRAVITY RECORDERS #
#=============================================================================#
# reset time and analysis
setTime 0.0
wipeAnalysis
remove recorders
# recorder time step
set recDT [expr 2*$motionDT]
# record nodal displacment, acceleration, and porepressure
# recorder Node -file acceleration.out -time -dT $recDT -nodeRange 1 $nNodeT -dof 1 accel
recorder Node -file displacement.out -time -dT $recDT -nodeRange 1 $nNodeT -dof 1 disp
# recorder Node -file velocity.out -time -dT $recDT -nodeRange 1 $nNodeT -dof 1 vel
recorder Node -file porePressure.out -time -dT $recDT -nodeRange 1 $nNodeT -dof 3 vel
# record stress and strain in the soil elements (Only 1 Gauss Point)
recorder Element -eleRange 1 $nElemT -time -dT $recDT -file Stress.out stress
recorder Element -eleRange 1 $nElemT -time -dT $recDT -file Strain.out strain
puts "Finished creating all recorders..."
#=============================================================================#
# 10. DYNAMIC ANALYSIS #
#=============================================================================#
# Reassign problem domain (2D with 3 DoF)
model BasicBuilder -ndm 2 -ndf 3
# define constant scaling factor for applied velocity
# Outcrop motion with elastic base
set scaleFactor 1.2
set cFactor [expr 2.0*$colArea*$dashpotCoeff*$scaleFactor]
# timeseries object for force history
set mSeries "Path -dt $motionDT -filePath $filePath -factor $cFactor"
# loading object
pattern Plain 10 $mSeries {
load 1 1.0 0.0 0.0
}
puts "Dynamic loading created..."
constraints Transformation
test NormDispIncr 1e-5 1000 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh $a0 $a1 0.0 0.0
analysis Transient
# RUN Non-linear Analysis with timestep reduction loop
set ok [analyze $nSteps $dT]
# if analysis fails, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
set mTime $curTime
puts "curTime: $curTime"
set curStep [expr $curTime/$dT]
puts "curStep: $curStep"
set rStep [expr ($nSteps-$curStep)*2.0]
set remStep [expr int(($nSteps-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
set ok [analyze $remStep $dT]
# if analysis fails again, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
puts "curTime: $curTime"
set curStep [expr ($curTime-$mTime)/$dT]
puts "curStep: $curStep"
set remStep [expr int(($rStep-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
analyze $remStep $dT
}
}
#=============================================================================#
# 11. DIFFUSSION ANALYSIS #
#=============================================================================#
wipeAnalysis
#CONSTRUCT RAYLEIGH DAMPING MATRIX (Match EQL in weak motion)
set pi 3.141592654
# natural frequency of soil column
set fn $fn
# critical damping ratio
set dampAft 0.9
# lower angular frequency
set omega1Aft [expr 2 * $pi * $fn]
puts "Lower Frequency = $omega1"
# upper angular frequency (5.0 is obtained from iteration)
set omega2Aft [expr 2 * $pi * (5.0*$fn)]
puts "Higher Frequency = $omega2"
# Rayleigh damping coefficients
set a0Aft [expr 2*$dampAft*$omega1Aft*$omega2Aft/($omega1Aft + $omega2Aft)]
set a1Aft [expr 2*$dampAft/($omega1Aft + $omega2Aft)]
puts "damping coefficients: a_0 = $a0; a_1 = $a1"
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"
# delta time for analysis
set dT 0.01
# number of steps in analysis
set nSteps [expr int(40 / $dT)]
# Reassign problem domain (2D with 3 DoF)
model BasicBuilder -ndm 2 -ndf 3
constraints Transformation
test NormDispIncr 1e-5 100 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh $a0Aft $a1Aft 0.0 0.0
analysis Transient
# RUN Non-linear Analysis with timestep reduction loop
set ok [analyze $nSteps $dT]
# if analysis fails, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
set mTime $curTime
puts "curTime: $curTime"
set curStep [expr $curTime/$dT]
puts "curStep: $curStep"
set rStep [expr ($nSteps-$curStep)*2.0]
set remStep [expr int(($nSteps-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
set ok [analyze $remStep $dT]
# if analysis fails again, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
puts "curTime: $curTime"
set curStep [expr ($curTime-$mTime)/$dT]
puts "curStep: $curStep"
set remStep [expr int(($rStep-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
analyze $remStep $dT
}
}
set endT [clock seconds]
puts "Finished with dynamic analysis..."